From 5e9e14948e3c250824f24ae4c33e6f5884505cf5 Mon Sep 17 00:00:00 2001
From: Giovanni Bussi <giovanni.bussi@gmail.com>
Date: Fri, 2 Dec 2011 18:08:43 +0100
Subject: [PATCH] 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 | 128 ++++++++++++++++++++++++++++++++++++-
src/GREX.cpp | 105 ++++++++++++++++++++++++++++++
src/GREX.h | 36 +++++++++++
src/PlumedMain.cpp | 3 +
4 files changed, 269 insertions(+), 3 deletions(-)
create mode 100644 src/GREX.cpp
create mode 100644 src/GREX.h
diff --git a/patches/gromacs-4.5.3.diff b/patches/gromacs-4.5.3.diff
index f5da3ceb6..73f1014d0 100644
--- a/patches/gromacs-4.5.3.diff
+++ b/patches/gromacs-4.5.3.diff
@@ -1,3 +1,83 @@
+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 */
+
diff --git a/src/GREX.cpp b/src/GREX.cpp
new file mode 100644
index 000000000..9536ef7b7
--- /dev/null
+++ b/src/GREX.cpp
@@ -0,0 +1,105 @@
+#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);
+}
diff --git a/src/GREX.h b/src/GREX.h
new file mode 100644
index 000000000..a240029be
--- /dev/null
+++ b/src/GREX.h
@@ -0,0 +1,36 @@
+#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
diff --git a/src/PlumedMain.cpp b/src/PlumedMain.cpp
index 2fa4350c1..c944c753c 100644
--- a/src/PlumedMain.cpp
+++ b/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{
--
GitLab