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.
39 #include "gromacs/commandline/pargs.h"
40 #include "gromacs/legacyheaders/typedefs.h"
41 #include "gromacs/legacyheaders/types/commrec.h"
42 #include "gromacs/utility/smalloc.h"
43 #include "gromacs/math/vec.h"
44 #include "gromacs/legacyheaders/copyrite.h"
45 #include "gromacs/fileio/tpxio.h"
46 #include "gromacs/legacyheaders/readinp.h"
47 #include "gromacs/legacyheaders/calcgrid.h"
48 #include "gromacs/legacyheaders/checkpoint.h"
50 #include "gromacs/random/random.h"
51 #include "gromacs/math/units.h"
52 #include "gromacs/legacyheaders/mdatoms.h"
53 #include "gromacs/legacyheaders/coulomb.h"
54 #include "gromacs/topology/mtop_util.h"
55 #include "gromacs/legacyheaders/network.h"
56 #include "gromacs/legacyheaders/main.h"
57 #include "gromacs/legacyheaders/macros.h"
59 #include "gromacs/utility/fatalerror.h"
63 /* We use the same defines as in mvdata.c here */
64 #define block_bc(cr, d) gmx_bcast( sizeof(d), &(d), (cr))
65 #define nblock_bc(cr, nr, d) gmx_bcast((nr)*sizeof((d)[0]), (d), (cr))
66 #define snew_bc(cr, d, nr) { if (!MASTER(cr)) {snew((d), (nr)); }}
67 /* #define TAKETIME */
70 ddnoSEL, ddnoINTERLEAVE, ddnoPP_PME, ddnoCARTESIAN, ddnoNR
73 /* Enum for situations that can occur during log file parsing */
79 eParselogResetProblem,
86 gmx_int64_t orig_sim_steps; /* Number of steps to be done in the real simulation */
87 int n_entries; /* Number of entries in arrays */
88 real volume; /* The volume of the box */
89 matrix recipbox; /* The reciprocal box */
90 int natoms; /* The number of atoms in the MD system */
91 real *fac; /* The scaling factor */
92 real *rcoulomb; /* The coulomb radii [0...nr_inputfiles] */
93 real *rvdw; /* The vdW radii */
94 int *nkx, *nky, *nkz; /* Number of k vectors in each spatial dimension */
95 real *fourier_sp; /* Fourierspacing */
96 real *ewald_rtol; /* Real space tolerance for Ewald, determines */
97 /* the real/reciprocal space relative weight */
98 real *ewald_beta; /* Splitting parameter [1/nm] */
99 real fracself; /* fraction of particles for SI error */
100 real q2all; /* sum ( q ^2 ) */
101 real q2allnr; /* nr of charges */
102 int *pme_order; /* Interpolation order for PME (bsplines) */
103 char **fn_out; /* Name of the output tpr file */
104 real *e_dir; /* Direct space part of PME error with these settings */
105 real *e_rec; /* Reciprocal space part of PME error */
106 gmx_bool bTUNE; /* flag for tuning */
110 /* Returns TRUE when atom is charged */
111 static gmx_bool is_charge(real charge)
113 if (charge*charge > GMX_REAL_EPS)
124 /* calculate charge density */
125 static void calc_q2all(
126 gmx_mtop_t *mtop, /* molecular topology */
127 real *q2all, real *q2allnr)
129 int imol, iatom; /* indices for loops */
130 real q2_all = 0; /* Sum of squared charges */
131 int nrq_mol; /* Number of charges in a single molecule */
132 int nrq_all; /* Total number of charges in the MD system */
134 gmx_moltype_t *molecule;
135 gmx_molblock_t *molblock;
138 fprintf(stderr, "\nCharge density:\n");
140 q2_all = 0.0; /* total q squared */
141 nrq_all = 0; /* total number of charges in the system */
142 for (imol = 0; imol < mtop->nmolblock; imol++) /* Loop over molecule types */
144 q2_mol = 0.0; /* q squared value of this molecule */
145 nrq_mol = 0; /* number of charges this molecule carries */
146 molecule = &(mtop->moltype[imol]);
147 molblock = &(mtop->molblock[imol]);
148 for (iatom = 0; iatom < molblock->natoms_mol; iatom++) /* Loop over atoms in this molecule */
150 qi = molecule->atoms.atom[iatom].q; /* Charge of this atom */
151 /* Is this charge worth to be considered? */
158 /* Multiply with the number of molecules present of this type and add */
159 q2_all += q2_mol*molblock->nmol;
160 nrq_all += nrq_mol*molblock->nmol;
162 fprintf(stderr, "Molecule %2d (%5d atoms) q2_mol=%10.3e nr.mol.charges=%5d (%6dx) q2_all=%10.3e tot.charges=%d\n",
163 imol, molblock->natoms_mol, q2_mol, nrq_mol, molblock->nmol, q2_all, nrq_all);
173 /* Estimate the direct space part error of the SPME Ewald sum */
174 static real estimate_direct(
178 real e_dir = 0; /* Error estimate */
179 real beta = 0; /* Splitting parameter (1/nm) */
180 real r_coulomb = 0; /* Cut-off in direct space */
183 beta = info->ewald_beta[0];
184 r_coulomb = info->rcoulomb[0];
186 e_dir = 2.0 * info->q2all * gmx_invsqrt( info->q2allnr * r_coulomb * info->volume );
187 e_dir *= exp (-beta*beta*r_coulomb*r_coulomb);
189 return ONE_4PI_EPS0*e_dir;
194 /* the following 4 functions determine polynomials required for the reciprocal error estimate */
196 static inline real eps_poly1(
197 real m, /* grid coordinate in certain direction */
198 real K, /* grid size in corresponding direction */
199 real n) /* spline interpolation order of the SPME */
202 real nom = 0; /* nominator */
203 real denom = 0; /* denominator */
211 for (i = -SUMORDER; i < 0; i++)
215 nom += pow( tmp, -n );
218 for (i = SUMORDER; i > 0; i--)
222 nom += pow( tmp, -n );
227 denom = pow( tmp, -n )+nom;
233 static inline real eps_poly2(
234 real m, /* grid coordinate in certain direction */
235 real K, /* grid size in corresponding direction */
236 real n) /* spline interpolation order of the SPME */
239 real nom = 0; /* nominator */
240 real denom = 0; /* denominator */
248 for (i = -SUMORDER; i < 0; i++)
252 nom += pow( tmp, -2*n );
255 for (i = SUMORDER; i > 0; i--)
259 nom += pow( tmp, -2*n );
262 for (i = -SUMORDER; i < SUMORDER+1; i++)
266 denom += pow( tmp, -n );
268 tmp = eps_poly1(m, K, n);
269 return nom / denom / denom + tmp*tmp;
273 static inline real eps_poly3(
274 real m, /* grid coordinate in certain direction */
275 real K, /* grid size in corresponding direction */
276 real n) /* spline interpolation order of the SPME */
279 real nom = 0; /* nominator */
280 real denom = 0; /* denominator */
288 for (i = -SUMORDER; i < 0; i++)
292 nom += i * pow( tmp, -2*n );
295 for (i = SUMORDER; i > 0; i--)
299 nom += i * pow( tmp, -2*n );
302 for (i = -SUMORDER; i < SUMORDER+1; i++)
306 denom += pow( tmp, -n );
309 return 2.0 * M_PI * nom / denom / denom;
313 static inline real eps_poly4(
314 real m, /* grid coordinate in certain direction */
315 real K, /* grid size in corresponding direction */
316 real n) /* spline interpolation order of the SPME */
319 real nom = 0; /* nominator */
320 real denom = 0; /* denominator */
328 for (i = -SUMORDER; i < 0; i++)
332 nom += i * i * pow( tmp, -2*n );
335 for (i = SUMORDER; i > 0; i--)
339 nom += i * i * pow( tmp, -2*n );
342 for (i = -SUMORDER; i < SUMORDER+1; i++)
346 denom += pow( tmp, -n );
349 return 4.0 * M_PI * M_PI * nom / denom / denom;
353 static inline real eps_self(
354 real m, /* grid coordinate in certain direction */
355 real K, /* grid size in corresponding direction */
356 rvec rboxv, /* reciprocal box vector */
357 real n, /* spline interpolation order of the SPME */
358 rvec x) /* coordinate of charge */
361 real tmp = 0; /* temporary variables for computations */
362 real tmp1 = 0; /* temporary variables for computations */
363 real tmp2 = 0; /* temporary variables for computations */
364 real rcoord = 0; /* coordinate in certain reciprocal space direction */
365 real nom = 0; /* nominator */
366 real denom = 0; /* denominator */
374 rcoord = iprod(rboxv, x);
377 for (i = -SUMORDER; i < 0; i++)
379 tmp = -sin(2.0 * M_PI * i * K * rcoord);
380 tmp1 = 2.0 * M_PI * m / K + 2.0 * M_PI * i;
381 tmp2 = pow(tmp1, -n);
382 nom += tmp * tmp2 * i;
386 for (i = SUMORDER; i > 0; i--)
388 tmp = -sin(2.0 * M_PI * i * K * rcoord);
389 tmp1 = 2.0 * M_PI * m / K + 2.0 * M_PI * i;
390 tmp2 = pow(tmp1, -n);
391 nom += tmp * tmp2 * i;
396 tmp = 2.0 * M_PI * m / K;
400 return 2.0 * M_PI * nom / denom * K;
406 /* The following routine is just a copy from pme.c */
408 static void calc_recipbox(matrix box, matrix recipbox)
410 /* Save some time by assuming upper right part is zero */
412 real tmp = 1.0/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
414 recipbox[XX][XX] = box[YY][YY]*box[ZZ][ZZ]*tmp;
415 recipbox[XX][YY] = 0;
416 recipbox[XX][ZZ] = 0;
417 recipbox[YY][XX] = -box[YY][XX]*box[ZZ][ZZ]*tmp;
418 recipbox[YY][YY] = box[XX][XX]*box[ZZ][ZZ]*tmp;
419 recipbox[YY][ZZ] = 0;
420 recipbox[ZZ][XX] = (box[YY][XX]*box[ZZ][YY]-box[YY][YY]*box[ZZ][XX])*tmp;
421 recipbox[ZZ][YY] = -box[ZZ][YY]*box[XX][XX]*tmp;
422 recipbox[ZZ][ZZ] = box[XX][XX]*box[YY][YY]*tmp;
426 /* Estimate the reciprocal space part error of the SPME Ewald sum. */
427 static real estimate_reciprocal(
429 rvec x[], /* array of particles */
430 real q[], /* array of charges */
431 int nr, /* number of charges = size of the charge array */
432 FILE gmx_unused *fp_out,
434 unsigned int seed, /* The seed for the random number generator */
435 int *nsamples, /* Return the number of samples used if Monte Carlo
436 * algorithm is used for self energy error estimate */
439 real e_rec = 0; /* reciprocal error estimate */
440 real e_rec1 = 0; /* Error estimate term 1*/
441 real e_rec2 = 0; /* Error estimate term 2*/
442 real e_rec3 = 0; /* Error estimate term 3 */
443 real e_rec3x = 0; /* part of Error estimate term 3 in x */
444 real e_rec3y = 0; /* part of Error estimate term 3 in y */
445 real e_rec3z = 0; /* part of Error estimate term 3 in z */
447 int nx, ny, nz; /* grid coordinates */
448 real q2_all = 0; /* sum of squared charges */
449 rvec gridpx; /* reciprocal grid point in x direction*/
450 rvec gridpxy; /* reciprocal grid point in x and y direction*/
451 rvec gridp; /* complete reciprocal grid point in 3 directions*/
452 rvec tmpvec; /* template to create points from basis vectors */
453 rvec tmpvec2; /* template to create points from basis vectors */
454 real coeff = 0; /* variable to compute coefficients of the error estimate */
455 real coeff2 = 0; /* variable to compute coefficients of the error estimate */
456 real tmp = 0; /* variables to compute different factors from vectors */
461 /* Random number generator */
462 gmx_rng_t rng = NULL;
465 /* Index variables for parallel work distribution */
466 int startglobal, stopglobal;
467 int startlocal, stoplocal;
476 rng = gmx_rng_init(seed);
484 for (i = 0; i < nr; i++)
489 /* Calculate indices for work distribution */
490 startglobal = -info->nkx[0]/2;
491 stopglobal = info->nkx[0]/2;
492 xtot = stopglobal*2+1;
495 x_per_core = static_cast<int>(ceil(static_cast<real>(xtot) / cr->nnodes));
496 startlocal = startglobal + x_per_core*cr->nodeid;
497 stoplocal = startlocal + x_per_core -1;
498 if (stoplocal > stopglobal)
500 stoplocal = stopglobal;
505 startlocal = startglobal;
506 stoplocal = stopglobal;
511 MPI_Barrier(MPI_COMM_WORLD);
527 fprintf(stderr, "Calculating reciprocal error part 1 ...");
531 for (nx = startlocal; nx <= stoplocal; nx++)
533 svmul(nx, info->recipbox[XX], gridpx);
534 for (ny = -info->nky[0]/2; ny < info->nky[0]/2+1; ny++)
536 svmul(ny, info->recipbox[YY], tmpvec);
537 rvec_add(gridpx, tmpvec, gridpxy);
538 for (nz = -info->nkz[0]/2; nz < info->nkz[0]/2+1; nz++)
540 if (0 == nx && 0 == ny && 0 == nz)
544 svmul(nz, info->recipbox[ZZ], tmpvec);
545 rvec_add(gridpxy, tmpvec, gridp);
547 coeff = exp(-1.0 * M_PI * M_PI * tmp / info->ewald_beta[0] / info->ewald_beta[0] );
548 coeff /= 2.0 * M_PI * info->volume * tmp;
552 tmp = eps_poly2(nx, info->nkx[0], info->pme_order[0]);
553 tmp += eps_poly2(ny, info->nkx[0], info->pme_order[0]);
554 tmp += eps_poly2(nz, info->nkx[0], info->pme_order[0]);
556 tmp1 = eps_poly1(nx, info->nkx[0], info->pme_order[0]);
557 tmp2 = eps_poly1(ny, info->nky[0], info->pme_order[0]);
559 tmp += 2.0 * tmp1 * tmp2;
561 tmp1 = eps_poly1(nz, info->nkz[0], info->pme_order[0]);
562 tmp2 = eps_poly1(ny, info->nky[0], info->pme_order[0]);
564 tmp += 2.0 * tmp1 * tmp2;
566 tmp1 = eps_poly1(nz, info->nkz[0], info->pme_order[0]);
567 tmp2 = eps_poly1(nx, info->nkx[0], info->pme_order[0]);
569 tmp += 2.0 * tmp1 * tmp2;
571 tmp1 = eps_poly1(nx, info->nkx[0], info->pme_order[0]);
572 tmp1 += eps_poly1(ny, info->nky[0], info->pme_order[0]);
573 tmp1 += eps_poly1(nz, info->nkz[0], info->pme_order[0]);
577 e_rec1 += 32.0 * M_PI * M_PI * coeff * coeff * coeff2 * tmp * q2_all * q2_all / nr;
579 tmp1 = eps_poly3(nx, info->nkx[0], info->pme_order[0]);
580 tmp1 *= info->nkx[0];
581 tmp2 = iprod(gridp, info->recipbox[XX]);
585 tmp1 = eps_poly3(ny, info->nky[0], info->pme_order[0]);
586 tmp1 *= info->nky[0];
587 tmp2 = iprod(gridp, info->recipbox[YY]);
591 tmp1 = eps_poly3(nz, info->nkz[0], info->pme_order[0]);
592 tmp1 *= info->nkz[0];
593 tmp2 = iprod(gridp, info->recipbox[ZZ]);
599 tmp1 = eps_poly4(nx, info->nkx[0], info->pme_order[0]);
600 tmp1 *= norm2(info->recipbox[XX]);
601 tmp1 *= info->nkx[0] * info->nkx[0];
605 tmp1 = eps_poly4(ny, info->nky[0], info->pme_order[0]);
606 tmp1 *= norm2(info->recipbox[YY]);
607 tmp1 *= info->nky[0] * info->nky[0];
611 tmp1 = eps_poly4(nz, info->nkz[0], info->pme_order[0]);
612 tmp1 *= norm2(info->recipbox[ZZ]);
613 tmp1 *= info->nkz[0] * info->nkz[0];
617 e_rec2 += 4.0 * coeff * coeff * tmp * q2_all * q2_all / nr;
623 fprintf(stderr, "\rCalculating reciprocal error part 1 ... %3.0f%%", 100.0*(nx-startlocal+1)/(x_per_core));
630 fprintf(stderr, "\n");
633 /* Use just a fraction of all charges to estimate the self energy error term? */
634 bFraction = (info->fracself > 0.0) && (info->fracself < 1.0);
638 /* Here xtot is the number of samples taken for the Monte Carlo calculation
639 * of the average of term IV of equation 35 in Wang2010. Round up to a
640 * number of samples that is divisible by the number of nodes */
641 x_per_core = static_cast<int>(ceil(info->fracself * nr / cr->nnodes));
642 xtot = x_per_core * cr->nnodes;
646 /* In this case we use all nr particle positions */
648 x_per_core = static_cast<int>(ceil(static_cast<real>(xtot) / cr->nnodes));
651 startlocal = x_per_core * cr->nodeid;
652 stoplocal = std::min(startlocal + x_per_core, xtot); /* min needed if xtot == nr */
656 /* Make shure we get identical results in serial and parallel. Therefore,
657 * take the sample indices from a single, global random number array that
658 * is constructed on the master node and that only depends on the seed */
662 for (i = 0; i < xtot; i++)
664 numbers[i] = static_cast<int>(floor(gmx_rng_uniform_real(rng) * nr));
667 /* Broadcast the random number array to the other nodes */
670 nblock_bc(cr, xtot, numbers);
673 if (bVerbose && MASTER(cr))
675 fprintf(stdout, "Using %d sample%s to approximate the self interaction error term",
676 xtot, xtot == 1 ? "" : "s");
679 fprintf(stdout, " (%d sample%s per rank)", x_per_core, x_per_core == 1 ? "" : "s");
681 fprintf(stdout, ".\n");
685 /* Return the number of positions used for the Monte Carlo algorithm */
688 for (i = startlocal; i < stoplocal; i++)
696 /* Randomly pick a charge */
701 /* Use all charges */
705 /* for(nx=startlocal; nx<=stoplocal; nx++)*/
706 for (nx = -info->nkx[0]/2; nx < info->nkx[0]/2+1; nx++)
708 svmul(nx, info->recipbox[XX], gridpx);
709 for (ny = -info->nky[0]/2; ny < info->nky[0]/2+1; ny++)
711 svmul(ny, info->recipbox[YY], tmpvec);
712 rvec_add(gridpx, tmpvec, gridpxy);
713 for (nz = -info->nkz[0]/2; nz < info->nkz[0]/2+1; nz++)
716 if (0 == nx && 0 == ny && 0 == nz)
721 svmul(nz, info->recipbox[ZZ], tmpvec);
722 rvec_add(gridpxy, tmpvec, gridp);
724 coeff = exp(-1.0 * M_PI * M_PI * tmp / info->ewald_beta[0] / info->ewald_beta[0] );
726 e_rec3x += coeff*eps_self(nx, info->nkx[0], info->recipbox[XX], info->pme_order[0], x[ci]);
727 e_rec3y += coeff*eps_self(ny, info->nky[0], info->recipbox[YY], info->pme_order[0], x[ci]);
728 e_rec3z += coeff*eps_self(nz, info->nkz[0], info->recipbox[ZZ], info->pme_order[0], x[ci]);
736 svmul(e_rec3x, info->recipbox[XX], tmpvec);
737 rvec_inc(tmpvec2, tmpvec);
738 svmul(e_rec3y, info->recipbox[YY], tmpvec);
739 rvec_inc(tmpvec2, tmpvec);
740 svmul(e_rec3z, info->recipbox[ZZ], tmpvec);
741 rvec_inc(tmpvec2, tmpvec);
743 e_rec3 += q[ci]*q[ci]*q[ci]*q[ci]*norm2(tmpvec2) / ( xtot * M_PI * info->volume * M_PI * info->volume);
746 fprintf(stderr, "\rCalculating reciprocal error part 2 ... %3.0f%%",
747 100.0*(i+1)/stoplocal);
754 fprintf(stderr, "\n");
761 t1 = MPI_Wtime() - t0;
762 fprintf(fp_out, "Recip. err. est. took : %lf s\n", t1);
770 fprintf(stderr, "Rank %3d: nx=[%3d...%3d] e_rec3=%e\n",
771 cr->nodeid, startlocal, stoplocal, e_rec3);
777 gmx_sum(1, &e_rec1, cr);
778 gmx_sum(1, &e_rec2, cr);
779 gmx_sum(1, &e_rec3, cr);
782 /* e_rec1*=8.0 * q2_all / info->volume / info->volume / nr ;
783 e_rec2*= q2_all / M_PI / M_PI / info->volume / info->volume / nr ;
784 e_rec3/= M_PI * M_PI * info->volume * info->volume * nr ;
786 e_rec = sqrt(e_rec1+e_rec2+e_rec3);
789 return ONE_4PI_EPS0 * e_rec;
793 /* Allocate memory for the inputinfo struct: */
794 static void create_info(t_inputinfo *info)
796 snew(info->fac, info->n_entries);
797 snew(info->rcoulomb, info->n_entries);
798 snew(info->rvdw, info->n_entries);
799 snew(info->nkx, info->n_entries);
800 snew(info->nky, info->n_entries);
801 snew(info->nkz, info->n_entries);
802 snew(info->fourier_sp, info->n_entries);
803 snew(info->ewald_rtol, info->n_entries);
804 snew(info->ewald_beta, info->n_entries);
805 snew(info->pme_order, info->n_entries);
806 snew(info->fn_out, info->n_entries);
807 snew(info->e_dir, info->n_entries);
808 snew(info->e_rec, info->n_entries);
812 /* Allocate and fill an array with coordinates and charges,
813 * returns the number of charges found
815 static int prepare_x_q(real *q[], rvec *x[], gmx_mtop_t *mtop, rvec x_orig[], t_commrec *cr)
818 int nq; /* number of charged particles */
819 gmx_mtop_atomloop_all_t aloop;
825 snew(*q, mtop->natoms);
826 snew(*x, mtop->natoms);
829 aloop = gmx_mtop_atomloop_all_init(mtop);
831 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
833 if (is_charge(atom->q))
836 (*x)[nq][XX] = x_orig[i][XX];
837 (*x)[nq][YY] = x_orig[i][YY];
838 (*x)[nq][ZZ] = x_orig[i][ZZ];
842 /* Give back some unneeded memory */
846 /* Broadcast x and q in the parallel case */
849 /* Transfer the number of charges */
853 nblock_bc(cr, nq, *x);
854 nblock_bc(cr, nq, *q);
862 /* Read in the tpr file and save information we need later in info */
863 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)
865 read_tpx_state(fn_sim_tpr, ir, state, NULL, mtop);
867 /* The values of the original tpr input file are save in the first
868 * place [0] of the arrays */
869 info->orig_sim_steps = ir->nsteps;
870 info->pme_order[0] = ir->pme_order;
871 info->rcoulomb[0] = ir->rcoulomb;
872 info->rvdw[0] = ir->rvdw;
873 info->nkx[0] = ir->nkx;
874 info->nky[0] = ir->nky;
875 info->nkz[0] = ir->nkz;
876 info->ewald_rtol[0] = ir->ewald_rtol;
877 info->fracself = fracself;
880 info->ewald_beta[0] = user_beta;
884 info->ewald_beta[0] = calc_ewaldcoeff_q(info->rcoulomb[0], info->ewald_rtol[0]);
887 /* Check if PME was chosen */
888 if (EEL_PME(ir->coulombtype) == FALSE)
890 gmx_fatal(FARGS, "Can only do optimizations for simulations with PME");
893 /* Check if rcoulomb == rlist, which is necessary for PME */
894 if (!(ir->rcoulomb == ir->rlist))
896 gmx_fatal(FARGS, "PME requires rcoulomb (%f) to be equal to rlist (%f).", ir->rcoulomb, ir->rlist);
901 /* Transfer what we need for parallelizing the reciprocal error estimate */
902 static void bcast_info(t_inputinfo *info, t_commrec *cr)
904 nblock_bc(cr, info->n_entries, info->nkx);
905 nblock_bc(cr, info->n_entries, info->nky);
906 nblock_bc(cr, info->n_entries, info->nkz);
907 nblock_bc(cr, info->n_entries, info->ewald_beta);
908 nblock_bc(cr, info->n_entries, info->pme_order);
909 nblock_bc(cr, info->n_entries, info->e_dir);
910 nblock_bc(cr, info->n_entries, info->e_rec);
911 block_bc(cr, info->volume);
912 block_bc(cr, info->recipbox);
913 block_bc(cr, info->natoms);
914 block_bc(cr, info->fracself);
915 block_bc(cr, info->bTUNE);
916 block_bc(cr, info->q2all);
917 block_bc(cr, info->q2allnr);
921 /* Estimate the error of the SPME Ewald sum. This estimate is based upon
922 * a) a homogeneous distribution of the charges
923 * b) a total charge of zero.
925 static void estimate_PME_error(t_inputinfo *info, t_state *state,
926 gmx_mtop_t *mtop, FILE *fp_out, gmx_bool bVerbose, unsigned int seed,
929 rvec *x = NULL; /* The coordinates */
930 real *q = NULL; /* The charges */
931 real edir = 0.0; /* real space error */
932 real erec = 0.0; /* reciprocal space error */
933 real derr = 0.0; /* difference of real and reciprocal space error */
934 real derr0 = 0.0; /* difference of real and reciprocal space error */
935 real beta = 0.0; /* splitting parameter beta */
936 real beta0 = 0.0; /* splitting parameter beta */
937 int ncharges; /* The number of atoms with charges */
938 int nsamples; /* The number of samples used for the calculation of the
939 * self-energy error term */
944 fprintf(fp_out, "\n--- PME ERROR ESTIMATE ---\n");
947 /* Prepare an x and q array with only the charged atoms */
948 ncharges = prepare_x_q(&q, &x, mtop, state->x, cr);
951 calc_q2all(mtop, &(info->q2all), &(info->q2allnr));
952 info->ewald_rtol[0] = gmx_erfc(info->rcoulomb[0]*info->ewald_beta[0]);
953 /* Write some info to log file */
954 fprintf(fp_out, "Box volume : %g nm^3\n", info->volume);
955 fprintf(fp_out, "Number of charged atoms : %d (total atoms %d)\n", ncharges, info->natoms);
956 fprintf(fp_out, "Coulomb radius : %g nm\n", info->rcoulomb[0]);
957 fprintf(fp_out, "Ewald_rtol : %g\n", info->ewald_rtol[0]);
958 fprintf(fp_out, "Ewald parameter beta : %g\n", info->ewald_beta[0]);
959 fprintf(fp_out, "Interpolation order : %d\n", info->pme_order[0]);
960 fprintf(fp_out, "Fourier grid (nx,ny,nz) : %d x %d x %d\n",
961 info->nkx[0], info->nky[0], info->nkz[0]);
968 bcast_info(info, cr);
972 /* Calculate direct space error */
973 info->e_dir[0] = estimate_direct(info);
975 /* Calculate reciprocal space error */
976 info->e_rec[0] = estimate_reciprocal(info, x, q, ncharges, fp_out, bVerbose,
977 seed, &nsamples, cr);
981 bcast_info(info, cr);
986 fprintf(fp_out, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
987 fprintf(fp_out, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
988 fprintf(fp_out, "Self-energy error term was estimated using %d samples\n", nsamples);
990 fprintf(stderr, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
991 fprintf(stderr, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
1000 fprintf(stderr, "Starting tuning ...\n");
1002 edir = info->e_dir[0];
1003 erec = info->e_rec[0];
1005 beta0 = info->ewald_beta[0];
1008 info->ewald_beta[0] += 0.1;
1012 info->ewald_beta[0] -= 0.1;
1014 info->e_dir[0] = estimate_direct(info);
1015 info->e_rec[0] = estimate_reciprocal(info, x, q, ncharges, fp_out, bVerbose,
1016 seed, &nsamples, cr);
1020 bcast_info(info, cr);
1024 edir = info->e_dir[0];
1025 erec = info->e_rec[0];
1027 while (fabs(derr/std::min(erec, edir)) > 1e-4)
1030 beta = info->ewald_beta[0];
1031 beta -= derr*(info->ewald_beta[0]-beta0)/(derr-derr0);
1032 beta0 = info->ewald_beta[0];
1033 info->ewald_beta[0] = beta;
1036 info->e_dir[0] = estimate_direct(info);
1037 info->e_rec[0] = estimate_reciprocal(info, x, q, ncharges, fp_out, bVerbose,
1038 seed, &nsamples, cr);
1042 bcast_info(info, cr);
1045 edir = info->e_dir[0];
1046 erec = info->e_rec[0];
1052 fprintf(stderr, "difference between real and rec. space error (step %d): %g\n", i, fabs(derr));
1053 fprintf(stderr, "old beta: %f\n", beta0);
1054 fprintf(stderr, "new beta: %f\n", beta);
1058 info->ewald_rtol[0] = gmx_erfc(info->rcoulomb[0]*info->ewald_beta[0]);
1062 /* Write some info to log file */
1064 fprintf(fp_out, "========= After tuning ========\n");
1065 fprintf(fp_out, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
1066 fprintf(fp_out, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
1067 fprintf(stderr, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
1068 fprintf(stderr, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
1069 fprintf(fp_out, "Ewald_rtol : %g\n", info->ewald_rtol[0]);
1070 fprintf(fp_out, "Ewald parameter beta : %g\n", info->ewald_beta[0]);
1080 int gmx_pme_error(int argc, char *argv[])
1082 const char *desc[] = {
1083 "[THISMODULE] estimates the error of the electrostatic forces",
1084 "if using the sPME algorithm. The flag [TT]-tune[tt] will determine",
1085 "the splitting parameter such that the error is equally",
1086 "distributed over the real and reciprocal space part.",
1087 "The part of the error that stems from self interaction of the particles "
1088 "is computationally demanding. However, a good a approximation is to",
1089 "just use a fraction of the particles for this term which can be",
1090 "indicated by the flag [TT]-self[tt].[PAR]",
1093 real fs = 0.0; /* 0 indicates: not set by the user */
1094 real user_beta = -1.0;
1095 real fracself = 1.0;
1097 t_state state; /* The state from the tpr input file */
1098 gmx_mtop_t mtop; /* The topology from the tpr input file */
1099 t_inputrec *ir = NULL; /* The inputrec from the tpr file */
1102 unsigned long PCA_Flags;
1103 gmx_bool bTUNE = FALSE;
1104 gmx_bool bVerbose = FALSE;
1108 static t_filenm fnm[] = {
1109 { efTPX, "-s", NULL, ffREAD },
1110 { efOUT, "-o", "error", ffWRITE },
1111 { efTPX, "-so", "tuned", ffOPTWR }
1114 output_env_t oenv = NULL;
1117 { "-beta", FALSE, etREAL, {&user_beta},
1118 "If positive, overwrite ewald_beta from [TT].tpr[tt] file with this value" },
1119 { "-tune", FALSE, etBOOL, {&bTUNE},
1120 "Tune the splitting parameter such that the error is equally distributed between real and reciprocal space" },
1121 { "-self", FALSE, etREAL, {&fracself},
1122 "If between 0.0 and 1.0, determine self interaction error from just this fraction of the charged particles" },
1123 { "-seed", FALSE, etINT, {&seed},
1124 "Random number seed used for Monte Carlo algorithm when [TT]-self[tt] is set to a value between 0.0 and 1.0" },
1125 { "-v", FALSE, etBOOL, {&bVerbose},
1126 "Be loud and noisy" }
1130 #define NFILE asize(fnm)
1132 cr = init_commrec();
1134 PCA_Flags = PCA_NOEXIT_ON_ARGS;
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");