From 1ff40130c200fc00f3d689a1627501c796b3468b Mon Sep 17 00:00:00 2001 From: Joe Jordan Date: Mon, 12 Apr 2021 13:36:26 +0000 Subject: [PATCH] Expose init_interaction_const and move to correct header/source --- src/gromacs/mdlib/forcerec.cpp | 337 +--------------------- src/gromacs/mdtypes/forcerec.h | 2 + src/gromacs/mdtypes/interaction_const.cpp | 326 ++++++++++++++++++++- src/gromacs/mdtypes/interaction_const.h | 15 + 4 files changed, 346 insertions(+), 334 deletions(-) diff --git a/src/gromacs/mdlib/forcerec.cpp b/src/gromacs/mdlib/forcerec.cpp index 7628ce37d1..ad147cf911 100644 --- a/src/gromacs/mdlib/forcerec.cpp +++ b/src/gromacs/mdlib/forcerec.cpp @@ -40,7 +40,6 @@ #include "config.h" -#include #include #include #include @@ -52,7 +51,6 @@ #include "gromacs/domdec/domdec.h" #include "gromacs/domdec/domdec_struct.h" #include "gromacs/ewald/ewald.h" -#include "gromacs/ewald/ewald_utils.h" #include "gromacs/ewald/pme_pp_comm_gpu.h" #include "gromacs/fileio/filetypes.h" #include "gromacs/gmxlib/network.h" @@ -63,7 +61,6 @@ #include "gromacs/listed_forces/listed_forces.h" #include "gromacs/listed_forces/pairs.h" #include "gromacs/math/functions.h" -#include "gromacs/math/units.h" #include "gromacs/math/utilities.h" #include "gromacs/math/vec.h" #include "gromacs/mdlib/dispersioncorrection.h" @@ -71,7 +68,6 @@ #include "gromacs/mdlib/forcerec_threading.h" #include "gromacs/mdlib/gmx_omp_nthreads.h" #include "gromacs/mdlib/md_support.h" -#include "gromacs/mdlib/rf_util.h" #include "gromacs/mdlib/wall.h" #include "gromacs/mdlib/wholemoleculetransform.h" #include "gromacs/mdtypes/commrec.h" @@ -434,62 +430,6 @@ static bool set_chargesum(FILE* log, t_forcerec* fr, const gmx_mtop_t& mtop) return (std::abs(fr->qsum[0]) > 1e-4 || std::abs(fr->qsum[1]) > 1e-4); } -static real calcBuckinghamBMax(FILE* fplog, const gmx_mtop_t& mtop) -{ - const t_atoms *at1, *at2; - int i, j, tpi, tpj, ntypes; - real b, bmin; - - if (fplog) - { - fprintf(fplog, "Determining largest Buckingham b parameter for table\n"); - } - ntypes = mtop.ffparams.atnr; - - bmin = -1; - real bham_b_max = 0; - for (size_t mt1 = 0; mt1 < mtop.moltype.size(); mt1++) - { - at1 = &mtop.moltype[mt1].atoms; - for (i = 0; (i < at1->nr); i++) - { - tpi = at1->atom[i].type; - if (tpi >= ntypes) - { - gmx_fatal(FARGS, "Atomtype[%d] = %d, maximum = %d", i, tpi, ntypes); - } - - for (size_t mt2 = mt1; mt2 < mtop.moltype.size(); mt2++) - { - at2 = &mtop.moltype[mt2].atoms; - for (j = 0; (j < at2->nr); j++) - { - tpj = at2->atom[j].type; - if (tpj >= ntypes) - { - gmx_fatal(FARGS, "Atomtype[%d] = %d, maximum = %d", j, tpj, ntypes); - } - b = mtop.ffparams.iparams[tpi * ntypes + tpj].bham.b; - if (b > bham_b_max) - { - bham_b_max = b; - } - if ((b < bmin) || (bmin == -1)) - { - bmin = b; - } - } - } - } - } - if (fplog) - { - fprintf(fplog, "Buckingham b parameters, min: %g, max: %g\n", bmin, bham_b_max); - } - - return bham_b_max; -} - /*!\brief If there's bonded interactions of type \c ftype1 or \c * ftype2 present in the topology, build an array of the number of * interactions present for each bonded interaction index found in the @@ -624,104 +564,6 @@ void forcerec_set_ranges(t_forcerec* fr, int natoms_force, int natoms_force_cons } } -static real cutoff_inf(real cutoff) -{ - if (cutoff == 0) - { - cutoff = GMX_CUTOFF_INF; - } - - return cutoff; -} - -/*! \brief Print Coulomb Ewald citations and set ewald coefficients */ -static void initCoulombEwaldParameters(FILE* fp, - const t_inputrec& ir, - bool systemHasNetCharge, - interaction_const_t* ic) -{ - if (!EEL_PME_EWALD(ir.coulombtype)) - { - return; - } - - if (fp) - { - fprintf(fp, "Will do PME sum in reciprocal space for electrostatic interactions.\n"); - - if (ir.coulombtype == CoulombInteractionType::P3mAD) - { - please_cite(fp, "Hockney1988"); - please_cite(fp, "Ballenegger2012"); - } - else - { - please_cite(fp, "Essmann95a"); - } - - if (ir.ewald_geometry == EwaldGeometry::ThreeDC) - { - if (fp) - { - fprintf(fp, - "Using the Ewald3DC correction for systems with a slab geometry%s.\n", - systemHasNetCharge ? " and net charge" : ""); - } - please_cite(fp, "In-Chul99a"); - if (systemHasNetCharge) - { - please_cite(fp, "Ballenegger2009"); - } - } - } - - ic->ewaldcoeff_q = calc_ewaldcoeff_q(ir.rcoulomb, ir.ewald_rtol); - if (fp) - { - fprintf(fp, "Using a Gaussian width (1/beta) of %g nm for Ewald\n", 1 / ic->ewaldcoeff_q); - } - - if (ic->coulomb_modifier == InteractionModifiers::PotShift) - { - GMX_RELEASE_ASSERT(ic->rcoulomb != 0, "Cutoff radius cannot be zero"); - ic->sh_ewald = std::erfc(ic->ewaldcoeff_q * ic->rcoulomb) / ic->rcoulomb; - } - else - { - ic->sh_ewald = 0; - } -} - -/*! \brief Print Van der Waals Ewald citations and set ewald coefficients */ -static void initVdwEwaldParameters(FILE* fp, const t_inputrec& ir, interaction_const_t* ic) -{ - if (!EVDW_PME(ir.vdwtype)) - { - return; - } - - if (fp) - { - fprintf(fp, "Will do PME sum in reciprocal space for LJ dispersion interactions.\n"); - please_cite(fp, "Essmann95a"); - } - ic->ewaldcoeff_lj = calc_ewaldcoeff_lj(ir.rvdw, ir.ewald_rtol_lj); - if (fp) - { - fprintf(fp, "Using a Gaussian width (1/beta) of %g nm for LJ Ewald\n", 1 / ic->ewaldcoeff_lj); - } - - if (ic->vdw_modifier == InteractionModifiers::PotShift) - { - real crc2 = gmx::square(ic->ewaldcoeff_lj * ic->rvdw); - ic->sh_lj_ewald = (std::exp(-crc2) * (1 + crc2 + 0.5 * crc2 * crc2) - 1) / gmx::power6(ic->rvdw); - } - else - { - ic->sh_lj_ewald = 0; - } -} - /* Generate Coulomb and/or Van der Waals Ewald long-range correction tables * * Tables are generated for one or both, depending on if the pointers @@ -785,185 +627,14 @@ void init_interaction_const_tables(FILE* fp, interaction_const_t* ic, const real } } -static void clear_force_switch_constants(shift_consts_t* sc) +real cutoff_inf(real cutoff) { - sc->c2 = 0; - sc->c3 = 0; - sc->cpot = 0; -} - -static void force_switch_constants(real p, real rsw, real rc, shift_consts_t* sc) -{ - /* Here we determine the coefficient for shifting the force to zero - * between distance rsw and the cut-off rc. - * For a potential of r^-p, we have force p*r^-(p+1). - * But to save flops we absorb p in the coefficient. - * Thus we get: - * force/p = r^-(p+1) + c2*r^2 + c3*r^3 - * potential = r^-p + c2/3*r^3 + c3/4*r^4 + cpot - */ - sc->c2 = ((p + 1) * rsw - (p + 4) * rc) / (pow(rc, p + 2) * gmx::square(rc - rsw)); - sc->c3 = -((p + 1) * rsw - (p + 3) * rc) / (pow(rc, p + 2) * gmx::power3(rc - rsw)); - sc->cpot = -pow(rc, -p) + p * sc->c2 / 3 * gmx::power3(rc - rsw) - + p * sc->c3 / 4 * gmx::power4(rc - rsw); -} - -static void potential_switch_constants(real rsw, real rc, switch_consts_t* sc) -{ - /* The switch function is 1 at rsw and 0 at rc. - * The derivative and second derivate are zero at both ends. - * rsw = max(r - r_switch, 0) - * sw = 1 + c3*rsw^3 + c4*rsw^4 + c5*rsw^5 - * dsw = 3*c3*rsw^2 + 4*c4*rsw^3 + 5*c5*rsw^4 - * force = force*dsw - potential*sw - * potential *= sw - */ - sc->c3 = -10 / gmx::power3(rc - rsw); - sc->c4 = 15 / gmx::power4(rc - rsw); - sc->c5 = -6 / gmx::power5(rc - rsw); -} - -/*! \brief Construct interaction constants - * - * This data is used (particularly) by search and force code for - * short-range interactions. Many of these are constant for the whole - * simulation; some are constant only after PME tuning completes. - */ -static interaction_const_t init_interaction_const(FILE* fp, - const t_inputrec& ir, - const gmx_mtop_t& mtop, - bool systemHasNetCharge) -{ - interaction_const_t interactionConst; - - interactionConst.coulombEwaldTables = std::make_unique(); - interactionConst.vdwEwaldTables = std::make_unique(); - - /* Lennard-Jones */ - interactionConst.vdwtype = ir.vdwtype; - interactionConst.vdw_modifier = ir.vdw_modifier; - interactionConst.reppow = mtop.ffparams.reppow; - interactionConst.rvdw = cutoff_inf(ir.rvdw); - interactionConst.rvdw_switch = ir.rvdw_switch; - interactionConst.ljpme_comb_rule = ir.ljpme_combination_rule; - interactionConst.useBuckingham = (mtop.ffparams.functype[0] == F_BHAM); - if (interactionConst.useBuckingham) - { - interactionConst.buckinghamBMax = calcBuckinghamBMax(fp, mtop); - } - - initVdwEwaldParameters(fp, ir, &interactionConst); - - clear_force_switch_constants(&interactionConst.dispersion_shift); - clear_force_switch_constants(&interactionConst.repulsion_shift); - - switch (interactionConst.vdw_modifier) - { - case InteractionModifiers::PotShift: - /* Only shift the potential, don't touch the force */ - interactionConst.dispersion_shift.cpot = -1.0 / gmx::power6(interactionConst.rvdw); - interactionConst.repulsion_shift.cpot = -1.0 / gmx::power12(interactionConst.rvdw); - break; - case InteractionModifiers::ForceSwitch: - /* Switch the force, switch and shift the potential */ - force_switch_constants( - 6.0, interactionConst.rvdw_switch, interactionConst.rvdw, &interactionConst.dispersion_shift); - force_switch_constants( - 12.0, interactionConst.rvdw_switch, interactionConst.rvdw, &interactionConst.repulsion_shift); - break; - case InteractionModifiers::PotSwitch: - /* Switch the potential and force */ - potential_switch_constants( - interactionConst.rvdw_switch, interactionConst.rvdw, &interactionConst.vdw_switch); - break; - case InteractionModifiers::None: - case InteractionModifiers::ExactCutoff: - /* Nothing to do here */ - break; - default: gmx_incons("unimplemented potential modifier"); - } - - /* Electrostatics */ - interactionConst.eeltype = ir.coulombtype; - interactionConst.coulomb_modifier = ir.coulomb_modifier; - interactionConst.rcoulomb = cutoff_inf(ir.rcoulomb); - interactionConst.rcoulomb_switch = ir.rcoulomb_switch; - interactionConst.epsilon_r = ir.epsilon_r; - - /* Set the Coulomb energy conversion factor */ - if (interactionConst.epsilon_r != 0) - { - interactionConst.epsfac = gmx::c_one4PiEps0 / interactionConst.epsilon_r; - } - else - { - /* eps = 0 is infinite dieletric: no Coulomb interactions */ - interactionConst.epsfac = 0; - } - - /* Reaction-field */ - if (EEL_RF(interactionConst.eeltype)) - { - GMX_RELEASE_ASSERT(interactionConst.eeltype != CoulombInteractionType::GRFNotused, - "GRF is no longer supported"); - interactionConst.reactionFieldPermitivity = ir.epsilon_rf; - calc_rffac(fp, - interactionConst.epsilon_r, - interactionConst.reactionFieldPermitivity, - interactionConst.rcoulomb, - &interactionConst.reactionFieldCoefficient, - &interactionConst.reactionFieldShift); - } - else - { - /* For plain cut-off we might use the reaction-field kernels */ - interactionConst.reactionFieldPermitivity = interactionConst.epsilon_r; - interactionConst.reactionFieldCoefficient = 0; - if (ir.coulomb_modifier == InteractionModifiers::PotShift) - { - interactionConst.reactionFieldShift = 1 / interactionConst.rcoulomb; - } - else - { - interactionConst.reactionFieldShift = 0; - } - } - - initCoulombEwaldParameters(fp, ir, systemHasNetCharge, &interactionConst); - - if (fp != nullptr) - { - real dispersion_shift; - - dispersion_shift = interactionConst.dispersion_shift.cpot; - if (EVDW_PME(interactionConst.vdwtype)) - { - dispersion_shift -= interactionConst.sh_lj_ewald; - } - fprintf(fp, - "Potential shift: LJ r^-12: %.3e r^-6: %.3e", - interactionConst.repulsion_shift.cpot, - dispersion_shift); - - if (interactionConst.eeltype == CoulombInteractionType::Cut) - { - fprintf(fp, ", Coulomb %.e", -interactionConst.reactionFieldShift); - } - else if (EEL_PME(interactionConst.eeltype)) - { - fprintf(fp, ", Ewald %.3e", -interactionConst.sh_ewald); - } - fprintf(fp, "\n"); - } - - if (ir.efep != FreeEnergyPerturbationType::No) + if (cutoff == 0) { - GMX_RELEASE_ASSERT(ir.fepvals, "ir.fepvals should be set wth free-energy"); - interactionConst.softCoreParameters = - std::make_unique(*ir.fepvals); + cutoff = GMX_CUTOFF_INF; } - return interactionConst; + return cutoff; } void init_forcerec(FILE* fplog, diff --git a/src/gromacs/mdtypes/forcerec.h b/src/gromacs/mdtypes/forcerec.h index 414e533a6f..218c7b26c9 100644 --- a/src/gromacs/mdtypes/forcerec.h +++ b/src/gromacs/mdtypes/forcerec.h @@ -106,6 +106,8 @@ class WholeMoleculeTransform; * this value should be slighlty smaller than sqrt(GMX_FLOAT_MAX). */ #define GMX_CUTOFF_INF 1E+18 +//! Check the cuttoff +real cutoff_inf(real cutoff); struct cginfo_mb_t { diff --git a/src/gromacs/mdtypes/interaction_const.cpp b/src/gromacs/mdtypes/interaction_const.cpp index 230dca2bff..9353a24047 100644 --- a/src/gromacs/mdtypes/interaction_const.cpp +++ b/src/gromacs/mdtypes/interaction_const.cpp @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2020, by the GROMACS development team, led by + * Copyright (c) 2020,2021, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -38,8 +38,15 @@ #include +#include "gromacs/ewald/ewald_utils.h" #include "gromacs/math/functions.h" +#include "gromacs/math/units.h" +#include "gromacs/mdlib/rf_util.h" +#include "gromacs/mdtypes/forcerec.h" #include "gromacs/mdtypes/inputrec.h" +#include "gromacs/topology/topology.h" +#include "gromacs/utility/fatalerror.h" +#include "gromacs/utility/pleasecite.h" interaction_const_t::SoftCoreParameters::SoftCoreParameters(const t_lambda& fepvals) : alphaVdw(fepvals.sc_alpha), @@ -51,3 +58,320 @@ interaction_const_t::SoftCoreParameters::SoftCoreParameters(const t_lambda& fepv // This is checked during tpr reading, so we can assert here GMX_RELEASE_ASSERT(fepvals.sc_r_power == 6.0, "We only support soft-core r-power 6"); } + +/*! \brief Print Coulomb Ewald citations and set ewald coefficients */ +static void initCoulombEwaldParameters(FILE* fp, + const t_inputrec& ir, + bool systemHasNetCharge, + interaction_const_t* ic) +{ + if (!EEL_PME_EWALD(ir.coulombtype)) + { + return; + } + + if (fp) + { + fprintf(fp, "Will do PME sum in reciprocal space for electrostatic interactions.\n"); + + if (ir.coulombtype == CoulombInteractionType::P3mAD) + { + please_cite(fp, "Hockney1988"); + please_cite(fp, "Ballenegger2012"); + } + else + { + please_cite(fp, "Essmann95a"); + } + + if (ir.ewald_geometry == EwaldGeometry::ThreeDC) + { + if (fp) + { + fprintf(fp, + "Using the Ewald3DC correction for systems with a slab geometry%s.\n", + systemHasNetCharge ? " and net charge" : ""); + } + please_cite(fp, "In-Chul99a"); + if (systemHasNetCharge) + { + please_cite(fp, "Ballenegger2009"); + } + } + } + + ic->ewaldcoeff_q = calc_ewaldcoeff_q(ir.rcoulomb, ir.ewald_rtol); + if (fp) + { + fprintf(fp, "Using a Gaussian width (1/beta) of %g nm for Ewald\n", 1 / ic->ewaldcoeff_q); + } + + if (ic->coulomb_modifier == InteractionModifiers::PotShift) + { + GMX_RELEASE_ASSERT(ic->rcoulomb != 0, "Cutoff radius cannot be zero"); + ic->sh_ewald = std::erfc(ic->ewaldcoeff_q * ic->rcoulomb) / ic->rcoulomb; + } + else + { + ic->sh_ewald = 0; + } +} + +/*! \brief Print Van der Waals Ewald citations and set ewald coefficients */ +static void initVdwEwaldParameters(FILE* fp, const t_inputrec& ir, interaction_const_t* ic) +{ + if (!EVDW_PME(ir.vdwtype)) + { + return; + } + + if (fp) + { + fprintf(fp, "Will do PME sum in reciprocal space for LJ dispersion interactions.\n"); + please_cite(fp, "Essmann95a"); + } + ic->ewaldcoeff_lj = calc_ewaldcoeff_lj(ir.rvdw, ir.ewald_rtol_lj); + if (fp) + { + fprintf(fp, "Using a Gaussian width (1/beta) of %g nm for LJ Ewald\n", 1 / ic->ewaldcoeff_lj); + } + + if (ic->vdw_modifier == InteractionModifiers::PotShift) + { + real crc2 = gmx::square(ic->ewaldcoeff_lj * ic->rvdw); + ic->sh_lj_ewald = (std::exp(-crc2) * (1 + crc2 + 0.5 * crc2 * crc2) - 1) / gmx::power6(ic->rvdw); + } + else + { + ic->sh_lj_ewald = 0; + } +} + +static real calcBuckinghamBMax(FILE* fplog, const gmx_mtop_t& mtop) +{ + const t_atoms *at1, *at2; + int i, j, tpi, tpj, ntypes; + real b, bmin; + + if (fplog) + { + fprintf(fplog, "Determining largest Buckingham b parameter for table\n"); + } + ntypes = mtop.ffparams.atnr; + + bmin = -1; + real bham_b_max = 0; + for (size_t mt1 = 0; mt1 < mtop.moltype.size(); mt1++) + { + at1 = &mtop.moltype[mt1].atoms; + for (i = 0; (i < at1->nr); i++) + { + tpi = at1->atom[i].type; + if (tpi >= ntypes) + { + gmx_fatal(FARGS, "Atomtype[%d] = %d, maximum = %d", i, tpi, ntypes); + } + + for (size_t mt2 = mt1; mt2 < mtop.moltype.size(); mt2++) + { + at2 = &mtop.moltype[mt2].atoms; + for (j = 0; (j < at2->nr); j++) + { + tpj = at2->atom[j].type; + if (tpj >= ntypes) + { + gmx_fatal(FARGS, "Atomtype[%d] = %d, maximum = %d", j, tpj, ntypes); + } + b = mtop.ffparams.iparams[tpi * ntypes + tpj].bham.b; + if (b > bham_b_max) + { + bham_b_max = b; + } + if ((b < bmin) || (bmin == -1)) + { + bmin = b; + } + } + } + } + } + if (fplog) + { + fprintf(fplog, "Buckingham b parameters, min: %g, max: %g\n", bmin, bham_b_max); + } + + return bham_b_max; +} + +static void clear_force_switch_constants(shift_consts_t* sc) +{ + sc->c2 = 0; + sc->c3 = 0; + sc->cpot = 0; +} + +static void force_switch_constants(real p, real rsw, real rc, shift_consts_t* sc) +{ + /* Here we determine the coefficient for shifting the force to zero + * between distance rsw and the cut-off rc. + * For a potential of r^-p, we have force p*r^-(p+1). + * But to save flops we absorb p in the coefficient. + * Thus we get: + * force/p = r^-(p+1) + c2*r^2 + c3*r^3 + * potential = r^-p + c2/3*r^3 + c3/4*r^4 + cpot + */ + sc->c2 = ((p + 1) * rsw - (p + 4) * rc) / (pow(rc, p + 2) * gmx::square(rc - rsw)); + sc->c3 = -((p + 1) * rsw - (p + 3) * rc) / (pow(rc, p + 2) * gmx::power3(rc - rsw)); + sc->cpot = -pow(rc, -p) + p * sc->c2 / 3 * gmx::power3(rc - rsw) + + p * sc->c3 / 4 * gmx::power4(rc - rsw); +} + +static void potential_switch_constants(real rsw, real rc, switch_consts_t* sc) +{ + /* The switch function is 1 at rsw and 0 at rc. + * The derivative and second derivate are zero at both ends. + * rsw = max(r - r_switch, 0) + * sw = 1 + c3*rsw^3 + c4*rsw^4 + c5*rsw^5 + * dsw = 3*c3*rsw^2 + 4*c4*rsw^3 + 5*c5*rsw^4 + * force = force*dsw - potential*sw + * potential *= sw + */ + sc->c3 = -10 / gmx::power3(rc - rsw); + sc->c4 = 15 / gmx::power4(rc - rsw); + sc->c5 = -6 / gmx::power5(rc - rsw); +} + + +interaction_const_t init_interaction_const(FILE* fp, const t_inputrec& ir, const gmx_mtop_t& mtop, bool systemHasNetCharge) +{ + interaction_const_t interactionConst; + + interactionConst.coulombEwaldTables = std::make_unique(); + interactionConst.vdwEwaldTables = std::make_unique(); + + /* Lennard-Jones */ + interactionConst.vdwtype = ir.vdwtype; + interactionConst.vdw_modifier = ir.vdw_modifier; + interactionConst.reppow = mtop.ffparams.reppow; + interactionConst.rvdw = cutoff_inf(ir.rvdw); + interactionConst.rvdw_switch = ir.rvdw_switch; + interactionConst.ljpme_comb_rule = ir.ljpme_combination_rule; + interactionConst.useBuckingham = (mtop.ffparams.functype[0] == F_BHAM); + if (interactionConst.useBuckingham) + { + interactionConst.buckinghamBMax = calcBuckinghamBMax(fp, mtop); + } + + initVdwEwaldParameters(fp, ir, &interactionConst); + + clear_force_switch_constants(&interactionConst.dispersion_shift); + clear_force_switch_constants(&interactionConst.repulsion_shift); + + switch (interactionConst.vdw_modifier) + { + case InteractionModifiers::PotShift: + /* Only shift the potential, don't touch the force */ + interactionConst.dispersion_shift.cpot = -1.0 / gmx::power6(interactionConst.rvdw); + interactionConst.repulsion_shift.cpot = -1.0 / gmx::power12(interactionConst.rvdw); + break; + case InteractionModifiers::ForceSwitch: + /* Switch the force, switch and shift the potential */ + force_switch_constants( + 6.0, interactionConst.rvdw_switch, interactionConst.rvdw, &interactionConst.dispersion_shift); + force_switch_constants( + 12.0, interactionConst.rvdw_switch, interactionConst.rvdw, &interactionConst.repulsion_shift); + break; + case InteractionModifiers::PotSwitch: + /* Switch the potential and force */ + potential_switch_constants( + interactionConst.rvdw_switch, interactionConst.rvdw, &interactionConst.vdw_switch); + break; + case InteractionModifiers::None: + case InteractionModifiers::ExactCutoff: + /* Nothing to do here */ + break; + default: gmx_incons("unimplemented potential modifier"); + } + + /* Electrostatics */ + interactionConst.eeltype = ir.coulombtype; + interactionConst.coulomb_modifier = ir.coulomb_modifier; + interactionConst.rcoulomb = cutoff_inf(ir.rcoulomb); + interactionConst.rcoulomb_switch = ir.rcoulomb_switch; + interactionConst.epsilon_r = ir.epsilon_r; + + /* Set the Coulomb energy conversion factor */ + if (interactionConst.epsilon_r != 0) + { + interactionConst.epsfac = gmx::c_one4PiEps0 / interactionConst.epsilon_r; + } + else + { + /* eps = 0 is infinite dieletric: no Coulomb interactions */ + interactionConst.epsfac = 0; + } + + /* Reaction-field */ + if (EEL_RF(interactionConst.eeltype)) + { + GMX_RELEASE_ASSERT(interactionConst.eeltype != CoulombInteractionType::GRFNotused, + "GRF is no longer supported"); + interactionConst.reactionFieldPermitivity = ir.epsilon_rf; + calc_rffac(fp, + interactionConst.epsilon_r, + interactionConst.reactionFieldPermitivity, + interactionConst.rcoulomb, + &interactionConst.reactionFieldCoefficient, + &interactionConst.reactionFieldShift); + } + else + { + /* For plain cut-off we might use the reaction-field kernels */ + interactionConst.reactionFieldPermitivity = interactionConst.epsilon_r; + interactionConst.reactionFieldCoefficient = 0; + if (ir.coulomb_modifier == InteractionModifiers::PotShift) + { + interactionConst.reactionFieldShift = 1 / interactionConst.rcoulomb; + } + else + { + interactionConst.reactionFieldShift = 0; + } + } + + initCoulombEwaldParameters(fp, ir, systemHasNetCharge, &interactionConst); + + if (fp != nullptr) + { + real dispersion_shift; + + dispersion_shift = interactionConst.dispersion_shift.cpot; + if (EVDW_PME(interactionConst.vdwtype)) + { + dispersion_shift -= interactionConst.sh_lj_ewald; + } + fprintf(fp, + "Potential shift: LJ r^-12: %.3e r^-6: %.3e", + interactionConst.repulsion_shift.cpot, + dispersion_shift); + + if (interactionConst.eeltype == CoulombInteractionType::Cut) + { + fprintf(fp, ", Coulomb %.e", -interactionConst.reactionFieldShift); + } + else if (EEL_PME(interactionConst.eeltype)) + { + fprintf(fp, ", Ewald %.3e", -interactionConst.sh_ewald); + } + fprintf(fp, "\n"); + } + + if (ir.efep != FreeEnergyPerturbationType::No) + { + GMX_RELEASE_ASSERT(ir.fepvals, "ir.fepvals should be set wth free-energy"); + interactionConst.softCoreParameters = + std::make_unique(*ir.fepvals); + } + + return interactionConst; +} diff --git a/src/gromacs/mdtypes/interaction_const.h b/src/gromacs/mdtypes/interaction_const.h index ce70634b8e..aee0a64ea9 100644 --- a/src/gromacs/mdtypes/interaction_const.h +++ b/src/gromacs/mdtypes/interaction_const.h @@ -36,6 +36,8 @@ #ifndef GMX_MDTYPES_INTERACTION_CONST_H #define GMX_MDTYPES_INTERACTION_CONST_H +#include + #include #include @@ -44,6 +46,8 @@ #include "gromacs/utility/real.h" struct t_lambda; +struct t_inputrec; +struct gmx_mtop_t; /* Used with force switching or a constant potential shift: * rsw = max(r - r_switch, 0) @@ -173,4 +177,15 @@ struct interaction_const_t std::unique_ptr softCoreParameters; }; +/*! \brief Construct interaction constants + * + * This data is used (particularly) by search and force code for + * short-range interactions. Many of these are constant for the whole + * simulation; some are constant only after PME tuning completes. + */ +interaction_const_t init_interaction_const(FILE* fp, + const t_inputrec& ir, + const gmx_mtop_t& mtop, + bool systemHasNetCharge); + #endif -- 2.22.0