include/gmx_simd_math_single.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_single_h_
  36 #define _gmx_simd_math_single_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 one   = gmx_set1_pr(1.0);
  45
  46     gmx_mm_pr       lu = gmx_rsqrt_pr(x);
  47
  48     return gmx_madd_pr(gmx_nmsub_pr(x, gmx_mul_pr(lu, lu), one), gmx_mul_pr(lu, half), lu);
  49 }
  50
  51
  52 /* 1.0/x */
  53 static gmx_inline gmx_mm_pr
  54 gmx_inv_pr(gmx_mm_pr x)
  55 {
  56     const gmx_mm_pr two = gmx_set1_pr(2.0);
  57
  58     gmx_mm_pr       lu = gmx_rcp_pr(x);
  59
  60     return gmx_mul_pr(lu, gmx_nmsub_pr(lu, x, two));
  61 }
  62
  63
  64 /* Calculate the force correction due to PME analytically.
  65  *
  66  * This routine is meant to enable analytical evaluation of the
  67  * direct-space PME electrostatic force to avoid tables.
  68  *
  69  * The direct-space potential should be Erfc(beta*r)/r, but there
  70  * are some problems evaluating that:
  71  *
  72  * First, the error function is difficult (read: expensive) to
  73  * approxmiate accurately for intermediate to large arguments, and
  74  * this happens already in ranges of beta*r that occur in simulations.
  75  * Second, we now try to avoid calculating potentials in Gromacs but
  76  * use forces directly.
  77  *
  78  * We can simply things slight by noting that the PME part is really
  79  * a correction to the normal Coulomb force since Erfc(z)=1-Erf(z), i.e.
  80  *
  81  * V= 1/r - Erf(beta*r)/r
  82  *
  83  * The first term we already have from the inverse square root, so
  84  * that we can leave out of this routine.
  85  *
  86  * For pme tolerances of 1e-3 to 1e-8 and cutoffs of 0.5nm to 1.8nm,
  87  * the argument beta*r will be in the range 0.15 to ~4. Use your
  88  * favorite plotting program to realize how well-behaved Erf(z)/z is
  89  * in this range!
  90  *
  91  * We approximate f(z)=erf(z)/z with a rational minimax polynomial.
  92  * However, it turns out it is more efficient to approximate f(z)/z and
  93  * then only use even powers. This is another minor optimization, since
  94  * we actually WANT f(z)/z, because it is going to be multiplied by
  95  * the vector between the two atoms to get the vectorial force. The
  96  * fastest flops are the ones we can avoid calculating!
  97  *
  98  * So, here's how it should be used:
  99  *
 100  * 1. Calculate r^2.
 101  * 2. Multiply by beta^2, so you get z^2=beta^2*r^2.
 102  * 3. Evaluate this routine with z^2 as the argument.
 103  * 4. The return value is the expression:
 104  *
 105  *
 106  *       2*exp(-z^2)     erf(z)
 107  *       ------------ - --------
 108  *       sqrt(Pi)*z^2      z^3
 109  *
 110  * 5. Multiply the entire expression by beta^3. This will get you
 111  *
 112  *       beta^3*2*exp(-z^2)     beta^3*erf(z)
 113  *       ------------------  - ---------------
 114  *          sqrt(Pi)*z^2            z^3
 115  *
 116  *    or, switching back to r (z=r*beta):
 117  *
 118  *       2*beta*exp(-r^2*beta^2)   erf(r*beta)
 119  *       ----------------------- - -----------
 120  *            sqrt(Pi)*r^2            r^3
 121  *
 122  *
 123  *    With a bit of math exercise you should be able to confirm that
 124  *    this is exactly D[Erf[beta*r]/r,r] divided by r another time.
 125  *
 126  * 6. Add the result to 1/r^3, multiply by the product of the charges,
 127  *    and you have your force (divided by r). A final multiplication
 128  *    with the vector connecting the two particles and you have your
 129  *    vectorial force to add to the particles.
 130  *
 131  */
 132 static gmx_mm_pr
 133 gmx_pmecorrF_pr(gmx_mm_pr z2)
 134 {
 135     const gmx_mm_pr  FN6      = gmx_set1_pr(-1.7357322914161492954e-8f);
 136     const gmx_mm_pr  FN5      = gmx_set1_pr(1.4703624142580877519e-6f);
 137     const gmx_mm_pr  FN4      = gmx_set1_pr(-0.000053401640219807709149f);
 138     const gmx_mm_pr  FN3      = gmx_set1_pr(0.0010054721316683106153f);
 139     const gmx_mm_pr  FN2      = gmx_set1_pr(-0.019278317264888380590f);
 140     const gmx_mm_pr  FN1      = gmx_set1_pr(0.069670166153766424023f);
 141     const gmx_mm_pr  FN0      = gmx_set1_pr(-0.75225204789749321333f);
 142
 143     const gmx_mm_pr  FD4      = gmx_set1_pr(0.0011193462567257629232f);
 144     const gmx_mm_pr  FD3      = gmx_set1_pr(0.014866955030185295499f);
 145     const gmx_mm_pr  FD2      = gmx_set1_pr(0.11583842382862377919f);
 146     const gmx_mm_pr  FD1      = gmx_set1_pr(0.50736591960530292870f);
 147     const gmx_mm_pr  FD0      = gmx_set1_pr(1.0f);
 148
 149     gmx_mm_pr        z4;
 150     gmx_mm_pr        polyFN0, polyFN1, polyFD0, polyFD1;
 151
 152     z4             = gmx_mul_pr(z2, z2);
 153
 154     polyFD0        = gmx_madd_pr(FD4, z4, FD2);
 155     polyFD1        = gmx_madd_pr(FD3, z4, FD1);
 156     polyFD0        = gmx_madd_pr(polyFD0, z4, FD0);
 157     polyFD0        = gmx_madd_pr(polyFD1, z2, polyFD0);
 158
 159     polyFD0        = gmx_inv_pr(polyFD0);
 160
 161     polyFN0        = gmx_madd_pr(FN6, z4, FN4);
 162     polyFN1        = gmx_madd_pr(FN5, z4, FN3);
 163     polyFN0        = gmx_madd_pr(polyFN0, z4, FN2);
 164     polyFN1        = gmx_madd_pr(polyFN1, z4, FN1);
 165     polyFN0        = gmx_madd_pr(polyFN0, z4, FN0);
 166     polyFN0        = gmx_madd_pr(polyFN1, z2, polyFN0);
 167
 168     return gmx_mul_pr(polyFN0, polyFD0);
 169 }
 170
 171
 172 /* Calculate the potential correction due to PME analytically.
 173  *
 174  * See gmx_pmecorrF_pr() for details about the approximation.
 175  *
 176  * This routine calculates Erf(z)/z, although you should provide z^2
 177  * as the input argument.
 178  *
 179  * Here's how it should be used:
 180  *
 181  * 1. Calculate r^2.
 182  * 2. Multiply by beta^2, so you get z^2=beta^2*r^2.
 183  * 3. Evaluate this routine with z^2 as the argument.
 184  * 4. The return value is the expression:
 185  *
 186  *
 187  *        erf(z)
 188  *       --------
 189  *          z
 190  *
 191  * 5. Multiply the entire expression by beta and switching back to r (z=r*beta):
 192  *
 193  *       erf(r*beta)
 194  *       -----------
 195  *           r
 196  *
 197  * 6. Add the result to 1/r, multiply by the product of the charges,
 198  *    and you have your potential.
 199  */
 200 static gmx_mm_pr
 201 gmx_pmecorrV_pr(gmx_mm_pr z2)
 202 {
 203     const gmx_mm_pr  VN6      = gmx_set1_pr(1.9296833005951166339e-8f);
 204     const gmx_mm_pr  VN5      = gmx_set1_pr(-1.4213390571557850962e-6f);
 205     const gmx_mm_pr  VN4      = gmx_set1_pr(0.000041603292906656984871f);
 206     const gmx_mm_pr  VN3      = gmx_set1_pr(-0.00013134036773265025626f);
 207     const gmx_mm_pr  VN2      = gmx_set1_pr(0.038657983986041781264f);
 208     const gmx_mm_pr  VN1      = gmx_set1_pr(0.11285044772717598220f);
 209     const gmx_mm_pr  VN0      = gmx_set1_pr(1.1283802385263030286f);
 210
 211     const gmx_mm_pr  VD3      = gmx_set1_pr(0.0066752224023576045451f);
 212     const gmx_mm_pr  VD2      = gmx_set1_pr(0.078647795836373922256f);
 213     const gmx_mm_pr  VD1      = gmx_set1_pr(0.43336185284710920150f);
 214     const gmx_mm_pr  VD0      = gmx_set1_pr(1.0f);
 215
 216     gmx_mm_pr        z4;
 217     gmx_mm_pr        polyVN0, polyVN1, polyVD0, polyVD1;
 218
 219     z4             = gmx_mul_pr(z2, z2);
 220
 221     polyVD1        = gmx_madd_pr(VD3, z4, VD1);
 222     polyVD0        = gmx_madd_pr(VD2, z4, VD0);
 223     polyVD0        = gmx_madd_pr(polyVD1, z2, polyVD0);
 224
 225     polyVD0        = gmx_inv_pr(polyVD0);
 226
 227     polyVN0        = gmx_madd_pr(VN6, z4, VN4);
 228     polyVN1        = gmx_madd_pr(VN5, z4, VN3);
 229     polyVN0        = gmx_madd_pr(polyVN0, z4, VN2);
 230     polyVN1        = gmx_madd_pr(polyVN1, z4, VN1);
 231     polyVN0        = gmx_madd_pr(polyVN0, z4, VN0);
 232     polyVN0        = gmx_madd_pr(polyVN1, z2, polyVN0);
 233
 234     return gmx_mul_pr(polyVN0, polyVD0);
 235 }
 236
 237
 238 #endif /* _gmx_simd_math_single_h_ */