From 11be00e412e2fe2375ab0da2726660469c8173b3 Mon Sep 17 00:00:00 2001 From: Berk Hess Date: Mon, 13 May 2019 14:39:16 +0200 Subject: [PATCH] Allow AWH geometry 'direction' to be periodic When an AWH dimension with pull geometry 'direction' has an interval length (nearly) matching the periodic unit cell dimension, it is now made periodic. Fixes #2946 Change-Id: I9a09b8eb7464e04bd3d72059324bceb590e531f5 --- docs/release-notes/2019/2019.3.rst | 9 ++++ docs/user-guide/mdp-options.rst | 6 ++- src/gromacs/awh/read-params.cpp | 61 ++++++++++++++++++++++++---- src/gromacs/awh/read-params.h | 4 +- src/gromacs/gmxpreprocess/grompp.cpp | 10 ++++- 5 files changed, 77 insertions(+), 13 deletions(-) diff --git a/docs/release-notes/2019/2019.3.rst b/docs/release-notes/2019/2019.3.rst index 6864d3729b..a689072367 100644 --- a/docs/release-notes/2019/2019.3.rst +++ b/docs/release-notes/2019/2019.3.rst @@ -37,6 +37,15 @@ Fix segmentation fault when using membrane embedding :issue:`2947` +Allow AWH with pull-geometry 'direction' to be periodic +""""""""""""""""""""""""""""""""""""""""""""""""""""""" + +When applying AWH to a pull coordinate with geometry 'direction' +with a AWH interval length of more than 95% of the box size, +the dimension is now made periodic. + +:issue:`2946` + Fixes for ``gmx`` tools ^^^^^^^^^^^^^^^^^^^^^^^ diff --git a/docs/user-guide/mdp-options.rst b/docs/user-guide/mdp-options.rst index 1f13bdd54c..137a2718de 100644 --- a/docs/user-guide/mdp-options.rst +++ b/docs/user-guide/mdp-options.rst @@ -2143,8 +2143,12 @@ AWH adaptive biasing (0.0) [nm] or [rad] Start value of the sampling interval along this dimension. The range of allowed values depends on the relevant pull geometry (see :mdp:`pull-coord1-geometry`). - For periodic geometries :mdp:`awh1-dim1-start` greater than :mdp:`awh1-dim1-end` + For dihedral geometries :mdp:`awh1-dim1-start` greater than :mdp:`awh1-dim1-end` is allowed. The interval will then wrap around from +period/2 to -period/2. + For the direction geometry, the dimension is made periodic when + the direction is along a box vector and covers more than 95% + of the box length. Note that one should not apply pressure coupling + along a periodic dimension. .. mdp:: awh1-dim1-end diff --git a/src/gromacs/awh/read-params.cpp b/src/gromacs/awh/read-params.cpp index a8db35eed7..b9843fb457 100644 --- a/src/gromacs/awh/read-params.cpp +++ b/src/gromacs/awh/read-params.cpp @@ -554,20 +554,43 @@ AwhParams *readAndCheckAwhParams(std::vector *inp, const t_inputrec * * Gets the period of a pull coordinate. * * \param[in] pullCoordParams The parameters for the pull coordinate. + * \param[in] pbc The PBC setup + * \param[in] intervalLength The length of the AWH interval for this pull coordinate * \returns the period (or 0 if not periodic). */ -static double get_pull_coord_period(const t_pull_coord &pullCoordParams) +static double get_pull_coord_period(const t_pull_coord &pullCoordParams, + const t_pbc &pbc, + const real intervalLength) { - double period; + double period = 0; - if (pullCoordParams.eGeom == epullgDIHEDRAL) + if (pullCoordParams.eGeom == epullgDIR) { - /* The dihedral angle is periodic in -180 to 180 deg */ - period = 360; + const real margin = 0.001; + // Make dims periodic when the interval covers > 95% + const real periodicFraction = 0.95; + + // Check if the pull direction is along a box vector + for (int dim = 0; dim < pbc.ndim_ePBC; dim++) + { + const real boxLength = norm(pbc.box[dim]); + const real innerProduct = iprod(pullCoordParams.vec, pbc.box[dim]); + if (innerProduct >= (1 - margin)*boxLength && + innerProduct <= (1 + margin)*boxLength) + { + GMX_RELEASE_ASSERT(intervalLength < (1 + margin)*boxLength, + "We have checked before that interval <= period"); + if (intervalLength > periodicFraction*boxLength) + { + period = boxLength; + } + } + } } - else + else if (pullCoordParams.eGeom == epullgDIHEDRAL) { - period = 0; + /* The dihedral angle is periodic in -180 to 180 deg */ + period = 360; } return period; @@ -678,7 +701,7 @@ static void checkInputConsistencyInterval(const AwhParams *awhParams, warninp_t void setStateDependentAwhParams(AwhParams *awhParams, const pull_params_t *pull_params, pull_t *pull_work, - const matrix box, int ePBC, + const matrix box, int ePBC, const tensor &compressibility, const t_grpopts *inputrecGroupOptions, warninp_t wi) { /* The temperature is not really state depenendent but is not known @@ -719,7 +742,27 @@ void setStateDependentAwhParams(AwhParams *awhParams, } - dimParams->period = get_pull_coord_period(pullCoordParams); + dimParams->period = get_pull_coord_period(pullCoordParams, pbc, dimParams->end - dimParams->origin); + // We would like to check for scaling, but we don't have the full inputrec available here + if (dimParams->period > 0 && !(pullCoordParams.eGeom == epullgANGLE || + pullCoordParams.eGeom == epullgDIHEDRAL)) + { + bool coordIsScaled = false; + for (int d2 = 0; d2 < DIM; d2++) + { + if (pullCoordParams.vec[d2] != 0 && norm2(compressibility[d2]) != 0) + { + coordIsScaled = true; + } + } + if (coordIsScaled) + { + std::string mesg = gmx::formatString("AWH dimension %d of bias %d is periodic with pull geometry '%s', " + "while you should are applying pressure scaling to the corresponding box vector, this is not supported.", + d + 1, k + 1, EPULLGEOM(pullCoordParams.eGeom)); + warning(wi, mesg.c_str()); + } + } /* The initial coordinate value, converted to external user units. */ dimParams->coordValueInit = diff --git a/src/gromacs/awh/read-params.h b/src/gromacs/awh/read-params.h index 5c5be7a4f9..ef2f8eab06 100644 --- a/src/gromacs/awh/read-params.h +++ b/src/gromacs/awh/read-params.h @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by + * Copyright (c) 2013,2014,2015,2016,2017,2018,2019, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -78,6 +78,7 @@ AwhParams *readAndCheckAwhParams(std::vector *inp, * \param[in,out] pull_work Pull working struct to register AWH bias in. * \param[in] box Box vectors. * \param[in] ePBC Periodic boundary conditions enum. + * \param[in] compressibility Compressibility matrix for pressure coupling, pass all 0 without pressure coupling * \param[in] inputrecGroupOptions Parameters for atom groups. * \param[in,out] wi Struct for bookeeping warnings. * @@ -88,6 +89,7 @@ void setStateDependentAwhParams(AwhParams *awhParams, pull_t *pull_work, const matrix box, int ePBC, + const tensor &compressibility, const t_grpopts *inputrecGroupOptions, warninp_t wi); diff --git a/src/gromacs/gmxpreprocess/grompp.cpp b/src/gromacs/gmxpreprocess/grompp.cpp index 2485d5a2e2..2a87cec0c4 100644 --- a/src/gromacs/gmxpreprocess/grompp.cpp +++ b/src/gromacs/gmxpreprocess/grompp.cpp @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by + * Copyright (c) 2013,2014,2015,2016,2017,2018,2019, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -2245,8 +2245,14 @@ int gmx_grompp(int argc, char *argv[]) if (ir->bDoAwh) { + tensor compressibility = { { 0 } }; + if (ir->epc != epcNO) + { + copy_mat(ir->compress, compressibility); + } setStateDependentAwhParams(ir->awhParams, ir->pull, pull, - state.box, ir->ePBC, &ir->opts, wi); + state.box, ir->ePBC, compressibility, + &ir->opts, wi); } if (ir->bPull) -- 2.22.0