Gromacs-like replica exchange is there!

I finally upload the GREX stuff. It seems to work, but I did not really
debug it (there might be wrong units or so). I also provide a patch
for gromacs 4.5.3. Interface is still a bit tricky, and probably it
will be necessary to adapt it a bit as we try to port it to other codes.

The overall flow is:
* at steps corresponding to attempts, positions/total energy are
  saved on a temporary buffer when normal plumed is called.
* when repl_ex is called in gromacs (which, in the case of gromacs
  is after moving atoms, so with different positions) and
  the partner replica is decided by gromacs, a fake exchange of all coordinates
  is done inside plumed and a new calculation is also done.
  Notice that since presently all "analysis" tools are acting in the
  calculate() method, they are called once more at this stage. This
  should be fixed.
* Then, change in local bias potential and foreign bias potential
  are collected by gromacs and added to the acceptance.

Notice that with this approach the plumed.dat files used for different
replicas could be completely different (but for UNITS! this probably
can be fixed, even if it seems a bit stupid to me), and, a part from
trivial errors which could be there, the acceptance
is evaluated consistently

patches/gromacs-4.5.3.diff

+patch -u -l -b -F 5 --suffix=.preplumed "./src/kernel/repl_ex.c" << \EOF_EOF
+--- ./src/kernel/repl_ex.c.preplumed
++++ ./src/kernel/repl_ex.c
+@@ -49,10 +49,16 @@
+ #include "names.h"
+ #include "mvdata.h"
+ #include "domdec.h"
+ #include "partdec.h"
+ 
++/* PLUMED */
++#include "../../Plumed.h"
++extern int    plumedswitch;
++extern plumed plumedmain;
++/* END PLUMED */
++
+ typedef struct gmx_repl_ex {
+   int  repl;
+   int  nrepl;
+   real temp;
+   int  type;
+@@ -538,10 +544,25 @@
+   snew(bEx,re->nrepl);
+   snew(prob,re->nrepl);
+ 
+   exchange = -1;
+   m = (step / re->nst) % 2;
++
++  if(plumedswitch){
++    int partner=re->repl;
++
++    for(i=1; i<re->nrepl; i++) {
++      if (i % 2 != m) continue;
++      a = re->ind[i-1];
++      b = re->ind[i];
++      if(re->repl==a) partner=b;
++      if(re->repl==b) partner=a;
++    }
++    plumed_cmd(plumedmain,"GREX setPartner",&partner);
++    plumed_cmd(plumedmain,"GREX calculate",NULL);
++  }
++
+   for(i=1; i<re->nrepl; i++) {
+     a = re->ind[i-1];
+     b = re->ind[i];
+     bPrint = (re->repl==a || re->repl==b);
+     if (i % 2 == m) {
+@@ -572,10 +593,19 @@
+ 	dpV = (betaA*re->pres[a]-betaB*re->pres[b])*(Vol[b]-Vol[a])/PRESFAC;
+ 	if (bPrint)
+ 	  fprintf(fplog,"  dpV = %10.3e  d = %10.3e",dpV,delta + dpV);
+ 	delta += dpV;
+       }
++      if(plumedswitch){
++        real ldb,fdb,dplumed;
++        plumed_cmd(plumedmain,"GREX getLocalDeltaBias",&ldb);
++        plumed_cmd(plumedmain,"GREX getForeignDeltaBias",&fdb);
++        dplumed=ldb*betaA+fdb*betaB;
++        delta+=dplumed;
++	if (bPrint)
++	  fprintf(fplog,"  dplumed = %10.3e  d = %10.3e",dplumed,delta);
++      }
+       if (bPrint)
+ 	fprintf(fplog,"\n");
+       if (delta <= 0) {
+ 	prob[i] = 1;
+ 	bEx[i] = TRUE;
+@@ -633,10 +663,12 @@
+     gmx_multisim_t *ms;
+     int  exchange=-1,shift;
+     gmx_bool bExchanged=FALSE;
+     
+     ms = cr->ms;
++
++    if(plumedswitch)plumed_cmd(plumedmain,"GREX prepare",NULL);
+  
+     if (MASTER(cr))
+     {
+         exchange = get_replica_exchange(fplog,ms,re,ener,det(state->box),
+                                         step,time);
+EOF_EOF
 patch -u -l -b -F 5 --suffix=.preplumed "./src/kernel/Makefile" << \EOF_EOF
 --- ./src/kernel/Makefile.preplumed
 +++ ./src/kernel/Makefile
@@ -47,7 +127,46 @@ patch -u -l -b -F 5 --suffix=.preplumed "./src/kernel/md.c" << \EOF_EOF
  #endif
  #ifdef GMX_THREADS
  #include "tmpi.h"
-@@ -1522,10 +1528,41 @@
+@@ -1375,13 +1381,37 @@
+         for(i=0; (i<gnx); i++) {
+             grpindex[i] = i;
+         }
+     }
+ 
+-    if (repl_ex_nst > 0 && MASTER(cr))
++    if (repl_ex_nst > 0 && MASTER(cr)){
+         repl_ex = init_replica_exchange(fplog,cr->ms,state_global,ir,
+                                         repl_ex_nst,repl_ex_seed);
++        if(plumedswitch) {
++          plumed_cmd(plumedmain,"GREX setMPIIntercomm",&cr->ms->mpi_comm_masters);
++          if(PAR(cr)){
++            if(DOMAINDECOMP(cr)) {
++              plumed_cmd(plumedmain,"GREX setMPIIntracomm",&cr->dd->mpi_comm_all);
++            }else{
++              plumed_cmd(plumedmain,"GREX setMPIIntracomm",&cr->mpi_comm_mysim);
++            }
++          }
++          plumed_cmd(plumedmain,"GREX init",NULL);
++        }
++    }
++
++    if (repl_ex_nst>0 && plumedswitch && ! MASTER(cr)){
++        if(PAR(cr)){
++          if(DOMAINDECOMP(cr)) {
++            plumed_cmd(plumedmain,"GREX setMPIIntracomm",&cr->dd->mpi_comm_all);
++          }else{
++            plumed_cmd(plumedmain,"GREX setMPIIntracomm",&cr->mpi_comm_mysim);
++          }
++        }
++        plumed_cmd(plumedmain,"GREX init",NULL);
++    }
++
+ 
+     if (!ir->bContinuation && !bRerunMD)
+     {
+         if (mdatoms->cFREEZE && (state->flags & (1<<estV)))
+         {
+@@ -1522,10 +1552,41 @@
              }
          }
          fprintf(fplog,"\n");
@@ -89,7 +208,7 @@ patch -u -l -b -F 5 --suffix=.preplumed "./src/kernel/md.c" << \EOF_EOF
      print_date_and_time(fplog,cr->nodeid,"Started mdrun",runtime);
      wallcycle_start(wcycle,ewcRUN);
      if (fplog)
-@@ -1802,10 +1839,17 @@
+@@ -1802,10 +1863,17 @@
                                      state,&f,mdatoms,top,fr,
                                      vsite,shellfc,constr,
                                      nrnb,wcycle,do_verbose);
@@ -107,7 +226,7 @@ patch -u -l -b -F 5 --suffix=.preplumed "./src/kernel/md.c" << \EOF_EOF
  
          if (MASTER(cr) && do_log && !bFFscan)
          {
-@@ -1928,16 +1972,37 @@
+@@ -1928,16 +1996,40 @@
               * in do_force.
               * This is parallellized as well, and does communication too. 
               * Check comments in sim_util.c
@@ -121,6 +240,7 @@ patch -u -l -b -F 5 --suffix=.preplumed "./src/kernel/md.c" << \EOF_EOF
 +         plumed_cmd(plumedmain,"setMasses",&mdatoms->massT[mdatoms->start]);
 +         plumed_cmd(plumedmain,"setCharges",&mdatoms->chargeA[mdatoms->start]);
 +         plumed_cmd(plumedmain,"setBox",&state->box[0][0]);
++
 +         plumed_cmd(plumedmain,"prepareCalc",NULL);
 +       }
 + /* END PLUMED */
@@ -137,6 +257,8 @@ patch -u -l -b -F 5 --suffix=.preplumed "./src/kernel/md.c" << \EOF_EOF
 +         plumed_cmd(plumedmain,"setForces",&f[mdatoms->start][0]);
 +         plumed_cmd(plumedmain,"setVirial",&force_vir[0][0]);
 +         plumed_cmd(plumedmain,"performCalc",NULL);
++         if ((repl_ex_nst > 0) && (step > 0) && !bLastStep &&
++            do_per_step(step,repl_ex_nst)) plumed_cmd(plumedmain,"GREX savePositions",NULL);
 +       }
 + /* END PLUMED */
 +

src/GREX.cpp

+#include "GREX.h"
+#include "sstream"
+#include "PlumedMain.h"
+
+using namespace std;
+using namespace PLMD;
+
+GREX::GREX(PlumedMain&p):
+  initialized(false),
+  plumedMain(p),
+  atoms(p.getAtoms()),
+  partner(-1), // = unset
+  myreplica(-1) // = unset
+{
+  p.setSuffix(".NA");
+}
+
+GREX::~GREX(){
+}
+
+void GREX::cmd(const string&key,void*val){
+  if(false){
+  }else if(key=="initialized"){
+    *static_cast<int*>(val)=initialized;
+  }else if(key=="setMPIIntracomm"){
+    assert(!initialized);
+    intracomm.Set_comm(val);
+  }else if(key=="setMPIIntercomm"){
+    assert(!initialized);
+    intercomm.Set_comm(val);
+  }else if(key=="init"){
+    assert(!initialized);
+    initialized=true;
+    std::string s;
+// note that for PEs!=root this is automatically 0 (comm defaults to MPI_COMM_SELF)
+    myreplica=intercomm.Get_rank();
+    intracomm.Sum(&myreplica,1);
+    Tools::convert(myreplica,s);
+    plumedMain.setSuffix("."+s);
+  }else if(key=="prepare"){
+    assert(initialized);
+    if(intracomm.Get_rank()==0) return;
+    intracomm.Bcast(&partner,1,0);
+    calculate();
+  }else if(key=="setPartner"){
+    assert(initialized);
+    partner=*static_cast<int*>(val);
+  }else if(key=="savePositions"){
+    assert(initialized);
+    savePositions();
+  }else if(key=="calculate"){
+    assert(initialized);
+    if(intracomm.Get_rank()!=0) return;
+    intracomm.Bcast(&partner,1,0);
+    calculate();
+  }else if(key=="getLocalDeltaBias"){
+    assert(initialized);
+    double x=localDeltaBias/(atoms.getMDUnits().energy/atoms.getUnits().energy);
+    atoms.double2MD(x,val);
+  }else if(key=="getForeignDeltaBias"){
+    assert(initialized);
+    double x=foreignDeltaBias/(atoms.getMDUnits().energy/atoms.getUnits().energy);
+    atoms.double2MD(x,val);
+  } else {
+    fprintf(stderr,"+++ PLUMED GREX ERROR\n");
+    fprintf(stderr,"+++ CANNOT INTERPRET CALL TO cmd() ROUTINE WITH ARG %s\n",key.c_str());
+    fprintf(stderr,"+++ There might be a mistake in the MD code\n");
+    fprintf(stderr,"+++ or you may be using an out-dated plumed version\n");
+    assert(0);
+  }
+}
+
+void GREX::savePositions(){
+  plumedMain.prepareDependencies();
+  atoms.shareAll();
+  plumedMain.waitData();
+  ostringstream o;
+  atoms.writeBinary(o);
+  buffer=o.str();
+}
+
+void GREX::calculate(){
+//fprintf(stderr,"CALCULATE %d %d\n",intercomm.Get_rank(),partner);
+  unsigned nn=buffer.size();
+  vector<char> rbuf(nn);
+  localDeltaBias=-plumedMain.getBias();
+  if(intracomm.Get_rank()==0){
+    PlumedCommunicator::Request req=intercomm.Isend(&buffer.c_str()[0],nn,partner,1066);
+    intercomm.Recv(&rbuf[0],rbuf.size(),partner,1066);
+    req.wait();
+  }
+  intracomm.Bcast(&rbuf[0],nn,0);
+  istringstream i(string(&rbuf[0],rbuf.size()));
+  atoms.readBinary(i);
+  plumedMain.prepareDependencies();
+  plumedMain.justCalculate();
+  localDeltaBias+=plumedMain.getBias();
+  if(intracomm.Get_rank()==0){
+    PlumedCommunicator::Request req=intercomm.Isend(&localDeltaBias,1,partner,1067);
+    intercomm.Recv(&foreignDeltaBias,1,partner,1067);
+    req.wait();
+//fprintf(stderr,">>> %d %d %20.12f %20.12f\n",intercomm.Get_rank(),partner,localDeltaBias,foreignDeltaBias);
+  }
+  intracomm.Bcast(&foreignDeltaBias,1,0);
+}

src/GREX.h

+#ifndef __PLUMED_Grex_h
+#define __PLUMED_Grex_h
+
+#include "WithCmd.h"
+#include "PlumedCommunicator.h"
+#include <string>
+
+namespace PLMD{
+
+class PlumedMain;
+class Atoms;
+
+class GREX:
+  public WithCmd
+{
+  bool initialized;
+  PlumedCommunicator intracomm;
+  PlumedCommunicator intercomm;
+  PlumedMain& plumedMain;
+  Atoms&      atoms;
+  int partner;
+  double localDeltaBias;
+  double foreignDeltaBias;
+  std::string buffer;
+  int myreplica;
+public:
+  GREX(PlumedMain&);
+  ~GREX();
+  void cmd(const std::string&key,void*val=NULL);
+  void calculate();
+  void savePositions();
+};
+
+}
+
+#endif

src/PlumedMain.cpp

@@ -14,6 +14,8 @@
 
 #include "ActionRegister.h"
 
+#include "GREX.h"
+
 using namespace PLMD;
 using namespace std;
 
@@ -265,6 +267,7 @@ void PlumedMain::cmd(const std::string & word,void*val){
        if(actionRegister().check(words[1])) check=1;
        *(static_cast<int*>(val))=check;
      } else if(nw==2 && words[0]=="GREX"){
+       if(!grex) grex=new GREX(*this);
        assert(grex);
        grex->cmd(words[1],val);
      } else{