2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2010,2011,2012,2013,2014, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
37 #include "gromacs/commandline/pargs.h"
39 #include "types/commrec.h"
40 #include "gromacs/utility/smalloc.h"
43 #include "gromacs/fileio/tpxio.h"
46 #include "checkpoint.h"
48 #include "gromacs/random/random.h"
52 #include "mtop_util.h"
57 #include "gromacs/legacyheaders/gmx_fatal.h"
59 /* We use the same defines as in mvdata.c here */
60 #define block_bc(cr, d) gmx_bcast( sizeof(d), &(d), (cr))
61 #define nblock_bc(cr, nr, d) gmx_bcast((nr)*sizeof((d)[0]), (d), (cr))
62 #define snew_bc(cr, d, nr) { if (!MASTER(cr)) {snew((d), (nr)); }}
63 /* #define TAKETIME */
66 ddnoSEL, ddnoINTERLEAVE, ddnoPP_PME, ddnoCARTESIAN, ddnoNR
69 /* Enum for situations that can occur during log file parsing */
75 eParselogResetProblem,
82 gmx_int64_t orig_sim_steps; /* Number of steps to be done in the real simulation */
83 int n_entries; /* Number of entries in arrays */
84 real volume; /* The volume of the box */
85 matrix recipbox; /* The reciprocal box */
86 int natoms; /* The number of atoms in the MD system */
87 real *fac; /* The scaling factor */
88 real *rcoulomb; /* The coulomb radii [0...nr_inputfiles] */
89 real *rvdw; /* The vdW radii */
90 int *nkx, *nky, *nkz; /* Number of k vectors in each spatial dimension */
91 real *fourier_sp; /* Fourierspacing */
92 real *ewald_rtol; /* Real space tolerance for Ewald, determines */
93 /* the real/reciprocal space relative weight */
94 real *ewald_beta; /* Splitting parameter [1/nm] */
95 real fracself; /* fraction of particles for SI error */
96 real q2all; /* sum ( q ^2 ) */
97 real q2allnr; /* nr of charges */
98 int *pme_order; /* Interpolation order for PME (bsplines) */
99 char **fn_out; /* Name of the output tpr file */
100 real *e_dir; /* Direct space part of PME error with these settings */
101 real *e_rec; /* Reciprocal space part of PME error */
102 gmx_bool bTUNE; /* flag for tuning */
106 /* Returns TRUE when atom is charged */
107 static gmx_bool is_charge(real charge)
109 if (charge*charge > GMX_REAL_EPS)
120 /* calculate charge density */
121 static void calc_q2all(
122 gmx_mtop_t *mtop, /* molecular topology */
123 real *q2all, real *q2allnr)
125 int imol, iatom; /* indices for loops */
126 real q2_all = 0; /* Sum of squared charges */
127 int nrq_mol; /* Number of charges in a single molecule */
128 int nrq_all; /* Total number of charges in the MD system */
130 gmx_moltype_t *molecule;
131 gmx_molblock_t *molblock;
134 fprintf(stderr, "\nCharge density:\n");
136 q2_all = 0.0; /* total q squared */
137 nrq_all = 0; /* total number of charges in the system */
138 for (imol = 0; imol < mtop->nmolblock; imol++) /* Loop over molecule types */
140 q2_mol = 0.0; /* q squared value of this molecule */
141 nrq_mol = 0; /* number of charges this molecule carries */
142 molecule = &(mtop->moltype[imol]);
143 molblock = &(mtop->molblock[imol]);
144 for (iatom = 0; iatom < molblock->natoms_mol; iatom++) /* Loop over atoms in this molecule */
146 qi = molecule->atoms.atom[iatom].q; /* Charge of this atom */
147 /* Is this charge worth to be considered? */
154 /* Multiply with the number of molecules present of this type and add */
155 q2_all += q2_mol*molblock->nmol;
156 nrq_all += nrq_mol*molblock->nmol;
158 fprintf(stderr, "Molecule %2d (%5d atoms) q2_mol=%10.3e nr.mol.charges=%5d (%6dx) q2_all=%10.3e tot.charges=%d\n",
159 imol, molblock->natoms_mol, q2_mol, nrq_mol, molblock->nmol, q2_all, nrq_all);
169 /* Estimate the direct space part error of the SPME Ewald sum */
170 static real estimate_direct(
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 *= exp (-beta*beta*r_coulomb*r_coulomb);
185 return ONE_4PI_EPS0*e_dir;
190 /* the following 4 functions determine polynomials required for the reciprocal error estimate */
192 static inline real eps_poly1(
193 real m, /* grid coordinate in certain direction */
194 real K, /* grid size in corresponding direction */
195 real n) /* spline interpolation order of the SPME */
198 real nom = 0; /* nominator */
199 real denom = 0; /* denominator */
207 for (i = -SUMORDER; i < 0; i++)
211 nom += pow( tmp, -n );
214 for (i = SUMORDER; i > 0; i--)
218 nom += pow( tmp, -n );
223 denom = pow( tmp, -n )+nom;
229 static inline real eps_poly2(
230 real m, /* grid coordinate in certain direction */
231 real K, /* grid size in corresponding direction */
232 real n) /* spline interpolation order of the SPME */
235 real nom = 0; /* nominator */
236 real denom = 0; /* denominator */
244 for (i = -SUMORDER; i < 0; i++)
248 nom += pow( tmp, -2*n );
251 for (i = SUMORDER; i > 0; i--)
255 nom += pow( tmp, -2*n );
258 for (i = -SUMORDER; i < SUMORDER+1; i++)
262 denom += pow( tmp, -n );
264 tmp = eps_poly1(m, K, n);
265 return nom / denom / denom + tmp*tmp;
269 static inline real eps_poly3(
270 real m, /* grid coordinate in certain direction */
271 real K, /* grid size in corresponding direction */
272 real n) /* spline interpolation order of the SPME */
275 real nom = 0; /* nominator */
276 real denom = 0; /* denominator */
284 for (i = -SUMORDER; i < 0; i++)
288 nom += i * pow( tmp, -2*n );
291 for (i = SUMORDER; i > 0; i--)
295 nom += i * pow( tmp, -2*n );
298 for (i = -SUMORDER; i < SUMORDER+1; i++)
302 denom += pow( tmp, -n );
305 return 2.0 * M_PI * nom / denom / denom;
309 static inline real eps_poly4(
310 real m, /* grid coordinate in certain direction */
311 real K, /* grid size in corresponding direction */
312 real n) /* spline interpolation order of the SPME */
315 real nom = 0; /* nominator */
316 real denom = 0; /* denominator */
324 for (i = -SUMORDER; i < 0; i++)
328 nom += i * i * pow( tmp, -2*n );
331 for (i = SUMORDER; i > 0; i--)
335 nom += i * i * pow( tmp, -2*n );
338 for (i = -SUMORDER; i < SUMORDER+1; i++)
342 denom += pow( tmp, -n );
345 return 4.0 * M_PI * M_PI * nom / denom / denom;
349 static inline real eps_self(
350 real m, /* grid coordinate in certain direction */
351 real K, /* grid size in corresponding direction */
352 rvec rboxv, /* reciprocal box vector */
353 real n, /* spline interpolation order of the SPME */
354 rvec x) /* coordinate of charge */
357 real tmp = 0; /* temporary variables for computations */
358 real tmp1 = 0; /* temporary variables for computations */
359 real tmp2 = 0; /* temporary variables for computations */
360 real rcoord = 0; /* coordinate in certain reciprocal space direction */
361 real nom = 0; /* nominator */
362 real denom = 0; /* denominator */
370 rcoord = iprod(rboxv, x);
373 for (i = -SUMORDER; i < 0; i++)
375 tmp = -sin(2.0 * M_PI * i * K * rcoord);
376 tmp1 = 2.0 * M_PI * m / K + 2.0 * M_PI * i;
377 tmp2 = pow(tmp1, -n);
378 nom += tmp * tmp2 * i;
382 for (i = SUMORDER; i > 0; i--)
384 tmp = -sin(2.0 * M_PI * i * K * rcoord);
385 tmp1 = 2.0 * M_PI * m / K + 2.0 * M_PI * i;
386 tmp2 = pow(tmp1, -n);
387 nom += tmp * tmp2 * i;
392 tmp = 2.0 * M_PI * m / K;
396 return 2.0 * M_PI * nom / denom * K;
402 /* The following routine is just a copy from pme.c */
404 static void calc_recipbox(matrix box, matrix recipbox)
406 /* Save some time by assuming upper right part is zero */
408 real tmp = 1.0/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
410 recipbox[XX][XX] = box[YY][YY]*box[ZZ][ZZ]*tmp;
411 recipbox[XX][YY] = 0;
412 recipbox[XX][ZZ] = 0;
413 recipbox[YY][XX] = -box[YY][XX]*box[ZZ][ZZ]*tmp;
414 recipbox[YY][YY] = box[XX][XX]*box[ZZ][ZZ]*tmp;
415 recipbox[YY][ZZ] = 0;
416 recipbox[ZZ][XX] = (box[YY][XX]*box[ZZ][YY]-box[YY][YY]*box[ZZ][XX])*tmp;
417 recipbox[ZZ][YY] = -box[ZZ][YY]*box[XX][XX]*tmp;
418 recipbox[ZZ][ZZ] = box[XX][XX]*box[YY][YY]*tmp;
422 /* Estimate the reciprocal space part error of the SPME Ewald sum. */
423 static real estimate_reciprocal(
425 rvec x[], /* array of particles */
426 real q[], /* array of charges */
427 int nr, /* number of charges = size of the charge array */
428 FILE gmx_unused *fp_out,
430 unsigned int seed, /* The seed for the random number generator */
431 int *nsamples, /* Return the number of samples used if Monte Carlo
432 * algorithm is used for self energy error estimate */
435 real e_rec = 0; /* reciprocal error estimate */
436 real e_rec1 = 0; /* Error estimate term 1*/
437 real e_rec2 = 0; /* Error estimate term 2*/
438 real e_rec3 = 0; /* Error estimate term 3 */
439 real e_rec3x = 0; /* part of Error estimate term 3 in x */
440 real e_rec3y = 0; /* part of Error estimate term 3 in y */
441 real e_rec3z = 0; /* part of Error estimate term 3 in z */
443 int nx, ny, nz; /* grid coordinates */
444 real q2_all = 0; /* sum of squared charges */
445 rvec gridpx; /* reciprocal grid point in x direction*/
446 rvec gridpxy; /* reciprocal grid point in x and y direction*/
447 rvec gridp; /* complete reciprocal grid point in 3 directions*/
448 rvec tmpvec; /* template to create points from basis vectors */
449 rvec tmpvec2; /* template to create points from basis vectors */
450 real coeff = 0; /* variable to compute coefficients of the error estimate */
451 real coeff2 = 0; /* variable to compute coefficients of the error estimate */
452 real tmp = 0; /* variables to compute different factors from vectors */
457 /* Random number generator */
458 gmx_rng_t rng = NULL;
461 /* Index variables for parallel work distribution */
462 int startglobal, stopglobal;
463 int startlocal, stoplocal;
472 rng = gmx_rng_init(seed);
480 for (i = 0; i < nr; i++)
485 /* Calculate indices for work distribution */
486 startglobal = -info->nkx[0]/2;
487 stopglobal = info->nkx[0]/2;
488 xtot = stopglobal*2+1;
491 x_per_core = static_cast<int>(ceil(static_cast<real>(xtot) / cr->nnodes));
492 startlocal = startglobal + x_per_core*cr->nodeid;
493 stoplocal = startlocal + x_per_core -1;
494 if (stoplocal > stopglobal)
496 stoplocal = stopglobal;
501 startlocal = startglobal;
502 stoplocal = stopglobal;
507 MPI_Barrier(MPI_COMM_WORLD);
523 fprintf(stderr, "Calculating reciprocal error part 1 ...");
527 for (nx = startlocal; nx <= stoplocal; nx++)
529 svmul(nx, info->recipbox[XX], gridpx);
530 for (ny = -info->nky[0]/2; ny < info->nky[0]/2+1; ny++)
532 svmul(ny, info->recipbox[YY], tmpvec);
533 rvec_add(gridpx, tmpvec, gridpxy);
534 for (nz = -info->nkz[0]/2; nz < info->nkz[0]/2+1; nz++)
536 if (0 == nx && 0 == ny && 0 == nz)
540 svmul(nz, info->recipbox[ZZ], tmpvec);
541 rvec_add(gridpxy, tmpvec, gridp);
543 coeff = exp(-1.0 * M_PI * M_PI * tmp / info->ewald_beta[0] / info->ewald_beta[0] );
544 coeff /= 2.0 * M_PI * info->volume * tmp;
548 tmp = eps_poly2(nx, info->nkx[0], info->pme_order[0]);
549 tmp += eps_poly2(ny, info->nkx[0], info->pme_order[0]);
550 tmp += eps_poly2(nz, info->nkx[0], info->pme_order[0]);
552 tmp1 = eps_poly1(nx, info->nkx[0], info->pme_order[0]);
553 tmp2 = eps_poly1(ny, info->nky[0], info->pme_order[0]);
555 tmp += 2.0 * tmp1 * tmp2;
557 tmp1 = eps_poly1(nz, info->nkz[0], info->pme_order[0]);
558 tmp2 = eps_poly1(ny, info->nky[0], info->pme_order[0]);
560 tmp += 2.0 * tmp1 * tmp2;
562 tmp1 = eps_poly1(nz, info->nkz[0], info->pme_order[0]);
563 tmp2 = eps_poly1(nx, info->nkx[0], info->pme_order[0]);
565 tmp += 2.0 * tmp1 * tmp2;
567 tmp1 = eps_poly1(nx, info->nkx[0], info->pme_order[0]);
568 tmp1 += eps_poly1(ny, info->nky[0], info->pme_order[0]);
569 tmp1 += eps_poly1(nz, info->nkz[0], info->pme_order[0]);
573 e_rec1 += 32.0 * M_PI * M_PI * coeff * coeff * coeff2 * tmp * q2_all * q2_all / nr;
575 tmp1 = eps_poly3(nx, info->nkx[0], info->pme_order[0]);
576 tmp1 *= info->nkx[0];
577 tmp2 = iprod(gridp, info->recipbox[XX]);
581 tmp1 = eps_poly3(ny, info->nky[0], info->pme_order[0]);
582 tmp1 *= info->nky[0];
583 tmp2 = iprod(gridp, info->recipbox[YY]);
587 tmp1 = eps_poly3(nz, info->nkz[0], info->pme_order[0]);
588 tmp1 *= info->nkz[0];
589 tmp2 = iprod(gridp, info->recipbox[ZZ]);
595 tmp1 = eps_poly4(nx, info->nkx[0], info->pme_order[0]);
596 tmp1 *= norm2(info->recipbox[XX]);
597 tmp1 *= info->nkx[0] * info->nkx[0];
601 tmp1 = eps_poly4(ny, info->nky[0], info->pme_order[0]);
602 tmp1 *= norm2(info->recipbox[YY]);
603 tmp1 *= info->nky[0] * info->nky[0];
607 tmp1 = eps_poly4(nz, info->nkz[0], info->pme_order[0]);
608 tmp1 *= norm2(info->recipbox[ZZ]);
609 tmp1 *= info->nkz[0] * info->nkz[0];
613 e_rec2 += 4.0 * coeff * coeff * tmp * q2_all * q2_all / nr;
619 fprintf(stderr, "\rCalculating reciprocal error part 1 ... %3.0f%%", 100.0*(nx-startlocal+1)/(x_per_core));
626 fprintf(stderr, "\n");
629 /* Use just a fraction of all charges to estimate the self energy error term? */
630 bFraction = (info->fracself > 0.0) && (info->fracself < 1.0);
634 /* Here xtot is the number of samples taken for the Monte Carlo calculation
635 * of the average of term IV of equation 35 in Wang2010. Round up to a
636 * number of samples that is divisible by the number of nodes */
637 x_per_core = static_cast<int>(ceil(info->fracself * nr / cr->nnodes));
638 xtot = x_per_core * cr->nnodes;
642 /* In this case we use all nr particle positions */
644 x_per_core = static_cast<int>(ceil(static_cast<real>(xtot) / cr->nnodes));
647 startlocal = x_per_core * cr->nodeid;
648 stoplocal = std::min(startlocal + x_per_core, xtot); /* min needed if xtot == nr */
652 /* Make shure we get identical results in serial and parallel. Therefore,
653 * take the sample indices from a single, global random number array that
654 * is constructed on the master node and that only depends on the seed */
658 for (i = 0; i < xtot; i++)
660 numbers[i] = static_cast<int>(floor(gmx_rng_uniform_real(rng) * nr));
663 /* Broadcast the random number array to the other nodes */
666 nblock_bc(cr, xtot, numbers);
669 if (bVerbose && MASTER(cr))
671 fprintf(stdout, "Using %d sample%s to approximate the self interaction error term",
672 xtot, xtot == 1 ? "" : "s");
675 fprintf(stdout, " (%d sample%s per node)", x_per_core, x_per_core == 1 ? "" : "s");
677 fprintf(stdout, ".\n");
681 /* Return the number of positions used for the Monte Carlo algorithm */
684 for (i = startlocal; i < stoplocal; i++)
692 /* Randomly pick a charge */
697 /* Use all charges */
701 /* for(nx=startlocal; nx<=stoplocal; nx++)*/
702 for (nx = -info->nkx[0]/2; nx < info->nkx[0]/2+1; nx++)
704 svmul(nx, info->recipbox[XX], gridpx);
705 for (ny = -info->nky[0]/2; ny < info->nky[0]/2+1; ny++)
707 svmul(ny, info->recipbox[YY], tmpvec);
708 rvec_add(gridpx, tmpvec, gridpxy);
709 for (nz = -info->nkz[0]/2; nz < info->nkz[0]/2+1; nz++)
712 if (0 == nx && 0 == ny && 0 == nz)
717 svmul(nz, info->recipbox[ZZ], tmpvec);
718 rvec_add(gridpxy, tmpvec, gridp);
720 coeff = exp(-1.0 * M_PI * M_PI * tmp / info->ewald_beta[0] / info->ewald_beta[0] );
722 e_rec3x += coeff*eps_self(nx, info->nkx[0], info->recipbox[XX], info->pme_order[0], x[ci]);
723 e_rec3y += coeff*eps_self(ny, info->nky[0], info->recipbox[YY], info->pme_order[0], x[ci]);
724 e_rec3z += coeff*eps_self(nz, info->nkz[0], info->recipbox[ZZ], info->pme_order[0], x[ci]);
732 svmul(e_rec3x, info->recipbox[XX], tmpvec);
733 rvec_inc(tmpvec2, tmpvec);
734 svmul(e_rec3y, info->recipbox[YY], tmpvec);
735 rvec_inc(tmpvec2, tmpvec);
736 svmul(e_rec3z, info->recipbox[ZZ], tmpvec);
737 rvec_inc(tmpvec2, tmpvec);
739 e_rec3 += q[ci]*q[ci]*q[ci]*q[ci]*norm2(tmpvec2) / ( xtot * M_PI * info->volume * M_PI * info->volume);
742 fprintf(stderr, "\rCalculating reciprocal error part 2 ... %3.0f%%",
743 100.0*(i+1)/stoplocal);
750 fprintf(stderr, "\n");
757 t1 = MPI_Wtime() - t0;
758 fprintf(fp_out, "Recip. err. est. took : %lf s\n", t1);
766 fprintf(stderr, "Node %3d: nx=[%3d...%3d] e_rec3=%e\n",
767 cr->nodeid, startlocal, stoplocal, e_rec3);
773 gmx_sum(1, &e_rec1, cr);
774 gmx_sum(1, &e_rec2, cr);
775 gmx_sum(1, &e_rec3, cr);
778 /* e_rec1*=8.0 * q2_all / info->volume / info->volume / nr ;
779 e_rec2*= q2_all / M_PI / M_PI / info->volume / info->volume / nr ;
780 e_rec3/= M_PI * M_PI * info->volume * info->volume * nr ;
782 e_rec = sqrt(e_rec1+e_rec2+e_rec3);
785 return ONE_4PI_EPS0 * e_rec;
789 /* Allocate memory for the inputinfo struct: */
790 static void create_info(t_inputinfo *info)
792 snew(info->fac, info->n_entries);
793 snew(info->rcoulomb, info->n_entries);
794 snew(info->rvdw, info->n_entries);
795 snew(info->nkx, info->n_entries);
796 snew(info->nky, info->n_entries);
797 snew(info->nkz, info->n_entries);
798 snew(info->fourier_sp, info->n_entries);
799 snew(info->ewald_rtol, info->n_entries);
800 snew(info->ewald_beta, info->n_entries);
801 snew(info->pme_order, info->n_entries);
802 snew(info->fn_out, info->n_entries);
803 snew(info->e_dir, info->n_entries);
804 snew(info->e_rec, info->n_entries);
808 /* Allocate and fill an array with coordinates and charges,
809 * returns the number of charges found
811 static int prepare_x_q(real *q[], rvec *x[], gmx_mtop_t *mtop, rvec x_orig[], t_commrec *cr)
814 int nq; /* number of charged particles */
815 gmx_mtop_atomloop_all_t aloop;
821 snew(*q, mtop->natoms);
822 snew(*x, mtop->natoms);
825 aloop = gmx_mtop_atomloop_all_init(mtop);
827 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
829 if (is_charge(atom->q))
832 (*x)[nq][XX] = x_orig[i][XX];
833 (*x)[nq][YY] = x_orig[i][YY];
834 (*x)[nq][ZZ] = x_orig[i][ZZ];
838 /* Give back some unneeded memory */
842 /* Broadcast x and q in the parallel case */
845 /* Transfer the number of charges */
849 nblock_bc(cr, nq, *x);
850 nblock_bc(cr, nq, *q);
858 /* Read in the tpr file and save information we need later in info */
859 static void read_tpr_file(const char *fn_sim_tpr, t_inputinfo *info, t_state *state, gmx_mtop_t *mtop, t_inputrec *ir, real user_beta, real fracself)
861 read_tpx_state(fn_sim_tpr, ir, state, NULL, mtop);
863 /* The values of the original tpr input file are save in the first
864 * place [0] of the arrays */
865 info->orig_sim_steps = ir->nsteps;
866 info->pme_order[0] = ir->pme_order;
867 info->rcoulomb[0] = ir->rcoulomb;
868 info->rvdw[0] = ir->rvdw;
869 info->nkx[0] = ir->nkx;
870 info->nky[0] = ir->nky;
871 info->nkz[0] = ir->nkz;
872 info->ewald_rtol[0] = ir->ewald_rtol;
873 info->fracself = fracself;
876 info->ewald_beta[0] = user_beta;
880 info->ewald_beta[0] = calc_ewaldcoeff_q(info->rcoulomb[0], info->ewald_rtol[0]);
883 /* Check if PME was chosen */
884 if (EEL_PME(ir->coulombtype) == FALSE)
886 gmx_fatal(FARGS, "Can only do optimizations for simulations with PME");
889 /* Check if rcoulomb == rlist, which is necessary for PME */
890 if (!(ir->rcoulomb == ir->rlist))
892 gmx_fatal(FARGS, "PME requires rcoulomb (%f) to be equal to rlist (%f).", ir->rcoulomb, ir->rlist);
897 /* Transfer what we need for parallelizing the reciprocal error estimate */
898 static void bcast_info(t_inputinfo *info, t_commrec *cr)
900 nblock_bc(cr, info->n_entries, info->nkx);
901 nblock_bc(cr, info->n_entries, info->nky);
902 nblock_bc(cr, info->n_entries, info->nkz);
903 nblock_bc(cr, info->n_entries, info->ewald_beta);
904 nblock_bc(cr, info->n_entries, info->pme_order);
905 nblock_bc(cr, info->n_entries, info->e_dir);
906 nblock_bc(cr, info->n_entries, info->e_rec);
907 block_bc(cr, info->volume);
908 block_bc(cr, info->recipbox);
909 block_bc(cr, info->natoms);
910 block_bc(cr, info->fracself);
911 block_bc(cr, info->bTUNE);
912 block_bc(cr, info->q2all);
913 block_bc(cr, info->q2allnr);
917 /* Estimate the error of the SPME Ewald sum. This estimate is based upon
918 * a) a homogeneous distribution of the charges
919 * b) a total charge of zero.
921 static void estimate_PME_error(t_inputinfo *info, t_state *state,
922 gmx_mtop_t *mtop, FILE *fp_out, gmx_bool bVerbose, unsigned int seed,
925 rvec *x = NULL; /* The coordinates */
926 real *q = NULL; /* 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, cr);
947 calc_q2all(mtop, &(info->q2all), &(info->q2allnr));
948 info->ewald_rtol[0] = gmx_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",
957 info->nkx[0], info->nky[0], info->nkz[0]);
964 bcast_info(info, cr);
968 /* Calculate direct space error */
969 info->e_dir[0] = estimate_direct(info);
971 /* Calculate reciprocal space error */
972 info->e_rec[0] = estimate_reciprocal(info, x, q, ncharges, fp_out, bVerbose,
973 seed, &nsamples, cr);
977 bcast_info(info, cr);
982 fprintf(fp_out, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
983 fprintf(fp_out, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
984 fprintf(fp_out, "Self-energy error term was estimated using %d samples\n", nsamples);
986 fprintf(stderr, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
987 fprintf(stderr, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
996 fprintf(stderr, "Starting tuning ...\n");
998 edir = info->e_dir[0];
999 erec = info->e_rec[0];
1001 beta0 = info->ewald_beta[0];
1004 info->ewald_beta[0] += 0.1;
1008 info->ewald_beta[0] -= 0.1;
1010 info->e_dir[0] = estimate_direct(info);
1011 info->e_rec[0] = estimate_reciprocal(info, x, q, ncharges, fp_out, bVerbose,
1012 seed, &nsamples, cr);
1016 bcast_info(info, cr);
1020 edir = info->e_dir[0];
1021 erec = info->e_rec[0];
1023 while (fabs(derr/std::min(erec, edir)) > 1e-4)
1026 beta = info->ewald_beta[0];
1027 beta -= derr*(info->ewald_beta[0]-beta0)/(derr-derr0);
1028 beta0 = info->ewald_beta[0];
1029 info->ewald_beta[0] = beta;
1032 info->e_dir[0] = estimate_direct(info);
1033 info->e_rec[0] = estimate_reciprocal(info, x, q, ncharges, fp_out, bVerbose,
1034 seed, &nsamples, cr);
1038 bcast_info(info, cr);
1041 edir = info->e_dir[0];
1042 erec = info->e_rec[0];
1048 fprintf(stderr, "difference between real and rec. space error (step %d): %g\n", i, fabs(derr));
1049 fprintf(stderr, "old beta: %f\n", beta0);
1050 fprintf(stderr, "new beta: %f\n", beta);
1054 info->ewald_rtol[0] = gmx_erfc(info->rcoulomb[0]*info->ewald_beta[0]);
1058 /* Write some info to log file */
1060 fprintf(fp_out, "========= After tuning ========\n");
1061 fprintf(fp_out, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
1062 fprintf(fp_out, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
1063 fprintf(stderr, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
1064 fprintf(stderr, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
1065 fprintf(fp_out, "Ewald_rtol : %g\n", info->ewald_rtol[0]);
1066 fprintf(fp_out, "Ewald parameter beta : %g\n", info->ewald_beta[0]);
1076 int gmx_pme_error(int argc, char *argv[])
1078 const char *desc[] = {
1079 "[THISMODULE] estimates the error of the electrostatic forces",
1080 "if using the sPME algorithm. The flag [TT]-tune[tt] will determine",
1081 "the splitting parameter such that the error is equally",
1082 "distributed over the real and reciprocal space part.",
1083 "The part of the error that stems from self interaction of the particles "
1084 "is computationally demanding. However, a good a approximation is to",
1085 "just use a fraction of the particles for this term which can be",
1086 "indicated by the flag [TT]-self[tt].[PAR]",
1089 real fs = 0.0; /* 0 indicates: not set by the user */
1090 real user_beta = -1.0;
1091 real fracself = 1.0;
1093 t_state state; /* The state from the tpr input file */
1094 gmx_mtop_t mtop; /* The topology from the tpr input file */
1095 t_inputrec *ir = NULL; /* The inputrec from the tpr file */
1098 unsigned long PCA_Flags;
1099 gmx_bool bTUNE = FALSE;
1100 gmx_bool bVerbose = FALSE;
1104 static t_filenm fnm[] = {
1105 { efTPX, "-s", NULL, ffREAD },
1106 { efOUT, "-o", "error", ffWRITE },
1107 { efTPX, "-so", "tuned", ffOPTWR }
1110 output_env_t oenv = NULL;
1113 { "-beta", FALSE, etREAL, {&user_beta},
1114 "If positive, overwrite ewald_beta from [TT].tpr[tt] file with this value" },
1115 { "-tune", FALSE, etBOOL, {&bTUNE},
1116 "Tune the splitting parameter such that the error is equally distributed between real and reciprocal space" },
1117 { "-self", FALSE, etREAL, {&fracself},
1118 "If between 0.0 and 1.0, determine self interaction error from just this fraction of the charged particles" },
1119 { "-seed", FALSE, etINT, {&seed},
1120 "Random number seed used for Monte Carlo algorithm when [TT]-self[tt] is set to a value between 0.0 and 1.0" },
1121 { "-v", FALSE, etBOOL, {&bVerbose},
1122 "Be loud and noisy" }
1126 #define NFILE asize(fnm)
1128 cr = init_commrec();
1130 MPI_Barrier(MPI_COMM_WORLD);
1133 PCA_Flags = PCA_NOEXIT_ON_ARGS;
1134 PCA_Flags |= (MASTER(cr) ? 0 : PCA_QUIET);
1136 if (!parse_common_args(&argc, argv, PCA_Flags,
1137 NFILE, fnm, asize(pa), pa, asize(desc), desc,
1145 bTUNE = opt2bSet("-so", NFILE, fnm);
1150 /* Allocate memory for the inputinfo struct: */
1152 info.fourier_sp[0] = fs;
1154 /* Read in the tpr file and open logfile for reading */
1158 read_tpr_file(opt2fn("-s", NFILE, fnm), &info, &state, &mtop, ir, user_beta, fracself);
1160 fp = fopen(opt2fn("-o", NFILE, fnm), "w");
1163 /* Check consistency if the user provided fourierspacing */
1164 if (fs > 0 && MASTER(cr))
1166 /* Recalculate the grid dimensions using fourierspacing from user input */
1170 calc_grid(stdout, state.box, info.fourier_sp[0], &(info.nkx[0]), &(info.nky[0]), &(info.nkz[0]));
1171 if ( (ir->nkx != info.nkx[0]) || (ir->nky != info.nky[0]) || (ir->nkz != info.nkz[0]) )
1173 gmx_fatal(FARGS, "Wrong fourierspacing %f nm, input file grid = %d x %d x %d, computed grid = %d x %d x %d",
1174 fs, ir->nkx, ir->nky, ir->nkz, info.nkx[0], info.nky[0], info.nkz[0]);
1178 /* Estimate (S)PME force error */
1180 /* Determine the volume of the simulation box */
1183 info.volume = det(state.box);
1184 calc_recipbox(state.box, info.recipbox);
1185 info.natoms = mtop.natoms;
1191 bcast_info(&info, cr);
1194 /* Get an error estimate of the input tpr file and do some tuning if requested */
1195 estimate_PME_error(&info, &state, &mtop, fp, bVerbose, seed, cr);
1199 /* Write out optimized tpr file if requested */
1200 if (opt2bSet("-so", NFILE, fnm) || bTUNE)
1202 ir->ewald_rtol = info.ewald_rtol[0];
1203 write_tpx_state(opt2fn("-so", NFILE, fnm), ir, &state, &mtop);
1205 please_cite(fp, "Wang2010");