Abrupt GC-AdResS: A new and more general implementation¶
The original idea of our proposal: to work on a general implementation of AdResS in class. MD packages. The current implementation of GC- AdResS in GROMACS has several performance problems. This module presents a very straight forward way to implement a new partitioning scheme, which solves two problems which affect the performance, the neighborlist search and the generic force kernel. Furthermore, we update the implementation to address this in a way that decouples the method directly from the core of any MD code, which does not hinder the performance and makes the scheme hardware independent. Theory, application and tests see https://aip.scitation.org/doi/10.1063/1.5031206 or https://arxiv.org/abs/1806.09870.
Purpose of Module¶
The main performance loss of AdResS simulations in GROMACS is in the neighboring list search and the generic serial force kernel, linking the atomistic (AT) and coarse grained (CG) forces together via a smooth weighting function. Thus, to get rid of the bottleneck with respect to performance and a hindrance regarding the easy/general implementation into other codes and thus get rid of the not optimized force kernel used in GROMACS we had to change the neighborlist search. This lead to a considerable speed up of the code. Furthermore it decouples the method directly from the core of any MD code, which does not hinder the performance and makes the scheme hardware independent. For the theory, application and tests see https://aip.scitation.org/doi/10.1063/1.5031206 or https://arxiv.org/abs/1806.09870.
Background Information¶
This module presents a very straight forward way to implement a new partitioning scheme. And this solves two problems which affect the performance, the neighborlist search and the generic force kernel.
In GROMACS the neighbor list is put together and organized in the file ‘ns.c’. In GROMACS 5.1
there are two functions which basically sort the incoming
particles into the different neighbor list. In its current official GROMACS release everything other than CG (with 
) or AT (with 
) is sorted into the neighbor lists. Any other particles are sorted into a special neighbor list only for AdResS.
We now changed this neighborlist sorting into: Everything is taken into account other than: (AT and ( 
or 
)) or (CG and ( 
and 
)). This leads to 5 distinct interactions: (1) AT-AT in the atomistic region, (2) CG-CG in the
CG region, (3) AT-AT between particles in the hybrid region, (4) AT-AT between particles of the
atomistic region with the hybrid region and (5) CG-CG between particles of the CG region with the
hybrid region. This if statement excludes the CG-CG interaction in the hybrid region.
This is a very straight forward way to implement a new partitioning scheme and utilize a constant weighting function. This solves both parts of the performance problem, the neighborlist search and the generic force kernel which can be simply switch off by switching to the standard interaction scheme implemented in GROMACS.
Building and Testing¶
We tested this new implementation on SPC water with varying system sizes. GROMACS is optimized especially for handling of bio-systems, i.e. GROMACS has the best performance in case of water simulations. We set up a couple of Abrupt GC-AdResS simulations ranging from small 6912 water molecules to 48k water molecules. We used a standard desktop machine (Intel Core i5-4590 CPU @ 3.30GHz x4) and run small 20 ps runs. We can see that our performance is much improved up to a factor of 2.5.
Here is a short manual on how to run the test and set up AdResS simulations in GROMACS:
mandatory requirement: a working full atomistic simulation (molecular configuration, force fields and optimal MD input parameter)
You need a very well converged NVT run, which can be used as starting point for the coarse grained (CG) and then later the AdResS simulations.
You have to generate a coarse grained (CG) potential. We use, for convenience, the inverse Boltzmann iteration provided in the VOTCA package (http://www.votca.org/home). The resulting tabulated CG potentials are used in the AdResS simulation. Alternatively you can use WCA potentials or standard Lennard-Jones potentials. The main requirement for the AdResS simulation is that the density in the CG region is the same as in the atomistic (AT) region.
NOTE: The method can be used with any potential, which preserves the correct density. If only a SPC/E CG potential is available it can be used for SPC/e water models as well as for a more advanced water model. It is possible to use a WCA potential, which is basically a Lennard-Jones potential. And it is possible to switch the CG ptential completely off. That will transform the CG region to a true thermodynamic reservoir with a non-interacting gas.
The next step is to create a double resolution configuration and adjust the dependencies (force field, topology, index file, GROMACS input file). Creating the configuration is straight forward (we use http://www.votca.org/home).
Example from VOTCA:
csg_map --top topol.tpr --cg cg_mapping_scheme --hybrid --trj conf.gro --out conf_hybrid.gro
Of course, if you want to use this configuration in a MD simulation you have to adjust the force field (see example file: spc.adress.itp). You have to define a virtual side:
[ virtual_sites3 ]
; Site from funct a d
; atom dependencies func a b
4 1 2 3 1 0.05595E+00 0.05595E+00
The next step is to adjust the status of the CG particle in the topology file (in our example: topol.top) from A for atom to V as virtual particle. And of course insert the new force field.
#include "spc.adress.itp"
Then you have to generate an index file with the different energy groups. In this example, we have 2 groups (EXW and WCG, the name of the CG particle):
gmx make_ndx -f conf_hybrid.gro
> a WCG
Found 3456 atoms with name WCG
3 WCG : 3456 atoms
> !3
Copied index group 3 'WCG'
Complemented group: 10368 atoms
4 !WCG : 10368 atoms
> name 4 EXW
> q
The next step is to adjust the GROMACS input file. AdResS needs the Langevin dynamics, so you have to choose:
integrator = sd
Since the system is double resolution, meaning we have the atomistci detals and the virtual particles, we have to define the energygroups:
; Selection of energy groups
energygrps = EXW WCG
energygrp_table = WCG WCG
GROMACS version 5.1.5 is using verlet as standard cutoff-scheme, so we have to change that to group:
; nblist update frequency
cutoff-scheme = group
Furthermore, in our simulations we use:
coulombtype = reaction-field
rcoulomb = 1.0
vdw-type = user
rvdw = 1.0
In case of local thermostat simulations (see https://aip.scitation.org/doi/10.1063/1.5031206 or https://arxiv.org/abs/1806.09870) we use:
coulombtype = reaction-field-zero
rcoulomb = 1.0
vdw-type = user
rvdw = 1.0
If you use the stochastic dynamics, we add the following entries to make sure we have only NVT and a thermalization via the Langevin dynamics.
; Temperature coupling
Tcoupl = no
Pcoupl = no
To switch the simulation to AdResS this is the key part. This starts the AdResS runs.
; AdResS parameters adress = yes ;no
Here you define the geometry of the atomistic region, either sphere (a spherical region anywhere in the simualtion box) or xsplit (a cuboid slice of the whole simulation box for the atomistic region, with the transition and coarse grained region on each side).
adress_type = sphere ;xsplit sphere or constant
This defines the width of the atomistic region, starting from the given reference coordinate (keyword adress_reference_coords, by simply using: tail conf_hybrid.gro | awk ‘(NF==3){print $1/2., $2/2., $3/2.}’). In the older versions of AdResS, with a smooth coupling between AT and CG the width of the hybrid region width (adress_hy_width) was also defined. In the Abrupt_AdResS setup it is not necessary any more, even if you put a number that region is counted (in the code) as AT.
adress_ex_width = 1.5
adress_hy_width = 1.5
adress_ex_forcecap = 2000
adress_interface_correction = thermoforce ;off
adress_site = com
adress_reference_coords = 3.7500 1.860355 1.860355
adress_tf_grp_names = WCG
adress_cg_grp_names = WCG
adress_do_hybridpairs = no
Another important aspect is the force capping. Abrupt AdResS works fine for small molecules like water, but for larger or more complex molecules the force capping is very important. We cap every force component (i.e. f(x),f(y),f(z)) acting on a particle and not the norm of the force, which reduces the computational time spend. This is described in another module.
adress_interface_correction defines if you use an external force to correct the density or not. In case of the old AdResS (smooth coupling) that correction simply refined the simulation, as the density difference was not significant. For the Abrupt AdResS, and the method development based on it, and more complex molecules (i.e. polymers) the thermodynamic force is essential. If it is not taken into account the risk to form interfaces between AT and CG is high. Also if particles coming too close (basically overlap) the run can crash. The role of the thermodynamic force, the force cap and the basic theory behind it see https://aip.scitation.org/doi/10.1063/1.5031206 or https://arxiv.org/abs/1806.09870. For this to work you must have a file e.g. in our example case: tabletf_WCG.xvg in the directory, otherwise you have to set:
adress_interface_correction = off
There is a number of properties you have to check. The first check is always the density and you see if the patch works from the density. If you have no thermodynamic force you have rather pronounced spikes in the density at the interfaces. If you have a converged thermodynamic force the density has to be within +/- 3% off from a comparable full atomistc simulation / experimental data. Then you need further properties to make sure you have an open system. The problem with the simulation is that an “artificial” interface is introduced and checks for the diffusion, the RDF’s… (full list see below) ensure that those regions mix and that you have proper particle transfer.
Source Code¶
The patch file for Abrupt GC-Adress is:
diff -ru /storage/mi/ck69giso/gromacs-5.1.5/src/gromacs/mdlib/adress.c /home/mi/ck69giso/gmx-515-hck/src/gromacs/mdlib/adress.c
--- /storage/mi/ck69giso/gromacs-5.1.5/src/gromacs/mdlib/adress.c 2016-07-13 14:56:04.000000000 +0200
+++ /home/mi/ck69giso/gmx-515-hck/src/gromacs/mdlib/adress.c 2018-08-15 12:39:32.000000000 +0200
@@ -101,17 +101,17 @@
return 0;
}
/* molecule is explicit */
- else if (sqr_dl < adressr*adressr)
+ else //if (sqr_dl < adressr*adressr)
{
return 1;
}
- /* hybrid region */
+ /* hybrid region
else
- {
- dl = sqrt(sqr_dl);
- tmp = cos((dl-adressr)*M_PI/2/adressw);
- return tmp*tmp;
- }
+ {
+// dl = sqrt(sqr_dl);
+// tmp = cos((dl-adressr)*M_PI/2/adressw);
+ return 0.5;
+ } */
}
void
diff -ru /storage/mi/ck69giso/gromacs-5.1.5/src/gromacs/mdlib/ns.c /home/mi/ck69giso/gmx-515-hck/src/gromacs/mdlib/ns.c
--- /storage/mi/ck69giso/gromacs-5.1.5/src/gromacs/mdlib/ns.c 2016-07-13 14:56:04.000000000 +0200
+++ /home/mi/ck69giso/gmx-515-hck/src/gromacs/mdlib/ns.c 2018-08-15 12:55:06.000000000 +0200
@@ -264,10 +264,12 @@
for (i = 0; i < fr->nnblists; i++)
{
nbl = &(fr->nblists[i]);
-
- if ((fr->adress_type != eAdressOff) && (i >= fr->nnblists/2))
+/* chk */
+// if ((fr->adress_type != eAdressOff) && (i >= fr->nnblists/2))
+ if ((fr->adress_type != eAdressOff))
{
- type = GMX_NBLIST_INTERACTION_ADRESS;
+ /* type = GMX_NBLIST_INTERACTION_ADRESS; */
+ type = GMX_NBLIST_INTERACTION_STANDARD;
}
init_nblist(log, &nbl->nlist_sr[eNL_VDWQQ], &nbl->nlist_lr[eNL_VDWQQ],
maxsr, maxlr, ivdw, ivdwmod, ielec, ielecmod, igeometry_def, type, bElecAndVdwSwitchDiffers);
@@ -601,11 +603,14 @@
int *cginfo;
int *type, *typeB;
real *charge, *chargeB;
+ real *wf;
real qi, qiB, qq, rlj;
gmx_bool bFreeEnergy, bFree, bFreeJ, bNotEx, *bPert;
gmx_bool bDoVdW_i, bDoCoul_i, bDoCoul_i_sol;
+ gmx_bool b_hybrid;
int iwater, jwater;
t_nblist *nlist;
+ gmx_bool bEnergyGroupCG;
/* Copy some pointers */
cginfo = fr->cginfo;
@@ -614,6 +619,7 @@
type = md->typeA;
typeB = md->typeB;
bPert = md->bPerturbed;
+ wf = md->wf;
/* Get atom range */
i0 = index[icg];
@@ -625,7 +631,7 @@
iwater = (solvent_opt != esolNO) ? GET_CGINFO_SOLOPT(cginfo[icg]) : esolNO;
bFreeEnergy = FALSE;
- if (md->nPerturbed)
+ if (md->nPerturbed )
{
/* Check if any of the particles involved are perturbed.
* If not we can do the cheaper normal put_in_list
@@ -683,6 +689,7 @@
}
if (!bFreeEnergy)
+/* if (!bFreeEnergy || (fr->adress_type != eAdressOff)) */
{
if (iwater != esolNO)
{
@@ -846,8 +853,13 @@
bDoVdW_i = (bDoVdW && bHaveVdW[type[i_atom]]);
bDoCoul_i = (bDoCoul && qi != 0);
+ /* chk */
+ bEnergyGroupCG = !egp_explicit(fr, igid);
+ /* chk */
+
if (bDoVdW_i || bDoCoul_i)
{
+
/* Loop over the j charge groups */
for (j = 0; (j < nj); j++)
{
@@ -867,7 +879,19 @@
/* Finally loop over the atoms in the j-charge group */
for (jj = jj0; jj < jj1; jj++)
{
+
bNotEx = NOTEXCL(bExcl, i, jj);
+ /* change 7.11.2017 chk*/
+ if ( fr->adress_type != eAdressOff )
+ {
+ if ( ( !bEnergyGroupCG && ( wf[i_atom] <= GMX_REAL_EPS || wf[jj] <= GMX_REAL_EPS ) ) ||
+ ( ( bEnergyGroupCG ) && ( wf[i_atom] > GMX_REAL_EPS || wf[jj] > GMX_REAL_EPS ) ))
+// abrupt-GC ( ( bEnergyGroupCG ) && ( wf[i_atom] > GMX_REAL_EPS && wf[jj] > GMX_REAL_EPS ) ))
+ {
+ continue;
+ }
+ }
+ /* change 7.11.2017 chk*/
if (bNotEx)
{
@@ -984,6 +1008,10 @@
jj1 = index[jcg+1];
/* Finally loop over the atoms in the j-charge group */
bFree = bPert[i_atom];
+
+ /* chk
+ bEnergyGroupCG = !egp_explicit(fr, igid); */
+
for (jj = jj0; (jj < jj1); jj++)
{
bFreeJ = bFree || bPert[jj];
@@ -994,6 +1022,16 @@
{
bNotEx = NOTEXCL(bExcl, i, jj);
+ /* chk
+
+ if ( ( !bEnergyGroupCG && ( wf[i_atom] <= GMX_REAL_EPS || wf[jj] <= GMX_REAL_EPS ) ) ||
+ ( ( bEnergyGroupCG ) && ( wf[i_atom] >= GMX_REAL_EPS && wf[jj] >= GMX_REAL_EPS ) )
+ )
+
+ {
+ continue;
+ } */
+
if (bNotEx)
{
if (bFreeJ)
@@ -1250,6 +1288,19 @@
* b_hybrid=true are placed into the _adress neighbour lists and
* processed by the generic AdResS kernel.
*/
+ /* change 7.11.2017 chk*/
+ /* if ( fr->adress_type != eAdressOff )
+ { */
+ if ( ( !bEnergyGroupCG && ( wf[i_atom] <= GMX_REAL_EPS || wf[jj] <= GMX_REAL_EPS ) ) ||
+ ( ( bEnergyGroupCG ) && ( wf[i_atom] > GMX_REAL_EPS || wf[jj] > GMX_REAL_EPS ) )
+// ( ( bEnergyGroupCG ) && ( wf[i_atom] > GMX_REAL_EPS && wf[jj] > GMX_REAL_EPS ) )
+ )
+
+ {
+ continue;
+ }
+ /* } */
+/* old version from normal GC-AdResS before october 7
if ( (bEnergyGroupCG &&
wf[i_atom] >= 1-GMX_REAL_EPS && wf[jj] >= 1-GMX_REAL_EPS ) ||
( !bEnergyGroupCG && wf[jj] <= GMX_REAL_EPS ) )
@@ -1259,6 +1310,7 @@
b_hybrid = !((wf[i_atom] >= 1-GMX_REAL_EPS && wf[jj] >= 1-GMX_REAL_EPS) ||
(wf[i_atom] <= GMX_REAL_EPS && wf[jj] <= GMX_REAL_EPS));
+*/
if (bNotEx)
{
@@ -1266,28 +1318,15 @@
{
if (charge[jj] != 0)
{
- if (!b_hybrid)
- {
- add_j_to_nblist(coul, jj, bLR);
- }
- else
- {
- add_j_to_nblist(coul_adress, jj, bLR);
- }
+ /* chk: removed the !b_hybrid if loops */
+ add_j_to_nblist(coul, jj, bLR);
}
}
else if (!bDoCoul_i)
{
if (bHaveVdW[type[jj]])
{
- if (!b_hybrid)
- {
add_j_to_nblist(vdw, jj, bLR);
- }
- else
- {
- add_j_to_nblist(vdw_adress, jj, bLR);
- }
}
}
else
@@ -1296,38 +1335,16 @@
{
if (charge[jj] != 0)
{
- if (!b_hybrid)
- {
add_j_to_nblist(vdwc, jj, bLR);
- }
- else
- {
- add_j_to_nblist(vdwc_adress, jj, bLR);
- }
}
else
{
- if (!b_hybrid)
- {
add_j_to_nblist(vdw, jj, bLR);
- }
- else
- {
- add_j_to_nblist(vdw_adress, jj, bLR);
- }
-
}
}
else if (charge[jj] != 0)
{
- if (!b_hybrid)
- {
add_j_to_nblist(coul, jj, bLR);
- }
- else
- {
- add_j_to_nblist(coul_adress, jj, bLR);
- }
}
}
@@ -2671,8 +2688,10 @@
rvec box_size, grid_x0, grid_x1;
int i, j, m, ngid;
real min_size, grid_dens;
+ real b_hybrid;
int nsearch;
gmx_bool bGrid;
+ gmx_bool bEnergyGroupCG;
char *ptr;
gmx_bool *i_egp_flags;
int cg_start, cg_end, start, end;
@@ -2774,8 +2793,9 @@
}
debug_gmx();
- if (fr->adress_type == eAdressOff)
- {
+/* chk */
+// if (fr->adress_type == eAdressOff)
+// {
if (!fr->ns.bCGlist)
{
put_in_list = put_in_list_at;
@@ -2784,11 +2804,12 @@
{
put_in_list = put_in_list_cg;
}
- }
- else
- {
- put_in_list = put_in_list_adress;
- }
+// }
+// else
+// {
+// put_in_list = put_in_list_adress;
+// }
+/* chk */
/* Do the core! */
if (bGrid)
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,
To apply the patch:
- copy into the main directory (gromacs/)
- patch < abrupt_adress.patch
In this module we also include a test scenario for GROMACS version 5.1.5 with a possible CG potential and all necessary input files. To run it simply run gmx grompp -f grompp.mdp -c conf.gro -p topol.top -n index.ndx -maxwarn 5; gmx mdrun using the patched version of GROMACS version 5.1.5 (see above).
When gmx mdrun finished normally (with the above mentioned setup), we have several mandatory checks to see if the simulation was successful or not.
- Easiest check: load the conf.gro and the trajectory file in vmd and check if you see particle diffusion or depleted areas.
- we check the density along the X-direction (xsplit: e.g. gmx density -f traj_compt.xtc -d X) or along the radius (sphere: e.g. via VOTCA: csg_density –axis r –rmax <value> –ref [x_ref,y_ref,z_ref] –trj traj_comp.xtc –top topol.tpr –out test.dens.comp), the density has to be less then 3% different from experimental data or the density from a full atomistic MD simulation. The density of the example is 1000 kg m^-3.
- static properties: crucial RDF’s (e.g. for water the oxygen-oxygen RDF)
- p(N): It describes the average number of particles in the AT region throughout the simulation.
- the density diffusion for each region (via a very helpful expansion for http://www.ks.uiuc.edu/Research/vmd/, the density profile tool see https://github.com/tonigi/vmd_density_profile).
- If we only thermalize the transition region, the AT region is NVE-like, which means it is even possible to determine the dynamics of the system.
The files for the water example can be found here:
spc-example.tar.gz