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,2021, 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 "
171 /* Estimate the direct space part error of the SPME Ewald sum */
172 static real estimate_direct(t_inputinfo* info)
174 real e_dir = 0; /* Error estimate */
175 real beta = 0; /* Splitting parameter (1/nm) */
176 real r_coulomb = 0; /* Cut-off in direct space */
179 beta = info->ewald_beta[0];
180 r_coulomb = info->rcoulomb[0];
182 e_dir = 2.0 * info->q2all * gmx::invsqrt(info->q2allnr * r_coulomb * info->volume);
183 e_dir *= std::exp(-beta * beta * r_coulomb * r_coulomb);
185 return gmx::c_one4PiEps0 * e_dir;
190 /* the following 4 functions determine polynomials required for the reciprocal error estimate */
192 static inline real eps_poly1(real m, /* grid coordinate in certain direction */
193 real K, /* grid size in corresponding direction */
194 real n) /* spline interpolation order of the SPME */
197 real nom = 0; /* nominator */
198 real denom = 0; /* denominator */
206 for (i = -SUMORDER; i < 0; i++)
210 nom += std::pow(tmp, -n);
213 for (i = SUMORDER; i > 0; i--)
217 nom += std::pow(tmp, -n);
222 denom = std::pow(tmp, -n) + nom;
227 static inline real eps_poly2(real m, /* grid coordinate in certain direction */
228 real K, /* grid size in corresponding direction */
229 real n) /* spline interpolation order of the SPME */
232 real nom = 0; /* nominator */
233 real denom = 0; /* denominator */
241 for (i = -SUMORDER; i < 0; i++)
245 nom += std::pow(tmp, -2 * n);
248 for (i = SUMORDER; i > 0; i--)
252 nom += std::pow(tmp, -2 * n);
255 for (i = -SUMORDER; i < SUMORDER + 1; i++)
259 denom += std::pow(tmp, -n);
261 tmp = eps_poly1(m, K, n);
262 return nom / denom / denom + tmp * tmp;
265 static inline real eps_poly3(real m, /* grid coordinate in certain direction */
266 real K, /* grid size in corresponding direction */
267 real n) /* spline interpolation order of the SPME */
270 real nom = 0; /* nominator */
271 real denom = 0; /* denominator */
279 for (i = -SUMORDER; i < 0; i++)
283 nom += i * std::pow(tmp, -2 * n);
286 for (i = SUMORDER; i > 0; i--)
290 nom += i * std::pow(tmp, -2 * n);
293 for (i = -SUMORDER; i < SUMORDER + 1; i++)
297 denom += std::pow(tmp, -n);
300 return 2.0 * M_PI * nom / denom / denom;
303 static inline real eps_poly4(real m, /* grid coordinate in certain direction */
304 real K, /* grid size in corresponding direction */
305 real n) /* spline interpolation order of the SPME */
308 real nom = 0; /* nominator */
309 real denom = 0; /* denominator */
317 for (i = -SUMORDER; i < 0; i++)
321 nom += i * i * std::pow(tmp, -2 * n);
324 for (i = SUMORDER; i > 0; i--)
328 nom += i * i * std::pow(tmp, -2 * n);
331 for (i = -SUMORDER; i < SUMORDER + 1; i++)
335 denom += std::pow(tmp, -n);
338 return 4.0 * M_PI * M_PI * nom / denom / denom;
341 static inline real eps_self(real m, /* grid coordinate in certain direction */
342 real K, /* grid size in corresponding direction */
343 rvec rboxv, /* reciprocal box vector */
344 real n, /* spline interpolation order of the SPME */
345 rvec x) /* coordinate of charge */
348 real tmp = 0; /* temporary variables for computations */
349 real tmp1 = 0; /* temporary variables for computations */
350 real tmp2 = 0; /* temporary variables for computations */
351 real rcoord = 0; /* coordinate in certain reciprocal space direction */
352 real nom = 0; /* nominator */
353 real denom = 0; /* denominator */
361 rcoord = iprod(rboxv, x);
364 for (i = -SUMORDER; i < 0; i++)
366 tmp = -std::sin(2.0 * M_PI * i * K * rcoord);
367 tmp1 = 2.0 * M_PI * m / K + 2.0 * M_PI * i;
368 tmp2 = std::pow(tmp1, -n);
369 nom += tmp * tmp2 * i;
373 for (i = SUMORDER; i > 0; i--)
375 tmp = -std::sin(2.0 * M_PI * i * K * rcoord);
376 tmp1 = 2.0 * M_PI * m / K + 2.0 * M_PI * i;
377 tmp2 = std::pow(tmp1, -n);
378 nom += tmp * tmp2 * i;
383 tmp = 2.0 * M_PI * m / K;
387 return 2.0 * M_PI * nom / denom * K;
392 /* The following routine is just a copy from pme.c */
394 static void calc_recipbox(matrix box, matrix recipbox)
396 /* Save some time by assuming upper right part is zero */
398 real tmp = 1.0 / (box[XX][XX] * box[YY][YY] * box[ZZ][ZZ]);
400 recipbox[XX][XX] = box[YY][YY] * box[ZZ][ZZ] * tmp;
401 recipbox[XX][YY] = 0;
402 recipbox[XX][ZZ] = 0;
403 recipbox[YY][XX] = -box[YY][XX] * box[ZZ][ZZ] * tmp;
404 recipbox[YY][YY] = box[XX][XX] * box[ZZ][ZZ] * tmp;
405 recipbox[YY][ZZ] = 0;
406 recipbox[ZZ][XX] = (box[YY][XX] * box[ZZ][YY] - box[YY][YY] * box[ZZ][XX]) * tmp;
407 recipbox[ZZ][YY] = -box[ZZ][YY] * box[XX][XX] * tmp;
408 recipbox[ZZ][ZZ] = box[XX][XX] * box[YY][YY] * tmp;
412 /* Estimate the reciprocal space part error of the SPME Ewald sum. */
413 static real estimate_reciprocal(t_inputinfo* info,
414 rvec x[], /* array of particles */
415 const real q[], /* array of charges */
416 int nr, /* number of charges = size of the charge array */
417 FILE gmx_unused* fp_out,
419 int seed, /* The seed for the random number generator */
420 int* nsamples, /* Return the number of samples used if Monte Carlo
421 * algorithm is used for self energy error estimate */
424 real e_rec = 0; /* reciprocal error estimate */
425 real e_rec1 = 0; /* Error estimate term 1*/
426 real e_rec2 = 0; /* Error estimate term 2*/
427 real e_rec3 = 0; /* Error estimate term 3 */
428 real e_rec3x = 0; /* part of Error estimate term 3 in x */
429 real e_rec3y = 0; /* part of Error estimate term 3 in y */
430 real e_rec3z = 0; /* part of Error estimate term 3 in z */
432 int nx, ny, nz; /* grid coordinates */
433 real q2_all = 0; /* sum of squared charges */
434 rvec gridpx; /* reciprocal grid point in x direction*/
435 rvec gridpxy; /* reciprocal grid point in x and y direction*/
436 rvec gridp; /* complete reciprocal grid point in 3 directions*/
437 rvec tmpvec; /* template to create points from basis vectors */
438 rvec tmpvec2; /* template to create points from basis vectors */
439 real coeff = 0; /* variable to compute coefficients of the error estimate */
440 real coeff2 = 0; /* variable to compute coefficients of the error estimate */
441 real tmp = 0; /* variables to compute different factors from vectors */
446 int* numbers = nullptr;
448 /* Index variables for parallel work distribution */
449 int startglobal, stopglobal;
450 int startlocal, stoplocal;
459 GMX_RELEASE_ASSERT(q != nullptr, "Must have charges");
463 seed = static_cast<int>(gmx::makeRandomSeed());
465 fprintf(stderr, "Using random seed %d.\n", seed);
467 gmx::DefaultRandomEngine rng(seed);
468 gmx::UniformIntDistribution<int> dist(0, nr - 1);
476 for (i = 0; i < nr; i++)
478 q2_all += q[i] * q[i];
481 /* Calculate indices for work distribution */
482 startglobal = -info->nkx[0] / 2;
483 stopglobal = info->nkx[0] / 2;
484 xtot = stopglobal * 2 + 1;
487 x_per_core = static_cast<int>(std::ceil(static_cast<real>(xtot) / cr->nnodes));
488 startlocal = startglobal + x_per_core * cr->nodeid;
489 stoplocal = startlocal + x_per_core - 1;
490 if (stoplocal > stopglobal)
492 stoplocal = stopglobal;
497 startlocal = startglobal;
498 stoplocal = stopglobal;
514 fprintf(stderr, "Calculating reciprocal error part 1 ...");
517 for (nx = startlocal; nx <= stoplocal; nx++)
519 svmul(nx, info->recipbox[XX], gridpx);
520 for (ny = -info->nky[0] / 2; ny < info->nky[0] / 2 + 1; ny++)
522 svmul(ny, info->recipbox[YY], tmpvec);
523 rvec_add(gridpx, tmpvec, gridpxy);
524 for (nz = -info->nkz[0] / 2; nz < info->nkz[0] / 2 + 1; nz++)
526 if (0 == nx && 0 == ny && 0 == nz)
530 svmul(nz, info->recipbox[ZZ], tmpvec);
531 rvec_add(gridpxy, tmpvec, gridp);
533 coeff = std::exp(-1.0 * M_PI * M_PI * tmp / info->ewald_beta[0] / info->ewald_beta[0]);
534 coeff /= 2.0 * M_PI * info->volume * tmp;
538 tmp = eps_poly2(nx, info->nkx[0], info->pme_order[0]);
539 tmp += eps_poly2(ny, info->nkx[0], info->pme_order[0]);
540 tmp += eps_poly2(nz, info->nkx[0], info->pme_order[0]);
542 tmp1 = eps_poly1(nx, info->nkx[0], info->pme_order[0]);
543 tmp2 = eps_poly1(ny, info->nky[0], info->pme_order[0]);
545 tmp += 2.0 * tmp1 * tmp2;
547 tmp1 = eps_poly1(nz, info->nkz[0], info->pme_order[0]);
548 tmp2 = eps_poly1(ny, info->nky[0], info->pme_order[0]);
550 tmp += 2.0 * tmp1 * tmp2;
552 tmp1 = eps_poly1(nz, info->nkz[0], info->pme_order[0]);
553 tmp2 = eps_poly1(nx, info->nkx[0], info->pme_order[0]);
555 tmp += 2.0 * tmp1 * tmp2;
557 tmp1 = eps_poly1(nx, info->nkx[0], info->pme_order[0]);
558 tmp1 += eps_poly1(ny, info->nky[0], info->pme_order[0]);
559 tmp1 += eps_poly1(nz, info->nkz[0], info->pme_order[0]);
563 e_rec1 += 32.0 * M_PI * M_PI * coeff * coeff * coeff2 * tmp * q2_all * q2_all / nr;
565 tmp1 = eps_poly3(nx, info->nkx[0], info->pme_order[0]);
566 tmp1 *= info->nkx[0];
567 tmp2 = iprod(gridp, info->recipbox[XX]);
571 tmp1 = eps_poly3(ny, info->nky[0], info->pme_order[0]);
572 tmp1 *= info->nky[0];
573 tmp2 = iprod(gridp, info->recipbox[YY]);
577 tmp1 = eps_poly3(nz, info->nkz[0], info->pme_order[0]);
578 tmp1 *= info->nkz[0];
579 tmp2 = iprod(gridp, info->recipbox[ZZ]);
585 tmp1 = eps_poly4(nx, info->nkx[0], info->pme_order[0]);
586 tmp1 *= norm2(info->recipbox[XX]);
587 tmp1 *= info->nkx[0] * info->nkx[0];
591 tmp1 = eps_poly4(ny, info->nky[0], info->pme_order[0]);
592 tmp1 *= norm2(info->recipbox[YY]);
593 tmp1 *= info->nky[0] * info->nky[0];
597 tmp1 = eps_poly4(nz, info->nkz[0], info->pme_order[0]);
598 tmp1 *= norm2(info->recipbox[ZZ]);
599 tmp1 *= info->nkz[0] * info->nkz[0];
603 e_rec2 += 4.0 * coeff * coeff * tmp * q2_all * q2_all / nr;
609 "\rCalculating reciprocal error part 1 ... %3.0f%%",
610 100.0 * (nx - startlocal + 1) / (x_per_core));
617 fprintf(stderr, "\n");
620 /* Use just a fraction of all charges to estimate the self energy error term? */
621 bFraction = (info->fracself > 0.0) && (info->fracself < 1.0);
625 /* Here xtot is the number of samples taken for the Monte Carlo calculation
626 * of the average of term IV of equation 35 in Wang2010. Round up to a
627 * number of samples that is divisible by the number of nodes */
628 x_per_core = static_cast<int>(std::ceil(info->fracself * nr / cr->nnodes));
629 xtot = x_per_core * cr->nnodes;
633 /* In this case we use all nr particle positions */
635 x_per_core = static_cast<int>(std::ceil(static_cast<real>(xtot) / cr->nnodes));
638 startlocal = x_per_core * cr->nodeid;
639 stoplocal = std::min(startlocal + x_per_core, xtot); /* min needed if xtot == nr */
643 /* Make shure we get identical results in serial and parallel. Therefore,
644 * take the sample indices from a single, global random number array that
645 * is constructed on the master node and that only depends on the seed */
649 for (i = 0; i < xtot; i++)
651 numbers[i] = dist(rng); // [0,nr-1]
654 /* Broadcast the random number array to the other nodes */
657 nblock_bc(cr->mpi_comm_mygroup, xtot, numbers);
660 if (bVerbose && MASTER(cr))
663 "Using %d sample%s to approximate the self interaction error term",
665 xtot == 1 ? "" : "s");
668 fprintf(stdout, " (%d sample%s per rank)", x_per_core, x_per_core == 1 ? "" : "s");
670 fprintf(stdout, ".\n");
674 /* Return the number of positions used for the Monte Carlo algorithm */
677 for (i = startlocal; i < stoplocal; i++)
685 /* Randomly pick a charge */
690 /* Use all charges */
694 /* for(nx=startlocal; nx<=stoplocal; nx++)*/
695 for (nx = -info->nkx[0] / 2; nx < info->nkx[0] / 2 + 1; nx++)
697 svmul(nx, info->recipbox[XX], gridpx);
698 for (ny = -info->nky[0] / 2; ny < info->nky[0] / 2 + 1; ny++)
700 svmul(ny, info->recipbox[YY], tmpvec);
701 rvec_add(gridpx, tmpvec, gridpxy);
702 for (nz = -info->nkz[0] / 2; nz < info->nkz[0] / 2 + 1; nz++)
705 if (0 == nx && 0 == ny && 0 == nz)
710 svmul(nz, info->recipbox[ZZ], tmpvec);
711 rvec_add(gridpxy, tmpvec, gridp);
713 coeff = std::exp(-1.0 * M_PI * M_PI * tmp / info->ewald_beta[0] / info->ewald_beta[0]);
716 * eps_self(nx, info->nkx[0], info->recipbox[XX], info->pme_order[0], x[ci]);
718 * eps_self(ny, info->nky[0], info->recipbox[YY], info->pme_order[0], x[ci]);
720 * eps_self(nz, info->nkz[0], info->recipbox[ZZ], info->pme_order[0], x[ci]);
727 svmul(e_rec3x, info->recipbox[XX], tmpvec);
728 rvec_inc(tmpvec2, tmpvec);
729 svmul(e_rec3y, info->recipbox[YY], tmpvec);
730 rvec_inc(tmpvec2, tmpvec);
731 svmul(e_rec3z, info->recipbox[ZZ], tmpvec);
732 rvec_inc(tmpvec2, tmpvec);
734 e_rec3 += q[ci] * q[ci] * q[ci] * q[ci] * norm2(tmpvec2)
735 / (xtot * M_PI * info->volume * M_PI * info->volume);
738 fprintf(stderr, "\rCalculating reciprocal error part 2 ... %3.0f%%", 100.0 * (i + 1) / stoplocal);
745 fprintf(stderr, "\n");
752 t1 = MPI_Wtime() - t0;
753 fprintf(fp_out, "Recip. err. est. took : %lf s\n", t1);
761 fprintf(stderr, "Rank %3d: nx=[%3d...%3d] e_rec3=%e\n", cr->nodeid, startlocal, stoplocal, e_rec3);
767 gmx_sum(1, &e_rec1, cr);
768 gmx_sum(1, &e_rec2, cr);
769 gmx_sum(1, &e_rec3, cr);
772 /* e_rec1*=8.0 * q2_all / info->volume / info->volume / nr ;
773 e_rec2*= q2_all / M_PI / M_PI / info->volume / info->volume / nr ;
774 e_rec3/= M_PI * M_PI * info->volume * info->volume * nr ;
776 e_rec = std::sqrt(e_rec1 + e_rec2 + e_rec3);
779 return gmx::c_one4PiEps0 * e_rec;
783 /* Allocate memory for the inputinfo struct: */
784 static void create_info(t_inputinfo* info)
786 snew(info->fac, info->n_entries);
787 snew(info->rcoulomb, info->n_entries);
788 snew(info->rvdw, info->n_entries);
789 snew(info->nkx, info->n_entries);
790 snew(info->nky, info->n_entries);
791 snew(info->nkz, info->n_entries);
792 snew(info->fourier_sp, info->n_entries);
793 snew(info->ewald_rtol, info->n_entries);
794 snew(info->ewald_beta, info->n_entries);
795 snew(info->pme_order, info->n_entries);
796 snew(info->fn_out, info->n_entries);
797 snew(info->e_dir, info->n_entries);
798 snew(info->e_rec, info->n_entries);
802 /* Allocate and fill an array with coordinates and charges,
803 * returns the number of charges found
805 static int prepare_x_q(real* q[], rvec* x[], const gmx_mtop_t* mtop, const rvec x_orig[], t_commrec* cr)
807 int nq; /* number of charged particles */
812 snew(*q, mtop->natoms);
813 snew(*x, mtop->natoms);
816 for (const AtomProxy atomP : AtomRange(*mtop))
818 const t_atom& local = atomP.atom();
819 int i = atomP.globalAtomNumber();
820 if (is_charge(local.q))
823 (*x)[nq][XX] = x_orig[i][XX];
824 (*x)[nq][YY] = x_orig[i][YY];
825 (*x)[nq][ZZ] = x_orig[i][ZZ];
829 /* Give back some unneeded memory */
833 /* Broadcast x and q in the parallel case */
836 /* Transfer the number of charges */
837 block_bc(cr->mpi_comm_mygroup, nq);
838 snew_bc(MASTER(cr), *x, nq);
839 snew_bc(MASTER(cr), *q, nq);
840 nblock_bc(cr->mpi_comm_mygroup, nq, *x);
841 nblock_bc(cr->mpi_comm_mygroup, nq, *q);
848 /* Read in the tpr file and save information we need later in info */
849 static void read_tpr_file(const char* fn_sim_tpr,
857 read_tpx_state(fn_sim_tpr, ir, state, mtop);
859 /* The values of the original tpr input file are save in the first
860 * place [0] of the arrays */
861 info->orig_sim_steps = ir->nsteps;
862 info->pme_order[0] = ir->pme_order;
863 info->rcoulomb[0] = ir->rcoulomb;
864 info->rvdw[0] = ir->rvdw;
865 info->nkx[0] = ir->nkx;
866 info->nky[0] = ir->nky;
867 info->nkz[0] = ir->nkz;
868 info->ewald_rtol[0] = ir->ewald_rtol;
869 info->fracself = fracself;
872 info->ewald_beta[0] = user_beta;
876 info->ewald_beta[0] = calc_ewaldcoeff_q(info->rcoulomb[0], info->ewald_rtol[0]);
879 /* Check if PME was chosen */
880 if (EEL_PME(ir->coulombtype) == FALSE)
882 gmx_fatal(FARGS, "Can only do optimizations for simulations with PME");
885 /* Check if rcoulomb == rlist, which is necessary for PME */
886 if (!(ir->rcoulomb == ir->rlist))
888 gmx_fatal(FARGS, "PME requires rcoulomb (%f) to be equal to rlist (%f).", ir->rcoulomb, ir->rlist);
893 /* Transfer what we need for parallelizing the reciprocal error estimate */
894 static void bcast_info(t_inputinfo* info, const t_commrec* cr)
896 nblock_bc(cr->mpi_comm_mygroup, info->n_entries, info->nkx);
897 nblock_bc(cr->mpi_comm_mygroup, info->n_entries, info->nky);
898 nblock_bc(cr->mpi_comm_mygroup, info->n_entries, info->nkz);
899 nblock_bc(cr->mpi_comm_mygroup, info->n_entries, info->ewald_beta);
900 nblock_bc(cr->mpi_comm_mygroup, info->n_entries, info->pme_order);
901 nblock_bc(cr->mpi_comm_mygroup, info->n_entries, info->e_dir);
902 nblock_bc(cr->mpi_comm_mygroup, info->n_entries, info->e_rec);
903 block_bc(cr->mpi_comm_mygroup, info->volume);
904 block_bc(cr->mpi_comm_mygroup, info->recipbox);
905 block_bc(cr->mpi_comm_mygroup, info->natoms);
906 block_bc(cr->mpi_comm_mygroup, info->fracself);
907 block_bc(cr->mpi_comm_mygroup, info->bTUNE);
908 block_bc(cr->mpi_comm_mygroup, info->q2all);
909 block_bc(cr->mpi_comm_mygroup, info->q2allnr);
913 /* Estimate the error of the SPME Ewald sum. This estimate is based upon
914 * a) a homogeneous distribution of the charges
915 * b) a total charge of zero.
917 static void estimate_PME_error(t_inputinfo* info,
918 const t_state* state,
919 const gmx_mtop_t* mtop,
925 rvec* x = nullptr; /* The coordinates */
926 real* q = nullptr; /* The charges */
927 real edir = 0.0; /* real space error */
928 real erec = 0.0; /* reciprocal space error */
929 real derr = 0.0; /* difference of real and reciprocal space error */
930 real derr0 = 0.0; /* difference of real and reciprocal space error */
931 real beta = 0.0; /* splitting parameter beta */
932 real beta0 = 0.0; /* splitting parameter beta */
933 int ncharges; /* The number of atoms with charges */
934 int nsamples; /* The number of samples used for the calculation of the
935 * self-energy error term */
940 fprintf(fp_out, "\n--- PME ERROR ESTIMATE ---\n");
943 /* Prepare an x and q array with only the charged atoms */
944 ncharges = prepare_x_q(&q, &x, mtop, state->x.rvec_array(), cr);
947 calc_q2all(mtop, &(info->q2all), &(info->q2allnr));
948 info->ewald_rtol[0] = std::erfc(info->rcoulomb[0] * info->ewald_beta[0]);
949 /* Write some info to log file */
950 fprintf(fp_out, "Box volume : %g nm^3\n", info->volume);
951 fprintf(fp_out, "Number of charged atoms : %d (total atoms %d)\n", ncharges, info->natoms);
952 fprintf(fp_out, "Coulomb radius : %g nm\n", info->rcoulomb[0]);
953 fprintf(fp_out, "Ewald_rtol : %g\n", info->ewald_rtol[0]);
954 fprintf(fp_out, "Ewald parameter beta : %g\n", info->ewald_beta[0]);
955 fprintf(fp_out, "Interpolation order : %d\n", info->pme_order[0]);
956 fprintf(fp_out, "Fourier grid (nx,ny,nz) : %d x %d x %d\n", info->nkx[0], info->nky[0], info->nkz[0]);
962 bcast_info(info, cr);
966 /* Calculate direct space error */
967 info->e_dir[0] = estimate_direct(info);
969 /* Calculate reciprocal space error */
970 info->e_rec[0] = estimate_reciprocal(info, x, q, ncharges, fp_out, bVerbose, seed, &nsamples, cr);
974 bcast_info(info, cr);
979 fprintf(fp_out, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
980 fprintf(fp_out, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
981 fprintf(fp_out, "Self-energy error term was estimated using %d samples\n", nsamples);
983 fprintf(stderr, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
984 fprintf(stderr, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
993 fprintf(stderr, "Starting tuning ...\n");
995 edir = info->e_dir[0];
996 erec = info->e_rec[0];
998 beta0 = info->ewald_beta[0];
1001 info->ewald_beta[0] += 0.1;
1005 info->ewald_beta[0] -= 0.1;
1007 info->e_dir[0] = estimate_direct(info);
1008 info->e_rec[0] = estimate_reciprocal(info, x, q, ncharges, fp_out, bVerbose, seed, &nsamples, cr);
1012 bcast_info(info, cr);
1016 edir = info->e_dir[0];
1017 erec = info->e_rec[0];
1019 while (std::abs(derr / std::min(erec, edir)) > 1e-4)
1022 beta = info->ewald_beta[0];
1023 beta -= derr * (info->ewald_beta[0] - beta0) / (derr - derr0);
1024 beta0 = info->ewald_beta[0];
1025 info->ewald_beta[0] = beta;
1028 info->e_dir[0] = estimate_direct(info);
1030 estimate_reciprocal(info, x, q, ncharges, fp_out, bVerbose, seed, &nsamples, cr);
1034 bcast_info(info, cr);
1037 edir = info->e_dir[0];
1038 erec = info->e_rec[0];
1045 "difference between real and rec. space error (step %d): %g\n",
1048 fprintf(stderr, "old beta: %f\n", beta0);
1049 fprintf(stderr, "new beta: %f\n", beta);
1053 info->ewald_rtol[0] = std::erfc(info->rcoulomb[0] * info->ewald_beta[0]);
1057 /* Write some info to log file */
1059 fprintf(fp_out, "========= After tuning ========\n");
1060 fprintf(fp_out, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
1061 fprintf(fp_out, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
1062 fprintf(stderr, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
1063 fprintf(stderr, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
1064 fprintf(fp_out, "Ewald_rtol : %g\n", info->ewald_rtol[0]);
1065 fprintf(fp_out, "Ewald parameter beta : %g\n", info->ewald_beta[0]);
1072 int gmx_pme_error(int argc, char* argv[])
1074 const char* desc[] = {
1075 "[THISMODULE] estimates the error of the electrostatic forces",
1076 "if using the sPME algorithm. The flag [TT]-tune[tt] will determine",
1077 "the splitting parameter such that the error is equally",
1078 "distributed over the real and reciprocal space part.",
1079 "The part of the error that stems from self interaction of the particles ",
1080 "is computationally demanding. However, a good a approximation is to",
1081 "just use a fraction of the particles for this term which can be",
1082 "indicated by the flag [TT]-self[tt].[PAR]",
1085 real fs = 0.0; /* 0 indicates: not set by the user */
1086 real user_beta = -1.0;
1087 real fracself = 1.0;
1089 t_state state; /* The state from the tpr input file */
1090 gmx_mtop_t mtop; /* The topology from the tpr input file */
1092 unsigned long PCA_Flags;
1093 gmx_bool bTUNE = FALSE;
1094 gmx_bool bVerbose = FALSE;
1098 static t_filenm fnm[] = { { efTPR, "-s", nullptr, ffREAD },
1099 { efOUT, "-o", "error", ffWRITE },
1100 { efTPR, "-so", "tuned", ffOPTWR } };
1102 gmx_output_env_t* oenv = nullptr;
1109 "If positive, overwrite ewald_beta from [REF].tpr[ref] file with this value" },
1114 "Tune the splitting parameter such that the error is equally distributed between "
1115 "real and reciprocal space" },
1120 "If between 0.0 and 1.0, determine self interaction error from just this "
1121 "fraction of the charged particles" },
1126 "Random number seed used for Monte Carlo algorithm when [TT]-self[tt] is set to "
1127 "a value between 0.0 and 1.0" },
1128 { "-v", FALSE, etBOOL, { &bVerbose }, "Be loud and noisy" }
1132 #define NFILE asize(fnm)
1134 CommrecHandle commrecHandle = init_commrec(MPI_COMM_WORLD);
1135 t_commrec* cr = commrecHandle.get();
1136 PCA_Flags = PCA_NOEXIT_ON_ARGS;
1138 if (!parse_common_args(
1139 &argc, argv, PCA_Flags, NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
1146 bTUNE = opt2bSet("-so", NFILE, fnm);
1151 /* Allocate memory for the inputinfo struct: */
1153 info.fourier_sp[0] = fs;
1158 read_tpr_file(opt2fn("-s", NFILE, fnm), &info, &state, &mtop, &ir, user_beta, fracself);
1159 /* Open logfile for reading */
1160 fp = fopen(opt2fn("-o", NFILE, fnm), "w");
1162 /* Determine the volume of the simulation box */
1163 info.volume = det(state.box);
1164 calc_recipbox(state.box, info.recipbox);
1165 info.natoms = mtop.natoms;
1169 /* Check consistency if the user provided fourierspacing */
1170 if (fs > 0 && MASTER(cr))
1172 /* Recalculate the grid dimensions using fourierspacing from user input */
1179 minimalPmeGridSize(info.pme_order[0]),
1183 if ((ir.nkx != info.nkx[0]) || (ir.nky != info.nky[0]) || (ir.nkz != info.nkz[0]))
1186 "Wrong fourierspacing %f nm, input file grid = %d x %d x %d, computed grid = "
1198 /* Estimate (S)PME force error */
1202 bcast_info(&info, cr);
1205 /* Get an error estimate of the input tpr file and do some tuning if requested */
1206 estimate_PME_error(&info, &state, &mtop, fp, bVerbose, seed, cr);
1210 /* Write out optimized tpr file if requested */
1211 if (opt2bSet("-so", NFILE, fnm) || bTUNE)
1213 ir.ewald_rtol = info.ewald_rtol[0];
1214 write_tpx_state(opt2fn("-so", NFILE, fnm), &ir, &state, &mtop);
1216 please_cite(fp, "Wang2010");