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 by the GROMACS development team.
7 * Copyright (c) 2019,2020, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
40 #include "calculate_spline_moduli.h"
46 #include "gromacs/math/utilities.h"
47 #include "gromacs/math/vec.h"
48 #include "gromacs/utility/exceptions.h"
49 #include "gromacs/utility/gmxassert.h"
50 #include "gromacs/utility/smalloc.h"
52 static void make_dft_mod(real
* mod
, const double* data
, int splineOrder
, int ndata
)
54 for (int i
= 0; i
< ndata
; i
++)
56 /* We use double precision, since this is only called once per grid.
57 * But for single precision bsp_mod, single precision also seems
58 * to give full accuracy.
62 for (int j
= 0; j
< splineOrder
; j
++)
64 double arg
= (2.0 * M_PI
* i
* (j
+ 1)) / ndata
;
65 sc
+= data
[j
] * cos(arg
);
66 ss
+= data
[j
] * sin(arg
);
68 mod
[i
] = sc
* sc
+ ss
* ss
;
70 if (splineOrder
% 2 == 0 && ndata
% 2 == 0)
72 /* Note that pme_order = splineOrder + 1 */
73 GMX_RELEASE_ASSERT(mod
[ndata
/ 2] < GMX_DOUBLE_EPS
,
74 "With even spline order and even grid size (ndata), dft_mod[ndata/2] "
75 "should first come out as zero");
76 /* This factor causes a division by zero. But since this occurs in
77 * the tail of the distribution, the term with this factor can
78 * be ignored (see Essmann et al. JCP 103, 8577).
79 * Using the average of the neighbors probably originates from
80 * Tom Darden's original PME code. It seems to give slighlty better
81 * accuracy than using a large value.
83 mod
[ndata
/ 2] = (mod
[ndata
/ 2 - 1] + mod
[ndata
/ 2 + 1]) * 0.5;
87 void make_bspline_moduli(splinevec bsp_mod
, int nx
, int ny
, int nz
, int pme_order
)
89 /* We use double precision, since this is only called once per grid.
90 * But for single precision bsp_mod, single precision also seems
91 * to give full accuracy.
95 /* In GROMACS we, confusingly, defined pme-order as the order
96 * of the cardinal B-spline + 1. This probably happened because
97 * the smooth PME paper only talks about "n" which is the number
98 * of points we spread to and that was chosen to be pme-order.
100 const int splineOrder
= pme_order
- 1;
102 snew(data
, splineOrder
);
105 for (int k
= 1; k
< splineOrder
; k
++)
110 for (int k
= 2; k
<= splineOrder
; k
++)
112 double div
= 1.0 / k
;
113 for (int m
= k
- 1; m
> 0; m
--)
115 data
[m
] = div
* ((k
- m
) * data
[m
- 1] + (m
+ 1) * data
[m
]);
117 data
[0] = div
* data
[0];
120 make_dft_mod(bsp_mod
[XX
], data
, splineOrder
, nx
);
121 make_dft_mod(bsp_mod
[YY
], data
, splineOrder
, ny
);
122 make_dft_mod(bsp_mod
[ZZ
], data
, splineOrder
, nz
);
127 /* Return the P3M optimal influence function */
128 static double do_p3m_influence(double z
, int order
)
135 /* The formula and most constants can be found in:
136 * Ballenegger et al., JCTC 8, 936 (2012)
140 case 2: return 1.0 - 2.0 * z2
/ 3.0;
141 case 3: return 1.0 - z2
+ 2.0 * z4
/ 15.0;
142 case 4: return 1.0 - 4.0 * z2
/ 3.0 + 2.0 * z4
/ 5.0 + 4.0 * z2
* z4
/ 315.0;
144 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;
146 return 1.0 - 2.0 * z2
+ 19.0 * z4
/ 15.0 - 256.0 * z2
* z4
/ 945.0
147 + 62.0 * z4
* z4
/ 4725.0 + 4.0 * z2
* z4
* z4
/ 155925.0;
149 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
150 - 2.0 * z2
* z4
* z4
/ 1485.0 + 4.0 * z4
* z4
* z4
/ 6081075.0;
152 return 1.0 - 8.0 * z2
/ 3.0 + 116.0 * z4
/ 45.0 - 344.0 * z2
* z4
/ 315.0
153 + 914.0 * z4
* z4
/ 4725.0 - 248.0 * z4
* z4
* z2
/ 22275.0
154 + 21844.0 * z4
* z4
* z4
/ 212837625.0 - 8.0 * z4
* z4
* z4
* z2
/ 638512875.0;
160 /* Calculate the P3M B-spline moduli for one dimension */
161 static void make_p3m_bspline_moduli_dim(real
* bsp_mod
, int n
, int order
)
163 double zarg
, zai
, sinzai
, infl
;
168 GMX_THROW(gmx::InconsistentInputError("The current P3M code only supports orders up to 8"));
175 for (i
= -maxk
; i
< 0; i
++)
179 infl
= do_p3m_influence(sinzai
, order
);
180 bsp_mod
[n
+ i
] = infl
* infl
* pow(sinzai
/ zai
, -2.0 * order
);
183 for (i
= 1; i
< maxk
; i
++)
187 infl
= do_p3m_influence(sinzai
, order
);
188 bsp_mod
[i
] = infl
* infl
* pow(sinzai
/ zai
, -2.0 * order
);
192 /* Calculate the P3M B-spline moduli */
193 void make_p3m_bspline_moduli(splinevec bsp_mod
, int nx
, int ny
, int nz
, int order
)
195 make_p3m_bspline_moduli_dim(bsp_mod
[XX
], nx
, order
);
196 make_p3m_bspline_moduli_dim(bsp_mod
[YY
], ny
, order
);
197 make_p3m_bspline_moduli_dim(bsp_mod
[ZZ
], nz
, order
);