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, 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.ener[egCOULSR]);
 110     epot[F_LJ]      = sum_v(grpp.nener, grpp.ener[egLJSR]);
 111     epot[F_LJ14]    = sum_v(grpp.nener, grpp.ener[egLJ14]);
 112     epot[F_COUL14]  = sum_v(grpp.nener, grpp.ener[egCOUL14]);
 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.ener[egBHAMSR]);
 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 (int i = 0; i < efptNR; i++)
 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 (efptMASS): index = F_DKDL; break;
 144                 case (efptCOUL): index = F_DVDL_COUL; break;
 145                 case (efptVDW): index = F_DVDL_VDW; break;
 146                 case (efptBONDED): index = F_DVDL_BONDED; break;
 147                 case (efptRESTRAINT): index = F_DVDL_RESTRAINT; break;
 148                 default: index = F_DVDL; break;
 149             }
 150             enerd->term[index] = enerd->dvdl_lin[i] + enerd->dvdl_nonlin[i];
 151             if (debug)
 152             {
 153                 fprintf(debug, "dvdl-%s[%2d]: %f: non-linear %f + linear %f\n", efpt_names[i], i,
 154                         enerd->term[index], enerd->dvdl_nonlin[i], enerd->dvdl_lin[i]);
 155             }
 156         }
 157         else
 158         {
 159             enerd->term[F_DVDL] += enerd->dvdl_lin[i] + enerd->dvdl_nonlin[i];
 160             if (debug)
 161             {
 162                 fprintf(debug, "dvd-%sl[%2d]: %f: non-linear %f + linear %f\n", efpt_names[0], i,
 163                         enerd->term[F_DVDL], enerd->dvdl_nonlin[i], enerd->dvdl_lin[i]);
 164             }
 165         }
 166     }
 167 }
 168
 169 void ForeignLambdaTerms::addConstantDhdl(const double dhdl)
 170 {
 171     for (double& foreignDhdl : dhdl_)
 172     {
 173         foreignDhdl += dhdl;
 174     }
 175 }
 176
 177 void ForeignLambdaTerms::finalizePotentialContributions(gmx::ArrayRef<const double> dvdlLinear,
 178                                                         gmx::ArrayRef<const real>   lambda,
 179                                                         const t_lambda&             fepvals)
 180 {
 181     if (finalizedPotentialContributions_)
 182     {
 183         return;
 184     }
 185
 186     double dvdl_lin = 0;
 187     for (int i = 0; i < efptNR; i++)
 188     {
 189         dvdl_lin += dvdlLinear[i];
 190     }
 191     addConstantDhdl(dvdl_lin);
 192
 193     /* Sum the foreign lambda energy difference contributions.
 194      * Note that here we only add the potential energy components.
 195      * The constraint and kinetic energy components are add after integration
 196      * by sum_dhdl().
 197      */
 198     for (int i = 0; i < fepvals.n_lambda; i++)
 199     {
 200         /* note we are iterating over fepvals here!
 201            For the current lam, dlam = 0 automatically,
 202            so we don't need to add anything to the
 203            enerd->enerpart_lambda[0] */
 204
 205         /* we don't need to worry about dvdl_lin contributions to dE at
 206            current lambda, because the contributions to the current
 207            lambda are automatically zeroed */
 208
 209         double enerpart_lambda = 0;
 210         for (gmx::index j = 0; j < lambda.ssize(); j++)
 211         {
 212             /* Note that this loop is over all dhdl components, not just the separated ones */
 213             const double dlam = fepvals.all_lambda[j][i] - lambda[j];
 214
 215             enerpart_lambda += dlam * dvdlLinear[j];
 216         }
 217         accumulate(1 + i, enerpart_lambda, 0);
 218     }
 219
 220     finalizedPotentialContributions_ = true;
 221 }
 222
 223 void accumulatePotentialEnergies(gmx_enerdata_t* enerd, gmx::ArrayRef<const real> lambda, const t_lambda* fepvals)
 224 {
 225     sum_epot(enerd->grpp, enerd->term);
 226
 227     if (fepvals)
 228     {
 229         set_dvdl_output(enerd, *fepvals);
 230
 231         enerd->foreignLambdaTerms.finalizePotentialContributions(enerd->dvdl_lin, lambda, *fepvals);
 232     }
 233 }
 234
 235 void ForeignLambdaTerms::accumulateKinetic(int listIndex, double energy, double dhdl)
 236 {
 237     energies_[listIndex] += energy;
 238     dhdl_[listIndex] += dhdl;
 239 }
 240
 241 void ForeignLambdaTerms::finalizeKineticContributions(gmx::ArrayRef<const real> energyTerms,
 242                                                       const double              dhdlMass,
 243                                                       gmx::ArrayRef<const real> lambda,
 244                                                       const t_lambda&           fepvals)
 245 {
 246     // Add perturbed mass contributions
 247     addConstantDhdl(dhdlMass);
 248
 249     // Treat current lambda, the deltaH contribution is 0 as delta-lambda=0 for the current lambda
 250     accumulateKinetic(0, 0.0, energyTerms[F_DVDL_CONSTR]);
 251     if (!fepvals.separate_dvdl[efptMASS])
 252     {
 253         accumulateKinetic(0, 0.0, energyTerms[F_DKDL]);
 254     }
 255
 256     for (int i = 0; i < fepvals.n_lambda; i++)
 257     {
 258         /* Note that potential energy terms have been added by sum_epot() -> sum_dvdl() */
 259
 260         /* Constraints can not be evaluated at foreign lambdas, so we add
 261          * a linear extrapolation. This is an approximation, but usually
 262          * quite accurate since constraints change little between lambdas.
 263          */
 264         const int    lambdaIndex = (fepvals.separate_dvdl[efptBONDED] ? efptBONDED : efptFEP);
 265         const double dlam        = fepvals.all_lambda[lambdaIndex][i] - lambda[lambdaIndex];
 266         accumulateKinetic(1 + i, dlam * energyTerms[F_DVDL_CONSTR], energyTerms[F_DVDL_CONSTR]);
 267
 268         if (!fepvals.separate_dvdl[efptMASS])
 269         {
 270             const double dlam = fepvals.all_lambda[efptMASS][i] - lambda[efptMASS];
 271             accumulateKinetic(1 + i, dlam * energyTerms[F_DKDL], energyTerms[F_DKDL]);
 272         }
 273     }
 274 }
 275
 276 void accumulateKineticLambdaComponents(gmx_enerdata_t*           enerd,
 277                                        gmx::ArrayRef<const real> lambda,
 278                                        const t_lambda&           fepvals)
 279 {
 280     if (fepvals.separate_dvdl[efptBONDED])
 281     {
 282         enerd->term[F_DVDL_BONDED] += enerd->term[F_DVDL_CONSTR];
 283     }
 284     else
 285     {
 286         enerd->term[F_DVDL] += enerd->term[F_DVDL_CONSTR];
 287     }
 288
 289     enerd->foreignLambdaTerms.finalizeKineticContributions(enerd->term, enerd->dvdl_lin[efptMASS],
 290                                                            lambda, fepvals);
 291
 292     /* The constrain contribution is now included in other terms, so clear it */
 293     enerd->term[F_DVDL_CONSTR] = 0;
 294 }
 295
 296 void reset_foreign_enerdata(gmx_enerdata_t* enerd)
 297 {
 298     int i, j;
 299
 300     /* First reset all foreign energy components.  Foreign energies always called on
 301        neighbor search steps */
 302     for (i = 0; (i < egNR); i++)
 303     {
 304         for (j = 0; (j < enerd->grpp.nener); j++)
 305         {
 306             enerd->foreign_grpp.ener[i][j] = 0.0;
 307         }
 308     }
 309
 310     /* potential energy components */
 311     for (i = 0; (i <= F_EPOT); i++)
 312     {
 313         enerd->foreign_term[i] = 0.0;
 314     }
 315 }
 316
 317 void reset_dvdl_enerdata(gmx_enerdata_t* enerd)
 318 {
 319     for (int i = 0; i < efptNR; i++)
 320     {
 321         enerd->dvdl_lin[i]    = 0.0;
 322         enerd->dvdl_nonlin[i] = 0.0;
 323     }
 324 }
 325
 326 void reset_enerdata(gmx_enerdata_t* enerd)
 327 {
 328     int i, j;
 329
 330     /* First reset all energy components. */
 331     for (i = 0; (i < egNR); i++)
 332     {
 333         for (j = 0; (j < enerd->grpp.nener); j++)
 334         {
 335             enerd->grpp.ener[i][j] = 0.0_real;
 336         }
 337     }
 338
 339     /* Normal potential energy components */
 340     for (i = 0; (i <= F_EPOT); i++)
 341     {
 342         enerd->term[i] = 0.0_real;
 343     }
 344     enerd->term[F_PDISPCORR]      = 0.0_real;
 345     enerd->term[F_DVDL]           = 0.0_real;
 346     enerd->term[F_DVDL_COUL]      = 0.0_real;
 347     enerd->term[F_DVDL_VDW]       = 0.0_real;
 348     enerd->term[F_DVDL_BONDED]    = 0.0_real;
 349     enerd->term[F_DVDL_RESTRAINT] = 0.0_real;
 350     enerd->term[F_DKDL]           = 0.0_real;
 351     enerd->foreignLambdaTerms.zeroAllTerms();
 352     /* reset foreign energy data and dvdl - separate functions since they are also called elsewhere */
 353     reset_foreign_enerdata(enerd);
 354     reset_dvdl_enerdata(enerd);
 355 }