src/gromacs/pulling/transformationcoordinate.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2021, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 #include "gmxpre.h"
  36
  37 #include "transformationcoordinate.h"
  38
  39 #include "gromacs/utility/arrayref.h"
  40 #include "gromacs/utility/exceptions.h"
  41 #include "gromacs/utility/fatalerror.h"
  42 #include "gromacs/utility/gmxassert.h"
  43 #include "gromacs/utility/stringutil.h"
  44
  45 #include "pull_internal.h"
  46 #include "pullcoordexpressionparser.h"
  47
  48 namespace gmx
  49 {
  50
  51 namespace
  52 {
  53
  54 //! Calculates the value a for transformation pull coordinate
  55 double getTransformationPullCoordinateValue(pull_coord_work_t* coord)
  56 {
  57     const int transformationPullCoordinateIndex = coord->params.coordIndex;
  58     GMX_ASSERT(ssize(coord->transformationVariables) == transformationPullCoordinateIndex,
  59                "We need as many variables as the transformation pull coordinate index");
  60     double result = 0;
  61     try
  62     {
  63         result = coord->expressionParser.evaluate(coord->transformationVariables);
  64     }
  65 #if HAVE_MUPARSER
  66     catch (mu::Parser::exception_type& e)
  67     {
  68         GMX_THROW(InternalError(
  69                 formatString("failed to evaluate expression for transformation pull-coord%d: %s\n",
  70                              transformationPullCoordinateIndex + 1,
  71                              e.GetMsg().c_str())));
  72     }
  73 #endif
  74     catch (std::exception& e)
  75     {
  76         GMX_THROW(InternalError(
  77                 formatString("failed to evaluate expression for transformation pull-coord%d.\n"
  78                              "Last variable pull-coord-index: %d.\n"
  79                              "Message:  %s\n",
  80                              transformationPullCoordinateIndex + 1,
  81                              transformationPullCoordinateIndex + 1,
  82                              e.what())));
  83     }
  84     return result;
  85 }
  86
  87 } // namespace
  88
  89 double getTransformationPullCoordinateValue(pull_coord_work_t*                coord,
  90                                             ArrayRef<const pull_coord_work_t> variableCoords)
  91 {
  92     GMX_ASSERT(ssize(variableCoords) == coord->params.coordIndex,
  93                "We need as many variables as the transformation pull coordinate index");
  94     int coordIndex = 0;
  95     for (const auto& variableCoord : variableCoords)
  96     {
  97         coord->transformationVariables[coordIndex++] = variableCoord.spatialData.value;
  98     }
  99
 100     return getTransformationPullCoordinateValue(coord);
 101 }
 102
 103 /*! \brief Calculates and returns the derivative of a transformation pull coordinate from a dependent coordinate
 104  *
 105  * Note #1: this requires that getTransformationPullCoordinateValue() has been called
 106  * before with the current coordinates.
 107  *
 108  * Note #2: this method will not compute inner derivates. That is taken care of in the regular pull code
 109  *
 110  * \param[in] coord  The (transformation) coordinate to compute the value for
 111  * \param[in] variablePcrdIndex Pull coordinate index of a variable.
 112  */
 113 static double computeDerivativeForTransformationPullCoord(pull_coord_work_t* coord, const int variablePcrdIndex)
 114 {
 115     GMX_ASSERT(variablePcrdIndex >= 0 && variablePcrdIndex < coord->params.coordIndex,
 116                "The variable index should be in range of the transformation coordinate");
 117
 118     // epsilon for numerical differentiation.
 119     const double transformationPcrdValue = coord->spatialData.value;
 120     // Perform numerical differentiation of 1st order
 121     const double valueBackup = coord->transformationVariables[variablePcrdIndex];
 122     double       dx          = coord->params.dx;
 123     coord->transformationVariables[variablePcrdIndex] += dx;
 124     double transformationPcrdValueEps = getTransformationPullCoordinateValue(coord);
 125     double derivative                 = (transformationPcrdValueEps - transformationPcrdValue) / dx;
 126     // reset pull coordinate value
 127     coord->transformationVariables[variablePcrdIndex] = valueBackup;
 128     return derivative;
 129 }
 130
 131
 132 /*!
 133  * \brief Distributes the force from a transformation pull coordiante to the dependent pull
 134  * coordinates by computing the inner derivatives
 135  *
 136  * \param pcrd The transformation pull coord
 137  * \param variableCoords The dependent pull coordinates
 138  * \param transformationCoordForce The force to distribute
 139  */
 140 static void distributeTransformationPullCoordForce(pull_coord_work_t*               pcrd,
 141                                                    gmx::ArrayRef<pull_coord_work_t> variableCoords,
 142                                                    const double transformationCoordForce)
 143 {
 144     if (std::abs(transformationCoordForce) < 1e-9)
 145     {
 146         // the force is effectively 0. Don't proceed and distribute it recursively
 147         return;
 148     }
 149     GMX_ASSERT(pcrd->params.eGeom == PullGroupGeometry::Transformation,
 150                "We shouldn't end up here when not using a transformation pull coordinate.");
 151     GMX_ASSERT(ssize(variableCoords) == pcrd->params.coordIndex,
 152                "We should have as many variable coords as the coord index of the transformation "
 153                "coordinate");
 154
 155     // pcrd->scalarForce += transformationCoordForce;
 156     for (auto& variableCoord : variableCoords)
 157     {
 158         const double derivative =
 159                 computeDerivativeForTransformationPullCoord(pcrd, variableCoord.params.coordIndex);
 160         const double variablePcrdForce = transformationCoordForce * derivative;
 161         /* Since we loop over all pull coordinates with smaller index, there can be ones
 162          * that are not referenced by the transformation coordinate. Avoid apply forces
 163          * on those by skipping application of zero force.
 164          */
 165         if (variablePcrdForce != 0)
 166         {
 167             if (debug)
 168             {
 169                 fprintf(debug,
 170                         "Distributing force %4.4f for transformation coordinate %d to coordinate "
 171                         "%d with "
 172                         "force "
 173                         "%4.4f\n",
 174                         transformationCoordForce,
 175                         pcrd->params.coordIndex,
 176                         variableCoord.params.coordIndex,
 177                         variablePcrdForce);
 178             }
 179             // Note that we add to the force here, in case multiple biases act on the same pull
 180             // coord (although that is not recommended it should still work)
 181             variableCoord.scalarForce += variablePcrdForce;
 182             if (variableCoord.params.eGeom == PullGroupGeometry::Transformation)
 183             {
 184                 /*
 185                  * We can have a transformation pull coordinate depend on another transformation pull coordinate
 186                  * which in turn leads to inner derivatives between pull coordinates.
 187                  * Here we redistribute the force via the inner product
 188                  *
 189                  * Note that this only works properly if the lower ranked transformation pull coordinate has it's scalarForce set to zero
 190                  */
 191                 distributeTransformationPullCoordForce(
 192                         &variableCoord,
 193                         variableCoords.subArray(0, variableCoord.params.coordIndex),
 194                         variablePcrdForce);
 195             }
 196         }
 197     }
 198 }
 199
 200 void applyTransformationPullCoordForce(pull_coord_work_t*               pcrd,
 201                                        gmx::ArrayRef<pull_coord_work_t> variableCoords,
 202                                        const double                     transformationCoordForce)
 203 {
 204     pcrd->scalarForce = transformationCoordForce;
 205     // Note on why we need to call another method here:
 206     // applyTransformationPullCoordForce is the method that should be called by the rest of the pull code.
 207     // It's non-recursive and called exactly once for every transformation coordinate for every timestep.
 208     // In it, we set the force on the transformation coordinate,
 209     // then pass the force on to the other pull coordinates via the method distributeTransformationPullCoordForce.
 210     // The latter method is recursive to account for inner derivatives.
 211     // Note that we don't set the force on the top-level transformation coordinate in distributeTransformationPullCoordForce,
 212     // we only add to the force, which is why it can be recursive.
 213     distributeTransformationPullCoordForce(pcrd, variableCoords, transformationCoordForce);
 214 }
 215
 216
 217 } // namespace gmx