Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2014-2023 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/PlumedMain.h"
32 : #include "core/ActionSet.h"
33 : #include "core/GenericMolInfo.h"
34 : #include "core/PbcAction.h"
35 : #include "tools/PDB.h"
36 : #include "tools/Pbc.h"
37 :
38 : #include <vector>
39 : #include <string>
40 : #include <memory>
41 :
42 : namespace PLMD {
43 : namespace generic {
44 :
45 : //+PLUMEDOC GENERIC FIT_TO_TEMPLATE
46 : /*
47 : This action is used to align a molecule to a template.
48 :
49 : This can be used to move the coordinates stored in plumed
50 : so as to be aligned with a provided template in PDB format. Pdb should contain
51 : also weights for alignment (see the format of PDB files used e.g. for [RMSD](RMSD.md)).
52 : Make sure your PDB file is correctly formatted as explained in the documentation for [PDB2CONSTANT](PDB2CONSTANT.md).
53 : Weights for displacement are ignored, since no displacement is computed here.
54 : Notice that all atoms (not only those in the template) are aligned.
55 : To see what effect try
56 : the [DUMPATOMS](DUMPATOMS.md) directive to output the atomic positions.
57 :
58 : Also notice that PLUMED propagate forces correctly so that you can add a bias on a CV computed
59 : after alignment. For many CVs this has no effect, but in some case the alignment can
60 : change the result. Examples are:
61 : - [POSITION](POSITION.md) CV since it is affected by a rigid shift of the system.
62 : - [DISTANCE](DISTANCE.md) CV with COMPONENTS. Since the alignment could involve a rotation (with TYPE=OPTIMAL) the actual components could be different
63 : from the original ones.
64 : - [CELL](CELL.md) components for a similar reason.
65 : - [DISTANCE](DISTANCE.md) from a [FIXEDATOM](FIXEDATOM.md), provided the fixed atom is introduced _after_ the FIT_TO_TEMPLATE action.
66 :
67 : !!! caution ""
68 :
69 : The implementation of TYPE=OPTIMAL is available but should be considered in testing phase. Please report any
70 : strange behavior.
71 :
72 : !!! caution "directive modifies stored positions"
73 :
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 [WHOLEMOLECULES](WHOLEMOLECULES.md) . Notice that
85 : rebuilding is local to this action. This is different from [WHOLEMOLECULES](WHOLEMOLECULES.md)
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 : ## 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 :
99 : ```plumed
100 : #SETTINGS INPUTFILES=regtest/basic/rt63/align.pdb
101 : # dump coordinates before fitting, to see the difference:
102 : DUMPATOMS FILE=dump-before.xyz ATOMS=1-20
103 :
104 : # fit coordinates to ref.pdb template
105 : # this is a "TYPE=SIMPLE" fit, so that only translations are used.
106 : FIT_TO_TEMPLATE STRIDE=1 REFERENCE=regtest/basic/rt63/align.pdb TYPE=SIMPLE
107 :
108 : # dump coordinates after fitting, to see the difference:
109 : DUMPATOMS FILE=dump-after.xyz ATOMS=1-20
110 : ```
111 :
112 : The following example instead performs a rototranslational fit.
113 :
114 : ```plumed
115 : #SETTINGS INPUTFILES=regtest/basic/rt63/align.pdb
116 : # dump coordinates before fitting, to see the difference:
117 : DUMPATOMS FILE=dump-before.xyz ATOMS=1-20
118 :
119 : # fit coordinates to ref.pdb template
120 : # this is a "TYPE=OPTIMAL" fit, so that rototranslations are used.
121 : FIT_TO_TEMPLATE STRIDE=1 REFERENCE=regtest/basic/rt63/align.pdb TYPE=OPTIMAL
122 :
123 : # dump coordinates after fitting, to see the difference:
124 : DUMPATOMS FILE=dump-after.xyz ATOMS=1-20
125 : ```
126 :
127 : In the following example you see two completely equivalent way
128 : to restrain an atom close to a position that is defined in the reference
129 : frame of an aligned molecule. You could for instance use this command to calculate the
130 : position of the center of mass of a ligand after having aligned the atoms to the reference
131 : frame of the protein that is determined by aligning the atoms in the protein to the coordinates
132 : provided in the file ref.pdb
133 :
134 : ```plumed
135 : #SETTINGS INPUTFILES=regtest/basic/rt63/align.pdb
136 : # center of the ligand:
137 : center: CENTER ATOMS=100-110
138 :
139 : FIT_TO_TEMPLATE REFERENCE=regtest/basic/rt63/align.pdb TYPE=OPTIMAL
140 :
141 : # place a fixed atom in the protein reference coordinates:
142 : fix: FIXEDATOM AT=1.0,1.1,1.0
143 :
144 : # take the distance between the fixed atom and the center of the ligand
145 : d: DISTANCE ATOMS=center,fix
146 :
147 : # apply a restraint
148 : RESTRAINT ARG=d AT=0.0 KAPPA=100.0
149 : ```
150 :
151 : Notice that you could have obtained an (almost) identical result by adding a fictitious
152 : atom to `ref.pdb` with the serial number corresponding to the atom labelled `center` (there is no automatic way
153 : to get it, but in this example it should be the number of atoms of the system plus one),
154 : and properly setting the weights for alignment and displacement in [RMSD](RMSD.md).
155 : There are two differences to be expected:
156 : (ab) FIT_TO_TEMPLATE might be slower since it has to rototranslate all the available atoms and
157 : (b) variables employing periodic boundary conditions (such as [DISTANCE](DISTANCE.md) without `NOPBC`, as in the example above)
158 : are allowed after FIT_TO_TEMPLATE, whereas [RMSD](RMSD.md) expects the issues related to the periodic boundary conditions to be already solved.
159 : The latter means that before the [RMSD](RMSD.md) statement one should use [WRAPAROUND](WRAPAROUND.md) or [WHOLEMOLECULES](WHOLEMOLECULES.md) to properly place
160 : the ligand.
161 :
162 :
163 : */
164 : //+ENDPLUMEDOC
165 :
166 :
167 : class FitToTemplate:
168 : public ActionPilot,
169 : public ActionAtomistic,
170 : public ActionWithValue {
171 : std::string type;
172 : bool nopbc;
173 : std::vector<double> weights;
174 : std::vector<std::pair<std::size_t,std::size_t> > p_aligned;
175 : Vector center;
176 : Vector shift;
177 : // optimal alignment related stuff
178 : std::unique_ptr<PLMD::RMSD> rmsd;
179 : Tensor rotation;
180 : Matrix< std::vector<Vector> > drotdpos;
181 : // not used anymore (see notes below at doNotRetrieve())
182 : // std::vector<Vector> positions;
183 : std::vector<Vector> DDistDRef;
184 : std::vector<Vector> ddistdpos;
185 : std::vector<Vector> centeredpositions;
186 : Vector center_positions;
187 : // Copy of the box value
188 : Value* boxValue;
189 : PbcAction* pbc_action;
190 : public:
191 : explicit FitToTemplate(const ActionOptions&ao);
192 : static void registerKeywords( Keywords& keys );
193 99 : bool actionHasForces() override {
194 99 : return true;
195 : }
196 : void calculate() override;
197 : void apply() override;
198 0 : unsigned getNumberOfDerivatives() override {
199 0 : plumed_merror("You should not call this function");
200 : };
201 : };
202 :
203 : PLUMED_REGISTER_ACTION(FitToTemplate,"FIT_TO_TEMPLATE")
204 :
205 13 : void FitToTemplate::registerKeywords( Keywords& keys ) {
206 13 : Action::registerKeywords( keys );
207 13 : ActionAtomistic::registerKeywords( keys );
208 13 : 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 13 : keys.add("compulsory","REFERENCE","a file in pdb format containing the reference structure and the atoms involved in the CV.");
210 13 : keys.add("compulsory","TYPE","SIMPLE","the manner in which RMSD alignment is performed. Should be OPTIMAL or SIMPLE.");
211 13 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
212 26 : keys.setValueDescription("scalar","the RMSD distance");
213 13 : }
214 :
215 9 : FitToTemplate::FitToTemplate(const ActionOptions&ao):
216 : Action(ao),
217 : ActionPilot(ao),
218 : ActionAtomistic(ao),
219 : ActionWithValue(ao),
220 9 : nopbc(false) {
221 : std::string reference;
222 9 : parse("REFERENCE",reference);
223 9 : type.assign("SIMPLE");
224 9 : parse("TYPE",type);
225 :
226 9 : parseFlag("NOPBC",nopbc);
227 9 : if(type!="SIMPLE" && type!="OPTIMAL" && type!="OPTIMAL-FAST") {
228 0 : error("unknown TYPE " + type + ": use SIMPLE, OPTIMAL or OPTIMAL-FAST");
229 : }
230 :
231 9 : checkRead();
232 :
233 9 : PDB pdb;
234 :
235 : // read everything in ang and transform to nm if we are not in natural units
236 9 : if( !pdb.read(reference,usingNaturalUnits(),0.1/getUnits().getLength()) ) {
237 0 : error("missing input file " + reference );
238 : }
239 :
240 9 : requestAtoms(pdb.getAtomNumbers());
241 9 : log.printf(" found %zu atoms in input \n",pdb.getAtomNumbers().size());
242 9 : log.printf(" with indices : ");
243 42 : for(unsigned i=0; i<pdb.getAtomNumbers().size(); ++i) {
244 33 : if(i%25==0) {
245 9 : log<<"\n";
246 : }
247 33 : log.printf("%d ",pdb.getAtomNumbers()[i].serial());
248 : }
249 9 : log.printf("\n");
250 :
251 9 : std::vector<Vector> myPositions=pdb.getPositions();
252 9 : weights=pdb.getOccupancy();
253 9 : std::vector<AtomNumber> aligned=pdb.getAtomNumbers();
254 9 : p_aligned.resize( aligned.size() );
255 42 : for(unsigned i=0; i<aligned.size(); ++i) {
256 33 : p_aligned[i] = getValueIndices( aligned[i] );
257 : }
258 :
259 :
260 : // normalize weights
261 : double n=0.0;
262 42 : for(unsigned i=0; i<weights.size(); ++i) {
263 33 : n+=weights[i];
264 : }
265 9 : if(n==0.0) {
266 0 : error("PDB file " + reference + " has zero weights. Please check the occupancy column.");
267 : }
268 9 : n=1.0/n;
269 42 : for(unsigned i=0; i<weights.size(); ++i) {
270 33 : weights[i]*=n;
271 : }
272 :
273 : // normalize weights for rmsd calculation
274 9 : std::vector<double> weights_measure=pdb.getBeta();
275 : n=0.0;
276 42 : for(unsigned i=0; i<weights_measure.size(); ++i) {
277 33 : n+=weights_measure[i];
278 : }
279 9 : if(n==0.0 && type!="SIMPLE") {
280 0 : error("PDB file " + reference + " has zero weights for the RMSD. Please check the beta column.");
281 : }
282 9 : n=1.0/n;
283 42 : for(unsigned i=0; i<weights_measure.size(); ++i) {
284 33 : weights_measure[i]*=n;
285 : }
286 :
287 : // subtract the center
288 42 : for(unsigned i=0; i<weights.size(); ++i) {
289 33 : center+=myPositions[i]*weights[i];
290 : }
291 42 : for(unsigned i=0; i<weights.size(); ++i) {
292 : myPositions[i]-=center;
293 : }
294 :
295 9 : if(type=="OPTIMAL" or type=="OPTIMAL-FAST" ) {
296 5 : rmsd=Tools::make_unique<RMSD>();
297 5 : rmsd->set(weights,weights_measure,myPositions,type,false,false);// note: the reference is shifted now with center in the origin
298 10 : log<<" Method chosen for fitting: "<<rmsd->getMethod()<<" \n";
299 : }
300 9 : if(nopbc) {
301 1 : log<<" Ignoring PBCs when doing alignment, make sure your molecule is whole!<n";
302 : }
303 : // register the value of rmsd (might be useful sometimes)
304 9 : addValue();
305 9 : setNotPeriodic();
306 :
307 : // I remove this optimization now in order to use makeWhole()
308 : // Notice that for FIT_TO_TEMPLATE TYPE=OPTIMAL a copy was made anyway
309 : // (due to the need to store position to propagate forces on rotational matrix later)
310 : // For FIT_TO_TEMPLATE TYPE=SIMPLE in principle we could use it and write an ad hoc
311 : // version of makeWhole that only computes the center. Too lazy to do it now.
312 : // In case we do it later, remember that uncommenting this line means that
313 : // getPositions will not work anymore! GB
314 : // doNotRetrieve();
315 :
316 : // this is required so as to allow modifyGlobalForce() to return correct
317 : // also for forces that are not owned (and thus not zeored) by all processors.
318 9 : pbc_action=plumed.getActionSet().selectWithLabel<PbcAction*>("Box");
319 9 : if( !pbc_action ) {
320 0 : error("cannot align box has not been set");
321 : }
322 9 : boxValue=pbc_action->copyOutput(0);
323 18 : }
324 :
325 :
326 108 : void FitToTemplate::calculate() {
327 :
328 108 : if(!nopbc) {
329 96 : makeWhole();
330 : }
331 :
332 108 : if (type=="SIMPLE") {
333 : Vector cc;
334 :
335 144 : for(unsigned i=0; i<p_aligned.size(); ++i) {
336 96 : cc+=weights[i]*getPosition(i);
337 : }
338 :
339 48 : shift=center-cc;
340 48 : setValue(shift.modulo());
341 48 : translateAllGlobalPositions(shift);
342 60 : } else if( type=="OPTIMAL" or type=="OPTIMAL-FAST") {
343 : // specific stuff that provides all that is needed
344 60 : double r=rmsd->calc_FitElements( getPositions(), rotation, drotdpos, centeredpositions, center_positions);
345 60 : setValue(r);
346 60 : rotateAllGlobalPositions(rotation,center_positions,center);
347 : // rotate box
348 60 : pbc_action->getPbc().setBox(matmul(pbc_action->getPbc().getBox(),transpose(rotation)));
349 : }
350 108 : }
351 :
352 108 : void FitToTemplate::apply() {
353 108 : if (type=="SIMPLE") {
354 48 : Vector totForce=getTotalForceOnAllAtoms();
355 48 : Tensor vv=Tensor(center,totForce);
356 192 : for(unsigned i=0; i<3; ++i)
357 576 : for(unsigned j=0; j<3; ++j) {
358 432 : boxValue->addForce( 3*i+j, vv(i,j) );
359 : }
360 144 : for(unsigned i=0; i<p_aligned.size(); ++i) {
361 96 : addForce( p_aligned[i], -totForce*weights[i]);
362 : }
363 60 : } else if ( type=="OPTIMAL" or type=="OPTIMAL-FAST") {
364 60 : const Tensor rotationT=transpose(rotation);
365 : // rotate back forces and accumulate rotated c.o.m. forces - this is already in the non rotated reference frame
366 60 : Vector totForce=rotateForcesOnAllAtoms(rotationT);
367 60 : Tensor virial;
368 240 : for(unsigned i=0; i<3; ++i)
369 720 : for(unsigned j=0; j<3; ++j) {
370 540 : virial[i][j] = boxValue->getForce( 3*i+j );
371 : }
372 : // notice that an extra Tensor(center,matmul(rotation,totForce)) is required to
373 : // compute the derivatives of the rotation with respect to center
374 60 : Tensor ww=matmul(rotationT,virial+Tensor(center,matmul(rotation,totForce)));
375 : // rotate back virial
376 60 : virial=matmul(rotationT,matmul(virial,rotation));
377 :
378 : // now we compute the force due to alignment
379 360 : for(unsigned i=0; i<p_aligned.size(); i++) {
380 300 : Vector g;
381 1200 : for(unsigned k=0; k<3; k++) {
382 : // trace of matmul(ww,RMSD::getMatrixFromDRot(drotdpos,i,k)): only the diagonal is computed,
383 : // summing in the same order as matmul
384 : double trace=0.0;
385 3600 : for(unsigned l=0; l<3; l++) {
386 : double dll=0.0;
387 10800 : for(unsigned m=0; m<3; m++) {
388 8100 : dll+=ww(l,m)*drotdpos(m,l)[i][k];
389 : }
390 2700 : trace+=dll;
391 : }
392 900 : g[k]=trace;
393 : }
394 : // here is the extra contribution
395 300 : addForce( p_aligned[i], -g-weights[i]*totForce );
396 : // here it the contribution to the virial
397 : // notice that here we can use absolute positions since, for the alignment to be defined,
398 : // positions should be in one well defined periodic image
399 600 : virial+=extProduct(getPosition(i),g);
400 : }
401 : // finally, correction to the virial
402 60 : boxValue->clearInputForce();
403 60 : virial+=extProduct(matmul(rotationT,center),totForce);
404 240 : for(unsigned i=0; i<3; ++i)
405 720 : for(unsigned j=0; j<3; ++j) {
406 540 : boxValue->addForce( 3*i+j, virial(i,j) );
407 : }
408 : }
409 108 : }
410 :
411 : }
412 : }
|
