src/gromacs/mdtypes/interaction_const.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2020,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 "interaction_const.h"
  38
  39 #include <cstdio>
  40
  41 #include "gromacs/ewald/ewald_utils.h"
  42 #include "gromacs/math/functions.h"
  43 #include "gromacs/math/units.h"
  44 #include "gromacs/mdlib/rf_util.h"
  45 #include "gromacs/mdtypes/forcerec.h"
  46 #include "gromacs/mdtypes/inputrec.h"
  47 #include "gromacs/topology/topology.h"
  48 #include "gromacs/utility/fatalerror.h"
  49 #include "gromacs/utility/pleasecite.h"
  50
  51 interaction_const_t::SoftCoreParameters::SoftCoreParameters(const t_lambda& fepvals) :
  52     alphaVdw(fepvals.sc_alpha),
  53     alphaCoulomb(fepvals.bScCoul ? fepvals.sc_alpha : 0),
  54     lambdaPower(fepvals.sc_power),
  55     sigma6WithInvalidSigma(gmx::power6(fepvals.sc_sigma)),
  56     sigma6Minimum(fepvals.bScCoul ? gmx::power6(fepvals.sc_sigma_min) : 0),
  57     softcoreType(fepvals.softcoreFunction),
  58     gapsysScaleLinpointVdW(fepvals.scGapsysScaleLinpointLJ),
  59     gapsysScaleLinpointCoul(fepvals.scGapsysScaleLinpointQ),
  60     sigma6VdWGapsys(gmx::power6(fepvals.scSigmaLJGapsys))
  61 {
  62     // This is checked during tpr reading, so we can assert here
  63     GMX_RELEASE_ASSERT(fepvals.sc_r_power == 6.0, "We only support soft-core r-power 6");
  64 }
  65
  66 /*! \brief Print Coulomb Ewald citations and set ewald coefficients */
  67 static void initCoulombEwaldParameters(FILE*                fp,
  68                                        const t_inputrec&    ir,
  69                                        bool                 systemHasNetCharge,
  70                                        interaction_const_t* ic)
  71 {
  72     if (!EEL_PME_EWALD(ir.coulombtype))
  73     {
  74         return;
  75     }
  76
  77     if (fp)
  78     {
  79         fprintf(fp, "Will do PME sum in reciprocal space for electrostatic interactions.\n");
  80
  81         if (ir.coulombtype == CoulombInteractionType::P3mAD)
  82         {
  83             please_cite(fp, "Hockney1988");
  84             please_cite(fp, "Ballenegger2012");
  85         }
  86         else
  87         {
  88             please_cite(fp, "Essmann95a");
  89         }
  90
  91         if (ir.ewald_geometry == EwaldGeometry::ThreeDC)
  92         {
  93             if (fp)
  94             {
  95                 fprintf(fp,
  96                         "Using the Ewald3DC correction for systems with a slab geometry%s.\n",
  97                         systemHasNetCharge ? " and net charge" : "");
  98             }
  99             please_cite(fp, "In-Chul99a");
 100             if (systemHasNetCharge)
 101             {
 102                 please_cite(fp, "Ballenegger2009");
 103             }
 104         }
 105     }
 106
 107     ic->ewaldcoeff_q = calc_ewaldcoeff_q(ir.rcoulomb, ir.ewald_rtol);
 108     if (fp)
 109     {
 110         fprintf(fp, "Using a Gaussian width (1/beta) of %g nm for Ewald\n", 1 / ic->ewaldcoeff_q);
 111     }
 112
 113     if (ic->coulomb_modifier == InteractionModifiers::PotShift)
 114     {
 115         GMX_RELEASE_ASSERT(ic->rcoulomb != 0, "Cutoff radius cannot be zero");
 116         ic->sh_ewald = std::erfc(ic->ewaldcoeff_q * ic->rcoulomb) / ic->rcoulomb;
 117     }
 118     else
 119     {
 120         ic->sh_ewald = 0;
 121     }
 122 }
 123
 124 /*! \brief Print Van der Waals Ewald citations and set ewald coefficients */
 125 static void initVdwEwaldParameters(FILE* fp, const t_inputrec& ir, interaction_const_t* ic)
 126 {
 127     if (!EVDW_PME(ir.vdwtype))
 128     {
 129         return;
 130     }
 131
 132     if (fp)
 133     {
 134         fprintf(fp, "Will do PME sum in reciprocal space for LJ dispersion interactions.\n");
 135         please_cite(fp, "Essmann95a");
 136     }
 137     ic->ewaldcoeff_lj = calc_ewaldcoeff_lj(ir.rvdw, ir.ewald_rtol_lj);
 138     if (fp)
 139     {
 140         fprintf(fp, "Using a Gaussian width (1/beta) of %g nm for LJ Ewald\n", 1 / ic->ewaldcoeff_lj);
 141     }
 142
 143     if (ic->vdw_modifier == InteractionModifiers::PotShift)
 144     {
 145         real crc2       = gmx::square(ic->ewaldcoeff_lj * ic->rvdw);
 146         ic->sh_lj_ewald = (std::exp(-crc2) * (1 + crc2 + 0.5 * crc2 * crc2) - 1) / gmx::power6(ic->rvdw);
 147     }
 148     else
 149     {
 150         ic->sh_lj_ewald = 0;
 151     }
 152 }
 153
 154 static real calcBuckinghamBMax(FILE* fplog, const gmx_mtop_t& mtop)
 155 {
 156     const t_atoms *at1, *at2;
 157     int            i, j, tpi, tpj, ntypes;
 158     real           b, bmin;
 159
 160     if (fplog)
 161     {
 162         fprintf(fplog, "Determining largest Buckingham b parameter for table\n");
 163     }
 164     ntypes = mtop.ffparams.atnr;
 165
 166     bmin            = -1;
 167     real bham_b_max = 0;
 168     for (size_t mt1 = 0; mt1 < mtop.moltype.size(); mt1++)
 169     {
 170         at1 = &mtop.moltype[mt1].atoms;
 171         for (i = 0; (i < at1->nr); i++)
 172         {
 173             tpi = at1->atom[i].type;
 174             if (tpi >= ntypes)
 175             {
 176                 gmx_fatal(FARGS, "Atomtype[%d] = %d, maximum = %d", i, tpi, ntypes);
 177             }
 178
 179             for (size_t mt2 = mt1; mt2 < mtop.moltype.size(); mt2++)
 180             {
 181                 at2 = &mtop.moltype[mt2].atoms;
 182                 for (j = 0; (j < at2->nr); j++)
 183                 {
 184                     tpj = at2->atom[j].type;
 185                     if (tpj >= ntypes)
 186                     {
 187                         gmx_fatal(FARGS, "Atomtype[%d] = %d, maximum = %d", j, tpj, ntypes);
 188                     }
 189                     b = mtop.ffparams.iparams[tpi * ntypes + tpj].bham.b;
 190                     if (b > bham_b_max)
 191                     {
 192                         bham_b_max = b;
 193                     }
 194                     if ((b < bmin) || (bmin == -1))
 195                     {
 196                         bmin = b;
 197                     }
 198                 }
 199             }
 200         }
 201     }
 202     if (fplog)
 203     {
 204         fprintf(fplog, "Buckingham b parameters, min: %g, max: %g\n", bmin, bham_b_max);
 205     }
 206
 207     return bham_b_max;
 208 }
 209
 210 static void clear_force_switch_constants(shift_consts_t* sc)
 211 {
 212     sc->c2   = 0;
 213     sc->c3   = 0;
 214     sc->cpot = 0;
 215 }
 216
 217 static void force_switch_constants(real p, real rsw, real rc, shift_consts_t* sc)
 218 {
 219     /* Here we determine the coefficient for shifting the force to zero
 220      * between distance rsw and the cut-off rc.
 221      * For a potential of r^-p, we have force p*r^-(p+1).
 222      * But to save flops we absorb p in the coefficient.
 223      * Thus we get:
 224      * force/p   = r^-(p+1) + c2*r^2 + c3*r^3
 225      * potential = r^-p + c2/3*r^3 + c3/4*r^4 + cpot
 226      */
 227     sc->c2   = ((p + 1) * rsw - (p + 4) * rc) / (pow(rc, p + 2) * gmx::square(rc - rsw));
 228     sc->c3   = -((p + 1) * rsw - (p + 3) * rc) / (pow(rc, p + 2) * gmx::power3(rc - rsw));
 229     sc->cpot = -pow(rc, -p) + p * sc->c2 / 3 * gmx::power3(rc - rsw)
 230                + p * sc->c3 / 4 * gmx::power4(rc - rsw);
 231 }
 232
 233 static void potential_switch_constants(real rsw, real rc, switch_consts_t* sc)
 234 {
 235     /* The switch function is 1 at rsw and 0 at rc.
 236      * The derivative and second derivate are zero at both ends.
 237      * rsw        = max(r - r_switch, 0)
 238      * sw         = 1 + c3*rsw^3 + c4*rsw^4 + c5*rsw^5
 239      * dsw        = 3*c3*rsw^2 + 4*c4*rsw^3 + 5*c5*rsw^4
 240      * force      = force*dsw - potential*sw
 241      * potential *= sw
 242      */
 243     sc->c3 = -10 / gmx::power3(rc - rsw);
 244     sc->c4 = 15 / gmx::power4(rc - rsw);
 245     sc->c5 = -6 / gmx::power5(rc - rsw);
 246 }
 247
 248
 249 interaction_const_t init_interaction_const(FILE* fp, const t_inputrec& ir, const gmx_mtop_t& mtop, bool systemHasNetCharge)
 250 {
 251     interaction_const_t interactionConst;
 252
 253     interactionConst.coulombEwaldTables = std::make_unique<EwaldCorrectionTables>();
 254     interactionConst.vdwEwaldTables     = std::make_unique<EwaldCorrectionTables>();
 255
 256     /* Lennard-Jones */
 257     interactionConst.vdwtype         = ir.vdwtype;
 258     interactionConst.vdw_modifier    = ir.vdw_modifier;
 259     interactionConst.reppow          = mtop.ffparams.reppow;
 260     interactionConst.rvdw            = cutoff_inf(ir.rvdw);
 261     interactionConst.rvdw_switch     = ir.rvdw_switch;
 262     interactionConst.ljpme_comb_rule = ir.ljpme_combination_rule;
 263     interactionConst.useBuckingham   = (mtop.ffparams.functype[0] == F_BHAM);
 264     if (interactionConst.useBuckingham)
 265     {
 266         interactionConst.buckinghamBMax = calcBuckinghamBMax(fp, mtop);
 267     }
 268
 269     initVdwEwaldParameters(fp, ir, &interactionConst);
 270
 271     clear_force_switch_constants(&interactionConst.dispersion_shift);
 272     clear_force_switch_constants(&interactionConst.repulsion_shift);
 273
 274     switch (interactionConst.vdw_modifier)
 275     {
 276         case InteractionModifiers::PotShift:
 277             /* Only shift the potential, don't touch the force */
 278             interactionConst.dispersion_shift.cpot = -1.0 / gmx::power6(interactionConst.rvdw);
 279             interactionConst.repulsion_shift.cpot  = -1.0 / gmx::power12(interactionConst.rvdw);
 280             break;
 281         case InteractionModifiers::ForceSwitch:
 282             /* Switch the force, switch and shift the potential */
 283             force_switch_constants(
 284                     6.0, interactionConst.rvdw_switch, interactionConst.rvdw, &interactionConst.dispersion_shift);
 285             force_switch_constants(
 286                     12.0, interactionConst.rvdw_switch, interactionConst.rvdw, &interactionConst.repulsion_shift);
 287             break;
 288         case InteractionModifiers::PotSwitch:
 289             /* Switch the potential and force */
 290             potential_switch_constants(
 291                     interactionConst.rvdw_switch, interactionConst.rvdw, &interactionConst.vdw_switch);
 292             break;
 293         case InteractionModifiers::None:
 294         case InteractionModifiers::ExactCutoff:
 295             /* Nothing to do here */
 296             break;
 297         default: gmx_incons("unimplemented potential modifier");
 298     }
 299
 300     /* Electrostatics */
 301     interactionConst.eeltype          = ir.coulombtype;
 302     interactionConst.coulomb_modifier = ir.coulomb_modifier;
 303     interactionConst.rcoulomb         = cutoff_inf(ir.rcoulomb);
 304     interactionConst.rcoulomb_switch  = ir.rcoulomb_switch;
 305     interactionConst.epsilon_r        = ir.epsilon_r;
 306
 307     /* Set the Coulomb energy conversion factor */
 308     if (interactionConst.epsilon_r != 0)
 309     {
 310         interactionConst.epsfac = gmx::c_one4PiEps0 / interactionConst.epsilon_r;
 311     }
 312     else
 313     {
 314         /* eps = 0 is infinite dieletric: no Coulomb interactions */
 315         interactionConst.epsfac = 0;
 316     }
 317
 318     /* Reaction-field */
 319     if (EEL_RF(interactionConst.eeltype))
 320     {
 321         GMX_RELEASE_ASSERT(interactionConst.eeltype != CoulombInteractionType::GRFNotused,
 322                            "GRF is no longer supported");
 323         interactionConst.reactionFieldPermitivity = ir.epsilon_rf;
 324         calc_rffac(fp,
 325                    interactionConst.epsilon_r,
 326                    interactionConst.reactionFieldPermitivity,
 327                    interactionConst.rcoulomb,
 328                    &interactionConst.reactionFieldCoefficient,
 329                    &interactionConst.reactionFieldShift);
 330     }
 331     else
 332     {
 333         /* For plain cut-off we might use the reaction-field kernels */
 334         interactionConst.reactionFieldPermitivity = interactionConst.epsilon_r;
 335         interactionConst.reactionFieldCoefficient = 0;
 336         if (ir.coulomb_modifier == InteractionModifiers::PotShift)
 337         {
 338             interactionConst.reactionFieldShift = 1 / interactionConst.rcoulomb;
 339         }
 340         else
 341         {
 342             interactionConst.reactionFieldShift = 0;
 343         }
 344     }
 345
 346     initCoulombEwaldParameters(fp, ir, systemHasNetCharge, &interactionConst);
 347
 348     if (fp != nullptr)
 349     {
 350         real dispersion_shift;
 351
 352         dispersion_shift = interactionConst.dispersion_shift.cpot;
 353         if (EVDW_PME(interactionConst.vdwtype))
 354         {
 355             dispersion_shift -= interactionConst.sh_lj_ewald;
 356         }
 357         fprintf(fp,
 358                 "Potential shift: LJ r^-12: %.3e r^-6: %.3e",
 359                 interactionConst.repulsion_shift.cpot,
 360                 dispersion_shift);
 361
 362         if (interactionConst.eeltype == CoulombInteractionType::Cut)
 363         {
 364             fprintf(fp, ", Coulomb %.e", -interactionConst.reactionFieldShift);
 365         }
 366         else if (EEL_PME(interactionConst.eeltype))
 367         {
 368             fprintf(fp, ", Ewald %.3e", -interactionConst.sh_ewald);
 369         }
 370         fprintf(fp, "\n");
 371     }
 372
 373     if (ir.efep != FreeEnergyPerturbationType::No)
 374     {
 375         GMX_RELEASE_ASSERT(ir.fepvals, "ir.fepvals should be set wth free-energy");
 376         interactionConst.softCoreParameters =
 377                 std::make_unique<interaction_const_t::SoftCoreParameters>(*ir.fepvals);
 378     }
 379
 380     return interactionConst;
 381 }