From 62e43da89cc33b842e51f6933a8f38e68feb1bec Mon Sep 17 00:00:00 2001 From: Berk Hess Date: Tue, 19 Dec 2017 11:55:25 +0100 Subject: [PATCH] Make acceleration correction VCM mode work The new acceleration correction VCM mode did not actually correct the coordinate for the acceleration, since a null pointer was passed. Introduced an extra CGLO flag to allow for correction of the coordinates, but leave the initial coordinates unaffected. Change-Id: I673793902df7241085fff20c63cf3ce88ef60313 --- src/gromacs/mdlib/md_support.cpp | 12 ++++++++---- src/gromacs/mdlib/md_support.h | 2 ++ src/programs/mdrun/md.cpp | 2 +- 3 files changed, 11 insertions(+), 5 deletions(-) diff --git a/src/gromacs/mdlib/md_support.cpp b/src/gromacs/mdlib/md_support.cpp index 9115a942a5..72d90484db 100644 --- a/src/gromacs/mdlib/md_support.cpp +++ b/src/gromacs/mdlib/md_support.cpp @@ -247,12 +247,16 @@ void compute_globals(FILE *fplog, gmx_global_stat *gstat, t_commrec *cr, t_input if (bStopCM) { check_cm_grp(fplog, vcm, ir, 1); - /* Don't pass x with linear modes to avoid correction of the initial - * coordinates for the initial COM velocity. + /* At initialization, do not pass x with acceleration-correction mode + * to avoid (incorrect) correction of the initial coordinates. */ + rvec *xPtr = nullptr; + if (vcm->mode == ecmANGULAR || (vcm->mode == ecmLINEAR_ACCELERATION_CORRECTION && !(flags & CGLO_INITIALIZATION))) + { + xPtr = as_rvec_array(state->x.data()); + } do_stopcm_grp(mdatoms->homenr, mdatoms->cVCM, - vcm->mode == ecmANGULAR ? as_rvec_array(state->x.data()) : nullptr, - as_rvec_array(state->v.data()), *vcm); + xPtr, as_rvec_array(state->v.data()), *vcm); inc_nrnb(nrnb, eNR_STOPCM, mdatoms->homenr); } diff --git a/src/gromacs/mdlib/md_support.h b/src/gromacs/mdlib/md_support.h index 20bbcc8b0b..ddee3c6792 100644 --- a/src/gromacs/mdlib/md_support.h +++ b/src/gromacs/mdlib/md_support.h @@ -64,6 +64,8 @@ class SimulationSignaller; * passed to compute_globals in md.c and global_stat. */ +/* we are initializing and not yet in the actual MD loop */ +#define CGLO_INITIALIZATION (1<<1) /* we are computing the kinetic energy from average velocities */ #define CGLO_EKINAVEVEL (1<<2) /* we are removing the center of mass momenta */ diff --git a/src/programs/mdrun/md.cpp b/src/programs/mdrun/md.cpp index 9604388c07..683460fa6b 100644 --- a/src/programs/mdrun/md.cpp +++ b/src/programs/mdrun/md.cpp @@ -704,7 +704,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, const gmx::MDLogger &mdlog, restore_ekinstate_from_state(cr, ekind, &state_global->ekinstate); } - cglo_flags = (CGLO_TEMPERATURE | CGLO_GSTAT + cglo_flags = (CGLO_INITIALIZATION | CGLO_TEMPERATURE | CGLO_GSTAT | (EI_VV(ir->eI) ? CGLO_PRESSURE : 0) | (EI_VV(ir->eI) ? CGLO_CONSTRAINT : 0) | (continuationOptions.haveReadEkin ? CGLO_READEKIN : 0)); -- 2.11.4.GIT