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 : #include <memory>
41 :
42 : using namespace std;
43 :
44 : namespace PLMD {
45 : namespace generic {
46 :
47 : //+PLUMEDOC GENERIC FIT_TO_TEMPLATE
48 : /*
49 : This action is used to align a molecule to a template.
50 :
51 : This can be used to move the coordinates stored in plumed
52 : so as to be aligned with a provided template in PDB format. Pdb should contain
53 : also weights for alignment (see the format of PDB files used e.g. for \ref RMSD).
54 : Make sure your PDB file is correctly formatted as explained \ref pdbreader "in this page".
55 : Weights for displacement are ignored, since no displacement is computed here.
56 : Notice that all atoms (not only those in the template) are aligned.
57 : To see what effect try
58 : the \ref DUMPATOMS directive to output the atomic positions.
59 :
60 : Also notice that PLUMED propagate forces correctly so that you can add a bias on a CV computed
61 : after alignment. For many CVs this has no effect, but in some case the alignment can
62 : change the result. Examples are:
63 : - \ref POSITION CV since it is affected by a rigid shift of the system.
64 : - \ref DISTANCE CV with COMPONENTS. Since the alignment could involve a rotation (with TYPE=OPTIMAL) the actual components could be different
65 : from the original ones.
66 : - \ref CELL components for a similar reason.
67 : - \ref DISTANCE from a \ref FIXEDATOM, provided the fixed atom is introduced _after_ the \ref FIT_TO_TEMPLATE action.
68 :
69 : \attention
70 : The implementation of TYPE=OPTIMAL is available but should be considered in testing phase. Please report any
71 : strange behavior.
72 :
73 : \attention
74 : This directive modifies the stored position at the precise moment
75 : it is executed. This means that only collective variables
76 : which are below it in the input script will see the corrected positions.
77 : As a general rule, put it at the top of the input file. Also, unless you
78 : know exactly what you are doing, leave the default stride (1), so that
79 : this action is performed at every MD step.
80 :
81 : When running with periodic boundary conditions, the atoms should be
82 : in the proper periodic image. This is done automatically since PLUMED 2.5,
83 : by considering the ordered list of atoms and rebuilding the molecules using a procedure
84 : that is equivalent to that done in \ref WHOLEMOLECULES . Notice that
85 : rebuilding is local to this action. This is different from \ref WHOLEMOLECULES
86 : which actually modifies the coordinates stored in PLUMED.
87 :
88 : In case you want to recover the old behavior you should use the NOPBC flag.
89 : In that case you need to take care that atoms are in the correct
90 : periodic image.
91 :
92 : \par Examples
93 :
94 : Align the atomic position to a template then print them.
95 : The following example is only translating the system so as
96 : to align the center of mass of a molecule to the one in the reference
97 : structure `ref.pdb`:
98 : \plumedfile
99 : # dump coordinates before fitting, to see the difference:
100 : DUMPATOMS FILE=dump-before.xyz ATOMS=1-20
101 :
102 : # fit coordinates to ref.pdb template
103 : # this is a "TYPE=SIMPLE" fit, so that only translations are used.
104 : FIT_TO_TEMPLATE STRIDE=1 REFERENCE=ref.pdb TYPE=SIMPLE
105 :
106 : # dump coordinates after fitting, to see the difference:
107 : DUMPATOMS FILE=dump-after.xyz ATOMS=1-20
108 : \endplumedfile
109 :
110 : The following example instead performs a rototranslational fit.
111 : \plumedfile
112 : # dump coordinates before fitting, to see the difference:
113 : DUMPATOMS FILE=dump-before.xyz ATOMS=1-20
114 :
115 : # fit coordinates to ref.pdb template
116 : # this is a "TYPE=OPTIMAL" fit, so that rototranslations are used.
117 : FIT_TO_TEMPLATE STRIDE=1 REFERENCE=ref.pdb TYPE=OPTIMAL
118 :
119 : # dump coordinates after fitting, to see the difference:
120 : DUMPATOMS FILE=dump-after.xyz ATOMS=1-20
121 : \endplumedfile
122 :
123 : In both these cases the reference structure should be provided in a reference pdb file such as the one below:
124 :
125 : \auxfile{ref.pdb}
126 : ATOM 8 HT3 ALA 2 -1.480 -1.560 1.212 1.00 1.00 DIA H
127 : ATOM 9 CAY ALA 2 -0.096 2.144 -0.669 1.00 1.00 DIA C
128 : ATOM 10 HY1 ALA 2 0.871 2.385 -0.588 1.00 1.00 DIA H
129 : ATOM 12 HY3 ALA 2 -0.520 2.679 -1.400 1.00 1.00 DIA H
130 : ATOM 14 OY ALA 2 -1.139 0.931 -0.973 1.00 1.00 DIA O
131 : END
132 : \endauxfile
133 :
134 : In the following example you see two completely equivalent way
135 : to restrain an atom close to a position that is defined in the reference
136 : frame of an aligned molecule. You could for instance use this command to calculate the
137 : position of the center of mass of a ligand after having aligned the atoms to the reference
138 : frame of the protein that is determined by aligning the atoms in the protein to the coordinates
139 : provided in the file ref.pdb
140 : \plumedfile
141 : # center of the ligand:
142 : center: CENTER ATOMS=100-110
143 :
144 : FIT_TO_TEMPLATE REFERENCE=ref.pdb TYPE=OPTIMAL
145 :
146 : # place a fixed atom in the protein reference coordinates:
147 : fix: FIXEDATOM AT=1.0,1.1,1.0
148 :
149 : # take the distance between the fixed atom and the center of the ligand
150 : d: DISTANCE ATOMS=center,fix
151 :
152 : # apply a restraint
153 : RESTRAINT ARG=d AT=0.0 KAPPA=100.0
154 : \endplumedfile
155 :
156 : Notice that you could have obtained an (almost) identical result by adding a fictitious
157 : atom to `ref.pdb` with the serial number corresponding to the atom labelled `center` (there is no automatic way
158 : to get it, but in this example it should be the number of atoms of the system plus one),
159 : and properly setting the weights for alignment and displacement in \ref RMSD.
160 : There are two differences to be expected:
161 : (ab) \ref FIT_TO_TEMPLATE might be slower since it has to rototranslate all the available atoms and
162 : (b) variables employing periodic boundary conditions (such as \ref DISTANCE without `NOPBC`, as in the example above)
163 : are allowed after \ref FIT_TO_TEMPLATE, whereas \ref RMSD expects the issues related to the periodic boundary conditions to be already solved.
164 : The latter means that before the \ref RMSD statement one should use \ref WRAPAROUND or \ref WHOLEMOLECULES to properly place
165 : the ligand.
166 :
167 :
168 : */
169 : //+ENDPLUMEDOC
170 :
171 :
172 45 : class FitToTemplate:
173 : public ActionPilot,
174 : public ActionAtomistic,
175 : public ActionWithValue
176 : {
177 : std::string type;
178 : bool nopbc;
179 : std::vector<double> weights;
180 : std::vector<AtomNumber> aligned;
181 : Vector center;
182 : Vector shift;
183 : // optimal alignment related stuff
184 : std::unique_ptr<PLMD::RMSD> rmsd;
185 : Tensor rotation;
186 : Matrix< std::vector<Vector> > drotdpos;
187 : // not used anymore (see notes below at doNotRetrieve())
188 : // std::vector<Vector> positions;
189 : std::vector<Vector> DDistDRef;
190 : std::vector<Vector> ddistdpos;
191 : std::vector<Vector> centeredpositions;
192 : Vector center_positions;
193 :
194 :
195 : public:
196 : explicit FitToTemplate(const ActionOptions&ao);
197 : static void registerKeywords( Keywords& keys );
198 : void calculate() override;
199 : void apply() override;
200 0 : unsigned getNumberOfDerivatives() override {plumed_merror("You should not call this function");};
201 : };
202 :
203 7850 : PLUMED_REGISTER_ACTION(FitToTemplate,"FIT_TO_TEMPLATE")
204 :
205 10 : void FitToTemplate::registerKeywords( Keywords& keys ) {
206 10 : Action::registerKeywords( keys );
207 10 : ActionAtomistic::registerKeywords( keys );
208 50 : 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!");
209 40 : keys.add("compulsory","REFERENCE","a file in pdb format containing the reference structure and the atoms involved in the CV.");
210 50 : keys.add("compulsory","TYPE","SIMPLE","the manner in which RMSD alignment is performed. Should be OPTIMAL or SIMPLE.");
211 30 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
212 10 : }
213 :
214 9 : FitToTemplate::FitToTemplate(const ActionOptions&ao):
215 : Action(ao),
216 : ActionPilot(ao),
217 : ActionAtomistic(ao),
218 : ActionWithValue(ao),
219 36 : nopbc(false)
220 : {
221 : string reference;
222 18 : parse("REFERENCE",reference);
223 9 : type.assign("SIMPLE");
224 18 : parse("TYPE",type);
225 :
226 18 : parseFlag("NOPBC",nopbc);
227 : // if(type!="SIMPLE") error("Only TYPE=SIMPLE is implemented in FIT_TO_TEMPLATE");
228 :
229 9 : checkRead();
230 :
231 18 : PDB pdb;
232 :
233 : // read everything in ang and transform to nm if we are not in natural units
234 18 : if( !pdb.read(reference,plumed.getAtoms().usingNaturalUnits(),0.1/atoms.getUnits().getLength()) )
235 0 : error("missing input file " + reference );
236 :
237 9 : requestAtoms(pdb.getAtomNumbers());
238 18 : log.printf(" found %z atoms in input \n",pdb.getAtomNumbers().size());
239 9 : log.printf(" with indices : ");
240 75 : for(unsigned i=0; i<pdb.getAtomNumbers().size(); ++i) {
241 33 : if(i%25==0) log<<"\n";
242 66 : log.printf("%d ",pdb.getAtomNumbers()[i].serial());
243 : }
244 9 : log.printf("\n");
245 :
246 9 : std::vector<Vector> positions=pdb.getPositions();
247 9 : weights=pdb.getOccupancy();
248 9 : aligned=pdb.getAtomNumbers();
249 :
250 :
251 : // normalize weights
252 75 : double n=0.0; for(unsigned i=0; i<weights.size(); ++i) n+=weights[i];
253 9 : if(n==0.0) {
254 0 : error("PDB file " + reference + " has zero weights. Please check the occupancy column.");
255 : }
256 9 : n=1.0/n;
257 84 : for(unsigned i=0; i<weights.size(); ++i) weights[i]*=n;
258 :
259 : // normalize weights for rmsd calculation
260 9 : vector<double> weights_measure=pdb.getBeta();
261 75 : n=0.0; for(unsigned i=0; i<weights_measure.size(); ++i) n+=weights_measure[i]; n=1.0/n;
262 84 : for(unsigned i=0; i<weights_measure.size(); ++i) weights_measure[i]*=n;
263 :
264 : // subtract the center
265 108 : for(unsigned i=0; i<weights.size(); ++i) center+=positions[i]*weights[i];
266 75 : for(unsigned i=0; i<weights.size(); ++i) positions[i]-=center;
267 :
268 13 : if(type=="OPTIMAL" or type=="OPTIMAL-FAST" ) {
269 5 : rmsd.reset(new RMSD());
270 5 : rmsd->set(weights,weights_measure,positions,type,false,false);// note: the reference is shifted now with center in the origin
271 10 : log<<" Method chosen for fitting: "<<rmsd->getMethod()<<" \n";
272 : }
273 9 : if(nopbc) {
274 1 : log<<" Ignoring PBCs when doing alignment, make sure your molecule is whole!<n";
275 : }
276 : // register the value of rmsd (might be useful sometimes)
277 9 : addValue(); setNotPeriodic();
278 :
279 : // I remove this optimization now in order to use makeWhole()
280 : // Notice that for FIT_TO_TEMPLATE TYPE=OPTIMAL a copy was made anyway
281 : // (due to the need to store position to propagate forces on rotational matrix later)
282 : // For FIT_TO_TEMPLATE TYPE=SIMPLE in principle we could use it and write an ad hoc
283 : // version of makeWhole that only computes the center. Too lazy to do it now.
284 : // In case we do it later, remember that uncommenting this line means that
285 : // getPositions will not work anymore! GB
286 : // doNotRetrieve();
287 :
288 : // this is required so as to allow modifyGlobalForce() to return correct
289 : // also for forces that are not owned (and thus not zeored) by all processors.
290 : allowToAccessGlobalForces();
291 9 : }
292 :
293 :
294 108 : void FitToTemplate::calculate() {
295 :
296 108 : if(!nopbc) makeWhole();
297 :
298 216 : if (type=="SIMPLE") {
299 48 : Vector cc;
300 :
301 288 : for(unsigned i=0; i<aligned.size(); ++i) {
302 192 : cc+=weights[i]*getPosition(i);
303 : }
304 :
305 48 : shift=center-cc;
306 48 : setValue(shift.modulo());
307 12840 : for(unsigned i=0; i<getTotAtoms(); i++) {
308 : Vector & ato (modifyGlobalPosition(AtomNumber::index(i)));
309 6372 : ato+=shift;
310 : }
311 : }
312 60 : else if( type=="OPTIMAL" or type=="OPTIMAL-FAST") {
313 : // specific stuff that provides all that is needed
314 120 : double r=rmsd->calc_FitElements( getPositions(), rotation, drotdpos, centeredpositions, center_positions);
315 60 : setValue(r);
316 15948 : for(unsigned i=0; i<getTotAtoms(); i++) {
317 : Vector & ato (modifyGlobalPosition(AtomNumber::index(i)));
318 7944 : ato=matmul(rotation,ato-center_positions)+center;
319 : }
320 : // rotate box
321 : Pbc & pbc(modifyGlobalPbc());
322 60 : pbc.setBox(matmul(pbc.getBox(),transpose(rotation)));
323 : }
324 :
325 108 : }
326 :
327 108 : void FitToTemplate::apply() {
328 216 : if (type=="SIMPLE") {
329 48 : Vector totForce;
330 12840 : for(unsigned i=0; i<getTotAtoms(); i++) {
331 6372 : totForce+=modifyGlobalForce(AtomNumber::index(i));
332 : }
333 48 : Tensor & vv(modifyGlobalVirial());
334 48 : vv+=Tensor(center,totForce);
335 288 : for(unsigned i=0; i<aligned.size(); ++i) {
336 : Vector & ff(modifyGlobalForce(aligned[i]));
337 96 : ff-=totForce*weights[i];
338 : }
339 60 : } else if ( type=="OPTIMAL" or type=="OPTIMAL-FAST") {
340 60 : Vector totForce;
341 16008 : for(unsigned i=0; i<getTotAtoms(); i++) {
342 : Vector & f(modifyGlobalForce(AtomNumber::index(i)));
343 : // rotate back forces
344 7944 : f=matmul(transpose(rotation),f);
345 : // accumulate rotated c.o.m. forces - this is already in the non rotated reference frame
346 7944 : totForce+=f;
347 : }
348 60 : Tensor& virial(modifyGlobalVirial());
349 : // notice that an extra Tensor(center,matmul(rotation,totForce)) is required to
350 : // compute the derivatives of the rotation with respect to center
351 60 : Tensor ww=matmul(transpose(rotation),virial+Tensor(center,matmul(rotation,totForce)));
352 : // rotate back virial
353 60 : virial=matmul(transpose(rotation),matmul(virial,rotation));
354 :
355 : // now we compute the force due to alignment
356 720 : for(unsigned i=0; i<aligned.size(); i++) {
357 300 : Vector g;
358 1200 : for(unsigned k=0; k<3; k++) {
359 : // this could be made faster computing only the diagonal of d
360 900 : Tensor d=matmul(ww,RMSD::getMatrixFromDRot(drotdpos,i,k));
361 900 : g[k]=(d(0,0)+d(1,1)+d(2,2));
362 : }
363 : // here is the extra contribution
364 900 : modifyGlobalForce(aligned[i])+=-g-weights[i]*totForce;
365 : // here it the contribution to the virial
366 : // notice that here we can use absolute positions since, for the alignment to be defined,
367 : // positions should be in one well defined periodic image
368 600 : virial+=extProduct(getPosition(i),g);
369 : }
370 : // finally, correction to the virial
371 60 : virial+=extProduct(matmul(transpose(rotation),center),totForce);
372 : }
373 108 : }
374 :
375 : }
376 5874 : }
|
