From 7fc99ab9884f771cea80072bae4f2aa70602e5c0 Mon Sep 17 00:00:00 2001 From: Paul Bauer Date: Fri, 3 Sep 2021 11:16:31 +0000 Subject: [PATCH] Add support for transformation pull coordinates --- docs/user-guide/mdp-options.rst | 23 ++ src/gromacs/fileio/tpxio.cpp | 24 +- src/gromacs/gmxpreprocess/readir.cpp | 3 +- src/gromacs/gmxpreprocess/readpull.cpp | 146 +++++++++--- src/gromacs/mdtypes/inputrec.cpp | 10 +- src/gromacs/mdtypes/md_enums.cpp | 4 +- src/gromacs/mdtypes/md_enums.h | 1 + src/gromacs/mdtypes/pull_params.h | 4 + src/gromacs/pulling/pull.cpp | 91 ++++++-- src/gromacs/pulling/pull_internal.h | 14 +- .../pulling/pullcoordexpressionparser.cpp | 94 ++++++++ .../pulling/pullcoordexpressionparser.h | 110 +++++++++ src/gromacs/pulling/tests/pull.cpp | 182 +++++++++++++++ .../pulling/transformationcoordinate.cpp | 215 ++++++++++++++++++ .../pulling/transformationcoordinate.h | 79 +++++++ src/testutils/TestMacros.cmake | 1 + 16 files changed, 929 insertions(+), 72 deletions(-) create mode 100644 src/gromacs/pulling/pullcoordexpressionparser.cpp create mode 100644 src/gromacs/pulling/pullcoordexpressionparser.h create mode 100644 src/gromacs/pulling/transformationcoordinate.cpp create mode 100644 src/gromacs/pulling/transformationcoordinate.h diff --git a/docs/user-guide/mdp-options.rst b/docs/user-guide/mdp-options.rst index beb3f2972b..590b5425e1 100644 --- a/docs/user-guide/mdp-options.rst +++ b/docs/user-guide/mdp-options.rst @@ -1792,6 +1792,29 @@ pull-coord2-vec, pull-coord2-k, and so on. 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 diff --git a/src/gromacs/fileio/tpxio.cpp b/src/gromacs/fileio/tpxio.cpp index d70585a9e4..f5bf5874eb 100644 --- a/src/gromacs/fileio/tpxio.cpp +++ b/src/gromacs/fileio/tpxio.cpp @@ -137,6 +137,7 @@ enum tpxv 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 */ }; @@ -301,17 +302,7 @@ static void do_pull_coord(gmx::ISerializer* serializer, { 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 { @@ -350,6 +341,17 @@ static void do_pull_coord(gmx::ISerializer* serializer, 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 { diff --git a/src/gromacs/gmxpreprocess/readir.cpp b/src/gromacs/gmxpreprocess/readir.cpp index adba76ef7d..6b8c64ea3a 100644 --- a/src/gromacs/gmxpreprocess/readir.cpp +++ b/src/gromacs/gmxpreprocess/readir.cpp @@ -4551,7 +4551,8 @@ void triple_check(const char* mdparin, t_inputrec* ir, gmx_mtop_t* sys, warninp_ 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); diff --git a/src/gromacs/gmxpreprocess/readpull.cpp b/src/gromacs/gmxpreprocess/readpull.cpp index e14cd2a534..b8062ee54c 100644 --- a/src/gromacs/gmxpreprocess/readpull.cpp +++ b/src/gromacs/gmxpreprocess/readpull.cpp @@ -40,6 +40,7 @@ #include #include #include +#include #include "gromacs/domdec/localatomsetmanager.h" #include "gromacs/fileio/readinp.h" @@ -57,9 +58,11 @@ #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) @@ -88,13 +91,12 @@ static std::vector setupPullGroupWeights(const char* wbuf) 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); @@ -132,13 +134,71 @@ static void process_pull_dim(char* dim_buf, ivec dim, const t_pull_coord* pcrd) } } -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]; @@ -297,6 +357,11 @@ static void init_pull_coord(t_pull_coord* pcrd, warning(wi, buf); } } + if (pcrd->eGeom == PullGroupGeometry::Transformation) + { + initTransformationPullCoord(pcrd, pull); + } + for (m = 0; m < DIM; m++) { pcrd->origin[m] = origin[m]; @@ -308,7 +373,7 @@ std::vector read_pullparams(std::vector* inp, pull_param { 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 */ @@ -375,6 +440,11 @@ std::vector read_pullparams(std::vector* inp, pull_param 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(inp, buf, wi); sprintf(buf, "pull-coord%d-groups", coordNum); @@ -385,6 +455,7 @@ std::vector read_pullparams(std::vector* inp, pull_param 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; } @@ -397,13 +468,21 @@ std::vector read_pullparams(std::vector* inp, pull_param &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); } @@ -439,7 +518,7 @@ std::vector read_pullparams(std::vector* inp, pull_param 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); @@ -520,12 +599,19 @@ void process_pull_groups(gmx::ArrayRef pullGroups, void checkPullCoords(gmx::ArrayRef pullGroups, gmx::ArrayRef 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())) { @@ -551,6 +637,7 @@ void checkPullCoords(gmx::ArrayRef pullGroups, gmx::ArrayRef "Weights are not supported for the reference group with cylinder pulling"); } } + c++; } } @@ -563,8 +650,6 @@ pull_t* set_pull_init(t_inputrec* ir, { pull_t* pull_work; t_pbc pbc; - int c; - double t_start; pull_params_t* pull = ir->pull.get(); gmx::LocalAtomSetManager atomSets; @@ -579,7 +664,7 @@ pull_t* set_pull_init(t_inputrec* ir, 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) { @@ -641,17 +726,14 @@ pull_t* set_pull_init(t_inputrec* ir, } 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); @@ -661,7 +743,7 @@ pull_t* set_pull_init(t_inputrec* ir, 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)); diff --git a/src/gromacs/mdtypes/inputrec.cpp b/src/gromacs/mdtypes/inputrec.cpp index e54e8c3a7c..09d85d86b1 100644 --- a/src/gromacs/mdtypes/inputrec.cpp +++ b/src/gromacs/mdtypes/inputrec.cpp @@ -472,8 +472,6 @@ static void pr_pull_group(FILE* fp, int indent, int g, const t_pull_group* pgrp) 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)); @@ -482,12 +480,10 @@ static void pr_pull_coord(FILE* fp, int indent, int c, const t_pull_coord* pcrd) 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); diff --git a/src/gromacs/mdtypes/md_enums.cpp b/src/gromacs/mdtypes/md_enums.cpp index 7221e51ab2..38647e60b0 100644 --- a/src/gromacs/mdtypes/md_enums.cpp +++ b/src/gromacs/mdtypes/md_enums.cpp @@ -340,8 +340,8 @@ const char* enumValueToString(PullingAlgorithm enumValue) const char* enumValueToString(PullGroupGeometry enumValue) { static constexpr gmx::EnumerationArray 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]; } diff --git a/src/gromacs/mdtypes/md_enums.h b/src/gromacs/mdtypes/md_enums.h index 57f2a34b72..168c0dc5e9 100644 --- a/src/gromacs/mdtypes/md_enums.h +++ b/src/gromacs/mdtypes/md_enums.h @@ -685,6 +685,7 @@ enum class PullGroupGeometry : int Angle, Dihedral, AngleAxis, + Transformation, Count, Default = Distance }; diff --git a/src/gromacs/mdtypes/pull_params.h b/src/gromacs/mdtypes/pull_params.h index 8e623df1bc..095b2484ac 100644 --- a/src/gromacs/mdtypes/pull_params.h +++ b/src/gromacs/mdtypes/pull_params.h @@ -80,6 +80,10 @@ struct t_pull_coord 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: diff --git a/src/gromacs/pulling/pull.cpp b/src/gromacs/pulling/pull.cpp index 80e81b44b4..106861bb68 100644 --- a/src/gromacs/pulling/pull.cpp +++ b/src/gromacs/pulling/pull.cpp @@ -49,6 +49,7 @@ #include #include #include +#include #include "gromacs/commandline/filenm.h" #include "gromacs/domdec/domdec_struct.h" @@ -84,6 +85,7 @@ #include "gromacs/utility/stringutil.h" #include "pull_internal.h" +#include "transformationcoordinate.h" namespace gmx { @@ -339,6 +341,9 @@ static void apply_forces_coord(const pull_coord_work_t& pcrd, * 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) @@ -354,6 +359,10 @@ static void apply_forces_coord(const pull_coord_work_t& pcrd, } 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) @@ -746,6 +755,11 @@ static void get_pull_coord_distance(const pull_t& pull, pull_coord_work_t* pcrd, 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(pull.coord).subArray(0, pcrd->params.coordIndex)); + break; default: gmx_incons("Unsupported pull type in get_pull_coord_distance"); } } @@ -999,6 +1013,9 @@ static void do_constraint(struct pull_t* pull, 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"); } @@ -1520,6 +1537,40 @@ static void check_external_potential_registration(const struct pull_t* pull) } } +/*! \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 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 @@ -1547,24 +1598,18 @@ void apply_external_pull_coord_force(struct pull_t* pull, 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).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--; } @@ -1631,12 +1676,21 @@ real pull_potential(struct pull_t* pull, { 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).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)) @@ -2039,6 +2093,7 @@ struct pull_t* init_pull(FILE* fplog, 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; diff --git a/src/gromacs/pulling/pull_internal.h b/src/gromacs/pulling/pull_internal.h index 14f19778db..3d095e3f69 100644 --- a/src/gromacs/pulling/pull_internal.h +++ b/src/gromacs/pulling/pull_internal.h @@ -58,6 +58,8 @@ #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. @@ -73,6 +75,8 @@ static const int c_pullMaxNumLocalAtomsSingleThreaded = 1; class PullHistory; enum class PbcType : int; +class t_state; + enum { epgrppbcNONE, @@ -154,7 +158,10 @@ struct pull_coord_work_t 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) { } @@ -174,6 +181,11 @@ struct pull_coord_work_t //! 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 transformationVariables; }; /* Struct for storing vectorial forces for a pull coordinate */ diff --git a/src/gromacs/pulling/pullcoordexpressionparser.cpp b/src/gromacs/pulling/pullcoordexpressionparser.cpp new file mode 100644 index 0000000000..27e94e0cc9 --- /dev/null +++ b/src/gromacs/pulling/pullcoordexpressionparser.cpp @@ -0,0 +1,94 @@ +/* + * 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(); + 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 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 diff --git a/src/gromacs/pulling/pullcoordexpressionparser.h b/src/gromacs/pulling/pullcoordexpressionparser.h new file mode 100644 index 0000000000..053b38473b --- /dev/null +++ b/src/gromacs/pulling/pullcoordexpressionparser.h @@ -0,0 +1,110 @@ +/* + * 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 + * \author Paul Bauer + * \author Joe Jordan + * \author Berk Hess + * + */ +#ifndef GMX_PULL_PULLCOORDEXPRESSIONPARSER_H +#define GMX_PULL_PULLCOORDEXPRESSIONPARSER_H + +#include "config.h" + +#include +#include +#include + +#if HAVE_MUPARSER +# include +#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 +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 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 variableValues_; + + /*! \brief The parser_ which compiles and evaluates the mathematical expression */ + std::unique_ptr parser_; +}; + +} // namespace gmx + +#endif // GMX_PULL_PULLCOORDEXPRESSIONPARSER_H diff --git a/src/gromacs/pulling/tests/pull.cpp b/src/gromacs/pulling/tests/pull.cpp index 07c2d545ee..046b5668b4 100644 --- a/src/gromacs/pulling/tests/pull.cpp +++ b/src/gromacs/pulling/tests/pull.cpp @@ -43,6 +43,8 @@ #include "gromacs/pulling/pull.h" +#include "config.h" + #include #include @@ -52,6 +54,7 @@ #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" @@ -195,6 +198,185 @@ TEST_F(PullTest, MaxPullDistanceXySkewedBox) 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).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).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).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{}); + EXPECT_REAL_EQ_TOL(value, 10, defaultRealTolerance()); +} +#endif // HAVE_MUPARSER + } // namespace } // namespace gmx diff --git a/src/gromacs/pulling/transformationcoordinate.cpp b/src/gromacs/pulling/transformationcoordinate.cpp new file mode 100644 index 0000000000..cc52053af5 --- /dev/null +++ b/src/gromacs/pulling/transformationcoordinate.cpp @@ -0,0 +1,215 @@ +/* + * 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 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 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 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 diff --git a/src/gromacs/pulling/transformationcoordinate.h b/src/gromacs/pulling/transformationcoordinate.h new file mode 100644 index 0000000000..c9b6c50119 --- /dev/null +++ b/src/gromacs/pulling/transformationcoordinate.h @@ -0,0 +1,79 @@ +/* + * 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 + * \author Paul Bauer + * \author Joe Jordan + * \author Berk Hess + * + */ +#ifndef GMX_PULL_TRANSFORMATIONCOORDINATE_H +#define GMX_PULL_TRANSFORMATIONCOORDINATE_H + +struct pull_t; +struct pull_coord_work_t; + +namespace gmx +{ +template +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 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 variableCoords, + double transformationCoordForce); + +} // namespace gmx + +#endif // GMX_PULL_TRANSFORMATIONCOORDINATE_H diff --git a/src/testutils/TestMacros.cmake b/src/testutils/TestMacros.cmake index 26794bb28a..8acc6f94b7 100644 --- a/src/testutils/TestMacros.cmake +++ b/src/testutils/TestMacros.cmake @@ -172,6 +172,7 @@ function (gmx_add_gtest_executable EXENAME) # 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) -- 2.22.0