PRE: removed mpi

RDC: cleaning

src/colvar/PRE.cpp

@@ -41,7 +41,6 @@ Calculates the PRE intensity ratio.
 class PRE : public Colvar {   
 private:
   bool             pbc;
-  bool             serial;
   double           constant, inept;
   vector<double>   rtwo;
   vector<unsigned> nga;
@@ -72,26 +71,21 @@ void PRE::registerKeywords( Keywords& keys ){
                    "Keywords like RTWO1, RTWO2, RTWO3,... should be listed.");
   keys.addFlag("ADDEXPVALUES",false,"Set to TRUE if you want to have fixed components with the experimetnal values.");  
   keys.add("numbered","PREINT","Add an experimental value for each PRE.");
-  keys.addFlag("SERIAL",false,"Perform the calculation in serial - for debug purpose");
   keys.addOutputComponent("pre","default","the # PRE");
   keys.addOutputComponent("exp","ADDEXPVALUES","the # PRE experimental intensity");
 }
 
 PRE::PRE(const ActionOptions&ao):
 PLUMED_COLVAR_INIT(ao),
-pbc(true),
-serial(false)
+pbc(true)
 {
-  parseFlag("SERIAL",serial);
-
   bool nopbc=!pbc;
   parseFlag("NOPBC",nopbc);
   pbc=!nopbc;  
 
   vector<AtomNumber> atom;
   parseAtomList("SPINLABEL",atom);
-  if(atom.size()!=1)
-    error("Number of specified atom should be 1");
+  if(atom.size()!=1) error("Number of specified atom should be 1");
 
   // Read in the atoms
   vector<AtomNumber> t, ga_lista, gb_lista;
@@ -167,8 +161,6 @@ serial(false)
     log.printf("\n");
   }
 
-  if(!serial)  log.printf("  The PREs are calculated in parallel\n");
-  else         log.printf("  The PREs are calculated in serial\n");
   if(pbc)      log.printf("  using periodic boundary conditions\n");
   else         log.printf("  without periodic boundary conditions\n");
 
@@ -195,78 +187,50 @@ PRE::~PRE(){
   delete nl;
 } 
 
-void PRE::calculate(){ 
-  unsigned sga = nga.size();
-  std::vector<Vector> deriv(getNumberOfAtoms());
-  std::vector<Tensor> dervir(sga);
-  vector<double> pre(sga), ratio(sga);
-
-  // internal parallelisation
-  unsigned stride=comm.Get_size();
-  unsigned rank=comm.Get_rank();
-  if(serial){
-    stride=1;
-    rank=0;
-  }
- 
-  for(unsigned i=0;i<sga;i++) ratio[i]=pre[i]=0.;
-
-  unsigned index=0; for(unsigned k=0;k<rank;k++) index += nga[k];
+void PRE::calculate()
+{
+  std::vector<Vector> deriv(getNumberOfAtoms()); 
+  unsigned index=0;
 
-  for(unsigned i=rank;i<sga;i+=stride) { //cycle over the number of pre 
+  for(unsigned i=0;i<nga.size();i++) { //cycle over the number of pre
+    Tensor dervir;
+    double pre=0;
     for(unsigned j=0;j<nga[i];j++) {
-      double aver=1./((double)nga[i]);
-      unsigned i0=nl->getClosePair(index).first;
-      unsigned i1=nl->getClosePair(index).second;
+      const double c_aver=constant/((double)nga[i]);
+      const unsigned i0=nl->getClosePair(index).first;
+      const unsigned i1=nl->getClosePair(index).second;
+
       Vector distance;
-      if(pbc){
-        distance=pbcDistance(getPosition(i0),getPosition(i1));
-      } else {
-        distance=delta(getPosition(i0),getPosition(i1));
-      }
-      double d=distance.modulo();
-      double r2=d*d;
-      double r6=r2*r2*r2;
-      double r8=r6*r2;
-      double tmpir6=constant*aver/r6;
-      double tmpir8=-6.*constant*aver/r8;
-
-      pre[i] += tmpir6;
+      if(pbc) distance=pbcDistance(getPosition(i0),getPosition(i1));
+      else    distance=delta(getPosition(i0),getPosition(i1));
 
-      deriv[i0] = -tmpir8*distance;
-      deriv[i1] = +tmpir8*distance;
+      const double r2=distance.modulo2();
+      const double r6=r2*r2*r2;
+      const double r8=r6*r2;
+      const double tmpir6=c_aver/r6;
+      const double tmpir8=-6.*c_aver/r8;
 
-      dervir[i] += Tensor(distance,deriv[i0]);
+      pre += tmpir6;
 
+      deriv[i0] = -tmpir8*distance;
+      deriv[i1] =  tmpir8*distance;
+      dervir   +=  Tensor(distance,deriv[i0]);
       index++;
     }
-    ratio[i] = rtwo[i]*exp(-pre[i]*inept) / (rtwo[i]+pre[i]);
-    for(unsigned k=i+1;k<i+stride;k++) index += nga[k];
-  }
-
-  if(!serial){
-    comm.Sum(&pre[0],sga);
-    comm.Sum(&ratio[0],sga);
-    comm.Sum(&deriv[0][0],3*deriv.size());
-    comm.Sum(&dervir[0][0][0],9*sga);
-  }
+    const double ratio = rtwo[i]*exp(-pre*inept) / (rtwo[i]+pre);
+    const double fact = -ratio*(inept+1./(rtwo[i]+pre));
 
-  index = 0;
-  for(unsigned i=0;i<sga;i++) { //cycle over the number of pre 
     Value* val=getPntrToComponent(i);
-    val->set(ratio[i]);
-    double fact = -ratio[i]*(inept+1./(rtwo[i]+pre[i]));
-    Tensor virial = fact*dervir[i];
-    setBoxDerivatives(val, virial);
+    val->set(ratio);
+    setBoxDerivatives(val, fact*dervir);
+
     for(unsigned j=0;j<nga[i];j++) {
-      unsigned i0=nl->getClosePair(index).first;
-      unsigned i1=nl->getClosePair(index).second;
+      const unsigned i0=nl->getClosePair(index).first;
+      const unsigned i1=nl->getClosePair(index).second;
       setAtomsDerivatives(val, i0, fact*deriv[i0]);
-      setAtomsDerivatives(val, i1, fact*deriv[i1]); 
-      index++;
-    }
+      setAtomsDerivatives(val, i1, fact*deriv[i1]);
+    } 
   }
-
 }
 
 }

src/colvar/RDC.cpp

@@ -275,16 +275,19 @@ void RDC::calculate()
       const double d    = distance.modulo();
       const double ind  = 1./d;
       const double id3  = ind*ind*ind; 
-      const double id7  = id3*id3*ind;
-      const double id9  = id7*ind*ind;
       const double max  = -Const*scale[index]*mu_s[index];
       const double dmax = id3*max;
       const double cos_theta = distance[2]*ind;
+
       const double rdc = 0.5*dmax*(3.*cos_theta*cos_theta-1.);
+
+      const double id7  = id3*id3*ind;
+      const double id9  = id7*ind*ind;
       const double x2=distance[0]*distance[0];
       const double y2=distance[1]*distance[1];
       const double z2=distance[2]*distance[2];
       const double prod = -max*id7*(1.5*x2 +1.5*y2 -6.*z2);
+
       Vector dRDC;
       dRDC[0] = prod*distance[0];
       dRDC[1] = prod*distance[1];