Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-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/MultiDomainRMSD.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 MultiRMSD : public Colvar {
37 :
38 : PLMD::MultiDomainRMSD* rmsd;
39 : bool squared;
40 : MultiValue myvals;
41 : ReferenceValuePack mypack;
42 :
43 : public:
44 : explicit MultiRMSD(const ActionOptions&);
45 : ~MultiRMSD();
46 : virtual void calculate();
47 : static void registerKeywords(Keywords& keys);
48 : };
49 :
50 :
51 : using namespace std;
52 :
53 : //+PLUMEDOC DCOLVAR MULTI-RMSD
54 : /*
55 : Calculate the RMSD distance moved by a number of separated domains from their positions in a reference structure.
56 :
57 :
58 : When you have large proteins the calculation of the root mean squared deviation between all the atoms in a reference
59 : structure and the instantaneous configuration becomes prohibitively expensive. You may thus instead want to calculate
60 : the RMSD between the atoms in a set of domains of your protein and your reference structure. That is to say:
61 :
62 : \f[
63 : d(X,X_r) = \sqrt{ \sum_{i} w_i\vert X_i - X_i' \vert^2 }
64 : \f]
65 :
66 : where here the sum is over the domains of the protein, \f$X_i\f$ represents the positions of the atoms in domain \f$i\f$
67 : in the instantaneous configuration and \f$X_i'\f$ is the positions of the atoms in domain \f$i\f$ in the reference
68 : configuration. \f$w_i\f$ is an optional weight.
69 :
70 : The distances for each of the domains in the above sum can be calculated using the \ref DRMSD or \ref RMSD measures or
71 : using a combination of these distance. The reference configuration is specified in a pdb file like the one below:
72 :
73 : \verbatim
74 : ATOM 2 O ALA 2 -0.926 -2.447 -0.497 1.00 1.00 DIA O
75 : ATOM 4 HNT ALA 2 0.533 -0.396 1.184 1.00 1.00 DIA H
76 : ATOM 6 HT1 ALA 2 -0.216 -2.590 1.371 1.00 1.00 DIA H
77 : ATOM 7 HT2 ALA 2 -0.309 -1.255 2.315 1.00 1.00 DIA H
78 : ATOM 8 HT3 ALA 2 -1.480 -1.560 1.212 1.00 1.00 DIA H
79 : ATOM 9 CAY ALA 2 -0.096 2.144 -0.669 1.00 1.00 DIA C
80 : ATOM 10 HY1 ALA 2 0.871 2.385 -0.588 1.00 1.00 DIA H
81 : TER
82 : ATOM 12 HY3 ALA 2 -0.520 2.679 -1.400 1.00 1.00 DIA H
83 : ATOM 14 OY ALA 2 -1.139 0.931 -0.973 1.00 1.00 DIA O
84 : ATOM 16 HN ALA 2 1.713 1.021 -0.873 1.00 1.00 DIA H
85 : ATOM 18 HA ALA 2 0.099 -0.774 -2.218 1.00 1.00 DIA H
86 : ATOM 19 CB ALA 2 2.063 -1.223 -1.276 1.00 1.00 DIA C
87 : ATOM 20 HB1 ALA 2 2.670 -0.716 -2.057 1.00 1.00 DIA H
88 : ATOM 21 HB2 ALA 2 2.556 -1.051 -0.295 1.00 1.00 DIA H
89 : ATOM 22 HB3 ALA 2 2.070 -2.314 -1.490 1.00 1.00 DIA H
90 : END
91 : \endverbatim
92 :
93 : with the TER keyword being used to separate the various domains in you protein.
94 :
95 :
96 : \par Examples
97 :
98 : The following tells plumed to calculate the RMSD distance between
99 : the positions of the atoms in the reference file and their instantaneous
100 : position. The Kearseley algorithm for each of the domains.
101 :
102 : \verbatim
103 : MULTI-RMSD REFERENCE=file.pdb TYPE=MULTI-OPTIMAL
104 : \endverbatim
105 :
106 : The following tells plumed to calculate the RMSD distance btween the positions of
107 : the atoms in the domains of reference the reference structure and their instantaneous
108 : positions. Here distances are calculated using the \ref DRMSD measure.
109 :
110 : \verbatim
111 : MULTI-RMSD REFERENCE=file.pdb TYPE=MULTI-DRMSD
112 : \endverbatim
113 :
114 : in this case it is possible to use the following DRMSD options in the pdb file using the REMARK syntax:
115 : \verbatim
116 : NOPBC to calculate distances without PBC
117 : LOWER_CUTOFF=# only pairs of atoms further than LOWER_CUTOFF are considered in the calculation
118 : UPPER_CUTOFF=# only pairs of atoms further than UPPER_CUTOFF are considered in the calculation
119 : \endverbatim
120 : as shown in the following example
121 :
122 : \verbatim
123 : REMARK NOPBC
124 : REMARK LOWER_CUTOFF=0.1
125 : REMARK UPPER_CUTOFF=0.8
126 : ATOM 2 O ALA 2 -0.926 -2.447 -0.497 1.00 1.00 DIA O
127 : ATOM 4 HNT ALA 2 0.533 -0.396 1.184 1.00 1.00 DIA H
128 : ATOM 6 HT1 ALA 2 -0.216 -2.590 1.371 1.00 1.00 DIA H
129 : ATOM 7 HT2 ALA 2 -0.309 -1.255 2.315 1.00 1.00 DIA H
130 : ATOM 8 HT3 ALA 2 -1.480 -1.560 1.212 1.00 1.00 DIA H
131 : ATOM 9 CAY ALA 2 -0.096 2.144 -0.669 1.00 1.00 DIA C
132 : ATOM 10 HY1 ALA 2 0.871 2.385 -0.588 1.00 1.00 DIA H
133 : TER
134 : ATOM 12 HY3 ALA 2 -0.520 2.679 -1.400 1.00 1.00 DIA H
135 : ATOM 14 OY ALA 2 -1.139 0.931 -0.973 1.00 1.00 DIA O
136 : ATOM 16 HN ALA 2 1.713 1.021 -0.873 1.00 1.00 DIA H
137 : ATOM 18 HA ALA 2 0.099 -0.774 -2.218 1.00 1.00 DIA H
138 : ATOM 19 CB ALA 2 2.063 -1.223 -1.276 1.00 1.00 DIA C
139 : ATOM 20 HB1 ALA 2 2.670 -0.716 -2.057 1.00 1.00 DIA H
140 : ATOM 21 HB2 ALA 2 2.556 -1.051 -0.295 1.00 1.00 DIA H
141 : ATOM 22 HB3 ALA 2 2.070 -2.314 -1.490 1.00 1.00 DIA H
142 : END
143 : \endverbatim
144 :
145 :
146 : */
147 : //+ENDPLUMEDOC
148 :
149 2526 : PLUMED_REGISTER_ACTION(MultiRMSD,"MULTI-RMSD")
150 :
151 4 : void MultiRMSD::registerKeywords(Keywords& keys) {
152 4 : Colvar::registerKeywords(keys);
153 4 : keys.add("compulsory","REFERENCE","a file in pdb format containing the reference structure and the atoms involved in the CV.");
154 4 : keys.add("compulsory","TYPE","MULTI-SIMPLE","the manner in which RMSD alignment is performed. Should be MULTI-OPTIMAL, MULTI-OPTIMAL-FAST, MULTI-SIMPLE or MULTI-DRMSD.");
155 4 : keys.addFlag("SQUARED",false," This should be setted if you want MSD instead of RMSD ");
156 4 : keys.remove("NOPBC");
157 4 : }
158 :
159 3 : MultiRMSD::MultiRMSD(const ActionOptions&ao):
160 3 : PLUMED_COLVAR_INIT(ao),squared(false),myvals(1,0), mypack(0,0,myvals)
161 : {
162 3 : string reference;
163 3 : parse("REFERENCE",reference);
164 6 : string type;
165 3 : type.assign("SIMPLE");
166 3 : parse("TYPE",type);
167 3 : parseFlag("SQUARED",squared);
168 :
169 3 : checkRead();
170 :
171 :
172 3 : addValueWithDerivatives(); setNotPeriodic();
173 6 : PDB pdb;
174 :
175 : // read everything in ang and transform to nm if we are not in natural units
176 3 : if( !pdb.read(reference,plumed.getAtoms().usingNaturalUnits(),0.1/atoms.getUnits().getLength()) )
177 0 : error("missing input file " + reference );
178 :
179 3 : rmsd = metricRegister().create<MultiDomainRMSD>(type,pdb);
180 :
181 6 : std::vector<AtomNumber> atoms;
182 3 : rmsd->getAtomRequests( atoms );
183 3 : requestAtoms( atoms );
184 :
185 3 : myvals.resize( 1, 3*atoms.size()+9 ); mypack.resize( 0, atoms.size() );
186 3 : for(unsigned i=0; i<atoms.size(); ++i) mypack.setAtomIndex( i, i );
187 :
188 3 : log.printf(" reference from file %s\n",reference.c_str());
189 3 : log.printf(" which contains %d atoms\n",getNumberOfAtoms());
190 3 : log.printf(" method for alignment : %s \n",type.c_str() );
191 6 : if(squared)log.printf(" chosen to use SQUARED option for MSD instead of RMSD\n");
192 3 : }
193 :
194 12 : MultiRMSD::~MultiRMSD() {
195 3 : delete rmsd;
196 9 : }
197 :
198 :
199 : // calculator
200 15 : void MultiRMSD::calculate() {
201 15 : double r=rmsd->calculate( getPositions(), getPbc(), mypack, squared );
202 :
203 15 : setValue(r);
204 15 : for(unsigned i=0; i<getNumberOfAtoms(); i++) setAtomsDerivatives( i, mypack.getAtomDerivative(i) );
205 :
206 15 : if( !mypack.virialWasSet() ) setBoxDerivativesNoPbc();
207 5 : else setBoxDerivatives( mypack.getBoxDerivatives() );
208 15 : }
209 :
210 : }
211 2523 : }
212 :
213 :
214 :
|
