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"
61 /* We use the same defines as in mvdata.c here */
62 #define block_bc(cr, d) gmx_bcast( sizeof(d), &(d), (cr))
63 #define nblock_bc(cr, nr, d) gmx_bcast((nr)*sizeof((d)[0]), (d), (cr))
64 #define snew_bc(cr, d, nr) { if (!MASTER(cr)) {snew((d), (nr)); }}
65 /* #define TAKETIME */
68 ddnoSEL, ddnoINTERLEAVE, ddnoPP_PME, ddnoCARTESIAN, ddnoNR
71 /* Enum for situations that can occur during log file parsing */
77 eParselogResetProblem,
84 gmx_int64_t orig_sim_steps; /* Number of steps to be done in the real simulation */
85 int n_entries; /* Number of entries in arrays */
86 real volume; /* The volume of the box */
87 matrix recipbox; /* The reciprocal box */
88 int natoms; /* The number of atoms in the MD system */
89 real *fac; /* The scaling factor */
90 real *rcoulomb; /* The coulomb radii [0...nr_inputfiles] */
91 real *rvdw; /* The vdW radii */
92 int *nkx, *nky, *nkz; /* Number of k vectors in each spatial dimension */
93 real *fourier_sp; /* Fourierspacing */
94 real *ewald_rtol; /* Real space tolerance for Ewald, determines */
95 /* the real/reciprocal space relative weight */
96 real *ewald_beta; /* Splitting parameter [1/nm] */
97 real fracself; /* fraction of particles for SI error */
98 real q2all; /* sum ( q ^2 ) */
99 real q2allnr; /* nr of charges */
100 int *pme_order; /* Interpolation order for PME (bsplines) */
101 char **fn_out; /* Name of the output tpr file */
102 real *e_dir; /* Direct space part of PME error with these settings */
103 real *e_rec; /* Reciprocal space part of PME error */
104 gmx_bool bTUNE; /* flag for tuning */
108 /* Returns TRUE when atom is charged */
109 static gmx_bool is_charge(real charge)
111 if (charge*charge > GMX_REAL_EPS)
122 /* calculate charge density */
123 static void calc_q2all(
124 gmx_mtop_t *mtop, /* molecular topology */
125 real *q2all, real *q2allnr)
127 int imol, iatom; /* indices for loops */
128 real q2_all = 0; /* Sum of squared charges */
129 int nrq_mol; /* Number of charges in a single molecule */
130 int nrq_all; /* Total number of charges in the MD system */
132 gmx_moltype_t *molecule;
133 gmx_molblock_t *molblock;
136 fprintf(stderr, "\nCharge density:\n");
138 q2_all = 0.0; /* total q squared */
139 nrq_all = 0; /* total number of charges in the system */
140 for (imol = 0; imol < mtop->nmolblock; imol++) /* Loop over molecule types */
142 q2_mol = 0.0; /* q squared value of this molecule */
143 nrq_mol = 0; /* number of charges this molecule carries */
144 molecule = &(mtop->moltype[imol]);
145 molblock = &(mtop->molblock[imol]);
146 for (iatom = 0; iatom < molblock->natoms_mol; iatom++) /* Loop over atoms in this molecule */
148 qi = molecule->atoms.atom[iatom].q; /* Charge of this atom */
149 /* Is this charge worth to be considered? */
156 /* Multiply with the number of molecules present of this type and add */
157 q2_all += q2_mol*molblock->nmol;
158 nrq_all += nrq_mol*molblock->nmol;
160 fprintf(stderr, "Molecule %2d (%5d atoms) q2_mol=%10.3e nr.mol.charges=%5d (%6dx) q2_all=%10.3e tot.charges=%d\n",
161 imol, molblock->natoms_mol, q2_mol, nrq_mol, molblock->nmol, q2_all, nrq_all);
171 /* Estimate the direct space part error of the SPME Ewald sum */
172 static real estimate_direct(
176 real e_dir = 0; /* Error estimate */
177 real beta = 0; /* Splitting parameter (1/nm) */
178 real r_coulomb = 0; /* Cut-off in direct space */
181 beta = info->ewald_beta[0];
182 r_coulomb = info->rcoulomb[0];
184 e_dir = 2.0 * info->q2all * gmx_invsqrt( info->q2allnr * r_coulomb * info->volume );
185 e_dir *= exp (-beta*beta*r_coulomb*r_coulomb);
187 return ONE_4PI_EPS0*e_dir;
192 /* the following 4 functions determine polynomials required for the reciprocal error estimate */
194 static inline real eps_poly1(
195 real m, /* grid coordinate in certain direction */
196 real K, /* grid size in corresponding direction */
197 real n) /* spline interpolation order of the SPME */
200 real nom = 0; /* nominator */
201 real denom = 0; /* denominator */
209 for (i = -SUMORDER; i < 0; i++)
213 nom += pow( tmp, -n );
216 for (i = SUMORDER; i > 0; i--)
220 nom += pow( tmp, -n );
225 denom = pow( tmp, -n )+nom;
231 static inline real eps_poly2(
232 real m, /* grid coordinate in certain direction */
233 real K, /* grid size in corresponding direction */
234 real n) /* spline interpolation order of the SPME */
237 real nom = 0; /* nominator */
238 real denom = 0; /* denominator */
246 for (i = -SUMORDER; i < 0; i++)
250 nom += pow( tmp, -2*n );
253 for (i = SUMORDER; i > 0; i--)
257 nom += pow( tmp, -2*n );
260 for (i = -SUMORDER; i < SUMORDER+1; i++)
264 denom += pow( tmp, -n );
266 tmp = eps_poly1(m, K, n);
267 return nom / denom / denom + tmp*tmp;
271 static inline real eps_poly3(
272 real m, /* grid coordinate in certain direction */
273 real K, /* grid size in corresponding direction */
274 real n) /* spline interpolation order of the SPME */
277 real nom = 0; /* nominator */
278 real denom = 0; /* denominator */
286 for (i = -SUMORDER; i < 0; i++)
290 nom += i * pow( tmp, -2*n );
293 for (i = SUMORDER; i > 0; i--)
297 nom += i * pow( tmp, -2*n );
300 for (i = -SUMORDER; i < SUMORDER+1; i++)
304 denom += pow( tmp, -n );
307 return 2.0 * M_PI * nom / denom / denom;
311 static inline real eps_poly4(
312 real m, /* grid coordinate in certain direction */
313 real K, /* grid size in corresponding direction */
314 real n) /* spline interpolation order of the SPME */
317 real nom = 0; /* nominator */
318 real denom = 0; /* denominator */
326 for (i = -SUMORDER; i < 0; i++)
330 nom += i * i * pow( tmp, -2*n );
333 for (i = SUMORDER; i > 0; i--)
337 nom += i * i * pow( tmp, -2*n );
340 for (i = -SUMORDER; i < SUMORDER+1; i++)
344 denom += pow( tmp, -n );
347 return 4.0 * M_PI * M_PI * nom / denom / denom;
351 static inline real eps_self(
352 real m, /* grid coordinate in certain direction */
353 real K, /* grid size in corresponding direction */
354 rvec rboxv, /* reciprocal box vector */
355 real n, /* spline interpolation order of the SPME */
356 rvec x) /* coordinate of charge */
359 real tmp = 0; /* temporary variables for computations */
360 real tmp1 = 0; /* temporary variables for computations */
361 real tmp2 = 0; /* temporary variables for computations */
362 real rcoord = 0; /* coordinate in certain reciprocal space direction */
363 real nom = 0; /* nominator */
364 real denom = 0; /* denominator */
372 rcoord = iprod(rboxv, x);
375 for (i = -SUMORDER; i < 0; i++)
377 tmp = -sin(2.0 * M_PI * i * K * rcoord);
378 tmp1 = 2.0 * M_PI * m / K + 2.0 * M_PI * i;
379 tmp2 = pow(tmp1, -n);
380 nom += tmp * tmp2 * i;
384 for (i = SUMORDER; i > 0; i--)
386 tmp = -sin(2.0 * M_PI * i * K * rcoord);
387 tmp1 = 2.0 * M_PI * m / K + 2.0 * M_PI * i;
388 tmp2 = pow(tmp1, -n);
389 nom += tmp * tmp2 * i;
394 tmp = 2.0 * M_PI * m / K;
398 return 2.0 * M_PI * nom / denom * K;
404 /* The following routine is just a copy from pme.c */
406 static void calc_recipbox(matrix box, matrix recipbox)
408 /* Save some time by assuming upper right part is zero */
410 real tmp = 1.0/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
412 recipbox[XX][XX] = box[YY][YY]*box[ZZ][ZZ]*tmp;
413 recipbox[XX][YY] = 0;
414 recipbox[XX][ZZ] = 0;
415 recipbox[YY][XX] = -box[YY][XX]*box[ZZ][ZZ]*tmp;
416 recipbox[YY][YY] = box[XX][XX]*box[ZZ][ZZ]*tmp;
417 recipbox[YY][ZZ] = 0;
418 recipbox[ZZ][XX] = (box[YY][XX]*box[ZZ][YY]-box[YY][YY]*box[ZZ][XX])*tmp;
419 recipbox[ZZ][YY] = -box[ZZ][YY]*box[XX][XX]*tmp;
420 recipbox[ZZ][ZZ] = box[XX][XX]*box[YY][YY]*tmp;
424 /* Estimate the reciprocal space part error of the SPME Ewald sum. */
425 static real estimate_reciprocal(
427 rvec x[], /* array of particles */
428 real q[], /* array of charges */
429 int nr, /* number of charges = size of the charge array */
430 FILE gmx_unused *fp_out,
432 unsigned int seed, /* The seed for the random number generator */
433 int *nsamples, /* Return the number of samples used if Monte Carlo
434 * algorithm is used for self energy error estimate */
437 real e_rec = 0; /* reciprocal error estimate */
438 real e_rec1 = 0; /* Error estimate term 1*/
439 real e_rec2 = 0; /* Error estimate term 2*/
440 real e_rec3 = 0; /* Error estimate term 3 */
441 real e_rec3x = 0; /* part of Error estimate term 3 in x */
442 real e_rec3y = 0; /* part of Error estimate term 3 in y */
443 real e_rec3z = 0; /* part of Error estimate term 3 in z */
445 int nx, ny, nz; /* grid coordinates */
446 real q2_all = 0; /* sum of squared charges */
447 rvec gridpx; /* reciprocal grid point in x direction*/
448 rvec gridpxy; /* reciprocal grid point in x and y direction*/
449 rvec gridp; /* complete reciprocal grid point in 3 directions*/
450 rvec tmpvec; /* template to create points from basis vectors */
451 rvec tmpvec2; /* template to create points from basis vectors */
452 real coeff = 0; /* variable to compute coefficients of the error estimate */
453 real coeff2 = 0; /* variable to compute coefficients of the error estimate */
454 real tmp = 0; /* variables to compute different factors from vectors */
459 /* Random number generator */
460 gmx_rng_t rng = NULL;
463 /* Index variables for parallel work distribution */
464 int startglobal, stopglobal;
465 int startlocal, stoplocal;
474 rng = gmx_rng_init(seed);
482 for (i = 0; i < nr; i++)
487 /* Calculate indices for work distribution */
488 startglobal = -info->nkx[0]/2;
489 stopglobal = info->nkx[0]/2;
490 xtot = stopglobal*2+1;
493 x_per_core = static_cast<int>(ceil(static_cast<real>(xtot) / cr->nnodes));
494 startlocal = startglobal + x_per_core*cr->nodeid;
495 stoplocal = startlocal + x_per_core -1;
496 if (stoplocal > stopglobal)
498 stoplocal = stopglobal;
503 startlocal = startglobal;
504 stoplocal = stopglobal;
509 MPI_Barrier(MPI_COMM_WORLD);
525 fprintf(stderr, "Calculating reciprocal error part 1 ...");
529 for (nx = startlocal; nx <= stoplocal; nx++)
531 svmul(nx, info->recipbox[XX], gridpx);
532 for (ny = -info->nky[0]/2; ny < info->nky[0]/2+1; ny++)
534 svmul(ny, info->recipbox[YY], tmpvec);
535 rvec_add(gridpx, tmpvec, gridpxy);
536 for (nz = -info->nkz[0]/2; nz < info->nkz[0]/2+1; nz++)
538 if (0 == nx && 0 == ny && 0 == nz)
542 svmul(nz, info->recipbox[ZZ], tmpvec);
543 rvec_add(gridpxy, tmpvec, gridp);
545 coeff = exp(-1.0 * M_PI * M_PI * tmp / info->ewald_beta[0] / info->ewald_beta[0] );
546 coeff /= 2.0 * M_PI * info->volume * tmp;
550 tmp = eps_poly2(nx, info->nkx[0], info->pme_order[0]);
551 tmp += eps_poly2(ny, info->nkx[0], info->pme_order[0]);
552 tmp += eps_poly2(nz, info->nkx[0], info->pme_order[0]);
554 tmp1 = eps_poly1(nx, info->nkx[0], info->pme_order[0]);
555 tmp2 = eps_poly1(ny, info->nky[0], info->pme_order[0]);
557 tmp += 2.0 * tmp1 * tmp2;
559 tmp1 = eps_poly1(nz, info->nkz[0], info->pme_order[0]);
560 tmp2 = eps_poly1(ny, info->nky[0], info->pme_order[0]);
562 tmp += 2.0 * tmp1 * tmp2;
564 tmp1 = eps_poly1(nz, info->nkz[0], info->pme_order[0]);
565 tmp2 = eps_poly1(nx, info->nkx[0], info->pme_order[0]);
567 tmp += 2.0 * tmp1 * tmp2;
569 tmp1 = eps_poly1(nx, info->nkx[0], info->pme_order[0]);
570 tmp1 += eps_poly1(ny, info->nky[0], info->pme_order[0]);
571 tmp1 += eps_poly1(nz, info->nkz[0], info->pme_order[0]);
575 e_rec1 += 32.0 * M_PI * M_PI * coeff * coeff * coeff2 * tmp * q2_all * q2_all / nr;
577 tmp1 = eps_poly3(nx, info->nkx[0], info->pme_order[0]);
578 tmp1 *= info->nkx[0];
579 tmp2 = iprod(gridp, info->recipbox[XX]);
583 tmp1 = eps_poly3(ny, info->nky[0], info->pme_order[0]);
584 tmp1 *= info->nky[0];
585 tmp2 = iprod(gridp, info->recipbox[YY]);
589 tmp1 = eps_poly3(nz, info->nkz[0], info->pme_order[0]);
590 tmp1 *= info->nkz[0];
591 tmp2 = iprod(gridp, info->recipbox[ZZ]);
597 tmp1 = eps_poly4(nx, info->nkx[0], info->pme_order[0]);
598 tmp1 *= norm2(info->recipbox[XX]);
599 tmp1 *= info->nkx[0] * info->nkx[0];
603 tmp1 = eps_poly4(ny, info->nky[0], info->pme_order[0]);
604 tmp1 *= norm2(info->recipbox[YY]);
605 tmp1 *= info->nky[0] * info->nky[0];
609 tmp1 = eps_poly4(nz, info->nkz[0], info->pme_order[0]);
610 tmp1 *= norm2(info->recipbox[ZZ]);
611 tmp1 *= info->nkz[0] * info->nkz[0];
615 e_rec2 += 4.0 * coeff * coeff * tmp * q2_all * q2_all / nr;
621 fprintf(stderr, "\rCalculating reciprocal error part 1 ... %3.0f%%", 100.0*(nx-startlocal+1)/(x_per_core));
628 fprintf(stderr, "\n");
631 /* Use just a fraction of all charges to estimate the self energy error term? */
632 bFraction = (info->fracself > 0.0) && (info->fracself < 1.0);
636 /* Here xtot is the number of samples taken for the Monte Carlo calculation
637 * of the average of term IV of equation 35 in Wang2010. Round up to a
638 * number of samples that is divisible by the number of nodes */
639 x_per_core = static_cast<int>(ceil(info->fracself * nr / cr->nnodes));
640 xtot = x_per_core * cr->nnodes;
644 /* In this case we use all nr particle positions */
646 x_per_core = static_cast<int>(ceil(static_cast<real>(xtot) / cr->nnodes));
649 startlocal = x_per_core * cr->nodeid;
650 stoplocal = std::min(startlocal + x_per_core, xtot); /* min needed if xtot == nr */
654 /* Make shure we get identical results in serial and parallel. Therefore,
655 * take the sample indices from a single, global random number array that
656 * is constructed on the master node and that only depends on the seed */
660 for (i = 0; i < xtot; i++)
662 numbers[i] = static_cast<int>(floor(gmx_rng_uniform_real(rng) * nr));
665 /* Broadcast the random number array to the other nodes */
668 nblock_bc(cr, xtot, numbers);
671 if (bVerbose && MASTER(cr))
673 fprintf(stdout, "Using %d sample%s to approximate the self interaction error term",
674 xtot, xtot == 1 ? "" : "s");
677 fprintf(stdout, " (%d sample%s per rank)", x_per_core, x_per_core == 1 ? "" : "s");
679 fprintf(stdout, ".\n");
683 /* Return the number of positions used for the Monte Carlo algorithm */
686 for (i = startlocal; i < stoplocal; i++)
694 /* Randomly pick a charge */
699 /* Use all charges */
703 /* for(nx=startlocal; nx<=stoplocal; nx++)*/
704 for (nx = -info->nkx[0]/2; nx < info->nkx[0]/2+1; nx++)
706 svmul(nx, info->recipbox[XX], gridpx);
707 for (ny = -info->nky[0]/2; ny < info->nky[0]/2+1; ny++)
709 svmul(ny, info->recipbox[YY], tmpvec);
710 rvec_add(gridpx, tmpvec, gridpxy);
711 for (nz = -info->nkz[0]/2; nz < info->nkz[0]/2+1; nz++)
714 if (0 == nx && 0 == ny && 0 == nz)
719 svmul(nz, info->recipbox[ZZ], tmpvec);
720 rvec_add(gridpxy, tmpvec, gridp);
722 coeff = exp(-1.0 * M_PI * M_PI * tmp / info->ewald_beta[0] / info->ewald_beta[0] );
724 e_rec3x += coeff*eps_self(nx, info->nkx[0], info->recipbox[XX], info->pme_order[0], x[ci]);
725 e_rec3y += coeff*eps_self(ny, info->nky[0], info->recipbox[YY], info->pme_order[0], x[ci]);
726 e_rec3z += coeff*eps_self(nz, info->nkz[0], info->recipbox[ZZ], info->pme_order[0], x[ci]);
734 svmul(e_rec3x, info->recipbox[XX], tmpvec);
735 rvec_inc(tmpvec2, tmpvec);
736 svmul(e_rec3y, info->recipbox[YY], tmpvec);
737 rvec_inc(tmpvec2, tmpvec);
738 svmul(e_rec3z, info->recipbox[ZZ], tmpvec);
739 rvec_inc(tmpvec2, tmpvec);
741 e_rec3 += q[ci]*q[ci]*q[ci]*q[ci]*norm2(tmpvec2) / ( xtot * M_PI * info->volume * M_PI * info->volume);
744 fprintf(stderr, "\rCalculating reciprocal error part 2 ... %3.0f%%",
745 100.0*(i+1)/stoplocal);
752 fprintf(stderr, "\n");
759 t1 = MPI_Wtime() - t0;
760 fprintf(fp_out, "Recip. err. est. took : %lf s\n", t1);
768 fprintf(stderr, "Rank %3d: nx=[%3d...%3d] e_rec3=%e\n",
769 cr->nodeid, startlocal, stoplocal, e_rec3);
775 gmx_sum(1, &e_rec1, cr);
776 gmx_sum(1, &e_rec2, cr);
777 gmx_sum(1, &e_rec3, cr);
780 /* e_rec1*=8.0 * q2_all / info->volume / info->volume / nr ;
781 e_rec2*= q2_all / M_PI / M_PI / info->volume / info->volume / nr ;
782 e_rec3/= M_PI * M_PI * info->volume * info->volume * nr ;
784 e_rec = sqrt(e_rec1+e_rec2+e_rec3);
787 return ONE_4PI_EPS0 * e_rec;
791 /* Allocate memory for the inputinfo struct: */
792 static void create_info(t_inputinfo *info)
794 snew(info->fac, info->n_entries);
795 snew(info->rcoulomb, info->n_entries);
796 snew(info->rvdw, info->n_entries);
797 snew(info->nkx, info->n_entries);
798 snew(info->nky, info->n_entries);
799 snew(info->nkz, info->n_entries);
800 snew(info->fourier_sp, info->n_entries);
801 snew(info->ewald_rtol, info->n_entries);
802 snew(info->ewald_beta, info->n_entries);
803 snew(info->pme_order, info->n_entries);
804 snew(info->fn_out, info->n_entries);
805 snew(info->e_dir, info->n_entries);
806 snew(info->e_rec, info->n_entries);
810 /* Allocate and fill an array with coordinates and charges,
811 * returns the number of charges found
813 static int prepare_x_q(real *q[], rvec *x[], gmx_mtop_t *mtop, rvec x_orig[], t_commrec *cr)
816 int nq; /* number of charged particles */
817 gmx_mtop_atomloop_all_t aloop;
823 snew(*q, mtop->natoms);
824 snew(*x, mtop->natoms);
827 aloop = gmx_mtop_atomloop_all_init(mtop);
829 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
831 if (is_charge(atom->q))
834 (*x)[nq][XX] = x_orig[i][XX];
835 (*x)[nq][YY] = x_orig[i][YY];
836 (*x)[nq][ZZ] = x_orig[i][ZZ];
840 /* Give back some unneeded memory */
844 /* Broadcast x and q in the parallel case */
847 /* Transfer the number of charges */
851 nblock_bc(cr, nq, *x);
852 nblock_bc(cr, nq, *q);
860 /* Read in the tpr file and save information we need later in info */
861 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)
863 read_tpx_state(fn_sim_tpr, ir, state, NULL, mtop);
865 /* The values of the original tpr input file are save in the first
866 * place [0] of the arrays */
867 info->orig_sim_steps = ir->nsteps;
868 info->pme_order[0] = ir->pme_order;
869 info->rcoulomb[0] = ir->rcoulomb;
870 info->rvdw[0] = ir->rvdw;
871 info->nkx[0] = ir->nkx;
872 info->nky[0] = ir->nky;
873 info->nkz[0] = ir->nkz;
874 info->ewald_rtol[0] = ir->ewald_rtol;
875 info->fracself = fracself;
878 info->ewald_beta[0] = user_beta;
882 info->ewald_beta[0] = calc_ewaldcoeff_q(info->rcoulomb[0], info->ewald_rtol[0]);
885 /* Check if PME was chosen */
886 if (EEL_PME(ir->coulombtype) == FALSE)
888 gmx_fatal(FARGS, "Can only do optimizations for simulations with PME");
891 /* Check if rcoulomb == rlist, which is necessary for PME */
892 if (!(ir->rcoulomb == ir->rlist))
894 gmx_fatal(FARGS, "PME requires rcoulomb (%f) to be equal to rlist (%f).", ir->rcoulomb, ir->rlist);
899 /* Transfer what we need for parallelizing the reciprocal error estimate */
900 static void bcast_info(t_inputinfo *info, t_commrec *cr)
902 nblock_bc(cr, info->n_entries, info->nkx);
903 nblock_bc(cr, info->n_entries, info->nky);
904 nblock_bc(cr, info->n_entries, info->nkz);
905 nblock_bc(cr, info->n_entries, info->ewald_beta);
906 nblock_bc(cr, info->n_entries, info->pme_order);
907 nblock_bc(cr, info->n_entries, info->e_dir);
908 nblock_bc(cr, info->n_entries, info->e_rec);
909 block_bc(cr, info->volume);
910 block_bc(cr, info->recipbox);
911 block_bc(cr, info->natoms);
912 block_bc(cr, info->fracself);
913 block_bc(cr, info->bTUNE);
914 block_bc(cr, info->q2all);
915 block_bc(cr, info->q2allnr);
919 /* Estimate the error of the SPME Ewald sum. This estimate is based upon
920 * a) a homogeneous distribution of the charges
921 * b) a total charge of zero.
923 static void estimate_PME_error(t_inputinfo *info, t_state *state,
924 gmx_mtop_t *mtop, FILE *fp_out, gmx_bool bVerbose, unsigned int seed,
927 rvec *x = NULL; /* The coordinates */
928 real *q = NULL; /* The charges */
929 real edir = 0.0; /* real space error */
930 real erec = 0.0; /* reciprocal space error */
931 real derr = 0.0; /* difference of real and reciprocal space error */
932 real derr0 = 0.0; /* difference of real and reciprocal space error */
933 real beta = 0.0; /* splitting parameter beta */
934 real beta0 = 0.0; /* splitting parameter beta */
935 int ncharges; /* The number of atoms with charges */
936 int nsamples; /* The number of samples used for the calculation of the
937 * self-energy error term */
942 fprintf(fp_out, "\n--- PME ERROR ESTIMATE ---\n");
945 /* Prepare an x and q array with only the charged atoms */
946 ncharges = prepare_x_q(&q, &x, mtop, state->x, cr);
949 calc_q2all(mtop, &(info->q2all), &(info->q2allnr));
950 info->ewald_rtol[0] = gmx_erfc(info->rcoulomb[0]*info->ewald_beta[0]);
951 /* Write some info to log file */
952 fprintf(fp_out, "Box volume : %g nm^3\n", info->volume);
953 fprintf(fp_out, "Number of charged atoms : %d (total atoms %d)\n", ncharges, info->natoms);
954 fprintf(fp_out, "Coulomb radius : %g nm\n", info->rcoulomb[0]);
955 fprintf(fp_out, "Ewald_rtol : %g\n", info->ewald_rtol[0]);
956 fprintf(fp_out, "Ewald parameter beta : %g\n", info->ewald_beta[0]);
957 fprintf(fp_out, "Interpolation order : %d\n", info->pme_order[0]);
958 fprintf(fp_out, "Fourier grid (nx,ny,nz) : %d x %d x %d\n",
959 info->nkx[0], info->nky[0], info->nkz[0]);
966 bcast_info(info, cr);
970 /* Calculate direct space error */
971 info->e_dir[0] = estimate_direct(info);
973 /* Calculate reciprocal space error */
974 info->e_rec[0] = estimate_reciprocal(info, x, q, ncharges, fp_out, bVerbose,
975 seed, &nsamples, cr);
979 bcast_info(info, cr);
984 fprintf(fp_out, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
985 fprintf(fp_out, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
986 fprintf(fp_out, "Self-energy error term was estimated using %d samples\n", nsamples);
988 fprintf(stderr, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
989 fprintf(stderr, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
998 fprintf(stderr, "Starting tuning ...\n");
1000 edir = info->e_dir[0];
1001 erec = info->e_rec[0];
1003 beta0 = info->ewald_beta[0];
1006 info->ewald_beta[0] += 0.1;
1010 info->ewald_beta[0] -= 0.1;
1012 info->e_dir[0] = estimate_direct(info);
1013 info->e_rec[0] = estimate_reciprocal(info, x, q, ncharges, fp_out, bVerbose,
1014 seed, &nsamples, cr);
1018 bcast_info(info, cr);
1022 edir = info->e_dir[0];
1023 erec = info->e_rec[0];
1025 while (fabs(derr/std::min(erec, edir)) > 1e-4)
1028 beta = info->ewald_beta[0];
1029 beta -= derr*(info->ewald_beta[0]-beta0)/(derr-derr0);
1030 beta0 = info->ewald_beta[0];
1031 info->ewald_beta[0] = beta;
1034 info->e_dir[0] = estimate_direct(info);
1035 info->e_rec[0] = estimate_reciprocal(info, x, q, ncharges, fp_out, bVerbose,
1036 seed, &nsamples, cr);
1040 bcast_info(info, cr);
1043 edir = info->e_dir[0];
1044 erec = info->e_rec[0];
1050 fprintf(stderr, "difference between real and rec. space error (step %d): %g\n", i, fabs(derr));
1051 fprintf(stderr, "old beta: %f\n", beta0);
1052 fprintf(stderr, "new beta: %f\n", beta);
1056 info->ewald_rtol[0] = gmx_erfc(info->rcoulomb[0]*info->ewald_beta[0]);
1060 /* Write some info to log file */
1062 fprintf(fp_out, "========= After tuning ========\n");
1063 fprintf(fp_out, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
1064 fprintf(fp_out, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
1065 fprintf(stderr, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
1066 fprintf(stderr, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
1067 fprintf(fp_out, "Ewald_rtol : %g\n", info->ewald_rtol[0]);
1068 fprintf(fp_out, "Ewald parameter beta : %g\n", info->ewald_beta[0]);
1078 int gmx_pme_error(int argc, char *argv[])
1080 const char *desc[] = {
1081 "[THISMODULE] estimates the error of the electrostatic forces",
1082 "if using the sPME algorithm. The flag [TT]-tune[tt] will determine",
1083 "the splitting parameter such that the error is equally",
1084 "distributed over the real and reciprocal space part.",
1085 "The part of the error that stems from self interaction of the particles "
1086 "is computationally demanding. However, a good a approximation is to",
1087 "just use a fraction of the particles for this term which can be",
1088 "indicated by the flag [TT]-self[tt].[PAR]",
1091 real fs = 0.0; /* 0 indicates: not set by the user */
1092 real user_beta = -1.0;
1093 real fracself = 1.0;
1095 t_state state; /* The state from the tpr input file */
1096 gmx_mtop_t mtop; /* The topology from the tpr input file */
1097 t_inputrec *ir = NULL; /* The inputrec from the tpr file */
1100 unsigned long PCA_Flags;
1101 gmx_bool bTUNE = FALSE;
1102 gmx_bool bVerbose = FALSE;
1106 static t_filenm fnm[] = {
1107 { efTPX, "-s", NULL, ffREAD },
1108 { efOUT, "-o", "error", ffWRITE },
1109 { efTPX, "-so", "tuned", ffOPTWR }
1112 output_env_t oenv = NULL;
1115 { "-beta", FALSE, etREAL, {&user_beta},
1116 "If positive, overwrite ewald_beta from [TT].tpr[tt] file with this value" },
1117 { "-tune", FALSE, etBOOL, {&bTUNE},
1118 "Tune the splitting parameter such that the error is equally distributed between real and reciprocal space" },
1119 { "-self", FALSE, etREAL, {&fracself},
1120 "If between 0.0 and 1.0, determine self interaction error from just this fraction of the charged particles" },
1121 { "-seed", FALSE, etINT, {&seed},
1122 "Random number seed used for Monte Carlo algorithm when [TT]-self[tt] is set to a value between 0.0 and 1.0" },
1123 { "-v", FALSE, etBOOL, {&bVerbose},
1124 "Be loud and noisy" }
1128 #define NFILE asize(fnm)
1130 cr = init_commrec();
1132 PCA_Flags = PCA_NOEXIT_ON_ARGS;
1134 if (!parse_common_args(&argc, argv, PCA_Flags,
1135 NFILE, fnm, asize(pa), pa, asize(desc), desc,
1143 bTUNE = opt2bSet("-so", NFILE, fnm);
1148 /* Allocate memory for the inputinfo struct: */
1150 info.fourier_sp[0] = fs;
1152 /* Read in the tpr file and open logfile for reading */
1156 read_tpr_file(opt2fn("-s", NFILE, fnm), &info, &state, &mtop, ir, user_beta, fracself);
1158 fp = fopen(opt2fn("-o", NFILE, fnm), "w");
1161 /* Check consistency if the user provided fourierspacing */
1162 if (fs > 0 && MASTER(cr))
1164 /* Recalculate the grid dimensions using fourierspacing from user input */
1168 calc_grid(stdout, state.box, info.fourier_sp[0], &(info.nkx[0]), &(info.nky[0]), &(info.nkz[0]));
1169 if ( (ir->nkx != info.nkx[0]) || (ir->nky != info.nky[0]) || (ir->nkz != info.nkz[0]) )
1171 gmx_fatal(FARGS, "Wrong fourierspacing %f nm, input file grid = %d x %d x %d, computed grid = %d x %d x %d",
1172 fs, ir->nkx, ir->nky, ir->nkz, info.nkx[0], info.nky[0], info.nkz[0]);
1176 /* Estimate (S)PME force error */
1178 /* Determine the volume of the simulation box */
1181 info.volume = det(state.box);
1182 calc_recipbox(state.box, info.recipbox);
1183 info.natoms = mtop.natoms;
1189 bcast_info(&info, cr);
1192 /* Get an error estimate of the input tpr file and do some tuning if requested */
1193 estimate_PME_error(&info, &state, &mtop, fp, bVerbose, seed, cr);
1197 /* Write out optimized tpr file if requested */
1198 if (opt2bSet("-so", NFILE, fnm) || bTUNE)
1200 ir->ewald_rtol = info.ewald_rtol[0];
1201 write_tpx_state(opt2fn("-so", NFILE, fnm), ir, &state, &mtop);
1203 please_cite(fp, "Wang2010");