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); keys.setDisplayName("RMSD");
44 162 : keys.addInputKeyword("compulsory","ARG","vector/matrix","the labels of two actions that you are calculating the RMSD between");
45 162 : keys.add("compulsory","TYPE","SIMPLE","the manner in which RMSD alignment is performed. Should be OPTIMAL or SIMPLE.");
46 162 : keys.add("compulsory","ALIGN","1.0","the weights to use when aligning to the reference structure");
47 162 : keys.add("compulsory","DISPLACE","1.0","the weights to use when calculating the displacement from the reference structure");
48 162 : keys.addFlag("SQUARED",false," This should be set if you want mean squared displacement instead of RMSD ");
49 162 : 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");
50 162 : keys.addFlag("DISPLACEMENT",false,"Calculate the vector of displacements instead of the length of this vector");
51 162 : keys.addOutputComponent("disp","DISPLACEMENT","vector/matrix","the vector of displacements for the atoms");
52 162 : keys.addOutputComponent("dist","DISPLACEMENT","scalar/vector","the RMSD distance the atoms have moved");
53 162 : keys.setValueDescription("scalar/vector","a vector containing the RMSD between the instantaneous structure and each of the reference structures that were input");
54 81 : }
55 :
56 49 : RMSDVector::RMSDVector(const ActionOptions&ao):
57 : Action(ao),
58 : ActionWithVector(ao),
59 49 : firststep(true)
60 : {
61 49 : if( getPntrToArgument(0)->getRank()!=1 ) error("first argument should be vector");
62 49 : if( getPntrToArgument(1)->getRank()<1 ) error("second argument should be matrix or a vector");
63 49 : if( getPntrToArgument(1)->getRank()==1 ) {
64 7 : if( getPntrToArgument(0)->getNumberOfValues()!=getPntrToArgument(1)->getNumberOfValues() ) error("mismatch between sizes of input vectors");
65 42 : } else if( getPntrToArgument(1)->getRank()==2 ) {
66 42 : if( getPntrToArgument(0)->getNumberOfValues()!=getPntrToArgument(1)->getShape()[1] ) error("mismatch between sizes of input vectors");
67 : }
68 49 : if( getPntrToArgument(0)->getNumberOfValues()%3!=0 ) error("number of components in input arguments should be multiple of three");
69 :
70 49 : unsigned natoms = getPntrToArgument(0)->getNumberOfValues() / 3;
71 98 : type.assign("SIMPLE"); parse("TYPE",type); parseFlag("SQUARED",squared);
72 49 : align.resize( natoms ); parseVector("ALIGN",align);
73 49 : displace.resize( natoms ); parseVector("DISPLACE",displace);
74 98 : bool unorm=false; parseFlag("UNORMALIZED",unorm); norm_weights=!unorm;
75 49 : double wa=0, wd=0; sqrtdisplace.resize( displace.size() );
76 630 : for(unsigned i=0; i<align.size(); ++i) { wa+=align[i]; wd+=displace[i]; }
77 :
78 49 : if( wa>epsilon ) {
79 49 : double iwa = 1. / wa;
80 630 : for(unsigned i=0; i<align.size(); ++i) align[i] *= iwa;
81 : } else {
82 0 : double iwa = 1. / natoms;
83 0 : for(unsigned i=0; i<align.size(); ++i) align[i] = iwa;
84 : }
85 49 : if( wd>epsilon ) {
86 49 : if( !norm_weights ) { wd = 1; } double iwd = 1. / wd;
87 630 : for(unsigned i=0; i<align.size(); ++i) displace[i] *= iwd;
88 : } else {
89 0 : double iwd = 1. / natoms;
90 0 : for(unsigned i=0; i<align.size(); ++i) displace[i] = iwd;
91 : }
92 630 : for(unsigned i=0; i<align.size(); ++i) sqrtdisplace[i] = sqrt(displace[i]);
93 :
94 49 : parseFlag("DISPLACEMENT",displacement);
95 49 : if( displacement && (getPntrToArgument(1)->getRank()==1 || getPntrToArgument(1)->getShape()[0]<=1) ) {
96 78 : addComponentWithDerivatives("dist"); componentIsNotPeriodic("dist");
97 26 : std::vector<unsigned> shape( 1, getPntrToArgument(0)->getNumberOfValues() );
98 78 : addComponent( "disp", shape ); getPntrToComponent(1)->buildDataStore(); componentIsNotPeriodic("disp");
99 23 : } else if( displacement ) {
100 2 : std::vector<unsigned> shape( 1, getPntrToArgument(1)->getShape()[0] );
101 4 : addComponent( "dist", shape ); getPntrToComponent(0)->buildDataStore();
102 2 : componentIsNotPeriodic("dist");
103 2 : shape.resize(2); shape[0] = getPntrToArgument(1)->getShape()[0]; shape[1] = getPntrToArgument(0)->getNumberOfValues();
104 4 : addComponent( "disp", shape ); getPntrToComponent(1)->buildDataStore(); getPntrToComponent(1)->reshapeMatrixStore( shape[1] );
105 4 : componentIsNotPeriodic("disp");
106 21 : } else if( (getPntrToArgument(1)->getRank()==1 || getPntrToArgument(1)->getShape()[0]==1) ) {
107 26 : addValue(); setNotPeriodic();
108 : } else {
109 8 : std::vector<unsigned> shape( 1, getPntrToArgument(1)->getShape()[0] );
110 8 : addValue( shape ); setNotPeriodic();
111 : }
112 49 : if( getPntrToArgument(1)->getRank()==1 || getPntrToArgument(1)->getNumberOfValues()==0 ) myrmsd.resize(1);
113 42 : else myrmsd.resize( getPntrToArgument(1)->getShape()[0] );
114 :
115 49 : if( getPntrToArgument(1)->getRank()==1 ) log.printf(" calculating RMSD distance between %d atoms. Distance between the avectors of atoms in %s and %s\n",
116 : natoms, getPntrToArgument(0)->getName().c_str(), getPntrToArgument(1)->getName().c_str() );
117 42 : else log.printf(" calculating RMSD distance between %d sets of %d atoms. Distance between vector %s of atoms and matrix of configurations in %s\n",
118 : getPntrToArgument(1)->getShape()[0], natoms, getPntrToArgument(0)->getName().c_str(), getPntrToArgument(1)->getName().c_str() );
119 49 : log.printf(" method for alignment : %s \n",type.c_str() );
120 49 : if(squared)log.printf(" chosen to use SQUARED option for MSD instead of RMSD\n");
121 21 : else log.printf(" using periodic boundary conditions\n");
122 49 : }
123 :
124 804 : unsigned RMSDVector::getNumberOfDerivatives() {
125 804 : return getPntrToArgument(0)->getNumberOfValues() + getPntrToArgument(1)->getNumberOfValues();
126 : }
127 :
128 10631 : void RMSDVector::setReferenceConfigurations() {
129 10631 : unsigned natoms = getPntrToArgument(0)->getShape()[0] / 3;
130 10631 : Vector center; std::vector<Vector> pos( natoms );
131 22578 : for(unsigned jconf=0; jconf<myrmsd.size(); ++jconf) {
132 11947 : center.zero();
133 172431 : for(unsigned i=0; i<pos.size(); ++i) {
134 641936 : for(unsigned j=0; j<3; ++j) pos[i][j] = getPntrToArgument(1)->get( (3*jconf+j)*pos.size() + i );
135 160484 : center+=pos[i]*align[i];
136 : }
137 172431 : for(unsigned i=0; i<pos.size(); ++i) pos[i] -= center;
138 11947 : myrmsd[jconf].clear(); myrmsd[jconf].set(align,displace,pos,type,true,norm_weights);
139 : }
140 10631 : }
141 :
142 192996 : double RMSDVector::calculateRMSD( const unsigned& current, std::vector<Vector>& pos, std::vector<Vector>& der, std::vector<Vector>& direction ) const {
143 192996 : unsigned natoms = pos.size();
144 2704147 : for(unsigned i=0; i<natoms; ++i) {
145 10044604 : for(unsigned j=0; j<3; ++j) pos[i][j] = getPntrToArgument(0)->get( j*natoms + i );
146 : }
147 :
148 192996 : if( displacement && type=="SIMPLE" ) {
149 646 : const Value* myval = getConstPntrToComponent(1);
150 646 : double r = myrmsd[current].simpleAlignment( align, displace, pos, myrmsd[current].getReference(), der, direction, squared );
151 646 : if( !doNotCalculateDerivatives() && myval->forcesWereAdded() ) {
152 33 : Vector comforce; comforce.zero();
153 187 : for(unsigned i=0; i<natoms; i++) {
154 616 : for(unsigned k=0; k<3; ++k) comforce[k] += align[i]*myval->getForce( (3*current+k)*natoms + i);
155 : }
156 187 : for(unsigned i=0; i<natoms; i++) {
157 616 : for(unsigned k=0; k<3; ++k) direction[i][k] = myval->getForce( (3*current+k)*natoms + i ) - comforce[k];
158 : }
159 : }
160 646 : return r;
161 192350 : } else if( displacement ) {
162 55589 : const Value* myval = getConstPntrToComponent(1);
163 111178 : Tensor rot; Matrix<std::vector<Vector> > DRotDPos(3,3); std::vector<Vector> centeredpos( natoms ), centeredreference( natoms );
164 55589 : double r = myrmsd[current].calc_PCAelements( pos, der, rot, DRotDPos, direction, centeredpos, centeredreference, squared ); std::vector<Vector> ref( myrmsd[current].getReference() );
165 55589 : if( !doNotCalculateDerivatives() && myval->forcesWereAdded() ) {
166 1671 : Tensor trot=rot.transpose(); double prefactor = 1 / static_cast<double>( natoms ); Vector v1; v1.zero();
167 22573 : for(unsigned n=0; n<natoms; n++) {
168 83608 : Vector ff; for(unsigned k=0; k<3; ++k ) ff[k] = myval->getForce( (3*current+k)*natoms + n );
169 20902 : v1+=prefactor*matmul(trot,ff);
170 : }
171 : // Notice that we use centreredreference here to accumulate the true forces
172 22573 : for(unsigned n=0; n<natoms; n++) {
173 83608 : Vector ff; for(unsigned k=0; k<3; ++k ) ff[k] = myval->getForce( (3*current+k)*natoms + n );
174 20902 : centeredreference[n] = sqrtdisplace[n]*( matmul(trot,ff) - v1 );
175 : }
176 6684 : for(unsigned a=0; a<3; a++) {
177 20052 : for(unsigned b=0; b<3; b++) {
178 203157 : double tmp1=0.; for(unsigned m=0; m<natoms; m++) tmp1+=centeredpos[m][b]*myval->getForce( (3*current+a)*natoms + m );
179 203157 : for(unsigned i=0; i<natoms; i++) centeredreference[i] += sqrtdisplace[i]*tmp1*DRotDPos[a][b][i];
180 : }
181 : }
182 : // Now subtract the current force and add on the true force
183 22573 : for(unsigned n=0; n<natoms; n++) {
184 83608 : for(unsigned k=0; k<3; ++k) direction[n][k] = centeredreference[n][k];
185 : }
186 : } else {
187 759226 : for(unsigned i=0; i<direction.size(); ++i) direction[i] = sqrtdisplace[i]*( direction[i] - ref[i] );
188 : }
189 : return r;
190 : }
191 136761 : return myrmsd[current].calculate( pos, der, squared );
192 : }
193 :
194 : // calculator
195 17182 : void RMSDVector::calculate() {
196 17182 : if( firststep || !getPntrToArgument(1)->isConstant() ) { setReferenceConfigurations(); firststep=false; }
197 :
198 17182 : if( getPntrToComponent(0)->getRank()==0 ) {
199 11796 : unsigned natoms = getPntrToArgument(0)->getShape()[0] / 3;
200 11796 : std::vector<Vector> pos( natoms ), der( natoms ), direction( natoms );
201 11796 : double r = calculateRMSD( 0, pos, der, direction );
202 :
203 11796 : getPntrToComponent(0)->set( r );
204 11796 : if( getNumberOfComponents()==2 ) {
205 11775 : Value* mydisp = getPntrToComponent(1);
206 168204 : for(unsigned i=0; i<natoms; i++) {
207 625716 : for(unsigned j=0; j<3; ++j ) mydisp->set( j*natoms+i, direction[i][j] );
208 : }
209 : }
210 11796 : if( doNotCalculateDerivatives() ) return;
211 :
212 612 : Value* myval = getPntrToComponent(0);
213 7472 : for(unsigned i=0; i<natoms; i++) {
214 27440 : for(unsigned j=0; j<3; ++j ) myval->setDerivative( j*natoms+i, der[i][j] );
215 : }
216 5386 : } else runAllTasks();
217 : }
218 :
219 160524 : bool RMSDVector::checkForTaskForce( const unsigned& itask, const Value* myval ) const {
220 160524 : if( myval->getRank()<2 ) return ActionWithVector::checkForTaskForce( itask, myval );
221 22932 : unsigned nelements = myval->getShape()[1], startr = itask*nelements;
222 874692 : for(unsigned j=0; j<nelements; ++j ) {
223 852852 : if( fabs( myval->getForce( startr + j ) )>epsilon ) return true;
224 : }
225 : return false;
226 : }
227 :
228 17162 : void RMSDVector::apply() {
229 17162 : if( doNotCalculateDerivatives() ) return;
230 :
231 4980 : if( getPntrToComponent(0)->getRank()==0 ) {
232 612 : std::vector<double> forces( getNumberOfDerivatives(), 0 );
233 612 : bool wasforced = getPntrToComponent(0)->applyForce( forces );
234 :
235 612 : if( getNumberOfComponents()==2 && getPntrToComponent(1)->forcesWereAdded() ) {
236 612 : unsigned natoms = getPntrToArgument(0)->getShape()[0] / 3;
237 612 : std::vector<Vector> pos( natoms ), der( natoms ), direction( natoms );
238 612 : double r = calculateRMSD( 0, pos, der, direction );
239 7472 : for(unsigned i=0; i<natoms; ++i) {
240 27440 : for(unsigned j=0; j<3; ++j ) forces[j*natoms+i] += direction[i][j];
241 : }
242 : wasforced=true;
243 : }
244 612 : if( wasforced ) { unsigned ss=0; addForcesOnArguments( 0, forces, ss, getLabel() ); }
245 4368 : } else ActionWithVector::apply();
246 : }
247 :
248 180588 : void RMSDVector::performTask( const unsigned& current, MultiValue& myvals ) const {
249 180588 : unsigned natoms = getPntrToArgument(0)->getShape()[0] / 3;
250 : std::vector<Vector>& pos( myvals.getFirstAtomVector() );
251 : std::vector<std::vector<Vector> > & allder( myvals.getFirstAtomDerivativeVector() );
252 180588 : if( allder.size()!=2 ) allder.resize(2);
253 : std::vector<Vector>& der( allder[0] ); std::vector<Vector>& direction( allder[1] );
254 180588 : if( pos.size()!=natoms ) { pos.resize( natoms ); der.resize( natoms ); direction.resize( natoms ); }
255 2528232 : for(unsigned i=0; i<pos.size(); ++i) {
256 9390576 : for(unsigned j=0; j<3; ++j) pos[i][j] = getPntrToArgument(0)->get( j*natoms + i );
257 : }
258 180588 : double r = calculateRMSD( current, pos, der, direction );
259 180588 : unsigned ostrn = getConstPntrToComponent(0)->getPositionInStream();
260 180588 : myvals.setValue( ostrn, r );
261 :
262 180588 : if( doNotCalculateDerivatives() ) return;
263 :
264 1929648 : for(unsigned i=0; i<natoms; i++) {
265 7167264 : for(unsigned j=0; j<3; ++j ) { myvals.addDerivative( ostrn, j*natoms+i, der[i][j] ); myvals.updateIndex( ostrn, j*natoms+i ); }
266 : }
267 : }
268 :
269 202902 : void RMSDVector::gatherStoredValue( const unsigned& valindex, const unsigned& code, const MultiValue& myvals,
270 : const unsigned& bufstart, std::vector<double>& buffer ) const {
271 202902 : if( getConstPntrToComponent(valindex)->getRank()==1 ) { ActionWithVector::gatherStoredValue( valindex, code, myvals, bufstart, buffer ); return; }
272 : const std::vector<Vector>& direction( myvals.getConstFirstAtomDerivativeVector()[1] );
273 42756 : unsigned natoms = direction.size(); unsigned vindex = bufstart + 3*code*natoms;
274 598584 : for(unsigned i=0; i<natoms; ++i) {
275 2223312 : for(unsigned j=0; j<3; ++j ) buffer[vindex + j*natoms + i] += direction[i][j];
276 : }
277 : }
278 :
279 20442 : void RMSDVector::gatherForcesOnStoredValue( const Value* myval, const unsigned& itask, const MultiValue& myvals, std::vector<double>& forces ) const {
280 20442 : if( myval->getRank()==1 ) { ActionWithVector::gatherForcesOnStoredValue( myval, itask, myvals, forces ); return; }
281 1092 : const std::vector<Vector>& direction( myvals.getConstFirstAtomDerivativeVector()[1] ); unsigned natoms = direction.size();
282 15288 : for(unsigned i=0; i<natoms; ++i) {
283 56784 : for(unsigned j=0; j<3; ++j ) forces[j*natoms+i] += direction[i][j];
284 : }
285 : }
286 :
287 :
288 : }
289 : }
290 :
291 :
292 :
|
