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