Fixed to ensure that regtests gridtools/rt-spherical-integral-2, symf…

…unc/rt-two-filter and crystdistrib/rt-bops all work with new force scheme

regtest/gridtools/rt-spherical-integral-2/config

@@ -3,10 +3,10 @@ plumed_modules=adjmat
 # this is to test a different name
 arg="--plumed plumed.dat --ixyz trajectory.xyz --dump-forces forces --dump-forces-fmt %8.4f" # --debug-forces forces.num"
 
-# function plumed_regtest_before() {
-#    export PLUMED_FORBID_CHAINS=yes
-# }
-# 
-# function plumed_regtest_after() {
-#    export PLUMED_FORBID_CHAINS=no
-# }
+function plumed_regtest_before() {
+   export PLUMED_FORBID_CHAINS=yes
+}
+
+function plumed_regtest_after() {
+   export PLUMED_FORBID_CHAINS=no
+}

regtest/symfunc/rt-two-filter/config

@@ -2,3 +2,11 @@ type=driver
 plumed_modules="symfunc adjmat"
 # this is to test a different name
 arg="--plumed plumed.dat --ixyz gas-one.xyz --dump-forces forces --dump-forces-fmt=%8.4f"  #--debug-forces forces.num"
+
+function plumed_regtest_before() {
+   export PLUMED_FORBID_CHAINS=yes
+}
+
+function plumed_regtest_after() {
+   export PLUMED_FORBID_CHAINS=no
+}

regtest/symfunc/rt-two-filter/plumed.dat

@@ -2,7 +2,7 @@ c1: COORDINATIONNUMBER SPECIES=1-150 SWITCH={EXP D_0=4.0 R_0=0.5 D_MAX=6.0}
 c1m: MORE_THAN ARG=c1 SWITCH={RATIONAL D_0=2.0 R_0=0.1}
 
 con: CONTACT_MATRIX GROUP=c1 SWITCH={EXP D_0=4.0 R_0=0.5 D_MAX=6.0}
-dot: OUTER_PRODUCT ARG=c1m,c1m
+dot: OUTER_PRODUCT ARG=c1m,c1m MASK=con
 mm: MATHEVAL ARG=con,dot FUNC=x*y PERIODIC=NO
 
 ones: ONES SIZE=150

src/crystdistrib/QuaternionBondProductMatrix.cpp

@@ -127,9 +127,9 @@ unsigned QuaternionBondProductMatrix::getNumberOfColumns() const {
 }
 
 void QuaternionBondProductMatrix::setupForTask( const unsigned& task_index, std::vector<unsigned>& indices, MultiValue& myvals ) const {
-  unsigned start_n = getPntrToArgument(4)->getShape()[0], size_v = getPntrToArgument(4)->getShape()[1];
+  Value* myarg = getPntrToArgument(4); unsigned start_n = myarg->getShape()[0], size_v = myarg->getRowLength(task_index);
   if( indices.size()!=size_v+1 ) indices.resize( size_v+1 );
-  for(unsigned i=0; i<size_v; ++i) indices[i+1] = start_n + i;
+  for(unsigned i=0; i<size_v; ++i) indices[i+1] = start_n + myarg->getRowIndex( task_index, i);
   myvals.setSplitIndex( size_v + 1 );
 }
 

src/function/FunctionOfMatrix.h

@@ -343,8 +343,7 @@ void FunctionOfMatrix<T>::performTask( const std::string& controller, const unsi
       if( getPntrToArgument(j)->getRank()==2 ) {
         for(int i=0; i<getNumberOfComponents(); ++i) {
           unsigned ostrn=getConstPntrToComponent(i)->getPositionInStream();
-          unsigned myind = base + getPntrToArgument(j)->getNumberOfColumns()*index1 + ind2;
-          if( getConstPntrToComponent(i)->valueIsStored() && getConstPntrToComponent(i)->getNumberOfColumns()<getConstPntrToComponent(i)->getShape()[1] ) myind = base + getPntrToArgument(j)->getNumberOfColumns()*index1 + getConstPntrToComponent(i)->getRowLength(index1);
+          unsigned myind = base + getPntrToArgument(j)->getIndexInStore( index1*getPntrToArgument(j)->getShape()[1] + ind2 );  
           myvals.addDerivative( ostrn, myind, derivatives(i,j) );
           myvals.updateIndex( ostrn, myind );
         }

src/matrixtools/OuterProduct.cpp

@@ -55,7 +55,7 @@ class OuterProduct : public ActionWithMatrix {
 PLUMED_REGISTER_ACTION(OuterProduct,"OUTER_PRODUCT")
 
 void OuterProduct::registerKeywords( Keywords& keys ) {
-  ActionWithMatrix::registerKeywords(keys); keys.use("ARG");
+  ActionWithMatrix::registerKeywords(keys); keys.use("ARG"); keys.use("MASK");
   keys.add("compulsory","FUNC","x*y","the function of the input vectors that should be put in the elements of the outer product");
   keys.addFlag("ELEMENTS_ON_DIAGONAL_ARE_ZERO",false,"set all diagonal elements to zero");
   keys.setValueDescription("a matrix containing the outer product of the two input vectors that was obtained using the function that was input");
@@ -67,9 +67,15 @@ OuterProduct::OuterProduct(const ActionOptions&ao):
   domin(false),
   domax(false)
 {
-  if( getNumberOfArguments()!=2 ) error("should be two arguments to this action, a matrix and a vector");
+  unsigned nargs=getNumberOfArguments(); if( getNumberOfMasks()>0 ) nargs = nargs - getNumberOfMasks();
+  if( nargs!=2 ) error("should be two arguments to this action, they should both be vectors");
   if( getPntrToArgument(0)->getRank()!=1 || getPntrToArgument(0)->hasDerivatives() ) error("first argument to this action should be a vector");
   if( getPntrToArgument(1)->getRank()!=1 || getPntrToArgument(1)->hasDerivatives() ) error("first argument to this action should be a vector");
+  if( getNumberOfMasks()==1 ) {
+      if( getPntrToArgument(2)->getRank()!=2 || getPntrToArgument(2)->hasDerivatives() ) error("mask argument should be a matrix");
+      if( getPntrToArgument(2)->getShape()[0]!=getPntrToArgument(0)->getShape()[0] ) error("mask argument has wrong size");
+      if( getPntrToArgument(2)->getShape()[1]!=getPntrToArgument(1)->getShape()[0] ) error("mask argument has wrong size");
+  }
 
   std::string func; parse("FUNC",func);
   if( func=="min") {
@@ -106,7 +112,15 @@ void OuterProduct::prepare() {
 }
 
 void OuterProduct::setupForTask( const unsigned& task_index, std::vector<unsigned>& indices, MultiValue& myvals ) const {
-  unsigned start_n = getPntrToArgument(0)->getShape()[0], size_v = getPntrToArgument(1)->getShape()[0];
+  unsigned start_n = getPntrToArgument(0)->getShape()[0];
+  if( getNumberOfMasks()>0 ) {
+      Value* maskarg = getPntrToArgument(2); unsigned size_v = maskarg->getRowLength(task_index);
+      if( indices.size()!=size_v+1 ) indices.resize( size_v+1 );
+      for(unsigned i=0; i<size_v; ++i) indices[i+1] = start_n + maskarg->getRowIndex( task_index, i );
+      myvals.setSplitIndex( 1 + size_v ); return;
+  }
+
+  unsigned size_v = getPntrToArgument(1)->getShape()[0];
   if( diagzero ) {
     if( indices.size()!=size_v ) indices.resize( size_v );
     unsigned k=1;