Abrupt-AdResS: Forcecap¶
The original idea of our proposal was to work on a general implementation of grand canonical adaptive resolution simulations (GC-AdResS) in classical MD packages. The current implementation of GC- AdResS in GROMACS (up to version 5.1.5) has performance problems. The Abrupt GC-AdResS implementation is avoiding those and make AdResS more interesting for other MD developers (especially since we could remove the force interpolation and weighting functions from the force kernel).
This module works in combination with abrupt_AdResS and is at the same time important for a successful simulation. It shows how to avoid particle overlap at the interface of the atomistic and coarse grained regions.
Purpose of Module¶
The implementation of Abrupt GC-AdResS is in itself only working for the smallest and simplest of molecules without problems. For larger and more complex molecules the simulation crashes. This module shows a way to avoid this.
Background Information¶
As studies of ionic liquids and polymer melts have shown for large and complicated molecules even the standard GC-AdResS is not working without additional force capping. The reason for that is when a molecule from the coarse grained region enters the hybrid region the atomistic representations, which are present due to the technically necessary double resolution, interact. It is possible for atoms to be too close together, which results in a too high force and thus in too high velocities of those particles. Since in Abrupt AdResS we can avoid the generic force kernel from GROMACS, the force capping (which was previously implemented at the end of the force calculation) had to be shifted and replaced. We finally looked at the integrator (in our case the stochastic dynamics integrator), which is the place where each force has to be read and handled.
This implementation of force capping is a rudimentary approach. The basic principle is when two particles are too close together, and thus the force are far higher than the average forces in the simulation, the force on the particles are re-scaled to a given value. That tactic makes sure that the insertion of particles in the atomistic region is introduced at reasonable velocities and temperature. As a side effect, the area of disturbance due to the introduction of a particle is limited.
Building and Testing¶
We have used this new addition to the code on two systems: water and ionic liquids. The results have been published in Ref. https://aip.scitation.org/doi/10.1063/1.5031206 or https://arxiv.org/abs/1806.09870. All the information about studied systems and the performance can be found there.
The patch provided can be applied alone without the Abrupt AdResS patch, in the main directory of GROMACS. The important part of the patch is that it adds an upper force limit in the stochastic dynamics integrator. If that upper limit is triggered the force is re-scaled to the given force cap. The rest is basically to make sure that the integrator is called correctly. There is a print command which is triggered once the force on a particle is higher than a given force cap value. The force capping simulation can in some cases cause lincs warnings. Since we take care of faulty configuration that way, we can disabling those warnings (export GMX_MAXBACKUP=-1 ; export GMX_MAXCONSTRWARN=-1). Otherwise it generates too many files and can crash as well.
Source Code¶
A note of caution: the chosen force cap trigger has to be a high enough value, otherwise normal interactions (interactions with forces around the average forces in a simulation) will trigger the force capping. That would change the dynamics and the structure of a system. Also it would decrease the performance of the code. If chosen too high it might run into impossible and unstable configuration, which will result in a program crash.
Recipe for hard coded force capping:
diff -ru /storage/mi/ck69giso/gromacs-5.1.5/src/gromacs/mdlib/update.cpp /home/mi/ck69giso/gmx-515-hck/src/gromacs/mdlib/update.cpp
--- /storage/mi/ck69giso/gromacs-5.1.5/src/gromacs/mdlib/update.cpp 2016-09-07 14:50:21.000000000 +0200
+++ /home/mi/ck69giso/gmx-515-hck/src/gromacs/mdlib/update.cpp 2018-07-24 16:07:27.000000000 +0200
@@ -67,6 +67,7 @@
#include "gromacs/utility/futil.h"
#include "gromacs/utility/gmxomp.h"
#include "gromacs/utility/smalloc.h"
+#include "adress.h"
/*For debugging, start at v(-dt/2) for velolcity verlet -- uncomment next line */
/*#define STARTFROMDT2*/
@@ -569,6 +570,8 @@
return upd;
}
+/* new */
+
static void do_update_sd1(gmx_stochd_t *sd,
int start, int nrend, double dt,
rvec accel[], ivec nFreeze[],
@@ -579,10 +582,11 @@
int ngtc, real ref_t[],
gmx_bool bDoConstr,
gmx_bool bFirstHalfConstr,
- gmx_int64_t step, int seed, int* gatindex)
+ gmx_int64_t step, int seed, int* gatindex, real fc)
{
gmx_sd_const_t *sdc;
gmx_sd_sigma_t *sig;
+
real kT;
int gf = 0, ga = 0, gt = 0;
real ism;
@@ -625,10 +629,21 @@
{
if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
{
- real sd_V, vn;
+// real sd_V, vn;
+ real sd_V, vn, fn;
+ fn = f[n][d];
+
+// fc = 10000.;
+
+ if (fabs(fn)>fc)
+ {
+ printf("SD (I) force-cap %e\n", fn);
+ fn = fc*fn/fabs(fn);
+ }
sd_V = ism*sig[gt].V*rnd[d];
- vn = v[n][d] + (invmass[n]*f[n][d] + accel[ga][d])*dt;
+ vn = v[n][d] + (invmass[n]*fn + accel[ga][d])*dt;
+// vn = v[n][d] + (invmass[n]*f[n][d] + accel[ga][d])*dt;
v[n][d] = vn*sdc[gt].em + sd_V;
/* Here we include half of the friction+noise
* update of v into the integration of x.
@@ -668,7 +683,20 @@
{
if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
{
- v[n][d] = v[n][d] + (im*f[n][d] + accel[ga][d])*dt;
+
+ real fn;
+
+// fc = 10000.;
+
+ fn = f[n][d];
+ if (fabs(fn)>fc)
+ {
+ printf("SD (II) force-cap %e\n", fn);
+ fn = fc*fn/fabs(fn);
+ }
+
+ v[n][d] = v[n][d] + (im*fn + accel[ga][d])*dt;
+// v[n][d] = v[n][d] + (im*f[n][d] + accel[ga][d])*dt;
xprime[n][d] = x[n][d] + v[n][d]*dt;
}
else
@@ -1644,6 +1672,8 @@
end_th = start + ((nrend-start)*(th+1))/nth;
/* The second part of the SD integration */
+ if (inputrec->bAdress)
+ {
do_update_sd1(upd->sd,
start_th, end_th, dt,
inputrec->opts.acc, inputrec->opts.nFreeze,
@@ -1653,7 +1683,23 @@
inputrec->opts.ngtc, inputrec->opts.ref_t,
bDoConstr, FALSE,
step, inputrec->ld_seed,
- DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL);
+ DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL,
+ inputrec->adress->ex_forcecap);
+ }
+ else
+ {
+ do_update_sd1(upd->sd,
+ start_th, end_th, dt,
+ inputrec->opts.acc, inputrec->opts.nFreeze,
+ md->invmass, md->ptype,
+ md->cFREEZE, md->cACC, md->cTC,
+ state->x, xprime, state->v, force,
+ inputrec->opts.ngtc, inputrec->opts.ref_t,
+ bDoConstr, FALSE,
+ step, inputrec->ld_seed,
+ DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL,
+ 5000.);
+ }
}
inc_nrnb(nrnb, eNR_UPDATE, homenr);
wallcycle_stop(wcycle, ewcUPDATE);
@@ -2031,6 +2077,21 @@
break;
case (eiSD1):
/* With constraints, the SD1 update is done in 2 parts */
+ if (inputrec->bAdress)
+ {
+ do_update_sd1(upd->sd,
+ start_th, end_th, dt,
+ inputrec->opts.acc, inputrec->opts.nFreeze,
+ md->invmass, md->ptype,
+ md->cFREEZE, md->cACC, md->cTC,
+ state->x, xprime, state->v, force,
+ inputrec->opts.ngtc, inputrec->opts.ref_t,
+ bDoConstr, TRUE,
+ step, inputrec->ld_seed, DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL,
+ inputrec->adress->ex_forcecap);
+ }
+ else
+ {
do_update_sd1(upd->sd,
start_th, end_th, dt,
inputrec->opts.acc, inputrec->opts.nFreeze,
@@ -2039,7 +2100,9 @@
state->x, xprime, state->v, force,
inputrec->opts.ngtc, inputrec->opts.ref_t,
bDoConstr, TRUE,
- step, inputrec->ld_seed, DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL);
+ step, inputrec->ld_seed, DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL,
+ 5000.);
+ }
break;
case (eiSD2):
/* The SD2 update is always done in 2 parts,
with:
fc = Chosen upper force limit for the ionic liquids simulations
fn = force acting on a particle
The patch provided can be applied in the main directory of GROMACS via:
patch < forcecap.patch