From 423d49bf03b6bd7f711444e48328f1c862c7d1fe Mon Sep 17 00:00:00 2001
From: Giovanni Bussi <giovanni.bussi@gmail.com>
Date: Mon, 23 Mar 2015 16:06:08 +0100
Subject: [PATCH] Fixes #132 for gmx 4.6.7
Workaround for #132.
Notice that gromacs wants something like -0.5*f*x,
whereas plumed provides something like -f*x.
Additionally, for not so clear reasons, when using plumed without computing
ENERGY the virial contribution from plumed was completely ignored.
There are two workarounds in this patch:
1. In case plumed does not ask gromacs to compute ENERGY, the plumed contribution
to virial is stored on a separate tensor and added to the gromacs one at the end
2. In case plumed asks gromacs to compute ENERGY, the a temporary tensor
is created containing twice the gromacs virial and passed to plumed. The virial
is then set to half of it.
This should solve both the issues discussed in #132
Needs gromacs to be repatched.
---
patches/gromacs-4.6.7.diff/src/kernel/md.c | 8 +++++++-
patches/gromacs-4.6.7.diff/src/mdlib/minimize.c | 8 +++++++-
2 files changed, 14 insertions(+), 2 deletions(-)
diff --git a/patches/gromacs-4.6.7.diff/src/kernel/md.c b/patches/gromacs-4.6.7.diff/src/kernel/md.c
index 741986ee9..8efc565a1 100644
--- a/patches/gromacs-4.6.7.diff/src/kernel/md.c
+++ b/patches/gromacs-4.6.7.diff/src/kernel/md.c
@@ -245,6 +245,7 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
/* PLUMED */
int plumedNeedsEnergy=0;
int plumedWantsToStop=0;
+ real plumed_vir[3][3];
/* END PLUMED */
#ifdef GMX_FAHCORE
@@ -1266,8 +1267,9 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
plumed_cmd(plumedmain,"prepareCalc",NULL);
plumed_cmd(plumedmain,"setStopFlag",&plumedWantsToStop);
plumed_cmd(plumedmain,"setForces",&f[mdatoms->start][0]);
- plumed_cmd(plumedmain,"setVirial",&force_vir[0][0]);
plumed_cmd(plumedmain,"isEnergyNeeded",&plumedNeedsEnergy);
+ for(unsigned i=0;i<3;i++) for(unsigned j=0;j<3;j++) plumed_vir[i][j]=0.0;
+ plumed_cmd(plumedmain,"setVirial",&plumed_vir[0][0]);
}
/* END PLUMED */
do_force(fplog, cr, ir, step, nrnb, wcycle, top, top_global, groups,
@@ -1279,8 +1281,12 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
/* PLUMED */
if(plumedswitch){
if(plumedNeedsEnergy){
+ for(unsigned i=0;i<3;i++) for(unsigned j=0;j<3;j++) plumed_vir[i][j]=force_vir[i][j]*2.0;
plumed_cmd(plumedmain,"setEnergy",&enerd->term[F_EPOT]);
plumed_cmd(plumedmain,"performCalc",NULL);
+ for(unsigned i=0;i<3;i++) for(unsigned j=0;j<3;j++) force_vir[i][j]=plumed_vir[i][j]*0.5;
+ } else {
+ for(unsigned i=0;i<3;i++) for(unsigned j=0;j<3;j++) force_vir[i][j]+=plumed_vir[i][j]*0.5;
}
if ((repl_ex_nst > 0) && (step > 0) && !bLastStep &&
do_per_step(step,repl_ex_nst)) plumed_cmd(plumedmain,"GREX savePositions",NULL);
diff --git a/patches/gromacs-4.6.7.diff/src/mdlib/minimize.c b/patches/gromacs-4.6.7.diff/src/mdlib/minimize.c
index f800e6225..89e705809 100644
--- a/patches/gromacs-4.6.7.diff/src/mdlib/minimize.c
+++ b/patches/gromacs-4.6.7.diff/src/mdlib/minimize.c
@@ -800,6 +800,7 @@ static void evaluate_energy(FILE *fplog, gmx_bool bVerbose, t_commrec *cr,
*/
/* PLUMED */
int plumedNeedsEnergy=0;
+ real plumed_vir[3][3];
if(plumedswitch){
long int lstep=count; (*plumedcmd)(plumedmain,"setStepLong",&count);
(*plumedcmd) (plumedmain,"setPositions",&ems->s.x[mdatoms->start][0]);
@@ -808,8 +809,9 @@ static void evaluate_energy(FILE *fplog, gmx_bool bVerbose, t_commrec *cr,
(*plumedcmd) (plumedmain,"setBox",&ems->s.box[0][0]);
(*plumedcmd) (plumedmain,"prepareCalc",NULL);
(*plumedcmd) (plumedmain,"setForces",&ems->f[mdatoms->start][0]);
- (*plumedcmd) (plumedmain,"setVirial",&force_vir[0][0]);
(*plumedcmd) (plumedmain,"isEnergyNeeded",&plumedNeedsEnergy);
+ for(unsigned i=0;i<3;i++) for(unsigned j=0;j<3;j++) plumed_vir[i][j]=0.0;
+ (*plumedcmd) (plumedmain,"setVirial",&plumed_vir[0][0]);
}
/* END PLUMED */
do_force(fplog, cr, inputrec,
@@ -823,8 +825,12 @@ static void evaluate_energy(FILE *fplog, gmx_bool bVerbose, t_commrec *cr,
/* PLUMED */
if(plumedswitch){
if(plumedNeedsEnergy) {
+ for(unsigned i=0;i<3;i++) for(unsigned j=0;j<3;j++) plumed_vir[i][j]=force_vir[i][j]*2.0;
(*plumedcmd) (plumedmain,"setEnergy",&enerd->term[F_EPOT]);
(*plumedcmd) (plumedmain,"performCalc",NULL);
+ for(unsigned i=0;i<3;i++) for(unsigned j=0;j<3;j++) force_vir[i][j]=plumed_vir[i][j]*0.5;
+ } else {
+ for(unsigned i=0;i<3;i++) for(unsigned j=0;j<3;j++) force_vir[i][j]+=plumed_vir[i][j]*0.5;
}
}
/* END PLUMED */
--
GitLab