2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
36 #ifndef GMX_GMXLIB_NONBONDED_SOFTCORE_H
37 #define GMX_GMXLIB_NONBONDED_SOFTCORE_H
41 #include "gromacs/math/functions.h"
42 #include "gromacs/simd/simd.h"
43 #include "gromacs/simd/simd_math.h"
45 /* linearized electrostatics */
46 template<class RealType, class BoolType>
47 static inline void quadraticApproximationCoulomb(const RealType qq,
51 const real dLambdaFac,
57 RealType constFac = qq * rInvQ;
58 RealType linFac = constFac * r * rInvQ;
59 RealType quadrFac = linFac * r * rInvQ;
61 /* Computing Coulomb force and potential energy */
62 *force = -2 * quadrFac + 3 * linFac;
64 *potential = quadrFac - 3 * (linFac - constFac);
66 RealType lambdaFacRevInv = gmx::maskzInv(1 - lambdaFac, dvdlMask);
67 *dvdl = dLambdaFac * 0.5_real * (lambdaFac * lambdaFacRevInv) * (quadrFac - 2 * linFac + constFac);
70 /* reaction-field linearized electrostatics */
71 template<bool computeForces, class RealType, class BoolType>
72 static inline void reactionFieldQuadraticPotential(const RealType qq,
77 const real dLambdaFac,
78 const RealType alphaEff,
80 const real potentialShift,
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))
90 RealType lambdaFacRev = gmx::selectByMask(1 - lambdaFac, computeValues);
92 RealType rQ = gmx::cbrt(lambdaFacRev);
93 rQ = gmx::sqrt(rQ) * (1 + gmx::abs(qq / facel));
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))
101 rQ = gmx::blend(rQ, rCutoff, beyondCutoff);
104 computeValues = computeValues && (r < rQ);
105 if (gmx::anyTrue(computeValues))
107 RealType rInvQ = gmx::maskzInv(rQ, computeValues);
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);
117 forceQuad = forceQuad - qq * 2 * krf * r * r;
118 potentialQuad = potentialQuad + qq * (krf * r * r - potentialShift);
121 if constexpr (computeForces)
123 *force = gmx::blend(*force, forceQuad, computeValues);
125 *potential = gmx::blend(*potential, potentialQuad, computeValues);
126 *dvdl = *dvdl + gmx::selectByMask(dvdlQuad, computeValues && withinCutoff);
131 /* ewald linearized electrostatics */
132 template<bool computeForces, class RealType, class BoolType>
133 static inline void ewaldQuadraticPotential(const RealType qq,
137 const real lambdaFac,
138 const real dLambdaFac,
139 const RealType alphaEff,
140 const real potentialShift,
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))
150 RealType lambdaFacRev = gmx::selectByMask(1 - lambdaFac, computeValues);
152 RealType rQ = gmx::cbrt(lambdaFacRev);
153 rQ = gmx::sqrt(rQ) * (1 + gmx::abs(qq / facel));
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))
161 rQ = gmx::blend(rQ, rCutoff, beyondCutoff);
164 computeValues = computeValues && (r < rQ);
165 if (gmx::anyTrue(computeValues))
167 RealType rInvQ = gmx::maskzInv(rQ, computeValues);
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);
176 // ewald modification
177 potentialQuad = potentialQuad - qq * potentialShift;
180 if constexpr (computeForces)
182 *force = gmx::blend(*force, forceQuad, computeValues);
184 *potential = gmx::blend(*potential, potentialQuad, computeValues);
185 *dvdl = *dvdl + gmx::selectByMask(dvdlQuad, computeValues && withinCutoff);
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,
196 const real lambdaFac,
197 const real dLambdaFac,
198 const RealType sigma6,
199 const RealType alphaEff,
200 const real repulsionShift,
201 const real dispersionShift,
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;
212 /* check if we have to use the hardcore values */
213 BoolType computeValues = mask && (lambdaFac < 1 && 0 < alphaEff);
214 if (gmx::anyTrue(computeValues))
216 RealType lambdaFacRev = gmx::selectByMask(1 - lambdaFac, computeValues);
217 RealType lambdaFacRevInv = gmx::maskzInv(1 - lambdaFac, computeValues);
219 RealType rQ = gmx::cbrt(c_twentySixSeventh * sigma6 * lambdaFacRev);
223 computeValues = (computeValues && r < rQ);
224 if (gmx::anyTrue(computeValues))
226 /* scaled values for c6 and c12 */
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;
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;
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));
258 potentialQuad = potentialQuad
259 + gmx::selectByMask(((c12s * repulsionShift) - (c6s * dispersionShift)),
261 if constexpr (computeForces)
263 *force = gmx::blend(*force, forceQuad, computeValues);
265 *potential = gmx::blend(*potential, potentialQuad, computeValues);
266 *dvdl = *dvdl + gmx::selectByMask(dvdlQuad, computeValues);