Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-2024 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "core/ActionShortcut.h"
23 : #include "core/ActionRegister.h"
24 : #include "core/PlumedMain.h"
25 : #include "core/ActionSet.h"
26 : #include "core/Group.h"
27 : #include "tools/PDB.h"
28 :
29 : namespace PLMD {
30 : namespace colvar {
31 :
32 : //+PLUMEDOC DCOLVAR DRMSD
33 : /*
34 : Calculate the distance RMSD with respect to a reference structure.
35 :
36 : To calculate the root-mean-square deviation between the atoms in two configurations
37 : you must first superimpose the two structures in some ways. Obviously, it is the internal vibrational
38 : motions of the structure - i.e. not the translations and rotations - that are interesting. However,
39 : aligning two structures by removing the translational and rotational motions is not easy. Furthermore,
40 : in some cases there can be alignment issues caused by so-called frame-fitting problems. It is thus
41 : often easier to calculate the distances between all the pairs of atoms. The distance
42 : between the two structures, ${\bf X}^a$ and ${\bf X}^b$ can then be measured as:
43 :
44 : $$
45 : d({\bf X}^A, {\bf X}^B) = \sqrt{\frac{1}{N(N-1)} \sum_{i \ne j} [ d({\bf x}_i^a,{\bf x}_j^a) - d({\bf x}_i^b,{\bf x}_j^b) ]^2}
46 : $$
47 :
48 : where $N$ is the number of atoms and $d({\bf x}_i,{\bf x}_j)$ represents the distance between
49 : atoms $i$ and $j$. Clearly, this representation of the configuration is invariant to translation and rotation.
50 : However, it can become expensive to calculate when the number of atoms is large. This can be resolved
51 : within the DRMSD colvar by setting `LOWER_CUTOFF` and `UPPER_CUTOFF`. These keywords ensure that only
52 : pairs of atoms that are within a certain range are incorporated into the above sum.
53 :
54 : In PDB files the atomic coordinates and box lengths should be in Angstroms unless
55 : you are working with natural units. If you are working with natural units then the coordinates
56 : should be in your natural length unit. [Click here](http://www.wwpdb.org/docs.html) for more details on the PDB file format.
57 :
58 : ## Examples
59 :
60 : The following tells plumed to calculate the distance RMSD between
61 : the positions of the atoms in the reference file and their instantaneous
62 : position. Only pairs of atoms whose distance in the reference structure is within
63 : 0.1 and 0.8 nm are considered.
64 :
65 : ```plumed
66 : #SETTINGS INPUTFILES=regtest/basic/rt-drmsd/test1.pdb
67 : DRMSD REFERENCE=regtest/basic/rt-drmsd/test1.pdb LOWER_CUTOFF=0.1 UPPER_CUTOFF=0.8
68 : ```
69 :
70 : The following tells plumed to calculate a DRMSD value for a pair of molecules.
71 :
72 : ```plumed
73 : #SETTINGS INPUTFILES=regtest/basic/rt-drmsd/test0.pdb
74 : DRMSD REFERENCE=regtest/basic/rt-drmsd/test0.pdb LOWER_CUTOFF=0.1 UPPER_CUTOFF=0.8 TYPE=INTER-DRMSD
75 : ```
76 :
77 : Notice that in the input reference file (which you can see by clicking on regtest/basic/rt-drmsd/test0.pdb )
78 : the atoms in each of the two molecules are separated by a TER
79 : command as shown below.
80 :
81 : In this example the INTER-DRMSD type ensures that the set of distances from which the final
82 : quantity is computed involve one atom from each of the two molecules. If this is replaced
83 : by INTRA-DRMSD then only those distances involving pairs of atoms that are both in the same
84 : molecule are computed.
85 :
86 : */
87 : //+ENDPLUMEDOC
88 :
89 : class DRMSD : public ActionShortcut {
90 : public:
91 : static void registerKeywords(Keywords& keys);
92 : explicit DRMSD(const ActionOptions&);
93 : };
94 :
95 : PLUMED_REGISTER_ACTION(DRMSD,"DRMSD")
96 :
97 16 : void DRMSD::registerKeywords( Keywords& keys ) {
98 16 : ActionShortcut::registerKeywords( keys );
99 16 : keys.add("compulsory","REFERENCE","a file in pdb format containing the reference structure and the atoms involved in the CV.");
100 16 : keys.add("optional","LOWER_CUTOFF","only pairs of atoms further than LOWER_CUTOFF are considered in the calculation.");
101 16 : keys.add("optional","UPPER_CUTOFF","only pairs of atoms closer than UPPER_CUTOFF are considered in the calculation.");
102 16 : keys.add("compulsory","TYPE","DRMSD","what kind of DRMSD would you like to calculate. You can use either the normal DRMSD involving all the distances between "
103 : "the atoms in your molecule. Alternatively, if you have multiple molecules you can use the type INTER-DRMSD "
104 : "to compute DRMSD values involving only those distances between the atoms at least two molecules or the type INTRA-DRMSD "
105 : "to compute DRMSD values involving only those distances between atoms in the same molecule");
106 16 : keys.addFlag("SQUARED",false,"This should be setted if you want MSD instead of RMSD ");
107 16 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
108 : // This is just ignored in reality which is probably bad
109 16 : keys.addFlag("NUMERICAL_DERIVATIVES",false,"calculate the derivatives for these quantities numerically");
110 32 : keys.setValueDescription("scalar/vector","the DRMSD distance between the instantaneous structure and the reference structure");
111 16 : keys.needsAction("SUM");
112 16 : keys.needsAction("DISTANCE");
113 16 : keys.needsAction("CONSTANT");
114 16 : keys.needsAction("EUCLIDEAN_DISTANCE");
115 16 : keys.needsAction("CUSTOM");
116 16 : }
117 :
118 13 : DRMSD::DRMSD( const ActionOptions& ao ):
119 : Action(ao),
120 13 : ActionShortcut(ao) {
121 : // Read in the reference configuration
122 : std::string reference;
123 13 : parse("REFERENCE",reference);
124 : // First bit of input for the instantaneous distances
125 : bool numder;
126 26 : parseFlag("NUMERICAL_DERIVATIVES",numder);
127 : double fake_unit=0.1;
128 13 : FILE* fp2=fopen(reference.c_str(),"r");
129 : bool do_read=true;
130 : unsigned nframes=0;
131 65 : while( do_read ) {
132 54 : PDB mypdb;
133 54 : do_read=mypdb.readFromFilepointer(fp2,false,fake_unit);
134 54 : if( !do_read && nframes>0 ) {
135 : break ;
136 : }
137 53 : nframes++;
138 54 : }
139 12 : fclose(fp2);
140 :
141 : // Get cutoff information
142 12 : double lcut=0;
143 25 : parse("LOWER_CUTOFF",lcut);
144 : std::string drmsd_type;
145 12 : parse("TYPE",drmsd_type);
146 12 : double ucut=std::numeric_limits<double>::max();
147 12 : parse("UPPER_CUTOFF",ucut);
148 : bool nopbc;
149 24 : parseFlag("NOPBC",nopbc);
150 : std::string pbc_str;
151 12 : if(nopbc) {
152 : pbc_str="NOPBC";
153 : }
154 : // Open the pdb file
155 12 : FILE* fp=fopen(reference.c_str(),"r");
156 : do_read=true;
157 12 : if(!fp) {
158 0 : error("could not open reference file " + reference );
159 : }
160 : unsigned n=0;
161 12 : std::string allpairs="";
162 : std::vector<std::pair<unsigned,unsigned> > upairs;
163 : std::vector<std::string> refvals;
164 65 : while ( do_read ) {
165 54 : PDB mypdb;
166 54 : do_read=mypdb.readFromFilepointer(fp,false,fake_unit);
167 54 : if( !do_read && n>0 ) {
168 : break ;
169 : }
170 53 : std::vector<Vector> pos( mypdb.getPositions() );
171 53 : unsigned nn=1;
172 53 : if( pos.size()==0 ) {
173 0 : error("read no atoms from file named " + reference );
174 : }
175 : // This is what we do for the first frame
176 53 : if( n==0 ) {
177 12 : std::vector<AtomNumber> atoms( mypdb.getAtomNumbers() );
178 12 : if( drmsd_type=="DRMSD" ) {
179 102 : for(unsigned i=0; i<atoms.size()-1; ++i) {
180 : std::string istr;
181 92 : Tools::convert( atoms[i].serial(), istr );
182 671 : for(unsigned j=i+1; j<atoms.size(); ++j) {
183 : std::string jstr;
184 579 : Tools::convert( atoms[j].serial(), jstr );
185 579 : double distance = delta( pos[i], pos[j] ).modulo();
186 579 : if( distance < ucut && distance > lcut ) {
187 : std::string num;
188 386 : Tools::convert( nn, num );
189 386 : nn++;
190 : // Add this pair to list of pairs
191 386 : upairs.push_back( std::pair<unsigned,unsigned>(i,j) );
192 : // Add this distance to list of reference values
193 : std::string dstr;
194 386 : Tools::convert( distance, dstr );
195 386 : refvals.push_back( dstr );
196 : // Calculate this distance
197 386 : if( nframes==1 ) {
198 616 : allpairs += " ATOMS" + num + "=" + istr + "," + jstr;
199 : } else {
200 156 : readInputLine( getShortcutLabel() + "_d" + num + ": DISTANCE ATOMS=" + istr + "," + jstr + " " + pbc_str );
201 : }
202 : }
203 : }
204 : }
205 : } else {
206 2 : unsigned nblocks = mypdb.getNumberOfAtomBlocks();
207 2 : std::vector<unsigned> blocks( 1 + nblocks );
208 2 : if( nblocks==1 ) {
209 0 : blocks[0]=0;
210 0 : blocks[1]=atoms.size();
211 : } else {
212 2 : blocks[0]=0;
213 6 : for(unsigned i=0; i<nblocks; ++i) {
214 4 : blocks[i+1]=mypdb.getAtomBlockEnds()[i];
215 : }
216 : }
217 2 : if( drmsd_type=="INTRA-DRMSD" ) {
218 3 : for(unsigned i=0; i<nblocks; ++i) {
219 7 : for(unsigned iatom=blocks[i]+1; iatom<blocks[i+1]; ++iatom) {
220 : std::string istr;
221 5 : Tools::convert( atoms[iatom].serial(), istr );
222 14 : for(unsigned jatom=blocks[i]; jatom<iatom; ++jatom) {
223 : std::string jstr;
224 9 : Tools::convert( atoms[jatom].serial(), jstr );
225 9 : double distance = delta( pos[iatom], pos[jatom] ).modulo();
226 9 : if(distance < ucut && distance > lcut ) {
227 : std::string num;
228 9 : Tools::convert( nn, num );
229 9 : nn++;
230 : // Add this pair to list of pairs
231 9 : upairs.push_back( std::pair<unsigned,unsigned>(iatom,jatom) );
232 : // Add this distance to list of reference values
233 : std::string dstr;
234 9 : Tools::convert( distance, dstr );
235 9 : refvals.push_back( dstr );
236 : // Calculate this distance
237 9 : if( nframes==1 ) {
238 18 : allpairs += " ATOMS" + num + "=" + istr + "," + jstr;
239 : } else {
240 0 : readInputLine( getShortcutLabel() + "_d" + num + ": DISTANCE ATOMS=" + istr + "," + jstr + " " + pbc_str );
241 : }
242 : }
243 : }
244 : }
245 : }
246 1 : } else if( drmsd_type=="INTER-DRMSD" ) {
247 2 : for(unsigned i=1; i<nblocks; ++i) {
248 2 : for(unsigned j=0; j<i; ++j) {
249 4 : for(unsigned iatom=blocks[i]; iatom<blocks[i+1]; ++iatom) {
250 : std::string istr;
251 3 : Tools::convert( atoms[iatom].serial(), istr );
252 15 : for(unsigned jatom=blocks[j]; jatom<blocks[j+1]; ++jatom) {
253 : std::string jstr;
254 12 : Tools::convert( atoms[jatom].serial(), jstr );
255 12 : double distance = delta( pos[iatom], pos[jatom] ).modulo();
256 12 : if(distance < ucut && distance > lcut ) {
257 : std::string num;
258 12 : Tools::convert( nn, num );
259 12 : nn++;
260 : // Add this pair to list of pairs
261 12 : upairs.push_back( std::pair<unsigned,unsigned>(iatom,jatom) );
262 : // Add this distance to list of reference values
263 : std::string dstr;
264 12 : Tools::convert( distance, dstr );
265 12 : refvals.push_back( dstr );
266 : // Calculate this distance
267 12 : if( nframes==1 ) {
268 24 : allpairs += " ATOMS" + num + "=" + istr + "," + jstr;
269 : } else {
270 0 : readInputLine( getShortcutLabel() + "_d" + num + ": DISTANCE ATOMS=" + istr + "," + jstr + " " + pbc_str );
271 : }
272 : }
273 : }
274 : }
275 : }
276 : }
277 : } else {
278 0 : plumed_merror( drmsd_type + " is not valid input to TYPE keyword");
279 : }
280 : }
281 : // This is for every subsequent frame
282 : } else {
283 3239 : for(unsigned i=0; i<refvals.size(); ++i) {
284 : std::string dstr;
285 3198 : Tools::convert( delta( pos[upairs[i].first], pos[upairs[i].second] ).modulo(), dstr );
286 6396 : refvals[i] += "," + dstr;
287 : }
288 : }
289 53 : n++;
290 54 : }
291 : // Now create values that hold all the reference distances
292 12 : fclose(fp);
293 :
294 12 : if( nframes==1 ) {
295 22 : readInputLine( getShortcutLabel() + "_d: DISTANCE" + allpairs + " " + pbc_str );
296 11 : std::string refstr = refvals[0];
297 329 : for(unsigned i=1; i<refvals.size(); ++i) {
298 636 : refstr += "," + refvals[i];
299 : }
300 22 : readInputLine( getShortcutLabel() + "_ref: CONSTANT VALUES=" + refstr );
301 22 : readInputLine( getShortcutLabel() + "_diffs: CUSTOM ARG=" + getShortcutLabel() + "_d," + getShortcutLabel() + "_ref FUNC=(x-y)*(x-y) PERIODIC=NO");
302 22 : readInputLine( getShortcutLabel() + "_u: SUM ARG=" + getShortcutLabel() + "_diffs PERIODIC=NO");
303 : } else {
304 : std::string arg_str1, arg_str2;
305 79 : for(unsigned i=0; i<refvals.size(); ++i ) {
306 : std::string inum;
307 78 : Tools::convert( i+1, inum );
308 156 : readInputLine( getShortcutLabel() + "_ref" + inum + ": CONSTANT VALUES=" + refvals[i] );
309 78 : if( i==0 ) {
310 2 : arg_str1 = getShortcutLabel() + "_d" + inum;
311 2 : arg_str2 = getShortcutLabel() + "_ref" + inum;
312 : } else {
313 154 : arg_str1 += "," + getShortcutLabel() + "_d" + inum;
314 154 : arg_str2 += "," + getShortcutLabel() + "_ref" + inum;
315 : }
316 : }
317 : // And calculate the euclidean distances between the true distances and the references
318 2 : readInputLine( getShortcutLabel() + "_u: EUCLIDEAN_DISTANCE SQUARED ARG1=" + arg_str1 + " ARG2=" + arg_str2 );
319 : }
320 : // And final value
321 : std::string nvals;
322 12 : Tools::convert( refvals.size(), nvals );
323 : bool squared;
324 12 : parseFlag("SQUARED",squared);
325 12 : if( squared ) {
326 6 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_u FUNC=x/" + nvals + " PERIODIC=NO");
327 : } else {
328 18 : readInputLine( getShortcutLabel() + "_2: CUSTOM ARG=" + getShortcutLabel() + "_u FUNC=(x/" + nvals + ") PERIODIC=NO");
329 18 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_2 FUNC=sqrt(x) PERIODIC=NO");
330 : }
331 26 : }
332 :
333 : }
334 : }
335 :
336 :
|
