then defined as the angle between two planes: the plane spanned by the
the two first vectors and the plane spanned the two last vectors.
+ .. mdp-value:: transformation
+
+ Transforms other pull coordinates using a mathematical expression defined by :mdp:`pull-coord1-expression`.
+ Pull coordinates of lower indices can be used as variables to this pull coordinate.
+ Thus, pull transformation coordinates should have a higher pull coordinate index
+ than all pull coordinates they transform.
+
+.. mdp:: pull-coord1-expression
+
+ Mathematical expression to transform pull coordinates of lower indices to a new one.
+ The pull coordinates are referred to as variables in the equation so that
+ pull-coord1's value becomes 'x1', pull-coord2 value becomes 'x2' etc.
+ The mathematical expression are evaluated using muParser.
+ Only relevant if :mdp:`pull-coord1-geometry` is set to :mdp-value:`transformation`.
+
+.. mdp:: pull-coord1-dx
+
+ (1e-9)
+ Size of finite difference to use in numerical derivation of the pull coordinate
+ with respect to other pull coordinates.
+ The current implementation uses a simple first order finite difference method to perform derivation so that
+ f'(x) = (f(x+dx)-f(x))/dx
+ Only relevant if :mdp:`pull-coord1-geometry` is set to :mdp-value:`transformation`.
.. mdp:: pull-coord1-groups
tpxv_VSite1, /**< Added 1 type virtual site */
tpxv_MTS, /**< Added multiple time stepping */
tpxv_RemovedConstantAcceleration, /**< Removed support for constant acceleration NEMD. */
+ tpxv_TransformationPullCoord, /**< Support for transformation pull coordinates */
tpxv_Count /**< the total number of tpxv versions */
};
{
if (pcrd->eType == PullingAlgorithm::External)
{
- std::string buf;
- if (serializer->reading())
- {
- serializer->doString(&buf);
- pcrd->externalPotentialProvider = gmx_strdup(buf.c_str());
- }
- else
- {
- buf = pcrd->externalPotentialProvider;
- serializer->doString(&buf);
- }
+ serializer->doString(&pcrd->externalPotentialProvider);
}
else
{
pcrd->ngroup = 0;
}
serializer->doIvec(&pcrd->dim.as_vec());
+ if (file_version >= tpxv_TransformationPullCoord)
+ {
+ serializer->doString(&pcrd->expression);
+ }
+ else
+ {
+ if (serializer->reading())
+ {
+ pcrd->expression.clear();
+ }
+ }
}
else
{
bWarned = FALSE;
for (i = 0; i < ir->pull->ncoord && !bWarned; i++)
{
- if (ir->pull->coord[i].group[0] == 0 || ir->pull->coord[i].group[1] == 0)
+ if (ir->pull->coord[i].eGeom != PullGroupGeometry::Transformation
+ && (ir->pull->coord[i].group[0] == 0 || ir->pull->coord[i].group[1] == 0))
{
const auto absRef = haveAbsoluteReference(*ir);
const auto havePosres = havePositionRestraints(*sys);
#include <cassert>
#include <cstdlib>
#include <cstring>
+#include <regex>
#include "gromacs/domdec/localatomsetmanager.h"
#include "gromacs/fileio/readinp.h"
#include "gromacs/utility/arrayref.h"
#include "gromacs/utility/basedefinitions.h"
#include "gromacs/utility/cstringutil.h"
+#include "gromacs/utility/exceptions.h"
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/futil.h"
#include "gromacs/utility/smalloc.h"
+#include "gromacs/utility/stringutil.h"
static void string2dvec(const char buf[], dvec nums)
static void process_pull_dim(char* dim_buf, ivec dim, const t_pull_coord* pcrd)
{
- int ndim, d, nchar;
- char *ptr, pulldim1[STRLEN];
-
- ptr = dim_buf;
- ndim = 0;
- for (d = 0; d < DIM; d++)
+ char* ptr = dim_buf;
+ int ndim = 0;
+ for (int d = 0; d < DIM; d++)
{
+ int nchar = 0;
+ char pulldim1[STRLEN];
if (sscanf(ptr, "%s%n", pulldim1, &nchar) != 1)
{
gmx_fatal(FARGS, "Less than 3 pull dimensions given in pull_dim: '%s'", dim_buf);
}
}
-static void init_pull_coord(t_pull_coord* pcrd,
- int coord_index_for_output,
- char* dim_buf,
- const char* origin_buf,
- const char* vec_buf,
- warninp_t wi)
+static void initTransformationPullCoord(t_pull_coord* pcrd, const pull_params_t& pull)
+{
+ const int coord_index_for_output = pull.coord.size() + 1;
+ if (pcrd->eType == PullingAlgorithm::Constraint)
+ {
+ GMX_THROW(gmx::InvalidInputError(gmx::formatString(
+ "pull-coord%d can not have type 'constraint' and geometry 'transformation'",
+ coord_index_for_output)));
+ }
+
+ /*Validate the mathematical expression to epullgTRANSFORMATION*/
+ if (pcrd->expression.empty())
+ {
+ GMX_THROW(gmx::InvalidInputError(
+ gmx::formatString("pull-coord%d-expression not set for pull coordinate of geometry "
+ "'transformation'",
+ coord_index_for_output)));
+ }
+ if (pcrd->dx == 0)
+ {
+ GMX_THROW(gmx::InvalidInputError(gmx::formatString(
+ "pull-coord%d-dx cannot be set to zero for pull coordinate of geometry "
+ "'transformation'",
+ coord_index_for_output)));
+ }
+ /* make sure that the kappa of all previous pull coords is 0*/
+ int previousCoordOutputIndex = 0;
+ for (const auto& previousPcrd : pull.coord)
+ {
+ previousCoordOutputIndex++;
+ // See if the previous variable is used by the transformatino coord
+ // Note that a simple std::string::find won't work since we don't want x1 to match x11 etc.
+ std::string previousPcrdName = gmx::formatString("x%d(\\D|$)", previousCoordOutputIndex);
+ std::regex rx(previousPcrdName);
+
+ std::ptrdiff_t number_of_matches = std::distance(
+ std::sregex_iterator(pcrd->expression.begin(), pcrd->expression.end(), rx),
+ std::sregex_iterator());
+
+ if (number_of_matches == 0)
+ {
+ // This previous coordinate is not used in this transformation, do not check it
+ continue;
+ }
+
+ if (previousPcrd.eType == PullingAlgorithm::Constraint)
+ {
+ GMX_THROW(gmx::InvalidInputError(
+ gmx::formatString("pull-coord%d can not use pull-coord%d in the "
+ "transformation since this is a constraint",
+ coord_index_for_output,
+ previousCoordOutputIndex)));
+ }
+ }
+}
+
+static void init_pull_coord(t_pull_coord* pcrd,
+ char* dim_buf,
+ const char* origin_buf,
+ const char* vec_buf,
+ const pull_params_t& pull,
+ warninp_t wi)
{
+ const int coord_index_for_output = pull.coord.size() + 1;
+
int m;
dvec origin, vec;
char buf[STRLEN];
warning(wi, buf);
}
}
+ if (pcrd->eGeom == PullGroupGeometry::Transformation)
+ {
+ initTransformationPullCoord(pcrd, pull);
+ }
+
for (m = 0; m < DIM; m++)
{
pcrd->origin[m] = origin[m];
{
int nscan, idum;
char buf[STRLEN];
- char provider[STRLEN], groups[STRLEN], dim_buf[STRLEN];
+ char provider[STRLEN], groups[STRLEN], dim_buf[STRLEN], expression[STRLEN];
char wbuf[STRLEN], origin_buf[STRLEN], vec_buf[STRLEN];
/* read pull parameters */
sprintf(buf, "pull-coord%d-potential-provider", coordNum);
setStringEntry(inp, buf, provider, "");
pullCoord.externalPotentialProvider = provider;
+ sprintf(buf, "pull-coord%d-expression", coordNum);
+ setStringEntry(inp, buf, expression, "");
+ pullCoord.expression = expression;
+ sprintf(buf, "pull-coord%d-dx", coordNum);
+ pullCoord.dx = get_ereal(inp, buf, 1e-9, wi);
sprintf(buf, "pull-coord%d-geometry", coordNum);
pullCoord.eGeom = getEnum<PullGroupGeometry>(inp, buf, wi);
sprintf(buf, "pull-coord%d-groups", coordNum);
case PullGroupGeometry::Dihedral: pullCoord.ngroup = 6; break;
case PullGroupGeometry::DirectionRelative:
case PullGroupGeometry::Angle: pullCoord.ngroup = 4; break;
+ case PullGroupGeometry::Transformation: pullCoord.ngroup = 0; break;
default: pullCoord.ngroup = 2; break;
}
&pullCoord.group[4],
&pullCoord.group[5],
&idum);
+ if (nscan < 0)
+ {
+ // If the groups are not defined we can get a negative value here.
+ // It makes more sense to change it to 0
+ nscan = 0;
+ }
if (nscan != pullCoord.ngroup)
{
- auto message =
- gmx::formatString("%s should contain %d pull group indices with geometry %s",
- buf,
- pullCoord.ngroup,
- enumValueToString(pullCoord.eGeom));
+ auto message = gmx::formatString(
+ "%s should contain %d pull group indices with geometry %s."
+ " Found %d groups.",
+ buf,
+ pullCoord.ngroup,
+ enumValueToString(pullCoord.eGeom),
+ nscan);
set_warning_line(wi, nullptr, -1);
warning_error(wi, message);
}
pullCoord.kB = get_ereal(inp, buf, pullCoord.k, wi);
/* Initialize the pull coordinate */
- init_pull_coord(&pullCoord, coordNum, dim_buf, origin_buf, vec_buf, wi);
+ init_pull_coord(&pullCoord, dim_buf, origin_buf, vec_buf, *pull, wi);
pullCoord.coordIndex = coordNum - 1;
pull->coord.emplace_back(pullCoord);
void checkPullCoords(gmx::ArrayRef<const t_pull_group> pullGroups, gmx::ArrayRef<const t_pull_coord> pullCoords)
{
- for (int c = 0; c < int(pullCoords.size()); c++)
+ for (int c = 0; c < pullCoords.ssize(); ++c)
{
const t_pull_coord& pcrd = pullCoords[c];
GMX_RELEASE_ASSERT(pcrd.coordIndex == c, "coordIndex should match the index in the vector");
+ if (pcrd.eGeom == PullGroupGeometry::Transformation)
+ {
+ GMX_RELEASE_ASSERT(pcrd.ngroup == 0,
+ "Transformation coordinates don't use groups and "
+ "should have 'ngroup' set to zero");
+ continue;
+ }
if (pcrd.group[0] < 0 || pcrd.group[0] >= int(pullGroups.size()) || pcrd.group[1] < 0
|| pcrd.group[1] >= int(pullGroups.size()))
{
"Weights are not supported for the reference group with cylinder pulling");
}
}
+ c++;
}
}
{
pull_t* pull_work;
t_pbc pbc;
- int c;
- double t_start;
pull_params_t* pull = ir->pull.get();
gmx::LocalAtomSetManager atomSets;
set_pbc(&pbc, ir->pbcType, box);
- t_start = ir->init_t + ir->init_step * ir->delta_t;
+ double t_start = ir->init_t + ir->init_step * ir->delta_t;
if (pull->bSetPbcRefToPrevStepCOM)
{
}
fprintf(stderr, "Pull group natoms pbc atom distance at start reference at t=0\n");
- for (c = 0; c < pull->ncoord; c++)
+ for (int c = 0; c < pull->ncoord; c++)
{
- t_pull_coord* pcrd;
- t_pull_group *pgrp0, *pgrp1;
- double value;
- real init = 0;
+ real init = 0;
- pcrd = &pull->coord[c];
+ t_pull_coord* pcrd = &pull->coord[c];
- pgrp0 = &pull->group[pcrd->group[0]];
- pgrp1 = &pull->group[pcrd->group[1]];
+ t_pull_group* pgrp0 = &pull->group[pcrd->group[0]];
+ t_pull_group* pgrp1 = &pull->group[pcrd->group[1]];
fprintf(stderr, "%8d %8zu %8d\n", pcrd->group[0], pgrp0->ind.size(), pgrp0->pbcatom + 1);
fprintf(stderr, "%8d %8zu %8d ", pcrd->group[1], pgrp1->ind.size(), pgrp1->pbcatom + 1);
pcrd->init = 0;
}
- value = get_pull_coord_value(pull_work, c, pbc);
+ double value = get_pull_coord_value(pull_work, c, pbc);
value *= pull_conversion_factor_internal2userinput(*pcrd);
fprintf(stderr, " %10.3f %s", value, pull_coordinate_units(*pcrd));
static void pr_pull_coord(FILE* fp, int indent, int c, const t_pull_coord* pcrd)
{
- int g;
-
pr_indent(fp, indent);
fprintf(fp, "pull-coord %d:\n", c);
PS("type", enumValueToString(pcrd->eType));
PS("potential-provider", pcrd->externalPotentialProvider.c_str());
}
PS("geometry", enumValueToString(pcrd->eGeom));
- for (g = 0; g < pcrd->ngroup; g++)
+ for (int g = 0; g < pcrd->ngroup; g++)
{
- char buf[STRLEN];
-
- sprintf(buf, "group[%d]", g);
- PI(buf, pcrd->group[g]);
+ std::string buffer = gmx::formatString("group[%d]", g);
+ PI(buffer.c_str(), pcrd->group[g]);
}
pr_ivec(fp, indent, "dim", pcrd->dim, DIM, TRUE);
pr_rvec(fp, indent, "origin", pcrd->origin, DIM, TRUE);
const char* enumValueToString(PullGroupGeometry enumValue)
{
static constexpr gmx::EnumerationArray<PullGroupGeometry, const char*> pullGroupControlNames = {
- "distance", "direction", "cylinder", "direction-periodic",
- "direction-relative", "angle", "dihedral", "angle-axis"
+ "distance", "direction", "cylinder", "direction-periodic", "direction-relative",
+ "angle", "dihedral", "angle-axis", "transformation"
};
return pullGroupControlNames[enumValue];
}
Angle,
Dihedral,
AngleAxis,
+ Transformation,
Count,
Default = Distance
};
std::string externalPotentialProvider;
//! The pull geometry
PullGroupGeometry eGeom = PullGroupGeometry::Distance;
+ //! Mathematical expression evaluated by the pull code for transformation coordinates.
+ std::string expression;
+ //! The finite difference to use in numerical derivation of mathematical expressions
+ double dx = 1e-9;
//! The number of groups, depends on eGeom
int ngroup = 0;
/*! \brief The pull groups:
#include <algorithm>
#include <memory>
#include <mutex>
+#include <string>
#include "gromacs/commandline/filenm.h"
#include "gromacs/domdec/domdec_struct.h"
#include "gromacs/utility/stringutil.h"
#include "pull_internal.h"
+#include "transformationcoordinate.h"
namespace gmx
{
* region instead of potential parallel regions within apply_forces_grp.
* But there could be overlap between pull groups and this would lead
* to data races.
+ *
+ * This may also lead to potential issues with force redistribution
+ * for transformation pull coordinates
*/
if (pcrd.params.eGeom == PullGroupGeometry::Cylinder)
}
apply_forces_grp(pullGroups[pcrd.params.group[1]], masses, f_tot, 1, f);
}
+ else if (pcrd.params.eGeom == PullGroupGeometry::Transformation)
+ {
+ return;
+ }
else
{
if (pcrd.params.eGeom == PullGroupGeometry::DirectionRelative)
case PullGroupGeometry::AngleAxis:
spatialData.value = gmx_angle_between_dvecs(spatialData.dr01, spatialData.vec);
break;
+ case PullGroupGeometry::Transformation:
+ // Note that we would only need to pass the part of coord up to coord_ind
+ spatialData.value = gmx::getTransformationPullCoordinateValue(
+ pcrd, ArrayRef<const pull_coord_work_t>(pull.coord).subArray(0, pcrd->params.coordIndex));
+ break;
default: gmx_incons("Unsupported pull type in get_pull_coord_distance");
}
}
dsvmul(lambda * rm * pgrp0->invtm, vec, dr0);
dr_tot[c] += -lambda;
break;
+ case PullGroupGeometry::Transformation:
+ GMX_RELEASE_ASSERT(false, "transformation with constraints should never occur");
+ break;
default: gmx_incons("Invalid enumeration value for eGeom");
}
}
}
+/*! \brief Applies a force of any non-transformation pull coordinate
+ *
+ * \param[in] pull The pull struct, we can't pass only the coord because of cylinder pulling
+ * \param[in] coord_index The index of the coord to apply to force for
+ * \param[in] coord_force The force working on this coord
+ * \param[in] masses Atom masses
+ * \param[in] forceWithVirial Atom force and virial object
+ */
+static void apply_default_pull_coord_force(pull_t* pull,
+ const int coord_index,
+ const double coord_force,
+ gmx::ArrayRef<const real> masses,
+ gmx::ForceWithVirial* forceWithVirial)
+{
+ pull_coord_work_t* pcrd = &pull->coord[coord_index];
+ /* Set the force */
+ pcrd->scalarForce = coord_force;
+ /* Calculate the forces on the pull groups */
+ PullCoordVectorForces pullCoordForces = calculateVectorForces(*pcrd);
+
+ /* Add the forces for this coordinate to the total virial and force */
+ if (forceWithVirial->computeVirial_ && pull->comm.isMasterRank)
+ {
+ matrix virial = { { 0 } };
+ add_virial_coord(virial, *pcrd, pullCoordForces);
+ forceWithVirial->addVirialContribution(virial);
+ }
+ apply_forces_coord(pull->coord[coord_index],
+ pull->group,
+ pullCoordForces,
+ masses,
+ as_rvec_array(forceWithVirial->force_.data()));
+}
+
/* Pull takes care of adding the forces of the external potential.
* The external potential module has to make sure that the corresponding
* potential energy is added either to the pull term or to a term
GMX_ASSERT(pcrd->bExternalPotentialProviderHasBeenRegistered,
"apply_external_pull_coord_force called for an unregistered pull coordinate");
- /* Set the force */
- pcrd->scalarForce = coord_force;
-
- /* Calculate the forces on the pull groups */
- PullCoordVectorForces pullCoordForces = calculateVectorForces(*pcrd);
-
- /* Add the forces for this coordinate to the total virial and force */
- if (forceWithVirial->computeVirial_ && pull->comm.isMasterRank)
+ if (pcrd->params.eGeom == PullGroupGeometry::Transformation)
{
- matrix virial = { { 0 } };
- add_virial_coord(virial, *pcrd, pullCoordForces);
- forceWithVirial->addVirialContribution(virial);
+ gmx::applyTransformationPullCoordForce(
+ pcrd,
+ gmx::ArrayRef<pull_coord_work_t>(pull->coord).subArray(0, pcrd->params.coordIndex),
+ coord_force);
+ }
+ else
+ {
+ apply_default_pull_coord_force(pull, coord_index, coord_force, masses, forceWithVirial);
}
-
- apply_forces_coord(
- *pcrd, pull->group, pullCoordForces, masses, as_rvec_array(forceWithVirial->force_.data()));
}
-
pull->numExternalPotentialsStillToBeAppliedThisStep--;
}
{
continue;
}
-
PullCoordVectorForces pullCoordForces = do_pull_pot_coord(
*pull, &pcrd, pbc, t, lambda, &V, computeVirial ? virial : nullptr, &dVdl);
- /* Distribute the force over the atoms in the pulled groups */
- apply_forces_coord(pcrd, pull->group, pullCoordForces, masses, f);
+ if (pcrd.params.eGeom == PullGroupGeometry::Transformation)
+ {
+ gmx::applyTransformationPullCoordForce(
+ &pcrd,
+ gmx::ArrayRef<pull_coord_work_t>(pull->coord).subArray(0, pcrd.params.coordIndex),
+ pcrd.scalarForce);
+ }
+ else
+ {
+ /* Distribute the force over the atoms in the pulled groups */
+ apply_forces_coord(pcrd, pull->group, pullCoordForces, masses, f);
+ }
}
if (MASTER(cr))
case PullGroupGeometry::Direction:
case PullGroupGeometry::DirectionPBC:
case PullGroupGeometry::Cylinder:
+ case PullGroupGeometry::Transformation:
case PullGroupGeometry::AngleAxis:
copy_rvec_to_dvec(pull_params->coord[c].vec, pcrd->spatialData.vec);
break;
#include "gromacs/mdtypes/pull_params.h"
#include "gromacs/utility/gmxmpi.h"
+#include "pullcoordexpressionparser.h"
+
/*! \brief Determines up to what local atom count a pull group gets processed single-threaded.
*
* We set this limit to 1 with debug to catch bugs.
class PullHistory;
enum class PbcType : int;
+class t_state;
+
enum
{
epgrppbcNONE,
value_ref(0),
spatialData(),
scalarForce(0),
- bExternalPotentialProviderHasBeenRegistered(false)
+ bExternalPotentialProviderHasBeenRegistered(false),
+ expressionParser(params.eGeom == PullGroupGeometry::Transformation ? params.expression : "",
+ params.coordIndex),
+ transformationVariables(params.eGeom == PullGroupGeometry::Transformation ? params.coordIndex : 0)
{
}
//! For external-potential coordinates only, for checking if a provider has been registered
bool bExternalPotentialProviderHasBeenRegistered;
+
+ //! The expression parser for a transformation coordinate
+ gmx::PullCoordExpressionParser expressionParser;
+ //! Variables from other pull coordinates for a transformation coordinate
+ std::vector<double> transformationVariables;
};
/* Struct for storing vectorial forces for a pull coordinate */
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2021, 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.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+#include "gmxpre.h"
+
+#include "pullcoordexpressionparser.h"
+
+#include "config.h"
+
+#include "gromacs/utility/arrayref.h"
+#include "gromacs/utility/exceptions.h"
+#include "gromacs/utility/gmxassert.h"
+#include "gromacs/utility/stringutil.h"
+
+#include "pull_internal.h"
+
+namespace gmx
+{
+
+PullCoordExpressionParser::PullCoordExpressionParser(const std::string& expression, const int numVariables) :
+ expression_(expression)
+{
+#if HAVE_MUPARSER
+ if (!expression.empty())
+ {
+ // Initialize the parser
+ parser_ = std::make_unique<mu::Parser>();
+ parser_->SetExpr(expression_);
+ variableValues_.resize(numVariables);
+ for (int n = 0; n < numVariables; n++)
+ {
+ variableValues_[n] = 0;
+ std::string name = "x" + std::to_string(n + 1);
+ parser_->DefineVar(name, &variableValues_[n]);
+ }
+ }
+#else
+ GMX_UNUSED_VALUE(numVariables);
+ GMX_RELEASE_ASSERT(expression.empty(),
+ "Can not use transformation pull coordinate without muparser");
+#endif
+}
+
+double PullCoordExpressionParser::evaluate(ArrayRef<const double> variables)
+{
+#if HAVE_MUPARSER
+ GMX_ASSERT(variables.size() == variableValues_.size(),
+ "The size of variables should match the size passed at the first call of this "
+ "method");
+ // Todo: consider if we can use variableValues_ directly without the extra variables buffer
+ std::copy(variables.begin(), variables.end(), variableValues_.begin());
+
+ return parser_->Eval();
+#else
+ GMX_UNUSED_VALUE(variables);
+
+ GMX_RELEASE_ASSERT(false, "evaluate() should not be called without muparser");
+
+ return 0;
+#endif
+}
+
+} // namespace gmx
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2021, 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.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ *
+ * \brief
+ * Contains classes and methods related to use of MuParser in pulling
+ *
+ * \author Oliver Fleetwood <oliver.fleetwood@gmail.com>
+ * \author Paul Bauer <paul.bauer.q@gmail.com>
+ * \author Joe Jordan <ejjordan@kth.se>
+ * \author Berk Hess <hess@kth.se>
+ *
+ */
+#ifndef GMX_PULL_PULLCOORDEXPRESSIONPARSER_H
+#define GMX_PULL_PULLCOORDEXPRESSIONPARSER_H
+
+#include "config.h"
+
+#include <memory>
+#include <string>
+#include <vector>
+
+#if HAVE_MUPARSER
+# include <muParser.h>
+#else
+namespace mu
+{
+//! Defines a dummy Parser type to reduce use of the preprocessor.
+using Parser = std::false_type;
+} // namespace mu
+#endif
+
+struct pull_coord_work_t;
+
+namespace gmx
+{
+template<typename>
+class ArrayRef;
+
+/*! \brief Class with a mathematical expression and parser.
+ * \internal
+ *
+ * The class handles parser instantiation from an mathematical expression, e.g. 'x1*x2',
+ * and evaluates the expression given the variables' numerical values.
+ *
+ * Note that for performance reasons you should not create a new PullCoordExpressionParser
+ * for every evaluation.
+ * Instead, instantiate one PullCoordExpressionParser per expression,
+ * then update the variables before the next evaluation.
+ *
+ */
+class PullCoordExpressionParser
+{
+public:
+ //! Constructor which takes a mathematical expression and the number of variables as arguments
+ PullCoordExpressionParser(const std::string& expression, int numVariables);
+
+ //! Evaluates the expression with the numerical values passed in \p variables.
+ double evaluate(ArrayRef<const double> variables);
+
+private:
+ /*! \brief The mathematical expression, e.g. 'x1*x2' */
+ std::string expression_;
+
+ /*! \brief A vector containing the numerical values of the variables before parser evaluation.
+ *
+ * muParser compiles the expression to bytecode, then binds to the memory address
+ * of these vector elements, making the evaluations fast and memory efficient.
+ * */
+ std::vector<double> variableValues_;
+
+ /*! \brief The parser_ which compiles and evaluates the mathematical expression */
+ std::unique_ptr<mu::Parser> parser_;
+};
+
+} // namespace gmx
+
+#endif // GMX_PULL_PULLCOORDEXPRESSIONPARSER_H
#include "gromacs/pulling/pull.h"
+#include "config.h"
+
#include <cmath>
#include <algorithm>
#include "gromacs/math/vec.h"
#include "gromacs/pbcutil/pbc.h"
#include "gromacs/pulling/pull_internal.h"
+#include "gromacs/pulling/transformationcoordinate.h"
#include "gromacs/utility/smalloc.h"
#include "testutils/refdata.h"
test(PbcType::XY, box);
}
+#if HAVE_MUPARSER
+
+/*
+ * Simple test case with one transformation coordinate and one pull coordinate.
+ *
+ */
+TEST_F(PullTest, TransformationCoordSimple)
+{
+ //-----_SETUP-------
+ pull_t pull;
+ // Create standard pull coordinates
+ const double dist = 10; // the pull coord value
+ t_pull_coord x1;
+ x1.coordIndex = 0;
+ pull.coord.emplace_back(x1);
+ pull.coord[0].spatialData.value = dist;
+
+ // Create transformation pull coordinates
+ t_pull_coord x2;
+ x2.eGeom = PullGroupGeometry::Transformation;
+ x2.expression = "-x1";
+ x2.coordIndex = 1;
+ pull.coord.emplace_back(x2);
+
+ //-----TESTS -------
+ // 1) check transformation pull coordinate values
+ pull.coord[1].spatialData.value = getTransformationPullCoordinateValue(
+ &pull.coord[1], constArrayRefFromArray(pull.coord.data(), 1));
+ EXPECT_REAL_EQ_TOL(pull.coord[1].spatialData.value, -dist, defaultRealTolerance());
+
+ // 2) check that a force on the transformation coord propagates to the original pull coordinate
+ const double force = 10; // The force on the transformation coordinate
+ applyTransformationPullCoordForce(
+ &pull.coord[1], gmx::ArrayRef<pull_coord_work_t>(pull.coord).subArray(0, 1), force);
+ EXPECT_REAL_EQ_TOL(force, pull.coord[1].scalarForce, defaultRealTolerance());
+ EXPECT_REAL_EQ_TOL(-force,
+ pull.coord[0].scalarForce,
+ // Since we do numerical differentiation we need to increase the tolerance
+ test::relativeToleranceAsFloatingPoint(force, 1e-5));
+}
+
+/*
+ * More complicated test case of transformation pull coordinates.
+ * To account for all the different ways transformation coordinates can depend on each other,
+ * we need a relatively complicated setup with multiple pull coordinates.
+ * We also test the code for multiple numerical values, which adds a for loop.
+ */
+TEST_F(PullTest, TransformationCoordAdvanced)
+{
+ pull_t pull;
+ // Create standard pull coordinates
+ t_pull_coord x1;
+ x1.eGeom = PullGroupGeometry::Distance;
+ x1.coordIndex = 0;
+ pull.coord.emplace_back(x1);
+ t_pull_coord x2;
+ x2.eGeom = PullGroupGeometry::Angle;
+ x2.coordIndex = 1;
+ pull.coord.emplace_back(x2);
+
+ // Create transformation pull coordinates
+ // x3, a pull coordinate that depends on another pull coordinate
+ t_pull_coord x3;
+ x3.eGeom = PullGroupGeometry::Transformation;
+ std::string expression1 = "x1^2";
+ x3.expression = expression1;
+ x3.coordIndex = 2;
+ pull.coord.emplace_back(x3);
+
+ // x4, the last transformation pull coordinate
+ t_pull_coord x4;
+ x4.eGeom = PullGroupGeometry::Transformation;
+ std::string expression2 = "x1 - 0.5*x2^3 + x3^2 + 3"; // note that x3^2 is equivalent to x1^4
+ x4.expression = expression2;
+ x4.coordIndex = 3;
+ x4.dx = 1e-4;
+ pull.coord.emplace_back(x4);
+
+ // below we set x1 and x2 to different values and make sure that
+ // 1) the transformation coordinates are correct, i.e. test getTransformationPullCoordinateValue
+ // 2) that the force is accurately distributed from the transformation coord to the normal
+ // pull coordinates, i.e. test applyTransformationPullCoordForce
+ for (double v1 = 0.5; v1 < 11; v1 += 1)
+ {
+ clear_pull_forces(&pull);
+
+ double v2 = -v1 * 10;
+ // transformation pull coord value
+ pull.coord[0].spatialData.value = v1;
+ pull.coord[1].spatialData.value = v2;
+ pull.coord[2].spatialData.value = getTransformationPullCoordinateValue(
+ &pull.coord[2], constArrayRefFromArray(pull.coord.data(), 2));
+ pull.coord[3].spatialData.value = getTransformationPullCoordinateValue(
+ &pull.coord[3], constArrayRefFromArray(pull.coord.data(), 3));
+
+ // 1) check transformation pull coordinate values
+ // Since we perform numerical differentiation and floating point operations
+ // we only expect the results below to be approximately equal
+ double expectedX3 = v1 * v1;
+ EXPECT_REAL_EQ_TOL(pull.coord[2].spatialData.value, expectedX3, defaultRealTolerance());
+
+ double expectedX4 = v1 - 0.5 * v2 * v2 * v2 + expectedX3 * expectedX3 + 3;
+ EXPECT_REAL_EQ_TOL(pull.coord[3].spatialData.value, expectedX4, defaultRealTolerance());
+
+ // 2.1) check derivatives and force on normal pull coordinates with a direct relationship
+ double transformationForcex3 = v2 / 9 + 1;
+ applyTransformationPullCoordForce(&pull.coord[2],
+ gmx::ArrayRef<pull_coord_work_t>(pull.coord).subArray(0, 2),
+ transformationForcex3);
+
+ double expectedFx1 = transformationForcex3 * 2 * v1;
+ // the theoretical error of first order numerical derivation is 0.5*f''(x)*h (not taking the numerical precision into account)
+ double tolX3 = 1e-4; // * x3.dx; // numerical error tolerance
+ EXPECT_REAL_EQ_TOL(expectedFx1,
+ pull.coord[0].scalarForce,
+ test::relativeToleranceAsFloatingPoint(expectedFx1, tolX3));
+
+ double expectedFx2 = 0;
+ EXPECT_REAL_EQ_TOL(expectedFx2, pull.coord[1].scalarForce, defaultRealTolerance());
+
+ double expectedFx3 = transformationForcex3;
+ EXPECT_REAL_EQ_TOL(expectedFx3, pull.coord[2].scalarForce, defaultRealTolerance());
+
+ double expectedFx4 = 0;
+ EXPECT_REAL_EQ_TOL(expectedFx4, pull.coord[3].scalarForce, defaultRealTolerance());
+
+ // 2.2) check derivatives and force on normal pull coordinates
+ // also taking inner derivatives accounts (remembering that x3 = x1^2)
+ // Only x4 has non-zero scalar force here
+ double transformationForcex4 = v1 + 4.5;
+ double tolX4 = 1e-2;
+ applyTransformationPullCoordForce(&pull.coord[3],
+ gmx::ArrayRef<pull_coord_work_t>(pull.coord).subArray(0, 3),
+ transformationForcex4);
+
+ expectedFx1 += transformationForcex4
+ * (1 + 4 * v1 * v1 * v1); // Note that we expect the inner derivative to work here
+ EXPECT_REAL_EQ_TOL(expectedFx1,
+ pull.coord[0].scalarForce,
+ test::relativeToleranceAsFloatingPoint(expectedFx1, tolX4));
+
+ expectedFx2 += -1.5 * v2 * v2 * transformationForcex4;
+ EXPECT_REAL_EQ_TOL(expectedFx2,
+ pull.coord[1].scalarForce,
+ test::relativeToleranceAsFloatingPoint(expectedFx2, tolX4));
+
+ expectedFx3 += 2 * expectedX3 * transformationForcex4;
+ EXPECT_REAL_EQ_TOL(expectedFx3,
+ pull.coord[2].scalarForce,
+ test::relativeToleranceAsFloatingPoint(expectedFx3, tolX4));
+
+ expectedFx4 += transformationForcex4;
+ EXPECT_REAL_EQ_TOL(expectedFx4, pull.coord[3].scalarForce, defaultRealTolerance());
+ }
+}
+
+/*
+ * Simple test case with one transformation coordinate set to a constant value
+ *
+ * The purpose of this test is just to make sure that the code works even in this case
+ *
+ */
+TEST_F(PullTest, TransformationCoordDummyExpression)
+{
+ //-----_SETUP-------
+ pull_t pull;
+ t_pull_coord x;
+ x.eGeom = PullGroupGeometry::Transformation;
+ x.coordIndex = 0;
+ x.expression = "10";
+ pull.coord.emplace_back(x);
+ //-----TESTS -------
+ // check transformation pull coordinate values
+ double value =
+ getTransformationPullCoordinateValue(&pull.coord[0], ArrayRef<const pull_coord_work_t>{});
+ EXPECT_REAL_EQ_TOL(value, 10, defaultRealTolerance());
+}
+#endif // HAVE_MUPARSER
+
} // namespace
} // namespace gmx
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2021, 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.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+#include "gmxpre.h"
+
+#include "transformationcoordinate.h"
+
+#include "gromacs/utility/arrayref.h"
+#include "gromacs/utility/exceptions.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/gmxassert.h"
+#include "gromacs/utility/stringutil.h"
+
+#include "pull_internal.h"
+#include "pullcoordexpressionparser.h"
+
+namespace gmx
+{
+
+namespace
+{
+
+//! Calculates the value a for transformation pull coordinate
+double getTransformationPullCoordinateValue(pull_coord_work_t* coord)
+{
+ const int transformationPullCoordinateIndex = coord->params.coordIndex;
+ GMX_ASSERT(ssize(coord->transformationVariables) == transformationPullCoordinateIndex,
+ "We need as many variables as the transformation pull coordinate index");
+ double result = 0;
+ try
+ {
+ result = coord->expressionParser.evaluate(coord->transformationVariables);
+ }
+ catch (mu::Parser::exception_type& e)
+ {
+ GMX_THROW(InternalError(
+ formatString("failed to evaluate expression for transformation pull-coord%d: %s\n",
+ transformationPullCoordinateIndex + 1,
+ e.GetMsg().c_str())));
+ }
+ catch (std::exception& e)
+ {
+ GMX_THROW(InternalError(
+ formatString("failed to evaluate expression for transformation pull-coord%d.\n"
+ "Last variable pull-coord-index: %d.\n"
+ "Message: %s\n",
+ transformationPullCoordinateIndex + 1,
+ transformationPullCoordinateIndex + 1,
+ e.what())));
+ }
+ return result;
+}
+
+} // namespace
+
+double getTransformationPullCoordinateValue(pull_coord_work_t* coord,
+ ArrayRef<const pull_coord_work_t> variableCoords)
+{
+ GMX_ASSERT(ssize(variableCoords) == coord->params.coordIndex,
+ "We need as many variables as the transformation pull coordinate index");
+ int coordIndex = 0;
+ for (const auto& variableCoord : variableCoords)
+ {
+ coord->transformationVariables[coordIndex++] = variableCoord.spatialData.value;
+ }
+
+ return getTransformationPullCoordinateValue(coord);
+}
+
+/*! \brief Calculates and returns the derivative of a transformation pull coordinate from a dependent coordinate
+ *
+ * Note #1: this requires that getTransformationPullCoordinateValue() has been called
+ * before with the current coordinates.
+ *
+ * Note #2: this method will not compute inner derivates. That is taken care of in the regular pull code
+ *
+ * \param[in] coord The (transformation) coordinate to compute the value for
+ * \param[in] variablePcrdIndex Pull coordinate index of a variable.
+ */
+static double computeDerivativeForTransformationPullCoord(pull_coord_work_t* coord, const int variablePcrdIndex)
+{
+ GMX_ASSERT(variablePcrdIndex >= 0 && variablePcrdIndex < coord->params.coordIndex,
+ "The variable index should be in range of the transformation coordinate");
+
+ // epsilon for numerical differentiation.
+ const double transformationPcrdValue = coord->spatialData.value;
+ // Perform numerical differentiation of 1st order
+ const double valueBackup = coord->transformationVariables[variablePcrdIndex];
+ double dx = coord->params.dx;
+ coord->transformationVariables[variablePcrdIndex] += dx;
+ double transformationPcrdValueEps = getTransformationPullCoordinateValue(coord);
+ double derivative = (transformationPcrdValueEps - transformationPcrdValue) / dx;
+ // reset pull coordinate value
+ coord->transformationVariables[variablePcrdIndex] = valueBackup;
+ return derivative;
+}
+
+
+/*!
+ * \brief Distributes the force from a transformation pull coordiante to the dependent pull
+ * coordinates by computing the inner derivatives
+ *
+ * \param pcrd The transformation pull coord
+ * \param variableCoords The dependent pull coordinates
+ * \param transformationCoordForce The force to distribute
+ */
+static void distributeTransformationPullCoordForce(pull_coord_work_t* pcrd,
+ gmx::ArrayRef<pull_coord_work_t> variableCoords,
+ const double transformationCoordForce)
+{
+ if (std::abs(transformationCoordForce) < 1e-9)
+ {
+ // the force is effectively 0. Don't proceed and distribute it recursively
+ return;
+ }
+ GMX_ASSERT(pcrd->params.eGeom == PullGroupGeometry::Transformation,
+ "We shouldn't end up here when not using a transformation pull coordinate.");
+ GMX_ASSERT(ssize(variableCoords) == pcrd->params.coordIndex,
+ "We should have as many variable coords as the coord index of the transformation "
+ "coordinate");
+
+ // pcrd->scalarForce += transformationCoordForce;
+ for (auto& variableCoord : variableCoords)
+ {
+ const double derivative =
+ computeDerivativeForTransformationPullCoord(pcrd, variableCoord.params.coordIndex);
+ const double variablePcrdForce = transformationCoordForce * derivative;
+ /* Since we loop over all pull coordinates with smaller index, there can be ones
+ * that are not referenced by the transformation coordinate. Avoid apply forces
+ * on those by skipping application of zero force.
+ */
+ if (variablePcrdForce != 0)
+ {
+ if (debug)
+ {
+ fprintf(debug,
+ "Distributing force %4.4f for transformation coordinate %d to coordinate "
+ "%d with "
+ "force "
+ "%4.4f\n",
+ transformationCoordForce,
+ pcrd->params.coordIndex,
+ variableCoord.params.coordIndex,
+ variablePcrdForce);
+ }
+ // Note that we add to the force here, in case multiple biases act on the same pull
+ // coord (although that is not recommended it should still work)
+ variableCoord.scalarForce += variablePcrdForce;
+ if (variableCoord.params.eGeom == PullGroupGeometry::Transformation)
+ {
+ /*
+ * We can have a transformation pull coordinate depend on another transformation pull coordinate
+ * which in turn leads to inner derivatives between pull coordinates.
+ * Here we redistribute the force via the inner product
+ *
+ * Note that this only works properly if the lower ranked transformation pull coordinate has it's scalarForce set to zero
+ */
+ distributeTransformationPullCoordForce(
+ &variableCoord,
+ variableCoords.subArray(0, variableCoord.params.coordIndex),
+ variablePcrdForce);
+ }
+ }
+ }
+}
+
+void applyTransformationPullCoordForce(pull_coord_work_t* pcrd,
+ gmx::ArrayRef<pull_coord_work_t> variableCoords,
+ const double transformationCoordForce)
+{
+ pcrd->scalarForce = transformationCoordForce;
+ // Note on why we need to call another method here:
+ // applyTransformationPullCoordForce is the method that should be called by the rest of the pull code.
+ // It's non-recursive and called exactly once for every transformation coordinate for every timestep.
+ // In it, we set the force on the transformation coordinate,
+ // then pass the force on to the other pull coordinates via the method distributeTransformationPullCoordForce.
+ // The latter method is recursive to account for inner derivatives.
+ // Note that we don't set the force on the top-level transformation coordinate in distributeTransformationPullCoordForce,
+ // we only add to the force, which is why it can be recursive.
+ distributeTransformationPullCoordForce(pcrd, variableCoords, transformationCoordForce);
+}
+
+
+} // namespace gmx
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2021, 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.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ *
+ * \brief
+ * Declares function for compute transformation coordinate values and forces
+ *
+ * \author Oliver Fleetwood <oliver.fleetwood@gmail.com>
+ * \author Paul Bauer <paul.bauer.q@gmail.com>
+ * \author Joe Jordan <ejjordan@kth.se>
+ * \author Berk Hess <hess@kth.se>
+ *
+ */
+#ifndef GMX_PULL_TRANSFORMATIONCOORDINATE_H
+#define GMX_PULL_TRANSFORMATIONCOORDINATE_H
+
+struct pull_t;
+struct pull_coord_work_t;
+
+namespace gmx
+{
+template<typename>
+class ArrayRef;
+
+/*! \brief Calculates pull->coord[coord_ind].spatialData.value for a transformation pull coordinate
+ *
+ * This requires the values of the pull coordinates of lower indices to be set
+ * \param[in] coord The (transformation) coordinate to compute the value for
+ * \param[in] variableCoords Pull coordinates used as variables, entries 0 to coord->coordIndex
+ * will be used \returns Transformation value for pull coordinate.
+ */
+double getTransformationPullCoordinateValue(pull_coord_work_t* coord,
+ ArrayRef<const pull_coord_work_t> variableCoords);
+
+/*! \brief Applies a force of a transformation pull coordinate and distributes it to pull coordinates of lower rank
+ *
+ * \param[in,out] pcrd The transformation pull coordinate to act on
+ * \param[in,out] variableCoords List of variable coords up to the coord index of \p pcrd
+ * \param[in] transformationCoordForce The force working on coord \p pcrd
+ */
+void applyTransformationPullCoordForce(pull_coord_work_t* pcrd,
+ gmx::ArrayRef<pull_coord_work_t> variableCoords,
+ double transformationCoordForce);
+
+} // namespace gmx
+
+#endif // GMX_PULL_TRANSFORMATIONCOORDINATE_H
# use for gmx::compat::optional. These are included as system
# headers so that no warnings are issued from them.
target_include_directories(${EXENAME} SYSTEM PRIVATE ${PROJECT_SOURCE_DIR}/src/external)
+ target_include_directories(${EXENAME} SYSTEM PRIVATE ${PROJECT_SOURCE_DIR}/src/external/muparser)
if(CYGWIN)
# Ensure GoogleTest headers can find POSIX things needed
target_compile_definitions(${EXENAME} PRIVATE _POSIX_C_SOURCE=200809L)
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff