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