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 : {
175 : private:
176 : /// The holders for the derivatives
177 : MultiValue myvals;
178 : ReferenceValuePack mypack;
179 : /// The position of the reference configuration (the one we align to)
180 : std::unique_ptr<ReferenceConfiguration> myref;
181 : /// The eigenvectors we are interested in
182 : std::vector<Direction> directions;
183 : /// Stuff for applying forces
184 : std::vector<double> forces, forcesToApply;
185 : bool nopbc;
186 : public:
187 : static void registerKeywords( Keywords& keys );
188 : explicit PCAVars(const ActionOptions&);
189 : unsigned getNumberOfDerivatives() override;
190 : void lockRequests() override;
191 : void unlockRequests() override;
192 : void calculateNumericalDerivatives( ActionWithValue* a ) override;
193 : void calculate() override;
194 : void apply() override;
195 : };
196 :
197 10429 : PLUMED_REGISTER_ACTION(PCAVars,"PCAVARS")
198 :
199 6 : void PCAVars::registerKeywords( Keywords& keys ) {
200 6 : Action::registerKeywords( keys );
201 6 : ActionWithValue::registerKeywords( keys );
202 6 : ActionAtomistic::registerKeywords( keys );
203 6 : ActionWithArguments::registerKeywords( keys );
204 6 : componentsAreNotOptional(keys); keys.use("ARG");
205 12 : keys.addOutputComponent("eig","default","the projections on each eigenvalue are stored on values labeled eig-1, eig-2, ...");
206 12 : keys.addOutputComponent("residual","default","the distance of the configuration from the linear subspace defined "
207 : "by the vectors, \\f$e_i\\f$, that are contained in the rows of \\f$A\\f$. In other words this is "
208 : "\\f$\\sqrt( r^2 - \\sum_i [\\mathbf{r}.\\mathbf{e_i}]^2)\\f$ where "
209 : "\\f$r\\f$ is the distance between the instantaneous position and the "
210 : "reference point.");
211 12 : keys.add("compulsory","REFERENCE","a pdb file containing the reference configuration and configurations that define the directions for each eigenvector");
212 12 : keys.add("compulsory","TYPE","OPTIMAL","The method we are using for alignment to the reference structure");
213 12 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
214 6 : }
215 :
216 5 : PCAVars::PCAVars(const ActionOptions& ao):
217 : Action(ao),
218 : ActionWithValue(ao),
219 : ActionAtomistic(ao),
220 : ActionWithArguments(ao),
221 10 : myvals(1,0),
222 5 : mypack(0,0,myvals),
223 10 : nopbc(false)
224 : {
225 :
226 : // What type of distance are we calculating
227 5 : std::string mtype; parse("TYPE",mtype);
228 :
229 10 : parseFlag("NOPBC",nopbc);
230 :
231 : // Open reference file
232 10 : std::string reference; parse("REFERENCE",reference);
233 5 : FILE* fp=fopen(reference.c_str(),"r");
234 5 : if(!fp) error("could not open reference file " + reference );
235 :
236 : // Read all reference configurations
237 : // MultiReferenceBase myframes( "", false );
238 : std::vector<std::unique_ptr<ReferenceConfiguration> > myframes;
239 : bool do_read=true; unsigned nfram=0;
240 20 : while (do_read) {
241 20 : PDB mypdb;
242 : // Read the pdb file
243 40 : do_read=mypdb.readFromFilepointer(fp,plumed.getAtoms().usingNaturalUnits(),0.1/atoms.getUnits().getLength());
244 : // Fix argument names
245 20 : if(do_read) {
246 15 : if( nfram==0 ) {
247 10 : myref=metricRegister().create<ReferenceConfiguration>( mtype, mypdb );
248 5 : Direction* tdir = dynamic_cast<Direction*>( myref.get() );
249 5 : if( tdir ) error("first frame should be reference configuration - not direction of vector");
250 5 : if( !myref->pcaIsEnabledForThisReference() ) error("can't do PCA with reference type " + mtype );
251 : // std::vector<std::string> remarks( mypdb.getRemark() ); std::string rtype;
252 : // bool found=Tools::parse( remarks, "TYPE", rtype );
253 : // if(!found){ std::vector<std::string> newrem(1); newrem[0]="TYPE="+mtype; mypdb.addRemark(newrem); }
254 : // myframes.push_back( metricRegister().create<ReferenceConfiguration>( "", mypdb ) );
255 : } else {
256 10 : auto mymsd = metricRegister().create<ReferenceConfiguration>( "", mypdb );
257 10 : myframes.emplace_back( std::move(mymsd) );
258 10 : }
259 15 : nfram++;
260 : } else {
261 : break;
262 : }
263 20 : }
264 5 : fclose(fp);
265 :
266 5 : if( nfram<=1 ) error("no eigenvectors were specified");
267 5 : log.printf(" found %u eigenvectors in file %s \n",nfram-1,reference.c_str() );
268 :
269 : // Finish the setup of the mapping object
270 : // Get the arguments and atoms that are required
271 5 : std::vector<AtomNumber> atoms; myref->getAtomRequests( atoms, false );
272 5 : std::vector<std::string> args; myref->getArgumentRequests( args, false );
273 5 : if( atoms.size()>0 ) {
274 5 : log.printf(" found %zu atoms in input \n",atoms.size());
275 5 : log.printf(" with indices : ");
276 40 : for(unsigned i=0; i<atoms.size(); ++i) {
277 35 : if(i%25==0) log<<"\n";
278 35 : log.printf("%d ",atoms[i].serial());
279 : }
280 5 : log.printf("\n");
281 : }
282 5 : requestAtoms( atoms ); std::vector<Value*> req_args;
283 5 : interpretArgumentList( args, req_args ); requestArguments( req_args );
284 :
285 : // And now check that the atoms/arguments are the same in all the eigenvalues
286 15 : for(unsigned i=0; i<myframes.size(); ++i) { myframes[i]->getAtomRequests( atoms, false ); myframes[i]->getArgumentRequests( args, false ); }
287 :
288 : // Setup the derivative pack
289 5 : if( atoms.size()>0 ) myvals.resize( 1, args.size() + 3*atoms.size() + 9 );
290 0 : else myvals.resize( 1, args.size() );
291 5 : mypack.resize( args.size(), atoms.size() );
292 40 : for(unsigned i=0; i<atoms.size(); ++i) mypack.setAtomIndex( i, i );
293 : /// This sets up all the storage data required by PCA in the pack
294 5 : myref->setupPCAStorage( mypack );
295 :
296 : // Check there are no periodic arguments
297 5 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
298 0 : if( getPntrToArgument(i)->isPeriodic() ) error("cannot use periodic variables in pca projections");
299 : }
300 : // Work out if the user wants to normalise the input vector
301 5 : checkRead();
302 :
303 5 : if(nopbc) log.printf(" without periodic boundary conditions\n");
304 4 : else log.printf(" using periodic boundary conditions\n");
305 :
306 : // Resize the matrices that will hold our eivenvectors
307 5 : PDB mypdb; mypdb.setAtomNumbers( atoms ); mypdb.addBlockEnd( atoms.size() );
308 5 : if( args.size()>0 ) mypdb.setArgumentNames( args );
309 : // Resize the matrices that will hold our eivenvectors
310 15 : for(unsigned i=0; i<myframes.size(); ++i) {
311 20 : directions.push_back( Direction(ReferenceConfigurationOptions("DIRECTION"))); directions[i].read( mypdb );
312 : }
313 :
314 : // Create fake periodic boundary condition (these would only be used for DRMSD which is not allowed)
315 : // Now calculate the eigenvectors
316 15 : for(unsigned i=0; i<myframes.size(); ++i) {
317 : // Calculate distance from reference configuration
318 10 : myframes[i]->extractDisplacementVector( myref->getReferencePositions(), getArguments(), myref->getReferenceArguments(), true, directions[i] );
319 : // Create a component to store the output
320 10 : std::string num; Tools::convert( i+1, num );
321 30 : addComponentWithDerivatives("eig-"+num); componentIsNotPeriodic("eig-"+num);
322 : }
323 15 : addComponentWithDerivatives("residual"); componentIsNotPeriodic("residual");
324 :
325 : // Get appropriate number of derivatives
326 : unsigned nder;
327 5 : if( getNumberOfAtoms()>0 ) {
328 5 : nder = 3*getNumberOfAtoms() + 9 + getNumberOfArguments();
329 : } else {
330 : nder = getNumberOfArguments();
331 : }
332 :
333 : // Resize all derivative arrays
334 5 : forces.resize( nder ); forcesToApply.resize( nder );
335 20 : for(int i=0; i<getNumberOfComponents(); ++i) getPntrToComponent(i)->resizeDerivatives(nder);
336 20 : }
337 :
338 36 : unsigned PCAVars::getNumberOfDerivatives() {
339 36 : if( getNumberOfAtoms()>0 ) {
340 36 : return 3*getNumberOfAtoms() + 9 + getNumberOfArguments();
341 : }
342 0 : return getNumberOfArguments();
343 : }
344 :
345 55 : void PCAVars::lockRequests() {
346 : ActionWithArguments::lockRequests();
347 : ActionAtomistic::lockRequests();
348 55 : }
349 :
350 55 : void PCAVars::unlockRequests() {
351 : ActionWithArguments::unlockRequests();
352 : ActionAtomistic::unlockRequests();
353 55 : }
354 :
355 55 : void PCAVars::calculate() {
356 :
357 55 : if(!nopbc && getNumberOfAtoms()>0) makeWhole();
358 :
359 : // Clear the reference value pack
360 55 : mypack.clear();
361 : // Calculate distance between instaneous configuration and reference
362 55 : double dist = myref->calculate( getPositions(), getPbc(), getArguments(), mypack, true );
363 :
364 : // Start accumulating residual by adding derivatives of distance
365 55 : Value* resid=getPntrToComponent( getNumberOfComponents()-1 ); unsigned nargs=getNumberOfArguments();
366 55 : for(unsigned j=0; j<getNumberOfArguments(); ++j) resid->addDerivative( j, mypack.getArgumentDerivative(j) );
367 440 : for(unsigned j=0; j<getNumberOfAtoms(); ++j) {
368 385 : Vector ader=mypack.getAtomDerivative( j );
369 1540 : for(unsigned k=0; k<3; ++k) resid->addDerivative( nargs +3*j+k, ader[k] );
370 : }
371 : // Retrieve the values of all arguments
372 55 : std::vector<double> args( getNumberOfArguments() ); for(unsigned i=0; i<getNumberOfArguments(); ++i) args[i]=getArgument(i);
373 :
374 : // Now calculate projections on pca vectors
375 55 : Vector adif, ader; Tensor fvir, tvir;
376 165 : for(int i=0; i<getNumberOfComponents()-1; ++i) { // One less component as we also have residual
377 110 : double proj=myref->projectDisplacementOnVector( directions[i], getArguments(), args, mypack );
378 : // And now accumulate derivatives
379 110 : Value* eid=getPntrToComponent(i);
380 110 : for(unsigned j=0; j<getNumberOfArguments(); ++j) eid->addDerivative( j, mypack.getArgumentDerivative(j) );
381 110 : if( getNumberOfAtoms()>0 ) {
382 110 : tvir.zero();
383 880 : for(unsigned j=0; j<getNumberOfAtoms(); ++j) {
384 770 : Vector myader=mypack.getAtomDerivative(j);
385 3080 : for(unsigned k=0; k<3; ++k) {
386 2310 : eid->addDerivative( nargs + 3*j+k, myader[k] );
387 2310 : resid->addDerivative( nargs + 3*j+k, -2*proj*myader[k] );
388 : }
389 770 : tvir += -1.0*Tensor( getPosition(j), myader );
390 : }
391 440 : for(unsigned j=0; j<3; ++j) {
392 1320 : for(unsigned k=0; k<3; ++k) eid->addDerivative( nargs + 3*getNumberOfAtoms() + 3*j + k, tvir(j,k) );
393 : }
394 : }
395 110 : dist -= proj*proj; // Subtract square from total squared distance to get residual squared
396 : // Derivatives of residual
397 110 : for(unsigned j=0; j<getNumberOfArguments(); ++j) resid->addDerivative( j, -2*proj*eid->getDerivative(j) );
398 : // for(unsigned j=0;j<getNumberOfArguments();++j) resid->addDerivative( j, -2*proj*arg_eigv(i,j) );
399 : // And set final value
400 110 : getPntrToComponent(i)->set( proj );
401 : }
402 55 : dist=std::sqrt(dist);
403 : resid->set( dist );
404 :
405 : // Take square root of residual derivatives
406 55 : double prefactor = 0.5 / dist;
407 55 : for(unsigned j=0; j<getNumberOfArguments(); ++j) resid->setDerivative( j, prefactor*resid->getDerivative(j) );
408 440 : for(unsigned j=0; j<getNumberOfAtoms(); ++j) {
409 1540 : for(unsigned k=0; k<3; ++k) resid->setDerivative( nargs + 3*j+k, prefactor*resid->getDerivative( nargs+3*j+k ) );
410 : }
411 :
412 : // And finally virial for residual
413 55 : if( getNumberOfAtoms()>0 ) {
414 55 : tvir.zero();
415 440 : for(unsigned j=0; j<getNumberOfAtoms(); ++j) {
416 1540 : Vector ader; for(unsigned k=0; k<3; ++k) ader[k]=resid->getDerivative( nargs + 3*j+k );
417 385 : tvir += -1.0*Tensor( getPosition(j), ader );
418 : }
419 220 : for(unsigned j=0; j<3; ++j) {
420 660 : for(unsigned k=0; k<3; ++k) resid->addDerivative( nargs + 3*getNumberOfAtoms() + 3*j + k, tvir(j,k) );
421 : }
422 : }
423 :
424 55 : }
425 :
426 0 : void PCAVars::calculateNumericalDerivatives( ActionWithValue* a ) {
427 0 : if( getNumberOfArguments()>0 ) {
428 0 : ActionWithArguments::calculateNumericalDerivatives( a );
429 : }
430 0 : if( getNumberOfAtoms()>0 ) {
431 0 : Matrix<double> save_derivatives( getNumberOfComponents(), getNumberOfArguments() );
432 0 : for(int j=0; j<getNumberOfComponents(); ++j) {
433 0 : for(unsigned i=0; i<getNumberOfArguments(); ++i) save_derivatives(j,i)=getPntrToComponent(j)->getDerivative(i);
434 : }
435 0 : calculateAtomicNumericalDerivatives( a, getNumberOfArguments() );
436 0 : for(int j=0; j<getNumberOfComponents(); ++j) {
437 0 : for(unsigned i=0; i<getNumberOfArguments(); ++i) getPntrToComponent(j)->addDerivative( i, save_derivatives(j,i) );
438 : }
439 : }
440 0 : }
441 :
442 55 : void PCAVars::apply() {
443 :
444 55 : bool wasforced=false; forcesToApply.assign(forcesToApply.size(),0.0);
445 220 : for(int i=0; i<getNumberOfComponents(); ++i) {
446 165 : if( getPntrToComponent(i)->applyForce( forces ) ) {
447 : wasforced=true;
448 0 : for(unsigned i=0; i<forces.size(); ++i) forcesToApply[i]+=forces[i];
449 : }
450 : }
451 55 : if( wasforced ) {
452 0 : addForcesOnArguments( forcesToApply );
453 0 : if( getNumberOfAtoms()>0 ) setForcesOnAtoms( forcesToApply, getNumberOfArguments() );
454 : }
455 :
456 55 : }
457 :
458 : }
459 : }
|
