src/gromacs/mdtypes/interaction_const.h

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012,2013,2014,2015,2017 by the GROMACS development team.
   5  * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
   6  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   7  * and including many others, as listed in the AUTHORS file in the
   8  * top-level source directory and at http://www.gromacs.org.
   9  *
  10  * GROMACS is free software; you can redistribute it and/or
  11  * modify it under the terms of the GNU Lesser General Public License
  12  * as published by the Free Software Foundation; either version 2.1
  13  * of the License, or (at your option) any later version.
  14  *
  15  * GROMACS is distributed in the hope that it will be useful,
  16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  18  * Lesser General Public License for more details.
  19  *
  20  * You should have received a copy of the GNU Lesser General Public
  21  * License along with GROMACS; if not, see
  22  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  23  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  24  *
  25  * If you want to redistribute modifications to GROMACS, please
  26  * consider that scientific software is very special. Version
  27  * control is crucial - bugs must be traceable. We will be happy to
  28  * consider code for inclusion in the official distribution, but
  29  * derived work must not be called official GROMACS. Details are found
  30  * in the README & COPYING files - if they are missing, get the
  31  * official version at http://www.gromacs.org.
  32  *
  33  * To help us fund GROMACS development, we humbly ask that you cite
  34  * the research papers on the package. Check out http://www.gromacs.org.
  35  */
  36 #ifndef GMX_MDTYPES_INTERACTION_CONST_H
  37 #define GMX_MDTYPES_INTERACTION_CONST_H
  38
  39 #include <cstdio>
  40
  41 #include <memory>
  42 #include <vector>
  43
  44 #include "gromacs/mdtypes/md_enums.h"
  45 #include "gromacs/utility/alignedallocator.h"
  46 #include "gromacs/utility/real.h"
  47
  48 struct t_lambda;
  49 struct t_inputrec;
  50 struct gmx_mtop_t;
  51
  52 /* Used with force switching or a constant potential shift:
  53  * rsw       = max(r - r_switch, 0)
  54  * force/p   = r^-(p+1) + c2*rsw^2 + c3*rsw^3
  55  * potential = r^-p + c2/3*rsw^3 + c3/4*rsw^4 + cpot
  56  * With a constant potential shift c2 and c3 are both 0.
  57  */
  58 struct shift_consts_t
  59 {
  60     real c2;
  61     real c3;
  62     real cpot;
  63 };
  64
  65 /* Used with potential switching:
  66  * rsw        = max(r - r_switch, 0)
  67  * sw         = 1 + c3*rsw^3 + c4*rsw^4 + c5*rsw^5
  68  * dsw        = 3*c3*rsw^2 + 4*c4*rsw^3 + 5*c5*rsw^4
  69  * force      = force*dsw - potential*sw
  70  * potential *= sw
  71  */
  72 struct switch_consts_t
  73 {
  74     real c3;
  75     real c4;
  76     real c5;
  77 };
  78
  79 /* Convenience type for vector with aligned memory */
  80 template<typename T>
  81 using AlignedVector = std::vector<T, gmx::AlignedAllocator<T>>;
  82
  83 /* Force/energy interpolation tables for Ewald long-range corrections
  84  *
  85  * Interpolation is linear for the force, quadratic for the potential.
  86  */
  87 struct EwaldCorrectionTables
  88 {
  89     // 1/table_spacing, units 1/nm
  90     real scale = 0;
  91     // Force table
  92     AlignedVector<real> tableF;
  93     // Energy table
  94     AlignedVector<real> tableV;
  95     // Coulomb force+energy table, size of array is tabq_size*4,
  96     // entry quadruplets are: F[i], F[i+1]-F[i], V[i], 0,
  97     // this is used with 4-wide SIMD for aligned loads
  98     AlignedVector<real> tableFDV0;
  99 };
 100
 101 /* The physical interaction parameters for non-bonded interaction calculations
 102  *
 103  * This struct contains copies of the physical interaction parameters
 104  * from the user input as well as processed values that are need in
 105  * non-bonded interaction kernels.
 106  *
 107  * The default constructor gives plain Coulomb and LJ interactions cut off
 108  * a 1 nm without potential shifting and a Coulomb pre-factor of 1.
 109  */
 110 struct interaction_const_t
 111 {
 112     /* This struct contains the soft-core parameters from t_lambda,
 113      * but processed for direct use in the kernels.
 114      */
 115     struct SoftCoreParameters
 116     {
 117         // Constructor
 118         SoftCoreParameters(const t_lambda& fepvals);
 119
 120         // Alpha parameter for Van der Waals interactions
 121         real alphaVdw;
 122         // Alpha parameter for Coulomb interactions
 123         real alphaCoulomb;
 124         // Exponent for the dependence of the soft-core on lambda
 125         int lambdaPower;
 126         // Value for sigma^6 for LJ interaction with C6<=0 and/or C12<=0
 127         real sigma6WithInvalidSigma;
 128         // Minimum value for sigma^6, used when soft-core is applied to Coulomb interactions
 129         real sigma6Minimum;
 130     };
 131
 132     /* VdW */
 133     VanDerWaalsType        vdwtype          = VanDerWaalsType::Cut;
 134     InteractionModifiers   vdw_modifier     = InteractionModifiers::None;
 135     double                 reppow           = 12;
 136     real                   rvdw             = 1;
 137     real                   rvdw_switch      = 0;
 138     struct shift_consts_t  dispersion_shift = { 0, 0, 0 };
 139     struct shift_consts_t  repulsion_shift  = { 0, 0, 0 };
 140     struct switch_consts_t vdw_switch       = { 0, 0, 0 };
 141     bool                   useBuckingham    = false;
 142     real                   buckinghamBMax   = 0;
 143
 144     /* type of electrostatics */
 145     CoulombInteractionType eeltype          = CoulombInteractionType::Cut;
 146     InteractionModifiers   coulomb_modifier = InteractionModifiers::None;
 147
 148     /* Coulomb */
 149     real rcoulomb        = 1;
 150     real rcoulomb_switch = 0;
 151
 152     /* PME/Ewald */
 153     real         ewaldcoeff_q  = 0;
 154     real         ewaldcoeff_lj = 0;
 155     LongRangeVdW ljpme_comb_rule = LongRangeVdW::Geom; /* LJ combination rule for the LJ PME mesh part */
 156     real         sh_ewald        = 0; /* -sh_ewald is added to the direct space potential */
 157     real         sh_lj_ewald = 0;     /* sh_lj_ewald is added to the correction potential */
 158
 159     /* Dielectric constant resp. multiplication factor for charges */
 160     real epsilon_r = 1;
 161     real epsfac    = 1;
 162
 163     /* Constants for reaction-field or plain cut-off */
 164     //! Dielectric constant for reaction field beyond the cutoff distance
 165     real reactionFieldPermitivity = 1;
 166     //! Coefficient for reaction field; scales relation between epsilon_r and reactionFieldPermitivity
 167     real reactionFieldCoefficient = 0;
 168     //! Constant shift to reaction field Coulomb interaction to make potential an integral of force
 169     real reactionFieldShift = 0;
 170
 171     // Coulomb Ewald correction table
 172     std::unique_ptr<EwaldCorrectionTables> coulombEwaldTables;
 173     // Van der Waals Ewald correction table
 174     std::unique_ptr<EwaldCorrectionTables> vdwEwaldTables;
 175
 176     // Free-energy parameters, only present when free-energy calculations are requested
 177     std::unique_ptr<SoftCoreParameters> softCoreParameters;
 178 };
 179
 180 /*! \brief Construct interaction constants
 181  *
 182  * This data is used (particularly) by search and force code for
 183  * short-range interactions. Many of these are constant for the whole
 184  * simulation; some are constant only after PME tuning completes.
 185  */
 186 interaction_const_t init_interaction_const(FILE*             fp,
 187                                            const t_inputrec& ir,
 188                                            const gmx_mtop_t& mtop,
 189                                            bool              systemHasNetCharge);
 190
 191 #endif