From ccd3ef1af844567cae97ef8779511a7f8f83961e Mon Sep 17 00:00:00 2001
From: Giovanni Bussi <giovanni.bussi@gmail.com>
Date: Thu, 26 Sep 2013 10:21:34 +0200
Subject: [PATCH] Merged patch into gromacs-4.6.3
As of this commit, to patch with Hamiltonian replica exchange
just use the gromacs-4.6.3 patch
This will simplify keeping up to date with upstream changes in the patch.
---
patches/gromacs-4.6.3-hrex.config | 29 -
patches/gromacs-4.6.3-hrex.diff | 585 ------------------
patches/gromacs-4.6.3.config | 23 +
patches/gromacs-4.6.3.diff/src/kernel/md.c | 91 +++
.../gromacs-4.6.3.diff/src/kernel/repl_ex.c | 20 +-
.../gromacs-4.6.3.diff/src/kernel/repl_ex.h | 82 +++
.../src/kernel/repl_ex.h.preplumed | 76 +++
7 files changed, 290 insertions(+), 616 deletions(-)
delete mode 100644 patches/gromacs-4.6.3-hrex.config
delete mode 100644 patches/gromacs-4.6.3-hrex.diff
create mode 100644 patches/gromacs-4.6.3.diff/src/kernel/repl_ex.h
create mode 100644 patches/gromacs-4.6.3.diff/src/kernel/repl_ex.h.preplumed
diff --git a/patches/gromacs-4.6.3-hrex.config b/patches/gromacs-4.6.3-hrex.config
deleted file mode 100644
index 6346748e7..000000000
--- a/patches/gromacs-4.6.3-hrex.config
+++ /dev/null
@@ -1,29 +0,0 @@
-
-
-function plumed_preliminary_test(){
-# check if the README contains the word GROMACS and if gromacs has been already configured
- grep -q GROMACS README 1>/dev/null 2>/dev/null
-}
-
-plumed_before_patch(){
-cat << EOF
-This is an experimental patch to perform Hamiltonian replica exchange with GROMACS
-Implementation described in:
- Bussi, Hamiltonian replica-exchange in GROMACS: a flexible implementation,
- Mol. Phys., DOI: 10.1080/00268976.2013.824126 (2013)
-
-To use it:
-* Prepare separate topologies (topol0.tpr, topol1.tpr, etc)
-* Set environment variable PLUMED_HREX (export PLUMED_HREX=1)
-* Run a normal replica exchange with gromacs
-
-Suggested checks:
-* Try with several identical force fields and different seed/starting point
- Acceptance should be 1.0
-
-Warnings:
-* Topologies should have the same number of atoms, same masses and same constraint topology
-* Choose neighbor list update that divides replex
-
-EOF
-}
diff --git a/patches/gromacs-4.6.3-hrex.diff b/patches/gromacs-4.6.3-hrex.diff
deleted file mode 100644
index 529123143..000000000
--- a/patches/gromacs-4.6.3-hrex.diff
+++ /dev/null
@@ -1,585 +0,0 @@
-patch -u -l -b -F 5 --suffix=.preplumed "./src/kernel/CMakeLists.txt" << \EOF_EOF
---- ./src/kernel/CMakeLists.txt.preplumed
-+++ ./src/kernel/CMakeLists.txt
-@@ -31,10 +31,12 @@
- #
- # To help us fund GROMACS development, we humbly ask that you cite
- # the research papers on the package. Check out http://www.gromacs.org.
- #
-
-+include(${CMAKE_SOURCE_DIR}/Plumed.cmake)
-+
- set(GMXPREPROCESS_SOURCES
- add_par.c
- calc_verletbuf.c
- compute_io.c
- convparm.c
-@@ -121,11 +123,11 @@
- set_target_properties(${PROGRAM} PROPERTIES OUTPUT_NAME "${PROGRAM}${GMX_BINARY_SUFFIX}")
- endforeach()
-
- add_executable(mdrun ${MDRUN_SOURCES} main.c)
- gmx_add_man_page(mdrun)
--target_link_libraries(mdrun gmxpreprocess md gmx ${OpenMP_LINKER_FLAGS})
-+target_link_libraries(mdrun gmxpreprocess md gmx ${OpenMP_LINKER_FLAGS} ${PLUMED_LOAD})
- set_target_properties(mdrun PROPERTIES OUTPUT_NAME "mdrun${GMX_BINARY_SUFFIX}" COMPILE_FLAGS "${OpenMP_C_FLAGS}")
-
- if(GMX_OPENMM)
- target_link_libraries(mdrun openmm_api_wrapper)
- endif()
-EOF_EOF
-patch -u -l -b -F 5 --suffix=.preplumed "./src/kernel/md.c" << \EOF_EOF
---- ./src/kernel/md.c.preplumed
-+++ ./src/kernel/md.c
-@@ -91,10 +91,16 @@
- #include "membed.h"
- #include "types/nlistheuristics.h"
- #include "types/iteratedconstraints.h"
- #include "nbnxn_cuda_data_mgmt.h"
-
-+/* PLUMED */
-+#include "../../Plumed.h"
-+extern int plumedswitch;
-+extern plumed plumedmain;
-+/* END PLUMED */
-+
- #ifdef GMX_LIB_MPI
- #include <mpi.h>
- #endif
- #ifdef GMX_THREAD_MPI
- #include "tmpi.h"
-@@ -235,10 +241,14 @@
- /* PME load balancing data for GPU kernels */
- pme_load_balancing_t pme_loadbal = NULL;
- double cycles_pmes;
- gmx_bool bPMETuneTry = FALSE, bPMETuneRunning = FALSE;
-
-+/* PLUMED */
-+ int plumedNeedsEnergy;
-+/* END PLUMED */
-+
- #ifdef GMX_FAHCORE
- /* Temporary addition for FAHCORE checkpointing */
- int chkpt_ret;
- #endif
-
-@@ -715,10 +725,50 @@
- }
- }
- fprintf(fplog, "\n");
- }
-
-+ /* PLUMED */
-+ if(plumedswitch){
-+ if(cr->ms && cr->ms->nsim>1) {
-+ if(MASTER(cr)) 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(PAR(cr)){
-+ if(DOMAINDECOMP(cr)) {
-+ plumed_cmd(plumedmain,"setMPIComm",&cr->dd->mpi_comm_all);
-+ }else{
-+ plumed_cmd(plumedmain,"setMPIComm",&cr->mpi_comm_mysim);
-+ }
-+ }
-+ plumed_cmd(plumedmain,"setNatoms",&top_global->natoms);
-+ plumed_cmd(plumedmain,"setMDEngine","gromacs");
-+ plumed_cmd(plumedmain,"setLog",fplog);
-+ real real_delta_t;
-+ real_delta_t=ir->delta_t;
-+ plumed_cmd(plumedmain,"setTimestep",&real_delta_t);
-+ plumed_cmd(plumedmain,"init",NULL);
-+
-+ if(PAR(cr)){
-+ if(DOMAINDECOMP(cr)) {
-+ plumed_cmd(plumedmain,"setAtomsNlocal",&cr->dd->nat_home);
-+ plumed_cmd(plumedmain,"setAtomsGatindex",cr->dd->gatindex);
-+ }else{
-+ plumed_cmd(plumedmain,"setAtomsNlocal",&mdatoms->homenr);
-+ plumed_cmd(plumedmain,"setAtomsContiguous",&mdatoms->start);
-+ }
-+ }
-+ }
-+ /* END PLUMED */
-+
- /* Set and write start time */
- runtime_start(runtime);
- print_date_and_time(fplog, cr->nodeid, "Started mdrun", runtime);
- wallcycle_start(wcycle, ewcRUN);
- if (fplog)
-@@ -1031,10 +1081,17 @@
- vsite, shellfc, constr,
- nrnb, wcycle,
- do_verbose && !bPMETuneRunning);
- wallcycle_stop(wcycle, ewcDOMDEC);
- /* If using an iterative integrator, reallocate space to match the decomposition */
-+
-+ /* PLUMED */
-+ if(plumedswitch){
-+ plumed_cmd(plumedmain,"setAtomsNlocal",&cr->dd->nat_home);
-+ plumed_cmd(plumedmain,"setAtomsGatindex",cr->dd->gatindex);
-+ }
-+ /* END PLUMED */
- }
- }
-
- if (MASTER(cr) && do_log && !bFFscan)
- {
-@@ -1080,10 +1137,96 @@
- }
- }
-
- GMX_MPE_LOG(ev_timestep2);
-
-+ gmx_bool bHREX;
-+ bHREX= repl_ex_nst > 0 && (step>0) && !bLastStep && do_per_step(step,repl_ex_nst)
-+ && getenv("PLUMED_HREX");
-+
-+/* Hamiltonian Replica Exchange */
-+ if(plumedswitch) if(bHREX){
-+ gmx_enerdata_t *hrex_enerd;
-+ snew(hrex_enerd,1);
-+ init_enerdata(top_global->groups.grps[egcENER].nr,ir->fepvals->n_lambda,hrex_enerd);
-+ int repl=-1;
-+ int nrepl=-1;
-+ if(MASTER(cr)){
-+ repl=replica_exchange_get_repl(repl_ex);
-+ nrepl=replica_exchange_get_nrepl(repl_ex);
-+ }
-+ if (PAR(cr)) {
-+ if (DOMAINDECOMP(cr))
-+ dd_collect_state(cr->dd,state,state_global);
-+ else
-+ pd_collect_state(cr,state_global);
-+ }
-+ if(MASTER(cr)){
-+ if(repl%2==step/repl_ex_nst%2){
-+ if(repl-1>=0) exchange_state(cr->ms,repl-1,state_global);
-+ }else{
-+ if(repl+1<nrepl) exchange_state(cr->ms,repl+1,state_global);
-+ }
-+ }
-+ if(PAR(cr)){
-+ if (DOMAINDECOMP(cr)) {
-+ dd_partition_system(fplog,step,cr,TRUE,1,
-+ state_global,top_global,ir,
-+ state,&f,mdatoms,top,fr,vsite,shellfc,constr,
-+ nrnb,wcycle,FALSE);
-+ } else {
-+ bcast_state(cr,state,FALSE);
-+ }
-+ }
-+ do_force(fplog, cr, ir, step, nrnb, wcycle, top, top_global, groups,
-+ state->box, state->x, &state->hist,
-+ f, force_vir, mdatoms, hrex_enerd, fcd,
-+ state->lambda, graph,
-+ fr, vsite, mu_tot, t, outf->fp_field, ed, bBornRadii,
-+ GMX_FORCE_STATECHANGED |
-+ ((DYNAMIC_BOX(*ir) || bRerunMD) ? GMX_FORCE_DYNAMICBOX : 0) |
-+ GMX_FORCE_ALLFORCES |
-+ GMX_FORCE_SEPLRF |
-+ GMX_FORCE_VIRIAL |
-+ GMX_FORCE_ENERGY |
-+ GMX_FORCE_DHDL |
-+ GMX_FORCE_NS);
-+ plumed_cmd(plumedmain,"GREX cacheLocalUSwap",&hrex_enerd->term[F_EPOT]);
-+ sfree(hrex_enerd);
-+
-+/* exchange back */
-+ if (PAR(cr)) {
-+ if (DOMAINDECOMP(cr))
-+ dd_collect_state(cr->dd,state,state_global);
-+ else
-+ pd_collect_state(cr,state_global);
-+ }
-+ if(MASTER(cr)){
-+ if(repl%2==step/repl_ex_nst%2){
-+ if(repl-1>=0) exchange_state(cr->ms,repl-1,state_global);
-+ }else{
-+ if(repl+1<nrepl) exchange_state(cr->ms,repl+1,state_global);
-+ }
-+ }
-+ if(PAR(cr)){
-+ if (DOMAINDECOMP(cr)) {
-+ dd_partition_system(fplog,step,cr,TRUE,1,
-+ state_global,top_global,ir,
-+ state,&f,mdatoms,top,fr,vsite,shellfc,constr,
-+ nrnb,wcycle,FALSE);
-+ if(plumedswitch){
-+ plumed_cmd(plumedmain,"setAtomsNlocal",&cr->dd->nat_home);
-+ plumed_cmd(plumedmain,"setAtomsGatindex",cr->dd->gatindex);
-+ }
-+ } else {
-+ bcast_state(cr,state,FALSE);
-+ }
-+ }
-+
-+ }
-+/* END Hamiltonian Replica Exchange */
-+
- /* We write a checkpoint at this MD step when:
- * either at an NS step when we signalled through gs,
- * or at the last step (but not when we do not want confout),
- * but never at the first step or with rerun.
- */
-@@ -1176,16 +1319,43 @@
- /* The coordinates (x) are shifted (to get whole molecules)
- * in do_force.
- * This is parallellized as well, and does communication too.
- * Check comments in sim_util.c
- */
-+
-+ /* PLUMED */
-+ if(plumedswitch){
-+ plumedNeedsEnergy=0;
-+ long int lstep=step; plumed_cmd(plumedmain,"setStepLong",&lstep);
-+ plumed_cmd(plumedmain,"setPositions",&state->x[mdatoms->start][0]);
-+ 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);
-+ plumed_cmd(plumedmain,"isEnergyNeeded",&plumedNeedsEnergy);
-+ }
-+ /* END PLUMED */
- do_force(fplog, cr, ir, step, nrnb, wcycle, top, top_global, groups,
- state->box, state->x, &state->hist,
- f, force_vir, mdatoms, enerd, fcd,
- state->lambda, graph,
- fr, vsite, mu_tot, t, outf->fp_field, ed, bBornRadii,
-- (bNS ? GMX_FORCE_NS : 0) | force_flags);
-+ (plumedNeedsEnergy? GMX_FORCE_ENERGY : 0) | (bNS ? GMX_FORCE_NS : 0) | force_flags);
-+ /* PLUMED */
-+ if(plumedswitch){
-+ if(plumedNeedsEnergy) plumed_cmd(plumedmain,"setEnergy",&enerd->term[F_EPOT]);
-+ 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);
-+
-+ if(bHREX)
-+ plumed_cmd(plumedmain,"GREX cacheLocalUNow",&enerd->term[F_EPOT]);
-+ }
-+ /* END PLUMED */
-+
- }
-
- GMX_BARRIER(cr->mpi_comm_mygroup);
-
- if (bTCR)
-EOF_EOF
-patch -u -l -b -F 5 --suffix=.preplumed "./src/kernel/mdrun.c" << \EOF_EOF
---- ./src/kernel/mdrun.c.preplumed
-+++ ./src/kernel/mdrun.c
-@@ -56,10 +56,16 @@
- #endif
-
- /* afm stuf */
- #include "pull.h"
-
-+/* PLUMED */
-+#include "../../Plumed.h"
-+int plumedswitch;
-+plumed plumedmain;
-+/* END PLUMED */
-+
- int cmain(int argc, char *argv[])
- {
- const char *desc[] = {
- "The [TT]mdrun[tt] program is the main computational chemistry engine",
- "within GROMACS. Obviously, it performs Molecular Dynamics simulations,",
-@@ -399,10 +405,11 @@
- { efLOG, "-rs", "rotslabs", ffOPTWR },
- { efLOG, "-rt", "rottorque", ffOPTWR },
- { efMTX, "-mtx", "nm", ffOPTWR },
- { efNDX, "-dn", "dipole", ffOPTWR },
- { efRND, "-multidir", NULL, ffOPTRDMULT},
-+ { efDAT, "-plumed", "plumed", ffOPTRD }, /* PLUMED */
- { efDAT, "-membed", "membed", ffOPTRD },
- { efTOP, "-mp", "membed", ffOPTRD },
- { efNDX, "-mn", "membed", ffOPTRD }
- };
- #define NFILE asize(fnm)
-@@ -731,19 +738,50 @@
- }
-
- ddxyz[XX] = (int)(realddxyz[XX] + 0.5);
- ddxyz[YY] = (int)(realddxyz[YY] + 0.5);
- ddxyz[ZZ] = (int)(realddxyz[ZZ] + 0.5);
-+ /* PLUMED */
-+ plumedswitch=0;
-+ if (opt2bSet("-plumed",NFILE,fnm)) plumedswitch=1;
-+ if(plumedswitch){
-+ int plumed_is_there=0;
-+ int real_precision=sizeof(real);
-+ real energyUnits=1.0;
-+ real lengthUnits=1.0;
-+ real timeUnits=1.0;
-+
-+
-+ if(!plumed_installed()){
-+ gmx_fatal(FARGS,"Plumed is not available. Check your PLUMED_KERNEL variable.");
-+ }
-+ plumedmain=plumed_create();
-+ plumed_cmd(plumedmain,"setRealPrecision",&real_precision);
-+ // this is not necessary for gromacs units:
-+ plumed_cmd(plumedmain,"setMDEnergyUnits",&energyUnits);
-+ plumed_cmd(plumedmain,"setMDLengthUnits",&lengthUnits);
-+ plumed_cmd(plumedmain,"setMDTimeUnits",&timeUnits);
-+ //
-+ plumed_cmd(plumedmain,"setPlumedDat",ftp2fn(efDAT,NFILE,fnm));
-+ plumedswitch=1;
-+ }
-+ /* END PLUMED */
-
- rc = mdrunner(&hw_opt, fplog, cr, NFILE, fnm, oenv, bVerbose, bCompact,
- nstglobalcomm, ddxyz, dd_node_order, rdd, rconstr,
- dddlb_opt[0], dlb_scale, ddcsx, ddcsy, ddcsz,
- nbpu_opt[0],
- nsteps, nstepout, resetstep,
- nmultisim, repl_ex_nst, repl_ex_nex, repl_ex_seed,
- pforce, cpt_period, max_hours, deviceOptions, Flags);
-
-+ /* PLUMED */
-+ if(plumedswitch){
-+ plumed_finalize(plumedmain);
-+ }
-+ /* END PLUMED */
-+
- gmx_finalize_par();
-
- if (MULTIMASTER(cr))
- {
- thanx(stderr);
-EOF_EOF
-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
-@@ -51,10 +51,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 */
-+
- #define PROBABILITYCUTOFF 100
- /* we don't bother evaluating if events are more rare than exp(-100) = 3.7x10^-44 */
-
- enum {
- ereTEMP, ereLAMBDA, ereENDSINGLE, ereTL, ereNR
-@@ -111,18 +117,20 @@
-
- snew(qall, ms->nsim);
- qall[re->repl] = q;
- gmx_sum_sim(ms->nsim, qall, ms);
-
-- bDiff = FALSE;
-- for (s = 1; s < ms->nsim; s++)
-- {
-- if (qall[s] != qall[0])
-- {
-+ /* PLUMED */
-+ //bDiff = FALSE;
-+ //for (s = 1; s < ms->nsim; s++)
-+ //{
-+ // if (qall[s] != qall[0])
-+ // {
- bDiff = TRUE;
-- }
-- }
-+ // }
-+ //}
-+ /* END PLUMED */
-
- if (bDiff)
- {
- /* Set the replica exchange type and quantities */
- re->type = ere;
-@@ -255,10 +263,14 @@
- for (i = 0; i < re->nrepl; i++)
- {
- re->ind[i] = i;
- }
-
-+ /* PLUMED */
-+ // plumed2: check if we want alternative patterns (i.e. for bias-exchange metaD)
-+ // in those cases replicas can share the same temperature.
-+ /*
- if (re->type < ereENDSINGLE)
- {
-
- for (i = 0; i < re->nrepl; i++)
- {
-@@ -275,10 +287,12 @@
- gmx_fatal(FARGS, "Two replicas have identical %ss", erename[re->type]);
- }
- }
- }
- }
-+ */
-+ /* END PLUMED */
-
- /* keep track of all the swaps, starting with the initial placement. */
- snew(re->allswaps, re->nrepl);
- for (i = 0; i < re->nrepl; i++)
- {
-@@ -524,11 +538,11 @@
- }
- sfree(buf);
- }
- }
-
--static void exchange_state(const gmx_multisim_t *ms, int b, t_state *state)
-+void exchange_state(const gmx_multisim_t *ms, int b, t_state *state)
- {
- /* When t_state changes, this code should be updated. */
- int ngtc, nnhpres;
- ngtc = state->ngtc * state->nhchainlength;
- nnhpres = state->nnhpres* state->nhchainlength;
-@@ -644,11 +658,11 @@
- svmul(fac, state->v[i], state->v[i]);
- }
- }
- }
-
--static void pd_collect_state(const t_commrec *cr, t_state *state)
-+void pd_collect_state(const t_commrec *cr, t_state *state)
- {
- int shift;
-
- if (debug)
- {
-@@ -869,10 +883,17 @@
- }
- if (bPrint)
- {
- fprintf(fplog, "Repl %d <-> %d dE_term = %10.3e (kT)\n", a, b, delta);
- }
-+
-+/* PLUMED */
-+/* this is necessary because with plumed HREX the energy contribution is
-+ already taken into account */
-+ if(getenv("PLUMED_HREX")) delta=0.0;
-+/* END PLUMED */
-+
- if (re->bNPT)
- {
- /* revist the calculation for 5.0. Might be some improvements. */
- dpV = (beta[ap]*re->pres[ap]-beta[bp]*re->pres[bp])*(Vol[b]-Vol[a])/PRESFAC;
- if (bPrint)
-@@ -1045,19 +1066,57 @@
- }
- else
- {
- /* standard nearest neighbor replica exchange */
- m = (step / re->nst) % 2;
-+ /* PLUMED */
-+ if(plumedswitch){
-+ int partner=re->repl;
-+ int test=0;
-+ plumed_cmd(plumedmain,"getExchangesFlag",&test);
-+ if(test>0){
-+ int *list;
-+ snew(list,re->nrepl);
-+ plumed_cmd(plumedmain,"setNumberOfReplicas",&(re->nrepl));
-+ plumed_cmd(plumedmain,"getExchangesList",list);
-+ for(i=0; i<re->nrepl; i++) re->ind[i]=list[i];
-+ sfree(list);
-+ }
-+
-+ 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);
-+ plumed_cmd(plumedmain,"GREX shareAllDeltaBias",NULL);
-+ }
-+ /* END PLUMED */
- 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)
- {
- delta = calc_delta(fplog, bPrint, re, a, b, a, b);
-+ /* PLUMED */
-+ if(plumedswitch){
-+ real adb,bdb,dplumed;
-+ char buf[300];
-+ sprintf(buf,"GREX getDeltaBias %d",a); plumed_cmd(plumedmain,buf,&adb);
-+ sprintf(buf,"GREX getDeltaBias %d",b); plumed_cmd(plumedmain,buf,&bdb);
-+ dplumed=adb*re->beta[a]+bdb*re->beta[b];
-+ delta+=dplumed;
-+ if (bPrint)
-+ fprintf(fplog,"dplumed = %10.3e dE_Term = %10.3e (kT)\n",dplumed,delta);
-+ }
-+ /* END PLUMED */
- if (delta <= 0)
- {
- /* accepted */
- prob[i] = 1;
- bEx[i] = TRUE;
-@@ -1306,10 +1365,14 @@
- * exchanges. */
- /* Where each replica ends up after the exchange attempt(s). */
- /* The order in which multiple exchanges will occur. */
- gmx_bool bThisReplicaExchanged = FALSE;
-
-+ /* PLUMED */
-+ if(plumedswitch)plumed_cmd(plumedmain,"GREX prepare",NULL);
-+ /* END PLUMED */
-+
- if (MASTER(cr))
- {
- replica_id = re->repl;
- test_for_replica_exchange(fplog, cr->ms, re, enerd, det(state_local->box), step, time);
- prepare_to_do_exchange(fplog, re->destinations, replica_id, re->nrepl, &maxswap,
-@@ -1437,5 +1500,12 @@
- fprintf(fplog, "\n");
- }
- /* print the transition matrix */
- print_transition_matrix(fplog, "", re->nrepl, re->nmoves, re->nattempt);
- }
-+
-+int replica_exchange_get_repl(const gmx_repl_ex_t re){
-+ return re->repl;
-+};
-+int replica_exchange_get_nrepl(const gmx_repl_ex_t re){
-+ return re->nrepl;
-+};
-EOF_EOF
-patch -u -l -b -F 5 --suffix=.preplumed "./src/kernel/repl_ex.h" << \EOF_EOF
---- ./src/kernel/repl_ex.h.preplumed
-+++ ./src/kernel/repl_ex.h
-@@ -71,6 +71,12 @@
- /* Should only be called on the master nodes */
-
- extern void pd_distribute_state(const t_commrec *cr, t_state *state);
- /* Distributes the state after exchange for particle decomposition */
-
-+extern int replica_exchange_get_repl(const gmx_repl_ex_t re);
-+extern int replica_exchange_get_nrepl(const gmx_repl_ex_t re);
-+extern void pd_collect_state(const t_commrec *cr, t_state *state);
-+extern void exchange_state(const gmx_multisim_t *ms, int b, t_state *state);
-+
-+
- #endif /* _repl_ex_h */
-EOF_EOF
diff --git a/patches/gromacs-4.6.3.config b/patches/gromacs-4.6.3.config
index a5d00fcea..6346748e7 100644
--- a/patches/gromacs-4.6.3.config
+++ b/patches/gromacs-4.6.3.config
@@ -4,3 +4,26 @@ function plumed_preliminary_test(){
# check if the README contains the word GROMACS and if gromacs has been already configured
grep -q GROMACS README 1>/dev/null 2>/dev/null
}
+
+plumed_before_patch(){
+cat << EOF
+This is an experimental patch to perform Hamiltonian replica exchange with GROMACS
+Implementation described in:
+ Bussi, Hamiltonian replica-exchange in GROMACS: a flexible implementation,
+ Mol. Phys., DOI: 10.1080/00268976.2013.824126 (2013)
+
+To use it:
+* Prepare separate topologies (topol0.tpr, topol1.tpr, etc)
+* Set environment variable PLUMED_HREX (export PLUMED_HREX=1)
+* Run a normal replica exchange with gromacs
+
+Suggested checks:
+* Try with several identical force fields and different seed/starting point
+ Acceptance should be 1.0
+
+Warnings:
+* Topologies should have the same number of atoms, same masses and same constraint topology
+* Choose neighbor list update that divides replex
+
+EOF
+}
diff --git a/patches/gromacs-4.6.3.diff/src/kernel/md.c b/patches/gromacs-4.6.3.diff/src/kernel/md.c
index 187b2796a..ec5ffd289 100644
--- a/patches/gromacs-4.6.3.diff/src/kernel/md.c
+++ b/patches/gromacs-4.6.3.diff/src/kernel/md.c
@@ -1139,6 +1139,94 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
GMX_MPE_LOG(ev_timestep2);
+ gmx_bool bHREX;
+ bHREX= repl_ex_nst > 0 && (step>0) && !bLastStep && do_per_step(step,repl_ex_nst)
+ && getenv("PLUMED_HREX");
+
+/* Hamiltonian Replica Exchange */
+ if(plumedswitch) if(bHREX){
+ gmx_enerdata_t *hrex_enerd;
+ snew(hrex_enerd,1);
+ init_enerdata(top_global->groups.grps[egcENER].nr,ir->fepvals->n_lambda,hrex_enerd);
+ int repl=-1;
+ int nrepl=-1;
+ if(MASTER(cr)){
+ repl=replica_exchange_get_repl(repl_ex);
+ nrepl=replica_exchange_get_nrepl(repl_ex);
+ }
+ if (PAR(cr)) {
+ if (DOMAINDECOMP(cr))
+ dd_collect_state(cr->dd,state,state_global);
+ else
+ pd_collect_state(cr,state_global);
+ }
+ if(MASTER(cr)){
+ if(repl%2==step/repl_ex_nst%2){
+ if(repl-1>=0) exchange_state(cr->ms,repl-1,state_global);
+ }else{
+ if(repl+1<nrepl) exchange_state(cr->ms,repl+1,state_global);
+ }
+ }
+ if(PAR(cr)){
+ if (DOMAINDECOMP(cr)) {
+ dd_partition_system(fplog,step,cr,TRUE,1,
+ state_global,top_global,ir,
+ state,&f,mdatoms,top,fr,vsite,shellfc,constr,
+ nrnb,wcycle,FALSE);
+ } else {
+ bcast_state(cr,state,FALSE);
+ }
+ }
+ do_force(fplog, cr, ir, step, nrnb, wcycle, top, top_global, groups,
+ state->box, state->x, &state->hist,
+ f, force_vir, mdatoms, hrex_enerd, fcd,
+ state->lambda, graph,
+ fr, vsite, mu_tot, t, outf->fp_field, ed, bBornRadii,
+ GMX_FORCE_STATECHANGED |
+ ((DYNAMIC_BOX(*ir) || bRerunMD) ? GMX_FORCE_DYNAMICBOX : 0) |
+ GMX_FORCE_ALLFORCES |
+ GMX_FORCE_SEPLRF |
+ GMX_FORCE_VIRIAL |
+ GMX_FORCE_ENERGY |
+ GMX_FORCE_DHDL |
+ GMX_FORCE_NS);
+ plumed_cmd(plumedmain,"GREX cacheLocalUSwap",&hrex_enerd->term[F_EPOT]);
+ sfree(hrex_enerd);
+
+/* exchange back */
+ if (PAR(cr)) {
+ if (DOMAINDECOMP(cr))
+ dd_collect_state(cr->dd,state,state_global);
+ else
+ pd_collect_state(cr,state_global);
+ }
+ if(MASTER(cr)){
+ if(repl%2==step/repl_ex_nst%2){
+ if(repl-1>=0) exchange_state(cr->ms,repl-1,state_global);
+ }else{
+ if(repl+1<nrepl) exchange_state(cr->ms,repl+1,state_global);
+ }
+ }
+ if(PAR(cr)){
+ if (DOMAINDECOMP(cr)) {
+ dd_partition_system(fplog,step,cr,TRUE,1,
+ state_global,top_global,ir,
+ state,&f,mdatoms,top,fr,vsite,shellfc,constr,
+ nrnb,wcycle,FALSE);
+ if(plumedswitch){
+ plumed_cmd(plumedmain,"setAtomsNlocal",&cr->dd->nat_home);
+ plumed_cmd(plumedmain,"setAtomsGatindex",cr->dd->gatindex);
+ }
+ } else {
+ bcast_state(cr,state,FALSE);
+ }
+ }
+
+ }
+/* END Hamiltonian Replica Exchange */
+
+
+
/* We write a checkpoint at this MD step when:
* either at an NS step when we signalled through gs,
* or at the last step (but not when we do not want confout),
@@ -1262,6 +1350,9 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
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);
+ if(bHREX)
+ plumed_cmd(plumedmain,"GREX cacheLocalUNow",&enerd->term[F_EPOT]);
+
}
/* END PLUMED */
}
diff --git a/patches/gromacs-4.6.3.diff/src/kernel/repl_ex.c b/patches/gromacs-4.6.3.diff/src/kernel/repl_ex.c
index 2e79c67ad..2f5f99159 100644
--- a/patches/gromacs-4.6.3.diff/src/kernel/repl_ex.c
+++ b/patches/gromacs-4.6.3.diff/src/kernel/repl_ex.c
@@ -540,7 +540,7 @@ static void exchange_rvecs(const gmx_multisim_t *ms, int b, rvec *v, int n)
}
}
-static void exchange_state(const gmx_multisim_t *ms, int b, t_state *state)
+void exchange_state(const gmx_multisim_t *ms, int b, t_state *state)
{
/* When t_state changes, this code should be updated. */
int ngtc, nnhpres;
@@ -660,7 +660,7 @@ static void scale_velocities(t_state *state, real fac)
}
}
-static void pd_collect_state(const t_commrec *cr, t_state *state)
+void pd_collect_state(const t_commrec *cr, t_state *state)
{
int shift;
@@ -885,6 +885,13 @@ static real calc_delta(FILE *fplog, gmx_bool bPrint, struct gmx_repl_ex *re, int
{
fprintf(fplog, "Repl %d <-> %d dE_term = %10.3e (kT)\n", a, b, delta);
}
+
+/* PLUMED */
+/* this is necessary because with plumed HREX the energy contribution is
+ already taken into account */
+ if(getenv("PLUMED_HREX")) delta=0.0;
+/* END PLUMED */
+
if (re->bNPT)
{
/* revist the calculation for 5.0. Might be some improvements. */
@@ -1518,3 +1525,12 @@ void print_replica_exchange_statistics(FILE *fplog, struct gmx_repl_ex *re)
/* print the transition matrix */
print_transition_matrix(fplog, "", re->nrepl, re->nmoves, re->nattempt);
}
+
+int replica_exchange_get_repl(const gmx_repl_ex_t re){
+ return re->repl;
+};
+
+int replica_exchange_get_nrepl(const gmx_repl_ex_t re){
+ return re->nrepl;
+};
+
diff --git a/patches/gromacs-4.6.3.diff/src/kernel/repl_ex.h b/patches/gromacs-4.6.3.diff/src/kernel/repl_ex.h
new file mode 100644
index 000000000..4a4003b96
--- /dev/null
+++ b/patches/gromacs-4.6.3.diff/src/kernel/repl_ex.h
@@ -0,0 +1,82 @@
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2004, The GROMACS development team,
+ * check out http://www.gromacs.org for more information.
+ * Copyright (c) 2012,2013, by the GROMACS development team, led by
+ * David van der Spoel, Berk Hess, Erik Lindahl, and including many
+ * others, as listed in the AUTHORS file in the top-level source
+ * directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+
+#ifndef _repl_ex_h
+#define _repl_ex_h
+
+#include "typedefs.h"
+#include "types/commrec.h"
+
+/* Abstract type for replica exchange */
+typedef struct gmx_repl_ex *gmx_repl_ex_t;
+
+extern gmx_repl_ex_t init_replica_exchange(FILE *fplog,
+ const gmx_multisim_t *ms,
+ const t_state *state,
+ const t_inputrec *ir,
+ int nst, int nmultiex, int init_seed);
+/* Should only be called on the master nodes */
+
+extern gmx_bool replica_exchange(FILE *fplog,
+ const t_commrec *cr,
+ gmx_repl_ex_t re,
+ t_state *state, gmx_enerdata_t *enerd,
+ t_state *state_local,
+ gmx_large_int_t step, real time);
+/* Attempts replica exchange, should be called on all nodes.
+ * Returns TRUE if this state has been exchanged.
+ * When running each replica in parallel,
+ * this routine collects the state on the master node before exchange.
+ * With particle the state is redistributed over the nodes after exchange.
+ * With domain decomposition the global state after exchanged in stored
+ * in state and still needs to be redistributed over the nodes.
+ */
+
+extern void print_replica_exchange_statistics(FILE *fplog, gmx_repl_ex_t re);
+/* Should only be called on the master nodes */
+
+extern void pd_distribute_state(const t_commrec *cr, t_state *state);
+/* Distributes the state after exchange for particle decomposition */
+
+extern int replica_exchange_get_repl(const gmx_repl_ex_t re);
+extern int replica_exchange_get_nrepl(const gmx_repl_ex_t re);
+extern void pd_collect_state(const t_commrec *cr, t_state *state);
+extern void exchange_state(const gmx_multisim_t *ms, int b, t_state *state);
+
+
+#endif /* _repl_ex_h */
diff --git a/patches/gromacs-4.6.3.diff/src/kernel/repl_ex.h.preplumed b/patches/gromacs-4.6.3.diff/src/kernel/repl_ex.h.preplumed
new file mode 100644
index 000000000..97c45c56b
--- /dev/null
+++ b/patches/gromacs-4.6.3.diff/src/kernel/repl_ex.h.preplumed
@@ -0,0 +1,76 @@
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2004, The GROMACS development team,
+ * check out http://www.gromacs.org for more information.
+ * Copyright (c) 2012,2013, by the GROMACS development team, led by
+ * David van der Spoel, Berk Hess, Erik Lindahl, and including many
+ * others, as listed in the AUTHORS file in the top-level source
+ * directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+
+#ifndef _repl_ex_h
+#define _repl_ex_h
+
+#include "typedefs.h"
+#include "types/commrec.h"
+
+/* Abstract type for replica exchange */
+typedef struct gmx_repl_ex *gmx_repl_ex_t;
+
+extern gmx_repl_ex_t init_replica_exchange(FILE *fplog,
+ const gmx_multisim_t *ms,
+ const t_state *state,
+ const t_inputrec *ir,
+ int nst, int nmultiex, int init_seed);
+/* Should only be called on the master nodes */
+
+extern gmx_bool replica_exchange(FILE *fplog,
+ const t_commrec *cr,
+ gmx_repl_ex_t re,
+ t_state *state, gmx_enerdata_t *enerd,
+ t_state *state_local,
+ gmx_large_int_t step, real time);
+/* Attempts replica exchange, should be called on all nodes.
+ * Returns TRUE if this state has been exchanged.
+ * When running each replica in parallel,
+ * this routine collects the state on the master node before exchange.
+ * With particle the state is redistributed over the nodes after exchange.
+ * With domain decomposition the global state after exchanged in stored
+ * in state and still needs to be redistributed over the nodes.
+ */
+
+extern void print_replica_exchange_statistics(FILE *fplog, gmx_repl_ex_t re);
+/* Should only be called on the master nodes */
+
+extern void pd_distribute_state(const t_commrec *cr, t_state *state);
+/* Distributes the state after exchange for particle decomposition */
+
+#endif /* _repl_ex_h */
--
GitLab