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