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