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