2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2012,2013,2014, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
35 #ifndef GMX_SIMD_MATH_DOUBLE_H_
36 #define GMX_SIMD_MATH_DOUBLE_H_
40 static gmx_inline gmx_simd_real_t
41 gmx_simd_invsqrt_r(gmx_simd_real_t x)
43 const gmx_simd_real_t half = gmx_simd_set1_r(0.5);
44 const gmx_simd_real_t three = gmx_simd_set1_r(3.0);
46 /* Lookup instruction only exists in single precision, convert back and forth... */
47 gmx_simd_real_t lu = gmx_simd_rsqrt_r(x);
49 lu = gmx_simd_mul_r(gmx_simd_mul_r(half, lu), gmx_simd_fnmadd_r(gmx_simd_mul_r(lu, lu), x, three));
50 return gmx_simd_mul_r(gmx_simd_mul_r(half, lu), gmx_simd_fnmadd_r(gmx_simd_mul_r(lu, lu), x, three));
55 static gmx_inline gmx_simd_real_t
56 gmx_simd_inv_r(gmx_simd_real_t x)
58 const gmx_simd_real_t two = gmx_simd_set1_r(2.0);
60 /* Lookup instruction only exists in single precision, convert back and forth... */
61 gmx_simd_real_t lu = gmx_simd_rcp_r(x);
63 /* Perform two N-R steps for double precision */
64 lu = gmx_simd_mul_r(lu, gmx_simd_fnmadd_r(lu, x, two));
65 return gmx_simd_mul_r(lu, gmx_simd_fnmadd_r(lu, x, two));
69 /* Calculate the force correction due to PME analytically.
71 * This routine is meant to enable analytical evaluation of the
72 * direct-space PME electrostatic force to avoid tables.
74 * The direct-space potential should be Erfc(beta*r)/r, but there
75 * are some problems evaluating that:
77 * First, the error function is difficult (read: expensive) to
78 * approxmiate accurately for intermediate to large arguments, and
79 * this happens already in ranges of beta*r that occur in simulations.
80 * Second, we now try to avoid calculating potentials in Gromacs but
81 * use forces directly.
83 * We can simply things slight by noting that the PME part is really
84 * a correction to the normal Coulomb force since Erfc(z)=1-Erf(z), i.e.
86 * V= 1/r - Erf(beta*r)/r
88 * The first term we already have from the inverse square root, so
89 * that we can leave out of this routine.
91 * For pme tolerances of 1e-3 to 1e-8 and cutoffs of 0.5nm to 1.8nm,
92 * the argument beta*r will be in the range 0.15 to ~4. Use your
93 * favorite plotting program to realize how well-behaved Erf(z)/z is
96 * We approximate f(z)=erf(z)/z with a rational minimax polynomial.
97 * However, it turns out it is more efficient to approximate f(z)/z and
98 * then only use even powers. This is another minor optimization, since
99 * we actually WANT f(z)/z, because it is going to be multiplied by
100 * the vector between the two atoms to get the vectorial force. The
101 * fastest flops are the ones we can avoid calculating!
103 * So, here's how it should be used:
106 * 2. Multiply by beta^2, so you get z^2=beta^2*r^2.
107 * 3. Evaluate this routine with z^2 as the argument.
108 * 4. The return value is the expression:
112 * ------------ - --------
115 * 5. Multiply the entire expression by beta^3. This will get you
117 * beta^3*2*exp(-z^2) beta^3*erf(z)
118 * ------------------ - ---------------
121 * or, switching back to r (z=r*beta):
123 * 2*beta*exp(-r^2*beta^2) erf(r*beta)
124 * ----------------------- - -----------
128 * With a bit of math exercise you should be able to confirm that
129 * this is exactly D[Erf[beta*r]/r,r] divided by r another time.
131 * 6. Add the result to 1/r^3, multiply by the product of the charges,
132 * and you have your force (divided by r). A final multiplication
133 * with the vector connecting the two particles and you have your
134 * vectorial force to add to the particles.
137 static gmx_simd_real_t
138 gmx_simd_pmecorrF_r(gmx_simd_real_t z2)
140 const gmx_simd_real_t FN10 = gmx_simd_set1_r(-8.0072854618360083154e-14);
141 const gmx_simd_real_t FN9 = gmx_simd_set1_r(1.1859116242260148027e-11);
142 const gmx_simd_real_t FN8 = gmx_simd_set1_r(-8.1490406329798423616e-10);
143 const gmx_simd_real_t FN7 = gmx_simd_set1_r(3.4404793543907847655e-8);
144 const gmx_simd_real_t FN6 = gmx_simd_set1_r(-9.9471420832602741006e-7);
145 const gmx_simd_real_t FN5 = gmx_simd_set1_r(0.000020740315999115847456);
146 const gmx_simd_real_t FN4 = gmx_simd_set1_r(-0.00031991745139313364005);
147 const gmx_simd_real_t FN3 = gmx_simd_set1_r(0.0035074449373659008203);
148 const gmx_simd_real_t FN2 = gmx_simd_set1_r(-0.031750380176100813405);
149 const gmx_simd_real_t FN1 = gmx_simd_set1_r(0.13884101728898463426);
150 const gmx_simd_real_t FN0 = gmx_simd_set1_r(-0.75225277815249618847);
152 const gmx_simd_real_t FD5 = gmx_simd_set1_r(0.000016009278224355026701);
153 const gmx_simd_real_t FD4 = gmx_simd_set1_r(0.00051055686934806966046);
154 const gmx_simd_real_t FD3 = gmx_simd_set1_r(0.0081803507497974289008);
155 const gmx_simd_real_t FD2 = gmx_simd_set1_r(0.077181146026670287235);
156 const gmx_simd_real_t FD1 = gmx_simd_set1_r(0.41543303143712535988);
157 const gmx_simd_real_t FD0 = gmx_simd_set1_r(1.0);
160 gmx_simd_real_t polyFN0, polyFN1, polyFD0, polyFD1;
162 z4 = gmx_simd_mul_r(z2, z2);
164 polyFD1 = gmx_simd_fmadd_r(FD5, z4, FD3);
165 polyFD1 = gmx_simd_fmadd_r(polyFD1, z4, FD1);
166 polyFD1 = gmx_simd_mul_r(polyFD1, z2);
167 polyFD0 = gmx_simd_fmadd_r(FD4, z4, FD2);
168 polyFD0 = gmx_simd_fmadd_r(polyFD0, z4, FD0);
169 polyFD0 = gmx_simd_add_r(polyFD0, polyFD1);
171 polyFD0 = gmx_simd_inv_r(polyFD0);
173 polyFN0 = gmx_simd_fmadd_r(FN10, z4, FN8);
174 polyFN0 = gmx_simd_fmadd_r(polyFN0, z4, FN6);
175 polyFN0 = gmx_simd_fmadd_r(polyFN0, z4, FN4);
176 polyFN0 = gmx_simd_fmadd_r(polyFN0, z4, FN2);
177 polyFN0 = gmx_simd_fmadd_r(polyFN0, z4, FN0);
178 polyFN1 = gmx_simd_fmadd_r(FN9, z4, FN7);
179 polyFN1 = gmx_simd_fmadd_r(polyFN1, z4, FN5);
180 polyFN1 = gmx_simd_fmadd_r(polyFN1, z4, FN3);
181 polyFN1 = gmx_simd_fmadd_r(polyFN1, z4, FN1);
182 polyFN0 = gmx_simd_fmadd_r(polyFN1, z2, polyFN0);
184 return gmx_simd_mul_r(polyFN0, polyFD0);
188 /* Calculate the potential correction due to PME analytically.
190 * This routine calculates Erf(z)/z, although you should provide z^2
191 * as the input argument.
193 * Here's how it should be used:
196 * 2. Multiply by beta^2, so you get z^2=beta^2*r^2.
197 * 3. Evaluate this routine with z^2 as the argument.
198 * 4. The return value is the expression:
205 * 5. Multiply the entire expression by beta and switching back to r (z=r*beta):
211 * 6. Subtract the result from 1/r, multiply by the product of the charges,
212 * and you have your potential.
215 static gmx_simd_real_t
216 gmx_simd_pmecorrV_r(gmx_simd_real_t z2)
218 const gmx_simd_real_t VN9 = gmx_simd_set1_r(-9.3723776169321855475e-13);
219 const gmx_simd_real_t VN8 = gmx_simd_set1_r(1.2280156762674215741e-10);
220 const gmx_simd_real_t VN7 = gmx_simd_set1_r(-7.3562157912251309487e-9);
221 const gmx_simd_real_t VN6 = gmx_simd_set1_r(2.6215886208032517509e-7);
222 const gmx_simd_real_t VN5 = gmx_simd_set1_r(-4.9532491651265819499e-6);
223 const gmx_simd_real_t VN4 = gmx_simd_set1_r(0.00025907400778966060389);
224 const gmx_simd_real_t VN3 = gmx_simd_set1_r(0.0010585044856156469792);
225 const gmx_simd_real_t VN2 = gmx_simd_set1_r(0.045247661136833092885);
226 const gmx_simd_real_t VN1 = gmx_simd_set1_r(0.11643931522926034421);
227 const gmx_simd_real_t VN0 = gmx_simd_set1_r(1.1283791671726767970);
229 const gmx_simd_real_t VD5 = gmx_simd_set1_r(0.000021784709867336150342);
230 const gmx_simd_real_t VD4 = gmx_simd_set1_r(0.00064293662010911388448);
231 const gmx_simd_real_t VD3 = gmx_simd_set1_r(0.0096311444822588683504);
232 const gmx_simd_real_t VD2 = gmx_simd_set1_r(0.085608012351550627051);
233 const gmx_simd_real_t VD1 = gmx_simd_set1_r(0.43652499166614811084);
234 const gmx_simd_real_t VD0 = gmx_simd_set1_r(1.0);
237 gmx_simd_real_t polyVN0, polyVN1, polyVD0, polyVD1;
239 z4 = gmx_simd_mul_r(z2, z2);
241 polyVD1 = gmx_simd_fmadd_r(VD5, z4, VD3);
242 polyVD0 = gmx_simd_fmadd_r(VD4, z4, VD2);
243 polyVD1 = gmx_simd_fmadd_r(polyVD1, z4, VD1);
244 polyVD0 = gmx_simd_fmadd_r(polyVD0, z4, VD0);
245 polyVD0 = gmx_simd_fmadd_r(polyVD1, z2, polyVD0);
247 polyVD0 = gmx_simd_inv_r(polyVD0);
249 polyVN1 = gmx_simd_fmadd_r(VN9, z4, VN7);
250 polyVN0 = gmx_simd_fmadd_r(VN8, z4, VN6);
251 polyVN1 = gmx_simd_fmadd_r(polyVN1, z4, VN5);
252 polyVN0 = gmx_simd_fmadd_r(polyVN0, z4, VN4);
253 polyVN1 = gmx_simd_fmadd_r(polyVN1, z4, VN3);
254 polyVN0 = gmx_simd_fmadd_r(polyVN0, z4, VN2);
255 polyVN1 = gmx_simd_fmadd_r(polyVN1, z4, VN1);
256 polyVN0 = gmx_simd_fmadd_r(polyVN0, z4, VN0);
257 polyVN0 = gmx_simd_fmadd_r(polyVN1, z2, polyVN0);
259 return gmx_simd_mul_r(polyVN0, polyVD0);