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