Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-2019 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/DRMSD.h"
27 : #include "reference/MetricRegister.h"
28 : #include "core/Atoms.h"
29 :
30 : using namespace std;
31 :
32 : namespace PLMD {
33 : namespace colvar {
34 :
35 : //+PLUMEDOC DCOLVAR DRMSD
36 : /*
37 : Calculate the distance RMSD with respect to a reference structure.
38 :
39 : To calculate the root-mean-square deviation between the atoms in two configurations
40 : you must first superimpose the two structures in some ways. Obviously, it is the internal vibrational
41 : motions of the structure - i.e. not the translations and rotations - that are interesting. However,
42 : aligning two structures by removing the translational and rotational motions is not easy. Furthermore,
43 : in some cases there can be alignment issues caused by so-called frame-fitting problems. It is thus
44 : often cheaper and easier to calculate the distances between all the pairs of atoms. The distance
45 : between the two structures, \f$\mathbf{X}^a\f$ and \f$\mathbf{X}^b\f$ can then be measured as:
46 :
47 : \f[
48 : 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}
49 : \f]
50 :
51 : where \f$N\f$ is the number of atoms and \f$d(\mathbf{x}_i,\mathbf{x}_j)\f$ represents the distance between
52 : atoms \f$i\f$ and \f$j\f$. Clearly, this representation of the configuration is invariant to translation and rotation.
53 : However, it can become expensive to calculate when the number of atoms is large. This can be resolved
54 : within the DRMSD colvar by setting LOWER_CUTOFF and UPPER_CUTOFF. These keywords ensure that only
55 : pairs of atoms that are within a certain range are incorporated into the above sum.
56 :
57 : In PDB files the atomic coordinates and box lengths should be in Angstroms unless
58 : you are working with natural units. If you are working with natural units then the coordinates
59 : should be in your natural length unit. For more details on the PDB file format visit http://www.wwpdb.org/docs.html
60 :
61 : \par Examples
62 :
63 : The following tells plumed to calculate the distance RMSD between
64 : the positions of the atoms in the reference file and their instantaneous
65 : position. Only pairs of atoms whose distance in the reference structure is within
66 : 0.1 and 0.8 nm are considered.
67 :
68 : \plumedfile
69 : DRMSD REFERENCE=file.pdb LOWER_CUTOFF=0.1 UPPER_CUTOFF=0.8
70 : \endplumedfile
71 :
72 : The following tells plumed to calculate a DRMSD value for a pair of molecules.
73 :
74 : \plumedfile
75 : DRMSD REFERENCE=file.pdb LOWER_CUTOFF=0.1 UPPER_CUTOFF=0.8 TYPE=INTER-DRMSD
76 : \endplumedfile
77 :
78 : In the input reference file (file.pdb) the atoms in each of the two molecules are separated by a TER
79 : command as shown below.
80 :
81 : \verbatim
82 : ATOM 8 HT3 ALA 2 -1.480 -1.560 1.212 1.00 1.00 DIA H
83 : ATOM 9 CAY ALA 2 -0.096 2.144 -0.669 1.00 1.00 DIA C
84 : ATOM 10 HY1 ALA 2 0.871 2.385 -0.588 1.00 1.00 DIA H
85 : TER
86 : ATOM 12 HY3 ALA 2 -0.520 2.679 -1.400 1.00 1.00 DIA H
87 : ATOM 14 OY ALA 2 -1.139 0.931 -0.973 1.00 1.00 DIA O
88 : END
89 : \endverbatim
90 :
91 : In this example the INTER-DRMSD type ensures that the set of distances from which the final
92 : quantity is computed involve one atom from each of the two molecules. If this is replaced
93 : by INTRA-DRMSD then only those distances involving pairs of atoms that are both in the same
94 : molecule are computed.
95 :
96 : */
97 : //+ENDPLUMEDOC
98 :
99 :
100 : class DRMSD : public Colvar {
101 :
102 : bool pbc_;
103 : MultiValue myvals;
104 : ReferenceValuePack mypack;
105 : PLMD::DRMSD* drmsd_;
106 :
107 : public:
108 : explicit DRMSD(const ActionOptions&);
109 : ~DRMSD();
110 : virtual void calculate();
111 : static void registerKeywords(Keywords& keys);
112 : };
113 :
114 6464 : PLUMED_REGISTER_ACTION(DRMSD,"DRMSD")
115 :
116 13 : void DRMSD::registerKeywords(Keywords& keys) {
117 13 : Colvar::registerKeywords(keys);
118 52 : keys.add("compulsory","REFERENCE","a file in pdb format containing the reference structure and the atoms involved in the CV.");
119 52 : keys.add("compulsory","LOWER_CUTOFF","only pairs of atoms further than LOWER_CUTOFF are considered in the calculation.");
120 52 : keys.add("compulsory","UPPER_CUTOFF","only pairs of atoms closer than UPPER_CUTOFF are considered in the calculation.");
121 65 : 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 "
122 : "the atoms in your molecule. Alternatively, if you have multiple molecules you can use the type INTER-DRMSD "
123 : "to compute DRMSD values involving only those distances between the atoms at least two molecules or the type INTRA-DRMSD "
124 : "to compute DRMSD values involving only those distances between atoms in the same molecule");
125 13 : }
126 :
127 12 : DRMSD::DRMSD(const ActionOptions&ao):
128 13 : PLUMED_COLVAR_INIT(ao), pbc_(true), myvals(1,0), mypack(0,0,myvals)
129 : {
130 : string reference;
131 24 : parse("REFERENCE",reference);
132 : double lcutoff;
133 24 : parse("LOWER_CUTOFF",lcutoff);
134 : double ucutoff;
135 24 : parse("UPPER_CUTOFF",ucutoff);
136 12 : bool nopbc(false);
137 24 : parseFlag("NOPBC",nopbc);
138 12 : pbc_=!nopbc;
139 :
140 12 : addValueWithDerivatives(); setNotPeriodic();
141 :
142 : // read everything in ang and transform to nm if we are not in natural units
143 24 : PDB pdb;
144 24 : if( !pdb.read(reference,plumed.getAtoms().usingNaturalUnits(),0.1/atoms.getUnits().getLength()) )
145 2 : error("missing input file " + reference );
146 :
147 : // store target_ distance
148 22 : std::string type; parse("TYPE",type);
149 11 : drmsd_= metricRegister().create<PLMD::DRMSD>( type );
150 11 : drmsd_->setBoundsOnDistances( !nopbc, lcutoff, ucutoff );
151 11 : drmsd_->set( pdb );
152 11 : checkRead();
153 :
154 : std::vector<AtomNumber> atoms;
155 11 : drmsd_->getAtomRequests( atoms );
156 : // drmsd_->setNumberOfAtoms( atoms.size() );
157 11 : requestAtoms( atoms );
158 :
159 : // Setup the derivative pack
160 22 : myvals.resize( 1, 3*atoms.size()+9 ); mypack.resize( 0, atoms.size() );
161 376 : for(unsigned i=0; i<atoms.size(); ++i) mypack.setAtomIndex( i, i );
162 :
163 22 : log.printf(" reference from file %s\n",reference.c_str());
164 22 : log.printf(" which contains %d atoms\n",getNumberOfAtoms());
165 :
166 11 : }
167 :
168 33 : DRMSD::~DRMSD() {
169 11 : delete drmsd_;
170 22 : }
171 :
172 : // calculator
173 595 : void DRMSD::calculate() {
174 :
175 595 : double drmsd; Tensor virial; mypack.clear();
176 1190 : drmsd=drmsd_->calculate(getPositions(), getPbc(), mypack, false);
177 :
178 595 : setValue(drmsd);
179 27185 : for(unsigned i=0; i<getNumberOfAtoms(); ++i) { if( myvals.isActive(3*i) ) setAtomsDerivatives( i, mypack.getAtomDerivative(i) ); }
180 1190 : setBoxDerivatives( mypack.getBoxDerivatives() );
181 595 : }
182 :
183 : }
184 4839 : }
|
