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/ActionWithValue.h"
23 : #include "core/ActionAtomistic.h"
24 : #include "core/ActionWithArguments.h"
25 : #include "reference/MetricRegister.h"
26 : #include "core/ActionRegister.h"
27 : #include "core/PlumedMain.h"
28 : #include "reference/Direction.h"
29 : #include "tools/Pbc.h"
30 :
31 : //+PLUMEDOC COLVAR PCAVARS
32 : /*
33 : Projection on principal component eigenvectors or other high dimensional linear subspace
34 :
35 : The collective variables described in \ref dists allow one to calculate the distance between the
36 : instantaneous structure adopted by the system and some high-dimensional, reference configuration. The
37 : problem with doing this is that, as one gets further and further from the reference configuration, the
38 : distance from it becomes a progressively poorer and poorer collective variable. This happens because
39 : the ``number" of structures at a distance \f$d\f$ from a reference configuration is proportional to \f$d^N\f$ in
40 : an \f$N\f$ dimensional space. Consequently, when \f$d\f$ is small the distance from the reference configuration
41 : may well be a good collective variable. However, when \f$d\f$ is large it is unlikely that the distance from the reference
42 : structure is a good CV. When the distance is large there will almost certainly be markedly different
43 : configuration that have the same CV value and hence barriers in transverse degrees of
44 : freedom.
45 :
46 : For these reasons dimensionality reduction is often employed so a projection \f$\mathbf{s}\f$ of a high-dimensional configuration
47 : \f$\mathbf{X}\f$ in a lower dimensionality space using a function:
48 :
49 : \f[
50 : \mathbf{s} = F(\mathbf{X}-\mathbf{X}^{ref})
51 : \f]
52 :
53 : where here we have introduced some high-dimensional reference configuration \f$\mathbf{X}^{ref}\f$. By far the simplest way to
54 : do this is to use some linear operator for \f$F\f$. That is to say we find a low-dimensional projection
55 : by rotating the basis vectors using some linear algebra:
56 :
57 : \f[
58 : \mathbf{s}_i = \sum_k A_{ik} ( X_{k} - X_{k}^{ref} )
59 : \f]
60 :
61 : Here \f$A\f$ is a \f$d\f$ by \f$D\f$ matrix where \f$D\f$ is the dimensionality of the high dimensional space and \f$d\f$ is
62 : the dimensionality of the lower dimensional subspace. In plumed when this kind of projection you can use the majority
63 : of the metrics detailed on \ref dists to calculate the displacement, \f$\mathbf{X}-\mathbf{X}^{ref}\f$, from the reference configuration.
64 : The matrix \f$A\f$ can be found by various means including principal component analysis and normal mode analysis. In both these methods the
65 : rows of \f$A\f$ would be the principle eigenvectors of a square matrix. For PCA the covariance while for normal modes the Hessian.
66 :
67 : \bug It is not possible to use the \ref DRMSD metric with this variable. You can get around this by listing the set of distances you wish to calculate for your DRMSD in the plumed file explicitly and using the EUCLIDEAN metric. MAHALONOBIS and NORM-EUCLIDEAN also do not work with this variable but using these options makes little sense when projecting on a linear subspace.
68 :
69 : \par Examples
70 :
71 : The following input calculates a projection on a linear subspace where the displacements
72 : from the reference configuration are calculated using the OPTIMAL metric. Consequently,
73 : both translation of the center of mass of the atoms and rotation of the reference
74 : frame are removed from these displacements. The matrix \f$A\f$ and the reference
75 : configuration \f$R^{ref}\f$ are specified in the pdb input file reference.pdb and the
76 : value of all projections (and the residual) are output to a file called colvar2.
77 :
78 : \plumedfile
79 : PCAVARS REFERENCE=reference.pdb TYPE=OPTIMAL LABEL=pca2
80 : PRINT ARG=pca2.* FILE=colvar2
81 : \endplumedfile
82 :
83 : The reference configurations can be specified using a pdb file. The first configuration that you provide is the reference configuration,
84 : which is referred to in the above as \f$X^{ref}\f$ subsequent configurations give the directions of row vectors that are contained in
85 : the matrix \f$A\f$ above. These directions can be specified by specifying a second configuration - in this case a vector will
86 : be constructed by calculating the displacement of this second configuration from the reference configuration. A pdb input prepared
87 : in this way would look as follows:
88 :
89 : \auxfile{reference.pdb}
90 : REMARK TYPE=OPTIMAL
91 : ATOM 2 CH3 ACE 1 12.932 -14.718 -6.016 1.00 1.00
92 : ATOM 5 C ACE 1 21.312 -9.928 -5.946 1.00 1.00
93 : ATOM 9 CA ALA 2 19.462 -11.088 -8.986 1.00 1.00
94 : ATOM 13 HB2 ALA 2 21.112 -10.688 -12.476 1.00 1.00
95 : ATOM 15 C ALA 2 19.422 7.978 -14.536 1.00 1.00
96 : ATOM 20 HH31 NME 3 20.122 -9.928 -17.746 1.00 1.00
97 : ATOM 21 HH32 NME 3 18.572 -13.148 -16.346 1.00 1.00
98 : END
99 : REMARK TYPE=OPTIMAL
100 : ATOM 2 CH3 ACE 1 13.932 -14.718 -6.016 1.00 1.00
101 : ATOM 5 C ACE 1 20.312 -9.928 -5.946 1.00 1.00
102 : ATOM 9 CA ALA 2 18.462 -11.088 -8.986 1.00 1.00
103 : ATOM 13 HB2 ALA 2 20.112 -11.688 -12.476 1.00 1.00
104 : ATOM 15 C ALA 2 19.422 7.978 -12.536 1.00 1.00
105 : ATOM 20 HH31 NME 3 20.122 -9.928 -17.746 1.00 1.00
106 : ATOM 21 HH32 NME 3 18.572 -13.148 -16.346 1.00 1.00
107 : END
108 : \endauxfile
109 :
110 : Alternatively, the second configuration can specify the components of \f$A\f$ explicitly. In this case you need to include the
111 : keyword TYPE=DIRECTION in the remarks to the pdb as shown below.
112 :
113 : \verbatim
114 : ATOM 2 CH3 ACE 1 12.932 -14.718 -6.016 1.00 1.00
115 : ATOM 5 C ACE 1 21.312 -9.928 -5.946 1.00 1.00
116 : ATOM 9 CA ALA 2 19.462 -11.088 -8.986 1.00 1.00
117 : ATOM 13 HB2 ALA 2 21.112 -10.688 -12.476 1.00 1.00
118 : ATOM 15 C ALA 2 19.422 7.978 -14.536 1.00 1.00
119 : ATOM 20 HH31 NME 3 20.122 -9.928 -17.746 1.00 1.00
120 : ATOM 21 HH32 NME 3 18.572 -13.148 -16.346 1.00 1.00
121 : END
122 : REMARK TYPE=DIRECTION
123 : ATOM 2 CH3 ACE 1 0.1414 0.3334 -0.0302 1.00 0.00
124 : ATOM 5 C ACE 1 0.0893 -0.1095 -0.1434 1.00 0.00
125 : ATOM 9 CA ALA 2 0.0207 -0.321 0.0321 1.00 0.00
126 : ATOM 13 HB2 ALA 2 0.0317 -0.6085 0.0783 1.00 0.00
127 : ATOM 15 C ALA 2 0.1282 -0.4792 0.0797 1.00 0.00
128 : ATOM 20 HH31 NME 3 0.0053 -0.465 0.0309 1.00 0.00
129 : ATOM 21 HH32 NME 3 -0.1019 -0.4261 -0.0082 1.00 0.00
130 : END
131 : \endverbatim
132 :
133 : If your metric involves arguments the labels of these arguments in your plumed input file should be specified in the REMARKS
134 : for each of the frames of your path. An input file in this case might look like this:
135 :
136 : \verbatim
137 : DESCRIPTION: a pca eigenvector specified using the start point and direction in the HD space.
138 : REMARK WEIGHT=1.0
139 : REMARK ARG=d1,d2
140 : REMARK d1=1.0 d2=1.0
141 : END
142 : REMARK TYPE=DIRECTION
143 : REMARK ARG=d1,d2
144 : REMARK d1=0.1 d2=0.25
145 : END
146 : \endverbatim
147 :
148 : Here we are working with the EUCLIDEAN metric and notice that we have specified the components of \f$A\f$ using DIRECTION.
149 : Consequently, the values of d1 and d2 in the second frame above do not specify a particular coordinate in the high-dimensional
150 : space as in they do in the first frame. Instead these values are the coefficients that can be used to construct a linear combination of d1 and d2.
151 : If we wanted to specify the direction in this metric using the start and end point of the vector we would write:
152 :
153 : \verbatim
154 : DESCRIPTION: a pca eigenvector specified using the start and end point of a vector in the HD space.
155 : REMARK WEIGHT=1.0
156 : REMARK ARG=d1,d2
157 : REMARK d1=1.0 d2=1.0
158 : END
159 : REMARK ARG=d1,d2
160 : REMARK d1=1.1 d2=1.25
161 : END
162 : \endverbatim
163 :
164 : */
165 : //+ENDPLUMEDOC
166 :
167 : namespace PLMD {
168 : namespace mapping {
169 :
170 : class PCAVars :
171 : public ActionWithValue,
172 : public ActionAtomistic,
173 : public ActionWithArguments {
174 : private:
175 : /// The holders for the derivatives
176 : MultiValue myvals;
177 : ReferenceValuePack mypack;
178 : /// The position of the reference configuration (the one we align to)
179 : std::unique_ptr<ReferenceConfiguration> myref;
180 : /// The eigenvectors we are interested in
181 : std::vector<Direction> directions;
182 : /// Stuff for applying forces
183 : std::vector<double> forces, forcesToApply;
184 : bool nopbc;
185 : public:
186 : static void registerKeywords( Keywords& keys );
187 : explicit PCAVars(const ActionOptions&);
188 : unsigned getNumberOfDerivatives() override;
189 : void lockRequests() override;
190 : void unlockRequests() override;
191 : void calculateNumericalDerivatives( ActionWithValue* a ) override;
192 : void calculate() override;
193 : void apply() override;
194 : };
195 :
196 13795 : PLUMED_REGISTER_ACTION(PCAVars,"PCAVARS")
197 :
198 9 : void PCAVars::registerKeywords( Keywords& keys ) {
199 9 : Action::registerKeywords( keys );
200 9 : ActionWithValue::registerKeywords( keys );
201 9 : ActionAtomistic::registerKeywords( keys );
202 9 : ActionWithArguments::registerKeywords( keys );
203 9 : componentsAreNotOptional(keys);
204 9 : keys.use("ARG");
205 18 : keys.addOutputComponent("eig","default","the projections on each eigenvalue are stored on values labeled eig-1, eig-2, ...");
206 18 : keys.addOutputComponent("residual","default","the distance of the configuration from the linear subspace defined "
207 : "by the vectors, eig-1, eig2, ... that are contained in the rows of A.");
208 18 : keys.add("compulsory","REFERENCE","a pdb file containing the reference configuration and configurations that define the directions for each eigenvector");
209 18 : keys.add("compulsory","TYPE","OPTIMAL","The method we are using for alignment to the reference structure");
210 18 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
211 9 : }
212 :
213 5 : PCAVars::PCAVars(const ActionOptions& ao):
214 : Action(ao),
215 : ActionWithValue(ao),
216 : ActionAtomistic(ao),
217 : ActionWithArguments(ao),
218 10 : myvals(1,0),
219 5 : mypack(0,0,myvals),
220 10 : nopbc(false) {
221 :
222 : // What type of distance are we calculating
223 : std::string mtype;
224 5 : parse("TYPE",mtype);
225 :
226 10 : parseFlag("NOPBC",nopbc);
227 :
228 : // Open reference file
229 : std::string reference;
230 10 : parse("REFERENCE",reference);
231 5 : FILE* fp=this->fopen(reference.c_str(),"r");
232 5 : if(!fp) {
233 0 : error("could not open reference file " + reference );
234 : }
235 :
236 : // call fclose when exiting this block
237 5 : auto deleter=[this](FILE* f) {
238 5 : this->fclose(f);
239 5 : };
240 : std::unique_ptr<FILE,decltype(deleter)> fp_deleter(fp,deleter);
241 :
242 : // Read all reference configurations
243 : // MultiReferenceBase myframes( "", false );
244 : std::vector<std::unique_ptr<ReferenceConfiguration> > myframes;
245 : bool do_read=true;
246 : unsigned nfram=0;
247 20 : while (do_read) {
248 20 : PDB mypdb;
249 : // Read the pdb file
250 40 : do_read=mypdb.readFromFilepointer(fp,plumed.getAtoms().usingNaturalUnits(),0.1/atoms.getUnits().getLength());
251 : // Fix argument names
252 20 : if(do_read) {
253 15 : if( nfram==0 ) {
254 10 : myref=metricRegister().create<ReferenceConfiguration>( mtype, mypdb );
255 5 : Direction* tdir = dynamic_cast<Direction*>( myref.get() );
256 5 : if( tdir ) {
257 0 : error("first frame should be reference configuration - not direction of vector");
258 : }
259 5 : if( !myref->pcaIsEnabledForThisReference() ) {
260 0 : error("can't do PCA with reference type " + mtype );
261 : }
262 : // std::vector<std::string> remarks( mypdb.getRemark() ); std::string rtype;
263 : // bool found=Tools::parse( remarks, "TYPE", rtype );
264 : // if(!found){ std::vector<std::string> newrem(1); newrem[0]="TYPE="+mtype; mypdb.addRemark(newrem); }
265 : // myframes.push_back( metricRegister().create<ReferenceConfiguration>( "", mypdb ) );
266 : } else {
267 10 : auto mymsd = metricRegister().create<ReferenceConfiguration>( "", mypdb );
268 10 : myframes.emplace_back( std::move(mymsd) );
269 10 : }
270 15 : nfram++;
271 : } else {
272 : break;
273 : }
274 20 : }
275 :
276 5 : if( nfram<=1 ) {
277 0 : error("no eigenvectors were specified");
278 : }
279 5 : log.printf(" found %u eigenvectors in file %s \n",nfram-1,reference.c_str() );
280 :
281 : // Finish the setup of the mapping object
282 : // Get the arguments and atoms that are required
283 : std::vector<AtomNumber> atoms;
284 5 : myref->getAtomRequests( atoms, false );
285 : std::vector<std::string> args;
286 5 : myref->getArgumentRequests( args, false );
287 5 : if( atoms.size()>0 ) {
288 5 : log.printf(" found %zu atoms in input \n",atoms.size());
289 5 : log.printf(" with indices : ");
290 40 : for(unsigned i=0; i<atoms.size(); ++i) {
291 35 : if(i%25==0) {
292 5 : log<<"\n";
293 : }
294 35 : log.printf("%d ",atoms[i].serial());
295 : }
296 5 : log.printf("\n");
297 : }
298 5 : requestAtoms( atoms );
299 : std::vector<Value*> req_args;
300 5 : interpretArgumentList( args, req_args );
301 5 : requestArguments( req_args );
302 :
303 : // And now check that the atoms/arguments are the same in all the eigenvalues
304 15 : for(unsigned i=0; i<myframes.size(); ++i) {
305 10 : myframes[i]->getAtomRequests( atoms, false );
306 10 : myframes[i]->getArgumentRequests( args, false );
307 : }
308 :
309 : // Setup the derivative pack
310 5 : if( atoms.size()>0 ) {
311 5 : myvals.resize( 1, args.size() + 3*atoms.size() + 9 );
312 : } else {
313 0 : myvals.resize( 1, args.size() );
314 : }
315 5 : mypack.resize( args.size(), atoms.size() );
316 40 : for(unsigned i=0; i<atoms.size(); ++i) {
317 : mypack.setAtomIndex( i, i );
318 : }
319 : /// This sets up all the storage data required by PCA in the pack
320 5 : myref->setupPCAStorage( mypack );
321 :
322 : // Check there are no periodic arguments
323 5 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
324 0 : if( getPntrToArgument(i)->isPeriodic() ) {
325 0 : error("cannot use periodic variables in pca projections");
326 : }
327 : }
328 : // Work out if the user wants to normalise the input vector
329 5 : checkRead();
330 :
331 5 : if(nopbc) {
332 1 : log.printf(" without periodic boundary conditions\n");
333 : } else {
334 4 : log.printf(" using periodic boundary conditions\n");
335 : }
336 :
337 : // Resize the matrices that will hold our eivenvectors
338 5 : PDB mypdb;
339 5 : mypdb.setAtomNumbers( atoms );
340 5 : mypdb.addBlockEnd( atoms.size() );
341 5 : if( args.size()>0 ) {
342 0 : mypdb.setArgumentNames( args );
343 : }
344 : // Resize the matrices that will hold our eivenvectors
345 15 : for(unsigned i=0; i<myframes.size(); ++i) {
346 20 : directions.push_back( Direction(ReferenceConfigurationOptions("DIRECTION")));
347 10 : directions[i].read( mypdb );
348 : }
349 :
350 : // Create fake periodic boundary condition (these would only be used for DRMSD which is not allowed)
351 : // Now calculate the eigenvectors
352 15 : for(unsigned i=0; i<myframes.size(); ++i) {
353 : // Calculate distance from reference configuration
354 10 : myframes[i]->extractDisplacementVector( myref->getReferencePositions(), getArguments(), myref->getReferenceArguments(), true, directions[i] );
355 : // Create a component to store the output
356 : std::string num;
357 10 : Tools::convert( i+1, num );
358 10 : addComponentWithDerivatives("eig-"+num);
359 20 : componentIsNotPeriodic("eig-"+num);
360 : }
361 5 : addComponentWithDerivatives("residual");
362 10 : componentIsNotPeriodic("residual");
363 :
364 : // Get appropriate number of derivatives
365 : unsigned nder;
366 5 : if( getNumberOfAtoms()>0 ) {
367 5 : nder = 3*getNumberOfAtoms() + 9 + getNumberOfArguments();
368 : } else {
369 : nder = getNumberOfArguments();
370 : }
371 :
372 : // Resize all derivative arrays
373 5 : forces.resize( nder );
374 5 : forcesToApply.resize( nder );
375 20 : for(int i=0; i<getNumberOfComponents(); ++i) {
376 15 : getPntrToComponent(i)->resizeDerivatives(nder);
377 : }
378 20 : }
379 :
380 36 : unsigned PCAVars::getNumberOfDerivatives() {
381 36 : if( getNumberOfAtoms()>0 ) {
382 36 : return 3*getNumberOfAtoms() + 9 + getNumberOfArguments();
383 : }
384 0 : return getNumberOfArguments();
385 : }
386 :
387 55 : void PCAVars::lockRequests() {
388 : ActionWithArguments::lockRequests();
389 : ActionAtomistic::lockRequests();
390 55 : }
391 :
392 55 : void PCAVars::unlockRequests() {
393 : ActionWithArguments::unlockRequests();
394 : ActionAtomistic::unlockRequests();
395 55 : }
396 :
397 55 : void PCAVars::calculate() {
398 :
399 55 : if(!nopbc && getNumberOfAtoms()>0) {
400 44 : makeWhole();
401 : }
402 :
403 : // Clear the reference value pack
404 55 : mypack.clear();
405 : // Calculate distance between instaneous configuration and reference
406 55 : double dist = myref->calculate( getPositions(), getPbc(), getArguments(), mypack, true );
407 :
408 : // Start accumulating residual by adding derivatives of distance
409 55 : Value* resid=getPntrToComponent( getNumberOfComponents()-1 );
410 55 : unsigned nargs=getNumberOfArguments();
411 55 : for(unsigned j=0; j<getNumberOfArguments(); ++j) {
412 : resid->addDerivative( j, mypack.getArgumentDerivative(j) );
413 : }
414 440 : for(unsigned j=0; j<getNumberOfAtoms(); ++j) {
415 385 : Vector ader=mypack.getAtomDerivative( j );
416 1540 : for(unsigned k=0; k<3; ++k) {
417 1155 : resid->addDerivative( nargs +3*j+k, ader[k] );
418 : }
419 : }
420 : // Retrieve the values of all arguments
421 55 : std::vector<double> args( getNumberOfArguments() );
422 55 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
423 0 : args[i]=getArgument(i);
424 : }
425 :
426 : // Now calculate projections on pca vectors
427 55 : Vector adif, ader;
428 55 : Tensor fvir, tvir;
429 165 : for(int i=0; i<getNumberOfComponents()-1; ++i) { // One less component as we also have residual
430 110 : double proj=myref->projectDisplacementOnVector( directions[i], getArguments(), args, mypack );
431 : // And now accumulate derivatives
432 110 : Value* eid=getPntrToComponent(i);
433 110 : for(unsigned j=0; j<getNumberOfArguments(); ++j) {
434 : eid->addDerivative( j, mypack.getArgumentDerivative(j) );
435 : }
436 110 : if( getNumberOfAtoms()>0 ) {
437 110 : tvir.zero();
438 880 : for(unsigned j=0; j<getNumberOfAtoms(); ++j) {
439 770 : Vector myader=mypack.getAtomDerivative(j);
440 3080 : for(unsigned k=0; k<3; ++k) {
441 2310 : eid->addDerivative( nargs + 3*j+k, myader[k] );
442 2310 : resid->addDerivative( nargs + 3*j+k, -2*proj*myader[k] );
443 : }
444 770 : tvir += -1.0*Tensor( getPosition(j), myader );
445 : }
446 440 : for(unsigned j=0; j<3; ++j) {
447 1320 : for(unsigned k=0; k<3; ++k) {
448 990 : eid->addDerivative( nargs + 3*getNumberOfAtoms() + 3*j + k, tvir(j,k) );
449 : }
450 : }
451 : }
452 110 : dist -= proj*proj; // Subtract square from total squared distance to get residual squared
453 : // Derivatives of residual
454 110 : for(unsigned j=0; j<getNumberOfArguments(); ++j) {
455 0 : resid->addDerivative( j, -2*proj*eid->getDerivative(j) );
456 : }
457 : // for(unsigned j=0;j<getNumberOfArguments();++j) resid->addDerivative( j, -2*proj*arg_eigv(i,j) );
458 : // And set final value
459 110 : getPntrToComponent(i)->set( proj );
460 : }
461 55 : dist=std::sqrt(dist);
462 : resid->set( dist );
463 :
464 : // Take square root of residual derivatives
465 55 : double prefactor = 0.5 / dist;
466 55 : for(unsigned j=0; j<getNumberOfArguments(); ++j) {
467 0 : resid->setDerivative( j, prefactor*resid->getDerivative(j) );
468 : }
469 440 : for(unsigned j=0; j<getNumberOfAtoms(); ++j) {
470 1540 : for(unsigned k=0; k<3; ++k) {
471 1155 : resid->setDerivative( nargs + 3*j+k, prefactor*resid->getDerivative( nargs+3*j+k ) );
472 : }
473 : }
474 :
475 : // And finally virial for residual
476 55 : if( getNumberOfAtoms()>0 ) {
477 55 : tvir.zero();
478 440 : for(unsigned j=0; j<getNumberOfAtoms(); ++j) {
479 385 : Vector ader;
480 1540 : for(unsigned k=0; k<3; ++k) {
481 1155 : ader[k]=resid->getDerivative( nargs + 3*j+k );
482 : }
483 385 : tvir += -1.0*Tensor( getPosition(j), ader );
484 : }
485 220 : for(unsigned j=0; j<3; ++j) {
486 660 : for(unsigned k=0; k<3; ++k) {
487 495 : resid->addDerivative( nargs + 3*getNumberOfAtoms() + 3*j + k, tvir(j,k) );
488 : }
489 : }
490 : }
491 :
492 55 : }
493 :
494 0 : void PCAVars::calculateNumericalDerivatives( ActionWithValue* a ) {
495 0 : if( getNumberOfArguments()>0 ) {
496 0 : ActionWithArguments::calculateNumericalDerivatives( a );
497 : }
498 0 : if( getNumberOfAtoms()>0 ) {
499 0 : Matrix<double> save_derivatives( getNumberOfComponents(), getNumberOfArguments() );
500 0 : for(int j=0; j<getNumberOfComponents(); ++j) {
501 0 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
502 0 : save_derivatives(j,i)=getPntrToComponent(j)->getDerivative(i);
503 : }
504 : }
505 0 : calculateAtomicNumericalDerivatives( a, getNumberOfArguments() );
506 0 : for(int j=0; j<getNumberOfComponents(); ++j) {
507 0 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
508 0 : getPntrToComponent(j)->addDerivative( i, save_derivatives(j,i) );
509 : }
510 : }
511 : }
512 0 : }
513 :
514 55 : void PCAVars::apply() {
515 :
516 : bool wasforced=false;
517 55 : forcesToApply.assign(forcesToApply.size(),0.0);
518 220 : for(int i=0; i<getNumberOfComponents(); ++i) {
519 165 : if( getPntrToComponent(i)->applyForce( forces ) ) {
520 : wasforced=true;
521 0 : for(unsigned i=0; i<forces.size(); ++i) {
522 0 : forcesToApply[i]+=forces[i];
523 : }
524 : }
525 : }
526 55 : if( wasforced ) {
527 0 : addForcesOnArguments( forcesToApply );
528 0 : if( getNumberOfAtoms()>0 ) {
529 0 : setForcesOnAtoms( forcesToApply, getNumberOfArguments() );
530 : }
531 : }
532 :
533 55 : }
534 :
535 : }
536 : }
|
