@@ -46,6 +46,8 @@ class MatrixTimesVector : public ActionWithVector {
46
46
void calculate () override ;
47
47
void performTask ( const unsigned & task_index, MultiValue& myvals ) const override ;
48
48
int checkTaskIsActive ( const unsigned & itask ) const override ;
49
+ void getNumberOfForceDerivatives ( unsigned & nforces, unsigned & nderiv ) const override ;
50
+ void gatherForces ( const unsigned & itask, const MultiValue& myvals, std::vector<double >& forces ) const override ;
49
51
};
50
52
51
53
PLUMED_REGISTER_ACTION (MatrixTimesVector," MATRIX_VECTOR_PRODUCT" )
@@ -173,7 +175,7 @@ int MatrixTimesVector::checkTaskIsActive( const unsigned& itask ) const {
173
175
174
176
void MatrixTimesVector::performTask ( const unsigned & task_index, MultiValue& myvals ) const {
175
177
if ( sumrows ) {
176
- unsigned base= 0 , n=getNumberOfArguments ()-1 ; Value* myvec = getPntrToArgument (n);
178
+ unsigned n=getNumberOfArguments ()-1 ; Value* myvec = getPntrToArgument (n);
177
179
for (unsigned i=0 ; i<n; ++i) {
178
180
Value* mymat = getPntrToArgument (i);
179
181
unsigned ncol = mymat->getNumberOfColumns ();
@@ -184,11 +186,10 @@ void MatrixTimesVector::performTask( const unsigned& task_index, MultiValue& myv
184
186
// And the derivatives
185
187
if ( doNotCalculateDerivatives () ) continue ;
186
188
187
- unsigned dloc = base + task_index*ncol;
189
+ unsigned dloc = task_index*ncol;
188
190
for (unsigned j=0 ; j<nmat; ++j) {
189
191
myvals.addDerivative ( i, dloc + j, 1.0 ); myvals.updateIndex ( i, dloc + j );
190
192
}
191
- base += mymat->getNumberOfStoredValues ();
192
193
}
193
194
} else if ( getPntrToArgument (1 )->getRank ()==1 ) {
194
195
Value* mymat = getPntrToArgument (0 );
@@ -238,5 +239,28 @@ void MatrixTimesVector::performTask( const unsigned& task_index, MultiValue& myv
238
239
}
239
240
}
240
241
242
+ void MatrixTimesVector::getNumberOfForceDerivatives ( unsigned & nforces, unsigned & nderiv ) const {
243
+ ActionWithVector::getNumberOfForceDerivatives ( nforces, nderiv );
244
+ if ( sumrows ) nderiv = getPntrToArgument (0 )->getNumberOfStoredValues () + getPntrToArgument (getNumberOfArguments ()-1 )->getNumberOfStoredValues ();
245
+ }
246
+
247
+ void MatrixTimesVector::gatherForces ( const unsigned & itask, const MultiValue& myvals, std::vector<double >& forces ) const {
248
+ if ( !sumrows ) { ActionWithVector::gatherForces ( itask, myvals, forces ); return ; }
249
+ if ( checkComponentsForForce () ) {
250
+ unsigned base = 0 ;
251
+ for (unsigned ival=0 ; ival<getNumberOfComponents (); ++ival) {
252
+ const Value* myval=getConstPntrToComponent (ival);
253
+ if ( myval->forcesWereAdded () ) {
254
+ double fforce = myval->getForce (itask);
255
+ for (unsigned j=0 ; j<myvals.getNumberActive (ival); ++j) {
256
+ unsigned jder=myvals.getActiveIndex (ival, j); plumed_dbg_assert ( jder<forces.size () );
257
+ forces[base+jder] += fforce*myvals.getDerivative ( ival, jder );
258
+ }
259
+ }
260
+ base += getPntrToArgument (ival)->getNumberOfStoredValues ();
261
+ }
262
+ }
263
+ }
264
+
241
265
}
242
266
}
0 commit comments