:issue:`2554`
+grompp now checks that pull groups are not close to half the box size
+""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+Pull groups that use a reference atom for periodic boundary treatment
+should have all their atoms well within half the box size of this reference.
+When this is not the case, grompp will issue a warning.
+
+:issue:`2397`
+
Fixes for ``gmx`` tools
^^^^^^^^^^^^^^^^^^^^^^^
if (ir->bPull)
{
- pull = set_pull_init(ir, sys, as_rvec_array(state.x.data()), state.box, state.lambda[efptMASS], oenv);
+ pull = set_pull_init(ir, sys, as_rvec_array(state.x.data()), state.box, state.lambda[efptMASS], oenv, wi);
}
/* Modules that supply external potential for pull coordinates
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2016,2017, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2016,2017,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.
pull_t *set_pull_init(t_inputrec *ir, const gmx_mtop_t *mtop,
rvec *x, matrix box, real lambda,
- const gmx_output_env_t *oenv);
+ const gmx_output_env_t *oenv,
+ warninp_t wi);
/* Prints the initial pull group distances in x.
* If requested, adds the current distance to the initial reference location.
* Returns the pull_t pull work struct. This should be passed to finish_pull()
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2016,2017, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2016,2017,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.
pull_t *set_pull_init(t_inputrec *ir, const gmx_mtop_t *mtop,
rvec *x, matrix box, real lambda,
- const gmx_output_env_t *oenv)
+ const gmx_output_env_t *oenv,
+ warninp_t wi)
{
pull_params_t *pull;
pull_t *pull_work;
pull_calc_coms(nullptr, pull_work, md, &pbc, t_start, x, nullptr);
+ int groupThatFailsPbc = pullCheckPbcWithinGroups(*pull_work, x, pbc);
+ if (groupThatFailsPbc >= 0)
+ {
+ 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);
+ }
+
fprintf(stderr, "Pull group natoms pbc atom distance at start reference at t=0\n");
for (c = 0; c < pull->ncoord; c++)
{
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2016,2017, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2016,2017,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.
rvec x[],
rvec *xp);
+/*! \brief Margin for checking PBC distances compared to half the box size in pullCheckPbcWithinGroups() */
+static constexpr real c_pullGroupPbcMargin = 0.9;
+
+/*! \brief Checks whether all groups that use a reference atom are within PBC restrictions
+ *
+ * Groups that use a reference atom for determining PBC should have all their
+ * atoms within half the box size from the PBC atom. The box size is used
+ * per dimension for rectangular boxes, but can be a combination of
+ * dimensions for triclinic boxes, depending on which dimensions are
+ * involved in the pull coordinates a group is involved in.
+ *
+ * Should be called without MPI parallelization and after pull_calc_coms()
+ * has been called at least once.
+ *
+ * \param[in] pull The pull data structure
+ * \param[in] x The coordinates
+ * \param[in] pbc Information struct about periodicity
+ * \returns -1 when all groups obey PBC or the first group index that fails PBC
+ */
+int pullCheckPbcWithinGroups(const pull_t &pull,
+ const rvec *x,
+ const t_pbc &pbc);
/*! \brief Returns if we have pull coordinates with potential pulling.
*
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2016,2017, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2016,2017,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.
make_cyl_refgrps(cr, pull, md, pbc, t, x);
}
}
+
+using BoolVec = gmx::BasicVector<bool>;
+
+/* Returns whether the pull group obeys the PBC restrictions */
+static bool pullGroupObeysPbcRestrictions(const pull_group_work_t &group,
+ const BoolVec &dimUsed,
+ const rvec *x,
+ const t_pbc &pbc,
+ const rvec &x_pbc)
+{
+ /* Determine which dimensions are relevant for PBC */
+ BoolVec dimUsesPbc = { false, false, false };
+ bool pbcIsRectangular = true;
+ for (int d = 0; d < pbc.ndim_ePBC; d++)
+ {
+ if (dimUsed[d])
+ {
+ dimUsesPbc[d] = true;
+ /* All non-zero dimensions of vector v are involved in PBC */
+ for (int d2 = d + 1; d2 < pbc.ndim_ePBC; d2++)
+ {
+ assert(d2 < DIM);
+ if (pbc.box[d2][d] != 0)
+ {
+ dimUsesPbc[d2] = true;
+ pbcIsRectangular = false;
+ }
+ }
+ }
+ }
+
+ rvec marginPerDim;
+ real marginDistance2 = 0;
+ if (pbcIsRectangular)
+ {
+ /* Use margins for dimensions independently */
+ for (int d = 0; d < pbc.ndim_ePBC; d++)
+ {
+ marginPerDim[d] = c_pullGroupPbcMargin*pbc.hbox_diag[d];
+ }
+ }
+ else
+ {
+ /* Check the total distance along the relevant dimensions */
+ for (int d = 0; d < pbc.ndim_ePBC; d++)
+ {
+ if (dimUsesPbc[d])
+ {
+ marginDistance2 += c_pullGroupPbcMargin*gmx::square(0.5)*norm2(pbc.box[d]);
+ }
+ }
+ }
+
+ for (int i = 0; i < group.nat_loc; i++)
+ {
+ rvec dx;
+ pbc_dx(&pbc, x[group.ind_loc[i]], x_pbc, dx);
+
+ bool atomIsTooFar = false;
+ if (pbcIsRectangular)
+ {
+ for (int d = 0; d < pbc.ndim_ePBC; d++)
+ {
+ if (dimUsesPbc[d] && (dx[d] < -marginPerDim[d] ||
+ dx[d] > marginPerDim[d]))
+ {
+ atomIsTooFar = true;
+ }
+ }
+ }
+ else
+ {
+ real pbcDistance2 = 0;
+ for (int d = 0; d < pbc.ndim_ePBC; d++)
+ {
+ if (dimUsesPbc[d])
+ {
+ pbcDistance2 += gmx::square(dx[d]);
+ }
+ }
+ atomIsTooFar = (pbcDistance2 > marginDistance2);
+ }
+ if (atomIsTooFar)
+ {
+ return false;
+ }
+ }
+
+ return true;
+}
+
+int pullCheckPbcWithinGroups(const pull_t &pull,
+ const rvec *x,
+ const t_pbc &pbc)
+{
+ if (pbc.ePBC == epbcNONE)
+ {
+ return -1;
+ }
+
+ /* Determine what dimensions are used for each group by pull coordinates */
+ std::vector<BoolVec> dimUsed(pull.ngroup, { false, false, false });
+ for (int c = 0; c < pull.ncoord; c++)
+ {
+ const t_pull_coord &coordParams = pull.coord[c].params;
+ for (int groupIndex = 0; groupIndex < coordParams.ngroup; groupIndex++)
+ {
+ for (int d = 0; d < DIM; d++)
+ {
+ if (coordParams.dim[d] &&
+ !(coordParams.eGeom == epullgCYL && groupIndex == 0))
+ {
+ dimUsed[coordParams.group[groupIndex]][d] = true;
+ }
+ }
+ }
+ }
+
+ /* Check PBC for every group that uses a PBC reference atom treatment */
+ for (int g = 0; g < pull.ngroup; g++)
+ {
+ const pull_group_work_t &group = pull.group[g];
+ if (group.epgrppbc == epgrppbcREFAT &&
+ !pullGroupObeysPbcRestrictions(group, dimUsed[g], x, pbc, pull.comm.rbuf[g]))
+ {
+ return g;
+ }
+ }
+
+ return -1;
+}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff