From: Magnus Lundborg Date: Thu, 5 Jul 2018 11:59:46 +0000 (+0200) Subject: Allow using COM of previous step as PBC reference X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=b59c3db50186e7d96e97c8c620d08caf628e6a84;p=alexxy%2Fgromacs.git Allow using COM of previous step as PBC reference Add an option, when pulling, to use the COM of the group of the previous step, to calculate PBC jumps, instead of a reference atom, which can sometimes move a lot during the simulation. When there is no previous step (when the COM of the previous step is not set) use the COM based on the reference atom. The COM of the previous step is written to the checkpoint. Fixes #2625 Change-Id: I6120b76a15f92167753463326125428f822848ab --- diff --git a/docs/reference-manual/special/pulling.rst b/docs/reference-manual/special/pulling.rst index 525da7a221..2a08303fbc 100644 --- a/docs/reference-manual/special/pulling.rst +++ b/docs/reference-manual/special/pulling.rst @@ -86,6 +86,16 @@ default, the middle (determined by the order in the topology) atom is used as a reference atom, but the user can also select any other atom if it would be closer to center of the group. +When there are large pull groups, such as a +lipid bilayer, ``pull-pbc-ref-prev-step-com`` can be used to avoid potential +large movements of the center of mass in case that atoms in the pull group +move so much that the reference atom is too far from the intended center of mass. +With this option enabled the center of mass from the previous step is used, +instead of the position of the reference atom, to determine the reference position. +The position of the reference atom is still used for the first step. For large pull +groups it is important to select a reference atom that is close to the intended +center of mass, i.e. do not use ``pull-group?-pbcatom = 0``. + For a layered system, for instance a lipid bilayer, it may be of interest to calculate the PMF of a lipid as function of its distance from the whole bilayer. The whole bilayer can be taken as reference diff --git a/docs/release-notes/features.rst b/docs/release-notes/features.rst index 866cf70350..6a432b0da7 100644 --- a/docs/release-notes/features.rst +++ b/docs/release-notes/features.rst @@ -26,3 +26,13 @@ On Intel CPUs with integrated GPUs, it is now possible to offload nonbonded task to the GPU the same way as offload is done to other GPU architectures. This can have performance benefits, in particular on modern desktop and mobile Intel CPUs this offload can give up to 20% higher simulation performance. + +Allow using COM of previous step as PBC reference +""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" + +Added an option (``pull-pbc-ref-from-prev-step-com``), when pulling, to use +the COM of the group of the previous step, to calculate PBC jumps instead of a +reference atom, which can sometimes move a lot during the simulation. +With this option the PBC reference atom is only used at initialization. +This can be of use when using large pull groups or groups with potentially +large relative movement of atoms. diff --git a/docs/user-guide/mdp-options.rst b/docs/user-guide/mdp-options.rst index cb93e10ad5..02907c6581 100644 --- a/docs/user-guide/mdp-options.rst +++ b/docs/user-guide/mdp-options.rst @@ -1612,6 +1612,20 @@ pull-coord2-vec, pull-coord2-k, and so on. frequency for writing out the force of all the pulled group (0 is never) +.. mdp:: pull-pbc-ref-prev-step-com + + .. mdp-value:: no + + Use the reference atom (:mdp:`pull-group1-pbcatom`) for the + treatment of periodic boundary conditions. + + .. mdp-value:: yes + + Use the COM of the previous step as reference for the treatment + of periodic boundary conditions. The reference is initialized + using the reference atom (:mdp:`pull-group1-pbcatom`), which should + be located centrally in the group. Using the COM from the + previous step can be useful if one or more pull groups are large. .. mdp:: pull-ngroups @@ -1649,7 +1663,10 @@ pull-coord2-vec, pull-coord2-k, and so on. vector. For determining the COM, all atoms in the group are put at their periodic image which is closest to :mdp:`pull-group1-pbcatom`. A value of 0 means that the middle - atom (number wise) is used. This parameter is not used with + atom (number wise) is used, which is only safe for small groups. + :ref:`gmx grompp` checks that the maximum distance from the reference + atom (specifically chosen, or not) to the other atoms in the group + is not too large. This parameter is not used with :mdp:`pull-coord1-geometry` cylinder. A value of -1 turns on cosine weighting, which is useful for a group of molecules in a periodic system, *e.g.* a water slab (see Engin et al. J. Chem. Phys. B diff --git a/src/gromacs/fileio/checkpoint.cpp b/src/gromacs/fileio/checkpoint.cpp index 620df7b6b9..38f76ee8ed 100644 --- a/src/gromacs/fileio/checkpoint.cpp +++ b/src/gromacs/fileio/checkpoint.cpp @@ -109,6 +109,7 @@ enum cptv { cptv_Unknown = 17, /**< Version before numbering scheme */ cptv_RemoveBuildMachineInformation, /**< remove functionality that makes mdrun builds non-reproducible */ + cptv_ComPrevStepAsPullGroupReference, /**< Allow using COM of previous step as pull group PBC reference */ cptv_Count /**< the total number of cptv versions */ }; @@ -213,6 +214,17 @@ static const char *eawhh_names[eawhhNR] = "awh_forceCorrelationGrid" }; +enum { + epullsPREVSTEPCOM, + epullsNR +}; + +static const char *epull_prev_step_com_names[epullsNR] = +{ + "Pull groups prev step COM" +}; + + //! Higher level vector element type, only used for formatting checkpoint dumps enum class CptElementType { @@ -229,7 +241,8 @@ enum class StatePart kineticEnergy, //!< Kinetic energy, needed for T/P-coupling state energyHistory, //!< Energy observable statistics freeEnergyHistory, //!< Free-energy state and observable statistics - accWeightHistogram //!< Accelerated weight histogram method state + accWeightHistogram, //!< Accelerated weight histogram method state + pullState //!< COM of previous step. }; //! \brief Return the name of a checkpoint entry based on part and part entry @@ -242,6 +255,7 @@ static const char *entryName(StatePart part, int ecpt) case StatePart::energyHistory: return eenh_names[ecpt]; case StatePart::freeEnergyHistory: return edfh_names[ecpt]; case StatePart::accWeightHistogram: return eawhh_names[ecpt]; + case StatePart::pullState: return epull_prev_step_com_names[ecpt]; } return nullptr; @@ -1177,13 +1191,13 @@ static int do_cpt_state(XDR *xd, case estDISRE_RM3TAV: ret = do_cpte_n_reals(xd, part, i, sflags, &state->hist.ndisrepairs, &state->hist.disre_rm3tav, list); break; case estORIRE_INITF: ret = do_cpte_real (xd, part, i, sflags, &state->hist.orire_initf, list); break; case estORIRE_DTAV: ret = do_cpte_n_reals(xd, part, i, sflags, &state->hist.norire_Dtav, &state->hist.orire_Dtav, list); break; + case estPREVSTEPCOM: ret = doVector(xd, part, i, sflags, &state->com_prev_step, list); break; default: gmx_fatal(FARGS, "Unknown state entry %d\n" "You are reading a checkpoint file written by different code, which is not supported", i); } } } - return ret; } diff --git a/src/gromacs/fileio/tpxio.cpp b/src/gromacs/fileio/tpxio.cpp index 4caaff0b27..2fcaaba4f6 100644 --- a/src/gromacs/fileio/tpxio.cpp +++ b/src/gromacs/fileio/tpxio.cpp @@ -120,6 +120,7 @@ enum tpxv { tpxv_GenericParamsForElectricField, /**< Introduced KeyValueTree and moved electric field parameters */ tpxv_AcceleratedWeightHistogram, /**< sampling with accelerated weight histogram method (AWH) */ tpxv_RemoveImplicitSolvation, /**< removed support for implicit solvation */ + tpxv_PullPrevStepCOMAsReference, /**< Enabled using the COM of the pull group of the last frame as reference for PBC */ tpxv_Count /**< the total number of tpxv versions */ }; @@ -746,6 +747,14 @@ static void do_pull(t_fileio *fio, pull_params_t *pull, gmx_bool bRead, } gmx_fio_do_int(fio, pull->nstxout); gmx_fio_do_int(fio, pull->nstfout); + if (file_version >= tpxv_PullPrevStepCOMAsReference) + { + gmx_fio_do_gmx_bool(fio, pull->bSetPbcRefToPrevStepCOM); + } + else + { + pull->bSetPbcRefToPrevStepCOM = FALSE; + } if (bRead) { snew(pull->group, pull->ngroup); diff --git a/src/gromacs/gmxpreprocess/readpull.cpp b/src/gromacs/gmxpreprocess/readpull.cpp index faa4c4af13..943560aaf4 100644 --- a/src/gromacs/gmxpreprocess/readpull.cpp +++ b/src/gromacs/gmxpreprocess/readpull.cpp @@ -279,13 +279,14 @@ char **read_pullparams(std::vector *inp, /* read pull parameters */ printStringNoNewline(inp, "Cylinder radius for dynamic reaction force groups (nm)"); - pull->cylinder_r = get_ereal(inp, "pull-cylinder-r", 1.5, wi); - pull->constr_tol = get_ereal(inp, "pull-constr-tol", 1E-6, wi); - pull->bPrintCOM = (get_eeenum(inp, "pull-print-com", yesno_names, wi) != 0); - pull->bPrintRefValue = (get_eeenum(inp, "pull-print-ref-value", yesno_names, wi) != 0); - pull->bPrintComp = (get_eeenum(inp, "pull-print-components", yesno_names, wi) != 0); - pull->nstxout = get_eint(inp, "pull-nstxout", 50, wi); - pull->nstfout = get_eint(inp, "pull-nstfout", 50, wi); + pull->cylinder_r = get_ereal(inp, "pull-cylinder-r", 1.5, wi); + pull->constr_tol = get_ereal(inp, "pull-constr-tol", 1E-6, wi); + pull->bPrintCOM = (get_eeenum(inp, "pull-print-com", yesno_names, wi) != 0); + pull->bPrintRefValue = (get_eeenum(inp, "pull-print-ref-value", yesno_names, wi) != 0); + pull->bPrintComp = (get_eeenum(inp, "pull-print-components", yesno_names, wi) != 0); + pull->nstxout = get_eint(inp, "pull-nstxout", 50, wi); + pull->nstfout = get_eint(inp, "pull-nstfout", 50, wi); + pull->bSetPbcRefToPrevStepCOM = (get_eeenum(inp, "pull-pbc-ref-prev-step-com", yesno_names, wi) != 0); printStringNoNewline(inp, "Number of pull groups"); pull->ngroup = get_eint(inp, "pull-ngroups", 1, wi); printStringNoNewline(inp, "Number of pull coordinates"); @@ -405,13 +406,15 @@ void make_pull_groups(pull_params_t *pull, t_pull_group *pgrp; /* Absolute reference group (might not be used) is special */ - pgrp = &pull->group[0]; - pgrp->nat = 0; - pgrp->pbcatom = -1; + pgrp = &pull->group[0]; + pgrp->nat = 0; + pgrp->pbcatom = -1; + pgrp->pbcatom_input = -1; for (g = 1; g < pull->ngroup; g++) { - pgrp = &pull->group[g]; + pgrp = &pull->group[g]; + pgrp->pbcatom_input = pgrp->pbcatom; if (strcmp(pgnames[g], "") == 0) { @@ -520,19 +523,57 @@ pull_t *set_pull_init(t_inputrec *ir, const gmx_mtop_t *mtop, t_start = ir->init_t + ir->init_step*ir->delta_t; + if (pull->bSetPbcRefToPrevStepCOM) + { + initPullComFromPrevStep(nullptr, pull_work, md, &pbc, x); + } pull_calc_coms(nullptr, pull_work, md, &pbc, t_start, x, nullptr); - int groupThatFailsPbc = pullCheckPbcWithinGroups(*pull_work, x, pbc, c_pullGroupPbcMargin); - if (groupThatFailsPbc >= 0) + for (int g = 0; g < pull->ngroup; g++) { - char buf[STRLEN]; - sprintf(buf, - "Pull group %d has atoms at a distance larger than %g times half the box size from the PBC atom (%d). If atoms are or will more beyond half the box size from the PBC atom, the COM will be ill defined.", - groupThatFailsPbc, - c_pullGroupPbcMargin, - pull->group[groupThatFailsPbc].pbcatom); - set_warning_line(wi, nullptr, -1); - warning(wi, buf); + bool groupObeysPbc = + pullCheckPbcWithinGroup(*pull_work, + gmx::arrayRefFromArray(reinterpret_cast(x), + mtop->natoms), + pbc, g, c_pullGroupSmallGroupThreshold); + if (!groupObeysPbc) + { + char buf[STRLEN]; + if (pull->group[g].pbcatom_input == 0) + { + sprintf(buf, "When the maximum distance from a pull group reference atom to other atoms in the " + "group is larger than %g times half the box size a centrally placed " + "atom should be chosen as pbcatom. Pull group %d is larger than that and does not have " + "a specific atom selected as reference atom.", c_pullGroupSmallGroupThreshold, g); + warning_error(wi, buf); + } + else if (!pull->bSetPbcRefToPrevStepCOM) + { + sprintf(buf, "The maximum distance from the chosen PBC atom (%d) of pull group %d to other " + "atoms in the group is larger than %g times half the box size. " + "Set the pull-pbc-ref-prev-step-com option to yes.", pull->group[g].pbcatom + 1, + g, c_pullGroupSmallGroupThreshold); + warning_error(wi, buf); + } + } + if (groupObeysPbc) + { + groupObeysPbc = + pullCheckPbcWithinGroup(*pull_work, + gmx::arrayRefFromArray(reinterpret_cast(x), + mtop->natoms), + pbc, g, c_pullGroupPbcMargin); + if (!groupObeysPbc) + { + char buf[STRLEN]; + sprintf(buf, + "Pull group %d has atoms at a distance larger than %g times half the box size from the PBC atom (%d). " + "If atoms are or will more beyond half the box size from the PBC atom, the COM will be ill defined.", + g, c_pullGroupPbcMargin, pull->group[g].pbcatom + 1); + set_warning_line(wi, nullptr, -1); + warning(wi, buf); + } + } } fprintf(stderr, "Pull group natoms pbc atom distance at start reference at t=0\n"); diff --git a/src/gromacs/mdlib/md_support.cpp b/src/gromacs/mdlib/md_support.cpp index b64b8185d7..1285af00fd 100644 --- a/src/gromacs/mdlib/md_support.cpp +++ b/src/gromacs/mdlib/md_support.cpp @@ -63,6 +63,7 @@ #include "gromacs/mdtypes/md_enums.h" #include "gromacs/mdtypes/state.h" #include "gromacs/pbcutil/pbc.h" +#include "gromacs/pulling/pull.h" #include "gromacs/timing/wallcycle.h" #include "gromacs/topology/mtop_util.h" #include "gromacs/trajectory/trajectoryframe.h" @@ -618,4 +619,9 @@ void set_state_entries(t_state *state, const t_inputrec *ir) snew(state->dfhist, 1); init_df_history(state->dfhist, ir->fepvals->n_lambda); } + + if (ir->pull && ir->pull->bSetPbcRefToPrevStepCOM) + { + state->flags |= (1<energyHistory = {}; } + if (ir->pull && ir->pull->bSetPbcRefToPrevStepCOM) + { + /* Copy the pull group COM of the previous step from the checkpoint state to the pull state */ + setPrevStepPullComFromState(ir->pull_work, state); + } + } + else if (ir->pull && ir->pull->bSetPbcRefToPrevStepCOM) + { + allocStatePrevStepPullCom(state, ir->pull_work); + t_pbc pbc; + set_pbc(&pbc, ir->ePBC, state->box); + initPullComFromPrevStep(cr, ir->pull_work, mdatoms, &pbc, as_rvec_array(state->x.data())); + updatePrevStepCom(ir->pull_work); + setStatePrevStepPullCom(ir->pull_work, state); } if (observablesHistory->energyHistory == nullptr) { @@ -1145,6 +1159,12 @@ void gmx::Integrator::do_md() state, graph, nrnb, wcycle, upd, constr); + if (MASTER(cr) && ir->bPull && ir->pull->bSetPbcRefToPrevStepCOM) + { + updatePrevStepCom(ir->pull_work); + setStatePrevStepPullCom(ir->pull_work, state); + } + if (ir->eI == eiVVAK) { /* erase F_EKIN and F_TEMP here? */ diff --git a/src/gromacs/mdtypes/inputrec.cpp b/src/gromacs/mdtypes/inputrec.cpp index 28aeecd0bd..7318d397b6 100644 --- a/src/gromacs/mdtypes/inputrec.cpp +++ b/src/gromacs/mdtypes/inputrec.cpp @@ -642,6 +642,7 @@ static void pr_pull(FILE *fp, int indent, const pull_params_t *pull) PS("pull-print-components", EBOOL(pull->bPrintComp)); PI("pull-nstxout", pull->nstxout); PI("pull-nstfout", pull->nstfout); + PS("pull-pbc-ref-prev-step-com", EBOOL(pull->bSetPbcRefToPrevStepCOM)); PI("pull-ngroups", pull->ngroup); for (g = 0; g < pull->ngroup; g++) { diff --git a/src/gromacs/mdtypes/pull-params.h b/src/gromacs/mdtypes/pull-params.h index 4de146229c..57d8aa1c47 100644 --- a/src/gromacs/mdtypes/pull-params.h +++ b/src/gromacs/mdtypes/pull-params.h @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2015,2016, by the GROMACS development team, led by + * Copyright (c) 2015,2016,2018, 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. @@ -57,11 +57,12 @@ /*! \brief Struct that defines a pull group */ typedef struct { - int nat; /**< Number of atoms in the pull group */ - int *ind; /**< The global atoms numbers */ - int nweight; /**< The number of weights (0 or nat) */ - real *weight; /**< Weights (use all 1 when weight==NULL) */ - int pbcatom; /**< The reference atom for pbc (global number) */ + int nat; /**< Number of atoms in the pull group */ + int *ind; /**< The global atoms numbers */ + int nweight; /**< The number of weights (0 or nat) */ + real *weight; /**< Weights (use all 1 when weight==NULL) */ + int pbcatom; /**< The reference atom for pbc (global number) */ + int pbcatom_input; /**< The reference atom for pbc (global number) as specified in the input parameters */ } t_pull_group; /*! Maximum number of pull groups that can be used in a pull coordinate */ @@ -86,18 +87,19 @@ typedef struct { /*! \brief Struct containing all pull parameters */ typedef struct pull_params_t { - int ngroup; /**< Number of pull groups */ - int ncoord; /**< Number of pull coordinates */ - real cylinder_r; /**< Radius of cylinder for dynamic COM (nm) */ - real constr_tol; /**< Absolute tolerance for constraints in (nm) */ - gmx_bool bPrintCOM; /**< Print coordinates of COM for each coord */ - gmx_bool bPrintRefValue; /**< Print the reference value for each coord */ - gmx_bool bPrintComp; /**< Print cartesian components for each coord with geometry=distance */ - int nstxout; /**< Output interval for pull x */ - int nstfout; /**< Output interval for pull f */ + int ngroup; /**< Number of pull groups */ + int ncoord; /**< Number of pull coordinates */ + real cylinder_r; /**< Radius of cylinder for dynamic COM (nm) */ + real constr_tol; /**< Absolute tolerance for constraints in (nm) */ + gmx_bool bPrintCOM; /**< Print coordinates of COM for each coord */ + gmx_bool bPrintRefValue; /**< Print the reference value for each coord */ + gmx_bool bPrintComp; /**< Print cartesian components for each coord with geometry=distance */ + gmx_bool bSetPbcRefToPrevStepCOM; /**< Use the COM of each group from the previous step as reference */ + int nstxout; /**< Output interval for pull x */ + int nstfout; /**< Output interval for pull f */ - t_pull_group *group; /**< groups to pull/restrain/etc/ */ - t_pull_coord *coord; /**< the pull coordinates */ + t_pull_group *group; /**< groups to pull/restrain/etc/ */ + t_pull_coord *coord; /**< the pull coordinates */ } pull_params_t; /*! \endcond */ diff --git a/src/gromacs/mdtypes/state.cpp b/src/gromacs/mdtypes/state.cpp index d8d1716ab0..d857729727 100644 --- a/src/gromacs/mdtypes/state.cpp +++ b/src/gromacs/mdtypes/state.cpp @@ -50,6 +50,7 @@ #include "gromacs/mdtypes/df_history.h" #include "gromacs/mdtypes/inputrec.h" #include "gromacs/mdtypes/md_enums.h" +#include "gromacs/mdtypes/pull-params.h" #include "gromacs/mdtypes/swaphistory.h" #include "gromacs/pbcutil/boxutilities.h" #include "gromacs/pbcutil/pbc.h" diff --git a/src/gromacs/mdtypes/state.h b/src/gromacs/mdtypes/state.h index 22a3c8c6b2..4e8705517f 100644 --- a/src/gromacs/mdtypes/state.h +++ b/src/gromacs/mdtypes/state.h @@ -94,7 +94,7 @@ enum { estORIRE_INITF, estORIRE_DTAV, estSVIR_PREV, estNH_VXI, estVETA, estVOL0, estNHPRES_XI, estNHPRES_VXI, estFVIR_PREV, estFEPSTATE, estMC_RNG_NOTSUPPORTED, estMC_RNGI_NOTSUPPORTED, - estBAROS_INT, + estBAROS_INT, estPREVSTEPCOM, estNR }; @@ -224,6 +224,8 @@ class t_state int ddp_count; //!< The DD partitioning count for this state int ddp_count_cg_gl; //!< The DD partitioning count for index_gl std::vector cg_gl; //!< The global cg number of the local cgs + + std::vector com_prev_step; //!< The COM of the previous step of each pull group }; #ifndef DOXYGEN diff --git a/src/gromacs/pulling/pull.cpp b/src/gromacs/pulling/pull.cpp index 250f92d3b3..f4134ca573 100644 --- a/src/gromacs/pulling/pull.cpp +++ b/src/gromacs/pulling/pull.cpp @@ -84,7 +84,7 @@ #include "pull_internal.h" -static int groupPbcFromParams(const t_pull_group ¶ms) +static int groupPbcFromParams(const t_pull_group ¶ms, bool setPbcRefToPrevStepCOM) { if (params.nat <= 1) { @@ -93,7 +93,14 @@ static int groupPbcFromParams(const t_pull_group ¶ms) } else if (params.pbcatom >= 0) { - return epgrppbcREFAT; + if (setPbcRefToPrevStepCOM) + { + return epgrppbcPREVSTEPCOM; + } + else + { + return epgrppbcREFAT; + } } else { @@ -107,9 +114,10 @@ static int groupPbcFromParams(const t_pull_group ¶ms) * This will be fixed when the pointers are replacted by std::vector. */ pull_group_work_t::pull_group_work_t(const t_pull_group ¶ms, - gmx::LocalAtomSet atomSet) : + gmx::LocalAtomSet atomSet, + bool bSetPbcRefToPrevStepCOM) : params(params), - epgrppbc(groupPbcFromParams(params)), + epgrppbc(groupPbcFromParams(params, bSetPbcRefToPrevStepCOM)), needToCalcCom(false), atomSet(atomSet), mwscale(0), @@ -1861,9 +1869,8 @@ void dd_make_local_pull_groups(const t_commrec *cr, struct pull_t *pull) /* We should participate if we have pull or pbc atoms */ if (!bMustParticipate && - (group.atomSet.numAtomsLocal() > 0 || - (group.epgrppbc == epgrppbcREFAT && - group.pbcAtomSet->numAtomsLocal() > 0))) + (group.atomSet.numAtomsLocal() > 0 || (group.epgrppbc == epgrppbcREFAT && + group.pbcAtomSet->numAtomsLocal() > 0))) { bMustParticipate = TRUE; } @@ -2120,7 +2127,8 @@ init_pull(FILE *fplog, const pull_params_t *pull_params, const t_inputrec *ir, for (int i = 0; i < pull_params->ngroup; ++i) { - pull->group.emplace_back(pull_params->group[i], atomSets->add({pull_params->group[i].ind, pull_params->group[i].ind+pull_params->group[i].nat})); + pull->group.emplace_back(pull_params->group[i], atomSets->add({pull_params->group[i].ind, pull_params->group[i].ind+pull_params->group[i].nat}), + pull_params->bSetPbcRefToPrevStepCOM); } if (cr != nullptr && DOMAINDECOMP(cr)) @@ -2142,6 +2150,7 @@ init_pull(FILE *fplog, const pull_params_t *pull_params, const t_inputrec *ir, pull->bAngle = FALSE; GMX_RELEASE_ASSERT(pull->group[0].params.nat == 0, "pull group 0 is an absolute reference group and should not contain atoms"); + pull->group[0].x_prev_step[XX] = NAN; pull->numCoordinatesWithExternalPotential = 0; @@ -2408,6 +2417,8 @@ init_pull(FILE *fplog, const pull_params_t *pull_params, const t_inputrec *ir, init_pull_group_index(fplog, cr, g, pgrp, bConstraint, pulldim_con, mtop, ir, lambda); + + pgrp->x_prev_step[XX] = NAN; } else { @@ -2432,7 +2443,7 @@ init_pull(FILE *fplog, const pull_params_t *pull_params, const t_inputrec *ir, } } const auto &referenceGroup = pull->group[coord.params.group[0]]; - pull->dyna.emplace_back(referenceGroup.params, referenceGroup.atomSet); + pull->dyna.emplace_back(referenceGroup.params, referenceGroup.atomSet, pull->params.bSetPbcRefToPrevStepCOM); } } diff --git a/src/gromacs/pulling/pull.h b/src/gromacs/pulling/pull.h index 89a9039bba..dc39c125de 100644 --- a/src/gromacs/pulling/pull.h +++ b/src/gromacs/pulling/pull.h @@ -68,6 +68,7 @@ struct t_filenm; struct t_inputrec; struct t_mdatoms; struct t_pbc; +class t_state; namespace gmx { @@ -291,6 +292,8 @@ void pull_calc_coms(const t_commrec *cr, /*! \brief Margin for checking pull group PBC distances compared to half the box size */ static constexpr real c_pullGroupPbcMargin = 0.9; +/*! \brief Threshold (as a factor of half the box size) for accepting pull groups without explicitly set refatom */ +static constexpr real c_pullGroupSmallGroupThreshold = 0.5; /*! \brief Checks whether all groups that use a reference atom are within PBC restrictions * @@ -368,4 +371,45 @@ gmx_bool pull_have_constraint(const struct pull_t *pull); real max_pull_distance2(const pull_coord_work_t *pcrd, const t_pbc *pbc); +/*! \brief Copies the COM from the previous step of all pull groups to the checkpoint state container + * + * \param[in] pull The COM pull force calculation data structure + * \param[in] state The global state container + */ +void setStatePrevStepPullCom(const struct pull_t *pull, t_state *state); + +/*! \brief Copies the pull group COM of the previous step from the checkpoint state to the pull state + * + * \param[in] pull The COM pull force calculation data structure + * \param[in] state The global state container + */ +void setPrevStepPullComFromState(struct pull_t *pull, const t_state *state); + +/*! \brief Sets the previous step COM to the current COM + * + * \param[in] pull The COM pull force calculation data structure + */ +void updatePrevStepCom(struct pull_t *pull); + +/*! \brief Resizes the vector, in the state container, containing the COMs from the previous step + * + * \param[in] state The global state container + * \param[in] pull The COM pull force calculation data structure + */ +void allocStatePrevStepPullCom(t_state *state, pull_t *pull); + +/*! \brief Initializes the COM of the previous step (set to initial COM) + * + * \param[in] cr Struct for communication info. + * \param[in] pull The pull data structure. + * \param[in] md All atoms. + * \param[in] pbc Information struct about periodicity. + * \param[in] x The local positions. + */ +void initPullComFromPrevStep(const t_commrec *cr, + pull_t *pull, + const t_mdatoms *md, + t_pbc *pbc, + const rvec x[]); + #endif diff --git a/src/gromacs/pulling/pull_internal.h b/src/gromacs/pulling/pull_internal.h index 9893601a4d..c56bb7c204 100644 --- a/src/gromacs/pulling/pull_internal.h +++ b/src/gromacs/pulling/pull_internal.h @@ -70,7 +70,7 @@ static const int c_pullMaxNumLocalAtomsSingleThreaded = 1; #endif enum { - epgrppbcNONE, epgrppbcREFAT, epgrppbcCOS + epgrppbcNONE, epgrppbcREFAT, epgrppbcCOS, epgrppbcPREVSTEPCOM }; /*! \internal @@ -80,11 +80,13 @@ struct pull_group_work_t { /*! \brief Constructor * - * \param[in] params The group parameters set by the user - * \param[in] atomSet The global to local atom set manager + * \param[in] params The group parameters set by the user + * \param[in] atomSet The global to local atom set manager + * \param[in] setPbcRefToPrevStepCOM Does this pull group use the COM from the previous step as reference position? */ pull_group_work_t(const t_pull_group ¶ms, - gmx::LocalAtomSet atomSet); + gmx::LocalAtomSet atomSet, + bool setPbcRefToPrevStepCOM); /* Data only modified at initialization */ const t_pull_group params; /**< The pull group parameters */ @@ -100,13 +102,14 @@ struct pull_group_work_t When no pbc refence atom is used, this pointer shall be null. */ /* Data, potentially, changed at every pull call */ - real mwscale; /**< mass*weight scaling factor 1/sum w m */ - real wscale; /**< scaling factor for the weights: sum w m/sum w w m */ - real invtm; /**< inverse total mass of the group: 1/wscale sum w m */ - std::vector < gmx::BasicVector < double>> mdw; /**< mass*gradient(weight) for atoms */ - std::vector dv; /**< distance to the other group(s) along vec */ - dvec x; /**< COM before update */ - dvec xp; /**< COM after update before constraining */ + real mwscale; /**< mass*weight scaling factor 1/sum w m */ + real wscale; /**< scaling factor for the weights: sum w m/sum w w m */ + real invtm; /**< inverse total mass of the group: 1/wscale sum w m */ + std::vector < gmx::BasicVector < double>> mdw; /**< mass*gradient(weight) for atoms */ + std::vector dv; /**< distance to the other group(s) along vec */ + dvec x; /**< COM before update */ + dvec xp; /**< COM after update before constraining */ + dvec x_prev_step; /**< center of mass of the previous step */ }; /* Struct describing the instantaneous spatial layout of a pull coordinate */ diff --git a/src/gromacs/pulling/pullutil.cpp b/src/gromacs/pulling/pullutil.cpp index b771a60b03..156458ddcf 100644 --- a/src/gromacs/pulling/pullutil.cpp +++ b/src/gromacs/pulling/pullutil.cpp @@ -47,8 +47,10 @@ #include "gromacs/math/utilities.h" #include "gromacs/math/vec.h" #include "gromacs/mdtypes/commrec.h" +#include "gromacs/mdtypes/inputrec.h" #include "gromacs/mdtypes/md_enums.h" #include "gromacs/mdtypes/mdatom.h" +#include "gromacs/mdtypes/state.h" #include "gromacs/pbcutil/pbc.h" #include "gromacs/pulling/pull.h" #include "gromacs/utility/fatalerror.h" @@ -163,7 +165,7 @@ static void pull_set_pbcatoms(const t_commrec *cr, struct pull_t *pull, for (size_t g = 0; g < pull->group.size(); g++) { const pull_group_work_t &group = pull->group[g]; - if (group.needToCalcCom && group.epgrppbc == epgrppbcREFAT) + if (group.needToCalcCom && (group.epgrppbc == epgrppbcREFAT || group.epgrppbc == epgrppbcPREVSTEPCOM)) { setPbcAtomCoords(pull->group[g], x, x_pbc[g]); numPbcAtoms++; @@ -582,12 +584,18 @@ void pull_calc_coms(const t_commrec *cr, { if (pgrp->epgrppbc != epgrppbcCOS) { - rvec x_pbc = { 0, 0, 0 }; + rvec x_pbc = { 0, 0, 0 }; - if (pgrp->epgrppbc == epgrppbcREFAT) + switch (pgrp->epgrppbc) { - /* Set the pbc atom */ - copy_rvec(comm->pbcAtomBuffer[g], x_pbc); + case epgrppbcREFAT: + /* Set the pbc atom */ + copy_rvec(comm->pbcAtomBuffer[g], x_pbc); + break; + case epgrppbcPREVSTEPCOM: + /* Set the pbc reference to the COM of the group of the last step */ + copy_dvec_to_rvec(pgrp->x_prev_step, comm->pbcAtomBuffer[g]); + copy_dvec_to_rvec(pgrp->x_prev_step, x_pbc); } /* The final sums should end up in comSums[0] */ @@ -742,7 +750,7 @@ void pull_calc_coms(const t_commrec *cr, { pgrp->xp[m] = dvecBuffer[1][m]*pgrp->mwscale; } - if (pgrp->epgrppbc == epgrppbcREFAT) + if (pgrp->epgrppbc == epgrppbcREFAT || pgrp->epgrppbc == epgrppbcPREVSTEPCOM) { pgrp->x[m] += comm->pbcAtomBuffer[g][m]; if (xp) @@ -925,7 +933,7 @@ int pullCheckPbcWithinGroups(const pull_t &pull, for (size_t g = 0; g < pull.group.size(); g++) { const pull_group_work_t &group = pull.group[g]; - if (group.epgrppbc == epgrppbcREFAT && + if ((group.epgrppbc == epgrppbcREFAT || group.epgrppbc == epgrppbcPREVSTEPCOM) && !pullGroupObeysPbcRestrictions(group, dimUsed[g], x, pbc, pull.comm.pbcAtomBuffer[g], pbcMargin)) { return g; @@ -949,7 +957,7 @@ bool pullCheckPbcWithinGroup(const pull_t &pull, /* Check PBC if the group uses a PBC reference atom treatment. */ const pull_group_work_t &group = pull.group[groupNr]; - if (group.epgrppbc != epgrppbcREFAT) + if (group.epgrppbc != epgrppbcREFAT && group.epgrppbc != epgrppbcPREVSTEPCOM) { return true; } @@ -977,3 +985,188 @@ bool pullCheckPbcWithinGroup(const pull_t &pull, return (pullGroupObeysPbcRestrictions(group, dimUsed, as_rvec_array(x.data()), pbc, pull.comm.pbcAtomBuffer[groupNr], pbcMargin)); } + +void setStatePrevStepPullCom(const struct pull_t *pull, t_state *state) +{ + for (size_t i = 0; i < state->com_prev_step.size()/DIM; i++) + { + for (int j = 0; j < DIM; j++) + { + state->com_prev_step[i*DIM+j] = pull->group[i].x_prev_step[j]; + } + } +} + +void setPrevStepPullComFromState(struct pull_t *pull, const t_state *state) +{ + for (size_t i = 0; i < state->com_prev_step.size()/DIM; i++) + { + for (int j = 0; j < DIM; j++) + { + pull->group[i].x_prev_step[j] = state->com_prev_step[i*DIM+j]; + } + } +} + +void updatePrevStepCom(struct pull_t *pull) +{ + for (size_t g = 0; g < pull->group.size(); g++) + { + if (pull->group[g].needToCalcCom) + { + for (int j = 0; j < DIM; j++) + { + pull->group[g].x_prev_step[j] = pull->group[g].x[j]; + } + } + } +} + +void allocStatePrevStepPullCom(t_state *state, pull_t *pull) +{ + if (!pull) + { + state->com_prev_step.clear(); + return; + } + size_t ngroup = pull->group.size(); + if (state->com_prev_step.size()/DIM != ngroup) + { + state->com_prev_step.resize(ngroup * DIM, NAN); + } +} + +void initPullComFromPrevStep(const t_commrec *cr, + pull_t *pull, + const t_mdatoms *md, + t_pbc *pbc, + const rvec x[]) +{ + pull_comm_t *comm = &pull->comm; + size_t ngroup = pull->group.size(); + + comm->pbcAtomBuffer.resize(ngroup); + comm->comBuffer.resize(ngroup*DIM); + + for (size_t g = 0; g < ngroup; g++) + { + pull_group_work_t *pgrp; + + pgrp = &pull->group[g]; + + if (pgrp->needToCalcCom && pgrp->epgrppbc == epgrppbcPREVSTEPCOM) + { + GMX_ASSERT(pgrp->params.nat > 1, "Groups with no atoms, or only one atom, should not " + "use the COM from the previous step as reference."); + + rvec x_pbc = { 0, 0, 0 }; + pull_set_pbcatoms(cr, pull, x, comm->pbcAtomBuffer); + copy_rvec(comm->pbcAtomBuffer[g], x_pbc); + + if (debug) + { + fprintf(debug, "Initialising prev step COM of pull group %zu. x_pbc =", g); + for (int m = 0; m < DIM; m++) + { + fprintf(debug, " %f", x_pbc[m]); + } + fprintf(debug, "\n"); + } + + /* The following is to a large extent similar to pull_calc_coms() */ + + /* The final sums should end up in sum_com[0] */ + ComSums &comSumsTotal = pull->comSums[0]; + + if (pgrp->atomSet.numAtomsLocal() <= c_pullMaxNumLocalAtomsSingleThreaded) + { + sum_com_part(pgrp, 0, pgrp->atomSet.numAtomsLocal(), + x, nullptr, md->massT, + pbc, x_pbc, + &comSumsTotal); + } + else + { +#pragma omp parallel for num_threads(pull->nthreads) schedule(static) + for (int t = 0; t < pull->nthreads; t++) + { + int ind_start = (pgrp->atomSet.numAtomsLocal()*(t + 0))/pull->nthreads; + int ind_end = (pgrp->atomSet.numAtomsLocal()*(t + 1))/pull->nthreads; + sum_com_part(pgrp, ind_start, ind_end, + x, nullptr, md->massT, + pbc, x_pbc, + &pull->comSums[t]); + } + + /* Reduce the thread contributions to sum_com[0] */ + for (int t = 1; t < pull->nthreads; t++) + { + comSumsTotal.sum_wm += pull->comSums[t].sum_wm; + comSumsTotal.sum_wwm += pull->comSums[t].sum_wwm; + dvec_inc(comSumsTotal.sum_wmx, pull->comSums[t].sum_wmx); + dvec_inc(comSumsTotal.sum_wmxp, pull->comSums[t].sum_wmxp); + } + } + + if (pgrp->localWeights.empty()) + { + comSumsTotal.sum_wwm = comSumsTotal.sum_wm; + } + + /* Copy local sums to a buffer for global summing */ + copy_dvec(comSumsTotal.sum_wmx, comm->comBuffer[g*3]); + copy_dvec(comSumsTotal.sum_wmxp, comm->comBuffer[g*3 + 1]); + comm->comBuffer[g*3 + 2][0] = comSumsTotal.sum_wm; + comm->comBuffer[g*3 + 2][1] = comSumsTotal.sum_wwm; + comm->comBuffer[g*3 + 2][2] = 0; + } + } + + pullAllReduce(cr, comm, ngroup*3*DIM, static_cast(comm->comBuffer[0])); + + for (size_t g = 0; g < ngroup; g++) + { + pull_group_work_t *pgrp; + + pgrp = &pull->group[g]; + if (pgrp->needToCalcCom) + { + if (pgrp->epgrppbc == epgrppbcPREVSTEPCOM) + { + double wmass, wwmass; + + /* Determine the inverse mass */ + wmass = comm->comBuffer[g*3+2][0]; + wwmass = comm->comBuffer[g*3+2][1]; + pgrp->mwscale = 1.0/wmass; + /* invtm==0 signals a frozen group, so then we should keep it zero */ + if (pgrp->invtm != 0) + { + pgrp->wscale = wmass/wwmass; + pgrp->invtm = wwmass/(wmass*wmass); + } + /* Divide by the total mass */ + for (int m = 0; m < DIM; m++) + { + pgrp->x[m] = comm->comBuffer[g*3 ][m]*pgrp->mwscale; + if (pgrp->epgrppbc == epgrppbcREFAT || pgrp->epgrppbc == epgrppbcPREVSTEPCOM) + { + pgrp->x[m] += comm->pbcAtomBuffer[g][m]; + } + } + if (debug) + { + fprintf(debug, "Pull group %zu wmass %f invtm %f\n", + g, 1.0/pgrp->mwscale, pgrp->invtm); + fprintf(debug, "Initialising prev step COM of pull group %zu to", g); + for (int m = 0; m < DIM; m++) + { + fprintf(debug, " %f", pgrp->x[m]); + } + fprintf(debug, "\n"); + } + copy_dvec(pgrp->x, pgrp->x_prev_step); + } + } + } +}