Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2014-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 "core/ActionAtomistic.h"
23 : #include "core/ActionPilot.h"
24 : #include "core/ActionRegister.h"
25 : #include "core/ActionWithValue.h"
26 : #include "tools/Vector.h"
27 : #include "tools/Matrix.h"
28 : #include "tools/AtomNumber.h"
29 : #include "tools/Tools.h"
30 : #include "tools/RMSD.h"
31 : #include "core/Atoms.h"
32 : #include "core/PlumedMain.h"
33 : #include "core/ActionSet.h"
34 : #include "core/SetupMolInfo.h"
35 : #include "tools/PDB.h"
36 : #include "tools/Pbc.h"
37 :
38 : #include <vector>
39 : #include <string>
40 :
41 : using namespace std;
42 :
43 : namespace PLMD {
44 : namespace generic {
45 :
46 : //+PLUMEDOC GENERIC FIT_TO_TEMPLATE
47 : /*
48 : This action is used to align a molecule to a template.
49 :
50 : This can be used to move the coordinates stored in plumed
51 : so as to be aligned with a provided template in PDB format. Pdb should contain
52 : also weights for alignment (see the format of PDB files used e.g. for \ref RMSD).
53 : Make sure your PDB file is correclty formatted as explained \ref pdbreader "in this page".
54 : Weights for displacement are ignored, since no displacement is computed here.
55 : Notice that all atoms (not only those in the template) are aligned.
56 : To see what effect try
57 : the \ref DUMPATOMS directive to output the atomic positions.
58 :
59 : Also notice that PLUMED propagate forces correctly so that you can add a bias on a CV computed
60 : after alignment. For many CVs this has no effect, but in some case the alignment can
61 : change the result. Examples are:
62 : - \ref POSITION CV since it is affected by a rigid shift of the system.
63 : - \ref DISTANCE CV with COMPONENTS. Since the alignment could involve a rotation (with TYPE=OPTIMAL) the actual components could be different
64 : from the original ones.
65 : - \ref CELL components for a similar reason.
66 : - \ref DISTANCE from a \ref FIXEDATOM, provided the fixed atom is introduced _after_ the \ref FIT_TO_TEMPLATE action.
67 :
68 : \attention
69 : The implementation of TYPE=OPTIMAL is available but should be considered in testing phase. Please report any
70 : strange behavior.
71 :
72 : \attention
73 : This directive modifies the stored position at the precise moment
74 : it is executed. This means that only collective variables
75 : which are below it in the input script will see the corrected positions.
76 : As a general rule, put it at the top of the input file. Also, unless you
77 : know exactly what you are doing, leave the default stride (1), so that
78 : this action is performed at every MD step.
79 :
80 : \warning
81 : The molecule used for alignment should be whole.
82 : In case it is broken by the host MD code, please use \ref WHOLEMOLECULES to reconstruct it before \ref FIT_TO_TEMPLATE .
83 :
84 :
85 : \par Examples
86 :
87 : Align the atomic position to a template then print them.
88 : The following example is only translating the system so as
89 : to align the center of mass of a molecule to the one in the reference
90 : structure `ref.pdb`:
91 : \plumedfile
92 : # dump coordinates before fitting, to see the difference:
93 : DUMPATOMS FILE=dump-before.xyz ATOMS=1-20
94 :
95 : # fit coordinates to ref.pdb template
96 : # this is a "TYPE=SIMPLE" fit, so that only translations are used.
97 : FIT_TO_TEMPLATE STRIDE=1 REFERENCE=ref.pdb TYPE=SIMPLE
98 :
99 : # dump coordinates after fitting, to see the difference:
100 : DUMPATOMS FILE=dump-after.xyz ATOMS=1-20
101 : \endplumedfile
102 :
103 : The following example instead performs a rototranslational fit.
104 : \plumedfile
105 : # dump coordinates before fitting, to see the difference:
106 : DUMPATOMS FILE=dump-before.xyz ATOMS=1-20
107 :
108 : # fit coordinates to ref.pdb template
109 : # this is a "TYPE=OPTIMAL" fit, so that rototranslations are used.
110 : FIT_TO_TEMPLATE STRIDE=1 REFERENCE=ref.pdb TYPE=OPTIMAL
111 :
112 : # dump coordinates after fitting, to see the difference:
113 : DUMPATOMS FILE=dump-after.xyz ATOMS=1-20
114 : \endplumedfile
115 :
116 : In the following example you see two completely equivalent way
117 : to restrain an atom close to a position that is defined in the reference
118 : frame of an aligned molecule. It could be for instance the center of mass
119 : of a ligand with respect to a protein
120 : \plumedfile
121 : # center of the ligand:
122 : ce: CENTER ATOMS=100-110
123 :
124 : FIT_TO_TEMPLATE REFERENCE=protein.pdb TYPE=OPTIMAL
125 :
126 : # place a fixed atom in the protein reference coordinates:
127 : fix: FIXEDATOM AT=1.0,1.1,1.0
128 :
129 : # take the distance between the fixed atom and the center of the ligand
130 : d: DISTANCE ATOMS=ce,fix
131 :
132 : # apply a restraint
133 : RESTRAINT ARG=d AT=0.0 KAPPA=100.0
134 : \endplumedfile
135 :
136 : Notice that you could have obtained an (almost) identical result adding a fictitious
137 : atom to `ref.pdb` with the serial number corresponding to the `ce` atom (there is no automatic way
138 : to get it, but in this example it should be the number of atoms of the system plus one),
139 : and properly setting the weights for alignment and displacement in \ref RMSD.
140 : There are two differences to be expected:
141 : (ab) \ref FIT_TO_TEMPLATE might be slower since it has to rototranslate all the available atoms and
142 : (b) variables employing PBCs (such as \ref DISTANCE without `NOPBC`, as in the example above)
143 : are allowed after \ref FIT_TO_TEMPLATE, whereas \ref RMSD expects PBCs to be already solved.
144 : The latter means that before the \ref RMSD statement one should use \ref WRAPAROUND or \ref WHOLEMOLECULES to properly place
145 : the ligand.
146 :
147 :
148 : */
149 : //+ENDPLUMEDOC
150 :
151 :
152 : class FitToTemplate:
153 : public ActionPilot,
154 : public ActionAtomistic,
155 : public ActionWithValue
156 : {
157 : std::string type;
158 : std::vector<double> weights;
159 : std::vector<AtomNumber> aligned;
160 : Vector center;
161 : Vector shift;
162 : // optimal alignment related stuff
163 : PLMD::RMSD* rmsd;
164 : Tensor rotation;
165 : Matrix< std::vector<Vector> > drotdpos;
166 : std::vector<Vector> positions;
167 : std::vector<Vector> DDistDRef;
168 : std::vector<Vector> ddistdpos;
169 : std::vector<Vector> centeredpositions;
170 : Vector center_positions;
171 :
172 :
173 : public:
174 : explicit FitToTemplate(const ActionOptions&ao);
175 : ~FitToTemplate();
176 : static void registerKeywords( Keywords& keys );
177 : void calculate();
178 : void apply();
179 0 : unsigned getNumberOfDerivatives() {plumed_merror("You should not call this function");};
180 : };
181 :
182 6460 : PLUMED_REGISTER_ACTION(FitToTemplate,"FIT_TO_TEMPLATE")
183 :
184 9 : void FitToTemplate::registerKeywords( Keywords& keys ) {
185 9 : Action::registerKeywords( keys );
186 9 : ActionAtomistic::registerKeywords( keys );
187 45 : keys.add("compulsory","STRIDE","1","the frequency with which molecules are reassembled. Unless you are completely certain about what you are doing leave this set equal to 1!");
188 36 : keys.add("compulsory","REFERENCE","a file in pdb format containing the reference structure and the atoms involved in the CV.");
189 45 : keys.add("compulsory","TYPE","SIMPLE","the manner in which RMSD alignment is performed. Should be OPTIMAL or SIMPLE.");
190 9 : }
191 :
192 8 : FitToTemplate::FitToTemplate(const ActionOptions&ao):
193 : Action(ao),
194 : ActionPilot(ao),
195 : ActionAtomistic(ao),
196 : ActionWithValue(ao),
197 40 : rmsd(NULL)
198 : {
199 : string reference;
200 16 : parse("REFERENCE",reference);
201 8 : type.assign("SIMPLE");
202 16 : parse("TYPE",type);
203 :
204 : // if(type!="SIMPLE") error("Only TYPE=SIMPLE is implemented in FIT_TO_TEMPLATE");
205 :
206 8 : checkRead();
207 :
208 16 : PDB pdb;
209 :
210 : // read everything in ang and transform to nm if we are not in natural units
211 16 : if( !pdb.read(reference,plumed.getAtoms().usingNaturalUnits(),0.1/atoms.getUnits().getLength()) )
212 0 : error("missing input file " + reference );
213 :
214 8 : requestAtoms(pdb.getAtomNumbers());
215 :
216 8 : std::vector<Vector> positions=pdb.getPositions();
217 8 : weights=pdb.getOccupancy();
218 8 : aligned=pdb.getAtomNumbers();
219 :
220 :
221 : // normalize weights
222 100 : double n=0.0; for(unsigned i=0; i<weights.size(); ++i) n+=weights[i];
223 8 : if(n==0.0) {
224 0 : error("PDB file " + reference + " has zero weights. Please check the occupancy column.");
225 : }
226 8 : n=1.0/n;
227 100 : for(unsigned i=0; i<weights.size(); ++i) weights[i]*=n;
228 :
229 : // normalize weights for rmsd calculation
230 8 : vector<double> weights_measure=pdb.getBeta();
231 100 : n=0.0; for(unsigned i=0; i<weights_measure.size(); ++i) n+=weights_measure[i]; n=1.0/n;
232 100 : for(unsigned i=0; i<weights_measure.size(); ++i) weights_measure[i]*=n;
233 :
234 : // subtract the center
235 128 : for(unsigned i=0; i<weights.size(); ++i) center+=positions[i]*weights[i];
236 100 : for(unsigned i=0; i<weights.size(); ++i) positions[i]-=center;
237 :
238 12 : if(type=="OPTIMAL" or type=="OPTIMAL-FAST" ) {
239 4 : rmsd=new RMSD();
240 4 : rmsd->set(weights,weights_measure,positions,type,false,false);// note: the reference is shifted now with center in the origin
241 8 : log<<" Method chosen for fitting: "<<rmsd->getMethod()<<" \n";
242 : }
243 : // register the value of rmsd (might be useful sometimes)
244 8 : addValue(); setNotPeriodic();
245 :
246 : doNotRetrieve();
247 :
248 : // this is required so as to allow modifyGlobalForce() to return correct
249 : // also for forces that are not owned (and thus not zeored) by all processors.
250 : allowToAccessGlobalForces();
251 8 : }
252 :
253 :
254 96 : void FitToTemplate::calculate() {
255 :
256 96 : Vector cc;
257 :
258 1200 : for(unsigned i=0; i<aligned.size(); ++i) {
259 336 : cc+=weights[i]*modifyPosition(aligned[i]);
260 : }
261 :
262 192 : if (type=="SIMPLE") {
263 48 : shift=center-cc;
264 48 : setValue(shift.modulo());
265 12792 : for(unsigned i=0; i<getTotAtoms(); i++) {
266 : Vector & ato (modifyPosition(AtomNumber::index(i)));
267 6372 : ato+=shift;
268 : }
269 : }
270 48 : else if( type=="OPTIMAL" or type=="OPTIMAL-FAST") {
271 : // we store positions here to be used in apply()
272 : // notice that in apply() it is not guaranteed that positions are still equal to their value here
273 : // since they could have been changed by a subsequent FIT_TO_TEMPLATE
274 96 : positions.resize(aligned.size());
275 816 : for (unsigned i=0; i<aligned.size(); i++) positions[i]=modifyPosition(aligned[i]);
276 :
277 : // specific stuff that provides all that is needed
278 48 : double r=rmsd->calc_FitElements( positions, rotation, drotdpos, centeredpositions, center_positions);
279 48 : setValue(r);
280 12720 : for(unsigned i=0; i<getTotAtoms(); i++) {
281 : Vector & ato (modifyPosition(AtomNumber::index(i)));
282 6336 : ato=matmul(rotation,ato-center_positions)+center;
283 : }
284 : // rotate box
285 : Pbc & pbc(modifyGlobalPbc());
286 48 : pbc.setBox(matmul(pbc.getBox(),transpose(rotation)));
287 : }
288 :
289 96 : }
290 :
291 96 : void FitToTemplate::apply() {
292 192 : if (type=="SIMPLE") {
293 48 : Vector totForce;
294 12792 : for(unsigned i=0; i<getTotAtoms(); i++) {
295 6372 : totForce+=modifyGlobalForce(AtomNumber::index(i));
296 : }
297 : Tensor & vv(modifyGlobalVirial());
298 48 : vv+=Tensor(center,totForce);
299 384 : for(unsigned i=0; i<aligned.size(); ++i) {
300 : Vector & ff(modifyGlobalForce(aligned[i]));
301 96 : ff-=totForce*weights[i];
302 : }
303 48 : } else if ( type=="OPTIMAL" or type=="OPTIMAL-FAST") {
304 48 : Vector totForce;
305 12720 : for(unsigned i=0; i<getTotAtoms(); i++) {
306 : Vector & f(modifyGlobalForce(AtomNumber::index(i)));
307 : // rotate back forces
308 6336 : f=matmul(transpose(rotation),f);
309 : // accumulate rotated c.o.m. forces - this is already in the non rotated reference frame
310 6336 : totForce+=f;
311 : }
312 : Tensor& virial(modifyGlobalVirial());
313 : // notice that an extra Tensor(center,matmul(rotation,totForce)) is required to
314 : // compute the derivatives of the rotation with respect to center
315 48 : Tensor ww=matmul(transpose(rotation),virial+Tensor(center,matmul(rotation,totForce)));
316 : // rotate back virial
317 48 : virial=matmul(transpose(rotation),matmul(virial,rotation));
318 :
319 : // now we compute the force due to alignment
320 576 : for(unsigned i=0; i<aligned.size(); i++) {
321 240 : Vector g;
322 960 : for(unsigned k=0; k<3; k++) {
323 : // this could be made faster computing only the diagonal of d
324 720 : Tensor d=matmul(ww,RMSD::getMatrixFromDRot(drotdpos,i,k));
325 720 : g[k]=(d(0,0)+d(1,1)+d(2,2));
326 : }
327 : // here is the extra contribution
328 720 : modifyGlobalForce(aligned[i])+=-g-weights[i]*totForce;
329 : // here it the contribution to the virial
330 : // notice that here we can use absolute positions since, for the alignment to be defined,
331 : // positions should be in one well defined periodic image
332 480 : virial+=extProduct(positions[i],g);
333 : }
334 : // finally, correction to the virial
335 48 : virial+=extProduct(matmul(transpose(rotation),center),totForce);
336 : }
337 96 : }
338 :
339 32 : FitToTemplate::~FitToTemplate() {
340 8 : if(rmsd) delete rmsd;
341 16 : }
342 :
343 : }
344 4839 : }
|
