src/gromacs/gmxlib/nonbonded/nb_softcore.h

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012,2013,2014,2016,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_GMXLIB_NONBONDED_SOFTCORE_H
  37 #define GMX_GMXLIB_NONBONDED_SOFTCORE_H
  38
  39 #include "config.h"
  40
  41 #include "gromacs/math/functions.h"
  42 #include "gromacs/simd/simd.h"
  43 #include "gromacs/simd/simd_math.h"
  44
  45 /* linearized electrostatics */
  46 template<bool computeForces, class RealType, class BoolType>
  47 static inline void quadraticApproximationCoulomb(const RealType qq,
  48                                                  const RealType rInvQ,
  49                                                  const RealType r,
  50                                                  const real     lambdaFac,
  51                                                  const real     dLambdaFac,
  52                                                  RealType gmx_unused* force,
  53                                                  RealType*            potential,
  54                                                  RealType*            dvdl,
  55                                                  BoolType             dvdlMask)
  56 {
  57     RealType constFac = qq * rInvQ;
  58     RealType linFac   = constFac * r * rInvQ;
  59     RealType quadrFac = linFac * r * rInvQ;
  60
  61     /* Computing Coulomb force and potential energy */
  62     if constexpr (computeForces)
  63     {
  64         *force = -2 * quadrFac + 3 * linFac;
  65     }
  66
  67     *potential = quadrFac - 3 * (linFac - constFac);
  68
  69     RealType lambdaFacRevInv = gmx::maskzInv(1 - lambdaFac, dvdlMask);
  70     *dvdl = dLambdaFac * 0.5_real * (lambdaFac * lambdaFacRevInv) * (quadrFac - 2 * linFac + constFac);
  71 }
  72
  73 /* reaction-field linearized electrostatics */
  74 template<bool computeForces, class RealType, class BoolType>
  75 static inline void reactionFieldQuadraticPotential(const RealType qq,
  76                                                    const real     facel,
  77                                                    const RealType r,
  78                                                    const real     rCutoff,
  79                                                    const real     lambdaFac,
  80                                                    const real     dLambdaFac,
  81                                                    const RealType alphaEff,
  82                                                    const real     krf,
  83                                                    const real     potentialShift,
  84                                                    RealType gmx_unused* force,
  85                                                    RealType*            potential,
  86                                                    RealType*            dvdl,
  87                                                    BoolType             mask)
  88 {
  89     /* check if we have to use the hardcore values */
  90     BoolType computeValues = mask && (lambdaFac < 1 && 0 < alphaEff && facel != 0);
  91     if (gmx::anyTrue(computeValues))
  92     {
  93         RealType lambdaFacRev = gmx::selectByMask(1 - lambdaFac, computeValues);
  94
  95         RealType rQ = gmx::cbrt(lambdaFacRev);
  96         rQ          = gmx::sqrt(rQ) * (1 + gmx::abs(qq / facel));
  97         rQ          = rQ * alphaEff;
  98
  99         // ensure that the linearization point doesn't go beyond rCutoff
 100         BoolType beyondCutoff = rCutoff < rQ;
 101         BoolType withinCutoff = rQ <= rCutoff;
 102         if (gmx::anyTrue(beyondCutoff))
 103         {
 104             rQ = gmx::blend(rQ, rCutoff, beyondCutoff);
 105         }
 106
 107         computeValues = computeValues && (r < rQ);
 108         if (gmx::anyTrue(computeValues))
 109         {
 110             RealType rInvQ = gmx::maskzInv(rQ, computeValues);
 111
 112             // Compute quadratic force, potential and dvdl
 113             RealType forceQuad(0);
 114             RealType potentialQuad(0);
 115             RealType dvdlQuad(0);
 116             quadraticApproximationCoulomb<computeForces>(
 117                     qq, rInvQ, r, lambdaFac, dLambdaFac, &forceQuad, &potentialQuad, &dvdlQuad, computeValues);
 118
 119             // rf modification
 120             forceQuad     = forceQuad - qq * 2 * krf * r * r;
 121             potentialQuad = potentialQuad + qq * (krf * r * r - potentialShift);
 122
 123             // update
 124             if constexpr (computeForces)
 125             {
 126                 *force = gmx::blend(*force, forceQuad, computeValues);
 127             }
 128             *potential = gmx::blend(*potential, potentialQuad, computeValues);
 129             *dvdl      = *dvdl + gmx::selectByMask(dvdlQuad, computeValues && withinCutoff);
 130         }
 131     }
 132 }
 133
 134 /* ewald linearized electrostatics */
 135 template<bool computeForces, class RealType, class BoolType>
 136 static inline void ewaldQuadraticPotential(const RealType qq,
 137                                            const real     facel,
 138                                            const RealType r,
 139                                            const real     rCutoff,
 140                                            const real     lambdaFac,
 141                                            const real     dLambdaFac,
 142                                            const RealType alphaEff,
 143                                            const real     potentialShift,
 144                                            RealType gmx_unused* force,
 145                                            RealType*            potential,
 146                                            RealType*            dvdl,
 147                                            BoolType             mask)
 148 {
 149     /* check if we have to use the hardcore values */
 150     BoolType computeValues = mask && (lambdaFac < 1 && 0 < alphaEff && facel != 0);
 151     if (gmx::anyTrue(computeValues))
 152     {
 153         RealType lambdaFacRev = gmx::selectByMask(1 - lambdaFac, computeValues);
 154
 155         RealType rQ = gmx::cbrt(lambdaFacRev);
 156         rQ          = gmx::sqrt(rQ) * (1 + gmx::abs(qq / facel));
 157         rQ          = rQ * alphaEff;
 158
 159         // ensure that the linearization point doesn't go beyond rCutoff
 160         BoolType beyondCutoff = rCutoff < rQ;
 161         BoolType withinCutoff = rQ <= rCutoff;
 162         if (gmx::anyTrue(beyondCutoff))
 163         {
 164             rQ = gmx::blend(rQ, rCutoff, beyondCutoff);
 165         }
 166
 167         computeValues = computeValues && (r < rQ);
 168         if (gmx::anyTrue(computeValues))
 169         {
 170             RealType rInvQ = gmx::maskzInv(rQ, computeValues);
 171
 172             // Compute quadratic force, potential and dvdl
 173             RealType forceQuad(0);
 174             RealType potentialQuad(0);
 175             RealType dvdlQuad(0);
 176             quadraticApproximationCoulomb<computeForces>(
 177                     qq, rInvQ, r, lambdaFac, dLambdaFac, &forceQuad, &potentialQuad, &dvdlQuad, computeValues);
 178
 179             // ewald modification
 180             potentialQuad = potentialQuad - qq * potentialShift;
 181
 182             // update
 183             if constexpr (computeForces)
 184             {
 185                 *force = gmx::blend(*force, forceQuad, computeValues);
 186             }
 187             *potential = gmx::blend(*potential, potentialQuad, computeValues);
 188             *dvdl      = *dvdl + gmx::selectByMask(dvdlQuad, computeValues && withinCutoff);
 189         }
 190     }
 191 }
 192
 193 /* cutoff LJ with quadratic appximation of lj-potential */
 194 template<bool computeForces, class RealType, class BoolType>
 195 static inline void lennardJonesQuadraticPotential(const RealType c6,
 196                                                   const RealType c12,
 197                                                   const RealType r,
 198                                                   const RealType rsq,
 199                                                   const real     lambdaFac,
 200                                                   const real     dLambdaFac,
 201                                                   const RealType sigma6,
 202                                                   const RealType alphaEff,
 203                                                   const real     repulsionShift,
 204                                                   const real     dispersionShift,
 205                                                   RealType gmx_unused* force,
 206                                                   RealType*            potential,
 207                                                   RealType*            dvdl,
 208                                                   BoolType             mask)
 209 {
 210     constexpr real c_twentySixSeventh = 26.0_real / 7.0_real;
 211     constexpr real c_oneSixth         = 1.0_real / 6.0_real;
 212     constexpr real c_oneTwelth        = 1.0_real / 12.0_real;
 213     constexpr real c_half             = 1.0_real / 2.0_real;
 214
 215     /* check if we have to use the hardcore values */
 216     BoolType computeValues = mask && (lambdaFac < 1 && 0 < alphaEff);
 217     if (gmx::anyTrue(computeValues))
 218     {
 219         RealType lambdaFacRev    = gmx::selectByMask(1 - lambdaFac, computeValues);
 220         RealType lambdaFacRevInv = gmx::maskzInv(1 - lambdaFac, computeValues);
 221
 222         RealType rQ = gmx::cbrt(c_twentySixSeventh * sigma6 * lambdaFacRev);
 223         rQ          = gmx::sqrt(rQ);
 224         rQ          = rQ * alphaEff;
 225
 226         computeValues = (computeValues && r < rQ);
 227         if (gmx::anyTrue(computeValues))
 228         {
 229             /* scaled values for c6 and c12 */
 230             RealType c6s, c12s;
 231             c6s  = c_oneSixth * c6;
 232             c12s = c_oneTwelth * c12;
 233             /* Temporary variables for inverted values */
 234             RealType rInvQ = gmx::maskzInv(rQ, computeValues);
 235             RealType rInv14C, rInv13C, rInv12C;
 236             RealType rInv8C, rInv7C, rInv6C;
 237             rInv6C  = rInvQ * rInvQ * rInvQ;
 238             rInv6C  = rInv6C * rInv6C;
 239             rInv7C  = rInv6C * rInvQ;
 240             rInv8C  = rInv7C * rInvQ;
 241             rInv14C = c12s * rInv7C * rInv7C * rsq;
 242             rInv13C = c12s * rInv7C * rInv6C * r;
 243             rInv12C = c12s * rInv6C * rInv6C;
 244             rInv8C  = rInv8C * c6s * rsq;
 245             rInv7C  = rInv7C * c6s * r;
 246             rInv6C  = rInv6C * c6s;
 247
 248             /* Temporary variables for A and B */
 249             RealType quadrFac, linearFac, constFac;
 250             quadrFac  = 156 * rInv14C - 42 * rInv8C;
 251             linearFac = 168 * rInv13C - 48 * rInv7C;
 252             constFac  = 91 * rInv12C - 28 * rInv6C;
 253
 254             /* Computing LJ force and potential energy */
 255             RealType gmx_unused forceQuad     = -quadrFac + linearFac;
 256             RealType            potentialQuad = c_half * quadrFac - linearFac + constFac;
 257             RealType            dvdlQuad      = dLambdaFac * 28 * (lambdaFac * lambdaFacRevInv)
 258                                 * ((6.5_real * rInv14C - rInv8C) - (13 * rInv13C - 2 * rInv7C)
 259                                    + (6.5_real * rInv12C - rInv6C));
 260
 261             potentialQuad = potentialQuad
 262                             + gmx::selectByMask(((c12s * repulsionShift) - (c6s * dispersionShift)),
 263                                                 computeValues);
 264             if constexpr (computeForces)
 265             {
 266                 *force = gmx::blend(*force, forceQuad, computeValues);
 267             }
 268             *potential = gmx::blend(*potential, potentialQuad, computeValues);
 269             *dvdl      = *dvdl + gmx::selectByMask(dvdlQuad, computeValues);
 270         }
 271     }
 272 }
 273
 274 #endif