2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2010,2011,2012,2013,2014 by the GROMACS development team.
5 * Copyright (c) 2015,2016,2017,2018,2019 by the GROMACS development team.
6 * Copyright (c) 2019,2020, 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 "pme_error.h"
47 #include "gromacs/commandline/pargs.h"
48 #include "gromacs/ewald/ewald_utils.h"
49 #include "gromacs/ewald/pme.h"
50 #include "gromacs/fft/calcgrid.h"
51 #include "gromacs/fileio/checkpoint.h"
52 #include "gromacs/fileio/tpxio.h"
53 #include "gromacs/gmxlib/network.h"
54 #include "gromacs/math/functions.h"
55 #include "gromacs/math/units.h"
56 #include "gromacs/math/vec.h"
57 #include "gromacs/mdlib/broadcaststructs.h"
58 #include "gromacs/mdtypes/commrec.h"
59 #include "gromacs/mdtypes/inputrec.h"
60 #include "gromacs/mdtypes/md_enums.h"
61 #include "gromacs/mdtypes/state.h"
62 #include "gromacs/random/threefry.h"
63 #include "gromacs/random/uniformintdistribution.h"
64 #include "gromacs/topology/mtop_util.h"
65 #include "gromacs/topology/topology.h"
66 #include "gromacs/utility/arraysize.h"
67 #include "gromacs/utility/fatalerror.h"
68 #include "gromacs/utility/pleasecite.h"
69 #include "gromacs/utility/smalloc.h"
71 /* #define TAKETIME */
74 /* Enum for situations that can occur during log file parsing */
81 eParselogResetProblem,
88 int64_t orig_sim_steps; /* Number of steps to be done in the real simulation */
89 int n_entries; /* Number of entries in arrays */
90 real volume; /* The volume of the box */
91 matrix recipbox; /* The reciprocal box */
92 int natoms; /* The number of atoms in the MD system */
93 real* fac; /* The scaling factor */
94 real* rcoulomb; /* The coulomb radii [0...nr_inputfiles] */
95 real* rvdw; /* The vdW radii */
96 int * nkx, *nky, *nkz; /* Number of k vectors in each spatial dimension */
97 real* fourier_sp; /* Fourierspacing */
98 real* ewald_rtol; /* Real space tolerance for Ewald, determines */
99 /* the real/reciprocal space relative weight */
100 real* ewald_beta; /* Splitting parameter [1/nm] */
101 real fracself; /* fraction of particles for SI error */
102 real q2all; /* sum ( q ^2 ) */
103 real q2allnr; /* nr of charges */
104 int* pme_order; /* Interpolation order for PME (bsplines) */
105 char** fn_out; /* Name of the output tpr file */
106 real* e_dir; /* Direct space part of PME error with these settings */
107 real* e_rec; /* Reciprocal space part of PME error */
108 gmx_bool bTUNE; /* flag for tuning */
112 /* Returns TRUE when atom is charged */
113 static gmx_bool is_charge(real charge)
115 return charge * charge > GMX_REAL_EPS;
119 /* calculate charge density */
120 static void calc_q2all(const gmx_mtop_t* mtop, /* molecular topology */
124 real q2_all = 0; /* Sum of squared charges */
125 int nrq_mol; /* Number of charges in a single molecule */
126 int nrq_all; /* Total number of charges in the MD system */
130 fprintf(stderr, "\nCharge density:\n");
132 q2_all = 0.0; /* total q squared */
133 nrq_all = 0; /* total number of charges in the system */
134 for (const gmx_molblock_t& molblock : mtop->molblock) /* Loop over molecule types */
136 q2_mol = 0.0; /* q squared value of this molecule */
137 nrq_mol = 0; /* number of charges this molecule carries */
138 const gmx_moltype_t& molecule = mtop->moltype[molblock.type];
139 for (int i = 0; i < molecule.atoms.nr; i++)
141 qi = molecule.atoms.atom[i].q;
142 /* Is this charge worth to be considered? */
149 /* Multiply with the number of molecules present of this type and add */
150 q2_all += q2_mol * molblock.nmol;
151 nrq_all += nrq_mol * molblock.nmol;
154 "Molecule %2d (%5d atoms) q2_mol=%10.3e nr.mol.charges=%5d (%6dx) q2_all=%10.3e "
156 imol, molecule.atoms.nr, q2_mol, nrq_mol, molblock.nmol, q2_all, nrq_all);
165 /* Estimate the direct space part error of the SPME Ewald sum */
166 static real estimate_direct(t_inputinfo* info)
168 real e_dir = 0; /* Error estimate */
169 real beta = 0; /* Splitting parameter (1/nm) */
170 real r_coulomb = 0; /* Cut-off in direct space */
173 beta = info->ewald_beta[0];
174 r_coulomb = info->rcoulomb[0];
176 e_dir = 2.0 * info->q2all * gmx::invsqrt(info->q2allnr * r_coulomb * info->volume);
177 e_dir *= std::exp(-beta * beta * r_coulomb * r_coulomb);
179 return ONE_4PI_EPS0 * e_dir;
184 /* the following 4 functions determine polynomials required for the reciprocal error estimate */
186 static inline real eps_poly1(real m, /* grid coordinate in certain direction */
187 real K, /* grid size in corresponding direction */
188 real n) /* spline interpolation order of the SPME */
191 real nom = 0; /* nominator */
192 real denom = 0; /* denominator */
200 for (i = -SUMORDER; i < 0; i++)
204 nom += std::pow(tmp, -n);
207 for (i = SUMORDER; i > 0; i--)
211 nom += std::pow(tmp, -n);
216 denom = std::pow(tmp, -n) + nom;
221 static inline real eps_poly2(real m, /* grid coordinate in certain direction */
222 real K, /* grid size in corresponding direction */
223 real n) /* spline interpolation order of the SPME */
226 real nom = 0; /* nominator */
227 real denom = 0; /* denominator */
235 for (i = -SUMORDER; i < 0; i++)
239 nom += std::pow(tmp, -2 * n);
242 for (i = SUMORDER; i > 0; i--)
246 nom += std::pow(tmp, -2 * n);
249 for (i = -SUMORDER; i < SUMORDER + 1; i++)
253 denom += std::pow(tmp, -n);
255 tmp = eps_poly1(m, K, n);
256 return nom / denom / denom + tmp * tmp;
259 static inline real eps_poly3(real m, /* grid coordinate in certain direction */
260 real K, /* grid size in corresponding direction */
261 real n) /* spline interpolation order of the SPME */
264 real nom = 0; /* nominator */
265 real denom = 0; /* denominator */
273 for (i = -SUMORDER; i < 0; i++)
277 nom += i * std::pow(tmp, -2 * n);
280 for (i = SUMORDER; i > 0; i--)
284 nom += i * std::pow(tmp, -2 * n);
287 for (i = -SUMORDER; i < SUMORDER + 1; i++)
291 denom += std::pow(tmp, -n);
294 return 2.0 * M_PI * nom / denom / denom;
297 static inline real eps_poly4(real m, /* grid coordinate in certain direction */
298 real K, /* grid size in corresponding direction */
299 real n) /* spline interpolation order of the SPME */
302 real nom = 0; /* nominator */
303 real denom = 0; /* denominator */
311 for (i = -SUMORDER; i < 0; i++)
315 nom += i * i * std::pow(tmp, -2 * n);
318 for (i = SUMORDER; i > 0; i--)
322 nom += i * i * std::pow(tmp, -2 * n);
325 for (i = -SUMORDER; i < SUMORDER + 1; i++)
329 denom += std::pow(tmp, -n);
332 return 4.0 * M_PI * M_PI * nom / denom / denom;
335 static inline real eps_self(real m, /* grid coordinate in certain direction */
336 real K, /* grid size in corresponding direction */
337 rvec rboxv, /* reciprocal box vector */
338 real n, /* spline interpolation order of the SPME */
339 rvec x) /* coordinate of charge */
342 real tmp = 0; /* temporary variables for computations */
343 real tmp1 = 0; /* temporary variables for computations */
344 real tmp2 = 0; /* temporary variables for computations */
345 real rcoord = 0; /* coordinate in certain reciprocal space direction */
346 real nom = 0; /* nominator */
347 real denom = 0; /* denominator */
355 rcoord = iprod(rboxv, x);
358 for (i = -SUMORDER; i < 0; i++)
360 tmp = -std::sin(2.0 * M_PI * i * K * rcoord);
361 tmp1 = 2.0 * M_PI * m / K + 2.0 * M_PI * i;
362 tmp2 = std::pow(tmp1, -n);
363 nom += tmp * tmp2 * i;
367 for (i = SUMORDER; i > 0; i--)
369 tmp = -std::sin(2.0 * M_PI * i * K * rcoord);
370 tmp1 = 2.0 * M_PI * m / K + 2.0 * M_PI * i;
371 tmp2 = std::pow(tmp1, -n);
372 nom += tmp * tmp2 * i;
377 tmp = 2.0 * M_PI * m / K;
381 return 2.0 * M_PI * nom / denom * K;
386 /* The following routine is just a copy from pme.c */
388 static void calc_recipbox(matrix box, matrix recipbox)
390 /* Save some time by assuming upper right part is zero */
392 real tmp = 1.0 / (box[XX][XX] * box[YY][YY] * box[ZZ][ZZ]);
394 recipbox[XX][XX] = box[YY][YY] * box[ZZ][ZZ] * tmp;
395 recipbox[XX][YY] = 0;
396 recipbox[XX][ZZ] = 0;
397 recipbox[YY][XX] = -box[YY][XX] * box[ZZ][ZZ] * tmp;
398 recipbox[YY][YY] = box[XX][XX] * box[ZZ][ZZ] * tmp;
399 recipbox[YY][ZZ] = 0;
400 recipbox[ZZ][XX] = (box[YY][XX] * box[ZZ][YY] - box[YY][YY] * box[ZZ][XX]) * tmp;
401 recipbox[ZZ][YY] = -box[ZZ][YY] * box[XX][XX] * tmp;
402 recipbox[ZZ][ZZ] = box[XX][XX] * box[YY][YY] * tmp;
406 /* Estimate the reciprocal space part error of the SPME Ewald sum. */
407 static real estimate_reciprocal(t_inputinfo* info,
408 rvec x[], /* array of particles */
409 const real q[], /* array of charges */
410 int nr, /* number of charges = size of the charge array */
411 FILE gmx_unused* fp_out,
413 int seed, /* The seed for the random number generator */
414 int* nsamples, /* Return the number of samples used if Monte Carlo
415 * algorithm is used for self energy error estimate */
418 real e_rec = 0; /* reciprocal error estimate */
419 real e_rec1 = 0; /* Error estimate term 1*/
420 real e_rec2 = 0; /* Error estimate term 2*/
421 real e_rec3 = 0; /* Error estimate term 3 */
422 real e_rec3x = 0; /* part of Error estimate term 3 in x */
423 real e_rec3y = 0; /* part of Error estimate term 3 in y */
424 real e_rec3z = 0; /* part of Error estimate term 3 in z */
426 int nx, ny, nz; /* grid coordinates */
427 real q2_all = 0; /* sum of squared charges */
428 rvec gridpx; /* reciprocal grid point in x direction*/
429 rvec gridpxy; /* reciprocal grid point in x and y direction*/
430 rvec gridp; /* complete reciprocal grid point in 3 directions*/
431 rvec tmpvec; /* template to create points from basis vectors */
432 rvec tmpvec2; /* template to create points from basis vectors */
433 real coeff = 0; /* variable to compute coefficients of the error estimate */
434 real coeff2 = 0; /* variable to compute coefficients of the error estimate */
435 real tmp = 0; /* variables to compute different factors from vectors */
440 int* numbers = nullptr;
442 /* Index variables for parallel work distribution */
443 int startglobal, stopglobal;
444 int startlocal, stoplocal;
453 GMX_RELEASE_ASSERT(q != nullptr, "Must have charges");
457 seed = static_cast<int>(gmx::makeRandomSeed());
459 fprintf(stderr, "Using random seed %d.\n", seed);
461 gmx::DefaultRandomEngine rng(seed);
462 gmx::UniformIntDistribution<int> dist(0, nr - 1);
470 for (i = 0; i < nr; i++)
472 q2_all += q[i] * q[i];
475 /* Calculate indices for work distribution */
476 startglobal = -info->nkx[0] / 2;
477 stopglobal = info->nkx[0] / 2;
478 xtot = stopglobal * 2 + 1;
481 x_per_core = static_cast<int>(std::ceil(static_cast<real>(xtot) / cr->nnodes));
482 startlocal = startglobal + x_per_core * cr->nodeid;
483 stoplocal = startlocal + x_per_core - 1;
484 if (stoplocal > stopglobal)
486 stoplocal = stopglobal;
491 startlocal = startglobal;
492 stoplocal = stopglobal;
508 fprintf(stderr, "Calculating reciprocal error part 1 ...");
511 for (nx = startlocal; nx <= stoplocal; nx++)
513 svmul(nx, info->recipbox[XX], gridpx);
514 for (ny = -info->nky[0] / 2; ny < info->nky[0] / 2 + 1; ny++)
516 svmul(ny, info->recipbox[YY], tmpvec);
517 rvec_add(gridpx, tmpvec, gridpxy);
518 for (nz = -info->nkz[0] / 2; nz < info->nkz[0] / 2 + 1; nz++)
520 if (0 == nx && 0 == ny && 0 == nz)
524 svmul(nz, info->recipbox[ZZ], tmpvec);
525 rvec_add(gridpxy, tmpvec, gridp);
527 coeff = std::exp(-1.0 * M_PI * M_PI * tmp / info->ewald_beta[0] / info->ewald_beta[0]);
528 coeff /= 2.0 * M_PI * info->volume * tmp;
532 tmp = eps_poly2(nx, info->nkx[0], info->pme_order[0]);
533 tmp += eps_poly2(ny, info->nkx[0], info->pme_order[0]);
534 tmp += eps_poly2(nz, info->nkx[0], info->pme_order[0]);
536 tmp1 = eps_poly1(nx, info->nkx[0], info->pme_order[0]);
537 tmp2 = eps_poly1(ny, info->nky[0], info->pme_order[0]);
539 tmp += 2.0 * tmp1 * tmp2;
541 tmp1 = eps_poly1(nz, info->nkz[0], info->pme_order[0]);
542 tmp2 = eps_poly1(ny, info->nky[0], info->pme_order[0]);
544 tmp += 2.0 * tmp1 * tmp2;
546 tmp1 = eps_poly1(nz, info->nkz[0], info->pme_order[0]);
547 tmp2 = eps_poly1(nx, info->nkx[0], info->pme_order[0]);
549 tmp += 2.0 * tmp1 * tmp2;
551 tmp1 = eps_poly1(nx, info->nkx[0], info->pme_order[0]);
552 tmp1 += eps_poly1(ny, info->nky[0], info->pme_order[0]);
553 tmp1 += eps_poly1(nz, info->nkz[0], info->pme_order[0]);
557 e_rec1 += 32.0 * M_PI * M_PI * coeff * coeff * coeff2 * tmp * q2_all * q2_all / nr;
559 tmp1 = eps_poly3(nx, info->nkx[0], info->pme_order[0]);
560 tmp1 *= info->nkx[0];
561 tmp2 = iprod(gridp, info->recipbox[XX]);
565 tmp1 = eps_poly3(ny, info->nky[0], info->pme_order[0]);
566 tmp1 *= info->nky[0];
567 tmp2 = iprod(gridp, info->recipbox[YY]);
571 tmp1 = eps_poly3(nz, info->nkz[0], info->pme_order[0]);
572 tmp1 *= info->nkz[0];
573 tmp2 = iprod(gridp, info->recipbox[ZZ]);
579 tmp1 = eps_poly4(nx, info->nkx[0], info->pme_order[0]);
580 tmp1 *= norm2(info->recipbox[XX]);
581 tmp1 *= info->nkx[0] * info->nkx[0];
585 tmp1 = eps_poly4(ny, info->nky[0], info->pme_order[0]);
586 tmp1 *= norm2(info->recipbox[YY]);
587 tmp1 *= info->nky[0] * info->nky[0];
591 tmp1 = eps_poly4(nz, info->nkz[0], info->pme_order[0]);
592 tmp1 *= norm2(info->recipbox[ZZ]);
593 tmp1 *= info->nkz[0] * info->nkz[0];
597 e_rec2 += 4.0 * coeff * coeff * tmp * q2_all * q2_all / nr;
602 fprintf(stderr, "\rCalculating reciprocal error part 1 ... %3.0f%%",
603 100.0 * (nx - startlocal + 1) / (x_per_core));
610 fprintf(stderr, "\n");
613 /* Use just a fraction of all charges to estimate the self energy error term? */
614 bFraction = (info->fracself > 0.0) && (info->fracself < 1.0);
618 /* Here xtot is the number of samples taken for the Monte Carlo calculation
619 * of the average of term IV of equation 35 in Wang2010. Round up to a
620 * number of samples that is divisible by the number of nodes */
621 x_per_core = static_cast<int>(std::ceil(info->fracself * nr / cr->nnodes));
622 xtot = x_per_core * cr->nnodes;
626 /* In this case we use all nr particle positions */
628 x_per_core = static_cast<int>(std::ceil(static_cast<real>(xtot) / cr->nnodes));
631 startlocal = x_per_core * cr->nodeid;
632 stoplocal = std::min(startlocal + x_per_core, xtot); /* min needed if xtot == nr */
636 /* Make shure we get identical results in serial and parallel. Therefore,
637 * take the sample indices from a single, global random number array that
638 * is constructed on the master node and that only depends on the seed */
642 for (i = 0; i < xtot; i++)
644 numbers[i] = dist(rng); // [0,nr-1]
647 /* Broadcast the random number array to the other nodes */
650 nblock_bc(cr, xtot, numbers);
653 if (bVerbose && MASTER(cr))
655 fprintf(stdout, "Using %d sample%s to approximate the self interaction error term",
656 xtot, xtot == 1 ? "" : "s");
659 fprintf(stdout, " (%d sample%s per rank)", x_per_core, x_per_core == 1 ? "" : "s");
661 fprintf(stdout, ".\n");
665 /* Return the number of positions used for the Monte Carlo algorithm */
668 for (i = startlocal; i < stoplocal; i++)
676 /* Randomly pick a charge */
681 /* Use all charges */
685 /* for(nx=startlocal; nx<=stoplocal; nx++)*/
686 for (nx = -info->nkx[0] / 2; nx < info->nkx[0] / 2 + 1; nx++)
688 svmul(nx, info->recipbox[XX], gridpx);
689 for (ny = -info->nky[0] / 2; ny < info->nky[0] / 2 + 1; ny++)
691 svmul(ny, info->recipbox[YY], tmpvec);
692 rvec_add(gridpx, tmpvec, gridpxy);
693 for (nz = -info->nkz[0] / 2; nz < info->nkz[0] / 2 + 1; nz++)
696 if (0 == nx && 0 == ny && 0 == nz)
701 svmul(nz, info->recipbox[ZZ], tmpvec);
702 rvec_add(gridpxy, tmpvec, gridp);
704 coeff = std::exp(-1.0 * M_PI * M_PI * tmp / info->ewald_beta[0] / info->ewald_beta[0]);
707 * eps_self(nx, info->nkx[0], info->recipbox[XX], info->pme_order[0], x[ci]);
709 * eps_self(ny, info->nky[0], info->recipbox[YY], info->pme_order[0], x[ci]);
711 * eps_self(nz, info->nkz[0], info->recipbox[ZZ], info->pme_order[0], x[ci]);
718 svmul(e_rec3x, info->recipbox[XX], tmpvec);
719 rvec_inc(tmpvec2, tmpvec);
720 svmul(e_rec3y, info->recipbox[YY], tmpvec);
721 rvec_inc(tmpvec2, tmpvec);
722 svmul(e_rec3z, info->recipbox[ZZ], tmpvec);
723 rvec_inc(tmpvec2, tmpvec);
725 e_rec3 += q[ci] * q[ci] * q[ci] * q[ci] * norm2(tmpvec2)
726 / (xtot * M_PI * info->volume * M_PI * info->volume);
729 fprintf(stderr, "\rCalculating reciprocal error part 2 ... %3.0f%%", 100.0 * (i + 1) / stoplocal);
736 fprintf(stderr, "\n");
743 t1 = MPI_Wtime() - t0;
744 fprintf(fp_out, "Recip. err. est. took : %lf s\n", t1);
752 fprintf(stderr, "Rank %3d: nx=[%3d...%3d] e_rec3=%e\n", cr->nodeid, startlocal, stoplocal, e_rec3);
758 gmx_sum(1, &e_rec1, cr);
759 gmx_sum(1, &e_rec2, cr);
760 gmx_sum(1, &e_rec3, cr);
763 /* e_rec1*=8.0 * q2_all / info->volume / info->volume / nr ;
764 e_rec2*= q2_all / M_PI / M_PI / info->volume / info->volume / nr ;
765 e_rec3/= M_PI * M_PI * info->volume * info->volume * nr ;
767 e_rec = std::sqrt(e_rec1 + e_rec2 + e_rec3);
770 return ONE_4PI_EPS0 * e_rec;
774 /* Allocate memory for the inputinfo struct: */
775 static void create_info(t_inputinfo* info)
777 snew(info->fac, info->n_entries);
778 snew(info->rcoulomb, info->n_entries);
779 snew(info->rvdw, info->n_entries);
780 snew(info->nkx, info->n_entries);
781 snew(info->nky, info->n_entries);
782 snew(info->nkz, info->n_entries);
783 snew(info->fourier_sp, info->n_entries);
784 snew(info->ewald_rtol, info->n_entries);
785 snew(info->ewald_beta, info->n_entries);
786 snew(info->pme_order, info->n_entries);
787 snew(info->fn_out, info->n_entries);
788 snew(info->e_dir, info->n_entries);
789 snew(info->e_rec, info->n_entries);
793 /* Allocate and fill an array with coordinates and charges,
794 * returns the number of charges found
796 static int prepare_x_q(real* q[], rvec* x[], const gmx_mtop_t* mtop, const rvec x_orig[], t_commrec* cr)
798 int nq; /* number of charged particles */
803 snew(*q, mtop->natoms);
804 snew(*x, mtop->natoms);
807 for (const AtomProxy atomP : AtomRange(*mtop))
809 const t_atom& local = atomP.atom();
810 int i = atomP.globalAtomNumber();
811 if (is_charge(local.q))
814 (*x)[nq][XX] = x_orig[i][XX];
815 (*x)[nq][YY] = x_orig[i][YY];
816 (*x)[nq][ZZ] = x_orig[i][ZZ];
820 /* Give back some unneeded memory */
824 /* Broadcast x and q in the parallel case */
827 /* Transfer the number of charges */
831 nblock_bc(cr, nq, *x);
832 nblock_bc(cr, nq, *q);
839 /* Read in the tpr file and save information we need later in info */
840 static void read_tpr_file(const char* fn_sim_tpr,
848 read_tpx_state(fn_sim_tpr, ir, state, mtop);
850 /* The values of the original tpr input file are save in the first
851 * place [0] of the arrays */
852 info->orig_sim_steps = ir->nsteps;
853 info->pme_order[0] = ir->pme_order;
854 info->rcoulomb[0] = ir->rcoulomb;
855 info->rvdw[0] = ir->rvdw;
856 info->nkx[0] = ir->nkx;
857 info->nky[0] = ir->nky;
858 info->nkz[0] = ir->nkz;
859 info->ewald_rtol[0] = ir->ewald_rtol;
860 info->fracself = fracself;
863 info->ewald_beta[0] = user_beta;
867 info->ewald_beta[0] = calc_ewaldcoeff_q(info->rcoulomb[0], info->ewald_rtol[0]);
870 /* Check if PME was chosen */
871 if (EEL_PME(ir->coulombtype) == FALSE)
873 gmx_fatal(FARGS, "Can only do optimizations for simulations with PME");
876 /* Check if rcoulomb == rlist, which is necessary for PME */
877 if (!(ir->rcoulomb == ir->rlist))
879 gmx_fatal(FARGS, "PME requires rcoulomb (%f) to be equal to rlist (%f).", ir->rcoulomb, ir->rlist);
884 /* Transfer what we need for parallelizing the reciprocal error estimate */
885 static void bcast_info(t_inputinfo* info, const t_commrec* cr)
887 nblock_bc(cr, info->n_entries, info->nkx);
888 nblock_bc(cr, info->n_entries, info->nky);
889 nblock_bc(cr, info->n_entries, info->nkz);
890 nblock_bc(cr, info->n_entries, info->ewald_beta);
891 nblock_bc(cr, info->n_entries, info->pme_order);
892 nblock_bc(cr, info->n_entries, info->e_dir);
893 nblock_bc(cr, info->n_entries, info->e_rec);
894 block_bc(cr, info->volume);
895 block_bc(cr, info->recipbox);
896 block_bc(cr, info->natoms);
897 block_bc(cr, info->fracself);
898 block_bc(cr, info->bTUNE);
899 block_bc(cr, info->q2all);
900 block_bc(cr, info->q2allnr);
904 /* Estimate the error of the SPME Ewald sum. This estimate is based upon
905 * a) a homogeneous distribution of the charges
906 * b) a total charge of zero.
908 static void estimate_PME_error(t_inputinfo* info,
909 const t_state* state,
910 const gmx_mtop_t* mtop,
916 rvec* x = nullptr; /* The coordinates */
917 real* q = nullptr; /* The charges */
918 real edir = 0.0; /* real space error */
919 real erec = 0.0; /* reciprocal space error */
920 real derr = 0.0; /* difference of real and reciprocal space error */
921 real derr0 = 0.0; /* difference of real and reciprocal space error */
922 real beta = 0.0; /* splitting parameter beta */
923 real beta0 = 0.0; /* splitting parameter beta */
924 int ncharges; /* The number of atoms with charges */
925 int nsamples; /* The number of samples used for the calculation of the
926 * self-energy error term */
931 fprintf(fp_out, "\n--- PME ERROR ESTIMATE ---\n");
934 /* Prepare an x and q array with only the charged atoms */
935 ncharges = prepare_x_q(&q, &x, mtop, state->x.rvec_array(), cr);
938 calc_q2all(mtop, &(info->q2all), &(info->q2allnr));
939 info->ewald_rtol[0] = std::erfc(info->rcoulomb[0] * info->ewald_beta[0]);
940 /* Write some info to log file */
941 fprintf(fp_out, "Box volume : %g nm^3\n", info->volume);
942 fprintf(fp_out, "Number of charged atoms : %d (total atoms %d)\n", ncharges, info->natoms);
943 fprintf(fp_out, "Coulomb radius : %g nm\n", info->rcoulomb[0]);
944 fprintf(fp_out, "Ewald_rtol : %g\n", info->ewald_rtol[0]);
945 fprintf(fp_out, "Ewald parameter beta : %g\n", info->ewald_beta[0]);
946 fprintf(fp_out, "Interpolation order : %d\n", info->pme_order[0]);
947 fprintf(fp_out, "Fourier grid (nx,ny,nz) : %d x %d x %d\n", info->nkx[0], info->nky[0],
954 bcast_info(info, cr);
958 /* Calculate direct space error */
959 info->e_dir[0] = estimate_direct(info);
961 /* Calculate reciprocal space error */
962 info->e_rec[0] = estimate_reciprocal(info, x, q, ncharges, fp_out, bVerbose, seed, &nsamples, cr);
966 bcast_info(info, cr);
971 fprintf(fp_out, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
972 fprintf(fp_out, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
973 fprintf(fp_out, "Self-energy error term was estimated using %d samples\n", nsamples);
975 fprintf(stderr, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
976 fprintf(stderr, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
985 fprintf(stderr, "Starting tuning ...\n");
987 edir = info->e_dir[0];
988 erec = info->e_rec[0];
990 beta0 = info->ewald_beta[0];
993 info->ewald_beta[0] += 0.1;
997 info->ewald_beta[0] -= 0.1;
999 info->e_dir[0] = estimate_direct(info);
1000 info->e_rec[0] = estimate_reciprocal(info, x, q, ncharges, fp_out, bVerbose, seed, &nsamples, cr);
1004 bcast_info(info, cr);
1008 edir = info->e_dir[0];
1009 erec = info->e_rec[0];
1011 while (std::abs(derr / std::min(erec, edir)) > 1e-4)
1014 beta = info->ewald_beta[0];
1015 beta -= derr * (info->ewald_beta[0] - beta0) / (derr - derr0);
1016 beta0 = info->ewald_beta[0];
1017 info->ewald_beta[0] = beta;
1020 info->e_dir[0] = estimate_direct(info);
1022 estimate_reciprocal(info, x, q, ncharges, fp_out, bVerbose, seed, &nsamples, cr);
1026 bcast_info(info, cr);
1029 edir = info->e_dir[0];
1030 erec = info->e_rec[0];
1036 fprintf(stderr, "difference between real and rec. space error (step %d): %g\n", i,
1038 fprintf(stderr, "old beta: %f\n", beta0);
1039 fprintf(stderr, "new beta: %f\n", beta);
1043 info->ewald_rtol[0] = std::erfc(info->rcoulomb[0] * info->ewald_beta[0]);
1047 /* Write some info to log file */
1049 fprintf(fp_out, "========= After tuning ========\n");
1050 fprintf(fp_out, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
1051 fprintf(fp_out, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
1052 fprintf(stderr, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
1053 fprintf(stderr, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
1054 fprintf(fp_out, "Ewald_rtol : %g\n", info->ewald_rtol[0]);
1055 fprintf(fp_out, "Ewald parameter beta : %g\n", info->ewald_beta[0]);
1062 int gmx_pme_error(int argc, char* argv[])
1064 const char* desc[] = {
1065 "[THISMODULE] estimates the error of the electrostatic forces",
1066 "if using the sPME algorithm. The flag [TT]-tune[tt] will determine",
1067 "the splitting parameter such that the error is equally",
1068 "distributed over the real and reciprocal space part.",
1069 "The part of the error that stems from self interaction of the particles ",
1070 "is computationally demanding. However, a good a approximation is to",
1071 "just use a fraction of the particles for this term which can be",
1072 "indicated by the flag [TT]-self[tt].[PAR]",
1075 real fs = 0.0; /* 0 indicates: not set by the user */
1076 real user_beta = -1.0;
1077 real fracself = 1.0;
1079 t_state state; /* The state from the tpr input file */
1080 gmx_mtop_t mtop; /* The topology from the tpr input file */
1082 unsigned long PCA_Flags;
1083 gmx_bool bTUNE = FALSE;
1084 gmx_bool bVerbose = FALSE;
1088 static t_filenm fnm[] = { { efTPR, "-s", nullptr, ffREAD },
1089 { efOUT, "-o", "error", ffWRITE },
1090 { efTPR, "-so", "tuned", ffOPTWR } };
1092 gmx_output_env_t* oenv = nullptr;
1099 "If positive, overwrite ewald_beta from [REF].tpr[ref] file with this value" },
1104 "Tune the splitting parameter such that the error is equally distributed between "
1105 "real and reciprocal space" },
1110 "If between 0.0 and 1.0, determine self interaction error from just this "
1111 "fraction of the charged particles" },
1116 "Random number seed used for Monte Carlo algorithm when [TT]-self[tt] is set to "
1117 "a value between 0.0 and 1.0" },
1118 { "-v", FALSE, etBOOL, { &bVerbose }, "Be loud and noisy" }
1122 #define NFILE asize(fnm)
1124 CommrecHandle commrecHandle = init_commrec(MPI_COMM_WORLD, nullptr);
1125 t_commrec* cr = commrecHandle.get();
1126 PCA_Flags = PCA_NOEXIT_ON_ARGS;
1128 if (!parse_common_args(&argc, argv, PCA_Flags, NFILE, fnm, asize(pa), pa, asize(desc), desc, 0,
1136 bTUNE = opt2bSet("-so", NFILE, fnm);
1141 /* Allocate memory for the inputinfo struct: */
1143 info.fourier_sp[0] = fs;
1148 read_tpr_file(opt2fn("-s", NFILE, fnm), &info, &state, &mtop, &ir, user_beta, fracself);
1149 /* Open logfile for reading */
1150 fp = fopen(opt2fn("-o", NFILE, fnm), "w");
1152 /* Determine the volume of the simulation box */
1153 info.volume = det(state.box);
1154 calc_recipbox(state.box, info.recipbox);
1155 info.natoms = mtop.natoms;
1159 /* Check consistency if the user provided fourierspacing */
1160 if (fs > 0 && MASTER(cr))
1162 /* Recalculate the grid dimensions using fourierspacing from user input */
1166 calcFftGrid(stdout, state.box, info.fourier_sp[0], minimalPmeGridSize(info.pme_order[0]),
1167 &(info.nkx[0]), &(info.nky[0]), &(info.nkz[0]));
1168 if ((ir.nkx != info.nkx[0]) || (ir.nky != info.nky[0]) || (ir.nkz != info.nkz[0]))
1171 "Wrong fourierspacing %f nm, input file grid = %d x %d x %d, computed grid = "
1173 fs, ir.nkx, ir.nky, ir.nkz, info.nkx[0], info.nky[0], info.nkz[0]);
1177 /* Estimate (S)PME force error */
1181 bcast_info(&info, cr);
1184 /* Get an error estimate of the input tpr file and do some tuning if requested */
1185 estimate_PME_error(&info, &state, &mtop, fp, bVerbose, seed, cr);
1189 /* Write out optimized tpr file if requested */
1190 if (opt2bSet("-so", NFILE, fnm) || bTUNE)
1192 ir.ewald_rtol = info.ewald_rtol[0];
1193 write_tpx_state(opt2fn("-so", NFILE, fnm), &ir, &state, &mtop);
1195 please_cite(fp, "Wang2010");