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 "ActionRegister.h"
25 : #include "tools/PDB.h"
26 : #include "reference/RMSDBase.h"
27 : #include "reference/MetricRegister.h"
28 : #include "core/Atoms.h"
29 :
30 : namespace PLMD {
31 : namespace colvar {
32 :
33 : class RMSD : public Colvar {
34 :
35 : MultiValue myvals;
36 : ReferenceValuePack mypack;
37 : std::unique_ptr<PLMD::RMSDBase> rmsd;
38 : bool squared;
39 : bool nopbc;
40 :
41 : public:
42 : explicit RMSD(const ActionOptions&);
43 : void calculate() override;
44 : static void registerKeywords(Keywords& keys);
45 : };
46 :
47 : //+PLUMEDOC DCOLVAR RMSD
48 : /*
49 : Calculate the RMSD with respect to a reference structure.
50 :
51 : The aim with this colvar it to calculate something like:
52 :
53 : \f[
54 : d(X,X') = \vert X-X' \vert
55 : \f]
56 :
57 : where \f$ X \f$ is the instantaneous position of all the atoms in the system and
58 : \f$ X' \f$ is the positions of the atoms in some reference structure provided as input.
59 : \f$ d(X,X') \f$ 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 the
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 \f$d(X,x')\f$ is to be calculated using:
68 :
69 : \f[
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 : \f]
72 : with
73 : \f[
74 : com_\alpha(X)= \sum_i \frac{w'_{i}}{\sum_j w'_j}X_{i,\alpha}
75 : \f]
76 : and
77 : \f[
78 : com_\alpha(X')= \sum_i \frac{w'_{i}}{\sum_j w'_j}X'_{i,\alpha}
79 : \f]
80 : Obviously, \f$ com_\alpha(X) \f$ and \f$ com_\alpha(X') \f$ represent the positions of the center of mass in the reference
81 : and instantaneous configurations if the weights $w'$ are set equal to the atomic masses. If the weights are all set equal to
82 : one, however, \f$com_\alpha(X) \f$ and \f$ com_\alpha(X') \f$ are the positions of the geometric centers.
83 : Notice that there are sets of weights: \f$ w' \f$ and \f$ w \f$. The first is used to calculate the position of the center of mass
84 : (so it determines how the atoms are \e aligned). Meanwhile, the second is used when calculating how far the atoms have actually been
85 : \e displaced. These weights are assigned in the reference configuration that you provide as input (i.e. the appear in the input file
86 : to this action that you set using REFERENCE=whatever.pdb). This input reference configuration consists of a simple pdb file
87 : containing the set of atoms for which you want to calculate the RMSD displacement and their positions in the reference configuration.
88 : 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
89 : instantaneous atomic positions that are used by PLUMED when calculating this colvar. As such if you want to calculate the RMSD distance
90 : 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
91 : file should be set equal to 1, 4, 6 and 28.
92 :
93 : The pdb input file should also contain the values of \f$w\f$ and \f$w'\f$. In particular, the OCCUPANCY column (the first column after the coordinates)
94 : is used provides the values of \f$ w'\f$ that are used to calculate the position of the center of mass. The BETA column (the second column
95 : after the Cartesian coordinates) is used to provide the \f$ w \f$ values which are used in the the calculation of the displacement.
96 : Please note that it is possible to use fractional values for beta and for the occupancy. However, we recommend you only do this when
97 : you really know what you are doing however as the results can be rather strange.
98 :
99 : In PDB files the atomic coordinates and box lengths should be in Angstroms unless
100 : you are working with natural units. If you are working with natural units then the coordinates
101 : should be in your natural length unit. For more details on the PDB file format visit http://www.wwpdb.org/docs.html.
102 : Make sure your PDB file is correctly formatted as explained \ref pdbreader "in this page".
103 :
104 : A different method is used to calculate the RMSD distance when you use TYPE=OPTIMAL on the input line. In this case the root mean square
105 : deviation is calculated after the positions of geometric centers in the reference and instantaneous configurations are aligned AND after
106 : an optimal alignment of the two frames is performed so that motion due to rotation of the reference frame between the two structures is
107 : removed. The equation for \f$d(X,X')\f$ in this case reads:
108 :
109 : \f[
110 : 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 }
111 : \f]
112 :
113 : where \f$ M(X,X',w') \f$ is the optimal alignment matrix which is calculated using the Kearsley \cite kearsley algorithm. Again different sets of
114 : weights are used for the alignment (\f$w'\f$) and for the displacement calculations (\f$w\f$).
115 : 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
116 : 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
117 : system and do no alignment of the ligand.
118 :
119 : (Note: when this form of RMSD is used to calculate the secondary structure variables (\ref ALPHARMSD, \ref ANTIBETARMSD and \ref PARABETARMSD
120 : 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)
121 :
122 : Please note that there are a number of other methods for calculating the distance between the instantaneous configuration and a reference configuration
123 : that are available in plumed. More information on these various methods can be found in the section of the manual on \ref dists.
124 :
125 : When running with periodic boundary conditions, the atoms should be
126 : in the proper periodic image. This is done automatically since PLUMED 2.5,
127 : by considering the ordered list of atoms and rebuilding molecules using a procedure
128 : that is equivalent to that done in \ref WHOLEMOLECULES . Notice that
129 : rebuilding is local to this action. This is different from \ref WHOLEMOLECULES
130 : which actually modifies the coordinates stored in PLUMED.
131 :
132 : In case you want to recover the old behavior you should use the NOPBC flag.
133 : In that case you need to take care that atoms are in the correct
134 : periodic image.
135 :
136 : \par Examples
137 :
138 : The following tells plumed to calculate the RMSD distance between
139 : the positions of the atoms in the reference file and their instantaneous
140 : position. The Kearsley algorithm is used so this is done optimally.
141 :
142 : \plumedfile
143 : RMSD REFERENCE=file.pdb TYPE=OPTIMAL
144 : \endplumedfile
145 :
146 : The reference configuration is specified in a pdb file that will have a format similar to the one shown below:
147 :
148 : \auxfile{file.pdb}
149 : ATOM 1 CL ALA 1 -3.171 0.295 2.045 1.00 1.00
150 : ATOM 5 CLP ALA 1 -1.819 -0.143 1.679 1.00 1.00
151 : ATOM 6 OL ALA 1 -1.177 -0.889 2.401 1.00 1.00
152 : ATOM 7 NL ALA 1 -1.313 0.341 0.529 1.00 1.00
153 : ATOM 8 HL ALA 1 -1.845 0.961 -0.011 1.00 1.00
154 : END
155 : \endauxfile
156 :
157 : ...
158 :
159 : */
160 : //+ENDPLUMEDOC
161 :
162 10561 : PLUMED_REGISTER_ACTION(RMSD,"RMSD")
163 :
164 73 : void RMSD::registerKeywords(Keywords& keys) {
165 73 : Colvar::registerKeywords(keys);
166 146 : keys.add("compulsory","REFERENCE","a file in pdb format containing the reference structure and the atoms involved in the CV.");
167 146 : keys.add("compulsory","TYPE","SIMPLE","the manner in which RMSD alignment is performed. Should be OPTIMAL or SIMPLE.");
168 146 : keys.addFlag("SQUARED",false," This should be set if you want mean squared displacement instead of RMSD ");
169 73 : }
170 :
171 72 : RMSD::RMSD(const ActionOptions&ao):
172 : PLUMED_COLVAR_INIT(ao),
173 146 : myvals(1,0),
174 72 : mypack(0,0,myvals),
175 72 : squared(false),
176 144 : nopbc(false)
177 : {
178 : std::string reference;
179 146 : parse("REFERENCE",reference);
180 : std::string type;
181 72 : type.assign("SIMPLE");
182 72 : parse("TYPE",type);
183 72 : parseFlag("SQUARED",squared);
184 74 : parseFlag("NOPBC",nopbc);
185 :
186 72 : checkRead();
187 :
188 :
189 72 : addValueWithDerivatives(); setNotPeriodic();
190 72 : PDB pdb;
191 :
192 : // read everything in ang and transform to nm if we are not in natural units
193 144 : if( !pdb.read(reference,plumed.getAtoms().usingNaturalUnits(),0.1/atoms.getUnits().getLength()) )
194 1 : error("missing input file " + reference );
195 :
196 142 : rmsd=metricRegister().create<RMSDBase>(type,pdb);
197 :
198 : std::vector<AtomNumber> atoms;
199 71 : rmsd->getAtomRequests( atoms );
200 : // rmsd->setNumberOfAtoms( atoms.size() );
201 71 : requestAtoms( atoms );
202 :
203 : // Setup the derivative pack
204 71 : myvals.resize( 1, 3*atoms.size()+9 ); mypack.resize( 0, atoms.size() );
205 3089 : for(unsigned i=0; i<atoms.size(); ++i) mypack.setAtomIndex( i, i );
206 :
207 70 : log.printf(" reference from file %s\n",reference.c_str());
208 70 : log.printf(" which contains %d atoms\n",getNumberOfAtoms());
209 70 : log.printf(" with indices : ");
210 3089 : for(unsigned i=0; i<atoms.size(); ++i) {
211 3019 : if(i%25==0) log<<"\n";
212 3019 : log.printf("%d ",atoms[i].serial());
213 : }
214 70 : log.printf("\n");
215 70 : log.printf(" method for alignment : %s \n",type.c_str() );
216 70 : if(squared)log.printf(" chosen to use SQUARED option for MSD instead of RMSD\n");
217 70 : if(nopbc) log.printf(" without periodic boundary conditions\n");
218 16 : else log.printf(" using periodic boundary conditions\n");
219 152 : }
220 :
221 :
222 : // calculator
223 40484 : void RMSD::calculate() {
224 40484 : if(!nopbc) makeWhole();
225 40484 : double r=rmsd->calculate( getPositions(), mypack, squared );
226 :
227 40484 : setValue(r);
228 14616239 : for(unsigned i=0; i<getNumberOfAtoms(); i++) setAtomsDerivatives( i, mypack.getAtomDerivative(i) );
229 :
230 40484 : Tensor virial; plumed_dbg_assert( !mypack.virialWasSet() );
231 40484 : setBoxDerivativesNoPbc();
232 40484 : }
233 :
234 : }
235 : }
236 :
237 :
238 :
|
