From 4d74b6a5c43cc391d1b18a7703ffd26ceee612e6 Mon Sep 17 00:00:00 2001
From: carlocamilloni <carlo.camilloni@gmail.com>
Date: Thu, 9 Nov 2017 22:37:44 +0100
Subject: [PATCH] hrex: gromacs-2016.4

---
 CHANGES/v2.4.md                               |  1 +
 .../src/programs/mdrun/md.cpp                 | 86 ++++++++++++++++++
 .../src/programs/mdrun/mdrun.cpp              | 15 ++++
 .../src/programs/mdrun/repl_ex.cpp            | 29 ++++++-
 .../src/programs/mdrun/repl_ex.h              | 87 +++++++++++++++++++
 .../src/programs/mdrun/repl_ex.h.preplumed    | 79 +++++++++++++++++
 6 files changed, 295 insertions(+), 2 deletions(-)
 create mode 100644 patches/gromacs-2016.4.diff/src/programs/mdrun/repl_ex.h
 create mode 100644 patches/gromacs-2016.4.diff/src/programs/mdrun/repl_ex.h.preplumed

diff --git a/CHANGES/v2.4.md b/CHANGES/v2.4.md
index 3d2673bee..d64b99fc5 100644
--- a/CHANGES/v2.4.md
+++ b/CHANGES/v2.4.md
@@ -154,6 +154,7 @@ Fixes after beta release:
   - Fixed a bug in EDS module for restarts, where accumulated rate was not stored to restart file, leading to convergence problems.
   - (developers) `plumed config makefile_conf` can be used to retrieve `Makefile.conf` file a posteriori.
   - Added patch for Quantum ESPRESSO 6.2 (thanks to Ralf Meyer).
+  - Implemented HREX for gromacs-2016.4.
   - (developers) Store `MPIEXEC` variable at configure time and use it later for running regtests. Notice that in case
     `MPIEXEC` is not specified regtests will be run using the command stored in env var `PLUMED_MPIRUN` or, if this is
     also not defined, using `mpirun`.
diff --git a/patches/gromacs-2016.4.diff/src/programs/mdrun/md.cpp b/patches/gromacs-2016.4.diff/src/programs/mdrun/md.cpp
index f3b524385..8e808e163 100644
--- a/patches/gromacs-2016.4.diff/src/programs/mdrun/md.cpp
+++ b/patches/gromacs-2016.4.diff/src/programs/mdrun/md.cpp
@@ -127,6 +127,10 @@ extern int    plumedswitch;
 extern plumed plumedmain;
 /* END PLUMED */
 
+/* PLUMED HREX */
+extern int plumed_hrex;
+/* END PLUMED HREX */
+
 #ifdef GMX_FAHCORE
 #include "corewrap.h"
 #endif
@@ -1087,6 +1091,86 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
         }
         clear_mat(force_vir);
 
+/* PLUMED HREX */
+        gmx_bool bHREX;
+        bHREX= repl_ex_nst > 0 && (step>0) && !bLastStep && do_per_step(step,repl_ex_nst) && plumed_hrex;
+
+        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 (DOMAINDECOMP(cr))
+            dd_collect_state(cr->dd,state,state_global);
+          else
+            copy_state_nonatomdata(state, 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 (!DOMAINDECOMP(cr))
+            copy_state_nonatomdata(state_global, state);
+          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,constr,
+                              nrnb,wcycle,FALSE);
+            }
+          }
+          do_force(fplog, cr, ir, step, nrnb, wcycle, top, groups,
+                     state->box, state->x, &state->hist,
+                     f, force_vir, mdatoms, hrex_enerd, fcd,
+                     state->lambda, graph,
+                     fr, vsite, mu_tot, t, mdoutf_get_fp_field(outf), ed, bBornRadii,
+                    GMX_FORCE_STATECHANGED |
+                       ((inputrecDynamicBox(ir) || bRerunMD) ? GMX_FORCE_DYNAMICBOX : 0) |
+                       GMX_FORCE_ALLFORCES |
+                       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 (DOMAINDECOMP(cr))
+            dd_collect_state(cr->dd,state,state_global);
+          else
+            copy_state_nonatomdata(state, 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 (!DOMAINDECOMP(cr))
+            copy_state_nonatomdata(state_global, state);
+          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,constr,
+                              nrnb,wcycle,FALSE);
+              if(plumedswitch){
+                plumed_cmd(plumedmain,"setAtomsNlocal",&cr->dd->nat_home);
+                plumed_cmd(plumedmain,"setAtomsGatindex",cr->dd->gatindex);
+              }
+            }
+          }
+
+        }
+/* END PLUMED HREX */
+
         /* 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),
@@ -1203,6 +1287,8 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
               if ((repl_ex_nst > 0) && (step > 0) && !bLastStep &&
                  do_per_step(step,repl_ex_nst)) plumed_cmd(plumedmain,"GREX savePositions",NULL);
               if(plumedWantsToStop) ir->nsteps=step_rel+1;
+              if(bHREX)
+                 plumed_cmd(plumedmain,"GREX cacheLocalUNow",&enerd->term[F_EPOT]);
             }
             /* END PLUMED */
         }
diff --git a/patches/gromacs-2016.4.diff/src/programs/mdrun/mdrun.cpp b/patches/gromacs-2016.4.diff/src/programs/mdrun/mdrun.cpp
index cf0f72cf8..2b131b65d 100644
--- a/patches/gromacs-2016.4.diff/src/programs/mdrun/mdrun.cpp
+++ b/patches/gromacs-2016.4.diff/src/programs/mdrun/mdrun.cpp
@@ -78,6 +78,10 @@ extern plumed plumedmain;
 extern void(*plumedcmd)(plumed,const char*,const void*);
 /* END PLUMED */
 
+/* PLUMED HREX */
+int plumed_hrex;
+/* END PLUMED HREX */
+
 /*! \brief Return whether either of the command-line parameters that
  *  will trigger a multi-simulation is set */
 static bool is_multisim_option_set(int argc, const char *const argv[])
@@ -403,6 +407,8 @@ int gmx_mdrun(int argc, char *argv[])
           "Number of random exchanges to carry out each exchange interval (N^3 is one suggestion).  -nex zero or not specified gives neighbor replica exchange." },
         { "-reseed",  FALSE, etINT, {&repl_ex_seed},
           "Seed for replica exchange, -1 is generate a seed" },
+        { "-hrex",  FALSE, etBOOL, {&plumed_hrex},
+          "Enable hamiltonian replica exchange" },
         { "-imdport",    FALSE, etINT, {&imdport},
           "HIDDENIMD listening port" },
         { "-imdwait",  FALSE, etBOOL, {&bIMDwait},
@@ -571,6 +577,15 @@ int gmx_mdrun(int argc, char *argv[])
       plumed_cmd(plumedmain,"setPlumedDat",ftp2fn(efDAT,NFILE,fnm));
       plumedswitch=1;
     }
+    /* PLUMED HREX*/
+    if(getenv("PLUMED_HREX")) plumed_hrex=1;
+    if(plumed_hrex){
+      if(!plumedswitch)  gmx_fatal(FARGS,"-hrex (or PLUMED_HREX) requires -plumed");
+      if(repl_ex_nst==0) gmx_fatal(FARGS,"-hrex (or PLUMED_HREX) replica exchange");
+      if(repl_ex_nex!=0) gmx_fatal(FARGS,"-hrex (or PLUMED_HREX) not compatible with -nex");
+    }
+    /* END PLUMED HREX */
+
     /* END PLUMED */
 
     rc = gmx::mdrunner(&hw_opt, fplog, cr, NFILE, fnm, oenv, bVerbose,
diff --git a/patches/gromacs-2016.4.diff/src/programs/mdrun/repl_ex.cpp b/patches/gromacs-2016.4.diff/src/programs/mdrun/repl_ex.cpp
index 70d9f30b8..2f35484ec 100644
--- a/patches/gromacs-2016.4.diff/src/programs/mdrun/repl_ex.cpp
+++ b/patches/gromacs-2016.4.diff/src/programs/mdrun/repl_ex.cpp
@@ -67,6 +67,10 @@ extern int    plumedswitch;
 extern plumed plumedmain;
 /* END PLUMED */
 
+/* PLUMED HREX */
+extern int plumed_hrex;
+/* END PLUMED HREX */
+
 #define PROBABILITYCUTOFF 100
 /* we don't bother evaluating if events are more rare than exp(-100) = 3.7x10^-44 */
 
@@ -543,7 +547,9 @@ static void exchange_rvecs(const gmx_multisim_t gmx_unused *ms, int gmx_unused b
     }
 }
 
-static void exchange_state(const gmx_multisim_t *ms, int b, t_state *state)
+/* PLUMED HREX */
+void exchange_state(const gmx_multisim_t *ms, int b, t_state *state)
+/* END PLUMED HREX */
 {
     /* When t_state changes, this code should be updated. */
     int ngtc, nnhpres;
@@ -623,7 +629,9 @@ static void copy_ints(const int *s, int *d, int n)
 #define scopy_reals(v, n) copy_reals(state->v, state_local->v, n);
 #define scopy_ints(v, n)   copy_ints(state->v, state_local->v, n);
 
-static void copy_state_nonatomdata(t_state *state, t_state *state_local)
+/* PLUMED HREX */
+void copy_state_nonatomdata(t_state *state, t_state *state_local)
+/* END PLUMED HREX */
 {
     /* When t_state changes, this code should be updated. */
     int ngtc, nnhpres;
@@ -866,6 +874,11 @@ 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 HREX */
+/* this is necessary because with plumed HREX the energy contribution is
+   already taken into account */
+    if(plumed_hrex) delta=0.0;
+/* END PLUMED HREX */
     if (re->bNPT)
     {
         /* revist the calculation for 5.0.  Might be some improvements. */
@@ -1488,3 +1501,15 @@ 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);
 }
+
+/* PLUMED HREX */
+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;
+};
+/* END PLUMED HREX */
+
+
diff --git a/patches/gromacs-2016.4.diff/src/programs/mdrun/repl_ex.h b/patches/gromacs-2016.4.diff/src/programs/mdrun/repl_ex.h
new file mode 100644
index 000000000..32499012a
--- /dev/null
+++ b/patches/gromacs-2016.4.diff/src/programs/mdrun/repl_ex.h
@@ -0,0 +1,87 @@
+/*
+ * 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.
+ * Copyright (c) 2011,2012,2013,2014,2015, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and 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 <cstdio>
+
+#include "gromacs/utility/basedefinitions.h"
+#include "gromacs/utility/real.h"
+
+struct gmx_enerdata_t;
+struct gmx_multisim_t;
+struct t_commrec;
+struct t_inputrec;
+struct t_state;
+
+/* Abstract type for replica exchange */
+typedef struct gmx_repl_ex *gmx_repl_ex_t;
+
+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 */
+
+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_int64_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 domain decomposition, the global state after exchange is stored
+ * in state and still needs to be redistributed over the nodes.
+ */
+
+void print_replica_exchange_statistics(FILE *fplog, gmx_repl_ex_t re);
+/* Should only be called on the master nodes */
+
+/* PLUMED HREX */
+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);
+extern void copy_state_nonatomdata(t_state *state, t_state *state_local);
+/* END PLUMED HREX */
+
+#endif  /* _repl_ex_h */
diff --git a/patches/gromacs-2016.4.diff/src/programs/mdrun/repl_ex.h.preplumed b/patches/gromacs-2016.4.diff/src/programs/mdrun/repl_ex.h.preplumed
new file mode 100644
index 000000000..dcb47d1fd
--- /dev/null
+++ b/patches/gromacs-2016.4.diff/src/programs/mdrun/repl_ex.h.preplumed
@@ -0,0 +1,79 @@
+/*
+ * 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.
+ * Copyright (c) 2011,2012,2013,2014,2015, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and 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 <cstdio>
+
+#include "gromacs/utility/basedefinitions.h"
+#include "gromacs/utility/real.h"
+
+struct gmx_enerdata_t;
+struct gmx_multisim_t;
+struct t_commrec;
+struct t_inputrec;
+struct t_state;
+
+/* Abstract type for replica exchange */
+typedef struct gmx_repl_ex *gmx_repl_ex_t;
+
+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 */
+
+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_int64_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 domain decomposition, the global state after exchange is stored
+ * in state and still needs to be redistributed over the nodes.
+ */
+
+void print_replica_exchange_statistics(FILE *fplog, gmx_repl_ex_t re);
+/* Should only be called on the master nodes */
+
+#endif  /* _repl_ex_h */
-- 
GitLab