Commit a84d9431 authored by Jony Castagna's avatar Jony Castagna
Browse files

Revert "Merge branch 'Abrupt_AdResS' into 'master'"

This reverts merge request !111
parent 4f13fddb
......@@ -113,16 +113,3 @@ The first Meso- and Multi-scale ESDW was held in Barcelona, Spain, in July 2017.
:maxdepth: 1
./modules/DL_MESO_DPD/sionlib_dlmeso_dpd/readme
GC-AdResS
---------
Adaptive Resolution Simulation: Implementation in GROMACS
.. toctree::
:glob:
:maxdepth: 1
./modules/GC-AdResS/Abrupt_AdResS/readme.rst
./modules/GC-AdResS/Abrupt_AdResS/abrupt_adress.rst
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,
.. _abrupt_adress_patch:
####################################
Patch file for module: Abrupt AdResS
####################################
The patch for the abrupt AdResS code is:
.. literalinclude:: ./abrupt_adress.patch
:linenos:
:download:`Downloadable version of patch file <abrupt_adress.patch>`
<
.. 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
Language
Implemented in GROMACS version 5.1.5
Licence
MD Simulation:
See GROMACS web page: `<http://www.gromacs.org>`_
Analysis tools and thermodynamic force calculation:
see VOTCA web page: `<http://www.votca.org/home>`_
Documentation Tool
Application Documentation
See GROMACS web page: `<http://www.gromacs.org>`_
See VOCTA web page: `<http://www.votca.org/Documentation>`_
Relevant Training Material
See GROMACS web page: `<http://www.gromacs.org>`_
See VOCTA web page: `<http://www.votca.org/tutorials>`_
.. 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).
.. _abrupt_adress:
#######################################################
Abrupt GC-AdResS: A new and more general implementation
#######################################################
.. Let's add a local table of contents to help people navigate the page
.. contents:: :local:
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>`_.
.. 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."
.. This is an example of what a *module* for E-CAM looks like. The original source of this page (:download:`readme.rst`) contains lots of additional comments to help you create your module (and understand ReST_ syntax) so please use this as a starting point. You are free add any level of complexity you wish (within the bounds of what ReST_ can do). More general instructions for making your contribution can be found in ":ref:`contributing`".
.. Remember that for a module to be accepted into the E-CAM repository, your source code changes in the target application must pass a number of acceptance criteria:
.. * Style *(use meaningful variable names, no global variables,...)*
.. * Source code documentation *(each function should be documented with each argument explained)*
.. * Tests *(everything you add should have either unit or regression tests)*
.. * Performance *(If what you introduce has a significant computational load you should make some performance optimisation effort using an appropriate tool. You should be able to verify that your changes have not introduced unexpected performance penalties, are threadsafe if needed,...)*
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
.. Give a brief overview of why the module is/was being created, explaining a little of the scientific background and how it fits into the larger picture of what you want to achieve.
.. If needed you can include latex mathematics like
.. :math:`\frac{ \sum_{t=0}^{N}f(t,k) }{N}`
.. which won't show up on GitLab/GitHub but will in final online documentation.
.. If you want to add a citation, such as [CIT2009]_. Note that citations may get rearranged, e.g., to the bottom of the "page".
.. : .. [CIT2009] A citation (as often used in journals).
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. We know that 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>`_.
.. The interface between the regions is more fluctuating and needs a more responsive thermodynamic force but it works reasonably well.
.. The second piece of the puzzle is the spatial partitioning as we showed at the ESDW8 in Berlin and as Guzman et al. (arXiv:1711.03290v1) published recently it is possible to use a spatial partitioning for GC-AdResS.
Background Information
______________________
.. Keep the helper text below around in your module by just adding ".. " in front of it, which turns it into a comment
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 :math:`w_i=w_j=1`) or AT (with :math:`w_i=w_j=0` ) 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 ( :math:`w_i=0` or :math:`w_j=0`)) or (CG and ( :math:`w_i>=0` and :math:`w_j>=0`)). 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
____________________