: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
^^^^^^^^^^^^^^^^^^^^^^^
(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
* 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;
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
}
- 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 =
/*
* 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.
* \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.
*
pull_t *pull_work,
const matrix box,
int ePBC,
+ const tensor &compressibility,
const t_grpopts *inputrecGroupOptions,
warninp_t wi);
*
* 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.
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)
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff