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