2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017,2019, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
39 #include "calculate_spline_moduli.h"
45 #include "gromacs/math/utilities.h"
46 #include "gromacs/math/vec.h"
47 #include "gromacs/utility/exceptions.h"
48 #include "gromacs/utility/gmxassert.h"
49 #include "gromacs/utility/smalloc.h"
51 #include "pme_internal.h"
53 static void make_dft_mod(real* mod, const double* data, int splineOrder, int ndata)
55 for (int i = 0; i < ndata; i++)
57 /* We use double precision, since this is only called once per grid.
58 * But for single precision bsp_mod, single precision also seems
59 * to give full accuracy.
63 for (int j = 0; j < splineOrder; j++)
65 double arg = (2.0 * M_PI * i * (j + 1)) / ndata;
66 sc += data[j] * cos(arg);
67 ss += data[j] * sin(arg);
69 mod[i] = sc * sc + ss * ss;
71 if (splineOrder % 2 == 0 && ndata % 2 == 0)
73 /* Note that pme_order = splineOrder + 1 */
74 GMX_RELEASE_ASSERT(mod[ndata / 2] < GMX_DOUBLE_EPS,
75 "With even spline order and even grid size (ndata), dft_mod[ndata/2] "
76 "should first come out as zero");
77 /* This factor causes a division by zero. But since this occurs in
78 * the tail of the distribution, the term with this factor can
79 * be ignored (see Essmann et al. JCP 103, 8577).
80 * Using the average of the neighbors probably originates from
81 * Tom Darden's original PME code. It seems to give slighlty better
82 * accuracy than using a large value.
84 mod[ndata / 2] = (mod[ndata / 2 - 1] + mod[ndata / 2 + 1]) * 0.5;
88 void make_bspline_moduli(splinevec bsp_mod, int nx, int ny, int nz, int pme_order)
90 /* We use double precision, since this is only called once per grid.
91 * But for single precision bsp_mod, single precision also seems
92 * to give full accuracy.
96 /* In GROMACS we, confusingly, defined pme-order as the order
97 * of the cardinal B-spline + 1. This probably happened because
98 * the smooth PME paper only talks about "n" which is the number
99 * of points we spread to and that was chosen to be pme-order.
101 const int splineOrder = pme_order - 1;
103 snew(data, splineOrder);
106 for (int k = 1; k < splineOrder; k++)
111 for (int k = 2; k <= splineOrder; k++)
113 double div = 1.0 / k;
114 for (int m = k - 1; m > 0; m--)
116 data[m] = div * ((k - m) * data[m - 1] + (m + 1) * data[m]);
118 data[0] = div * data[0];
121 make_dft_mod(bsp_mod[XX], data, splineOrder, nx);
122 make_dft_mod(bsp_mod[YY], data, splineOrder, ny);
123 make_dft_mod(bsp_mod[ZZ], data, splineOrder, nz);
128 /* Return the P3M optimal influence function */
129 static double do_p3m_influence(double z, int order)
136 /* The formula and most constants can be found in:
137 * Ballenegger et al., JCTC 8, 936 (2012)
141 case 2: return 1.0 - 2.0 * z2 / 3.0;
142 case 3: return 1.0 - z2 + 2.0 * z4 / 15.0;
143 case 4: return 1.0 - 4.0 * z2 / 3.0 + 2.0 * z4 / 5.0 + 4.0 * z2 * z4 / 315.0;
145 return 1.0 - 5.0 * z2 / 3.0 + 7.0 * z4 / 9.0 - 17.0 * z2 * z4 / 189.0 + 2.0 * z4 * z4 / 2835.0;
147 return 1.0 - 2.0 * z2 + 19.0 * z4 / 15.0 - 256.0 * z2 * z4 / 945.0
148 + 62.0 * z4 * z4 / 4725.0 + 4.0 * z2 * z4 * z4 / 155925.0;
150 return 1.0 - 7.0 * z2 / 3.0 + 28.0 * z4 / 15.0 - 16.0 * z2 * z4 / 27.0 + 26.0 * z4 * z4 / 405.0
151 - 2.0 * z2 * z4 * z4 / 1485.0 + 4.0 * z4 * z4 * z4 / 6081075.0;
153 return 1.0 - 8.0 * z2 / 3.0 + 116.0 * z4 / 45.0 - 344.0 * z2 * z4 / 315.0
154 + 914.0 * z4 * z4 / 4725.0 - 248.0 * z4 * z4 * z2 / 22275.0
155 + 21844.0 * z4 * z4 * z4 / 212837625.0 - 8.0 * z4 * z4 * z4 * z2 / 638512875.0;
161 /* Calculate the P3M B-spline moduli for one dimension */
162 static void make_p3m_bspline_moduli_dim(real* bsp_mod, int n, int order)
164 double zarg, zai, sinzai, infl;
169 GMX_THROW(gmx::InconsistentInputError("The current P3M code only supports orders up to 8"));
176 for (i = -maxk; i < 0; i++)
180 infl = do_p3m_influence(sinzai, order);
181 bsp_mod[n + i] = infl * infl * pow(sinzai / zai, -2.0 * order);
184 for (i = 1; i < maxk; i++)
188 infl = do_p3m_influence(sinzai, order);
189 bsp_mod[i] = infl * infl * pow(sinzai / zai, -2.0 * order);
193 /* Calculate the P3M B-spline moduli */
194 void make_p3m_bspline_moduli(splinevec bsp_mod, int nx, int ny, int nz, int order)
196 make_p3m_bspline_moduli_dim(bsp_mod[XX], nx, order);
197 make_p3m_bspline_moduli_dim(bsp_mod[YY], ny, order);
198 make_p3m_bspline_moduli_dim(bsp_mod[ZZ], nz, order);