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/ActionWithValue.h"
23 : #include "core/ActionAtomistic.h"
24 : #include "core/ActionWithArguments.h"
25 : #include "reference/MultiReferenceBase.h"
26 : #include "reference/MetricRegister.h"
27 : #include "core/ActionRegister.h"
28 : #include "core/PlumedMain.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 : instaneous 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 explicitally 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 : \verbatim
79 : PCAVARS REFERENCE=reference.pdb TYPE=OPTIMAL LABEL=pca2
80 : PRINT ARG=pca2.* FILE=colvar2
81 : \endverbatim
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 refered 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 : \verbatim
90 : ATOM 2 CH3 ACE 1 12.932 -14.718 -6.016 1.00 1.00
91 : ATOM 5 C ACE 1 21.312 -9.928 -5.946 1.00 1.00
92 : ATOM 9 CA ALA 2 19.462 -11.088 -8.986 1.00 1.00
93 : ATOM 13 HB2 ALA 2 21.112 -10.688 -12.476 1.00 1.00
94 : ATOM 15 C ALA 2 19.422 7.978 -14.536 1.00 1.00
95 : ATOM 20 HH31 NME 3 20.122 -9.928 -17.746 1.00 1.00
96 : ATOM 21 HH32 NME 3 18.572 -13.148 -16.346 1.00 1.00
97 : END
98 : ATOM 2 CH3 ACE 1 13.932 -14.718 -6.016 1.00 1.00
99 : ATOM 5 C ACE 1 20.312 -9.928 -5.946 1.00 1.00
100 : ATOM 9 CA ALA 2 18.462 -11.088 -8.986 1.00 1.00
101 : ATOM 13 HB2 ALA 2 20.112 -11.688 -12.476 1.00 1.00
102 : ATOM 15 C ALA 2 19.422 7.978 -12.536 1.00 1.00
103 : ATOM 20 HH31 NME 3 20.122 -9.928 -17.746 1.00 1.00
104 : ATOM 21 HH32 NME 3 18.572 -13.148 -16.346 1.00 1.00
105 : END
106 : \endverbatim
107 :
108 : Alternatively, the second configuration can specify the components of \f$A\f$ explicitally. In this case you need to include the
109 : keyword TYPE=DIRECTION in the remarks to the pdb as shown below.
110 :
111 : \verbatim
112 : ATOM 2 CH3 ACE 1 12.932 -14.718 -6.016 1.00 1.00
113 : ATOM 5 C ACE 1 21.312 -9.928 -5.946 1.00 1.00
114 : ATOM 9 CA ALA 2 19.462 -11.088 -8.986 1.00 1.00
115 : ATOM 13 HB2 ALA 2 21.112 -10.688 -12.476 1.00 1.00
116 : ATOM 15 C ALA 2 19.422 7.978 -14.536 1.00 1.00
117 : ATOM 20 HH31 NME 3 20.122 -9.928 -17.746 1.00 1.00
118 : ATOM 21 HH32 NME 3 18.572 -13.148 -16.346 1.00 1.00
119 : END
120 : REMARK TYPE=DIRECTION
121 : ATOM 2 CH3 ACE 1 0.1414 0.3334 -0.0302 1.00 0.00
122 : ATOM 5 C ACE 1 0.0893 -0.1095 -0.1434 1.00 0.00
123 : ATOM 9 CA ALA 2 0.0207 -0.321 0.0321 1.00 0.00
124 : ATOM 13 HB2 ALA 2 0.0317 -0.6085 0.0783 1.00 0.00
125 : ATOM 15 C ALA 2 0.1282 -0.4792 0.0797 1.00 0.00
126 : ATOM 20 HH31 NME 3 0.0053 -0.465 0.0309 1.00 0.00
127 : ATOM 21 HH32 NME 3 -0.1019 -0.4261 -0.0082 1.00 0.00
128 : END
129 : \endverbatim
130 :
131 : If your metric involves arguments the labels of these arguments in your plumed input file should be specified in the REMARKS
132 : for each of the frames of your path. An input file in this case might look like this:
133 :
134 : \verbatim
135 : DESCRIPTION: a pca eigenvector specified using the start point and direction in the HD space.
136 : REMARK WEIGHT=1.0
137 : REMARK ARG=d1,d2
138 : REMARK d1=1.0 d2=1.0
139 : END
140 : REMARK TYPE=DIRECTION
141 : REMARK ARG=d1,d2
142 : REMARK d1=0.1 d2=0.25
143 : END
144 : \endverbatim
145 :
146 : Here we are working with the EUCLIDEAN metric and notice that we have specified the components of \f$A\f$ using DIRECTION.
147 : Consequently, the values of d1 and d2 in the second frame above do not specify a particular coordinate in the high-dimensional
148 : 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.
149 : If we wanted to specify the direction in this metric using the start and end point of the vector we would write:
150 :
151 : \verbatim
152 : DESCRIPTION: a pca eigenvector specified using the start and end point of a vector in the HD space.
153 : REMARK WEIGHT=1.0
154 : REMARK ARG=d1,d2
155 : REMARK d1=1.0 d2=1.0
156 : END
157 : REMARK ARG=d1,d2
158 : REMARK d1=1.1 d2=1.25
159 : END
160 : \endverbatim
161 :
162 : */
163 : //+ENDPLUMEDOC
164 :
165 : namespace PLMD {
166 : namespace mapping {
167 :
168 : class PCAVars :
169 : public ActionWithValue,
170 : public ActionAtomistic,
171 : public ActionWithArguments
172 : {
173 : private:
174 : /// The holders for the derivatives
175 : MultiValue myvals;
176 : ReferenceValuePack mypack;
177 : /// The position of the reference configuration (the one we align to)
178 : ReferenceConfiguration* myref;
179 : /// The eigenvectors for the atomic displacements
180 : Matrix<Vector> atom_eigv;
181 : /// The eigenvectors for the displacements in argument space
182 : Matrix<double> arg_eigv;
183 : /// Stuff for applying forces
184 : std::vector<double> forces, forcesToApply;
185 : public:
186 : static void registerKeywords( Keywords& keys );
187 : explicit PCAVars(const ActionOptions&);
188 : ~PCAVars();
189 : unsigned getNumberOfDerivatives();
190 : void lockRequests();
191 : void unlockRequests();
192 : void calculateNumericalDerivatives( ActionWithValue* a );
193 : void calculate();
194 : void apply();
195 : };
196 :
197 2527 : PLUMED_REGISTER_ACTION(PCAVars,"PCAVARS")
198 :
199 5 : void PCAVars::registerKeywords( Keywords& keys ) {
200 5 : Action::registerKeywords( keys );
201 5 : ActionWithValue::registerKeywords( keys );
202 5 : ActionAtomistic::registerKeywords( keys );
203 5 : ActionWithArguments::registerKeywords( keys );
204 5 : componentsAreNotOptional(keys);
205 5 : keys.addOutputComponent("eig","default","the projections on each eigenvalue are stored on values labeled eig-1, eig-2, ...");
206 : 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 5 : "reference point.");
211 5 : keys.add("compulsory","REFERENCE","a pdb file containing the reference configuration and configurations that define the directions for each eigenvector");
212 5 : keys.add("compulsory","TYPE","OPTIMAL","The method we are using for alignment to the reference structure");
213 5 : keys.addFlag("NORMALIZE",false,"calculate the length of the eigenvector input and divide the components by it so as to have a normalised vector");
214 5 : }
215 :
216 4 : PCAVars::PCAVars(const ActionOptions& ao):
217 : Action(ao),
218 : ActionWithValue(ao),
219 : ActionAtomistic(ao),
220 : ActionWithArguments(ao),
221 : myvals(1,0),
222 4 : mypack(0,0,myvals)
223 : {
224 :
225 : // What type of distance are we calculating
226 4 : std::string mtype; parse("TYPE",mtype);
227 :
228 : // Open reference file
229 8 : std::string reference; parse("REFERENCE",reference);
230 4 : FILE* fp=fopen(reference.c_str(),"r");
231 4 : if(!fp) error("could not open reference file " + reference );
232 :
233 : // Read all reference configurations
234 8 : MultiReferenceBase myframes( "", false );
235 4 : bool do_read=true; unsigned nfram=0;
236 20 : while (do_read) {
237 16 : PDB mypdb;
238 : // Read the pdb file
239 16 : do_read=mypdb.readFromFilepointer(fp,plumed.getAtoms().usingNaturalUnits(),0.1/atoms.getUnits().getLength());
240 : // Fix argument names
241 16 : expandArgKeywordInPDB( mypdb );
242 16 : if(do_read) {
243 12 : if( nfram==0 ) {
244 4 : myref = metricRegister().create<ReferenceConfiguration>( mtype, mypdb );
245 4 : if( myref->isDirection() ) error("first frame should be reference configuration - not direction of vector");
246 4 : if( !myref->pcaIsEnabledForThisReference() ) error("can't do PCA with reference type " + mtype );
247 8 : std::vector<std::string> remarks( mypdb.getRemark() ); std::string rtype;
248 4 : bool found=Tools::parse( remarks, "TYPE", rtype );
249 4 : if(!found) { std::vector<std::string> newrem(1); newrem[0]="TYPE="+mtype; mypdb.addRemark(newrem); }
250 8 : myframes.readFrame( mypdb );
251 8 : } else myframes.readFrame( mypdb );
252 12 : nfram++;
253 : } else {
254 4 : break;
255 : }
256 12 : }
257 4 : fclose(fp);
258 :
259 4 : if( nfram<=2 ) error("no eigenvectors were specified");
260 4 : log.printf(" found %u eigenvectors in file %s \n",nfram-1,reference.c_str() );
261 :
262 : // Finish the setup of the mapping object
263 : // Get the arguments and atoms that are required
264 8 : std::vector<AtomNumber> atoms; std::vector<std::string> args;
265 4 : myframes.getAtomAndArgumentRequirements( atoms, args );
266 8 : requestAtoms( atoms ); std::vector<Value*> req_args;
267 4 : interpretArgumentList( args, req_args ); requestArguments( req_args );
268 :
269 : // Setup the derivative pack
270 4 : if( atoms.size()>0 ) myvals.resize( 1, args.size() + 3*atoms.size() + 9 );
271 0 : else myvals.resize( 1, args.size() );
272 4 : mypack.resize( args.size(), atoms.size() );
273 4 : for(unsigned i=0; i<atoms.size(); ++i) mypack.setAtomIndex( i, i );
274 : /// This sets up all the storage data required by PCA in the pack
275 4 : myframes.getFrame(0)->setupPCAStorage( mypack );
276 :
277 : // Retrieve the position of the first frame, as we use this for alignment
278 4 : myref->setNamesAndAtomNumbers( atoms, args );
279 : // Check there are no periodic arguments
280 4 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
281 0 : if( getPntrToArgument(i)->isPeriodic() ) error("cannot use periodic variables in pca projections");
282 : }
283 : // Work out if the user wants to normalise the input vector
284 4 : bool nflag; parseFlag("NORMALIZE",nflag);
285 4 : checkRead();
286 :
287 : // Resize the matrices that will hold our eivenvectors
288 4 : if( getNumberOfAtoms()>0 ) atom_eigv.resize( nfram-1, getNumberOfAtoms() );
289 4 : if( getNumberOfArguments()>0 ) arg_eigv.resize( nfram-1, getNumberOfArguments() );
290 :
291 : // Create fake periodic boundary condition (these would only be used for DRMSD which is not allowed)
292 8 : Pbc fake_pbc;
293 : // Now calculate the eigenvectors
294 12 : for(unsigned i=1; i<nfram; ++i) {
295 : // Calculate distance from reference configuration
296 8 : double dist=myframes.getFrame(i)->calc( myref->getReferencePositions(), fake_pbc, getArguments(), myref->getReferenceArguments(), mypack, true );
297 :
298 : // Calculate the length of the vector for normalization
299 8 : double tmp, norm=0.0;
300 64 : for(unsigned j=0; j<getNumberOfAtoms(); ++j) {
301 56 : for(unsigned k=0; k<3; ++k) { tmp = mypack.getAtomsDisplacementVector()[j][k]; norm+=tmp*tmp; }
302 : }
303 8 : for(unsigned j=0; j<getNumberOfArguments(); ++j) { tmp = 0.5*mypack.getArgumentDerivative(j); norm+=tmp*tmp; }
304 :
305 : // Normalize the eigevector
306 8 : if(nflag) { norm = 1.0 / sqrt(norm); } else { norm = 1.0; }
307 8 : for(unsigned j=0; j<getNumberOfAtoms(); ++j) atom_eigv(i-1,j) = norm*mypack.getAtomsDisplacementVector()[j];
308 8 : for(unsigned j=0; j<getNumberOfArguments(); ++j) arg_eigv(i-1,j) = -0.5*norm*mypack.getArgumentDerivative(j);
309 :
310 : // Create a component to store the output
311 8 : std::string num; Tools::convert( i, num );
312 8 : addComponentWithDerivatives("eig-"+num); componentIsNotPeriodic("eig-"+num);
313 8 : }
314 4 : addComponentWithDerivatives("residual"); componentIsNotPeriodic("residual");
315 :
316 : // Get appropriate number of derivatives
317 : unsigned nder;
318 4 : if( getNumberOfAtoms()>0 ) {
319 4 : nder = 3*getNumberOfAtoms() + 9 + getNumberOfArguments();
320 : } else {
321 0 : nder = getNumberOfArguments();
322 : }
323 :
324 : // Resize all derivative arrays
325 4 : forces.resize( nder ); forcesToApply.resize( nder );
326 8 : for(unsigned i=0; i<getNumberOfComponents(); ++i) getPntrToComponent(i)->resizeDerivatives(nder);
327 4 : }
328 :
329 16 : PCAVars::~PCAVars() {
330 4 : delete myref;
331 12 : }
332 :
333 36 : unsigned PCAVars::getNumberOfDerivatives() {
334 36 : if( getNumberOfAtoms()>0 ) {
335 36 : return 3*getNumberOfAtoms() + 9 + getNumberOfArguments();
336 : }
337 0 : return getNumberOfArguments();
338 : }
339 :
340 44 : void PCAVars::lockRequests() {
341 44 : ActionWithArguments::lockRequests();
342 44 : ActionAtomistic::lockRequests();
343 44 : }
344 :
345 44 : void PCAVars::unlockRequests() {
346 44 : ActionWithArguments::unlockRequests();
347 44 : ActionAtomistic::unlockRequests();
348 44 : }
349 :
350 44 : void PCAVars::calculate() {
351 : // Clear the reference value pack
352 44 : mypack.clear();
353 : // Calculate distance between instaneous configuration and reference
354 44 : double dist = myref->calculate( getPositions(), getPbc(), getArguments(), mypack, true );
355 :
356 : // Start accumulating residual by adding derivatives of distance
357 44 : Value* resid=getPntrToComponent( getNumberOfComponents()-1 ); unsigned nargs=getNumberOfArguments();
358 44 : for(unsigned j=0; j<getNumberOfArguments(); ++j) resid->addDerivative( j, mypack.getArgumentDerivative(j) );
359 352 : for(unsigned j=0; j<getNumberOfAtoms(); ++j) {
360 308 : Vector ader=mypack.getAtomDerivative( j );
361 308 : for(unsigned k=0; k<3; ++k) resid->addDerivative( nargs +3*j+k, ader[k] );
362 : }
363 :
364 : // Now calculate projections on pca vectors
365 44 : Vector adif, ader; Tensor fvir, tvir;
366 132 : for(unsigned i=0; i<getNumberOfComponents()-1; ++i) { // One less component as we also have residual
367 88 : double proj=0; tvir.zero(); Value* eid=getPntrToComponent(i);
368 88 : for(unsigned j=0; j<getNumberOfArguments(); ++j) {
369 0 : proj+=arg_eigv(i,j)*0.5*mypack.getArgumentDerivative(j);
370 0 : eid->addDerivative( j, arg_eigv(i,j) );
371 : }
372 88 : if( getNumberOfAtoms()>0 ) {
373 88 : proj += myref->projectAtomicDisplacementOnVector( i, atom_eigv, getPositions(), mypack );
374 704 : for(unsigned j=0; j<getNumberOfAtoms(); ++j) {
375 616 : Vector myader=mypack.getAtomDerivative(j);
376 2464 : for(unsigned k=0; k<3; ++k) {
377 1848 : eid->addDerivative( nargs + 3*j+k, myader[k] );
378 1848 : resid->addDerivative( nargs + 3*j+k, -2*proj*myader[k] );
379 : }
380 616 : tvir += -1.0*Tensor( getPosition(j), myader );
381 : }
382 352 : for(unsigned j=0; j<3; ++j) {
383 264 : for(unsigned k=0; k<3; ++k) eid->addDerivative( nargs + 3*getNumberOfAtoms() + 3*j + k, tvir(j,k) );
384 : }
385 : }
386 88 : dist -= proj*proj; // Subtract square from total squared distance to get residual squared
387 : // Derivatives of residual
388 88 : for(unsigned j=0; j<getNumberOfArguments(); ++j) resid->addDerivative( j, -2*proj*arg_eigv(i,j) );
389 : // And set final value
390 88 : getPntrToComponent(i)->set( proj );
391 : }
392 44 : dist=sqrt(dist);
393 44 : resid->set( dist );
394 :
395 : // Take square root of residual derivatives
396 44 : double prefactor = 0.5 / dist;
397 44 : for(unsigned j=0; j<getNumberOfArguments(); ++j) resid->setDerivative( j, prefactor*resid->getDerivative(j) );
398 352 : for(unsigned j=0; j<getNumberOfAtoms(); ++j) {
399 308 : for(unsigned k=0; k<3; ++k) resid->setDerivative( nargs + 3*j+k, prefactor*resid->getDerivative( nargs+3*j+k ) );
400 : }
401 :
402 : // And finally virial for residual
403 44 : if( getNumberOfAtoms()>0 ) {
404 44 : tvir.zero();
405 352 : for(unsigned j=0; j<getNumberOfAtoms(); ++j) {
406 308 : Vector ader; for(unsigned k=0; k<3; ++k) ader[k]=resid->getDerivative( nargs + 3*j+k );
407 308 : tvir += -1.0*Tensor( getPosition(j), ader );
408 : }
409 176 : for(unsigned j=0; j<3; ++j) {
410 132 : for(unsigned k=0; k<3; ++k) resid->addDerivative( nargs + 3*getNumberOfAtoms() + 3*j + k, tvir(j,k) );
411 : }
412 : }
413 :
414 44 : }
415 :
416 0 : void PCAVars::calculateNumericalDerivatives( ActionWithValue* a ) {
417 0 : if( getNumberOfArguments()>0 ) {
418 0 : ActionWithArguments::calculateNumericalDerivatives( a );
419 : }
420 0 : if( getNumberOfAtoms()>0 ) {
421 0 : Matrix<double> save_derivatives( getNumberOfComponents(), getNumberOfArguments() );
422 0 : for(unsigned j=0; j<getNumberOfComponents(); ++j) {
423 0 : for(unsigned i=0; i<getNumberOfArguments(); ++i) save_derivatives(j,i)=getPntrToComponent(j)->getDerivative(i);
424 : }
425 0 : calculateAtomicNumericalDerivatives( a, getNumberOfArguments() );
426 0 : for(unsigned j=0; j<getNumberOfComponents(); ++j) {
427 0 : for(unsigned i=0; i<getNumberOfArguments(); ++i) getPntrToComponent(j)->addDerivative( i, save_derivatives(j,i) );
428 0 : }
429 : }
430 0 : }
431 :
432 44 : void PCAVars::apply() {
433 :
434 44 : bool wasforced=false; forcesToApply.assign(forcesToApply.size(),0.0);
435 176 : for(unsigned i=0; i<getNumberOfComponents(); ++i) {
436 132 : if( getPntrToComponent(i)->applyForce( forces ) ) {
437 0 : wasforced=true;
438 0 : for(unsigned i=0; i<forces.size(); ++i) forcesToApply[i]+=forces[i];
439 : }
440 : }
441 44 : if( wasforced ) {
442 0 : addForcesOnArguments( forcesToApply );
443 0 : if( getNumberOfAtoms()>0 ) setForcesOnAtoms( forcesToApply, getNumberOfArguments() );
444 : }
445 :
446 44 : }
447 :
448 : }
449 2523 : }
|
