include/gmx_simd_math_double.h

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012,2013, by the GROMACS development team, led by
   5  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   6  * others, as listed in the AUTHORS file in the top-level source
   7  * directory and at http://www.gromacs.org.
   8  *
   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.
  13  *
  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.
  18  *
  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.
  23  *
  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.
  31  *
  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.
  34  */
  35 #ifndef _gmx_simd_math_double_h_
  36 #define _gmx_simd_math_double_h_
  37
  38
  39 /* 1.0/sqrt(x) */
  40 static gmx_inline gmx_mm_pr
  41 gmx_invsqrt_pr(gmx_mm_pr x)
  42 {
  43     const gmx_mm_pr half  = gmx_set1_pr(0.5);
  44     const gmx_mm_pr three = gmx_set1_pr(3.0);
  45
  46     /* Lookup instruction only exists in single precision, convert back and forth... */
  47     gmx_mm_pr lu = gmx_rsqrt_pr(x);
  48
  49     lu = gmx_mul_pr(gmx_mul_pr(half, lu), gmx_nmsub_pr(gmx_mul_pr(lu, lu), x, three));
  50     return gmx_mul_pr(gmx_mul_pr(half, lu), gmx_nmsub_pr(gmx_mul_pr(lu, lu), x, three));
  51 }
  52
  53
  54 /* 1.0/x */
  55 static gmx_inline gmx_mm_pr
  56 gmx_inv_pr(gmx_mm_pr x)
  57 {
  58     const gmx_mm_pr two  = gmx_set1_pr(2.0);
  59
  60     /* Lookup instruction only exists in single precision, convert back and forth... */
  61     gmx_mm_pr lu = gmx_rcp_pr(x);
  62
  63     /* Perform two N-R steps for double precision */
  64     lu         = gmx_mul_pr(lu, gmx_nmsub_pr(lu, x, two));
  65     return gmx_mul_pr(lu, gmx_nmsub_pr(lu, x, two));
  66 }
  67
  68
  69 /* Calculate the force correction due to PME analytically.
  70  *
  71  * This routine is meant to enable analytical evaluation of the
  72  * direct-space PME electrostatic force to avoid tables.
  73  *
  74  * The direct-space potential should be Erfc(beta*r)/r, but there
  75  * are some problems evaluating that:
  76  *
  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.
  82  *
  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.
  85  *
  86  * V= 1/r - Erf(beta*r)/r
  87  *
  88  * The first term we already have from the inverse square root, so
  89  * that we can leave out of this routine.
  90  *
  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
  94  * in this range!
  95  *
  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!
 102  *
 103  * So, here's how it should be used:
 104  *
 105  * 1. Calculate r^2.
 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:
 109  *
 110  *
 111  *       2*exp(-z^2)     erf(z)
 112  *       ------------ - --------
 113  *       sqrt(Pi)*z^2      z^3
 114  *
 115  * 5. Multiply the entire expression by beta^3. This will get you
 116  *
 117  *       beta^3*2*exp(-z^2)     beta^3*erf(z)
 118  *       ------------------  - ---------------
 119  *          sqrt(Pi)*z^2            z^3
 120  *
 121  *    or, switching back to r (z=r*beta):
 122  *
 123  *       2*beta*exp(-r^2*beta^2)   erf(r*beta)
 124  *       ----------------------- - -----------
 125  *            sqrt(Pi)*r^2            r^3
 126  *
 127  *
 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.
 130  *
 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.
 135  *
 136  */
 137 static gmx_mm_pr
 138 gmx_pmecorrF_pr(gmx_mm_pr z2)
 139 {
 140     const gmx_mm_pr  FN10     = gmx_set1_pr(-8.0072854618360083154e-14);
 141     const gmx_mm_pr  FN9      = gmx_set1_pr(1.1859116242260148027e-11);
 142     const gmx_mm_pr  FN8      = gmx_set1_pr(-8.1490406329798423616e-10);
 143     const gmx_mm_pr  FN7      = gmx_set1_pr(3.4404793543907847655e-8);
 144     const gmx_mm_pr  FN6      = gmx_set1_pr(-9.9471420832602741006e-7);
 145     const gmx_mm_pr  FN5      = gmx_set1_pr(0.000020740315999115847456);
 146     const gmx_mm_pr  FN4      = gmx_set1_pr(-0.00031991745139313364005);
 147     const gmx_mm_pr  FN3      = gmx_set1_pr(0.0035074449373659008203);
 148     const gmx_mm_pr  FN2      = gmx_set1_pr(-0.031750380176100813405);
 149     const gmx_mm_pr  FN1      = gmx_set1_pr(0.13884101728898463426);
 150     const gmx_mm_pr  FN0      = gmx_set1_pr(-0.75225277815249618847);
 151
 152     const gmx_mm_pr  FD5      = gmx_set1_pr(0.000016009278224355026701);
 153     const gmx_mm_pr  FD4      = gmx_set1_pr(0.00051055686934806966046);
 154     const gmx_mm_pr  FD3      = gmx_set1_pr(0.0081803507497974289008);
 155     const gmx_mm_pr  FD2      = gmx_set1_pr(0.077181146026670287235);
 156     const gmx_mm_pr  FD1      = gmx_set1_pr(0.41543303143712535988);
 157     const gmx_mm_pr  FD0      = gmx_set1_pr(1.0);
 158
 159     gmx_mm_pr        z4;
 160     gmx_mm_pr        polyFN0, polyFN1, polyFD0, polyFD1;
 161
 162     z4             = gmx_mul_pr(z2, z2);
 163
 164     polyFD1        = gmx_madd_pr(FD5, z4, FD3);
 165     polyFD1        = gmx_madd_pr(polyFD1, z4, FD1);
 166     polyFD1        = gmx_mul_pr(polyFD1, z2);
 167     polyFD0        = gmx_madd_pr(FD4, z4, FD2);
 168     polyFD0        = gmx_madd_pr(polyFD0, z4, FD0);
 169     polyFD0        = gmx_add_pr(polyFD0, polyFD1);
 170
 171     polyFD0        = gmx_inv_pr(polyFD0);
 172
 173     polyFN0        = gmx_madd_pr(FN10, z4, FN8);
 174     polyFN0        = gmx_madd_pr(polyFN0, z4, FN6);
 175     polyFN0        = gmx_madd_pr(polyFN0, z4, FN4);
 176     polyFN0        = gmx_madd_pr(polyFN0, z4, FN2);
 177     polyFN0        = gmx_madd_pr(polyFN0, z4, FN0);
 178     polyFN1        = gmx_madd_pr(FN9, z4, FN7);
 179     polyFN1        = gmx_madd_pr(polyFN1, z4, FN5);
 180     polyFN1        = gmx_madd_pr(polyFN1, z4, FN3);
 181     polyFN1        = gmx_madd_pr(polyFN1, z4, FN1);
 182     polyFN0        = gmx_madd_pr(polyFN1, z2, polyFN0);
 183
 184     return gmx_mul_pr(polyFN0, polyFD0);
 185 }
 186
 187
 188 /* Calculate the potential correction due to PME analytically.
 189  *
 190  * This routine calculates Erf(z)/z, although you should provide z^2
 191  * as the input argument.
 192  *
 193  * Here's how it should be used:
 194  *
 195  * 1. Calculate r^2.
 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:
 199  *
 200  *
 201  *        erf(z)
 202  *       --------
 203  *          z
 204  *
 205  * 5. Multiply the entire expression by beta and switching back to r (z=r*beta):
 206  *
 207  *       erf(r*beta)
 208  *       -----------
 209  *           r
 210  *
 211  * 6. Subtract the result from 1/r, multiply by the product of the charges,
 212  *    and you have your potential.
 213  *
 214  */
 215 static gmx_mm_pr
 216 gmx_pmecorrV_pr(gmx_mm_pr z2)
 217 {
 218     const gmx_mm_pr  VN9      = gmx_set1_pr(-9.3723776169321855475e-13);
 219     const gmx_mm_pr  VN8      = gmx_set1_pr(1.2280156762674215741e-10);
 220     const gmx_mm_pr  VN7      = gmx_set1_pr(-7.3562157912251309487e-9);
 221     const gmx_mm_pr  VN6      = gmx_set1_pr(2.6215886208032517509e-7);
 222     const gmx_mm_pr  VN5      = gmx_set1_pr(-4.9532491651265819499e-6);
 223     const gmx_mm_pr  VN4      = gmx_set1_pr(0.00025907400778966060389);
 224     const gmx_mm_pr  VN3      = gmx_set1_pr(0.0010585044856156469792);
 225     const gmx_mm_pr  VN2      = gmx_set1_pr(0.045247661136833092885);
 226     const gmx_mm_pr  VN1      = gmx_set1_pr(0.11643931522926034421);
 227     const gmx_mm_pr  VN0      = gmx_set1_pr(1.1283791671726767970);
 228
 229     const gmx_mm_pr  VD5      = gmx_set1_pr(0.000021784709867336150342);
 230     const gmx_mm_pr  VD4      = gmx_set1_pr(0.00064293662010911388448);
 231     const gmx_mm_pr  VD3      = gmx_set1_pr(0.0096311444822588683504);
 232     const gmx_mm_pr  VD2      = gmx_set1_pr(0.085608012351550627051);
 233     const gmx_mm_pr  VD1      = gmx_set1_pr(0.43652499166614811084);
 234     const gmx_mm_pr  VD0      = gmx_set1_pr(1.0);
 235
 236     gmx_mm_pr        z4;
 237     gmx_mm_pr        polyVN0, polyVN1, polyVD0, polyVD1;
 238
 239     z4             = gmx_mul_pr(z2, z2);
 240
 241     polyVD1        = gmx_madd_pr(VD5, z4, VD3);
 242     polyVD0        = gmx_madd_pr(VD4, z4, VD2);
 243     polyVD1        = gmx_madd_pr(polyVD1, z4, VD1);
 244     polyVD0        = gmx_madd_pr(polyVD0, z4, VD0);
 245     polyVD0        = gmx_madd_pr(polyVD1, z2, polyVD0);
 246
 247     polyVD0        = gmx_inv_pr(polyVD0);
 248
 249     polyVN1        = gmx_madd_pr(VN9, z4, VN7);
 250     polyVN0        = gmx_madd_pr(VN8, z4, VN6);
 251     polyVN1        = gmx_madd_pr(polyVN1, z4, VN5);
 252     polyVN0        = gmx_madd_pr(polyVN0, z4, VN4);
 253     polyVN1        = gmx_madd_pr(polyVN1, z4, VN3);
 254     polyVN0        = gmx_madd_pr(polyVN0, z4, VN2);
 255     polyVN1        = gmx_madd_pr(polyVN1, z4, VN1);
 256     polyVN0        = gmx_madd_pr(polyVN0, z4, VN0);
 257     polyVN0        = gmx_madd_pr(polyVN1, z2, polyVN0);
 258
 259     return gmx_mul_pr(polyVN0, polyVD0);
 260 }
 261
 262
 263 #endif /*_gmx_simd_math_double_h_ */