Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-2024 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 "core/ActionShortcut.h"
23 : #include "core/ActionRegister.h"
24 : #include "core/PlumedMain.h"
25 : #include "core/ActionSet.h"
26 : #include "core/Group.h"
27 : #include "tools/PDB.h"
28 :
29 : namespace PLMD {
30 : namespace colvar {
31 :
32 : //+PLUMEDOC DCOLVAR DRMSD
33 : /*
34 : Calculate the distance RMSD with respect to a reference structure.
35 :
36 : To calculate the root-mean-square deviation between the atoms in two configurations
37 : you must first superimpose the two structures in some ways. Obviously, it is the internal vibrational
38 : motions of the structure - i.e. not the translations and rotations - that are interesting. However,
39 : aligning two structures by removing the translational and rotational motions is not easy. Furthermore,
40 : in some cases there can be alignment issues caused by so-called frame-fitting problems. It is thus
41 : often cheaper and easier to calculate the distances between all the pairs of atoms. The distance
42 : between the two structures, \f$\mathbf{X}^a\f$ and \f$\mathbf{X}^b\f$ can then be measured as:
43 :
44 : \f[
45 : d(\mathbf{X}^A, \mathbf{X}^B) = \sqrt{\frac{1}{N(N-1)} \sum_{i \ne j} [ d(\mathbf{x}_i^a,\mathbf{x}_j^a) - d(\mathbf{x}_i^b,\mathbf{x}_j^b) ]^2}
46 : \f]
47 :
48 : where \f$N\f$ is the number of atoms and \f$d(\mathbf{x}_i,\mathbf{x}_j)\f$ represents the distance between
49 : atoms \f$i\f$ and \f$j\f$. Clearly, this representation of the configuration is invariant to translation and rotation.
50 : However, it can become expensive to calculate when the number of atoms is large. This can be resolved
51 : within the DRMSD colvar by setting LOWER_CUTOFF and UPPER_CUTOFF. These keywords ensure that only
52 : pairs of atoms that are within a certain range are incorporated into the above sum.
53 :
54 : In PDB files the atomic coordinates and box lengths should be in Angstroms unless
55 : you are working with natural units. If you are working with natural units then the coordinates
56 : should be in your natural length unit. For more details on the PDB file format visit http://www.wwpdb.org/docs.html
57 :
58 : \par Examples
59 :
60 : The following tells plumed to calculate the distance RMSD between
61 : the positions of the atoms in the reference file and their instantaneous
62 : position. Only pairs of atoms whose distance in the reference structure is within
63 : 0.1 and 0.8 nm are considered.
64 :
65 : \plumedfile
66 : DRMSD REFERENCE=file1.pdb LOWER_CUTOFF=0.1 UPPER_CUTOFF=0.8
67 : \endplumedfile
68 :
69 : The reference file is a PDB file that looks like this
70 :
71 : \auxfile{file1.pdb}
72 : ATOM 8 HT3 ALA 2 -1.480 -1.560 1.212 1.00 1.00 DIA H
73 : ATOM 9 CAY ALA 2 -0.096 2.144 -0.669 1.00 1.00 DIA C
74 : ATOM 10 HY1 ALA 2 0.871 2.385 -0.588 1.00 1.00 DIA H
75 : ATOM 12 HY3 ALA 2 -0.520 2.679 -1.400 1.00 1.00 DIA H
76 : ATOM 14 OY ALA 2 -1.139 0.931 -0.973 1.00 1.00 DIA O
77 : END
78 : \endauxfile
79 :
80 : The following tells plumed to calculate a DRMSD value for a pair of molecules.
81 :
82 : \plumedfile
83 : DRMSD REFERENCE=file2.pdb LOWER_CUTOFF=0.1 UPPER_CUTOFF=0.8 TYPE=INTER-DRMSD
84 : \endplumedfile
85 :
86 : In the input reference file (file.pdb) the atoms in each of the two molecules are separated by a TER
87 : command as shown below.
88 :
89 : \auxfile{file2.pdb}
90 : ATOM 8 HT3 ALA 2 -1.480 -1.560 1.212 1.00 1.00 DIA H
91 : ATOM 9 CAY ALA 2 -0.096 2.144 -0.669 1.00 1.00 DIA C
92 : ATOM 10 HY1 ALA 2 0.871 2.385 -0.588 1.00 1.00 DIA H
93 : TER
94 : ATOM 12 HY3 ALA 2 -0.520 2.679 -1.400 1.00 1.00 DIA H
95 : ATOM 14 OY ALA 2 -1.139 0.931 -0.973 1.00 1.00 DIA O
96 : END
97 : \endauxfile
98 :
99 : In this example the INTER-DRMSD type ensures that the set of distances from which the final
100 : quantity is computed involve one atom from each of the two molecules. If this is replaced
101 : by INTRA-DRMSD then only those distances involving pairs of atoms that are both in the same
102 : molecule are computed.
103 :
104 : */
105 : //+ENDPLUMEDOC
106 :
107 : class DRMSD : public ActionShortcut {
108 : public:
109 : static void registerKeywords(Keywords& keys);
110 : explicit DRMSD(const ActionOptions&);
111 : };
112 :
113 : PLUMED_REGISTER_ACTION(DRMSD,"DRMSD")
114 :
115 16 : void DRMSD::registerKeywords( Keywords& keys ) {
116 16 : ActionShortcut::registerKeywords( keys );
117 32 : keys.add("compulsory","REFERENCE","a file in pdb format containing the reference structure and the atoms involved in the CV.");
118 32 : keys.add("optional","LOWER_CUTOFF","only pairs of atoms further than LOWER_CUTOFF are considered in the calculation.");
119 32 : keys.add("optional","UPPER_CUTOFF","only pairs of atoms closer than UPPER_CUTOFF are considered in the calculation.");
120 32 : keys.add("compulsory","TYPE","DRMSD","what kind of DRMSD would you like to calculate. You can use either the normal DRMSD involving all the distances between "
121 : "the atoms in your molecule. Alternatively, if you have multiple molecules you can use the type INTER-DRMSD "
122 : "to compute DRMSD values involving only those distances between the atoms at least two molecules or the type INTRA-DRMSD "
123 : "to compute DRMSD values involving only those distances between atoms in the same molecule");
124 32 : keys.addFlag("SQUARED",false,"This should be setted if you want MSD instead of RMSD ");
125 32 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
126 : // This is just ignored in reality which is probably bad
127 32 : keys.addFlag("NUMERICAL_DERIVATIVES",false,"calculate the derivatives for these quantities numerically");
128 32 : keys.setValueDescription("scalar/vector","the DRMSD distance between the instantaneous structure and the reference structure");
129 80 : keys.needsAction("SUM"); keys.needsAction("DISTANCE"); keys.needsAction("CONSTANT"); keys.needsAction("EUCLIDEAN_DISTANCE"); keys.needsAction("CUSTOM");
130 16 : }
131 :
132 13 : DRMSD::DRMSD( const ActionOptions& ao ):
133 : Action(ao),
134 13 : ActionShortcut(ao)
135 : {
136 : // Read in the reference configuration
137 13 : std::string reference; parse("REFERENCE",reference);
138 : // First bit of input for the instantaneous distances
139 26 : bool numder; parseFlag("NUMERICAL_DERIVATIVES",numder); double fake_unit=0.1;
140 13 : FILE* fp2=fopen(reference.c_str(),"r"); bool do_read=true; unsigned nframes=0;
141 65 : while( do_read ) {
142 54 : PDB mypdb; do_read=mypdb.readFromFilepointer(fp2,false,fake_unit);
143 54 : if( !do_read && nframes>0 ) break ;
144 53 : nframes++;
145 54 : }
146 12 : fclose(fp2);
147 :
148 : // Get cutoff information
149 25 : double lcut=0; parse("LOWER_CUTOFF",lcut); std::string drmsd_type; parse("TYPE",drmsd_type);
150 12 : double ucut=std::numeric_limits<double>::max(); parse("UPPER_CUTOFF",ucut);
151 24 : bool nopbc; parseFlag("NOPBC",nopbc); std::string pbc_str; if(nopbc) pbc_str="NOPBC";
152 : // Open the pdb file
153 12 : FILE* fp=fopen(reference.c_str(),"r"); do_read=true;
154 12 : if(!fp) error("could not open reference file " + reference ); unsigned n=0; std::string allpairs="";
155 : std::vector<std::pair<unsigned,unsigned> > upairs; std::vector<std::string> refvals;
156 65 : while ( do_read ) {
157 54 : PDB mypdb; do_read=mypdb.readFromFilepointer(fp,false,fake_unit);
158 54 : if( !do_read && n>0 ) break ;
159 53 : std::vector<Vector> pos( mypdb.getPositions() ); unsigned nn=1;
160 53 : if( pos.size()==0 ) error("read no atoms from file named " + reference );
161 : // This is what we do for the first frame
162 53 : if( n==0 ) {
163 12 : std::vector<AtomNumber> atoms( mypdb.getAtomNumbers() );
164 12 : if( drmsd_type=="DRMSD" ) {
165 102 : for(unsigned i=0; i<atoms.size()-1; ++i) {
166 92 : std::string istr; Tools::convert( atoms[i].serial(), istr );
167 671 : for(unsigned j=i+1; j<atoms.size(); ++j) {
168 579 : std::string jstr; Tools::convert( atoms[j].serial(), jstr );
169 579 : double distance = delta( pos[i], pos[j] ).modulo();
170 579 : if( distance < ucut && distance > lcut ) {
171 386 : std::string num; Tools::convert( nn, num ); nn++;
172 : // Add this pair to list of pairs
173 386 : upairs.push_back( std::pair<unsigned,unsigned>(i,j) );
174 : // Add this distance to list of reference values
175 386 : std::string dstr; Tools::convert( distance, dstr ); refvals.push_back( dstr );
176 : // Calculate this distance
177 694 : if( nframes==1 ) allpairs += " ATOMS" + num + "=" + istr + "," + jstr;
178 156 : else readInputLine( getShortcutLabel() + "_d" + num + ": DISTANCE ATOMS=" + istr + "," + jstr + " " + pbc_str );
179 : }
180 : }
181 : }
182 : } else {
183 2 : unsigned nblocks = mypdb.getNumberOfAtomBlocks(); std::vector<unsigned> blocks( 1 + nblocks );
184 2 : if( nblocks==1 ) { blocks[0]=0; blocks[1]=atoms.size(); }
185 6 : else { blocks[0]=0; for(unsigned i=0; i<nblocks; ++i) blocks[i+1]=mypdb.getAtomBlockEnds()[i]; }
186 2 : if( drmsd_type=="INTRA-DRMSD" ) {
187 3 : for(unsigned i=0; i<nblocks; ++i) {
188 7 : for(unsigned iatom=blocks[i]+1; iatom<blocks[i+1]; ++iatom) {
189 5 : std::string istr; Tools::convert( atoms[iatom].serial(), istr );
190 14 : for(unsigned jatom=blocks[i]; jatom<iatom; ++jatom) {
191 9 : std::string jstr; Tools::convert( atoms[jatom].serial(), jstr );
192 9 : double distance = delta( pos[iatom], pos[jatom] ).modulo();
193 9 : if(distance < ucut && distance > lcut ) {
194 9 : std::string num; Tools::convert( nn, num ); nn++;
195 : // Add this pair to list of pairs
196 9 : upairs.push_back( std::pair<unsigned,unsigned>(iatom,jatom) );
197 : // Add this distance to list of reference values
198 9 : std::string dstr; Tools::convert( distance, dstr ); refvals.push_back( dstr );
199 : // Calculate this distance
200 18 : if( nframes==1 ) allpairs += " ATOMS" + num + "=" + istr + "," + jstr;
201 0 : else readInputLine( getShortcutLabel() + "_d" + num + ": DISTANCE ATOMS=" + istr + "," + jstr + " " + pbc_str );
202 : }
203 : }
204 : }
205 : }
206 1 : } else if( drmsd_type=="INTER-DRMSD" ) {
207 2 : for(unsigned i=1; i<nblocks; ++i) {
208 2 : for(unsigned j=0; j<i; ++j) {
209 4 : for(unsigned iatom=blocks[i]; iatom<blocks[i+1]; ++iatom) {
210 3 : std::string istr; Tools::convert( atoms[iatom].serial(), istr );
211 15 : for(unsigned jatom=blocks[j]; jatom<blocks[j+1]; ++jatom) {
212 12 : std::string jstr; Tools::convert( atoms[jatom].serial(), jstr );
213 12 : double distance = delta( pos[iatom], pos[jatom] ).modulo();
214 12 : if(distance < ucut && distance > lcut ) {
215 12 : std::string num; Tools::convert( nn, num ); nn++;
216 : // Add this pair to list of pairs
217 12 : upairs.push_back( std::pair<unsigned,unsigned>(iatom,jatom) );
218 : // Add this distance to list of reference values
219 12 : std::string dstr; Tools::convert( distance, dstr ); refvals.push_back( dstr );
220 : // Calculate this distance
221 24 : if( nframes==1 ) allpairs += " ATOMS" + num + "=" + istr + "," + jstr;
222 0 : else readInputLine( getShortcutLabel() + "_d" + num + ": DISTANCE ATOMS=" + istr + "," + jstr + " " + pbc_str );
223 : }
224 : }
225 : }
226 : }
227 : }
228 0 : } else plumed_merror( drmsd_type + " is not valid input to TYPE keyword");
229 : }
230 : // This is for every subsequent frame
231 : } else {
232 3239 : for(unsigned i=0; i<refvals.size(); ++i) {
233 3198 : std::string dstr; Tools::convert( delta( pos[upairs[i].first], pos[upairs[i].second] ).modulo(), dstr );
234 6396 : refvals[i] += "," + dstr;
235 : }
236 : }
237 53 : n++;
238 54 : }
239 : // Now create values that hold all the reference distances
240 12 : fclose(fp);
241 :
242 12 : if( nframes==1 ) {
243 22 : readInputLine( getShortcutLabel() + "_d: DISTANCE" + allpairs + " " + pbc_str );
244 340 : std::string refstr = refvals[0]; for(unsigned i=1; i<refvals.size(); ++i) refstr += "," + refvals[i];
245 22 : readInputLine( getShortcutLabel() + "_ref: CONSTANT VALUES=" + refstr );
246 22 : readInputLine( getShortcutLabel() + "_diffs: CUSTOM ARG=" + getShortcutLabel() + "_d," + getShortcutLabel() + "_ref FUNC=(x-y)*(x-y) PERIODIC=NO");
247 22 : readInputLine( getShortcutLabel() + "_u: SUM ARG=" + getShortcutLabel() + "_diffs PERIODIC=NO");
248 : } else {
249 : std::string arg_str1, arg_str2;
250 79 : for(unsigned i=0; i<refvals.size(); ++i ) {
251 78 : std::string inum; Tools::convert( i+1, inum );
252 156 : readInputLine( getShortcutLabel() + "_ref" + inum + ": CONSTANT VALUES=" + refvals[i] );
253 80 : if( i==0 ) { arg_str1 = getShortcutLabel() + "_d" + inum; arg_str2 = getShortcutLabel() + "_ref" + inum; }
254 231 : else { arg_str1 += "," + getShortcutLabel() + "_d" + inum; arg_str2 += "," + getShortcutLabel() + "_ref" + inum; }
255 : }
256 : // And calculate the euclidean distances between the true distances and the references
257 2 : readInputLine( getShortcutLabel() + "_u: EUCLIDEAN_DISTANCE SQUARED ARG1=" + arg_str1 + " ARG2=" + arg_str2 );
258 : }
259 : // And final value
260 12 : std::string nvals; Tools::convert( refvals.size(), nvals ); bool squared; parseFlag("SQUARED",squared);
261 15 : if( squared ) readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_u FUNC=x/" + nvals + " PERIODIC=NO");
262 : else {
263 18 : readInputLine( getShortcutLabel() + "_2: CUSTOM ARG=" + getShortcutLabel() + "_u FUNC=(x/" + nvals + ") PERIODIC=NO");
264 18 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_2 FUNC=sqrt(x) PERIODIC=NO");
265 : }
266 26 : }
267 :
268 : }
269 : }
270 :
271 :
|
