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