Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2018 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 : class RMSD : public Colvar {
37 :
38 : MultiValue myvals;
39 : ReferenceValuePack mypack;
40 : PLMD::RMSDBase* rmsd;
41 : bool squared;
42 :
43 : public:
44 : explicit RMSD(const ActionOptions&);
45 : ~RMSD();
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 1st, 4th, 6th and 28th atoms in the MD codes input file then the indices of the corresponding refernece 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 centre 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 correclty 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 calcuations (\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 : \warning
132 : The molecule used for \ref RMSD calculation should be whole (both atoms used in alignment and in displacement calculation).
133 : In case it is broken by the host MD code, please use \ref WHOLEMOLECULES to reconstruct it before \ref RMSD calculation.
134 :
135 : \par Examples
136 :
137 : The following tells plumed to calculate the RMSD distance between
138 : the positions of the atoms in the reference file and their instantaneous
139 : position. The Kearseley algorithm is used so this is done optimally.
140 :
141 : \verbatim
142 : RMSD REFERENCE=file.pdb TYPE=OPTIMAL
143 : \endverbatim
144 :
145 : ...
146 :
147 : */
148 : //+ENDPLUMEDOC
149 :
150 2592 : PLUMED_REGISTER_ACTION(RMSD,"RMSD")
151 :
152 70 : void RMSD::registerKeywords(Keywords& keys) {
153 70 : Colvar::registerKeywords(keys);
154 70 : keys.add("compulsory","REFERENCE","a file in pdb format containing the reference structure and the atoms involved in the CV.");
155 70 : keys.add("compulsory","TYPE","SIMPLE","the manner in which RMSD alignment is performed. Should be OPTIMAL or SIMPLE.");
156 70 : keys.addFlag("SQUARED",false," This should be setted if you want MSD instead of RMSD ");
157 70 : keys.remove("NOPBC");
158 70 : }
159 :
160 69 : RMSD::RMSD(const ActionOptions&ao):
161 69 : PLUMED_COLVAR_INIT(ao),myvals(1,0), mypack(0,0,myvals),squared(false)
162 : {
163 69 : string reference;
164 69 : parse("REFERENCE",reference);
165 138 : string type;
166 69 : type.assign("SIMPLE");
167 69 : parse("TYPE",type);
168 69 : parseFlag("SQUARED",squared);
169 :
170 69 : checkRead();
171 :
172 :
173 69 : addValueWithDerivatives(); setNotPeriodic();
174 138 : PDB pdb;
175 :
176 : // read everything in ang and transform to nm if we are not in natural units
177 69 : if( !pdb.read(reference,plumed.getAtoms().usingNaturalUnits(),0.1/atoms.getUnits().getLength()) )
178 0 : error("missing input file " + reference );
179 :
180 69 : rmsd = metricRegister().create<RMSDBase>(type,pdb);
181 :
182 138 : std::vector<AtomNumber> atoms;
183 69 : rmsd->getAtomRequests( atoms );
184 : // rmsd->setNumberOfAtoms( atoms.size() );
185 69 : requestAtoms( atoms );
186 :
187 : // Setup the derivative pack
188 69 : myvals.resize( 1, 3*atoms.size()+9 ); mypack.resize( 0, atoms.size() );
189 69 : for(unsigned i=0; i<atoms.size(); ++i) mypack.setAtomIndex( i, i );
190 :
191 69 : log.printf(" reference from file %s\n",reference.c_str());
192 69 : log.printf(" which contains %d atoms\n",getNumberOfAtoms());
193 69 : log.printf(" method for alignment : %s \n",type.c_str() );
194 138 : if(squared)log.printf(" chosen to use SQUARED option for MSD instead of RMSD\n");
195 69 : }
196 :
197 276 : RMSD::~RMSD() {
198 69 : delete rmsd;
199 207 : }
200 :
201 :
202 : // calculator
203 40483 : void RMSD::calculate() {
204 40483 : double r=rmsd->calculate( getPositions(), mypack, squared );
205 :
206 40483 : setValue(r);
207 40483 : for(unsigned i=0; i<getNumberOfAtoms(); i++) setAtomsDerivatives( i, mypack.getAtomDerivative(i) );
208 :
209 40483 : Tensor virial; plumed_dbg_assert( !mypack.virialWasSet() );
210 40483 : setBoxDerivativesNoPbc();
211 40483 : }
212 :
213 : }
214 2523 : }
215 :
216 :
217 :
|
