From e38b9130069194bb0d52db1c30ac122110677f80 Mon Sep 17 00:00:00 2001 From: Berk Hess Date: Tue, 3 Oct 2017 13:44:37 +0200 Subject: [PATCH] Add acceleration correction VCM mode Minor refactoring in vcm.cpp using templating to reduce code duplication and improve performance. Change-Id: I560027fe6f315eede0d7aa573bc22bf12eba4c5d --- docs/user-guide/mdp-options.rst | 14 ++- src/gromacs/gmxpreprocess/readir.cpp | 1 + src/gromacs/mdlib/md_support.cpp | 8 +- src/gromacs/mdlib/vcm.cpp | 190 ++++++++++++++++++++++++++--------- src/gromacs/mdlib/vcm.h | 15 ++- src/gromacs/mdtypes/md_enums.cpp | 4 +- src/gromacs/mdtypes/md_enums.h | 4 +- 7 files changed, 174 insertions(+), 62 deletions(-) diff --git a/docs/user-guide/mdp-options.rst b/docs/user-guide/mdp-options.rst index 3feb1d2737..8666f64d34 100644 --- a/docs/user-guide/mdp-options.rst +++ b/docs/user-guide/mdp-options.rst @@ -222,11 +222,21 @@ Run control .. mdp-value:: Linear - Remove center of mass translation + Remove center of mass translational velocity .. mdp-value:: Angular - Remove center of mass translation and rotation around the center of mass + Remove center of mass translational and rotational velocity around + the center of mass + + .. mdp-value:: Linear-acceleration-correction + + Remove center of mass translational velocity. Correct the center of + mass position assuming linear acceleration over :mdp:`nstcomm` steps. + This is useful for cases where an acceleration is expected on the + center of mass which is nearly constant over mdp:`nstcomm` steps. + This can occur for example when pulling on a group using an absolute + reference. .. mdp-value:: None diff --git a/src/gromacs/gmxpreprocess/readir.cpp b/src/gromacs/gmxpreprocess/readir.cpp index c213b455a9..1f0f2a054c 100644 --- a/src/gromacs/gmxpreprocess/readir.cpp +++ b/src/gromacs/gmxpreprocess/readir.cpp @@ -3018,6 +3018,7 @@ static void calc_nrdf(gmx_mtop_t *mtop, t_inputrec *ir, char **gnames) switch (ir->comm_mode) { case ecmLINEAR: + case ecmLINEAR_ACCELERATION_CORRECTION: nrdf_vcm_sub[j] = 0; for (d = 0; d < ndim_rm_vcm; d++) { diff --git a/src/gromacs/mdlib/md_support.cpp b/src/gromacs/mdlib/md_support.cpp index 3c2f23b81d..9115a942a5 100644 --- a/src/gromacs/mdlib/md_support.cpp +++ b/src/gromacs/mdlib/md_support.cpp @@ -247,8 +247,12 @@ void compute_globals(FILE *fplog, gmx_global_stat *gstat, t_commrec *cr, t_input if (bStopCM) { check_cm_grp(fplog, vcm, ir, 1); - do_stopcm_grp(0, mdatoms->homenr, mdatoms->cVCM, - as_rvec_array(state->x.data()), as_rvec_array(state->v.data()), vcm); + /* Don't pass x with linear modes to avoid correction of the initial + * coordinates for the initial COM velocity. + */ + do_stopcm_grp(mdatoms->homenr, mdatoms->cVCM, + vcm->mode == ecmANGULAR ? as_rvec_array(state->x.data()) : nullptr, + as_rvec_array(state->v.data()), *vcm); inc_nrnb(nrnb, eNR_STOPCM, mdatoms->homenr); } diff --git a/src/gromacs/mdlib/vcm.cpp b/src/gromacs/mdlib/vcm.cpp index 079c0342f7..16a0c2b412 100644 --- a/src/gromacs/mdlib/vcm.cpp +++ b/src/gromacs/mdlib/vcm.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, by the GROMACS development team, led by + * Copyright (c) 2013,2014,2015,2016,2017, 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. @@ -50,18 +50,20 @@ #include "gromacs/pbcutil/pbc.h" #include "gromacs/topology/topology.h" #include "gromacs/utility/fatalerror.h" +#include "gromacs/utility/gmxassert.h" #include "gromacs/utility/gmxomp.h" #include "gromacs/utility/smalloc.h" -t_vcm *init_vcm(FILE *fp, gmx_groups_t *groups, t_inputrec *ir) +t_vcm *init_vcm(FILE *fp, gmx_groups_t *groups, const t_inputrec *ir) { t_vcm *vcm; int g; snew(vcm, 1); - vcm->mode = (ir->nstcomm > 0) ? ir->comm_mode : ecmNO; - vcm->ndim = ndof_com(ir); + vcm->mode = (ir->nstcomm > 0) ? ir->comm_mode : ecmNO; + vcm->ndim = ndof_com(ir); + vcm->timeStep = ir->nstcomm*ir->delta_t; if (vcm->mode == ecmANGULAR && vcm->ndim < 3) { @@ -227,68 +229,158 @@ void calc_vcm_grp(int start, int homenr, t_mdatoms *md, } } -void do_stopcm_grp(int start, int homenr, unsigned short *group_id, - rvec x[], rvec v[], t_vcm *vcm) +/*! \brief Remove the COM motion velocity from the velocities + * + * \note This routine should be called from within an OpenMP parallel region. + * + * \tparam numDimensions Correct dimensions 0 to \p numDimensions-1 + * \param[in] homenr The number of atoms to correct + * \param[in] group_id List of VCM group ids, when nullptr is passed all atoms are assumed to be in group 0 + * \param[in,out] v The velocities to correct + * \param[in] vcm VCM data + */ +template +static void +doStopComMotionLinear(int homenr, + const unsigned short *group_id, + rvec *v, + const t_vcm &vcm) { - if (vcm->mode != ecmNO) + if (group_id == nullptr) { - // cppcheck-suppress unreadVariable - int gmx_unused nth = gmx_omp_nthreads_get(emntDefault); -#pragma omp parallel num_threads(nth) +#pragma omp for schedule(static) + for (int i = 0; i < homenr; i++) { - int i, g = 0; - rvec dv, dx; - /* Subtract linear momentum */ - switch (vcm->ndim) + for (int d = 0; d < numDimensions; d++) { - case 1: + v[i][d] -= vcm.group_v[0][d]; + } + } + } + else + { #pragma omp for schedule(static) - for (i = start; i < start+homenr; i++) - { - if (group_id) - { - g = group_id[i]; - } - v[i][XX] -= vcm->group_v[g][XX]; - } - break; - case 2: + for (int i = 0; i < homenr; i++) + { + const int g = group_id[i]; + for (int d = 0; d < numDimensions; d++) + { + v[i][d] -= vcm.group_v[g][d]; + } + } + } +} + +/*! \brief Remove the COM motion velocity from the velocities, correct the coordinates assuming constant acceleration + * + * \note This routine should be called from within an OpenMP parallel region. + * + * \tparam numDimensions Correct dimensions 0 to \p numDimensions-1 + * \param[in] homenr The number of atoms to correct + * \param[in] group_id List of VCM group ids, when nullptr is passed all atoms are assumed to be in group 0 + * \param[in,out] x The coordinates to correct + * \param[in,out] v The velocities to correct + * \param[in] vcm VCM data + */ +template +static void +doStopComMotionAccelerationCorrection(int homenr, + const unsigned short *group_id, + rvec * gmx_restrict x, + rvec * gmx_restrict v, + const t_vcm &vcm) +{ + const real xCorrectionFactor = 0.5*vcm.timeStep; + + if (group_id == nullptr) + { #pragma omp for schedule(static) - for (i = start; i < start+homenr; i++) - { - if (group_id) - { - g = group_id[i]; - } - v[i][XX] -= vcm->group_v[g][XX]; - v[i][YY] -= vcm->group_v[g][YY]; - } - break; - case 3: + for (int i = 0; i < homenr; i++) + { + for (int d = 0; d < numDimensions; d++) + { + x[i][d] -= vcm.group_v[0][d]*xCorrectionFactor; + v[i][d] -= vcm.group_v[0][d]; + } + } + } + else + { #pragma omp for schedule(static) - for (i = start; i < start+homenr; i++) - { - if (group_id) - { - g = group_id[i]; - } - rvec_dec(v[i], vcm->group_v[g]); - } - break; + for (int i = 0; i < homenr; i++) + { + const int g = group_id[i]; + for (int d = 0; d < numDimensions; d++) + { + x[i][d] -= vcm.group_v[g][d]*xCorrectionFactor; + v[i][d] -= vcm.group_v[g][d]; } - if (vcm->mode == ecmANGULAR) + } + } +} + +void do_stopcm_grp(int homenr, const unsigned short *group_id, + rvec x[], rvec v[], const t_vcm &vcm) +{ + if (vcm.mode != ecmNO) + { + // cppcheck-suppress unreadVariable + int gmx_unused nth = gmx_omp_nthreads_get(emntDefault); +#pragma omp parallel num_threads(nth) + { + if (vcm.mode == ecmLINEAR || + vcm.mode == ecmANGULAR || + (vcm.mode == ecmLINEAR_ACCELERATION_CORRECTION && x == nullptr)) + { + /* Subtract linear momentum for v */ + switch (vcm.ndim) + { + case 1: + doStopComMotionLinear<1>(homenr, group_id, v, vcm); + break; + case 2: + doStopComMotionLinear<2>(homenr, group_id, v, vcm); + break; + case 3: + doStopComMotionLinear<3>(homenr, group_id, v, vcm); + break; + } + } + else + { + GMX_ASSERT(vcm.mode == ecmLINEAR_ACCELERATION_CORRECTION, "When the mode is not linear or angular, it should be acceleration correction"); + /* Subtract linear momentum for v and x*/ + switch (vcm.ndim) + { + case 1: + doStopComMotionAccelerationCorrection<1>(homenr, group_id, x, v, vcm); + break; + case 2: + doStopComMotionAccelerationCorrection<2>(homenr, group_id, x, v, vcm); + break; + case 3: + doStopComMotionAccelerationCorrection<3>(homenr, group_id, x, v, vcm); + break; + } + + } + if (vcm.mode == ecmANGULAR) { /* Subtract angular momentum */ + GMX_ASSERT(x, "Need x to compute angular momentum correction"); + + int g = 0; #pragma omp for schedule(static) - for (i = start; i < start+homenr; i++) + for (int i = 0; i < homenr; i++) { if (group_id) { g = group_id[i]; } /* Compute the correction to the velocity for each atom */ - rvec_sub(x[i], vcm->group_x[g], dx); - cprod(vcm->group_w[g], dx, dv); + rvec dv, dx; + rvec_sub(x[i], vcm.group_x[g], dx); + cprod(vcm.group_w[g], dx, dv); rvec_dec(v[i], dv); } } diff --git a/src/gromacs/mdlib/vcm.h b/src/gromacs/mdlib/vcm.h index 4eebf427b6..7ffe477d82 100644 --- a/src/gromacs/mdlib/vcm.h +++ b/src/gromacs/mdlib/vcm.h @@ -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, by the GROMACS development team, led by + * Copyright (c) 2013,2014,2015,2016,2017, 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. @@ -61,6 +61,7 @@ typedef struct { int stride; /* Stride for thread data */ int mode; /* One of the enums above */ gmx_bool ndim; /* The number of dimensions for corr. */ + real timeStep; /* The time step for COMM removal */ real *group_ndf; /* Number of degrees of freedom */ rvec *group_p; /* Linear momentum per group */ rvec *group_v; /* Linear velocity per group */ @@ -73,15 +74,19 @@ typedef struct { t_vcm_thread *thread_vcm; /* Temporary data per thread and group */ } t_vcm; -t_vcm *init_vcm(FILE *fp, gmx_groups_t *groups, t_inputrec *ir); +t_vcm *init_vcm(FILE *fp, gmx_groups_t *groups, const t_inputrec *ir); /* Do a per group center of mass things */ void calc_vcm_grp(int start, int homenr, t_mdatoms *md, rvec x[], rvec v[], t_vcm *vcm); -void do_stopcm_grp(int start, int homenr, - unsigned short *group_id, - rvec x[], rvec v[], t_vcm *vcm); +/* Set the COM velocity to zero and potentially correct the COM position. + * + * With linear modes nullptr can be passed for x. + */ +void do_stopcm_grp(int homenr, + const unsigned short *group_id, + rvec x[], rvec v[], const t_vcm &vcm); void check_cm_grp(FILE *fp, t_vcm *vcm, t_inputrec *ir, real Temp_Max); diff --git a/src/gromacs/mdtypes/md_enums.cpp b/src/gromacs/mdtypes/md_enums.cpp index 2b6eb86c62..e47f1e0a0f 100644 --- a/src/gromacs/mdtypes/md_enums.cpp +++ b/src/gromacs/mdtypes/md_enums.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, by the GROMACS development team, led by + * Copyright (c) 2013,2014,2015,2016,2017, 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. @@ -182,7 +182,7 @@ const char *edispc_names[edispcNR+1] = { }; const char *ecm_names[ecmNR+1] = { - "Linear", "Angular", "None", nullptr + "Linear", "Angular", "None", "Linear-acceleration-correction", nullptr }; const char *eann_names[eannNR+1] = { diff --git a/src/gromacs/mdtypes/md_enums.h b/src/gromacs/mdtypes/md_enums.h index e2a016e48d..a0c1565499 100644 --- a/src/gromacs/mdtypes/md_enums.h +++ b/src/gromacs/mdtypes/md_enums.h @@ -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, by the GROMACS development team, led by + * Copyright (c) 2013,2014,2015,2016,2017, 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. @@ -453,7 +453,7 @@ extern const char *edispc_names[edispcNR+1]; //! Center of mass motion removal algorithm. enum { - ecmLINEAR, ecmANGULAR, ecmNO, ecmNR + ecmLINEAR, ecmANGULAR, ecmNO, ecmLINEAR_ACCELERATION_CORRECTION, ecmNR }; //! String corresponding to COM removal extern const char *ecm_names[ecmNR+1]; -- 2.11.4.GIT