Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2014-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 "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 :
67 : \attention
68 : The implementation of TYPE=OPTIMAL is available but should be considered in testing phase. Please report any
69 : strange behavior.
70 :
71 : \attention
72 : This directive modifies the stored position at the precise moment
73 : it is executed. This means that only collective variables
74 : which are below it in the input script will see the corrected positions.
75 : As a general rule, put it at the top of the input file. Also, unless you
76 : know exactly what you are doing, leave the default stride (1), so that
77 : this action is performed at every MD step.
78 :
79 : \warning
80 : The molecule used for alignment should be whole.
81 : In case it is broken by the host MD code, please use \ref WHOLEMOLECULES to reconstruct it before \ref FIT_TO_TEMPLATE .
82 :
83 :
84 : \par Examples
85 :
86 : Align the atomic position to a template then print them
87 : \verbatim
88 : # to see the effect, one could dump the atoms before alignment
89 : DUMPATOMS FILE=dump-before.xyz ATOMS=1-20
90 : FIT_TO_TEMPLATE STRIDE=1 REFERENCE=ref.pdb TYPE=SIMPLE
91 : DUMPATOMS FILE=dump-after.xyz ATOMS=1-20
92 : \endverbatim
93 : (see also \ref DUMPATOMS)
94 :
95 :
96 :
97 :
98 : */
99 : //+ENDPLUMEDOC
100 :
101 :
102 : class FitToTemplate:
103 : public ActionPilot,
104 : public ActionAtomistic,
105 : public ActionWithValue
106 : {
107 : std::string type;
108 : std::vector<double> weights;
109 : std::vector<AtomNumber> aligned;
110 : Vector center;
111 : Vector shift;
112 : // optimal alignment related stuff
113 : PLMD::RMSD* rmsd;
114 : Tensor rotation;
115 : Matrix< std::vector<Vector> > drotdpos;
116 : std::vector<Vector> positions;
117 : std::vector<Vector> DDistDRef;
118 : std::vector<Vector> ddistdpos;
119 : std::vector<Vector> centeredpositions;
120 : Vector center_positions;
121 :
122 :
123 : public:
124 : explicit FitToTemplate(const ActionOptions&ao);
125 : ~FitToTemplate();
126 : static void registerKeywords( Keywords& keys );
127 : void calculate();
128 : void apply();
129 0 : unsigned getNumberOfDerivatives() {plumed_merror("You should not call this function");};
130 : };
131 :
132 2531 : PLUMED_REGISTER_ACTION(FitToTemplate,"FIT_TO_TEMPLATE")
133 :
134 9 : void FitToTemplate::registerKeywords( Keywords& keys ) {
135 9 : Action::registerKeywords( keys );
136 9 : ActionAtomistic::registerKeywords( keys );
137 9 : 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!");
138 9 : keys.add("compulsory","REFERENCE","a file in pdb format containing the reference structure and the atoms involved in the CV.");
139 9 : keys.add("compulsory","TYPE","SIMPLE","the manner in which RMSD alignment is performed. Should be OPTIMAL or SIMPLE.");
140 9 : }
141 :
142 8 : FitToTemplate::FitToTemplate(const ActionOptions&ao):
143 : Action(ao),
144 : ActionPilot(ao),
145 : ActionAtomistic(ao),
146 : ActionWithValue(ao),
147 8 : rmsd(NULL)
148 : {
149 8 : string reference;
150 8 : parse("REFERENCE",reference);
151 8 : type.assign("SIMPLE");
152 8 : parse("TYPE",type);
153 :
154 : // if(type!="SIMPLE") error("Only TYPE=SIMPLE is implemented in FIT_TO_TEMPLATE");
155 :
156 8 : checkRead();
157 :
158 16 : PDB pdb;
159 :
160 : // read everything in ang and transform to nm if we are not in natural units
161 8 : if( !pdb.read(reference,plumed.getAtoms().usingNaturalUnits(),0.1/atoms.getUnits().getLength()) )
162 0 : error("missing input file " + reference );
163 :
164 8 : requestAtoms(pdb.getAtomNumbers());
165 :
166 16 : std::vector<Vector> positions=pdb.getPositions();
167 8 : weights=pdb.getOccupancy();
168 8 : aligned=pdb.getAtomNumbers();
169 :
170 :
171 : // normalize weights
172 8 : double n=0.0; for(unsigned i=0; i<weights.size(); ++i) n+=weights[i];
173 8 : if(n==0.0) {
174 0 : error("PDB file " + reference + " has zero weights. Please check the occupancy column.");
175 : }
176 8 : n=1.0/n;
177 8 : for(unsigned i=0; i<weights.size(); ++i) weights[i]*=n;
178 :
179 : // normalize weights for rmsd calculation
180 16 : vector<double> weights_measure=pdb.getBeta();
181 8 : n=0.0; for(unsigned i=0; i<weights_measure.size(); ++i) n+=weights_measure[i]; n=1.0/n;
182 8 : for(unsigned i=0; i<weights_measure.size(); ++i) weights_measure[i]*=n;
183 :
184 : // subtract the center
185 8 : for(unsigned i=0; i<weights.size(); ++i) center+=positions[i]*weights[i];
186 8 : for(unsigned i=0; i<weights.size(); ++i) positions[i]-=center;
187 :
188 8 : if(type=="OPTIMAL" or type=="OPTIMAL-FAST" ) {
189 4 : rmsd=new RMSD();
190 4 : rmsd->set(weights,weights_measure,positions,type,false,false);// note: the reference is shifted now with center in the origin
191 4 : log<<" Method chosen for fitting: "<<rmsd->getMethod()<<" \n";
192 : }
193 : // register the value of rmsd (might be useful sometimes)
194 8 : addValue(); setNotPeriodic();
195 :
196 8 : doNotRetrieve();
197 :
198 : // this is required so as to allow modifyGlobalForce() to return correct
199 : // also for forces that are not owned (and thus not zeored) by all processors.
200 16 : allowToAccessGlobalForces();
201 8 : }
202 :
203 :
204 96 : void FitToTemplate::calculate() {
205 :
206 96 : Vector cc;
207 :
208 432 : for(unsigned i=0; i<aligned.size(); ++i) {
209 336 : cc+=weights[i]*modifyPosition(aligned[i]);
210 : }
211 :
212 96 : if (type=="SIMPLE") {
213 48 : shift=center-cc;
214 48 : setValue(shift.modulo());
215 6420 : for(unsigned i=0; i<getTotAtoms(); i++) {
216 6372 : Vector & ato (modifyPosition(AtomNumber::index(i)));
217 6372 : ato+=shift;
218 : }
219 : }
220 48 : else if( type=="OPTIMAL" or type=="OPTIMAL-FAST") {
221 : // we store positions here to be used in apply()
222 : // notice that in apply() it is not guaranteed that positions are still equal to their value here
223 : // since they could have been changed by a subsequent FIT_TO_TEMPLATE
224 48 : positions.resize(aligned.size());
225 48 : for (unsigned i=0; i<aligned.size(); i++) positions[i]=modifyPosition(aligned[i]);
226 :
227 : // specific stuff that provides all that is needed
228 48 : double r=rmsd->calc_FitElements( positions, rotation, drotdpos, centeredpositions, center_positions);
229 48 : setValue(r);
230 6384 : for(unsigned i=0; i<getTotAtoms(); i++) {
231 6336 : Vector & ato (modifyPosition(AtomNumber::index(i)));
232 6336 : ato=matmul(rotation,ato-center_positions)+center;
233 : }
234 : // rotate box
235 48 : Pbc & pbc(modifyGlobalPbc());
236 48 : pbc.setBox(matmul(pbc.getBox(),transpose(rotation)));
237 : }
238 :
239 96 : }
240 :
241 96 : void FitToTemplate::apply() {
242 96 : if (type=="SIMPLE") {
243 48 : Vector totForce;
244 6420 : for(unsigned i=0; i<getTotAtoms(); i++) {
245 6372 : totForce+=modifyGlobalForce(AtomNumber::index(i));
246 : }
247 48 : Tensor & vv(modifyGlobalVirial());
248 48 : vv+=Tensor(center,totForce);
249 144 : for(unsigned i=0; i<aligned.size(); ++i) {
250 96 : Vector & ff(modifyGlobalForce(aligned[i]));
251 96 : ff-=totForce*weights[i];
252 : }
253 48 : } else if ( type=="OPTIMAL" or type=="OPTIMAL-FAST") {
254 48 : Vector totForce;
255 6384 : for(unsigned i=0; i<getTotAtoms(); i++) {
256 6336 : Vector & f(modifyGlobalForce(AtomNumber::index(i)));
257 : // rotate back forces
258 6336 : f=matmul(transpose(rotation),f);
259 : // accumulate rotated c.o.m. forces - this is already in the non rotated reference frame
260 6336 : totForce+=f;
261 : }
262 48 : Tensor& virial(modifyGlobalVirial());
263 : // notice that an extra Tensor(center,matmul(rotation,totForce)) is required to
264 : // compute the derivatives of the rotation with respect to center
265 48 : Tensor ww=matmul(transpose(rotation),virial+Tensor(center,matmul(rotation,totForce)));
266 : // rotate back virial
267 48 : virial=matmul(transpose(rotation),matmul(virial,rotation));
268 :
269 : // now we compute the force due to alignment
270 288 : for(unsigned i=0; i<aligned.size(); i++) {
271 240 : Vector g;
272 960 : for(unsigned k=0; k<3; k++) {
273 : // this could be made faster computing only the diagonal of d
274 720 : Tensor d=matmul(ww,RMSD::getMatrixFromDRot(drotdpos,i,k));
275 720 : g[k]=(d(0,0)+d(1,1)+d(2,2));
276 : }
277 : // here is the extra contribution
278 240 : modifyGlobalForce(aligned[i])+=-g-weights[i]*totForce;
279 : // here it the contribution to the virial
280 : // notice that here we can use absolute positions since, for the alignment to be defined,
281 : // positions should be in one well defined periodic image
282 240 : virial+=extProduct(positions[i],g);
283 : }
284 : // finally, correction to the virial
285 48 : virial+=extProduct(matmul(transpose(rotation),center),totForce);
286 : }
287 96 : }
288 :
289 32 : FitToTemplate::~FitToTemplate() {
290 8 : if(rmsd) delete rmsd;
291 24 : }
292 :
293 : }
294 2523 : }
|
