Fixed code that say regtest adjmat/rt45 works as well

plumed · Sep 6, 2024 · aaa48f6 · aaa48f6 · PlumedBot · Sep 6, 2024
1 parent 8080dcb
commit aaa48f6
Show file tree

Hide file tree

Showing 4 changed files with 28 additions and 7 deletions.
diff --git a/regtest/adjmat/rt45/config b/regtest/adjmat/rt45/config
@@ -3,3 +3,11 @@ plumed_modules=adjmat
 # this is to test a different name
 arg="--plumed plumed.dat --trajectory-stride 10 --timestep 0.005 --ixyz trajectory.xyz --dump-forces forces --dump-forces-fmt=%8.4f" # --debug-forces=forces.num"
 extra_files="../../trajectories/trajectory.xyz"
+
+function plumed_regtest_before() {
+   export PLUMED_FORBID_CHAINS=yes
+}
+
+function plumed_regtest_after() {
+   export PLUMED_FORBID_CHAINS=no
+}
diff --git a/regtest/adjmat/rt45/plumed.dat b/regtest/adjmat/rt45/plumed.dat
@@ -11,9 +11,6 @@ a2: COORD_ANGLES ...
   BETWEEN1={GAUSSIAN LOWER=0.25pi UPPER=0.75pi}
 ... 
 
-PRINT ARG=a1_wsum FILE=matrix
-
-DUMPDERIVATIVES ARG=a1_sum FILE=deriv FMT=%8.4f
 a2b: CUSTOM ARG=a2_between-1,a2_denom FUNC=x/y PERIODIC=NO
 PRINT ARG=a1.*,a2.*,a2b FILE=colvar FMT=%8.5f
 BIASVALUE ARG=a1_mean

diff --git a/src/matrixtools/MatrixTimesMatrix.cpp b/src/matrixtools/MatrixTimesMatrix.cpp
@@ -46,6 +46,7 @@ namespace matrixtools {
 class MatrixTimesMatrix : public ActionWithMatrix {
 private:
   bool squared;
+  bool diagzero;
   unsigned nderivatives;
   bool stored_matrix1, stored_matrix2;
 public:
@@ -66,6 +67,7 @@ PLUMED_REGISTER_ACTION(MatrixTimesMatrix,"DISSIMILARITIES")
 void MatrixTimesMatrix::registerKeywords( Keywords& keys ) {
   ActionWithMatrix::registerKeywords(keys); keys.use("ARG"); keys.use("MASK");
   keys.addFlag("SQUARED",false,"calculate the squares of the dissimilarities (this option cannot be used with MATRIX_PRODUCT)");
+  keys.addFlag("ELEMENTS_ON_DIAGONAL_ARE_ZERO",false,"set all diagonal elements to zero");
   keys.setValueDescription("the product of the two input matrices");
 }
 
@@ -83,6 +85,9 @@ MatrixTimesMatrix::MatrixTimesMatrix(const ActionOptions&ao):
   std::string headstr=getFirstActionInChain()->getLabel();
   stored_matrix1 = getPntrToArgument(0)->ignoreStoredValue( headstr );
   stored_matrix2 = getPntrToArgument(1)->ignoreStoredValue( headstr );
+  parseFlag("ELEMENTS_ON_DIAGONAL_ARE_ZERO",diagzero);
+  if( diagzero ) log.printf("  setting diagonal elements equal to zero\n");
+
   if( getName()=="DISSIMILARITIES" ) {
     parseFlag("SQUARED",squared);
     if( squared ) log.printf("  calculating the squares of the dissimilarities \n");
@@ -119,13 +124,24 @@ void MatrixTimesMatrix::getAdditionalTasksRequired( ActionWithVector* action, st
 }
 
 void MatrixTimesMatrix::setupForTask( const unsigned& task_index, std::vector<unsigned>& indices, MultiValue& myvals ) const {
+  unsigned start_n = getPntrToArgument(0)->getShape()[0]; 
   if( getNumberOfArguments()>2 ) {
-    unsigned start_n = getPntrToArgument(0)->getShape()[0], size_v = getPntrToArgument(2)->getRowLength(task_index);
+    unsigned size_v = getPntrToArgument(2)->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 + getPntrToArgument(2)->getRowIndex(task_index, i);
-    myvals.setSplitIndex( size_v + 1 );
+    myvals.setSplitIndex( size_v + 1 ); return;
+  } 
+
+  unsigned size_v = getPntrToArgument(1)->getShape()[1];
+  if( diagzero ) {
+    if( indices.size()!=size_v ) indices.resize( size_v );
+    unsigned k=1;
+    for(unsigned i=0; i<size_v; ++i) {
+      if( task_index==i ) continue ;
+      indices[k] = size_v + i; k++;
+    }
+    myvals.setSplitIndex( size_v );
   } else {
-    unsigned start_n = getPntrToArgument(0)->getShape()[0], size_v = getPntrToArgument(1)->getShape()[1];
     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;
     myvals.setSplitIndex( size_v + 1 );

diff --git a/src/multicolvar/CoordAngles.cpp b/src/multicolvar/CoordAngles.cpp
@@ -92,7 +92,7 @@ CoordAngles::CoordAngles(const ActionOptions& ao):
   // Avoid double counting
   readInputLine( getShortcutLabel() + "_wmat: CUSTOM ARG=" + getShortcutLabel() + "_swd FUNC=0.5*x PERIODIC=NO");
   // And the matrix of dot products and the angles
-  readInputLine( getShortcutLabel() + "_dpmat: MATRIX_PRODUCT ARG=" + getShortcutLabel() + "_stack," + getShortcutLabel() + "_stackT");
+  readInputLine( getShortcutLabel() + "_dpmat: MATRIX_PRODUCT ELEMENTS_ON_DIAGONAL_ARE_ZERO ARG=" + getShortcutLabel() + "_stack," + getShortcutLabel() + "_stackT");
   readInputLine( getShortcutLabel() + "_angles: CUSTOM ARG=" + getShortcutLabel() + "_dpmat FUNC=acos(x) PERIODIC=NO");
   // Read the input
   Keywords keys; MultiColvarShortcuts::shortcutKeywords( keys ); pruneShortcuts( keys ); bool do_mean; parseFlag("MEAN",do_mean);