src/gromacs/mdlib/enerdata_utils.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include "enerdata_utils.h"
  41
  42 #include "gromacs/gmxlib/network.h"
  43 #include "gromacs/mdtypes/commrec.h"
  44 #include "gromacs/mdtypes/enerdata.h"
  45 #include "gromacs/mdtypes/inputrec.h"
  46 #include "gromacs/utility/enumerationhelpers.h"
  47 #include "gromacs/utility/fatalerror.h"
  48 #include "gromacs/utility/smalloc.h"
  49
  50 ForeignLambdaTerms::ForeignLambdaTerms(int numLambdas) :
  51     numLambdas_(numLambdas), energies_(1 + numLambdas), dhdl_(1 + numLambdas)
  52 {
  53 }
  54
  55 std::pair<std::vector<double>, std::vector<double>> ForeignLambdaTerms::getTerms(const t_commrec* cr) const
  56 {
  57     GMX_RELEASE_ASSERT(finalizedPotentialContributions_,
  58                        "The object needs to be finalized before calling getTerms");
  59
  60     std::vector<double> data(2 * numLambdas_);
  61     for (int i = 0; i < numLambdas_; i++)
  62     {
  63         data[i]               = energies_[1 + i] - energies_[0];
  64         data[numLambdas_ + i] = dhdl_[1 + i];
  65     }
  66     if (cr && cr->nnodes > 1)
  67     {
  68         gmx_sumd(data.size(), data.data(), cr);
  69     }
  70     auto dataMid = data.begin() + numLambdas_;
  71
  72     return { { data.begin(), dataMid }, { dataMid, data.end() } };
  73 }
  74
  75 void ForeignLambdaTerms::zeroAllTerms()
  76 {
  77     std::fill(energies_.begin(), energies_.end(), 0.0);
  78     std::fill(dhdl_.begin(), dhdl_.end(), 0.0);
  79     finalizedPotentialContributions_ = false;
  80 }
  81
  82 gmx_enerdata_t::gmx_enerdata_t(int numEnergyGroups, int numFepLambdas) :
  83     grpp(numEnergyGroups), foreignLambdaTerms(numFepLambdas)
  84 {
  85 }
  86
  87 static real sum_v(int n, gmx::ArrayRef<const real> v)
  88 {
  89     real t;
  90     int  i;
  91
  92     t = 0.0;
  93     for (i = 0; (i < n); i++)
  94     {
  95         t = t + v[i];
  96     }
  97
  98     return t;
  99 }
 100
 101 void sum_epot(const gmx_grppairener_t& grpp, real* epot)
 102 {
 103     int i;
 104
 105     /* Accumulate energies */
 106     epot[F_COUL_SR] = sum_v(grpp.nener, grpp.energyGroupPairTerms[NonBondedEnergyTerms::CoulombSR]);
 107     epot[F_LJ]      = sum_v(grpp.nener, grpp.energyGroupPairTerms[NonBondedEnergyTerms::LJSR]);
 108     epot[F_LJ14]    = sum_v(grpp.nener, grpp.energyGroupPairTerms[NonBondedEnergyTerms::LJ14]);
 109     epot[F_COUL14]  = sum_v(grpp.nener, grpp.energyGroupPairTerms[NonBondedEnergyTerms::Coulomb14]);
 110
 111     /* lattice part of LR doesnt belong to any group
 112      * and has been added earlier
 113      */
 114     epot[F_BHAM] = sum_v(grpp.nener, grpp.energyGroupPairTerms[NonBondedEnergyTerms::BuckinghamSR]);
 115
 116     epot[F_EPOT] = 0;
 117     for (i = 0; (i < F_EPOT); i++)
 118     {
 119         if (i != F_DISRESVIOL && i != F_ORIRESDEV)
 120         {
 121             epot[F_EPOT] += epot[i];
 122         }
 123     }
 124 }
 125
 126 // Adds computed dH/dlambda contribution i to the requested output
 127 static void set_dhdl_output(gmx_enerdata_t* enerd, FreeEnergyPerturbationCouplingType i, const t_lambda& fepvals)
 128 {
 129     if (fepvals.separate_dvdl[i])
 130     {
 131         /* Translate free-energy term indices to idef term indices.
 132          * Could this be done more readably/compactly?
 133          */
 134         int index;
 135         switch (i)
 136         {
 137             case (FreeEnergyPerturbationCouplingType::Mass): index = F_DKDL; break;
 138             case (FreeEnergyPerturbationCouplingType::Coul): index = F_DVDL_COUL; break;
 139             case (FreeEnergyPerturbationCouplingType::Vdw): index = F_DVDL_VDW; break;
 140             case (FreeEnergyPerturbationCouplingType::Bonded): index = F_DVDL_BONDED; break;
 141             case (FreeEnergyPerturbationCouplingType::Restraint): index = F_DVDL_RESTRAINT; break;
 142             default: index = F_DVDL; break;
 143         }
 144         enerd->term[index] = enerd->dvdl_lin[i] + enerd->dvdl_nonlin[i];
 145         if (debug)
 146         {
 147             fprintf(debug,
 148                     "dvdl-%s[%2d]: %f: non-linear %f + linear %f\n",
 149                     enumValueToString(i),
 150                     static_cast<int>(i),
 151                     enerd->term[index],
 152                     enerd->dvdl_nonlin[i],
 153                     enerd->dvdl_lin[i]);
 154         }
 155     }
 156     else
 157     {
 158         enerd->term[F_DVDL] += enerd->dvdl_lin[i] + enerd->dvdl_nonlin[i];
 159         if (debug)
 160         {
 161             fprintf(debug,
 162                     "dvd-%sl[%2d]: %f: non-linear %f + linear %f\n",
 163                     enumValueToString(FreeEnergyPerturbationCouplingType::Fep),
 164                     static_cast<int>(i),
 165                     enerd->term[F_DVDL],
 166                     enerd->dvdl_nonlin[i],
 167                     enerd->dvdl_lin[i]);
 168         }
 169     }
 170 }
 171
 172 void ForeignLambdaTerms::addConstantDhdl(const double dhdl)
 173 {
 174     for (double& foreignDhdl : dhdl_)
 175     {
 176         foreignDhdl += dhdl;
 177     }
 178 }
 179
 180 void ForeignLambdaTerms::finalizePotentialContributions(gmx::ArrayRef<const double> dvdlLinear,
 181                                                         gmx::ArrayRef<const real>   lambda,
 182                                                         const t_lambda&             fepvals)
 183 {
 184     if (finalizedPotentialContributions_)
 185     {
 186         return;
 187     }
 188
 189     double dvdl_lin = 0;
 190     for (int i = 0; i < static_cast<int>(FreeEnergyPerturbationCouplingType::Count); i++)
 191     {
 192         dvdl_lin += dvdlLinear[i];
 193     }
 194     addConstantDhdl(dvdl_lin);
 195
 196     /* Sum the foreign lambda energy difference contributions.
 197      * Note that here we only add the potential energy components.
 198      * The constraint and kinetic energy components are add after integration
 199      * by sum_dhdl().
 200      */
 201     for (int i = 0; i < fepvals.n_lambda; i++)
 202     {
 203         /* note we are iterating over fepvals here!
 204            For the current lam, dlam = 0 automatically,
 205            so we don't need to add anything to the
 206            enerd->enerpart_lambda[0] */
 207
 208         /* we don't need to worry about dvdl_lin contributions to dE at
 209            current lambda, because the contributions to the current
 210            lambda are automatically zeroed */
 211
 212         double enerpart_lambda = 0;
 213         for (gmx::index j = 0; j < lambda.ssize(); j++)
 214         {
 215             /* Note that this loop is over all dhdl components, not just the separated ones */
 216             const double dlam = fepvals.all_lambda[j][i] - lambda[j];
 217
 218             enerpart_lambda += dlam * dvdlLinear[j];
 219         }
 220         accumulate(1 + i, enerpart_lambda, 0);
 221     }
 222
 223     finalizedPotentialContributions_ = true;
 224 }
 225
 226 void accumulatePotentialEnergies(gmx_enerdata_t* enerd, gmx::ArrayRef<const real> lambda, const t_lambda* fepvals)
 227 {
 228     sum_epot(enerd->grpp, enerd->term.data());
 229
 230     if (fepvals)
 231     {
 232         enerd->term[F_DVDL] = 0.0;
 233         for (auto i : gmx::EnumerationWrapper<FreeEnergyPerturbationCouplingType>{})
 234         {
 235             // Skip kinetic terms here, as those are not available here yet
 236             if (i != FreeEnergyPerturbationCouplingType::Mass)
 237             {
 238                 set_dhdl_output(enerd, i, *fepvals);
 239             }
 240         }
 241
 242         enerd->foreignLambdaTerms.finalizePotentialContributions(enerd->dvdl_lin, lambda, *fepvals);
 243     }
 244 }
 245
 246 void ForeignLambdaTerms::accumulateKinetic(int listIndex, double energy, double dhdl)
 247 {
 248     energies_[listIndex] += energy;
 249     dhdl_[listIndex] += dhdl;
 250 }
 251
 252 void ForeignLambdaTerms::finalizeKineticContributions(gmx::ArrayRef<const real> energyTerms,
 253                                                       const double              dhdlMass,
 254                                                       gmx::ArrayRef<const real> lambda,
 255                                                       const t_lambda&           fepvals)
 256 {
 257     // Add perturbed mass contributions
 258     addConstantDhdl(dhdlMass);
 259
 260     // Treat current lambda, the deltaH contribution is 0 as delta-lambda=0 for the current lambda
 261     accumulateKinetic(0, 0.0, energyTerms[F_DVDL_CONSTR]);
 262     if (!fepvals.separate_dvdl[FreeEnergyPerturbationCouplingType::Mass])
 263     {
 264         accumulateKinetic(0, 0.0, energyTerms[F_DKDL]);
 265     }
 266
 267     for (int i = 0; i < fepvals.n_lambda; i++)
 268     {
 269         /* Note that potential energy terms have been added by sum_epot() -> sum_dvdl() */
 270
 271         /* Constraints can not be evaluated at foreign lambdas, so we add
 272          * a linear extrapolation. This is an approximation, but usually
 273          * quite accurate since constraints change little between lambdas.
 274          */
 275         const FreeEnergyPerturbationCouplingType lambdaIndex =
 276                 (fepvals.separate_dvdl[FreeEnergyPerturbationCouplingType::Bonded]
 277                          ? FreeEnergyPerturbationCouplingType::Bonded
 278                          : FreeEnergyPerturbationCouplingType::Fep);
 279         const double dlam = fepvals.all_lambda[lambdaIndex][i] - lambda[static_cast<int>(lambdaIndex)];
 280         accumulateKinetic(1 + i, dlam * energyTerms[F_DVDL_CONSTR], energyTerms[F_DVDL_CONSTR]);
 281
 282         if (!fepvals.separate_dvdl[FreeEnergyPerturbationCouplingType::Mass])
 283         {
 284             const double dlam = fepvals.all_lambda[FreeEnergyPerturbationCouplingType::Mass][i]
 285                                 - lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Mass)];
 286             accumulateKinetic(1 + i, dlam * energyTerms[F_DKDL], energyTerms[F_DKDL]);
 287         }
 288     }
 289 }
 290
 291 void accumulateKineticLambdaComponents(gmx_enerdata_t*           enerd,
 292                                        gmx::ArrayRef<const real> lambda,
 293                                        const t_lambda&           fepvals)
 294 {
 295     if (fepvals.separate_dvdl[FreeEnergyPerturbationCouplingType::Bonded])
 296     {
 297         enerd->term[F_DVDL_BONDED] += enerd->term[F_DVDL_CONSTR];
 298     }
 299     else
 300     {
 301         enerd->term[F_DVDL] += enerd->term[F_DVDL_CONSTR];
 302     }
 303
 304     // Add computed mass dH/dlambda contribution to the requested output
 305     set_dhdl_output(enerd, FreeEnergyPerturbationCouplingType::Mass, fepvals);
 306
 307     enerd->foreignLambdaTerms.finalizeKineticContributions(
 308             enerd->term, enerd->dvdl_lin[FreeEnergyPerturbationCouplingType::Mass], lambda, fepvals);
 309
 310     /* The constrain contribution is now included in other terms, so clear it */
 311     enerd->term[F_DVDL_CONSTR] = 0;
 312 }
 313
 314 void gmx_grppairener_t::clear()
 315 {
 316     for (int i = 0; (i < static_cast<int>(NonBondedEnergyTerms::Count)); i++)
 317     {
 318         for (int j = 0; (j < nener); j++)
 319         {
 320             energyGroupPairTerms[i][j] = 0.0;
 321         }
 322     }
 323 }
 324
 325 void reset_dvdl_enerdata(gmx_enerdata_t* enerd)
 326 {
 327     for (auto i : keysOf(enerd->dvdl_lin))
 328     {
 329         enerd->dvdl_lin[i]    = 0.0;
 330         enerd->dvdl_nonlin[i] = 0.0;
 331     }
 332 }
 333
 334 void reset_enerdata(gmx_enerdata_t* enerd)
 335 {
 336     int i, j;
 337
 338     /* First reset all energy components. */
 339     for (i = 0; (i < static_cast<int>(NonBondedEnergyTerms::Count)); i++)
 340     {
 341         for (j = 0; (j < enerd->grpp.nener); j++)
 342         {
 343             enerd->grpp.energyGroupPairTerms[i][j] = 0.0_real;
 344         }
 345     }
 346
 347     /* Normal potential energy components */
 348     for (i = 0; (i <= F_EPOT); i++)
 349     {
 350         enerd->term[i] = 0.0_real;
 351     }
 352     enerd->term[F_PDISPCORR]      = 0.0_real;
 353     enerd->term[F_DVDL]           = 0.0_real;
 354     enerd->term[F_DVDL_COUL]      = 0.0_real;
 355     enerd->term[F_DVDL_VDW]       = 0.0_real;
 356     enerd->term[F_DVDL_BONDED]    = 0.0_real;
 357     enerd->term[F_DVDL_RESTRAINT] = 0.0_real;
 358     enerd->term[F_DKDL]           = 0.0_real;
 359     enerd->foreignLambdaTerms.zeroAllTerms();
 360     /* reset dvdl - separate function since it is also called elsewhere */
 361     reset_dvdl_enerdata(enerd);
 362 }