GADGET-4
ewald_test.cc
Go to the documentation of this file.
1/*******************************************************************************
2 * \copyright This file is part of the GADGET4 N-body/SPH code developed
3 * \copyright by Volker Springel. Copyright (C) 2014-2020 by Volker Springel
4 * \copyright (vspringel@mpa-garching.mpg.de) and all contributing authors.
5 *******************************************************************************/
6
12#include "gadgetconfig.h"
13
14#include <math.h>
15#include <mpi.h>
16#include <stdio.h>
17#include <stdlib.h>
18#include <string.h>
19
20#include "../data/allvars.h"
21#include "../data/dtypes.h"
22#include "../data/mymalloc.h"
23#include "../gravity/ewald.h"
24#include "../gravity/ewaldtensors.h"
25#include "../gravtree/gravtree.h"
26#include "../io/io.h"
27#include "../main/simulation.h"
28#include "../sort/cxxsort.h"
29#include "../system/system.h"
30
31#ifdef EWALD_TEST
32
45void ewald::ewald_corr(double dx, double dy, double dz, enum interpolate_options flag, ewald_data &fper)
46{
47 if(!ewald_is_initialized)
48 Terminate("How come that Ewald tables are not initialized?");
49
50 int signx, signy, signz;
51
52 if(dx < 0)
53 {
54 dx = -dx;
55 signx = -1;
56 }
57 else
58 signx = +1;
59
60 if(dy < 0)
61 {
62 dy = -dy;
63 signy = -1;
64 }
65 else
66 signy = +1;
67
68 if(dz < 0)
69 {
70 dz = -dz;
71 signz = -1;
72 }
73 else
74 signz = +1;
75
76 ewald_interpolate(dx, dy, dz, flag, fper);
77
78 /* change signs as needed */
79
80 fper.D1phi[0] *= signx;
81 fper.D1phi[1] *= signy;
82 fper.D1phi[2] *= signz;
83
84 if(flag == POINTMASS)
85 return;
86
87 fper.D2phi[qXY] *= signx * signy;
88 fper.D2phi[qXZ] *= signx * signz;
89 fper.D2phi[qYZ] *= signy * signz;
90
91 fper.D3phi[dXXX] *= signx;
92 fper.D3phi[dXXY] *= signy;
93 fper.D3phi[dXXZ] *= signz;
94 fper.D3phi[dXYY] *= signx;
95 fper.D3phi[dXYZ] *= signx * signy * signz;
96 fper.D3phi[dXZZ] *= signx;
97 fper.D3phi[dYYY] *= signy;
98 fper.D3phi[dYYZ] *= signz;
99 fper.D3phi[dYZZ] *= signy;
100 fper.D3phi[dZZZ] *= signz;
101
102 fper.D4phi[sXXXY] *= signx * signy;
103 fper.D4phi[sXYYY] *= signx * signy;
104 fper.D4phi[sXXXZ] *= signx * signz;
105 fper.D4phi[sXZZZ] *= signx * signz;
106 fper.D4phi[sYYYZ] *= signy * signz;
107 fper.D4phi[sYZZZ] *= signy * signz;
108 fper.D4phi[sXXYZ] *= signy * signz;
109 fper.D4phi[sXYYZ] *= signx * signz;
110 fper.D4phi[sXYZZ] *= signx * signy;
111
112 fper.D5phi[rXXXXX] *= signx;
113 fper.D5phi[rYYYYY] *= signy;
114 fper.D5phi[rZZZZZ] *= signz;
115
116 fper.D5phi[rXXXXY] *= signy;
117 fper.D5phi[rXXXXZ] *= signz;
118 fper.D5phi[rXYYYY] *= signx;
119 fper.D5phi[rXZZZZ] *= signx;
120 fper.D5phi[rYYYYZ] *= signz;
121 fper.D5phi[rYZZZZ] *= signy;
122
123 fper.D5phi[rXXXYY] *= signx;
124 fper.D5phi[rXXXZZ] *= signx;
125 fper.D5phi[rXXYYY] *= signy;
126 fper.D5phi[rXXZZZ] *= signz;
127 fper.D5phi[rYYYZZ] *= signy;
128 fper.D5phi[rYYZZZ] *= signz;
129
130 fper.D5phi[rXXYZZ] *= signy;
131 fper.D5phi[rXXYYZ] *= signz;
132 fper.D5phi[rXYYZZ] *= signx;
133
134 fper.D5phi[rXXXYZ] *= signx * signy * signz;
135 fper.D5phi[rXYYYZ] *= signx * signy * signz;
136 fper.D5phi[rXYZZZ] *= signx * signy * signz;
137}
138
139void ewald::ewald_corr_exact(double dx, double dy, double dz, enum interpolate_options flag, ewald_data &fper)
140{
141 double fac = 1.0 / All.BoxSize;
142 double x = dx * fac;
143 double y = dy * fac;
144 double z = dz * fac;
145
146 fper.D0phi = pow(fac, 1) * ewald_D0(x, y, z);
147 fper.D1phi = pow(fac, 2) * ewald_D1(x, y, z);
148
149 if(flag == POINTMASS)
150 return;
151
152 fper.D2phi = pow(fac, 3) * ewald_D2(x, y, z);
153 fper.D3phi = pow(fac, 4) * ewald_D3(x, y, z);
154 fper.D4phi = pow(fac, 5) * ewald_D4(x, y, z);
155 fper.D5phi = pow(fac, 6) * ewald_D5(x, y, z);
156}
157
158void ewald::ewald_interpolate(double dx, double dy, double dz, enum interpolate_options flag, ewald_data &fper)
159{
160 double u = dx * Ewd_fac_intp[0];
161 int i = (int)u;
162 if(i >= ENX)
163 i = ENX - 1;
164 u -= i;
165
166 double v = dy * Ewd_fac_intp[1];
167 int j = (int)v;
168 if(j >= ENY)
169 j = ENY - 1;
170 v -= j;
171
172 double w = dz * Ewd_fac_intp[2];
173 int k = (int)w;
174 if(k >= ENZ)
175 k = ENZ - 1;
176 w -= k;
177
178 double f1 = (1 - u) * (1 - v) * (1 - w);
179 double f2 = (1 - u) * (1 - v) * (w);
180 double f3 = (1 - u) * (v) * (1 - w);
181 double f4 = (1 - u) * (v) * (w);
182 double f5 = (u) * (1 - v) * (1 - w);
183 double f6 = (u) * (1 - v) * (w);
184 double f7 = (u) * (v) * (1 - w);
185 double f8 = (u) * (v) * (w);
186
187 ewald_data &C1 = Ewd[ewd_offset(i, j, k)];
188 ewald_data &C2 = Ewd[ewd_offset(i, j, k + 1)];
189 ewald_data &C3 = Ewd[ewd_offset(i, j + 1, k)];
190 ewald_data &C4 = Ewd[ewd_offset(i, j + 1, k + 1)];
191 ewald_data &C5 = Ewd[ewd_offset(i + 1, j, k)];
192 ewald_data &C6 = Ewd[ewd_offset(i + 1, j, k + 1)];
193 ewald_data &C7 = Ewd[ewd_offset(i + 1, j + 1, k)];
194 ewald_data &C8 = Ewd[ewd_offset(i + 1, j + 1, k + 1)];
195
196#ifdef EVALPOTENTIAL
197 fper.D0phi =
198 f1 * C1.D0phi + f2 * C2.D0phi + f3 * C3.D0phi + f4 * C4.D0phi + f5 * C5.D0phi + f6 * C6.D0phi + f7 * C7.D0phi + f8 * C8.D0phi;
199#endif
200 fper.D1phi =
201 f1 * C1.D1phi + f2 * C2.D1phi + f3 * C3.D1phi + f4 * C4.D1phi + f5 * C5.D1phi + f6 * C6.D1phi + f7 * C7.D1phi + f8 * C8.D1phi;
202
203 if(flag == POINTMASS)
204 return;
205
206 fper.D2phi =
207 f1 * C1.D2phi + f2 * C2.D2phi + f3 * C3.D2phi + f4 * C4.D2phi + f5 * C5.D2phi + f6 * C6.D2phi + f7 * C7.D2phi + f8 * C8.D2phi;
208
209 fper.D3phi =
210 f1 * C1.D3phi + f2 * C2.D3phi + f3 * C3.D3phi + f4 * C4.D3phi + f5 * C5.D3phi + f6 * C6.D3phi + f7 * C7.D3phi + f8 * C8.D3phi;
211
212 fper.D4phi =
213 f1 * C1.D4phi + f2 * C2.D4phi + f3 * C3.D4phi + f4 * C4.D4phi + f5 * C5.D4phi + f6 * C6.D4phi + f7 * C7.D4phi + f8 * C8.D4phi;
214
215 fper.D5phi =
216 f1 * C1.D5phi + f2 * C2.D5phi + f3 * C3.D5phi + f4 * C4.D5phi + f5 * C5.D5phi + f6 * C6.D5phi + f7 * C7.D5phi + f8 * C8.D5phi;
217}
218
219void ewald::test_interpolation_accuracy(void)
220{
221 init_rng(42);
222
223 printf("\n\n");
224
225 double errsum0 = 0;
226 double errmax0 = 0;
227 double count0 = 0;
228 for(int i = 0; i < 1000; i++)
229 {
230 double x = get_random_number() - 0.5;
231 double y = get_random_number() - 0.5;
232 double z = get_random_number() - 0.5;
233
234 double r = sqrt(x * x + y * y + z * z);
235
236 double D0phi_exact = (1.0 / All.BoxSize) * ewald_D0(x, y, z);
237
238 ewald_data ew;
239 ewald_corr(x * All.BoxSize, y * All.BoxSize, z * All.BoxSize, ewald::MULTIPOLES, ew);
240
241 double err = fabs((D0phi_exact - ew.D0phi) / D0phi_exact);
242
243 errsum0 += err;
244 if(err > errmax0)
245 errmax0 = err;
246 count0++;
247
248 if(err > 0.1)
249 printf("%4d r=%g D0_exact=%g D0_interpol=%g rel_error=%g\n", i, r, D0phi_exact, ew.D0phi, err);
250 }
251
252 double errsum1 = 0;
253 double errmax1 = 0;
254 double count1 = 0;
255 for(int i = 0; i < 1000; i++)
256 {
257 double x = get_random_number() - 0.5;
258 double y = get_random_number() - 0.5;
259 double z = get_random_number() - 0.5;
260
261 double r = sqrt(x * x + y * y + z * z);
262
263 vector<double> D1phi_exact = pow(All.BoxSize, -2) * ewald_D1(x, y, z);
264
265 ewald_data ew;
266 ewald_corr(x * All.BoxSize, y * All.BoxSize, z * All.BoxSize, ewald::MULTIPOLES, ew);
267
268 double norm = D1phi_exact.norm();
269
270 for(int j = 0; j < 3; j++)
271 {
272 double err = fabs(D1phi_exact[j] - ew.D1phi[j]) / norm;
273 errsum1 += err;
274 if(err > errmax1)
275 errmax1 = err;
276 count1++;
277
278 if(err > 0.1)
279 {
280 printf("%4d r=%g bin=%d D1_exact[%d]=%g D1_interpol[%d]=%g rel_error=%g\n", i, r,
281 (int)(r * All.BoxSize * Ewd_fac_intp[0]), j, D1phi_exact[j], j, ew.D1phi[j], err);
282 }
283 }
284 }
285
286 double errsum2 = 0;
287 double errmax2 = 0;
288 double count2 = 0;
289 for(int i = 0; i < 1000; i++)
290 {
291 double x = get_random_number() - 0.5;
292 double y = get_random_number() - 0.5;
293 double z = get_random_number() - 0.5;
294
295 double r = sqrt(x * x + y * y + z * z);
296
297 symtensor2<double> D2phi_exact = pow(All.BoxSize, -3) * ewald_D2(x, y, z);
298
299 ewald_data ew;
300 ewald_corr(x * All.BoxSize, y * All.BoxSize, z * All.BoxSize, ewald::MULTIPOLES, ew);
301
302 double norm = D2phi_exact.norm();
303
304 for(int j = 0; j < 6; j++)
305 {
306 double err = fabs((D2phi_exact[j] - ew.D2phi[j]) / norm);
307 errsum2 += err;
308 if(err > errmax2)
309 errmax2 = err;
310 count2++;
311
312 if(err > 0.1)
313 printf("%4d r=%g D2_exact[%d]=%g D2_interpol[%d]=%g rel_error=%g\n", i, r, j, D2phi_exact[j], j, ew.D2phi[j], err);
314 }
315 }
316
317 double errsum3 = 0;
318 double errmax3 = 0;
319 double count3 = 0;
320 for(int i = 0; i < 1000; i++)
321 {
322 double x = get_random_number() - 0.5;
323 double y = get_random_number() - 0.5;
324 double z = get_random_number() - 0.5;
325
326 double r = sqrt(x * x + y * y + z * z);
327
328 symtensor3<double> D3phi_exact = pow(All.BoxSize, -4) * ewald_D3(x, y, z);
329
330 ewald_data ew;
331 ewald_corr(x * All.BoxSize, y * All.BoxSize, z * All.BoxSize, ewald::MULTIPOLES, ew);
332
333 double norm = D3phi_exact.norm();
334
335 for(int j = 0; j < 10; j++)
336 {
337 double err = fabs((D3phi_exact[j] - ew.D3phi[j]) / norm);
338 errsum3 += err;
339 if(err > errmax3)
340 errmax3 = err;
341 count3++;
342
343 if(err > 0.1)
344 printf("%4d r=%g D3_exact[%d]=%g D3_interpol[%d]=%g rel_error=%g\n", i, r, j, D3phi_exact[j], j, ew.D3phi[j], err);
345 }
346 }
347
348 double errsum4 = 0;
349 double errmax4 = 0;
350 double count4 = 0;
351 for(int i = 0; i < 1000; i++)
352 {
353 double x = get_random_number() - 0.5;
354 double y = get_random_number() - 0.5;
355 double z = get_random_number() - 0.5;
356
357 double r = sqrt(x * x + y * y + z * z);
358
359 symtensor4<double> D4phi_exact = pow(All.BoxSize, -5) * ewald_D4(x, y, z);
360
361 ewald_data ew;
362 ewald_corr(x * All.BoxSize, y * All.BoxSize, z * All.BoxSize, ewald::MULTIPOLES, ew);
363
364 double norm = D4phi_exact.norm();
365
366 for(int j = 0; j < 15; j++)
367 {
368 double err = fabs((D4phi_exact[j] - ew.D4phi[j]) / norm);
369 errsum4 += err;
370 if(err > errmax4)
371 errmax4 = err;
372 count4++;
373
374 if(err > 0.1)
375 printf("%4d r=%g D4_exact[%d]=%g D4_interpol[%d]=%g rel_error=%g\n", i, r, j, D4phi_exact[j], j, ew.D4phi[j], err);
376 }
377 }
378
379 double errsum5 = 0;
380 double errmax5 = 0;
381 double count5 = 0;
382 for(int i = 0; i < 1000; i++)
383 {
384 double x = get_random_number() - 0.5;
385 double y = get_random_number() - 0.5;
386 double z = get_random_number() - 0.5;
387
388 symtensor5<double> D5phi_exact = pow(All.BoxSize, -6) * ewald_D5(x, y, z);
389
390 ewald_data ew;
391 ewald_corr(x * All.BoxSize, y * All.BoxSize, z * All.BoxSize, ewald::MULTIPOLES, ew);
392
393 double norm = D5phi_exact.norm();
394
395 for(int j = 0; j < 21; j++)
396 {
397 double err = fabs((D5phi_exact[j] - ew.D5phi[j]) / norm);
398 errsum5 += err;
399 if(err > errmax5)
400 errmax5 = err;
401 count5++;
402 }
403 }
404
405 printf("\n\n");
406
407 printf("D0: max error = %g mean error=%g\n", errmax0, errsum0 / count0);
408 printf("D1: max error = %g mean error=%g\n", errmax1, errsum1 / count1);
409 printf("D2: max error = %g mean error=%g\n", errmax2, errsum2 / count2);
410 printf("D3: max error = %g mean error=%g\n", errmax3, errsum3 / count3);
411 printf("D4: max error = %g mean error=%g\n", errmax4, errsum4 / count4);
412 printf("D5: max error = %g mean error=%g\n", errmax5, errsum5 / count5);
413
414 printf("\n\n");
415
416 {
417 double errsum0 = 0;
418 double errmax0 = 0;
419 double count0 = 0;
420 for(int i = 0; i < 1000; i++)
421 {
422 double x = get_random_number() - 0.5;
423 double y = get_random_number() - 0.5;
424 double z = get_random_number() - 0.5;
425
426 double D0phi_exact = (1.0 / All.BoxSize) * ewald_D0(x, y, z);
427
428 ewald_data ew;
429 double posd[3] = {x * All.BoxSize, y * All.BoxSize, z * All.BoxSize};
430 MyIntPosType pos[3];
431 pos_to_signedintpos(posd, (MySignedIntPosType *)pos);
432 MyIntPosType ref[3] = {0, 0, 0};
433 ewald_gridlookup(pos, ref, ewald::MULTIPOLES, ew);
434
435 double err = fabs((D0phi_exact - ew.D0phi) / D0phi_exact);
436
437 errsum0 += err;
438 if(err > errmax0)
439 errmax0 = err;
440 count0++;
441 }
442
443 double errsum1 = 0;
444 double errmax1 = 0;
445 double count1 = 0;
446 for(int i = 0; i < 1000; i++)
447 {
448 double x = get_random_number() - 0.5;
449 double y = get_random_number() - 0.5;
450 double z = get_random_number() - 0.5;
451
452 vector<double> D1phi_exact = pow(All.BoxSize, -2) * ewald_D1(x, y, z);
453
454 ewald_data ew;
455 double posd[3] = {x * All.BoxSize, y * All.BoxSize, z * All.BoxSize};
456 MyIntPosType pos[3];
457 pos_to_signedintpos(posd, (MySignedIntPosType *)pos);
458 MyIntPosType ref[3] = {0, 0, 0};
459 ewald_gridlookup(pos, ref, ewald::MULTIPOLES, ew);
460
461 double norm = D1phi_exact.norm();
462
463 for(int j = 0; j < 3; j++)
464 {
465 double err = fabs(D1phi_exact[j] - ew.D1phi[j]) / norm;
466
467 errsum1 += err;
468 if(err > errmax1)
469 errmax1 = err;
470 count1++;
471 }
472 }
473
474 double errsum2 = 0;
475 double errmax2 = 0;
476 double count2 = 0;
477 for(int i = 0; i < 1000; i++)
478 {
479 double x = get_random_number() - 0.5;
480 double y = get_random_number() - 0.5;
481 double z = get_random_number() - 0.5;
482
483 symtensor2<double> D2phi_exact = pow(All.BoxSize, -3) * ewald_D2(x, y, z);
484
485 ewald_data ew;
486 double posd[3] = {x * All.BoxSize, y * All.BoxSize, z * All.BoxSize};
487 MyIntPosType pos[3];
488 pos_to_signedintpos(posd, (MySignedIntPosType *)pos);
489 MyIntPosType ref[3] = {0, 0, 0};
490 ewald_gridlookup(pos, ref, ewald::MULTIPOLES, ew);
491
492 double norm = D2phi_exact.norm();
493
494 for(int j = 0; j < 6; j++)
495 {
496 double err = fabs((D2phi_exact[j] - ew.D2phi[j]) / norm);
497 errsum2 += err;
498 if(err > errmax2)
499 errmax2 = err;
500 count2++;
501 }
502 }
503
504 double errsum3 = 0;
505 double errmax3 = 0;
506 double count3 = 0;
507 for(int i = 0; i < 1000; i++)
508 {
509 double x = get_random_number() - 0.5;
510 double y = get_random_number() - 0.5;
511 double z = get_random_number() - 0.5;
512
513 symtensor3<double> D3phi_exact = pow(All.BoxSize, -4) * ewald_D3(x, y, z);
514
515 ewald_data ew;
516 double posd[3] = {x * All.BoxSize, y * All.BoxSize, z * All.BoxSize};
517 MyIntPosType pos[3];
518 pos_to_signedintpos(posd, (MySignedIntPosType *)pos);
519 MyIntPosType ref[3] = {0, 0, 0};
520 ewald_gridlookup(pos, ref, ewald::MULTIPOLES, ew);
521
522 double norm = D3phi_exact.norm();
523
524 for(int j = 0; j < 10; j++)
525 {
526 double err = fabs((D3phi_exact[j] - ew.D3phi[j]) / norm);
527 errsum3 += err;
528 if(err > errmax3)
529 errmax3 = err;
530 count3++;
531 }
532 }
533
534 double errsum4 = 0;
535 double errmax4 = 0;
536 double count4 = 0;
537 for(int i = 0; i < 1000; i++)
538 {
539 double x = get_random_number() - 0.5;
540 double y = get_random_number() - 0.5;
541 double z = get_random_number() - 0.5;
542
543 symtensor4<double> D4phi_exact = pow(All.BoxSize, -5) * ewald_D4(x, y, z);
544
545 ewald_data ew;
546 double posd[3] = {x * All.BoxSize, y * All.BoxSize, z * All.BoxSize};
547 MyIntPosType pos[3];
548 pos_to_signedintpos(posd, (MySignedIntPosType *)pos);
549 MyIntPosType ref[3] = {0, 0, 0};
550 ewald_gridlookup(pos, ref, ewald::MULTIPOLES, ew);
551
552 double norm = D4phi_exact.norm();
553
554 for(int j = 0; j < 15; j++)
555 {
556 double err = fabs((D4phi_exact[j] - ew.D4phi[j]) / norm);
557 errsum4 += err;
558 if(err > errmax4)
559 errmax4 = err;
560 count4++;
561 }
562 }
563
564 double errsum5 = 0;
565 double errmax5 = 0;
566 double count5 = 0;
567 for(int i = 0; i < 1000; i++)
568 {
569 double x = get_random_number() - 0.5;
570 double y = get_random_number() - 0.5;
571 double z = get_random_number() - 0.5;
572
573 symtensor5<double> D5phi_exact = pow(All.BoxSize, -6) * ewald_D5(x, y, z);
574
575 ewald_data ew;
576 double posd[3] = {x * All.BoxSize, y * All.BoxSize, z * All.BoxSize};
577 MyIntPosType pos[3];
578 pos_to_signedintpos(posd, (MySignedIntPosType *)pos);
579 MyIntPosType ref[3] = {0, 0, 0};
580 ewald_gridlookup(pos, ref, ewald::MULTIPOLES, ew);
581
582 double norm = D5phi_exact.norm();
583
584 for(int j = 0; j < 21; j++)
585 {
586 double err = fabs((D5phi_exact[j] - ew.D5phi[j]) / norm);
587 errsum5 += err;
588 if(err > errmax5)
589 errmax5 = err;
590 count5++;
591 }
592 }
593
594 printf("Grid look-up: \n\n");
595
596 printf("D0: max error = %g mean error=%g\n", errmax0, errsum0 / count0);
597 printf("D1: max error = %g mean error=%g\n", errmax1, errsum1 / count1);
598 printf("D2: max error = %g mean error=%g\n", errmax2, errsum2 / count2);
599 printf("D3: max error = %g mean error=%g\n", errmax3, errsum3 / count3);
600 printf("D4: max error = %g mean error=%g\n", errmax4, errsum4 / count4);
601 printf("D5: max error = %g mean error=%g\n", errmax5, errsum5 / count5);
602 }
603
604 Terminate("stop");
605}
606
607ewaldtensor6<double> ewald::ewald_P6(void)
608{
609#ifdef GRAVITY_TALLBOX
610 Terminate("GRAVITY_TALLBOX is not implemented");
611#endif
612
613 ewaldtensor6<double> P6 = 0.0;
614
615 double leff = pow((1.0 / LONG_X) * (1.0 / LONG_Y) * (1.0 / LONG_Z), 1.0 / 3);
616 double alpha = 2.0 / leff;
617 double alpha2 = alpha * alpha;
618
619 int qxmax = (int)(8.0 * LONG_X / alpha + 0.5);
620 int qymax = (int)(8.0 * LONG_Y / alpha + 0.5);
621 int qzmax = (int)(8.0 * LONG_Z / alpha + 0.5);
622
623 int nxmax = (int)(2.0 * alpha / LONG_X + 0.5);
624 int nymax = (int)(2.0 * alpha / LONG_Y + 0.5);
625 int nzmax = (int)(2.0 * alpha / LONG_Z + 0.5);
626
627 for(int nx = -qxmax; nx <= qxmax; nx++)
628 for(int ny = -qymax; ny <= qymax; ny++)
629 for(int nz = -qzmax; nz <= qzmax; nz++)
630 {
631 double dx = -nx * (1.0 / LONG_X);
632 double dy = -ny * (1.0 / LONG_Y);
633 double dz = -nz * (1.0 / LONG_Z);
634
635 vector<double> dxyz(dx, dy, dz);
636
637 double r2 = dx * dx + dy * dy + dz * dz;
638 double r = sqrt(r2);
639
640 double rinv = (r > 0) ? 1.0 / r : 0.0;
641 double r2inv = rinv * rinv;
642 double r3inv = r2inv * rinv;
643 double r4inv = r2inv * r2inv;
644 double r7inv = r3inv * r4inv;
645
646 if(nx != 0 || ny != 0 || nz != 0)
647 {
648 // derivatives of f(r) = -Erfc[alpha * r] /r
649
650 double ar = alpha * r;
651 double ar2 = ar * ar;
652 double ar3 = ar * ar * ar;
653 double ar5 = ar3 * ar * ar;
654 double ar7 = ar5 * ar * ar;
655 double ar9 = ar7 * ar * ar;
656 double ar11 = ar9 * ar * ar;
657 double xir2 = pow(dx * rinv, 2);
658 double xir4 = pow(dx * rinv, 4);
659 double xir6 = pow(dx * rinv, 6);
660
661 double yir2 = pow(dy * rinv, 2);
662 double zir2 = pow(dz * rinv, 2);
663
664 P6.XXXXXX += (450.0 * ar + 300.0 * ar3 + 120.0 * ar5 - xir2 * (9450.0 * ar + 6300.0 * ar3 + 2520.0 * ar5 + 720.0 * ar7) +
665 xir4 * (28350.0 * ar + 18900.0 * ar3 + 7560.0 * ar5 + 2160.0 * ar7 + 480.0 * ar9) -
666 xir6 * (20790.0 * ar + 13860.0 * ar3 + 5544.0 * ar5 + 1584.0 * ar7 + 352.0 * ar9 + 64.0 * ar11)) /
667 sqrt(M_PI) * exp(-ar2) * r7inv +
668 erfc(ar) * (225.0 - 4725.0 * xir2 + 14175.0 * xir4 - 10395.0 * xir6) * r7inv;
669
670 P6.XXXXYY += (90.0 * ar + 60.0 * ar3 + 24.0 * ar5 - xir2 * (1260.0 * ar + 840.0 * ar3 + 336.0 * ar5 + 96.0 * ar7) +
671 xir4 * (1890.0 * ar + 1260.0 * ar3 + 504.0 * ar5 + 144.0 * ar7 + 32.0 * ar9) -
672 yir2 * (630.0 * ar + 420.0 * ar3 + 168.0 * ar5 + 48.0 * ar7) +
673 xir2 * yir2 * (11340.0 * ar + 7560.0 * ar3 + 3024.0 * ar5 + 864.0 * ar7 + 192.0 * ar9) -
674 xir4 * yir2 * (20790.0 * ar + 13860.0 * ar3 + 5544.0 * ar5 + 1584.0 * ar7 + 352.0 * ar9 + 64.0 * ar11)) /
675 sqrt(M_PI) * exp(-ar2) * r7inv +
676 erfc(ar) *
677 (45.0 - 630.0 * xir2 + 945.0 * xir4 - 315.0 * yir2 + 5670.0 * xir2 * yir2 - 10395.0 * xir4 * yir2) *
678 r7inv;
679
680 P6.XXYYZZ +=
681 (30.0 * ar + 20.0 * ar3 + 8.0 * ar5 - (xir2 + yir2 + zir2) * (210.0 * ar + 140.0 * ar3 + 56.0 * ar5 + 16.0 * ar7) +
682 +(xir2 * yir2 + xir2 * zir2 + yir2 * zir2) * (1890.0 * ar + 1260.0 * ar3 + 504.0 * ar5 + 144.0 * ar7 + 32.0 * ar9) -
683 xir2 * yir2 * zir2 * (20790.0 * ar + 13860.0 * ar3 + 5544.0 * ar5 + 1584.0 * ar7 + 352.0 * ar9 + 64.0 * ar11)) /
684 sqrt(M_PI) * exp(-ar2) * r7inv +
685 erfc(ar) *
686 (15.0 - 105.0 * (xir2 + yir2 + zir2) + 945.0 * (xir2 * yir2 + xir2 * zir2 + yir2 * zir2) -
687 10395.0 * xir2 * yir2 * zir2) *
688 r7inv;
689 }
690 else
691 {
692 /* we add the 1/r term here to the (0|0|0) entry, followed by differentiation, and the limit r->0 to obtain accurate
693 * results at the origin
694 */
695
696 /* Note, for small r:
697 *
698 * [1/- erfc(a r)]/r = 2 a/sqrt(pi) * [ 1 - (a r)^2/3 + (a r)^4 / 10 - (a r)^6 / 42 + (a r)^8 / 216 - ...]
699 *
700 */
701
702 P6.XXXXXX += (-240.0) * pow(alpha, 7) / (7.0 * sqrt(M_PI));
703 P6.XXXXYY += (-48.0) * pow(alpha, 7) / (7.0 * sqrt(M_PI));
704 P6.XXYYZZ += 0;
705 }
706 }
707
708 for(int nx = -nxmax; nx <= nxmax; nx++)
709 for(int ny = -nymax; ny <= nymax; ny++)
710 for(int nz = -nzmax; nz <= nzmax; nz++)
711 {
712 if(nx != 0 || ny != 0 || nz != 0)
713 {
714 double kx = (2.0 * M_PI * LONG_X) * nx;
715 double ky = (2.0 * M_PI * LONG_Y) * ny;
716 double kz = (2.0 * M_PI * LONG_Z) * nz;
717 double k2 = kx * kx + ky * ky + kz * kz;
718
719 double val = 4.0 * M_PI * (LONG_X * LONG_Y * LONG_Z) / k2 * exp(-k2 / (4.0 * alpha2));
720
721 P6.XXXXXX += val * pow(kx, 6);
722 P6.XXXXYY += val * pow(kx, 4) * pow(ky, 2);
723 P6.XXYYZZ += val * pow(kx, 2) * pow(ky, 2) * pow(kz, 2);
724 }
725 }
726
727 return P6;
728}
729
730ewaldtensor8<double> ewald::ewald_P8(void)
731{
732#ifdef GRAVITY_TALLBOX
733 Terminate("GRAVITY_TALLBOX is not implemented");
734#endif
735
736 ewaldtensor8<double> P8 = 0.0;
737
738 double leff = pow((1.0 / LONG_X) * (1.0 / LONG_Y) * (1.0 / LONG_Z), 1.0 / 3);
739 double alpha = 2.0 / leff;
740 double alpha2 = alpha * alpha;
741
742 int qxmax = (int)(8.0 * LONG_X / alpha + 0.5);
743 int qymax = (int)(8.0 * LONG_Y / alpha + 0.5);
744 int qzmax = (int)(8.0 * LONG_Z / alpha + 0.5);
745
746 int nxmax = (int)(2.0 * alpha / LONG_X + 0.5);
747 int nymax = (int)(2.0 * alpha / LONG_Y + 0.5);
748 int nzmax = (int)(2.0 * alpha / LONG_Z + 0.5);
749
750 for(int nx = -qxmax; nx <= qxmax; nx++)
751 for(int ny = -qymax; ny <= qymax; ny++)
752 for(int nz = -qzmax; nz <= qzmax; nz++)
753 {
754 double dx = -nx * (1.0 / LONG_X);
755 double dy = -ny * (1.0 / LONG_Y);
756 double dz = -nz * (1.0 / LONG_Z);
757
758 vector<double> dxyz(dx, dy, dz);
759
760 double r2 = dx * dx + dy * dy + dz * dz;
761 double r = sqrt(r2);
762
763 double rinv = (r > 0) ? 1.0 / r : 0.0;
764 double r2inv = rinv * rinv;
765 double r3inv = r2inv * rinv;
766 double r4inv = r2inv * r2inv;
767 double r5inv = r2inv * r3inv;
768 double r9inv = r4inv * r5inv;
769
770 if(nx != 0 || ny != 0 || nz != 0)
771 {
772 // derivatives of f(r) = -Erfc[alpha * r] /r
773
774 double ar = alpha * r;
775 double ar2 = ar * ar;
776 double ar3 = ar * ar * ar;
777 double ar5 = ar3 * ar * ar;
778 double ar7 = ar5 * ar * ar;
779 double ar9 = ar7 * ar * ar;
780 double ar11 = ar9 * ar * ar;
781 double ar13 = ar11 * ar * ar;
782 double ar15 = ar13 * ar * ar;
783 double xir2 = pow(dx * rinv, 2);
784 double xir4 = pow(dx * rinv, 4);
785 double xir6 = pow(dx * rinv, 6);
786 double xir8 = pow(dx * rinv, 8);
787
788 double yir2 = pow(dy * rinv, 2);
789 double yir4 = pow(dy * rinv, 4);
790 double zir2 = pow(dz * rinv, 2);
791
792 P8.XXXXXXXX +=
793 (-(22050.0 * ar + 14700.0 * ar3 + 5880.0 * ar5 + 1680.0 * ar7) +
794 xir2 * (793800.0 * ar + 529200.0 * ar3 + 211680.0 * ar5 + 60480.0 * ar7 + 13440.0 * ar9) -
795 xir4 * (4365900.0 * ar + 2910600.0 * ar3 + 1164240.0 * ar5 + 332640.0 * ar7 + 73920.0 * ar9 + 13440.0 * ar11) +
796 xir6 * (7567560.0 * ar + 5045040.0 * ar3 + 2018016.0 * ar5 + 576576.0 * ar7 + 128128.0 * ar9 + 23296.0 * ar11 +
797 3584 * ar13) -
798 xir8 * (4054050.0 * ar + 2702700.0 * ar3 + 1081080.0 * ar5 + 308880.0 * ar7 + 68640.0 * ar9 + 12480.0 * ar11 +
799 1920.0 * ar13 + 256.0 * ar15)) /
800 sqrt(M_PI) * exp(-ar2) * r9inv +
801 erfc(ar) * (-11025.0 + 396900.0 * xir2 - 2182950.0 * xir4 + 3783780.0 * xir6 - 2027025.0 * xir8) * r9inv;
802
803 P8.XXXXXXYY =
804 (-(3150.0 * ar + 2100.0 * ar3 + 840.0 * ar5 + 240.0 * ar7) +
805 xir2 * (85050.0 * ar + 56700.0 * ar3 + 22680.0 * ar5 + 6480.0 * ar7 + 1440.0 * ar9) -
806 xir4 * (311850.0 * ar + 207900.0 * ar3 + 83160.0 * ar5 + 23760.0 * ar7 + 5280.0 * ar9 + 960.0 * ar11) +
807 xir6 *
808 (270270.0 * ar + 180180.0 * ar3 + 72072.0 * ar5 + 20592.0 * ar7 + 4576.0 * ar9 + 832.0 * ar11 + 128.0 * ar13) +
809 yir2 * (28350.0 * ar + 18900 * ar3 + 7560.0 * ar5 + 2160 * ar7 + 480 * ar9) -
810 xir2 * yir2 * (935550 * ar + 623700.0 * ar3 + 249480.0 * ar5 + 71280.0 * ar7 + 15840.0 * ar9 + 2880.0 * ar11) +
811 xir4 * yir2 *
812 (4054050.0 * ar + 2702700.0 * ar3 + 1081080.0 * ar5 + 308880.0 * ar7 + 68640.0 * ar9 + 12480.0 * ar11 +
813 1920.0 * ar13) -
814 xir6 * yir2 *
815 (4054050.0 * ar + 2702700.0 * ar3 + 1081080.0 * ar5 + 308880.0 * ar7 + 68640.0 * ar9 + 12480.0 * ar11 +
816 1920.0 * ar13 + 256.0 * ar15)) /
817 sqrt(M_PI) * exp(-ar2) * r9inv +
818 erfc(ar) *
819 (-1575.0 + 42525.0 * xir2 - 155925.0 * xir4 + 135135.0 * xir6 + 14175.0 * yir2 - 467775.0 * xir2 * yir2 +
820 2027025.0 * xir4 * yir2 - 2027025 * xir6 * yir2) *
821 r9inv;
822
823 P8.XXXXYYYY =
824 (-(1890.0 * ar + 1260.0 * ar3 + 504.0 * ar5 + 144.0 * ar7) +
825 (xir2 + yir2) * (34020.0 * ar + 22680.0 * ar3 + 9072.0 * ar5 + 2592.0 * ar7 + 576.0 * ar9) -
826 (xir4 + yir4) * (62370.0 * ar + 41580.0 * ar3 + 16632.0 * ar5 + 4752.0 * ar7 + 1056.0 * ar9 + 192.0 * ar11) +
827 -xir2 * yir2 * (748440.0 * ar + 498960.0 * ar3 + 199584.0 * ar5 + 57024.0 * ar7 + 12672.0 * ar9 + 2304.0 * ar11) +
828
829 (xir4 * yir2 + xir2 * yir4) * (1621620.0 * ar + 1081080.0 * ar3 + 432432.0 * ar5 + 123552.0 * ar7 + 27456.0 * ar9 +
830 4992.0 * ar11 + 768.0 * ar13) -
831 xir4 * yir4 *
832 (4054050.0 * ar + 2702700.0 * ar3 + 1081080.0 * ar5 + 308880.0 * ar7 + 68640.0 * ar9 + 12480.0 * ar11 +
833 1920.0 * ar13 + 256.0 * ar15)) /
834 sqrt(M_PI) * exp(-ar2) * r9inv +
835 erfc(ar) *
836 (-945.0 + 17010.0 * (xir2 + yir2) - 31185.0 * (xir4 + yir4) - 374220.0 * xir2 * yir2 +
837 810810.0 * (xir4 * yir2 + xir2 * yir4) - 2027025.0 * xir4 * yir4) *
838 r9inv;
839
840 P8.XXXXYYZZ = (-(630.0 * ar + 420.0 * ar3 + 168.0 * ar5 + 48.0 * ar7) +
841 (xir2) * (11340.0 * ar + 7560.0 * ar3 + 3024.0 * ar5 + 864.0 * ar7 + 192.0 * ar9) -
842 (xir4) * (20790.0 * ar + 13860.0 * ar3 + 5544.0 * ar5 + 1584.0 * ar7 + 352.0 * ar9 + 64.0 * ar11) +
843 +(yir2 + zir2) * (5670.0 * ar + 3780.0 * ar3 + 1512.0 * ar5 + 432.0 * ar7 + 96.0 * ar9) -
844 (xir2 * yir2 + xir2 * zir2) *
845 (124740 * ar + 83160.0 * ar3 + 33264.0 * ar5 + 9504 * ar7 + 2112.0 * ar9 + 384.0 * ar11) +
846 (xir4 * yir2 + xir4 * zir2) * (270270.0 * ar + 180180.0 * ar3 + 72072 * ar5 + 20592.0 * ar7 +
847 4576.0 * ar9 + 832.0 * ar11 + 128.0 * ar13) -
848 yir2 * zir2 * (62370.0 * ar + 41580 * ar3 + 16632.0 * ar5 + 4752 * ar7 + 1056.0 * ar9 + 192.0 * ar11) +
849 (xir2 * yir2 * zir2) * (1621620.0 * ar + 1081080.0 * ar3 + 432432.0 * ar5 + 123552.0 * ar7 +
850 27456.0 * ar9 + 4992.0 * ar11 + 768.0 * ar13) -
851 xir4 * yir2 * zir2 *
852 (4054050.0 * ar + 2702700.0 * ar3 + 1081080.0 * ar5 + 308880.0 * ar7 + 68640.0 * ar9 +
853 12480.0 * ar11 + 1920.0 * ar13 + 256.0 * ar15)) /
854 sqrt(M_PI) * exp(-ar2) * r9inv +
855 erfc(ar) *
856 (-315.0 + 5670.0 * xir2 - 10395.0 * xir4 + 2835.0 * (yir2 + zir2) -
857 62370.0 * (xir2 * yir2 + xir2 * zir2) + 135135.0 * xir4 * (yir2 + zir2) - 31185.0 * yir2 * zir2 +
858 810810.0 * xir2 * yir2 * zir2 - 2027025 * xir4 * yir2 * zir2) *
859 r9inv;
860 }
861 else
862 {
863 /* we add the 1/r term here to the (0|0|0) entry, followed by differentiation, and the limit r->0 to obtain accurate
864 * results at the origin
865 */
866
867 /* Note, for small r:
868 *
869 * [1/- erfc(a r)]/r = 2 a/sqrt(pi) * [ 1 - (a r)^2/3 + (a r)^4 / 10 - (a r)^6 / 42 + (a r)^8 / 216 - ...]
870 *
871 */
872
873 P8.XXXXXXXX += 1120.0 * pow(alpha, 9) / (3.0 * sqrt(M_PI));
874 P8.XXXXXXYY += 160.0 * pow(alpha, 9) / (3.0 * sqrt(M_PI));
875 P8.XXXXYYYY += 0;
876 P8.XXXXYYZZ += 32.0 * pow(alpha, 9) / (3.0 * sqrt(M_PI));
877 }
878 }
879
880 for(int nx = -nxmax; nx <= nxmax; nx++)
881 for(int ny = -nymax; ny <= nymax; ny++)
882 for(int nz = -nzmax; nz <= nzmax; nz++)
883 {
884 if(nx != 0 || ny != 0 || nz != 0)
885 {
886 double kx = (2.0 * M_PI * LONG_X) * nx;
887 double ky = (2.0 * M_PI * LONG_Y) * ny;
888 double kz = (2.0 * M_PI * LONG_Z) * nz;
889 double k2 = kx * kx + ky * ky + kz * kz;
890
891 double val = -4.0 * M_PI * (LONG_X * LONG_Y * LONG_Z) / k2 * exp(-k2 / (4.0 * alpha2));
892
893 P8.XXXXXXXX += val * pow(kx, 8);
894 P8.XXXXXXYY += val * pow(kx, 6) * pow(ky, 2);
895 P8.XXXXYYYY += val * pow(kx, 4) * pow(ky, 4);
896 P8.XXXXYYZZ += val * pow(kx, 4) * pow(ky, 2) * pow(kz, 2);
897 }
898 }
899
900 return P8;
901}
902
903ewaldtensor10<double> ewald::ewald_P10(void)
904{
905#ifdef GRAVITY_TALLBOX
906 Terminate("GRAVITY_TALLBOX is not implemented");
907#endif
908
909 ewaldtensor10<double> P10 = 0.0;
910
911 double leff = pow((1.0 / LONG_X) * (1.0 / LONG_Y) * (1.0 / LONG_Z), 1.0 / 3);
912 double alpha = 2.0 / leff;
913 double alpha2 = alpha * alpha;
914
915 int qxmax = (int)(8.0 * LONG_X / alpha + 0.5);
916 int qymax = (int)(8.0 * LONG_Y / alpha + 0.5);
917 int qzmax = (int)(8.0 * LONG_Z / alpha + 0.5);
918
919 int nxmax = (int)(2.0 * alpha / LONG_X + 0.5);
920 int nymax = (int)(2.0 * alpha / LONG_Y + 0.5);
921 int nzmax = (int)(2.0 * alpha / LONG_Z + 0.5);
922
923 for(int nx = -qxmax; nx <= qxmax; nx++)
924 for(int ny = -qymax; ny <= qymax; ny++)
925 for(int nz = -qzmax; nz <= qzmax; nz++)
926 {
927 double dx = -nx * (1.0 / LONG_X);
928 double dy = -ny * (1.0 / LONG_Y);
929 double dz = -nz * (1.0 / LONG_Z);
930
931 vector<double> dxyz(dx, dy, dz);
932
933 double r2 = dx * dx + dy * dy + dz * dz;
934 double r = sqrt(r2);
935
936 double rinv = (r > 0) ? 1.0 / r : 0.0;
937 double r2inv = rinv * rinv;
938 double r3inv = r2inv * rinv;
939 double r4inv = r2inv * r2inv;
940 double r7inv = r3inv * r4inv;
941 double r11inv = r4inv * r7inv;
942
943 if(nx != 0 || ny != 0 || nz != 0)
944 {
945 double ar = alpha * r;
946 double ar2 = ar * ar;
947 double ar3 = ar * ar * ar;
948 double ar5 = ar3 * ar * ar;
949 double ar7 = ar5 * ar * ar;
950 double ar9 = ar7 * ar * ar;
951 double ar11 = ar9 * ar * ar;
952 double ar13 = ar11 * ar * ar;
953 double ar15 = ar13 * ar * ar;
954 double ar17 = ar15 * ar * ar;
955 double ar19 = ar17 * ar * ar;
956
957 double xir2 = pow(dx * rinv, 2);
958 double xir4 = pow(dx * rinv, 4);
959 double xir6 = pow(dx * rinv, 6);
960 double xir8 = pow(dx * rinv, 8);
961 double xir10 = pow(dx * rinv, 10);
962
963 double yir2 = pow(dy * rinv, 2);
964 double yir4 = pow(dy * rinv, 4);
965 double zir2 = pow(dz * rinv, 2);
966
967 P10.XXXXXXXXXX +=
968 ((1786050.0 * ar + 1190700.0 * ar3 + 476280.0 * ar5 + 136080.0 * ar7 + 30240 * ar9) -
969 xir2 * (98232750.0 * ar + 65488500.0 * ar3 + 26195400.0 * ar5 + 7484400.0 * ar7 + 1663200.0 * ar9 + 302400 * ar11) +
970 xir4 * (851350500.0 * ar + 567567000.0 * ar3 + 227026800.0 * ar5 + 64864800.0 * ar7 + 14414400.0 * ar9 +
971 2620800.0 * ar11 + 403200 * ar13) -
972 xir6 * (2554051500.0 * ar + 1702701000.0 * ar3 + 681080400.0 * ar5 + 194594400.0 * ar7 + 43243200.0 * ar9 +
973 7862400.0 * ar11 + 1209600.0 * ar13 + 161280 * ar15) +
974 xir8 * (3101348250.0 * ar + 2067565500.0 * ar3 + 827026200.0 * ar5 + 236293200.0 * ar7 + 52509600.0 * ar9 +
975 9547200.0 * ar11 + 1468800.0 * ar13 + 195840.0 * ar15 + 23040 * ar17) -
976 xir10 * (1309458150 * ar + 872972100 * ar3 + 349188840 * ar5 + 99768240 * ar7 + 22170720 * ar9 + 4031040 * ar11 +
977 620160 * ar13 + 82688 * ar15 + 9728 * ar17 + 1024.0 * ar19)) /
978 sqrt(M_PI) * exp(-ar2) * r11inv +
979 erfc(ar) *
980 (893025.0 - 49116375.0 * xir2 + 425675250.0 * xir4 - 1277025750 * xir6 + 1550674125.0 * xir8 -
981 654729075.0 * xir10) *
982 r11inv;
983
984 P10.XXXXXXXXYY +=
985 ((198450.0 * ar + 132300.0 * ar3 + 52920.0 * ar5 + 15120.0 * ar7 + 3360.0 * ar9) -
986 xir2 * (8731800.0 * ar + 5821200.0 * ar3 + 2328480.0 * ar5 + 665280.0 * ar7 + 147840.0 * ar9 + 26880.0 * ar11) +
987 xir4 * (56756700.0 * ar + 37837800.0 * ar3 + 15135120.0 * ar5 + 4324320.0 * ar7 + 960960.0 * ar9 + 174720.0 * ar11 +
988 26880.0 * ar13) -
989 xir6 * (113513400.0 * ar + 75675600.0 * ar3 + 30270240.0 * ar5 + 8648640.0 * ar7 + 1921920.0 * ar9 +
990 349440.0 * ar11 + 53760.0 * ar13 + 7168.0 * ar15) +
991 xir8 * (68918850.0 * ar + 45945900.0 * ar3 + 18378360.0 * ar5 + 5250960.0 * ar7 + 1166880.0 * ar9 +
992 212160.0 * ar11 + 32640.0 * ar13 + 4352.0 * ar15 + 512.0 * ar17) -
993 yir2 * (2182950.0 * ar + 1455300.0 * ar3 + 582120.0 * ar5 + 166320.0 * ar7 + 36960.0 * ar9 + 6720.0 * ar11) +
994 xir2 * yir2 *
995 (113513400.0 * ar + 75675600.0 * ar3 + 30270240.0 * ar5 + 8648640.0 * ar7 + 1921920.0 * ar9 + 349440.0 * ar11 +
996 53760.0 * ar13) -
997 xir4 * yir2 *
998 (851350500.0 * ar + 567567000.0 * ar3 + 227026800.0 * ar5 + 64864800.0 * ar7 + 14414400.0 * ar9 +
999 2620800.0 * ar11 + 403200.0 * ar13 + 53760.0 * ar15) +
1000 xir6 * yir2 *
1001 (1929727800.0 * ar + 1286485200.0 * ar3 + 514594080.0 * ar5 + 147026880.0 * ar7 + 32672640.0 * ar9 +
1002 5940480.0 * ar11 + 913920.0 * ar13 + 121856.0 * ar15 + 14336.0 * ar17) -
1003 xir8 * yir2 *
1004 (1309458150.0 * ar + 872972100.0 * ar3 + 349188840.0 * ar5 + 99768240 * ar7 + 22170720 * ar9 + 4031040 * ar11 +
1005 620160 * ar13 + 82688 * ar15 + 9728 * ar17 + 1024.0 * ar19)) /
1006 sqrt(M_PI) * exp(-ar2) * r11inv +
1007 erfc(ar) *
1008 (99225.0 - 4365900.0 * xir2 + 28378350.0 * xir4 - 56756700 * xir6 + 34459425.0 * xir8 - 1091475.0 * yir2 +
1009 56756700.0 * xir2 * yir2 - 425675250.0 * xir4 * yir2 + 964863900.0 * xir6 * yir2 - 654729075.0 * xir8 * yir2) *
1010 r11inv;
1011
1012 P10.XXXXXXYYYY +=
1013 ((85050.0 * ar + 56700.0 * ar3 + 22680.0 * ar5 + 6480.0 * ar7 + 1440.0 * ar9) -
1014 xir2 * (2806650.0 * ar + 1871100.0 * ar3 + 748440.0 * ar5 + 213840.0 * ar7 + 47520.0 * ar9 + 8640.0 * ar11) +
1015 xir4 * (12162150.0 * ar + 8108100.0 * ar3 + 3243240.0 * ar5 + 926640.0 * ar7 + 205920.0 * ar9 + 37440.0 * ar11 +
1016 5760.0 * ar13) -
1017 xir6 * (12162150.0 * ar + 8108100.0 * ar3 + 3243240.0 * ar5 + 926640.0 * ar7 + 205920.0 * ar9 + 37440.0 * ar11 +
1018 5760.0 * ar13 + 768.0 * ar15) -
1019 yir2 * (1871100.0 * ar + 1247400.0 * ar3 + 498960.0 * ar5 + 142560.0 * ar7 + 31680.0 * ar9 + 5760.0 * ar11) +
1020 xir2 * yir2 *
1021 (72972900.0 * ar + 48648600 * ar3 + 19459440.0 * ar5 + 5559840.0 * ar7 + 1235520.0 * ar9 + 224640.0 * ar11 +
1022 34560.0 * ar13) -
1023 xir4 * yir2 *
1024 (364864500.0 * ar + 243243000.0 * ar3 + 97297200.0 * ar5 + 27799200.0 * ar7 + 6177600.0 * ar9 +
1025 1123200.0 * ar11 + 172800.0 * ar13 + 23040.0 * ar15) +
1026 xir6 * yir2 *
1027 (413513100.0 * ar + 275675400.0 * ar3 + 110270160.0 * ar5 + 31505760.0 * ar7 + 7001280.0 * ar9 +
1028 1272960.0 * ar11 + 195840.0 * ar13 + 26112.0 * ar15 + 3072.0 * ar17) +
1029 yir4 * (4054050.0 * ar + 2702700.0 * ar3 + 1081080.0 * ar5 + 308880.0 * ar7 + 68640.0 * ar9 + 12480.0 * ar11 +
1030 1920 * ar13) -
1031 xir2 * yir4 *
1032 (182432250.0 * ar + 121621500.0 * ar3 + 48648600.0 * ar5 + 13899600 * ar7 + 3088800 * ar9 + 561600 * ar11 +
1033 86400 * ar13 + 11520 * ar15) +
1034 xir4 * yir4 *
1035 (1033782750.0 * ar + 689188500.0 * ar3 + 275675400.0 * ar5 + 78764400 * ar7 + 17503200 * ar9 + 3182400 * ar11 +
1036 489600 * ar13 + 65280 * ar15 + 7680 * ar17) -
1037 xir6 * yir4 *
1038 (1309458150.0 * ar + 872972100.0 * ar3 + 349188840.0 * ar5 + 99768240.0 * ar7 + 22170720.0 * ar9 +
1039 4031040.0 * ar11 + 620160.0 * ar13 + 82688.0 * ar15 + 9728.0 * ar17 + 1024.0 * ar19)) /
1040 sqrt(M_PI) * exp(-ar2) * r11inv +
1041 erfc(ar) *
1042 (42525.0 - 1403325 * xir2 + 6081075.0 * xir4 - 6081075 * xir6 - 935550.0 * yir2 + 36486450.0 * xir2 * yir2 -
1043 182432250.0 * xir4 * yir2 + 206756550.0 * xir6 * yir2 + 2027025.0 * yir4 - 91216125.0 * xir2 * yir4 +
1044 516891375.0 * xir4 * yir4 - 654729075.0 * xir6 * yir4) *
1045 r11inv;
1046
1047 P10.XXXXXXYYZZ +=
1048 ((28350 * ar + 18900 * ar3 + 7560 * ar5 + 2160 * ar7 + 480 * ar9) -
1049 xir2 * (935550 * ar + 623700 * ar3 + 249480 * ar5 + 71280 * ar7 + 15840 * ar9 + 2880 * ar11) +
1050 xir4 * (4054050 * ar + 2702700 * ar3 + 1081080 * ar5 + 308880 * ar7 + 68640 * ar9 + 12480 * ar11 + 1920 * ar13) -
1051 xir6 * (4054050 * ar + 2702700 * ar3 + 1081080 * ar5 + 308880 * ar7 + 68640 * ar9 + 12480 * ar11 + 1920 * ar13 +
1052 256 * ar15) -
1053 (yir2 + zir2) * (311850 * ar + 207900 * ar3 + 83160 * ar5 + 23760 * ar7 + 5280 * ar9 + 960 * ar11) +
1054 (xir2 * yir2 + xir2 * zir2) *
1055 (12162150 * ar + 8108100 * ar3 + 3243240 * ar5 + 926640 * ar7 + 205920 * ar9 + 37440 * ar11 + 5760 * ar13) -
1056 (xir4 * yir2 + xir4 * zir2) * (60810750 * ar + 40540500 * ar3 + 16216200 * ar5 + 4633200 * ar7 + 1029600 * ar9 +
1057 187200 * ar11 + 28800 * ar13 + 3840 * ar15) +
1058 (xir6 * yir2 + xir6 * zir2) * (68918850 * ar + 45945900 * ar3 + 18378360 * ar5 + 5250960 * ar7 + 1166880 * ar9 +
1059 212160 * ar11 + 32640 * ar13 + 4352 * ar15 + 512 * ar17) +
1060 yir2 * zir2 *
1061 (4054050 * ar + 2702700 * ar3 + 1081080.0 * ar5 + 308880.0 * ar7 + 68640.0 * ar9 + 12480.0 * ar11 +
1062 1920 * ar13) -
1063 xir2 * yir2 * zir2 *
1064 (182432250.0 * ar + 121621500.0 * ar3 + 48648600.0 * ar5 + 13899600 * ar7 + 3088800 * ar9 + 561600 * ar11 +
1065 86400 * ar13 + 11520 * ar15) +
1066 xir4 * yir2 * zir2 *
1067 (1033782750.0 * ar + 689188500.0 * ar3 + 275675400.0 * ar5 + 78764400 * ar7 + 17503200 * ar9 + 3182400 * ar11 +
1068 489600 * ar13 + 65280 * ar15 + 7680 * ar17) -
1069 xir6 * yir2 * zir2 *
1070 (1309458150.0 * ar + 872972100.0 * ar3 + 349188840.0 * ar5 + 99768240.0 * ar7 + 22170720.0 * ar9 +
1071 4031040.0 * ar11 + 620160.0 * ar13 + 82688.0 * ar15 + 9728.0 * ar17 + 1024.0 * ar19)) /
1072 sqrt(M_PI) * exp(-ar2) * r11inv +
1073 erfc(ar) *
1074 (14175.0 - 467775 * xir2 + 2027025 * xir4 - 2027025 * xir6 - 155925.0 * (yir2 + zir2) +
1075 6081075.0 * xir2 * (yir2 + zir2) - 30405375.0 * xir4 * (yir2 + zir2) + 34459425.0 * xir6 * (yir2 + zir2) +
1076 2027025.0 * yir2 * zir2 - 91216125.0 * xir2 * yir2 * zir2 + 516891375.0 * xir4 * yir2 * zir2 -
1077 654729075.0 * xir6 * yir2 * zir2) *
1078 r11inv;
1079
1080 P10.XXXXYYYYZZ +=
1081 ((17010 * ar + 11340 * ar3 + 4536 * ar5 + 1296 * ar7 + 288 * ar9) -
1082 (xir2 + yir2) * (374220 * ar + 249480 * ar3 + 99792 * ar5 + 28512 * ar7 + 6336 * ar9 + 1152 * ar11) +
1083 (xir4 + yir4) * (810810 * ar + 540540 * ar3 + 216216 * ar5 + 61776 * ar7 + 13728 * ar9 + 2496 * ar11 + 384 * ar13) +
1084 xir2 * yir2 *
1085 (9729720 * ar + 6486480 * ar3 + 2594592 * ar5 + 741312 * ar7 + 164736 * ar9 + 29952 * ar11 + 4608 * ar13) -
1086 (xir4 * yir2 + xir2 * yir4) * (24324300 * ar + 16216200 * ar3 + 6486480 * ar5 + 1853280 * ar7 + 411840 * ar9 +
1087 74880 * ar11 + 11520 * ar13 + 1536 * ar15) +
1088 xir4 * yir4 *
1089 (68918850 * ar + 45945900 * ar3 + 18378360 * ar5 + 5250960 * ar7 + 1166880 * ar9 + 212160 * ar11 +
1090 32640 * ar13 + 4352 * ar15 + 512.0 * ar17) -
1091 zir2 * (187110 * ar + 124740 * ar3 + 49896 * ar5 + 14256 * ar7 + 3168 * ar9 + 576 * ar11) +
1092 (xir2 + yir2) * zir2 *
1093 (4864860 * ar + 3243240 * ar3 + 1297296 * ar5 + 370656 * ar7 + 82368 * ar9 + 14976 * ar11 + 2304 * ar13) -
1094 (xir4 + yir4) * zir2 *
1095 (12162150 * ar + 8108100 * ar3 + 3243240 * ar5 + 926640 * ar7 + 205920 * ar9 + 37440 * ar11 + 5760 * ar13 +
1096 768 * ar15) -
1097 xir2 * yir2 * zir2 *
1098 (145945800 * ar + 97297200 * ar3 + 38918880 * ar5 + 11119680 * ar7 + 2471040 * ar9 + 449280 * ar11 +
1099 69120 * ar13 + 9216 * ar15) +
1100 (xir4 * yir2 + xir2 * yir4) * zir2 *
1101 (413513100 * ar + 275675400 * ar3 + 110270160 * ar5 + 31505760 * ar7 + 7001280 * ar9 + 1272960 * ar11 +
1102 195840 * ar13 + 26112 * ar15 + 3072 * ar17) -
1103 xir4 * yir4 * zir2 *
1104 (1309458150.0 * ar + 872972100.0 * ar3 + 349188840.0 * ar5 + 99768240.0 * ar7 + 22170720.0 * ar9 +
1105 4031040.0 * ar11 + 620160.0 * ar13 + 82688.0 * ar15 + 9728.0 * ar17 + 1024.0 * ar19)) /
1106 sqrt(M_PI) * exp(-ar2) * r11inv +
1107 erfc(ar) *
1108 (8505.0 - 187110 * (xir2 + yir2) + 405405 * (xir4 + yir4) + 4864860 * xir2 * yir2 -
1109 12162150 * (xir4 * yir2 + xir2 * yir4) + 34459425 * xir4 * yir4 - 93555 * zir2 +
1110 2432430 * (xir2 + yir2) * zir2 - 6081075 * (xir4 + yir4) * zir2 - 72972900 * xir2 * yir2 * zir2 +
1111 206756550 * (xir4 * yir2 + xir2 * yir4) * zir2 - 654729075.0 * xir4 * yir4 * zir2) *
1112 r11inv;
1113 }
1114 else
1115 {
1116 /* we add the 1/r term here to the (0|0|0) entry, followed by differentiation, and the limit r->0 to obtain accurate
1117 * results at the origin
1118 */
1119
1120 /* Note, for small r:
1121 *
1122 * [1/- erfc(a r)]/r = 2 a/sqrt(pi) * [ 1 - (a r)^2/3 + (a r)^4 / 10 - (a r)^6 / 42 + (a r)^8 / 216 - ...]
1123 *
1124 */
1125
1126 P10.XXXXXXXXXX += (-60480.0) * pow(alpha, 11) / (11.0 * sqrt(M_PI));
1127 P10.XXXXXXXXYY += (-6720.0) * pow(alpha, 11) / (11.0 * sqrt(M_PI));
1128 P10.XXXXXXYYYY += (-2880.0) * pow(alpha, 11) / (11.0 * sqrt(M_PI));
1129 P10.XXXXXXYYZZ += (-960.0) * pow(alpha, 11) / (11.0 * sqrt(M_PI));
1130 P10.XXXXYYYYZZ += (-576.0) * pow(alpha, 11) / (11.0 * sqrt(M_PI));
1131 }
1132 }
1133
1134 for(int nx = -nxmax; nx <= nxmax; nx++)
1135 for(int ny = -nymax; ny <= nymax; ny++)
1136 for(int nz = -nzmax; nz <= nzmax; nz++)
1137 {
1138 if(nx != 0 || ny != 0 || nz != 0)
1139 {
1140 double kx = (2.0 * M_PI * LONG_X) * nx;
1141 double ky = (2.0 * M_PI * LONG_Y) * ny;
1142 double kz = (2.0 * M_PI * LONG_Z) * nz;
1143 double k2 = kx * kx + ky * ky + kz * kz;
1144
1145 double val = 4.0 * M_PI * (LONG_X * LONG_Y * LONG_Z) / k2 * exp(-k2 / (4.0 * alpha2));
1146
1147 P10.XXXXXXXXXX += val * pow(kx, 10);
1148 P10.XXXXXXXXYY += val * pow(kx, 8) * pow(ky, 2);
1149 P10.XXXXXXYYYY += val * pow(kx, 6) * pow(ky, 4);
1150 P10.XXXXXXYYZZ += val * pow(kx, 6) * pow(ky, 2) * pow(kz, 2);
1151 P10.XXXXYYYYZZ += val * pow(kx, 4) * pow(ky, 4) * pow(kz, 2);
1152 }
1153 }
1154
1155 return P10;
1156}
1157
1158#endif
All
global_data_all_processes All
Definition: main.cc:40
ewald::ewald_D3
symtensor3< double > ewald_D3(double x, double y, double z)
Definition: ewald.cc:1312
ewald::ewald_P6
ewaldtensor6< double > ewald_P6(void)
ewald::ewald_D0
double ewald_D0(double x, double y, double z)
This function computes the potential correction term by means of Ewald summation.
Definition: ewald.cc:705
ewald::ewald_D4
symtensor4< double > ewald_D4(double x, double y, double z)
Definition: ewald.cc:1550
ewald::ewald_corr
void ewald_corr(double dx, double dy, double dz, enum interpolate_options, ewald_data &fper)
ewald::ewald_P10
ewaldtensor10< double > ewald_P10(void)
ewald::POINTMASS
@ POINTMASS
Definition: ewald.h:67
ewald::MULTIPOLES
@ MULTIPOLES
Definition: ewald.h:68
ewald::ewald_P8
ewaldtensor8< double > ewald_P8(void)
ewald::ewald_D5
symtensor5< double > ewald_D5(double x, double y, double z)
Definition: ewald.cc:1694
ewald::ewald_D1
vector< double > ewald_D1(double x, double y, double z)
This function computes the force correction term (difference between full force of infinite lattice a...
Definition: ewald.cc:898
ewald::ewald_corr_exact
void ewald_corr_exact(double dx, double dy, double dz, enum interpolate_options flag, ewald_data &fper)
ewald::ewald_D2
symtensor2< double > ewald_D2(double x, double y, double z)
Definition: ewald.cc:1090
ewald::ewald_gridlookup
void ewald_gridlookup(const MyIntPosType *p_intpos, const MyIntPosType *target_intpos, enum interpolate_options flag, ewald_data &fper)
Definition: ewald.cc:355
intposconvert::pos_to_signedintpos
MySignedIntPosType pos_to_signedintpos(T posdiff)
Definition: intposconvert.h:325
symtensor2::norm
double norm(void)
Definition: symtensors.h:287
symtensor3::norm
double norm(void)
Definition: symtensors.h:417
symtensor4::norm
double norm(void)
Definition: symtensors.h:504
symtensor5::norm
double norm(void)
Definition: symtensors.h:597
vector::norm
double norm(void)
Definition: symtensors.h:189
M_PI
#define M_PI
Definition: constants.h:56
LONG_Y
#define LONG_Y
Definition: dtypes.h:369
MyIntPosType
uint32_t MyIntPosType
Definition: dtypes.h:35
MySignedIntPosType
int32_t MySignedIntPosType
Definition: dtypes.h:36
LONG_X
#define LONG_X
Definition: dtypes.h:361
LONG_Z
#define LONG_Z
Definition: dtypes.h:377
ENX
#define ENX
Definition: ewald.h:45
ENZ
#define ENZ
Definition: ewald.h:47
ENY
#define ENY
Definition: ewald.h:46
Terminate
#define Terminate(...)
Definition: macros.h:15
half_float::detail::exp
expr exp(half arg)
Definition: half.hpp:2724
half_float::detail::fabs
half fabs(half arg)
Definition: half.hpp:2627
half_float::detail::pow
expr pow(half base, half exp)
Definition: half.hpp:2803
half_float::detail::erfc
expr erfc(half arg)
Definition: half.hpp:2925
half_float::detail::sqrt
expr sqrt(half arg)
Definition: half.hpp:2777
ewald_data
Definition: gravtree.h:184
ewald_data::D1phi
vector< MyReal > D1phi
Definition: gravtree.h:186
ewald_data::D3phi
symtensor3< MyReal > D3phi
Definition: gravtree.h:188
ewald_data::D0phi
MyReal D0phi
Definition: gravtree.h:185
ewald_data::D2phi
symtensor2< MyReal > D2phi
Definition: gravtree.h:187
global_data_all_processes::BoxSize
double BoxSize
Definition: allvars.h:144
rYZZZZ
#define rYZZZZ
Definition: symtensor_indices.h:339
rXXYZZ
#define rXXYZZ
Definition: symtensor_indices.h:195
dYZZ
#define dYZZ
Definition: symtensor_indices.h:54
rYYZZZ
#define rYYZZZ
Definition: symtensor_indices.h:312
sXXXY
#define sXXXY
Definition: symtensor_indices.h:70
rXXZZZ
#define rXXZZZ
Definition: symtensor_indices.h:204
dXYY
#define dXYY
Definition: symtensor_indices.h:37
rXXXYY
#define rXXXYY
Definition: symtensor_indices.h:182
sXXYZ
#define sXXYZ
Definition: symtensor_indices.h:75
sXYYZ
#define sXYYZ
Definition: symtensor_indices.h:87
dZZZ
#define dZZZ
Definition: symtensor_indices.h:66
rXXXYZ
#define rXXXYZ
Definition: symtensor_indices.h:183
dYYY
#define dYYY
Definition: symtensor_indices.h:49
dXYZ
#define dXYZ
Definition: symtensor_indices.h:38
sYYYZ
#define sYYYZ
Definition: symtensor_indices.h:123
rXXYYY
#define rXXYYY
Definition: symtensor_indices.h:191
rXXXXZ
#define rXXXXZ
Definition: symtensor_indices.h:180
sXZZZ
#define sXZZZ
Definition: symtensor_indices.h:103
rXXXXY
#define rXXXXY
Definition: symtensor_indices.h:179
rYYYYZ
#define rYYYYZ
Definition: symtensor_indices.h:300
rXYYYZ
#define rXYYYZ
Definition: symtensor_indices.h:219
qYZ
#define qYZ
Definition: symtensor_indices.h:26
rXYYYY
#define rXYYYY
Definition: symtensor_indices.h:218
sXYYY
#define sXYYY
Definition: symtensor_indices.h:86
dXXY
#define dXXY
Definition: symtensor_indices.h:33
dXXX
#define dXXX
Definition: symtensor_indices.h:32
dXXZ
#define dXXZ
Definition: symtensor_indices.h:34
rYYYZZ
#define rYYYZZ
Definition: symtensor_indices.h:303
dYYZ
#define dYYZ
Definition: symtensor_indices.h:50
rXXXXX
#define rXXXXX
Definition: symtensor_indices.h:178
rXXYYZ
#define rXXYYZ
Definition: symtensor_indices.h:192
sXYZZ
#define sXYZZ
Definition: symtensor_indices.h:91
rZZZZZ
#define rZZZZZ
Definition: symtensor_indices.h:420
rXXXZZ
#define rXXXZZ
Definition: symtensor_indices.h:186
rXYYZZ
#define rXYYZZ
Definition: symtensor_indices.h:222
sXXXZ
#define sXXXZ
Definition: symtensor_indices.h:71
sYZZZ
#define sYZZZ
Definition: symtensor_indices.h:139
qXY
#define qXY
Definition: symtensor_indices.h:22
rXYZZZ
#define rXYZZZ
Definition: symtensor_indices.h:231
rYYYYY
#define rYYYYY
Definition: symtensor_indices.h:299
dXZZ
#define dXZZ
Definition: symtensor_indices.h:42
rXZZZZ
#define rXZZZZ
Definition: symtensor_indices.h:258
qXZ
#define qXZ
Definition: symtensor_indices.h:23
init_rng
void init_rng(int thistask)
Definition: system.cc:42
get_random_number
double get_random_number(void)
Definition: system.cc:49