Commit 2766bed5 authored by Jony Castagna's avatar Jony Castagna

Merge branch 'Abrupt_AdResS_forcecap' into 'master'

Abrupt AdResS Forcecap

See merge request e-cam/E-CAM-Library!113
parents 0133a8d0 f3874f34
......@@ -101,8 +101,8 @@ The following modules connected to the ParaDiS code have been produced so far:
:glob:
:maxdepth: 1
./modules/paradis_precipitate/paradis_precipitate_GC/readme.rst
./modules/paradis_precipitate/paradis_precipitate_HPC/readme.rst
./modules/paradis_precipitate/paradis_precipitate_GC/readme
./modules/paradis_precipitate/paradis_precipitate_HPC/readme
GC-AdResS
......@@ -114,9 +114,10 @@ Adaptive Resolution Simulation: Implementation in GROMACS
:glob:
:maxdepth: 1
./modules/GC-AdResS/Abrupt_AdResS/readme.rst
./modules/GC-AdResS/Abrupt_AdResS/abrupt_adress.rst
./modules/GC-AdResS/AdResS_RDF/readme.rst
./modules/GC-AdResS/Abrupt_AdResS/readme
./modules/GC-AdResS/Abrupt_AdResS/abrupt_adress
./modules/GC-AdResS/AdResS_RDF/readme
./modules/GC-AdResS/Abrupt_Adress_forcecap/readme
.. _ALL_background:
......
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,
.. In ReStructured Text (ReST) indentation and spacing are very important (it is how ReST knows what to do with your
document). For ReST to understand what you intend and to render it correctly please to keep the structure of this
template. Make sure that any time you use ReST syntax (such as for ".. sidebar::" below), it needs to be preceded
and followed by white space (if you see warnings when this file is built they this is a common origin for problems).
.. Firstly, let's add technical info as a sidebar and allow text below to wrap around it. This list is a work in
progress, please help us improve it. We use *definition lists* of ReST_ to make this readable.
.. sidebar:: Software Technical Information
Name
GC-AdResS -Abrupt scheme- Force Capping
Language
Implemented in GROMACS version 5.1.5
Licence
See GROMACS web page: `<http://www.gromacs.org/>`_
Documentation Tool
Application Documentation
See GROMACS web page: `<http://www.gromacs.org/>`_
Relevant Training Material
See GROMACS web page: `<http://www.gromacs.org/>`_
.. In the next line you have the name of how this module will be referenced in the main documentation (which you can
reference, in this case, as ":ref:`example`"). You *MUST* change the reference below from "example" to something
unique otherwise you will cause cross-referencing errors. The reference must come right before the heading for the
reference to work (so don't insert a comment between).
_example:
#######################
Abrupt-AdResS: Forcecap
#######################
.. Let's add a local table of contents to help people navigate the page
.. contents:: :local:
.. Add an abstract for a *general* audience here. Write a few lines that explains the "helicopter view" of why you are
creating this module. For example, you might say that "This module is a stepping stone to incorporating XXXX effects
into YYYY process, which in turn should allow ZZZZ to be simulated. If successful, this could make it possible to
produce compound AAAA while avoiding expensive process BBBB and CCCC."
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
_________________
.. Keep the helper text below around in your module by just adding ".. " in front of it, which turns it into a comment
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
______________________
.. Keep the helper text below around in your module by just adding ".. " in front of it, which turns it into a comment
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
____________________
.. Keep the helper ttxt below around in your module by just adding ".. " in front of it, which turns it into a comment
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
___________
.. Notice the syntax of a URL reference below `Text <URL>`_
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:
.. literalinclude:: ./forcecap.patch
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
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment