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<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);
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 *force = gmx::blend(*force, forceQuad, computeValues);
122 *potential = gmx::blend(*potential, potentialQuad, computeValues);
123 *dvdl = *dvdl + gmx::selectByMask(dvdlQuad, computeValues && withinCutoff);
128 /* ewald linearized electrostatics */
129 template<class RealType, class BoolType>
130 static inline void ewaldQuadraticPotential(const RealType qq,
134 const real lambdaFac,
135 const real dLambdaFac,
136 const RealType alphaEff,
137 const real potentialShift,
144 /* check if we have to use the hardcore values */
145 BoolType computeValues = mask && (lambdaFac < 1 && 0 < alphaEff);
146 if (gmx::anyTrue(computeValues))
148 RealType lambdaFacRev = gmx::selectByMask(1 - lambdaFac, computeValues);
150 RealType rQ = gmx::cbrt(lambdaFacRev);
151 rQ = gmx::sqrt(rQ) * (1 + gmx::abs(qq / facel));
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))
159 rQ = gmx::blend(rQ, rCutoff, beyondCutoff);
162 computeValues = computeValues && (r < rQ);
163 if (gmx::anyTrue(computeValues))
165 RealType rInvQ = gmx::maskzInv(rQ, computeValues);
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);
174 // ewald modification
175 potentialQuad = potentialQuad - qq * potentialShift;
178 *force = gmx::blend(*force, forceQuad, computeValues);
179 *potential = gmx::blend(*potential, potentialQuad, computeValues);
180 *dvdl = *dvdl + gmx::selectByMask(dvdlQuad, computeValues && withinCutoff);
185 /* cutoff LJ with quadratic appximation of lj-potential */
186 template<class RealType, class BoolType>
187 static inline void lennardJonesQuadraticPotential(const RealType c6,
191 const real lambdaFac,
192 const real dLambdaFac,
193 const RealType sigma6,
194 const RealType alphaEff,
195 const real repulsionShift,
196 const real dispersionShift,
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;
207 /* check if we have to use the hardcore values */
208 BoolType computeValues = mask && (lambdaFac < 1 && 0 < alphaEff);
209 if (gmx::anyTrue(computeValues))
211 RealType lambdaFacRev = gmx::selectByMask(1 - lambdaFac, computeValues);
212 RealType lambdaFacRevInv = gmx::maskzInv(1 - lambdaFac, computeValues);
214 RealType rQ = gmx::cbrt(c_twentySixSeventh * sigma6 * lambdaFacRev);
218 computeValues = (computeValues && r < rQ);
219 if (gmx::anyTrue(computeValues))
221 /* scaled values for c6 and c12 */
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;
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;
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));
253 potentialQuad = potentialQuad
254 + gmx::selectByMask(((c12s * repulsionShift) - (c6s * dispersionShift)),
256 *force = gmx::blend(*force, forceQuad, computeValues);
257 *potential = gmx::blend(*potential, potentialQuad, computeValues);
258 *dvdl = *dvdl + gmx::selectByMask(dvdlQuad, computeValues);