Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-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 "RMSDVector.h"
23 : #include "core/PlumedMain.h"
24 : #include "core/ActionRegister.h"
25 : #include "tools/PDB.h"
26 :
27 : //+PLUMEDOC DCOLVAR RMSD_VECTOR
28 : /*
29 : Calculate the RMSD distance between the instaneous configuration and multiple reference configurations
30 :
31 :
32 : \par Examples
33 :
34 : */
35 : //+ENDPLUMEDOC
36 :
37 : namespace PLMD {
38 : namespace colvar {
39 :
40 : PLUMED_REGISTER_ACTION(RMSDVector,"RMSD_VECTOR")
41 :
42 81 : void RMSDVector::registerKeywords(Keywords& keys) {
43 81 : ActionWithVector::registerKeywords(keys);
44 81 : keys.setDisplayName("RMSD");
45 162 : keys.addInputKeyword("compulsory","ARG","vector/matrix","the labels of two actions that you are calculating the RMSD between");
46 81 : keys.add("compulsory","TYPE","SIMPLE","the manner in which RMSD alignment is performed. Should be OPTIMAL or SIMPLE.");
47 81 : keys.add("compulsory","ALIGN","1.0","the weights to use when aligning to the reference structure");
48 81 : keys.add("compulsory","DISPLACE","1.0","the weights to use when calculating the displacement from the reference structure");
49 81 : keys.addFlag("SQUARED",false," This should be set if you want mean squared displacement instead of RMSD ");
50 81 : keys.addFlag("UNORMALIZED",false,"by default the mean sequare deviation or root mean square deviation is calculated. If this option is given no averaging is done");
51 81 : keys.addFlag("DISPLACEMENT",false,"Calculate the vector of displacements instead of the length of this vector");
52 162 : keys.addOutputComponent("disp","DISPLACEMENT","vector/matrix","the vector of displacements for the atoms");
53 162 : keys.addOutputComponent("dist","DISPLACEMENT","scalar/vector","the RMSD distance the atoms have moved");
54 162 : keys.setValueDescription("scalar/vector","a vector containing the RMSD between the instantaneous structure and each of the reference structures that were input");
55 81 : }
56 :
57 49 : RMSDVector::RMSDVector(const ActionOptions&ao):
58 : Action(ao),
59 : ActionWithVector(ao),
60 49 : firststep(true) {
61 49 : if( getPntrToArgument(0)->getRank()!=1 ) {
62 0 : error("first argument should be vector");
63 : }
64 49 : if( getPntrToArgument(1)->getRank()<1 ) {
65 0 : error("second argument should be matrix or a vector");
66 : }
67 49 : if( getPntrToArgument(1)->getRank()==1 ) {
68 7 : if( getPntrToArgument(0)->getNumberOfValues()!=getPntrToArgument(1)->getNumberOfValues() ) {
69 0 : error("mismatch between sizes of input vectors");
70 : }
71 42 : } else if( getPntrToArgument(1)->getRank()==2 ) {
72 42 : if( getPntrToArgument(0)->getNumberOfValues()!=getPntrToArgument(1)->getShape()[1] ) {
73 0 : error("mismatch between sizes of input vectors");
74 : }
75 : }
76 49 : if( getPntrToArgument(0)->getNumberOfValues()%3!=0 ) {
77 0 : error("number of components in input arguments should be multiple of three");
78 : }
79 :
80 49 : unsigned natoms = getPntrToArgument(0)->getNumberOfValues() / 3;
81 49 : type.assign("SIMPLE");
82 49 : parse("TYPE",type);
83 49 : parseFlag("SQUARED",squared);
84 49 : align.resize( natoms );
85 49 : parseVector("ALIGN",align);
86 49 : displace.resize( natoms );
87 49 : parseVector("DISPLACE",displace);
88 49 : bool unorm=false;
89 49 : parseFlag("UNORMALIZED",unorm);
90 49 : norm_weights=!unorm;
91 : double wa=0, wd=0;
92 49 : sqrtdisplace.resize( displace.size() );
93 630 : for(unsigned i=0; i<align.size(); ++i) {
94 581 : wa+=align[i];
95 581 : wd+=displace[i];
96 : }
97 :
98 49 : if( wa>epsilon ) {
99 49 : double iwa = 1. / wa;
100 630 : for(unsigned i=0; i<align.size(); ++i) {
101 581 : align[i] *= iwa;
102 : }
103 : } else {
104 0 : double iwa = 1. / natoms;
105 0 : for(unsigned i=0; i<align.size(); ++i) {
106 0 : align[i] = iwa;
107 : }
108 : }
109 49 : if( wd>epsilon ) {
110 49 : if( !norm_weights ) {
111 : wd = 1;
112 : }
113 49 : double iwd = 1. / wd;
114 630 : for(unsigned i=0; i<align.size(); ++i) {
115 581 : displace[i] *= iwd;
116 : }
117 : } else {
118 0 : double iwd = 1. / natoms;
119 0 : for(unsigned i=0; i<align.size(); ++i) {
120 0 : displace[i] = iwd;
121 : }
122 : }
123 630 : for(unsigned i=0; i<align.size(); ++i) {
124 581 : sqrtdisplace[i] = sqrt(displace[i]);
125 : }
126 :
127 49 : parseFlag("DISPLACEMENT",displacement);
128 49 : if( displacement && (getPntrToArgument(1)->getRank()==1 || getPntrToArgument(1)->getShape()[0]<=1) ) {
129 52 : addComponentWithDerivatives("dist");
130 52 : componentIsNotPeriodic("dist");
131 26 : std::vector<unsigned> shape( 1, getPntrToArgument(0)->getNumberOfValues() );
132 26 : addComponent( "disp", shape );
133 26 : getPntrToComponent(1)->buildDataStore();
134 52 : componentIsNotPeriodic("disp");
135 23 : } else if( displacement ) {
136 2 : std::vector<unsigned> shape( 1, getPntrToArgument(1)->getShape()[0] );
137 2 : addComponent( "dist", shape );
138 2 : getPntrToComponent(0)->buildDataStore();
139 2 : componentIsNotPeriodic("dist");
140 2 : shape.resize(2);
141 2 : shape[0] = getPntrToArgument(1)->getShape()[0];
142 2 : shape[1] = getPntrToArgument(0)->getNumberOfValues();
143 2 : addComponent( "disp", shape );
144 2 : getPntrToComponent(1)->buildDataStore();
145 2 : getPntrToComponent(1)->reshapeMatrixStore( shape[1] );
146 4 : componentIsNotPeriodic("disp");
147 21 : } else if( (getPntrToArgument(1)->getRank()==1 || getPntrToArgument(1)->getShape()[0]==1) ) {
148 13 : addValue();
149 13 : setNotPeriodic();
150 : } else {
151 8 : std::vector<unsigned> shape( 1, getPntrToArgument(1)->getShape()[0] );
152 8 : addValue( shape );
153 8 : setNotPeriodic();
154 : }
155 49 : if( getPntrToArgument(1)->getRank()==1 || getPntrToArgument(1)->getNumberOfValues()==0 ) {
156 7 : myrmsd.resize(1);
157 : } else {
158 42 : myrmsd.resize( getPntrToArgument(1)->getShape()[0] );
159 : }
160 :
161 49 : if( getPntrToArgument(1)->getRank()==1 )
162 7 : log.printf(" calculating RMSD distance between %d atoms. Distance between the avectors of atoms in %s and %s\n",
163 : natoms, getPntrToArgument(0)->getName().c_str(), getPntrToArgument(1)->getName().c_str() );
164 : else
165 42 : log.printf(" calculating RMSD distance between %d sets of %d atoms. Distance between vector %s of atoms and matrix of configurations in %s\n",
166 : getPntrToArgument(1)->getShape()[0], natoms, getPntrToArgument(0)->getName().c_str(), getPntrToArgument(1)->getName().c_str() );
167 49 : log.printf(" method for alignment : %s \n",type.c_str() );
168 49 : if(squared) {
169 28 : log.printf(" chosen to use SQUARED option for MSD instead of RMSD\n");
170 : } else {
171 21 : log.printf(" using periodic boundary conditions\n");
172 : }
173 49 : }
174 :
175 804 : unsigned RMSDVector::getNumberOfDerivatives() {
176 804 : return getPntrToArgument(0)->getNumberOfValues() + getPntrToArgument(1)->getNumberOfValues();
177 : }
178 :
179 10631 : void RMSDVector::setReferenceConfigurations() {
180 10631 : unsigned natoms = getPntrToArgument(0)->getShape()[0] / 3;
181 10631 : Vector center;
182 10631 : std::vector<Vector> pos( natoms );
183 22578 : for(unsigned jconf=0; jconf<myrmsd.size(); ++jconf) {
184 11947 : center.zero();
185 172431 : for(unsigned i=0; i<pos.size(); ++i) {
186 641936 : for(unsigned j=0; j<3; ++j) {
187 481452 : pos[i][j] = getPntrToArgument(1)->get( (3*jconf+j)*pos.size() + i );
188 : }
189 160484 : center+=pos[i]*align[i];
190 : }
191 172431 : for(unsigned i=0; i<pos.size(); ++i) {
192 160484 : pos[i] -= center;
193 : }
194 11947 : myrmsd[jconf].clear();
195 11947 : myrmsd[jconf].set(align,displace,pos,type,true,norm_weights);
196 : }
197 10631 : }
198 :
199 192996 : double RMSDVector::calculateRMSD( const unsigned& current, std::vector<Vector>& pos, std::vector<Vector>& der, std::vector<Vector>& direction ) const {
200 192996 : unsigned natoms = pos.size();
201 2704147 : for(unsigned i=0; i<natoms; ++i) {
202 10044604 : for(unsigned j=0; j<3; ++j) {
203 7533453 : pos[i][j] = getPntrToArgument(0)->get( j*natoms + i );
204 : }
205 : }
206 :
207 192996 : if( displacement && type=="SIMPLE" ) {
208 646 : const Value* myval = getConstPntrToComponent(1);
209 646 : double r = myrmsd[current].simpleAlignment( align, displace, pos, myrmsd[current].getReference(), der, direction, squared );
210 646 : if( !doNotCalculateDerivatives() && myval->forcesWereAdded() ) {
211 33 : Vector comforce;
212 33 : comforce.zero();
213 187 : for(unsigned i=0; i<natoms; i++) {
214 616 : for(unsigned k=0; k<3; ++k) {
215 462 : comforce[k] += align[i]*myval->getForce( (3*current+k)*natoms + i);
216 : }
217 : }
218 187 : for(unsigned i=0; i<natoms; i++) {
219 616 : for(unsigned k=0; k<3; ++k) {
220 462 : direction[i][k] = myval->getForce( (3*current+k)*natoms + i ) - comforce[k];
221 : }
222 : }
223 : }
224 646 : return r;
225 192350 : } else if( displacement ) {
226 55589 : const Value* myval = getConstPntrToComponent(1);
227 55589 : Tensor rot;
228 : Matrix<std::vector<Vector> > DRotDPos(3,3);
229 55589 : std::vector<Vector> centeredpos( natoms ), centeredreference( natoms );
230 55589 : double r = myrmsd[current].calc_PCAelements( pos, der, rot, DRotDPos, direction, centeredpos, centeredreference, squared );
231 55589 : std::vector<Vector> ref( myrmsd[current].getReference() );
232 55589 : if( !doNotCalculateDerivatives() && myval->forcesWereAdded() ) {
233 1671 : Tensor trot=rot.transpose();
234 1671 : double prefactor = 1 / static_cast<double>( natoms );
235 1671 : Vector v1;
236 1671 : v1.zero();
237 22573 : for(unsigned n=0; n<natoms; n++) {
238 20902 : Vector ff;
239 83608 : for(unsigned k=0; k<3; ++k ) {
240 62706 : ff[k] = myval->getForce( (3*current+k)*natoms + n );
241 : }
242 20902 : v1+=prefactor*matmul(trot,ff);
243 : }
244 : // Notice that we use centreredreference here to accumulate the true forces
245 22573 : for(unsigned n=0; n<natoms; n++) {
246 20902 : Vector ff;
247 83608 : for(unsigned k=0; k<3; ++k ) {
248 62706 : ff[k] = myval->getForce( (3*current+k)*natoms + n );
249 : }
250 20902 : centeredreference[n] = sqrtdisplace[n]*( matmul(trot,ff) - v1 );
251 : }
252 6684 : for(unsigned a=0; a<3; a++) {
253 20052 : for(unsigned b=0; b<3; b++) {
254 : double tmp1=0.;
255 203157 : for(unsigned m=0; m<natoms; m++) {
256 188118 : tmp1+=centeredpos[m][b]*myval->getForce( (3*current+a)*natoms + m );
257 : }
258 203157 : for(unsigned i=0; i<natoms; i++) {
259 188118 : centeredreference[i] += sqrtdisplace[i]*tmp1*DRotDPos[a][b][i];
260 : }
261 : }
262 : }
263 : // Now subtract the current force and add on the true force
264 22573 : for(unsigned n=0; n<natoms; n++) {
265 83608 : for(unsigned k=0; k<3; ++k) {
266 62706 : direction[n][k] = centeredreference[n][k];
267 : }
268 : }
269 : } else {
270 759226 : for(unsigned i=0; i<direction.size(); ++i) {
271 705308 : direction[i] = sqrtdisplace[i]*( direction[i] - ref[i] );
272 : }
273 : }
274 : return r;
275 : }
276 136761 : return myrmsd[current].calculate( pos, der, squared );
277 : }
278 :
279 : // calculator
280 17182 : void RMSDVector::calculate() {
281 17182 : if( firststep || !getPntrToArgument(1)->isConstant() ) {
282 10608 : setReferenceConfigurations();
283 10608 : firststep=false;
284 : }
285 :
286 17182 : if( getPntrToComponent(0)->getRank()==0 ) {
287 11796 : unsigned natoms = getPntrToArgument(0)->getShape()[0] / 3;
288 11796 : std::vector<Vector> pos( natoms ), der( natoms ), direction( natoms );
289 11796 : double r = calculateRMSD( 0, pos, der, direction );
290 :
291 11796 : getPntrToComponent(0)->set( r );
292 11796 : if( getNumberOfComponents()==2 ) {
293 11775 : Value* mydisp = getPntrToComponent(1);
294 168204 : for(unsigned i=0; i<natoms; i++) {
295 625716 : for(unsigned j=0; j<3; ++j ) {
296 469287 : mydisp->set( j*natoms+i, direction[i][j] );
297 : }
298 : }
299 : }
300 11796 : if( doNotCalculateDerivatives() ) {
301 : return;
302 : }
303 :
304 612 : Value* myval = getPntrToComponent(0);
305 7472 : for(unsigned i=0; i<natoms; i++) {
306 27440 : for(unsigned j=0; j<3; ++j ) {
307 20580 : myval->setDerivative( j*natoms+i, der[i][j] );
308 : }
309 : }
310 : } else {
311 5386 : runAllTasks();
312 : }
313 : }
314 :
315 160524 : bool RMSDVector::checkForTaskForce( const unsigned& itask, const Value* myval ) const {
316 160524 : if( myval->getRank()<2 ) {
317 137592 : return ActionWithVector::checkForTaskForce( itask, myval );
318 : }
319 22932 : unsigned nelements = myval->getShape()[1], startr = itask*nelements;
320 874692 : for(unsigned j=0; j<nelements; ++j ) {
321 852852 : if( fabs( myval->getForce( startr + j ) )>epsilon ) {
322 : return true;
323 : }
324 : }
325 : return false;
326 : }
327 :
328 17162 : void RMSDVector::apply() {
329 17162 : if( doNotCalculateDerivatives() ) {
330 : return;
331 : }
332 :
333 4980 : if( getPntrToComponent(0)->getRank()==0 ) {
334 612 : std::vector<double> forces( getNumberOfDerivatives(), 0 );
335 612 : bool wasforced = getPntrToComponent(0)->applyForce( forces );
336 :
337 612 : if( getNumberOfComponents()==2 && getPntrToComponent(1)->forcesWereAdded() ) {
338 612 : unsigned natoms = getPntrToArgument(0)->getShape()[0] / 3;
339 612 : std::vector<Vector> pos( natoms ), der( natoms ), direction( natoms );
340 612 : double r = calculateRMSD( 0, pos, der, direction );
341 7472 : for(unsigned i=0; i<natoms; ++i) {
342 27440 : for(unsigned j=0; j<3; ++j ) {
343 20580 : forces[j*natoms+i] += direction[i][j];
344 : }
345 : }
346 : wasforced=true;
347 : }
348 0 : if( wasforced ) {
349 612 : unsigned ss=0;
350 612 : addForcesOnArguments( 0, forces, ss, getLabel() );
351 : }
352 : } else {
353 4368 : ActionWithVector::apply();
354 : }
355 : }
356 :
357 180588 : void RMSDVector::performTask( const unsigned& current, MultiValue& myvals ) const {
358 180588 : unsigned natoms = getPntrToArgument(0)->getShape()[0] / 3;
359 : std::vector<Vector>& pos( myvals.getFirstAtomVector() );
360 : std::vector<std::vector<Vector> > & allder( myvals.getFirstAtomDerivativeVector() );
361 180588 : if( allder.size()!=2 ) {
362 14048 : allder.resize(2);
363 : }
364 : std::vector<Vector>& der( allder[0] );
365 : std::vector<Vector>& direction( allder[1] );
366 180588 : if( pos.size()!=natoms ) {
367 14048 : pos.resize( natoms );
368 14048 : der.resize( natoms );
369 14048 : direction.resize( natoms );
370 : }
371 2528232 : for(unsigned i=0; i<pos.size(); ++i) {
372 9390576 : for(unsigned j=0; j<3; ++j) {
373 7042932 : pos[i][j] = getPntrToArgument(0)->get( j*natoms + i );
374 : }
375 : }
376 180588 : double r = calculateRMSD( current, pos, der, direction );
377 180588 : unsigned ostrn = getConstPntrToComponent(0)->getPositionInStream();
378 180588 : myvals.setValue( ostrn, r );
379 :
380 180588 : if( doNotCalculateDerivatives() ) {
381 : return;
382 : }
383 :
384 1929648 : for(unsigned i=0; i<natoms; i++) {
385 7167264 : for(unsigned j=0; j<3; ++j ) {
386 5375448 : myvals.addDerivative( ostrn, j*natoms+i, der[i][j] );
387 5375448 : myvals.updateIndex( ostrn, j*natoms+i );
388 : }
389 : }
390 : }
391 :
392 202902 : void RMSDVector::gatherStoredValue( const unsigned& valindex, const unsigned& code, const MultiValue& myvals,
393 : const unsigned& bufstart, std::vector<double>& buffer ) const {
394 202902 : if( getConstPntrToComponent(valindex)->getRank()==1 ) {
395 160146 : ActionWithVector::gatherStoredValue( valindex, code, myvals, bufstart, buffer );
396 160146 : return;
397 : }
398 : const std::vector<Vector>& direction( myvals.getConstFirstAtomDerivativeVector()[1] );
399 42756 : unsigned natoms = direction.size();
400 42756 : unsigned vindex = bufstart + 3*code*natoms;
401 598584 : for(unsigned i=0; i<natoms; ++i) {
402 2223312 : for(unsigned j=0; j<3; ++j ) {
403 1667484 : buffer[vindex + j*natoms + i] += direction[i][j];
404 : }
405 : }
406 : }
407 :
408 20442 : void RMSDVector::gatherForcesOnStoredValue( const Value* myval, const unsigned& itask, const MultiValue& myvals, std::vector<double>& forces ) const {
409 20442 : if( myval->getRank()==1 ) {
410 19350 : ActionWithVector::gatherForcesOnStoredValue( myval, itask, myvals, forces );
411 19350 : return;
412 : }
413 : const std::vector<Vector>& direction( myvals.getConstFirstAtomDerivativeVector()[1] );
414 1092 : unsigned natoms = direction.size();
415 15288 : for(unsigned i=0; i<natoms; ++i) {
416 56784 : for(unsigned j=0; j<3; ++j ) {
417 42588 : forces[j*natoms+i] += direction[i][j];
418 : }
419 : }
420 : }
421 :
422 :
423 : }
424 : }
425 :
426 :
427 :
|
