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 : d: DRMSD ...
68 : REFERENCE=regtest/basic/rt-drmsd/test1.pdb
69 : LOWER_CUTOFF=0.1 UPPER_CUTOFF=0.8
70 : ...
71 : ```
72 :
73 : The following tells plumed to calculate the square of the DRMSD value for a pair of molecules.
74 :
75 : ```plumed
76 : #SETTINGS INPUTFILES=regtest/basic/rt-drmsd/test0.pdb
77 : d: DRMSD ...
78 : REFERENCE=regtest/basic/rt-drmsd/test0.pdb
79 : LOWER_CUTOFF=0.1 UPPER_CUTOFF=0.8 TYPE=INTER-DRMSD
80 : SQUARED
81 : ...
82 : ```
83 :
84 : Notice that in the input reference file (which you can see by clicking on regtest/basic/rt-drmsd/test0.pdb )
85 : the atoms in each of the two molecules are separated by a TER
86 : command as shown below.
87 :
88 : In this example the INTER-DRMSD type ensures that the set of distances from which the final
89 : quantity is computed involve one atom from each of the two molecules. If this is replaced
90 : by INTRA-DRMSD then only those distances involving pairs of atoms that are both in the same
91 : molecule are computed.
92 :
93 : ## Periodic boundary conditions
94 :
95 : Notice that PLUMED does not use periodic boundary conditions when computing the reference distances. However, the instantaneous
96 : values of the distances are computed using periodic boundary conditions. If you would like not to use periodic boundary conditions
97 : when computing these instantaneous distances you use the NOPBC flag as illustrated below:
98 :
99 : ```plumed
100 : #SETTINGS INPUTFILES=regtest/basic/rt-drmsd/test1.pdb
101 : d: DRMSD ...
102 : REFERENCE=regtest/basic/rt-drmsd/test1.pdb
103 : LOWER_CUTOFF=0.1 UPPER_CUTOFF=0.8 NOPBC
104 : ...
105 : ```
106 :
107 : */
108 : //+ENDPLUMEDOC
109 :
110 : class DRMSD : public ActionShortcut {
111 : public:
112 : static void registerKeywords(Keywords& keys);
113 : explicit DRMSD(const ActionOptions&);
114 : };
115 :
116 : PLUMED_REGISTER_ACTION(DRMSD,"DRMSD")
117 :
118 19 : void DRMSD::registerKeywords( Keywords& keys ) {
119 19 : ActionShortcut::registerKeywords( keys );
120 19 : keys.add("compulsory","REFERENCE","a file in pdb format containing the reference structure and the atoms involved in the CV.");
121 19 : keys.add("optional","LOWER_CUTOFF","only pairs of atoms further than LOWER_CUTOFF are considered in the calculation.");
122 19 : keys.add("optional","UPPER_CUTOFF","only pairs of atoms closer than UPPER_CUTOFF are considered in the calculation.");
123 19 : 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 "
124 : "the atoms in your molecule. Alternatively, if you have multiple molecules you can use the type INTER-DRMSD "
125 : "to compute DRMSD values involving only those distances between the atoms at least two molecules or the type INTRA-DRMSD "
126 : "to compute DRMSD values involving only those distances between atoms in the same molecule");
127 19 : keys.addFlag("SQUARED",false,"This should be set if you want MSD instead of RMSD ");
128 19 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
129 : // This is just ignored in reality which is probably bad
130 38 : keys.addDeprecatedFlag("NUMERICAL_DERIVATIVES","");
131 38 : keys.setValueDescription("scalar/vector","the DRMSD distance between the instantaneous structure and the reference structure");
132 19 : keys.needsAction("SUM");
133 19 : keys.needsAction("DISTANCE");
134 19 : keys.needsAction("CONSTANT");
135 19 : keys.needsAction("EUCLIDEAN_DISTANCE");
136 19 : keys.needsAction("CUSTOM");
137 19 : }
138 :
139 14 : DRMSD::DRMSD( const ActionOptions& ao ):
140 : Action(ao),
141 14 : ActionShortcut(ao) {
142 : // Read in the reference configuration
143 : std::string reference;
144 14 : parse("REFERENCE",reference);
145 : // First bit of input for the instantaneous distances
146 : bool numder;
147 28 : parseFlag("NUMERICAL_DERIVATIVES",numder);
148 : double fake_unit=0.1;
149 14 : FILE* fp2=fopen(reference.c_str(),"r");
150 : bool do_read=true;
151 : unsigned nframes=0;
152 67 : while( do_read ) {
153 56 : PDB mypdb;
154 56 : do_read=mypdb.readFromFilepointer(fp2,false,fake_unit);
155 56 : if( !do_read && nframes>0 ) {
156 : break ;
157 : }
158 54 : nframes++;
159 56 : }
160 13 : fclose(fp2);
161 :
162 : // Get cutoff information
163 13 : double lcut=0;
164 28 : parse("LOWER_CUTOFF",lcut);
165 : std::string drmsd_type;
166 13 : parse("TYPE",drmsd_type);
167 13 : double ucut=std::numeric_limits<double>::max();
168 13 : parse("UPPER_CUTOFF",ucut);
169 : bool nopbc;
170 27 : parseFlag("NOPBC",nopbc);
171 : std::string pbc_str;
172 13 : if(nopbc) {
173 : pbc_str="NOPBC";
174 : }
175 : // Open the pdb file
176 13 : FILE* fp=fopen(reference.c_str(),"r");
177 : do_read=true;
178 13 : if(!fp) {
179 0 : error("could not open reference file " + reference );
180 : }
181 : unsigned n=0;
182 14 : std::string allpairs="";
183 13 : std::vector<std::pair<unsigned,unsigned> > upairs;
184 13 : std::vector<std::string> refvals;
185 67 : while ( do_read ) {
186 56 : PDB mypdb;
187 56 : do_read=mypdb.readFromFilepointer(fp,false,fake_unit);
188 56 : if( !do_read && n>0 ) {
189 : break ;
190 : }
191 54 : std::vector<Vector> pos( mypdb.getPositions() );
192 54 : unsigned nn=1;
193 54 : if( pos.size()==0 ) {
194 0 : error("read no atoms from file named " + reference );
195 : }
196 : // This is what we do for the first frame
197 54 : if( n==0 ) {
198 13 : std::vector<AtomNumber> atoms( mypdb.getAtomNumbers() );
199 13 : if( drmsd_type=="DRMSD" ) {
200 114 : for(unsigned i=0; i<atoms.size()-1; ++i) {
201 : std::string istr;
202 103 : Tools::convert( atoms[i].serial(), istr );
203 748 : for(unsigned j=i+1; j<atoms.size(); ++j) {
204 : std::string jstr;
205 645 : Tools::convert( atoms[j].serial(), jstr );
206 645 : double distance = delta( pos[i], pos[j] ).modulo();
207 645 : if( distance < ucut && distance > lcut ) {
208 : std::string num;
209 386 : Tools::convert( nn, num );
210 386 : nn++;
211 : // Add this pair to list of pairs
212 386 : upairs.push_back( std::pair<unsigned,unsigned>(i,j) );
213 : // Add this distance to list of reference values
214 : std::string dstr;
215 386 : Tools::convert( distance, dstr );
216 386 : refvals.push_back( dstr );
217 : // Calculate this distance
218 386 : if( nframes==1 ) {
219 924 : allpairs += " ATOMS" + num + "=" + istr + "," + jstr;
220 : } else {
221 390 : readInputLine( getShortcutLabel() + "_d" + num + ": DISTANCE ATOMS=" + istr + "," + jstr + " " + pbc_str );
222 : }
223 : }
224 : }
225 : }
226 : } else {
227 2 : unsigned nblocks = mypdb.getNumberOfAtomBlocks();
228 2 : std::vector<unsigned> blocks( 1 + nblocks );
229 2 : if( nblocks==1 ) {
230 0 : blocks[0]=0;
231 0 : blocks[1]=atoms.size();
232 : } else {
233 2 : blocks[0]=0;
234 6 : for(unsigned i=0; i<nblocks; ++i) {
235 4 : blocks[i+1]=mypdb.getAtomBlockEnds()[i];
236 : }
237 : }
238 2 : if( drmsd_type=="INTRA-DRMSD" ) {
239 3 : for(unsigned i=0; i<nblocks; ++i) {
240 7 : for(unsigned iatom=blocks[i]+1; iatom<blocks[i+1]; ++iatom) {
241 : std::string istr;
242 5 : Tools::convert( atoms[iatom].serial(), istr );
243 14 : for(unsigned jatom=blocks[i]; jatom<iatom; ++jatom) {
244 : std::string jstr;
245 9 : Tools::convert( atoms[jatom].serial(), jstr );
246 9 : double distance = delta( pos[iatom], pos[jatom] ).modulo();
247 9 : if(distance < ucut && distance > lcut ) {
248 : std::string num;
249 9 : Tools::convert( nn, num );
250 9 : nn++;
251 : // Add this pair to list of pairs
252 9 : upairs.push_back( std::pair<unsigned,unsigned>(iatom,jatom) );
253 : // Add this distance to list of reference values
254 : std::string dstr;
255 9 : Tools::convert( distance, dstr );
256 9 : refvals.push_back( dstr );
257 : // Calculate this distance
258 9 : if( nframes==1 ) {
259 27 : allpairs += " ATOMS" + num + "=" + istr + "," + jstr;
260 : } else {
261 0 : readInputLine( getShortcutLabel() + "_d" + num + ": DISTANCE ATOMS=" + istr + "," + jstr + " " + pbc_str );
262 : }
263 : }
264 : }
265 : }
266 : }
267 1 : } else if( drmsd_type=="INTER-DRMSD" ) {
268 2 : for(unsigned i=1; i<nblocks; ++i) {
269 2 : for(unsigned j=0; j<i; ++j) {
270 4 : for(unsigned iatom=blocks[i]; iatom<blocks[i+1]; ++iatom) {
271 : std::string istr;
272 3 : Tools::convert( atoms[iatom].serial(), istr );
273 15 : for(unsigned jatom=blocks[j]; jatom<blocks[j+1]; ++jatom) {
274 : std::string jstr;
275 12 : Tools::convert( atoms[jatom].serial(), jstr );
276 12 : double distance = delta( pos[iatom], pos[jatom] ).modulo();
277 12 : if(distance < ucut && distance > lcut ) {
278 : std::string num;
279 12 : Tools::convert( nn, num );
280 12 : nn++;
281 : // Add this pair to list of pairs
282 12 : upairs.push_back( std::pair<unsigned,unsigned>(iatom,jatom) );
283 : // Add this distance to list of reference values
284 : std::string dstr;
285 12 : Tools::convert( distance, dstr );
286 12 : refvals.push_back( dstr );
287 : // Calculate this distance
288 12 : if( nframes==1 ) {
289 36 : allpairs += " ATOMS" + num + "=" + istr + "," + jstr;
290 : } else {
291 0 : readInputLine( getShortcutLabel() + "_d" + num + ": DISTANCE ATOMS=" + istr + "," + jstr + " " + pbc_str );
292 : }
293 : }
294 : }
295 : }
296 : }
297 : }
298 : } else {
299 0 : plumed_merror( drmsd_type + " is not valid input to TYPE keyword");
300 : }
301 2 : }
302 : // This is for every subsequent frame
303 13 : } else {
304 3239 : for(unsigned i=0; i<refvals.size(); ++i) {
305 : std::string dstr;
306 3198 : Tools::convert( delta( pos[upairs[i].first], pos[upairs[i].second] ).modulo(), dstr );
307 6396 : refvals[i] += "," + dstr;
308 : }
309 : }
310 54 : n++;
311 56 : }
312 : // Now create values that hold all the reference distances
313 13 : fclose(fp);
314 :
315 13 : if( refvals.empty() ) {
316 2 : plumed_merror("no pairs of atoms are within the cutoffs for the reference structure, no distances can be computed");
317 : }
318 12 : if( nframes==1 ) {
319 33 : readInputLine( getShortcutLabel() + "_d: DISTANCE" + allpairs + " " + pbc_str );
320 11 : std::string refstr = refvals[0];
321 329 : for(unsigned i=1; i<refvals.size(); ++i) {
322 636 : refstr += "," + refvals[i];
323 : }
324 11 : readInputLine( getShortcutLabel() + "_ref: CONSTANT VALUES=" + refstr );
325 33 : readInputLine( getShortcutLabel() + "_diffs: CUSTOM ARG=" + getShortcutLabel() + "_d," + getShortcutLabel() + "_ref FUNC=(x-y)*(x-y) PERIODIC=NO");
326 33 : readInputLine( getShortcutLabel() + "_u: SUM ARG=" + getShortcutLabel() + "_diffs PERIODIC=NO");
327 : } else {
328 : std::string arg_str1, arg_str2;
329 79 : for(unsigned i=0; i<refvals.size(); ++i ) {
330 : std::string inum;
331 78 : Tools::convert( i+1, inum );
332 156 : readInputLine( getShortcutLabel() + "_ref" + inum + ": CONSTANT VALUES=" + refvals[i] );
333 78 : if( i==0 ) {
334 1 : arg_str1 = getShortcutLabel() + "_d" + inum;
335 2 : arg_str2 = getShortcutLabel() + "_ref" + inum;
336 : } else {
337 154 : arg_str1 += "," + getShortcutLabel() + "_d" + inum;
338 231 : arg_str2 += "," + getShortcutLabel() + "_ref" + inum;
339 : }
340 : }
341 : // And calculate the euclidean distances between the true distances and the references
342 3 : readInputLine( getShortcutLabel() + "_u: EUCLIDEAN_DISTANCE SQUARED ARG2=" + arg_str1 + " ARG1=" + arg_str2 );
343 : }
344 : // And final value
345 : std::string nvals;
346 12 : Tools::convert( refvals.size(), nvals );
347 : bool squared;
348 12 : parseFlag("SQUARED",squared);
349 12 : if( squared ) {
350 12 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_u FUNC=x/" + nvals + " PERIODIC=NO");
351 : } else {
352 27 : readInputLine( getShortcutLabel() + "_2: CUSTOM ARG=" + getShortcutLabel() + "_u FUNC=(x/" + nvals + ") PERIODIC=NO");
353 27 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_2 FUNC=sqrt(x) PERIODIC=NO");
354 : }
355 359 : }
356 :
357 : }
358 : }
359 :
360 :
|
