Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2023 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 "Colvar.h"
23 : #include "core/PlumedMain.h"
24 : #include "core/ActionRegister.h"
25 : #include "tools/RMSD.h"
26 : #include "tools/PDB.h"
27 :
28 : namespace PLMD {
29 : namespace colvar {
30 :
31 : class RMSD : public Colvar {
32 :
33 : bool squared;
34 : bool nopbc;
35 : PLMD::RMSD myrmsd;
36 : std::vector<Vector> der;
37 : public:
38 : explicit RMSD(const ActionOptions&);
39 : void calculate() override;
40 : static void registerKeywords(Keywords& keys);
41 : };
42 :
43 : //+PLUMEDOC DCOLVAR RMSD
44 : /*
45 : Calculate the RMSD with respect to a reference structure.
46 :
47 : One colvar that has been shown to be very successful in studying protein folding is the distance between the instantaneous configuration
48 : and a reference configuration - often the structure of the folded state. When the free energy of a protein is shown as a function
49 : of this collective variable there is a minima for low values of the CV, which is due to the folded state of the protein. There is
50 : then a second minima at higher values of the CV, which is the minima corresponding to the unfolded state. The aim of this colvar is
51 : thus to calculate something like:
52 :
53 : $$
54 : d(X,X') = \vert X-X' \vert
55 : $$
56 :
57 : where $X$ is the instantaneous position of all the atoms in the system and
58 : $X'$ is the positions of the atoms in some reference structure that is provided as input.
59 : $d(X,X')$ thus measures the distance all the atoms have moved away from this reference configuration.
60 : Oftentimes, it is only the internal motions of the structure - i.e. not the translations of the center of
61 : mass or the rotations of the reference frame - that are interesting. Hence, when calculating
62 : the root-mean-square deviation between the atoms in two configurations
63 : you must first superimpose the two structures in some way. At present PLUMED provides two distinct ways
64 : of performing this superposition. The first method is applied when you use TYPE=SIMPLE in the input
65 : line. This instruction tells PLUMED that the root mean square deviation is to be calculated after the
66 : positions of the geometric centers in the reference and instantaneous configurations are aligned. In
67 : other words $d(X,x')$ is to be calculated using:
68 :
69 : $$
70 : d(X,X') = \sqrt{ \sum_i \sum_\alpha^{x,y,z} \frac{w_i}{\sum_j w_j} ( X _{i,\alpha}-com _\alpha(X)-{X'} _{i,\alpha}+com _\alpha(X') )^2 }
71 : $$
72 :
73 : with
74 :
75 : $$
76 : com_\alpha (X) = \sum_i \frac{w'_{i}}{\sum_j w'_j} X _{i,\alpha}
77 : $$
78 :
79 : and
80 :
81 : $$
82 : com_\alpha(X')= \sum_i \frac{w'_{i}}{\sum_j w'_j}X' _{i,\alpha}
83 : $$
84 :
85 : Obviously, $com_\alpha(X)$ and $com_\alpha(X')$ represent the positions of the center of mass in the reference
86 : and instantaneous configurations if the weights $w'$ are set equal to the atomic masses. If the weights are all set equal to
87 : one, however, $com_\alpha(X)$ and $com_\alpha(X')$ are the positions of the geometric centers.
88 :
89 : An example input that can be used to calculate and print this RMSD distance is shown below:
90 :
91 : ```plumed
92 : #SETTINGS INPUTFILES=regtest/basic/rt19/test0.pdb
93 : rmsd0: RMSD TYPE=SIMPLE REFERENCE=regtest/basic/rt19/test0.pdb
94 : PRINT ARG=rmsd0 FILE=colvar
95 : ```
96 :
97 : Notice that there are sets of weights: $w'$ and $w$ in the formulas above. The first of these weights is used to calculate the position of the center of mass
98 : (so it determines how the atoms are _aligned_). The second set of weights, $w$ is used when calculating how far the atoms have been
99 : _displaced_. These weights are assigned in the reference configuration that you provide as input (i.e. the appear in the input file
100 : to this action that you set using REFERENCE=whatever.pdb). As you can see in the input above, this input consists of a simple pdb file
101 : containing the set of atoms for which you want to calculate the RMSD displacement and their positions in the reference configuration.
102 : It is important to note that the indices in this pdb need to be set correctly. The indices in this file determine the indices of the
103 : instantaneous atomic positions that are used by PLUMED when calculating this colvar. As such if you want to calculate the RMSD distance
104 : moved by the first, fourth, sixth and twenty eighth atoms in the MD codes input file then the indices of the corresponding reference positions in this pdb
105 : file should be set equal to 1, 4, 6 and 28.
106 :
107 : The pdb input file should also contain the weights $w$ and $w'$. These second of these sets of weights, $w'$, appears in the OCCUPANCY column (the first column after the coordinates).
108 : These are the weights that are used to calculate the position of the center of mass. The BETA column (the second column
109 : after the Cartesian coordinates) is used to provide the $w$ values which are used in the calculation of the displacement.
110 : Please note that it is possible to use fractional values for beta and for the occupancy. However, we recommend you only do this when
111 : you really know what you are doing however as the results can be rather strange. A more common practise is to use different sets of atoms
112 : for the alignment and the displacement. You can do this by setting the $w$ and $w'$ values for all the atoms you wish to use for alignment only equal to zero and
113 : one respectively and by setting the $w$ and $w'$ values for all the atoms you wish to use for displacement only to one and zero respectively.
114 :
115 : In the PDB input files that you use for RMSD the atomic coordinates and box lengths should be in Angstroms unless
116 : you are working with natural units. If you are working with natural units then the coordinates
117 : should be in your natural length unit. You can find more details on the PDB file format [here](http://www.wwpdb.org/docs.html).
118 : Please make sure your PDB file is correctly formatted. More detail on the format for PDB files can be found in the documentation for the [PDB2CONSTANT](PDB2CONSTANT.md) action.
119 :
120 : The following input uses a different method to calculate the RMSD distance as you can tell from the TYPE=OPTIMAL on the input line. In addition, because we have added the
121 : SQURED flag on the input line, we are calculating $d(X,X')^2$ rather than $d(X,X')$. Calculating $d(X,X')^2$ is slightly less computationally expensive than computing $d(X,X')$ as
122 : you avoid a square root operation.
123 :
124 : ```plumed
125 : #SETTINGS INPUTFILES=regtest/basic/rt19/test0.pdb
126 : rmsd0: RMSD TYPE=OPTIMAL SQUARED REFERENCE=regtest/basic/rt19/test0.pdb
127 : PRINT ARG=rmsd0 FILE=colvar
128 : ```
129 :
130 : In this case the root mean square deviation is calculated after the positions of geometric centers in the reference and instantaneous configurations are aligned AND after
131 : an optimal alignment of the two frames is performed so that motion due to rotation of the reference frame between the two structures is
132 : removed. The equation for $d(X,X')$ in this case is:
133 :
134 : $$
135 : d(X,X') = \sqrt{ \sum_i \sum_\alpha^{x,y,z} \frac{w_i}{\sum_j w_j}[ X_{i,\alpha}-com_\alpha(X)- \sum_\beta M(X,X',w')_{\alpha,\beta}({X'}_{i,\beta}-com_\beta(X')) ]^2 }
136 : $$
137 :
138 : where $M(X,X',w')$ is the optimal alignment matrix which is calculated using the Kearsley algorithm that is described in the paper from the bibliography below. Again different sets of
139 : weights are used for the alignment ($w'$) and for the displacement calculations ($w$).
140 : This gives a great deal of flexibility as it allows you to use a different sets of atoms (which may or may not overlap) for the alignment and displacement
141 : parts of the calculation. This may be very useful when you want to calculate how a ligand moves about in a protein cavity as you can use the protein as a reference
142 : system and do no alignment of the ligand.
143 :
144 : (Note: when this form of RMSD is used to calculate the secondary structure variables ([ALPHARMSD](ALPHARMSD.md), [ANTIBETARMSD](ANTIBETARMSD.md) and [PARABETARMSD](PARABETARMSD.md)
145 : all the atoms in the segment are assumed to be part of both the alignment and displacement sets and all weights are set equal to one)
146 :
147 : ## Numerical derivatives
148 :
149 : If you are calculating the RMSD distance from a single RMSD reference frame you calculate numerical derivatives of the RMSD distance as shown below for the simple metric:
150 :
151 : ```plumed
152 : #SETTINGS INPUTFILES=regtest/basic/rt19/test0.pdb
153 : rmsd0: RMSD TYPE=SIMPLE REFERENCE=regtest/basic/rt19/test0.pdb
154 : rmsd0n: RMSD TYPE=SIMPLE NUMERICAL_DERIVATIVES REFERENCE=regtest/basic/rt19/test0.pdb
155 : DUMPDERIVATIVES ARG=rmsd0,rmsd0n FILE=deriv
156 : ```
157 :
158 : and as shown below for the OPTIMAL metric:
159 :
160 : ```plumed
161 : #SETTINGS INPUTFILES=regtest/basic/rt19/test0.pdb
162 : rmsd0: RMSD TYPE=OPTIMAL REFERENCE=regtest/basic/rt19/test0.pdb
163 : rmsd0n: RMSD TYPE=OPTIMAL NUMERICAL_DERIVATIVES REFERENCE=regtest/basic/rt19/test0.pdb
164 : DUMPDERIVATIVES ARG=rmsd0,rmsd0n FILE=deriv
165 : ```
166 :
167 : In practice there is no reason to use the NUMERICAL_DERIVATIVES option unless you are testing the RMSD implementation in the way the above two inputs are doing by printing
168 : a file that allows us to compare the analytic and numerical derivatives. Evaluating the derivtives numerically is __much__ more expensive than calculating them analytically.
169 :
170 : ## Computing RMSD displacements
171 :
172 : The $d(X,X')$ values that are calculated when you use the TYPE=SIMPLE and TYPE=OPTIMAL variants of RMSD are scalars. These scalars tell you the length of the vector of displacements,
173 : $X - X'$, between the instantaneous and reference positions. If you would like to access this vector of displacements instead of its length you can by using the DISPLACEMENT keyword
174 : as shown below for TYPE=SIMPLE
175 :
176 : ```plumed
177 : #SETTINGS INPUTFILES=regtest/basic/rt19/test0.pdb
178 : rmsd0: RMSD TYPE=SIMPLE DISPLACEMENT REFERENCE=regtest/basic/rt19/test0.pdb
179 : PRINT ARG=rmsd0.disp,rmsd0.dist FILE=colvar
180 : ```
181 :
182 : or as shown below for TYPE=OPTIMAL
183 :
184 : ```plumed
185 : #SETTINGS INPUTFILES=regtest/basic/rt19/test0.pdb
186 : rmsd0: RMSD TYPE=OPTIMAL DISPLACEMENT REFERENCE=regtest/basic/rt19/test0.pdb
187 : PRINT ARG=rmsd0.disp,rmsd0.dist FILE=colvar
188 : ```
189 :
190 : The RMSD command for these inputs output two components
191 :
192 : - `dist` - the length of displacement vector that is output when you don't use the DISPLACEMENT keyword
193 : - `disp` - the 3N dimensional vector of atomic dispacements, where N is the number of atoms.
194 :
195 : These vectors of displacements are used if you use the [PCAVARS](PCAVARS.md) action to compute the projection of the displacement on a particular reference vector.
196 :
197 : ## Computing multiple RMSD values
198 :
199 : You can also define multiple reference configurations in the reference input as is done in the following example:
200 :
201 : ```plumed
202 : #SETTINGS INPUTFILES=regtest/mapping/rt39/all1.pdb
203 : rmsd: RMSD TYPE=OPTIMAL REFERENCE=regtest/mapping/rt39/all1.pdb
204 : PRINT ARG=rmsd FILE=colvar
205 : ```
206 :
207 : The output from RMSD in this case is a vector that contains the RMSD distances from each of the reference configurations in your input file. This feature is used in the
208 : [PATH](PATH.md) shortcut. Furthermore, you can use the DISPLACEMENT keyword when there are multiple reference configurations in the input file as shown below:
209 :
210 : ```plumed
211 : #SETTINGS INPUTFILES=regtest/mapping/rt39/all1.pdb
212 : rmsd: RMSD TYPE=OPTIMAL DISPLACEMENT REFERENCE=regtest/mapping/rt39/all1.pdb
213 : PRINT ARG=rmsd.disp,rmsd.dist FILE=colvar
214 : ```
215 :
216 : The RMSD command here still outputs the `dist` and `disp` components but now `dist` is a vector and `disp` is a matrix. This type of command is used to implement
217 : the [GEOMETRIC_PATH](GEOMETRIC_PATH.md) shortcut. For this command you need information on the distances from a set of reference configurations that can be found
218 : in the `dist` component as well as the information on the displacement vectors between the instantaneous position and each of the reference configurations that is contained in the `dist` matrix.
219 :
220 : Please note that there are a number of other methods for calculating the distance between the instantaneous configuration and a reference configuration that are available in plumed.
221 :
222 : ## A note on periodic boundary conditions
223 :
224 : When periodic boundary conditions are used, the atoms should be
225 : in the proper periodic image. This has been done automatically since PLUMED 2.5,
226 : by considering the ordered list of atoms and rebuilding molecules using a procedure
227 : that is equivalent to that done in [WHOLEMOLECULES](WHOLEMOLECULES.md). Notice that
228 : rebuilding is local to this action. This is different from [WHOLEMOLECULES](WHOLEMOLECULES.md)
229 : which actually modifies the coordinates stored in PLUMED.
230 :
231 : In case you want to recover the old behavior you should use the NOPBC flag as shown below:
232 :
233 : ```plumed
234 : #SETTINGS INPUTFILES=regtest/basic/rt19/test0.pdb
235 : rmsd0: RMSD TYPE=SIMPLE NOPBC REFERENCE=regtest/basic/rt19/test0.pdb
236 : PRINT ARG=rmsd0 FILE=colvar
237 : ```
238 :
239 : In that case you need to take care that atoms are in the correct
240 : periodic image.
241 :
242 : ## A final thought
243 :
244 : Notice that there are many other ways of calculating the distance from a particular reference structure. These normally work by computing the square root of the sum of the squares
245 : of the differences in collective variable values. To compute these distances you use the functionality in the [refdist](module_refdist.md) module.
246 :
247 : */
248 : //+ENDPLUMEDOC
249 :
250 : PLUMED_REGISTER_ACTION(RMSD,"RMSD_SCALAR")
251 :
252 169 : void RMSD::registerKeywords(Keywords& keys) {
253 169 : Colvar::registerKeywords(keys);
254 338 : keys.setDisplayName("RMSD");
255 169 : keys.add("compulsory","REFERENCE","a file in pdb format containing the reference structure and the atoms involved in the CV.");
256 169 : keys.add("compulsory","TYPE","SIMPLE","the manner in which RMSD alignment is performed. Should be OPTIMAL or SIMPLE.");
257 169 : keys.addFlag("SQUARED",false," This should be set if you want mean squared displacement instead of RMSD ");
258 338 : keys.setValueDescription("scalar","the RMSD between the instantaneous structure and the reference structure that was input");
259 169 : }
260 :
261 83 : RMSD::RMSD(const ActionOptions&ao):
262 : PLUMED_COLVAR_INIT(ao),
263 83 : squared(false),
264 83 : nopbc(false) {
265 : std::string reference;
266 167 : parse("REFERENCE",reference);
267 : std::string type;
268 83 : type.assign("SIMPLE");
269 83 : parse("TYPE",type);
270 83 : parseFlag("SQUARED",squared);
271 83 : parseFlag("NOPBC",nopbc);
272 83 : checkRead();
273 :
274 84 : addValueWithDerivatives();
275 83 : setNotPeriodic();
276 83 : PDB pdb;
277 :
278 : // read everything in ang and transform to nm if we are not in natural units
279 83 : if( !pdb.read(reference,usingNaturalUnits(),0.1/getUnits().getLength()) ) {
280 0 : error("missing input file " + reference );
281 : }
282 83 : myrmsd.set( pdb, type, true, true );
283 :
284 83 : std::vector<AtomNumber> atoms( pdb.getAtomNumbers() );
285 83 : requestAtoms( atoms );
286 82 : der.resize( atoms.size() );
287 :
288 82 : log.printf(" reference from file %s\n",reference.c_str());
289 82 : log.printf(" which contains %d atoms\n",getNumberOfAtoms());
290 82 : log.printf(" with indices : ");
291 3996 : for(unsigned i=0; i<atoms.size(); ++i) {
292 3914 : if(i%25==0) {
293 189 : log<<"\n";
294 : }
295 3914 : log.printf("%d ",atoms[i].serial());
296 : }
297 82 : log.printf("\n");
298 82 : log.printf(" method for alignment : %s \n",type.c_str() );
299 82 : if(squared) {
300 11 : log.printf(" chosen to use SQUARED option for MSD instead of RMSD\n");
301 : }
302 82 : if(nopbc) {
303 61 : log.printf(" without periodic boundary conditions\n");
304 : } else {
305 21 : log.printf(" using periodic boundary conditions\n");
306 : }
307 170 : }
308 :
309 :
310 : // calculator
311 40521 : void RMSD::calculate() {
312 40521 : if(!nopbc) {
313 3959 : makeWhole();
314 : }
315 40521 : double r=myrmsd.calculate( getPositions(), der, squared );
316 :
317 40521 : setValue(r);
318 14618160 : for(unsigned i=0; i<getNumberOfAtoms(); i++) {
319 14577639 : setAtomsDerivatives( i, der[i] );
320 : }
321 40521 : setBoxDerivativesNoPbc();
322 40521 : }
323 :
324 : }
325 : }
326 :
327 :
328 :
|
