DRMSD and SECONDARY STRUCTURES: small optimisations

src/reference/DRMSD.cpp

@@ -76,38 +76,38 @@ double DRMSD::calc( const std::vector<Vector>& pos, const Pbc& pbc, ReferenceVal
   plumed_dbg_assert(!targets.empty());
 
   Vector distance; 
-  myder.clear(); double drmsd=0.; 
+  myder.clear(); 
+  double drmsd=0.; 
   for(std::map< std::pair <unsigned,unsigned> , double>::const_iterator it=targets.begin();it!=targets.end();++it){
       
-      unsigned i=getAtomIndex( it->first.first );
-      unsigned j=getAtomIndex( it->first.second );
+      const unsigned i=getAtomIndex( it->first.first );
+      const unsigned j=getAtomIndex( it->first.second );
 
-      if(nopbc){ distance=delta( pos[i] , pos[j] ); }
-      else{ distance=pbc.distance( pos[i] , pos[j] ); }
+      if(nopbc) distance=delta( pos[i] , pos[j] ); 
+      else      distance=pbc.distance( pos[i] , pos[j] );
 
-      double len = distance.modulo();
-      double diff = len - it->second;
+      const double len = distance.modulo();
+      const double diff = len - it->second;
+      const double der = diff / len; 
 
       drmsd += diff * diff;
-      myder.addAtomDerivatives( i, -( diff / len ) * distance );
-      myder.addAtomDerivatives( j, ( diff / len ) * distance );
-      myder.addBoxDerivatives( -( diff / len ) * Tensor(distance,distance) );
+      myder.addAtomDerivatives( i, -der * distance );
+      myder.addAtomDerivatives( j,  der * distance );
+      myder.addBoxDerivatives( - der * Tensor(distance,distance) );
   }
 
-  double npairs = static_cast<double>(targets.size());
+  const double inpairs = 1./static_cast<double>(targets.size());
   double idrmsd;
 
   if(squared){
-     drmsd = drmsd / npairs;
-     idrmsd = 2.0 / npairs;
+     drmsd = drmsd * inpairs;
+     idrmsd = 2.0 * inpairs;
   } else {
-     drmsd = sqrt( drmsd / npairs );
-     idrmsd = 1.0/( drmsd * npairs );
+     drmsd = sqrt( drmsd * inpairs );
+     idrmsd = inpairs / drmsd ;
   }
 
   myder.scaleAllDerivatives( idrmsd );
-  // virial *= idrmsd; 
-  // for(unsigned i=0;i<getNumberOfAtoms();++i){atom_ders[i] *= idrmsd;}
 
   return drmsd;
 }

src/secondarystructure/SecondaryStructureRMSD.cpp

@@ -76,7 +76,7 @@ ActionWithValue(ao),
 ActionWithVessel(ao),
 updateFreq(0),
 align_strands(false),
-s_cutoff(0),
+s_cutoff2(0),
 align_atom_1(0),
 align_atom_2(0)
 {
@@ -95,8 +95,10 @@ align_atom_2(0)
   }
 
   if( keywords.exists("STRANDS_CUTOFF") ){
+    double s_cutoff = 0;
     parse("STRANDS_CUTOFF",s_cutoff); align_strands=true;
     if( s_cutoff>0) log.printf("  ignoring contributions from strands that are more than %f apart\n",s_cutoff);
+    s_cutoff2=s_cutoff*s_cutoff;
   }
 }
 
@@ -201,12 +203,13 @@ void SecondaryStructureRMSD::calculate(){
 void SecondaryStructureRMSD::performTask( const unsigned& task_index, const unsigned& current, MultiValue& myvals ) const {
   // Retrieve the positions
   std::vector<Vector> pos( references[0]->getNumberOfAtoms() );
-  for(unsigned i=0;i<pos.size();++i) pos[i]=ActionAtomistic::getPosition( getAtomIndex(current,i) );
+  const unsigned n=pos.size();
+  for(unsigned i=0;i<n;++i) pos[i]=ActionAtomistic::getPosition( getAtomIndex(current,i) );
 
   // This does strands cutoff
   Vector distance=pbcDistance( pos[align_atom_1],pos[align_atom_2] ); 
-  if( s_cutoff>0 ){
-     if( distance.modulo()>s_cutoff ){
+  if( s_cutoff2>0 ){
+     if( distance.modulo2()>s_cutoff2 ){
        myvals.setValue( 0, 0.0 );
        return;
      }
@@ -222,15 +225,17 @@ void SecondaryStructureRMSD::performTask( const unsigned& task_index, const unsi
   }
   // Create a holder for the derivatives
   ReferenceValuePack mypack( 0, pos.size(), myvals ); mypack.setValIndex( 1 );
-  for(unsigned i=0;i<pos.size();++i) mypack.setAtomIndex( i, getAtomIndex(current,i) );
+  for(unsigned i=0;i<n;++i) mypack.setAtomIndex( i, getAtomIndex(current,i) );
 
   // And now calculate the RMSD
-  double r,nr; const Pbc& pbc=getPbc(); 
-  unsigned closest=0; r=references[0]->calculate( pos, pbc, mypack, false );
-  for(unsigned i=1;i<references.size();++i){
-      mypack.setValIndex( i+1 );
-      nr=references[i]->calculate( pos, pbc, mypack, false );
-      if( nr<r ){ closest=i; r=nr; }
+  const Pbc& pbc=getPbc(); 
+  unsigned closest=0; 
+  double r = references[0]->calculate( pos, pbc, mypack, false );
+  const unsigned rs = references.size();
+  for(unsigned i=1;i<rs;++i){
+    mypack.setValIndex( i+1 );
+    double nr=references[i]->calculate( pos, pbc, mypack, false );
+    if( nr<r ){ closest=i; r=nr; }
   }
 
   // Transfer everything to the value
@@ -238,11 +243,12 @@ void SecondaryStructureRMSD::performTask( const unsigned& task_index, const unsi
   if( closest>0 ) mypack.moveDerivatives( closest+1, 1 );
 
   if( !mypack.virialWasSet() ){
-      Tensor vir; vir.zero();
-      for(unsigned i=0;i<colvar_atoms[current].size();++i){
-         vir+=(-1.0*Tensor( pos[i], mypack.getAtomDerivative(i) ));
-      } 
-      mypack.setValIndex(1); mypack.addBoxDerivatives( vir );
+    Tensor vir;
+    const unsigned cacs = colvar_atoms[current].size();
+    for(unsigned i=0;i<cacs;++i){
+       vir+=(-1.0*Tensor( pos[i], mypack.getAtomDerivative(i) ));
+    } 
+    mypack.setValIndex(1); mypack.addBoxDerivatives( vir );
   }
 
   return;

src/secondarystructure/SecondaryStructureRMSD.h

@@ -54,7 +54,7 @@ private:
   bool firsttime;
 /// Variables for strands cutoff
   bool align_strands;
-  double s_cutoff;
+  double s_cutoff2;
   unsigned align_atom_1, align_atom_2;
   bool verbose_output;
 /// Tempory variables for getting positions of atoms and applying forces