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/fatalerror.h"
  47 #include "gromacs/utility/smalloc.h"
  48
  49 ForeignLambdaTerms::ForeignLambdaTerms(int numLambdas) :
  50     numLambdas_(numLambdas),
  51     energies_(1 + numLambdas),
  52     dhdl_(1 + numLambdas)
  53 {
  54 }
  55
  56 std::pair<std::vector<double>, std::vector<double>> ForeignLambdaTerms::getTerms(const t_commrec* cr) const
  57 {
  58     GMX_RELEASE_ASSERT(finalizedPotentialContributions_,
  59                        "The object needs to be finalized before calling getTerms");
  60
  61     std::vector<double> data(2 * numLambdas_);
  62     for (int i = 0; i < numLambdas_; i++)
  63     {
  64         data[i]               = energies_[1 + i] - energies_[0];
  65         data[numLambdas_ + i] = dhdl_[1 + i];
  66     }
  67     if (cr && cr->nnodes > 1)
  68     {
  69         gmx_sumd(data.size(), data.data(), cr);
  70     }
  71     auto dataMid = data.begin() + numLambdas_;
  72
  73     return { { data.begin(), dataMid }, { dataMid, data.end() } };
  74 }
  75
  76 void ForeignLambdaTerms::zeroAllTerms()
  77 {
  78     std::fill(energies_.begin(), energies_.end(), 0.0);
  79     std::fill(dhdl_.begin(), dhdl_.end(), 0.0);
  80     finalizedPotentialContributions_ = false;
  81 }
  82
  83 gmx_enerdata_t::gmx_enerdata_t(int numEnergyGroups, int numFepLambdas) :
  84     grpp(numEnergyGroups),
  85     foreignLambdaTerms(numFepLambdas),
  86     foreign_grpp(numEnergyGroups)
  87 {
  88 }
  89
  90 static real sum_v(int n, gmx::ArrayRef<const real> v)
  91 {
  92     real t;
  93     int  i;
  94
  95     t = 0.0;
  96     for (i = 0; (i < n); i++)
  97     {
  98         t = t + v[i];
  99     }
 100
 101     return t;
 102 }
 103
 104 void sum_epot(const gmx_grppairener_t& grpp, real* epot)
 105 {
 106     int i;
 107
 108     /* Accumulate energies */
 109     epot[F_COUL_SR] = sum_v(grpp.nener, grpp.energyGroupPairTerms[NonBondedEnergyTerms::CoulombSR]);
 110     epot[F_LJ]      = sum_v(grpp.nener, grpp.energyGroupPairTerms[NonBondedEnergyTerms::LJSR]);
 111     epot[F_LJ14]    = sum_v(grpp.nener, grpp.energyGroupPairTerms[NonBondedEnergyTerms::LJ14]);
 112     epot[F_COUL14]  = sum_v(grpp.nener, grpp.energyGroupPairTerms[NonBondedEnergyTerms::Coulomb14]);
 113
 114     /* lattice part of LR doesnt belong to any group
 115      * and has been added earlier
 116      */
 117     epot[F_BHAM] = sum_v(grpp.nener, grpp.energyGroupPairTerms[NonBondedEnergyTerms::BuckinghamSR]);
 118
 119     epot[F_EPOT] = 0;
 120     for (i = 0; (i < F_EPOT); i++)
 121     {
 122         if (i != F_DISRESVIOL && i != F_ORIRESDEV)
 123         {
 124             epot[F_EPOT] += epot[i];
 125         }
 126     }
 127 }
 128
 129 // Adds computed dV/dlambda contributions to the requested outputs
 130 static void set_dvdl_output(gmx_enerdata_t* enerd, const t_lambda& fepvals)
 131 {
 132     enerd->term[F_DVDL] = 0.0;
 133     for (auto i : keysOf(fepvals.separate_dvdl))
 134     {
 135         if (fepvals.separate_dvdl[i])
 136         {
 137             /* Translate free-energy term indices to idef term indices.
 138              * Could this be done more readably/compactly?
 139              */
 140             int index;
 141             switch (i)
 142             {
 143                 case (FreeEnergyPerturbationCouplingType::Mass): index = F_DKDL; break;
 144                 case (FreeEnergyPerturbationCouplingType::Coul): index = F_DVDL_COUL; break;
 145                 case (FreeEnergyPerturbationCouplingType::Vdw): index = F_DVDL_VDW; break;
 146                 case (FreeEnergyPerturbationCouplingType::Bonded): index = F_DVDL_BONDED; break;
 147                 case (FreeEnergyPerturbationCouplingType::Restraint):
 148                     index = F_DVDL_RESTRAINT;
 149                     break;
 150                 default: index = F_DVDL; break;
 151             }
 152             enerd->term[index] = enerd->dvdl_lin[i] + enerd->dvdl_nonlin[i];
 153             if (debug)
 154             {
 155                 fprintf(debug,
 156                         "dvdl-%s[%2d]: %f: non-linear %f + linear %f\n",
 157                         enumValueToString(i),
 158                         static_cast<int>(i),
 159                         enerd->term[index],
 160                         enerd->dvdl_nonlin[i],
 161                         enerd->dvdl_lin[i]);
 162             }
 163         }
 164         else
 165         {
 166             enerd->term[F_DVDL] += enerd->dvdl_lin[i] + enerd->dvdl_nonlin[i];
 167             if (debug)
 168             {
 169                 fprintf(debug,
 170                         "dvd-%sl[%2d]: %f: non-linear %f + linear %f\n",
 171                         enumValueToString(FreeEnergyPerturbationCouplingType::Fep),
 172                         static_cast<int>(i),
 173                         enerd->term[F_DVDL],
 174                         enerd->dvdl_nonlin[i],
 175                         enerd->dvdl_lin[i]);
 176             }
 177         }
 178     }
 179 }
 180
 181 void ForeignLambdaTerms::addConstantDhdl(const double dhdl)
 182 {
 183     for (double& foreignDhdl : dhdl_)
 184     {
 185         foreignDhdl += dhdl;
 186     }
 187 }
 188
 189 void ForeignLambdaTerms::finalizePotentialContributions(gmx::ArrayRef<const double> dvdlLinear,
 190                                                         gmx::ArrayRef<const real>   lambda,
 191                                                         const t_lambda&             fepvals)
 192 {
 193     if (finalizedPotentialContributions_)
 194     {
 195         return;
 196     }
 197
 198     double dvdl_lin = 0;
 199     for (int i = 0; i < static_cast<int>(FreeEnergyPerturbationCouplingType::Count); i++)
 200     {
 201         dvdl_lin += dvdlLinear[i];
 202     }
 203     addConstantDhdl(dvdl_lin);
 204
 205     /* Sum the foreign lambda energy difference contributions.
 206      * Note that here we only add the potential energy components.
 207      * The constraint and kinetic energy components are add after integration
 208      * by sum_dhdl().
 209      */
 210     for (int i = 0; i < fepvals.n_lambda; i++)
 211     {
 212         /* note we are iterating over fepvals here!
 213            For the current lam, dlam = 0 automatically,
 214            so we don't need to add anything to the
 215            enerd->enerpart_lambda[0] */
 216
 217         /* we don't need to worry about dvdl_lin contributions to dE at
 218            current lambda, because the contributions to the current
 219            lambda are automatically zeroed */
 220
 221         double enerpart_lambda = 0;
 222         for (gmx::index j = 0; j < lambda.ssize(); j++)
 223         {
 224             /* Note that this loop is over all dhdl components, not just the separated ones */
 225             const double dlam = fepvals.all_lambda[j][i] - lambda[j];
 226
 227             enerpart_lambda += dlam * dvdlLinear[j];
 228         }
 229         accumulate(1 + i, enerpart_lambda, 0);
 230     }
 231
 232     finalizedPotentialContributions_ = true;
 233 }
 234
 235 void accumulatePotentialEnergies(gmx_enerdata_t* enerd, gmx::ArrayRef<const real> lambda, const t_lambda* fepvals)
 236 {
 237     sum_epot(enerd->grpp, enerd->term.data());
 238
 239     if (fepvals)
 240     {
 241         set_dvdl_output(enerd, *fepvals);
 242
 243         enerd->foreignLambdaTerms.finalizePotentialContributions(enerd->dvdl_lin, lambda, *fepvals);
 244     }
 245 }
 246
 247 void ForeignLambdaTerms::accumulateKinetic(int listIndex, double energy, double dhdl)
 248 {
 249     energies_[listIndex] += energy;
 250     dhdl_[listIndex] += dhdl;
 251 }
 252
 253 void ForeignLambdaTerms::finalizeKineticContributions(gmx::ArrayRef<const real> energyTerms,
 254                                                       const double              dhdlMass,
 255                                                       gmx::ArrayRef<const real> lambda,
 256                                                       const t_lambda&           fepvals)
 257 {
 258     // Add perturbed mass contributions
 259     addConstantDhdl(dhdlMass);
 260
 261     // Treat current lambda, the deltaH contribution is 0 as delta-lambda=0 for the current lambda
 262     accumulateKinetic(0, 0.0, energyTerms[F_DVDL_CONSTR]);
 263     if (!fepvals.separate_dvdl[FreeEnergyPerturbationCouplingType::Mass])
 264     {
 265         accumulateKinetic(0, 0.0, energyTerms[F_DKDL]);
 266     }
 267
 268     for (int i = 0; i < fepvals.n_lambda; i++)
 269     {
 270         /* Note that potential energy terms have been added by sum_epot() -> sum_dvdl() */
 271
 272         /* Constraints can not be evaluated at foreign lambdas, so we add
 273          * a linear extrapolation. This is an approximation, but usually
 274          * quite accurate since constraints change little between lambdas.
 275          */
 276         const FreeEnergyPerturbationCouplingType lambdaIndex =
 277                 (fepvals.separate_dvdl[FreeEnergyPerturbationCouplingType::Bonded]
 278                          ? FreeEnergyPerturbationCouplingType::Bonded
 279                          : FreeEnergyPerturbationCouplingType::Fep);
 280         const double dlam = fepvals.all_lambda[lambdaIndex][i] - lambda[static_cast<int>(lambdaIndex)];
 281         accumulateKinetic(1 + i, dlam * energyTerms[F_DVDL_CONSTR], energyTerms[F_DVDL_CONSTR]);
 282
 283         if (!fepvals.separate_dvdl[FreeEnergyPerturbationCouplingType::Mass])
 284         {
 285             const double dlam = fepvals.all_lambda[FreeEnergyPerturbationCouplingType::Mass][i]
 286                                 - lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Mass)];
 287             accumulateKinetic(1 + i, dlam * energyTerms[F_DKDL], energyTerms[F_DKDL]);
 288         }
 289     }
 290 }
 291
 292 void accumulateKineticLambdaComponents(gmx_enerdata_t*           enerd,
 293                                        gmx::ArrayRef<const real> lambda,
 294                                        const t_lambda&           fepvals)
 295 {
 296     if (fepvals.separate_dvdl[FreeEnergyPerturbationCouplingType::Bonded])
 297     {
 298         enerd->term[F_DVDL_BONDED] += enerd->term[F_DVDL_CONSTR];
 299     }
 300     else
 301     {
 302         enerd->term[F_DVDL] += enerd->term[F_DVDL_CONSTR];
 303     }
 304
 305     enerd->foreignLambdaTerms.finalizeKineticContributions(
 306             enerd->term, enerd->dvdl_lin[FreeEnergyPerturbationCouplingType::Mass], lambda, fepvals);
 307
 308     /* The constrain contribution is now included in other terms, so clear it */
 309     enerd->term[F_DVDL_CONSTR] = 0;
 310 }
 311
 312 void reset_foreign_enerdata(gmx_enerdata_t* enerd)
 313 {
 314     int i, j;
 315
 316     /* First reset all foreign energy components.  Foreign energies always called on
 317        neighbor search steps */
 318     for (i = 0; (i < static_cast<int>(NonBondedEnergyTerms::Count)); i++)
 319     {
 320         for (j = 0; (j < enerd->grpp.nener); j++)
 321         {
 322             enerd->foreign_grpp.energyGroupPairTerms[i][j] = 0.0;
 323         }
 324     }
 325
 326     /* potential energy components */
 327     for (i = 0; (i <= F_EPOT); i++)
 328     {
 329         enerd->foreign_term[i] = 0.0;
 330     }
 331 }
 332
 333 void reset_dvdl_enerdata(gmx_enerdata_t* enerd)
 334 {
 335     for (auto i : keysOf(enerd->dvdl_lin))
 336     {
 337         enerd->dvdl_lin[i]    = 0.0;
 338         enerd->dvdl_nonlin[i] = 0.0;
 339     }
 340 }
 341
 342 void reset_enerdata(gmx_enerdata_t* enerd)
 343 {
 344     int i, j;
 345
 346     /* First reset all energy components. */
 347     for (i = 0; (i < static_cast<int>(NonBondedEnergyTerms::Count)); i++)
 348     {
 349         for (j = 0; (j < enerd->grpp.nener); j++)
 350         {
 351             enerd->grpp.energyGroupPairTerms[i][j] = 0.0_real;
 352         }
 353     }
 354
 355     /* Normal potential energy components */
 356     for (i = 0; (i <= F_EPOT); i++)
 357     {
 358         enerd->term[i] = 0.0_real;
 359     }
 360     enerd->term[F_PDISPCORR]      = 0.0_real;
 361     enerd->term[F_DVDL]           = 0.0_real;
 362     enerd->term[F_DVDL_COUL]      = 0.0_real;
 363     enerd->term[F_DVDL_VDW]       = 0.0_real;
 364     enerd->term[F_DVDL_BONDED]    = 0.0_real;
 365     enerd->term[F_DVDL_RESTRAINT] = 0.0_real;
 366     enerd->term[F_DKDL]           = 0.0_real;
 367     enerd->foreignLambdaTerms.zeroAllTerms();
 368     /* reset foreign energy data and dvdl - separate functions since they are also called elsewhere */
 369     reset_foreign_enerdata(enerd);
 370     reset_dvdl_enerdata(enerd);
 371 }