2 * This source code is part of
6 * GROningen MAchine for Chemical Simulations
8 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
9 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
10 * Copyright (c) 2001-2008, The GROMACS development team,
11 * check out http://www.gromacs.org for more information.
13 * This program is free software; you can redistribute it and/or
14 * modify it under the terms of the GNU General Public License
15 * as published by the Free Software Foundation; either version 2
16 * of the License, or (at your option) any later version.
18 * If you want to redistribute modifications, please consider that
19 * scientific software is very special. Version control is crucial -
20 * bugs must be traceable. We will be happy to consider code for
21 * inclusion in the official distribution, but derived work must not
22 * be called official GROMACS. Details are found in the README & COPYING
23 * files - if they are missing, get the official version at www.gromacs.org.
25 * To help us fund GROMACS development, we humbly ask that you cite
26 * the papers on the package - you can find them in the top README file.
28 * For more info, check our website at http://www.gromacs.org
31 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
44 #include "checkpoint.h"
46 #include "gmx_random.h"
50 #include "mtop_util.h"
55 /* We use the same defines as in mvdata.c here */
56 #define block_bc(cr, d) gmx_bcast( sizeof(d), &(d), (cr))
57 #define nblock_bc(cr, nr, d) gmx_bcast((nr)*sizeof((d)[0]), (d), (cr))
58 #define snew_bc(cr, d, nr) { if (!MASTER(cr)) {snew((d), (nr)); }}
59 /* #define TAKETIME */
62 ddnoSEL, ddnoINTERLEAVE, ddnoPP_PME, ddnoCARTESIAN, ddnoNR
65 /* Enum for situations that can occur during log file parsing */
71 eParselogResetProblem,
78 gmx_large_int_t orig_sim_steps; /* Number of steps to be done in the real simulation */
79 int n_entries; /* Number of entries in arrays */
80 real volume; /* The volume of the box */
81 matrix recipbox; /* The reciprocal box */
82 int natoms; /* The number of atoms in the MD system */
83 real *fac; /* The scaling factor */
84 real *rcoulomb; /* The coulomb radii [0...nr_inputfiles] */
85 real *rvdw; /* The vdW radii */
86 int *nkx, *nky, *nkz; /* Number of k vectors in each spatial dimension */
87 real *fourier_sp; /* Fourierspacing */
88 real *ewald_rtol; /* Real space tolerance for Ewald, determines */
89 /* the real/reciprocal space relative weight */
90 real *ewald_beta; /* Splitting parameter [1/nm] */
91 real fracself; /* fraction of particles for SI error */
92 real q2all; /* sum ( q ^2 ) */
93 real q2allnr; /* nr of charges */
94 int *pme_order; /* Interpolation order for PME (bsplines) */
95 char **fn_out; /* Name of the output tpr file */
96 real *e_dir; /* Direct space part of PME error with these settings */
97 real *e_rec; /* Reciprocal space part of PME error */
98 gmx_bool bTUNE; /* flag for tuning */
102 /* Returns TRUE when atom is charged */
103 static gmx_bool is_charge(real charge)
105 if (charge*charge > GMX_REAL_EPS)
116 /* calculate charge density */
117 static void calc_q2all(
118 gmx_mtop_t *mtop, /* molecular topology */
119 real *q2all, real *q2allnr)
121 int imol, iatom; /* indices for loops */
122 real q2_all = 0; /* Sum of squared charges */
123 int nrq_mol; /* Number of charges in a single molecule */
124 int nrq_all; /* Total number of charges in the MD system */
126 gmx_moltype_t *molecule;
127 gmx_molblock_t *molblock;
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 (imol = 0; imol < mtop->nmolblock; imol++) /* 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 molecule = &(mtop->moltype[imol]);
139 molblock = &(mtop->molblock[imol]);
140 for (iatom = 0; iatom < molblock->natoms_mol; iatom++) /* Loop over atoms in this molecule */
142 qi = molecule->atoms.atom[iatom].q; /* Charge of this atom */
143 /* Is this charge worth to be considered? */
150 /* Multiply with the number of molecules present of this type and add */
151 q2_all += q2_mol*molblock->nmol;
152 nrq_all += nrq_mol*molblock->nmol;
154 fprintf(stderr, "Molecule %2d (%5d atoms) q2_mol=%10.3e nr.mol.charges=%5d (%6dx) q2_all=%10.3e tot.charges=%d\n",
155 imol, molblock->natoms_mol, q2_mol, nrq_mol, molblock->nmol, q2_all, nrq_all);
165 /* Estimate the direct space part error of the SPME Ewald sum */
166 static real estimate_direct(
170 real e_dir = 0; /* Error estimate */
171 real beta = 0; /* Splitting parameter (1/nm) */
172 real r_coulomb = 0; /* Cut-off in direct space */
175 beta = info->ewald_beta[0];
176 r_coulomb = info->rcoulomb[0];
178 e_dir = 2.0 * info->q2all * gmx_invsqrt( info->q2allnr * r_coulomb * info->volume );
179 e_dir *= exp (-beta*beta*r_coulomb*r_coulomb);
181 return ONE_4PI_EPS0*e_dir;
186 /* the following 4 functions determine polynomials required for the reciprocal error estimate */
188 static inline real eps_poly1(
189 real m, /* grid coordinate in certain direction */
190 real K, /* grid size in corresponding direction */
191 real n) /* spline interpolation order of the SPME */
194 real nom = 0; /* nominator */
195 real denom = 0; /* denominator */
203 for (i = -SUMORDER; i < 0; i++)
207 nom += pow( tmp, -n );
210 for (i = SUMORDER; i > 0; i--)
214 nom += pow( tmp, -n );
219 denom = pow( tmp, -n )+nom;
225 static inline real eps_poly2(
226 real m, /* grid coordinate in certain direction */
227 real K, /* grid size in corresponding direction */
228 real n) /* spline interpolation order of the SPME */
231 real nom = 0; /* nominator */
232 real denom = 0; /* denominator */
240 for (i = -SUMORDER; i < 0; i++)
244 nom += pow( tmp, -2*n );
247 for (i = SUMORDER; i > 0; i--)
251 nom += pow( tmp, -2*n );
254 for (i = -SUMORDER; i < SUMORDER+1; i++)
258 denom += pow( tmp, -n );
260 tmp = eps_poly1(m, K, n);
261 return nom / denom / denom + tmp*tmp;
265 static inline real eps_poly3(
266 real m, /* grid coordinate in certain direction */
267 real K, /* grid size in corresponding direction */
268 real n) /* spline interpolation order of the SPME */
271 real nom = 0; /* nominator */
272 real denom = 0; /* denominator */
280 for (i = -SUMORDER; i < 0; i++)
284 nom += i * pow( tmp, -2*n );
287 for (i = SUMORDER; i > 0; i--)
291 nom += i * pow( tmp, -2*n );
294 for (i = -SUMORDER; i < SUMORDER+1; i++)
298 denom += pow( tmp, -n );
301 return 2.0 * M_PI * nom / denom / denom;
305 static inline real eps_poly4(
306 real m, /* grid coordinate in certain direction */
307 real K, /* grid size in corresponding direction */
308 real n) /* spline interpolation order of the SPME */
311 real nom = 0; /* nominator */
312 real denom = 0; /* denominator */
320 for (i = -SUMORDER; i < 0; i++)
324 nom += i * i * pow( tmp, -2*n );
327 for (i = SUMORDER; i > 0; i--)
331 nom += i * i * pow( tmp, -2*n );
334 for (i = -SUMORDER; i < SUMORDER+1; i++)
338 denom += pow( tmp, -n );
341 return 4.0 * M_PI * M_PI * nom / denom / denom;
345 static inline real eps_self(
346 real m, /* grid coordinate in certain direction */
347 real K, /* grid size in corresponding direction */
348 rvec rboxv, /* reciprocal box vector */
349 real n, /* spline interpolation order of the SPME */
350 rvec x) /* coordinate of charge */
353 real tmp = 0; /* temporary variables for computations */
354 real tmp1 = 0; /* temporary variables for computations */
355 real tmp2 = 0; /* temporary variables for computations */
356 real rcoord = 0; /* coordinate in certain reciprocal space direction */
357 real nom = 0; /* nominator */
358 real denom = 0; /* denominator */
366 rcoord = iprod(rboxv, x);
369 for (i = -SUMORDER; i < 0; i++)
371 tmp = -sin(2.0 * M_PI * i * K * rcoord);
372 tmp1 = 2.0 * M_PI * m / K + 2.0 * M_PI * i;
373 tmp2 = pow(tmp1, -n);
374 nom += tmp * tmp2 * i;
378 for (i = SUMORDER; i > 0; i--)
380 tmp = -sin(2.0 * M_PI * i * K * rcoord);
381 tmp1 = 2.0 * M_PI * m / K + 2.0 * M_PI * i;
382 tmp2 = pow(tmp1, -n);
383 nom += tmp * tmp2 * i;
388 tmp = 2.0 * M_PI * m / K;
392 return 2.0 * M_PI * nom / denom * K;
398 /* The following routine is just a copy from pme.c */
400 static void calc_recipbox(matrix box, matrix recipbox)
402 /* Save some time by assuming upper right part is zero */
404 real tmp = 1.0/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
406 recipbox[XX][XX] = box[YY][YY]*box[ZZ][ZZ]*tmp;
407 recipbox[XX][YY] = 0;
408 recipbox[XX][ZZ] = 0;
409 recipbox[YY][XX] = -box[YY][XX]*box[ZZ][ZZ]*tmp;
410 recipbox[YY][YY] = box[XX][XX]*box[ZZ][ZZ]*tmp;
411 recipbox[YY][ZZ] = 0;
412 recipbox[ZZ][XX] = (box[YY][XX]*box[ZZ][YY]-box[YY][YY]*box[ZZ][XX])*tmp;
413 recipbox[ZZ][YY] = -box[ZZ][YY]*box[XX][XX]*tmp;
414 recipbox[ZZ][ZZ] = box[XX][XX]*box[YY][YY]*tmp;
418 /* Estimate the reciprocal space part error of the SPME Ewald sum. */
419 static real estimate_reciprocal(
421 rvec x[], /* array of particles */
422 real q[], /* array of charges */
423 int nr, /* number of charges = size of the charge array */
424 FILE gmx_unused *fp_out,
426 unsigned int seed, /* The seed for the random number generator */
427 int *nsamples, /* Return the number of samples used if Monte Carlo
428 * algorithm is used for self energy error estimate */
431 real e_rec = 0; /* reciprocal error estimate */
432 real e_rec1 = 0; /* Error estimate term 1*/
433 real e_rec2 = 0; /* Error estimate term 2*/
434 real e_rec3 = 0; /* Error estimate term 3 */
435 real e_rec3x = 0; /* part of Error estimate term 3 in x */
436 real e_rec3y = 0; /* part of Error estimate term 3 in y */
437 real e_rec3z = 0; /* part of Error estimate term 3 in z */
439 int nx, ny, nz; /* grid coordinates */
440 real q2_all = 0; /* sum of squared charges */
441 rvec gridpx; /* reciprocal grid point in x direction*/
442 rvec gridpxy; /* reciprocal grid point in x and y direction*/
443 rvec gridp; /* complete reciprocal grid point in 3 directions*/
444 rvec tmpvec; /* template to create points from basis vectors */
445 rvec tmpvec2; /* template to create points from basis vectors */
446 real coeff = 0; /* variable to compute coefficients of the error estimate */
447 real coeff2 = 0; /* variable to compute coefficients of the error estimate */
448 real tmp = 0; /* variables to compute different factors from vectors */
453 /* Random number generator */
454 gmx_rng_t rng = NULL;
457 /* Index variables for parallel work distribution */
458 int startglobal, stopglobal;
459 int startlocal, stoplocal;
468 rng = gmx_rng_init(seed);
476 for (i = 0; i < nr; 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>(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;
503 MPI_Barrier(MPI_COMM_WORLD);
519 fprintf(stderr, "Calculating reciprocal error part 1 ...");
523 for (nx = startlocal; nx <= stoplocal; nx++)
525 svmul(nx, info->recipbox[XX], gridpx);
526 for (ny = -info->nky[0]/2; ny < info->nky[0]/2+1; ny++)
528 svmul(ny, info->recipbox[YY], tmpvec);
529 rvec_add(gridpx, tmpvec, gridpxy);
530 for (nz = -info->nkz[0]/2; nz < info->nkz[0]/2+1; nz++)
532 if (0 == nx && 0 == ny && 0 == nz)
536 svmul(nz, info->recipbox[ZZ], tmpvec);
537 rvec_add(gridpxy, tmpvec, gridp);
539 coeff = exp(-1.0 * M_PI * M_PI * tmp / info->ewald_beta[0] / info->ewald_beta[0] );
540 coeff /= 2.0 * M_PI * info->volume * tmp;
544 tmp = eps_poly2(nx, info->nkx[0], info->pme_order[0]);
545 tmp += eps_poly2(ny, info->nkx[0], info->pme_order[0]);
546 tmp += eps_poly2(nz, info->nkx[0], info->pme_order[0]);
548 tmp1 = eps_poly1(nx, info->nkx[0], info->pme_order[0]);
549 tmp2 = eps_poly1(ny, info->nky[0], info->pme_order[0]);
551 tmp += 2.0 * tmp1 * tmp2;
553 tmp1 = eps_poly1(nz, info->nkz[0], info->pme_order[0]);
554 tmp2 = eps_poly1(ny, info->nky[0], info->pme_order[0]);
556 tmp += 2.0 * tmp1 * tmp2;
558 tmp1 = eps_poly1(nz, info->nkz[0], info->pme_order[0]);
559 tmp2 = eps_poly1(nx, info->nkx[0], info->pme_order[0]);
561 tmp += 2.0 * tmp1 * tmp2;
563 tmp1 = eps_poly1(nx, info->nkx[0], info->pme_order[0]);
564 tmp1 += eps_poly1(ny, info->nky[0], info->pme_order[0]);
565 tmp1 += eps_poly1(nz, info->nkz[0], info->pme_order[0]);
569 e_rec1 += 32.0 * M_PI * M_PI * coeff * coeff * coeff2 * tmp * q2_all * q2_all / nr;
571 tmp1 = eps_poly3(nx, info->nkx[0], info->pme_order[0]);
572 tmp1 *= info->nkx[0];
573 tmp2 = iprod(gridp, info->recipbox[XX]);
577 tmp1 = eps_poly3(ny, info->nky[0], info->pme_order[0]);
578 tmp1 *= info->nky[0];
579 tmp2 = iprod(gridp, info->recipbox[YY]);
583 tmp1 = eps_poly3(nz, info->nkz[0], info->pme_order[0]);
584 tmp1 *= info->nkz[0];
585 tmp2 = iprod(gridp, info->recipbox[ZZ]);
591 tmp1 = eps_poly4(nx, info->nkx[0], info->pme_order[0]);
592 tmp1 *= norm2(info->recipbox[XX]);
593 tmp1 *= info->nkx[0] * info->nkx[0];
597 tmp1 = eps_poly4(ny, info->nky[0], info->pme_order[0]);
598 tmp1 *= norm2(info->recipbox[YY]);
599 tmp1 *= info->nky[0] * info->nky[0];
603 tmp1 = eps_poly4(nz, info->nkz[0], info->pme_order[0]);
604 tmp1 *= norm2(info->recipbox[ZZ]);
605 tmp1 *= info->nkz[0] * info->nkz[0];
609 e_rec2 += 4.0 * coeff * coeff * tmp * q2_all * q2_all / nr;
615 fprintf(stderr, "\rCalculating reciprocal error part 1 ... %3.0f%%", 100.0*(nx-startlocal+1)/(x_per_core));
622 fprintf(stderr, "\n");
625 /* Use just a fraction of all charges to estimate the self energy error term? */
626 bFraction = (info->fracself > 0.0) && (info->fracself < 1.0);
630 /* Here xtot is the number of samples taken for the Monte Carlo calculation
631 * of the average of term IV of equation 35 in Wang2010. Round up to a
632 * number of samples that is divisible by the number of nodes */
633 x_per_core = static_cast<int>(ceil(info->fracself * nr / cr->nnodes));
634 xtot = x_per_core * cr->nnodes;
638 /* In this case we use all nr particle positions */
640 x_per_core = static_cast<int>(ceil(static_cast<real>(xtot) / cr->nnodes));
643 startlocal = x_per_core * cr->nodeid;
644 stoplocal = std::min(startlocal + x_per_core, xtot); /* min needed if xtot == nr */
648 /* Make shure we get identical results in serial and parallel. Therefore,
649 * take the sample indices from a single, global random number array that
650 * is constructed on the master node and that only depends on the seed */
654 for (i = 0; i < xtot; i++)
656 numbers[i] = static_cast<int>(floor(gmx_rng_uniform_real(rng) * nr));
659 /* Broadcast the random number array to the other nodes */
662 nblock_bc(cr, xtot, numbers);
665 if (bVerbose && MASTER(cr))
667 fprintf(stdout, "Using %d sample%s to approximate the self interaction error term",
668 xtot, xtot == 1 ? "" : "s");
671 fprintf(stdout, " (%d sample%s per node)", x_per_core, x_per_core == 1 ? "" : "s");
673 fprintf(stdout, ".\n");
677 /* Return the number of positions used for the Monte Carlo algorithm */
680 for (i = startlocal; i < stoplocal; i++)
688 /* Randomly pick a charge */
693 /* Use all charges */
697 /* for(nx=startlocal; nx<=stoplocal; nx++)*/
698 for (nx = -info->nkx[0]/2; nx < info->nkx[0]/2+1; nx++)
700 svmul(nx, info->recipbox[XX], gridpx);
701 for (ny = -info->nky[0]/2; ny < info->nky[0]/2+1; ny++)
703 svmul(ny, info->recipbox[YY], tmpvec);
704 rvec_add(gridpx, tmpvec, gridpxy);
705 for (nz = -info->nkz[0]/2; nz < info->nkz[0]/2+1; nz++)
708 if (0 == nx && 0 == ny && 0 == nz)
713 svmul(nz, info->recipbox[ZZ], tmpvec);
714 rvec_add(gridpxy, tmpvec, gridp);
716 coeff = exp(-1.0 * M_PI * M_PI * tmp / info->ewald_beta[0] / info->ewald_beta[0] );
718 e_rec3x += coeff*eps_self(nx, info->nkx[0], info->recipbox[XX], info->pme_order[0], x[ci]);
719 e_rec3y += coeff*eps_self(ny, info->nky[0], info->recipbox[YY], info->pme_order[0], x[ci]);
720 e_rec3z += coeff*eps_self(nz, info->nkz[0], info->recipbox[ZZ], info->pme_order[0], x[ci]);
728 svmul(e_rec3x, info->recipbox[XX], tmpvec);
729 rvec_inc(tmpvec2, tmpvec);
730 svmul(e_rec3y, info->recipbox[YY], tmpvec);
731 rvec_inc(tmpvec2, tmpvec);
732 svmul(e_rec3z, info->recipbox[ZZ], tmpvec);
733 rvec_inc(tmpvec2, tmpvec);
735 e_rec3 += q[ci]*q[ci]*q[ci]*q[ci]*norm2(tmpvec2) / ( xtot * M_PI * info->volume * M_PI * info->volume);
738 fprintf(stderr, "\rCalculating reciprocal error part 2 ... %3.0f%%",
739 100.0*(i+1)/stoplocal);
746 fprintf(stderr, "\n");
753 t1 = MPI_Wtime() - t0;
754 fprintf(fp_out, "Recip. err. est. took : %lf s\n", t1);
762 fprintf(stderr, "Node %3d: nx=[%3d...%3d] e_rec3=%e\n",
763 cr->nodeid, startlocal, stoplocal, e_rec3);
769 gmx_sum(1, &e_rec1, cr);
770 gmx_sum(1, &e_rec2, cr);
771 gmx_sum(1, &e_rec3, cr);
774 /* e_rec1*=8.0 * q2_all / info->volume / info->volume / nr ;
775 e_rec2*= q2_all / M_PI / M_PI / info->volume / info->volume / nr ;
776 e_rec3/= M_PI * M_PI * info->volume * info->volume * nr ;
778 e_rec = sqrt(e_rec1+e_rec2+e_rec3);
781 return ONE_4PI_EPS0 * e_rec;
785 /* Allocate memory for the inputinfo struct: */
786 static void create_info(t_inputinfo *info)
788 snew(info->fac, info->n_entries);
789 snew(info->rcoulomb, info->n_entries);
790 snew(info->rvdw, info->n_entries);
791 snew(info->nkx, info->n_entries);
792 snew(info->nky, info->n_entries);
793 snew(info->nkz, info->n_entries);
794 snew(info->fourier_sp, info->n_entries);
795 snew(info->ewald_rtol, info->n_entries);
796 snew(info->ewald_beta, info->n_entries);
797 snew(info->pme_order, info->n_entries);
798 snew(info->fn_out, info->n_entries);
799 snew(info->e_dir, info->n_entries);
800 snew(info->e_rec, info->n_entries);
804 /* Allocate and fill an array with coordinates and charges,
805 * returns the number of charges found
807 static int prepare_x_q(real *q[], rvec *x[], gmx_mtop_t *mtop, rvec x_orig[], t_commrec *cr)
810 int nq; /* number of charged particles */
811 gmx_mtop_atomloop_all_t aloop;
817 snew(*q, mtop->natoms);
818 snew(*x, mtop->natoms);
821 aloop = gmx_mtop_atomloop_all_init(mtop);
823 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
825 if (is_charge(atom->q))
828 (*x)[nq][XX] = x_orig[i][XX];
829 (*x)[nq][YY] = x_orig[i][YY];
830 (*x)[nq][ZZ] = x_orig[i][ZZ];
834 /* Give back some unneeded memory */
838 /* Broadcast x and q in the parallel case */
841 /* Transfer the number of charges */
845 nblock_bc(cr, nq, *x);
846 nblock_bc(cr, nq, *q);
854 /* Read in the tpr file and save information we need later in info */
855 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)
857 read_tpx_state(fn_sim_tpr, ir, state, NULL, 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(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, t_commrec *cr)
896 nblock_bc(cr, info->n_entries, info->nkx);
897 nblock_bc(cr, info->n_entries, info->nky);
898 nblock_bc(cr, info->n_entries, info->nkz);
899 nblock_bc(cr, info->n_entries, info->ewald_beta);
900 nblock_bc(cr, info->n_entries, info->pme_order);
901 nblock_bc(cr, info->n_entries, info->e_dir);
902 nblock_bc(cr, info->n_entries, info->e_rec);
903 block_bc(cr, info->volume);
904 block_bc(cr, info->recipbox);
905 block_bc(cr, info->natoms);
906 block_bc(cr, info->fracself);
907 block_bc(cr, info->bTUNE);
908 block_bc(cr, info->q2all);
909 block_bc(cr, 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, t_state *state,
918 gmx_mtop_t *mtop, FILE *fp_out, gmx_bool bVerbose, unsigned int seed,
921 rvec *x = NULL; /* The coordinates */
922 real *q = NULL; /* The charges */
923 real edir = 0.0; /* real space error */
924 real erec = 0.0; /* reciprocal space error */
925 real derr = 0.0; /* difference of real and reciprocal space error */
926 real derr0 = 0.0; /* difference of real and reciprocal space error */
927 real beta = 0.0; /* splitting parameter beta */
928 real beta0 = 0.0; /* splitting parameter beta */
929 int ncharges; /* The number of atoms with charges */
930 int nsamples; /* The number of samples used for the calculation of the
931 * self-energy error term */
936 fprintf(fp_out, "\n--- PME ERROR ESTIMATE ---\n");
939 /* Prepare an x and q array with only the charged atoms */
940 ncharges = prepare_x_q(&q, &x, mtop, state->x, cr);
943 calc_q2all(mtop, &(info->q2all), &(info->q2allnr));
944 info->ewald_rtol[0] = gmx_erfc(info->rcoulomb[0]*info->ewald_beta[0]);
945 /* Write some info to log file */
946 fprintf(fp_out, "Box volume : %g nm^3\n", info->volume);
947 fprintf(fp_out, "Number of charged atoms : %d (total atoms %d)\n", ncharges, info->natoms);
948 fprintf(fp_out, "Coulomb radius : %g nm\n", info->rcoulomb[0]);
949 fprintf(fp_out, "Ewald_rtol : %g\n", info->ewald_rtol[0]);
950 fprintf(fp_out, "Ewald parameter beta : %g\n", info->ewald_beta[0]);
951 fprintf(fp_out, "Interpolation order : %d\n", info->pme_order[0]);
952 fprintf(fp_out, "Fourier grid (nx,ny,nz) : %d x %d x %d\n",
953 info->nkx[0], info->nky[0], info->nkz[0]);
960 bcast_info(info, cr);
964 /* Calculate direct space error */
965 info->e_dir[0] = estimate_direct(info);
967 /* Calculate reciprocal space error */
968 info->e_rec[0] = estimate_reciprocal(info, x, q, ncharges, fp_out, bVerbose,
969 seed, &nsamples, cr);
973 bcast_info(info, cr);
978 fprintf(fp_out, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
979 fprintf(fp_out, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
980 fprintf(fp_out, "Self-energy error term was estimated using %d samples\n", nsamples);
982 fprintf(stderr, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
983 fprintf(stderr, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
992 fprintf(stderr, "Starting tuning ...\n");
994 edir = info->e_dir[0];
995 erec = info->e_rec[0];
997 beta0 = info->ewald_beta[0];
1000 info->ewald_beta[0] += 0.1;
1004 info->ewald_beta[0] -= 0.1;
1006 info->e_dir[0] = estimate_direct(info);
1007 info->e_rec[0] = estimate_reciprocal(info, x, q, ncharges, fp_out, bVerbose,
1008 seed, &nsamples, cr);
1012 bcast_info(info, cr);
1016 edir = info->e_dir[0];
1017 erec = info->e_rec[0];
1019 while (fabs(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);
1029 info->e_rec[0] = estimate_reciprocal(info, x, q, ncharges, fp_out, bVerbose,
1030 seed, &nsamples, cr);
1034 bcast_info(info, cr);
1037 edir = info->e_dir[0];
1038 erec = info->e_rec[0];
1044 fprintf(stderr, "difference between real and rec. space error (step %d): %g\n", i, fabs(derr));
1045 fprintf(stderr, "old beta: %f\n", beta0);
1046 fprintf(stderr, "new beta: %f\n", beta);
1050 info->ewald_rtol[0] = gmx_erfc(info->rcoulomb[0]*info->ewald_beta[0]);
1054 /* Write some info to log file */
1056 fprintf(fp_out, "========= After tuning ========\n");
1057 fprintf(fp_out, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
1058 fprintf(fp_out, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
1059 fprintf(stderr, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
1060 fprintf(stderr, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
1061 fprintf(fp_out, "Ewald_rtol : %g\n", info->ewald_rtol[0]);
1062 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 "[TT]g_pme_error[tt] 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 */
1091 t_inputrec *ir = NULL; /* The inputrec from the tpr file */
1094 unsigned long PCA_Flags;
1095 gmx_bool bTUNE = FALSE;
1096 gmx_bool bVerbose = FALSE;
1100 static t_filenm fnm[] = {
1101 { efTPX, "-s", NULL, ffREAD },
1102 { efOUT, "-o", "error", ffWRITE },
1103 { efTPX, "-so", "tuned", ffOPTWR }
1106 output_env_t oenv = NULL;
1109 { "-beta", FALSE, etREAL, {&user_beta},
1110 "If positive, overwrite ewald_beta from [TT].tpr[tt] file with this value" },
1111 { "-tune", FALSE, etBOOL, {&bTUNE},
1112 "Tune the splitting parameter such that the error is equally distributed between real and reciprocal space" },
1113 { "-self", FALSE, etREAL, {&fracself},
1114 "If between 0.0 and 1.0, determine self interaction error from just this fraction of the charged particles" },
1115 { "-seed", FALSE, etINT, {&seed},
1116 "Random number seed used for Monte Carlo algorithm when [TT]-self[tt] is set to a value between 0.0 and 1.0" },
1117 { "-v", FALSE, etBOOL, {&bVerbose},
1118 "Be loud and noisy" }
1122 #define NFILE asize(fnm)
1126 MPI_Barrier(MPI_COMM_WORLD);
1129 PCA_Flags = PCA_NOEXIT_ON_ARGS;
1130 PCA_Flags |= (MASTER(cr) ? 0 : PCA_QUIET);
1132 if (!parse_common_args(&argc, argv, PCA_Flags,
1133 NFILE, fnm, asize(pa), pa, asize(desc), desc,
1141 bTUNE = opt2bSet("-so", NFILE, fnm);
1146 /* Allocate memory for the inputinfo struct: */
1148 info.fourier_sp[0] = fs;
1150 /* Read in the tpr file and open logfile for reading */
1154 read_tpr_file(opt2fn("-s", NFILE, fnm), &info, &state, &mtop, ir, user_beta, fracself);
1156 fp = fopen(opt2fn("-o", NFILE, fnm), "w");
1159 /* Check consistency if the user provided fourierspacing */
1160 if (fs > 0 && MASTER(cr))
1162 /* Recalculate the grid dimensions using fourierspacing from user input */
1166 calc_grid(stdout, state.box, info.fourier_sp[0], &(info.nkx[0]), &(info.nky[0]), &(info.nkz[0]));
1167 if ( (ir->nkx != info.nkx[0]) || (ir->nky != info.nky[0]) || (ir->nkz != info.nkz[0]) )
1169 gmx_fatal(FARGS, "Wrong fourierspacing %f nm, input file grid = %d x %d x %d, computed grid = %d x %d x %d",
1170 fs, ir->nkx, ir->nky, ir->nkz, info.nkx[0], info.nky[0], info.nkz[0]);
1174 /* Estimate (S)PME force error */
1176 /* Determine the volume of the simulation box */
1179 info.volume = det(state.box);
1180 calc_recipbox(state.box, info.recipbox);
1181 info.natoms = mtop.natoms;
1187 bcast_info(&info, cr);
1190 /* Get an error estimate of the input tpr file and do some tuning if requested */
1191 estimate_PME_error(&info, &state, &mtop, fp, bVerbose, seed, cr);
1195 /* Write out optimized tpr file if requested */
1196 if (opt2bSet("-so", NFILE, fnm) || bTUNE)
1198 ir->ewald_rtol = info.ewald_rtol[0];
1199 write_tpx_state(opt2fn("-so", NFILE, fnm), ir, &state, &mtop);
1201 please_cite(fp, "Wang2010");