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.
41 #include "gromacs/commandline/pargs.h"
42 #include "gromacs/fileio/tpxio.h"
43 #include "gromacs/gmxana/gmx_ana.h"
44 #include "gromacs/legacyheaders/calcgrid.h"
45 #include "gromacs/legacyheaders/checkpoint.h"
46 #include "gromacs/legacyheaders/copyrite.h"
47 #include "gromacs/legacyheaders/coulomb.h"
48 #include "gromacs/legacyheaders/macros.h"
49 #include "gromacs/legacyheaders/main.h"
50 #include "gromacs/legacyheaders/mdatoms.h"
51 #include "gromacs/legacyheaders/network.h"
52 #include "gromacs/legacyheaders/readinp.h"
53 #include "gromacs/legacyheaders/typedefs.h"
54 #include "gromacs/legacyheaders/types/commrec.h"
55 #include "gromacs/math/units.h"
56 #include "gromacs/math/vec.h"
57 #include "gromacs/random/random.h"
58 #include "gromacs/topology/mtop_util.h"
59 #include "gromacs/utility/fatalerror.h"
60 #include "gromacs/utility/smalloc.h"
62 /* We use the same defines as in mvdata.c here */
63 #define block_bc(cr, d) gmx_bcast( sizeof(d), &(d), (cr))
64 #define nblock_bc(cr, nr, d) gmx_bcast((nr)*sizeof((d)[0]), (d), (cr))
65 #define snew_bc(cr, d, nr) { if (!MASTER(cr)) {snew((d), (nr)); }}
66 /* #define TAKETIME */
69 ddnoSEL, ddnoINTERLEAVE, ddnoPP_PME, ddnoCARTESIAN, ddnoNR
72 /* Enum for situations that can occur during log file parsing */
78 eParselogResetProblem,
85 gmx_int64_t orig_sim_steps; /* Number of steps to be done in the real simulation */
86 int n_entries; /* Number of entries in arrays */
87 real volume; /* The volume of the box */
88 matrix recipbox; /* The reciprocal box */
89 int natoms; /* The number of atoms in the MD system */
90 real *fac; /* The scaling factor */
91 real *rcoulomb; /* The coulomb radii [0...nr_inputfiles] */
92 real *rvdw; /* The vdW radii */
93 int *nkx, *nky, *nkz; /* Number of k vectors in each spatial dimension */
94 real *fourier_sp; /* Fourierspacing */
95 real *ewald_rtol; /* Real space tolerance for Ewald, determines */
96 /* the real/reciprocal space relative weight */
97 real *ewald_beta; /* Splitting parameter [1/nm] */
98 real fracself; /* fraction of particles for SI error */
99 real q2all; /* sum ( q ^2 ) */
100 real q2allnr; /* nr of charges */
101 int *pme_order; /* Interpolation order for PME (bsplines) */
102 char **fn_out; /* Name of the output tpr file */
103 real *e_dir; /* Direct space part of PME error with these settings */
104 real *e_rec; /* Reciprocal space part of PME error */
105 gmx_bool bTUNE; /* flag for tuning */
109 /* Returns TRUE when atom is charged */
110 static gmx_bool is_charge(real charge)
112 if (charge*charge > GMX_REAL_EPS)
123 /* calculate charge density */
124 static void calc_q2all(
125 gmx_mtop_t *mtop, /* molecular topology */
126 real *q2all, real *q2allnr)
128 int imol, iatom; /* indices for loops */
129 real q2_all = 0; /* Sum of squared charges */
130 int nrq_mol; /* Number of charges in a single molecule */
131 int nrq_all; /* Total number of charges in the MD system */
133 gmx_moltype_t *molecule;
134 gmx_molblock_t *molblock;
137 fprintf(stderr, "\nCharge density:\n");
139 q2_all = 0.0; /* total q squared */
140 nrq_all = 0; /* total number of charges in the system */
141 for (imol = 0; imol < mtop->nmolblock; imol++) /* Loop over molecule types */
143 q2_mol = 0.0; /* q squared value of this molecule */
144 nrq_mol = 0; /* number of charges this molecule carries */
145 molecule = &(mtop->moltype[imol]);
146 molblock = &(mtop->molblock[imol]);
147 for (iatom = 0; iatom < molblock->natoms_mol; iatom++) /* Loop over atoms in this molecule */
149 qi = molecule->atoms.atom[iatom].q; /* Charge of this atom */
150 /* Is this charge worth to be considered? */
157 /* Multiply with the number of molecules present of this type and add */
158 q2_all += q2_mol*molblock->nmol;
159 nrq_all += nrq_mol*molblock->nmol;
161 fprintf(stderr, "Molecule %2d (%5d atoms) q2_mol=%10.3e nr.mol.charges=%5d (%6dx) q2_all=%10.3e tot.charges=%d\n",
162 imol, molblock->natoms_mol, q2_mol, nrq_mol, molblock->nmol, q2_all, nrq_all);
172 /* Estimate the direct space part error of the SPME Ewald sum */
173 static real estimate_direct(
177 real e_dir = 0; /* Error estimate */
178 real beta = 0; /* Splitting parameter (1/nm) */
179 real r_coulomb = 0; /* Cut-off in direct space */
182 beta = info->ewald_beta[0];
183 r_coulomb = info->rcoulomb[0];
185 e_dir = 2.0 * info->q2all * gmx_invsqrt( info->q2allnr * r_coulomb * info->volume );
186 e_dir *= exp (-beta*beta*r_coulomb*r_coulomb);
188 return ONE_4PI_EPS0*e_dir;
193 /* the following 4 functions determine polynomials required for the reciprocal error estimate */
195 static inline real eps_poly1(
196 real m, /* grid coordinate in certain direction */
197 real K, /* grid size in corresponding direction */
198 real n) /* spline interpolation order of the SPME */
201 real nom = 0; /* nominator */
202 real denom = 0; /* denominator */
210 for (i = -SUMORDER; i < 0; i++)
214 nom += pow( tmp, -n );
217 for (i = SUMORDER; i > 0; i--)
221 nom += pow( tmp, -n );
226 denom = pow( tmp, -n )+nom;
232 static inline real eps_poly2(
233 real m, /* grid coordinate in certain direction */
234 real K, /* grid size in corresponding direction */
235 real n) /* spline interpolation order of the SPME */
238 real nom = 0; /* nominator */
239 real denom = 0; /* denominator */
247 for (i = -SUMORDER; i < 0; i++)
251 nom += pow( tmp, -2*n );
254 for (i = SUMORDER; i > 0; i--)
258 nom += pow( tmp, -2*n );
261 for (i = -SUMORDER; i < SUMORDER+1; i++)
265 denom += pow( tmp, -n );
267 tmp = eps_poly1(m, K, n);
268 return nom / denom / denom + tmp*tmp;
272 static inline real eps_poly3(
273 real m, /* grid coordinate in certain direction */
274 real K, /* grid size in corresponding direction */
275 real n) /* spline interpolation order of the SPME */
278 real nom = 0; /* nominator */
279 real denom = 0; /* denominator */
287 for (i = -SUMORDER; i < 0; i++)
291 nom += i * pow( tmp, -2*n );
294 for (i = SUMORDER; i > 0; i--)
298 nom += i * pow( tmp, -2*n );
301 for (i = -SUMORDER; i < SUMORDER+1; i++)
305 denom += pow( tmp, -n );
308 return 2.0 * M_PI * nom / denom / denom;
312 static inline real eps_poly4(
313 real m, /* grid coordinate in certain direction */
314 real K, /* grid size in corresponding direction */
315 real n) /* spline interpolation order of the SPME */
318 real nom = 0; /* nominator */
319 real denom = 0; /* denominator */
327 for (i = -SUMORDER; i < 0; i++)
331 nom += i * i * pow( tmp, -2*n );
334 for (i = SUMORDER; i > 0; i--)
338 nom += i * i * pow( tmp, -2*n );
341 for (i = -SUMORDER; i < SUMORDER+1; i++)
345 denom += pow( tmp, -n );
348 return 4.0 * M_PI * M_PI * nom / denom / denom;
352 static inline real eps_self(
353 real m, /* grid coordinate in certain direction */
354 real K, /* grid size in corresponding direction */
355 rvec rboxv, /* reciprocal box vector */
356 real n, /* spline interpolation order of the SPME */
357 rvec x) /* coordinate of charge */
360 real tmp = 0; /* temporary variables for computations */
361 real tmp1 = 0; /* temporary variables for computations */
362 real tmp2 = 0; /* temporary variables for computations */
363 real rcoord = 0; /* coordinate in certain reciprocal space direction */
364 real nom = 0; /* nominator */
365 real denom = 0; /* denominator */
373 rcoord = iprod(rboxv, x);
376 for (i = -SUMORDER; i < 0; i++)
378 tmp = -sin(2.0 * M_PI * i * K * rcoord);
379 tmp1 = 2.0 * M_PI * m / K + 2.0 * M_PI * i;
380 tmp2 = pow(tmp1, -n);
381 nom += tmp * tmp2 * i;
385 for (i = SUMORDER; i > 0; i--)
387 tmp = -sin(2.0 * M_PI * i * K * rcoord);
388 tmp1 = 2.0 * M_PI * m / K + 2.0 * M_PI * i;
389 tmp2 = pow(tmp1, -n);
390 nom += tmp * tmp2 * i;
395 tmp = 2.0 * M_PI * m / K;
399 return 2.0 * M_PI * nom / denom * K;
405 /* The following routine is just a copy from pme.c */
407 static void calc_recipbox(matrix box, matrix recipbox)
409 /* Save some time by assuming upper right part is zero */
411 real tmp = 1.0/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
413 recipbox[XX][XX] = box[YY][YY]*box[ZZ][ZZ]*tmp;
414 recipbox[XX][YY] = 0;
415 recipbox[XX][ZZ] = 0;
416 recipbox[YY][XX] = -box[YY][XX]*box[ZZ][ZZ]*tmp;
417 recipbox[YY][YY] = box[XX][XX]*box[ZZ][ZZ]*tmp;
418 recipbox[YY][ZZ] = 0;
419 recipbox[ZZ][XX] = (box[YY][XX]*box[ZZ][YY]-box[YY][YY]*box[ZZ][XX])*tmp;
420 recipbox[ZZ][YY] = -box[ZZ][YY]*box[XX][XX]*tmp;
421 recipbox[ZZ][ZZ] = box[XX][XX]*box[YY][YY]*tmp;
425 /* Estimate the reciprocal space part error of the SPME Ewald sum. */
426 static real estimate_reciprocal(
428 rvec x[], /* array of particles */
429 real q[], /* array of charges */
430 int nr, /* number of charges = size of the charge array */
431 FILE gmx_unused *fp_out,
433 unsigned int seed, /* The seed for the random number generator */
434 int *nsamples, /* Return the number of samples used if Monte Carlo
435 * algorithm is used for self energy error estimate */
438 real e_rec = 0; /* reciprocal error estimate */
439 real e_rec1 = 0; /* Error estimate term 1*/
440 real e_rec2 = 0; /* Error estimate term 2*/
441 real e_rec3 = 0; /* Error estimate term 3 */
442 real e_rec3x = 0; /* part of Error estimate term 3 in x */
443 real e_rec3y = 0; /* part of Error estimate term 3 in y */
444 real e_rec3z = 0; /* part of Error estimate term 3 in z */
446 int nx, ny, nz; /* grid coordinates */
447 real q2_all = 0; /* sum of squared charges */
448 rvec gridpx; /* reciprocal grid point in x direction*/
449 rvec gridpxy; /* reciprocal grid point in x and y direction*/
450 rvec gridp; /* complete reciprocal grid point in 3 directions*/
451 rvec tmpvec; /* template to create points from basis vectors */
452 rvec tmpvec2; /* template to create points from basis vectors */
453 real coeff = 0; /* variable to compute coefficients of the error estimate */
454 real coeff2 = 0; /* variable to compute coefficients of the error estimate */
455 real tmp = 0; /* variables to compute different factors from vectors */
460 /* Random number generator */
461 gmx_rng_t rng = NULL;
464 /* Index variables for parallel work distribution */
465 int startglobal, stopglobal;
466 int startlocal, stoplocal;
475 rng = gmx_rng_init(seed);
483 for (i = 0; i < nr; i++)
488 /* Calculate indices for work distribution */
489 startglobal = -info->nkx[0]/2;
490 stopglobal = info->nkx[0]/2;
491 xtot = stopglobal*2+1;
494 x_per_core = static_cast<int>(ceil(static_cast<real>(xtot) / cr->nnodes));
495 startlocal = startglobal + x_per_core*cr->nodeid;
496 stoplocal = startlocal + x_per_core -1;
497 if (stoplocal > stopglobal)
499 stoplocal = stopglobal;
504 startlocal = startglobal;
505 stoplocal = stopglobal;
510 MPI_Barrier(MPI_COMM_WORLD);
526 fprintf(stderr, "Calculating reciprocal error part 1 ...");
530 for (nx = startlocal; nx <= stoplocal; nx++)
532 svmul(nx, info->recipbox[XX], gridpx);
533 for (ny = -info->nky[0]/2; ny < info->nky[0]/2+1; ny++)
535 svmul(ny, info->recipbox[YY], tmpvec);
536 rvec_add(gridpx, tmpvec, gridpxy);
537 for (nz = -info->nkz[0]/2; nz < info->nkz[0]/2+1; nz++)
539 if (0 == nx && 0 == ny && 0 == nz)
543 svmul(nz, info->recipbox[ZZ], tmpvec);
544 rvec_add(gridpxy, tmpvec, gridp);
546 coeff = exp(-1.0 * M_PI * M_PI * tmp / info->ewald_beta[0] / info->ewald_beta[0] );
547 coeff /= 2.0 * M_PI * info->volume * tmp;
551 tmp = eps_poly2(nx, info->nkx[0], info->pme_order[0]);
552 tmp += eps_poly2(ny, info->nkx[0], info->pme_order[0]);
553 tmp += eps_poly2(nz, info->nkx[0], info->pme_order[0]);
555 tmp1 = eps_poly1(nx, info->nkx[0], info->pme_order[0]);
556 tmp2 = eps_poly1(ny, info->nky[0], info->pme_order[0]);
558 tmp += 2.0 * tmp1 * tmp2;
560 tmp1 = eps_poly1(nz, info->nkz[0], info->pme_order[0]);
561 tmp2 = eps_poly1(ny, info->nky[0], info->pme_order[0]);
563 tmp += 2.0 * tmp1 * tmp2;
565 tmp1 = eps_poly1(nz, info->nkz[0], info->pme_order[0]);
566 tmp2 = eps_poly1(nx, info->nkx[0], info->pme_order[0]);
568 tmp += 2.0 * tmp1 * tmp2;
570 tmp1 = eps_poly1(nx, info->nkx[0], info->pme_order[0]);
571 tmp1 += eps_poly1(ny, info->nky[0], info->pme_order[0]);
572 tmp1 += eps_poly1(nz, info->nkz[0], info->pme_order[0]);
576 e_rec1 += 32.0 * M_PI * M_PI * coeff * coeff * coeff2 * tmp * q2_all * q2_all / nr;
578 tmp1 = eps_poly3(nx, info->nkx[0], info->pme_order[0]);
579 tmp1 *= info->nkx[0];
580 tmp2 = iprod(gridp, info->recipbox[XX]);
584 tmp1 = eps_poly3(ny, info->nky[0], info->pme_order[0]);
585 tmp1 *= info->nky[0];
586 tmp2 = iprod(gridp, info->recipbox[YY]);
590 tmp1 = eps_poly3(nz, info->nkz[0], info->pme_order[0]);
591 tmp1 *= info->nkz[0];
592 tmp2 = iprod(gridp, info->recipbox[ZZ]);
598 tmp1 = eps_poly4(nx, info->nkx[0], info->pme_order[0]);
599 tmp1 *= norm2(info->recipbox[XX]);
600 tmp1 *= info->nkx[0] * info->nkx[0];
604 tmp1 = eps_poly4(ny, info->nky[0], info->pme_order[0]);
605 tmp1 *= norm2(info->recipbox[YY]);
606 tmp1 *= info->nky[0] * info->nky[0];
610 tmp1 = eps_poly4(nz, info->nkz[0], info->pme_order[0]);
611 tmp1 *= norm2(info->recipbox[ZZ]);
612 tmp1 *= info->nkz[0] * info->nkz[0];
616 e_rec2 += 4.0 * coeff * coeff * tmp * q2_all * q2_all / nr;
622 fprintf(stderr, "\rCalculating reciprocal error part 1 ... %3.0f%%", 100.0*(nx-startlocal+1)/(x_per_core));
629 fprintf(stderr, "\n");
632 /* Use just a fraction of all charges to estimate the self energy error term? */
633 bFraction = (info->fracself > 0.0) && (info->fracself < 1.0);
637 /* Here xtot is the number of samples taken for the Monte Carlo calculation
638 * of the average of term IV of equation 35 in Wang2010. Round up to a
639 * number of samples that is divisible by the number of nodes */
640 x_per_core = static_cast<int>(ceil(info->fracself * nr / cr->nnodes));
641 xtot = x_per_core * cr->nnodes;
645 /* In this case we use all nr particle positions */
647 x_per_core = static_cast<int>(ceil(static_cast<real>(xtot) / cr->nnodes));
650 startlocal = x_per_core * cr->nodeid;
651 stoplocal = std::min(startlocal + x_per_core, xtot); /* min needed if xtot == nr */
655 /* Make shure we get identical results in serial and parallel. Therefore,
656 * take the sample indices from a single, global random number array that
657 * is constructed on the master node and that only depends on the seed */
661 for (i = 0; i < xtot; i++)
663 numbers[i] = static_cast<int>(floor(gmx_rng_uniform_real(rng) * nr));
666 /* Broadcast the random number array to the other nodes */
669 nblock_bc(cr, xtot, numbers);
672 if (bVerbose && MASTER(cr))
674 fprintf(stdout, "Using %d sample%s to approximate the self interaction error term",
675 xtot, xtot == 1 ? "" : "s");
678 fprintf(stdout, " (%d sample%s per rank)", x_per_core, x_per_core == 1 ? "" : "s");
680 fprintf(stdout, ".\n");
684 /* Return the number of positions used for the Monte Carlo algorithm */
687 for (i = startlocal; i < stoplocal; i++)
695 /* Randomly pick a charge */
700 /* Use all charges */
704 /* for(nx=startlocal; nx<=stoplocal; nx++)*/
705 for (nx = -info->nkx[0]/2; nx < info->nkx[0]/2+1; nx++)
707 svmul(nx, info->recipbox[XX], gridpx);
708 for (ny = -info->nky[0]/2; ny < info->nky[0]/2+1; ny++)
710 svmul(ny, info->recipbox[YY], tmpvec);
711 rvec_add(gridpx, tmpvec, gridpxy);
712 for (nz = -info->nkz[0]/2; nz < info->nkz[0]/2+1; nz++)
715 if (0 == nx && 0 == ny && 0 == nz)
720 svmul(nz, info->recipbox[ZZ], tmpvec);
721 rvec_add(gridpxy, tmpvec, gridp);
723 coeff = exp(-1.0 * M_PI * M_PI * tmp / info->ewald_beta[0] / info->ewald_beta[0] );
725 e_rec3x += coeff*eps_self(nx, info->nkx[0], info->recipbox[XX], info->pme_order[0], x[ci]);
726 e_rec3y += coeff*eps_self(ny, info->nky[0], info->recipbox[YY], info->pme_order[0], x[ci]);
727 e_rec3z += coeff*eps_self(nz, info->nkz[0], info->recipbox[ZZ], info->pme_order[0], x[ci]);
735 svmul(e_rec3x, info->recipbox[XX], tmpvec);
736 rvec_inc(tmpvec2, tmpvec);
737 svmul(e_rec3y, info->recipbox[YY], tmpvec);
738 rvec_inc(tmpvec2, tmpvec);
739 svmul(e_rec3z, info->recipbox[ZZ], tmpvec);
740 rvec_inc(tmpvec2, tmpvec);
742 e_rec3 += q[ci]*q[ci]*q[ci]*q[ci]*norm2(tmpvec2) / ( xtot * M_PI * info->volume * M_PI * info->volume);
745 fprintf(stderr, "\rCalculating reciprocal error part 2 ... %3.0f%%",
746 100.0*(i+1)/stoplocal);
753 fprintf(stderr, "\n");
760 t1 = MPI_Wtime() - t0;
761 fprintf(fp_out, "Recip. err. est. took : %lf s\n", t1);
769 fprintf(stderr, "Rank %3d: nx=[%3d...%3d] e_rec3=%e\n",
770 cr->nodeid, startlocal, stoplocal, e_rec3);
776 gmx_sum(1, &e_rec1, cr);
777 gmx_sum(1, &e_rec2, cr);
778 gmx_sum(1, &e_rec3, cr);
781 /* e_rec1*=8.0 * q2_all / info->volume / info->volume / nr ;
782 e_rec2*= q2_all / M_PI / M_PI / info->volume / info->volume / nr ;
783 e_rec3/= M_PI * M_PI * info->volume * info->volume * nr ;
785 e_rec = sqrt(e_rec1+e_rec2+e_rec3);
788 return ONE_4PI_EPS0 * e_rec;
792 /* Allocate memory for the inputinfo struct: */
793 static void create_info(t_inputinfo *info)
795 snew(info->fac, info->n_entries);
796 snew(info->rcoulomb, info->n_entries);
797 snew(info->rvdw, info->n_entries);
798 snew(info->nkx, info->n_entries);
799 snew(info->nky, info->n_entries);
800 snew(info->nkz, info->n_entries);
801 snew(info->fourier_sp, info->n_entries);
802 snew(info->ewald_rtol, info->n_entries);
803 snew(info->ewald_beta, info->n_entries);
804 snew(info->pme_order, info->n_entries);
805 snew(info->fn_out, info->n_entries);
806 snew(info->e_dir, info->n_entries);
807 snew(info->e_rec, info->n_entries);
811 /* Allocate and fill an array with coordinates and charges,
812 * returns the number of charges found
814 static int prepare_x_q(real *q[], rvec *x[], gmx_mtop_t *mtop, rvec x_orig[], t_commrec *cr)
817 int nq; /* number of charged particles */
818 gmx_mtop_atomloop_all_t aloop;
824 snew(*q, mtop->natoms);
825 snew(*x, mtop->natoms);
828 aloop = gmx_mtop_atomloop_all_init(mtop);
830 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
832 if (is_charge(atom->q))
835 (*x)[nq][XX] = x_orig[i][XX];
836 (*x)[nq][YY] = x_orig[i][YY];
837 (*x)[nq][ZZ] = x_orig[i][ZZ];
841 /* Give back some unneeded memory */
845 /* Broadcast x and q in the parallel case */
848 /* Transfer the number of charges */
852 nblock_bc(cr, nq, *x);
853 nblock_bc(cr, nq, *q);
861 /* Read in the tpr file and save information we need later in info */
862 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)
864 read_tpx_state(fn_sim_tpr, ir, state, NULL, mtop);
866 /* The values of the original tpr input file are save in the first
867 * place [0] of the arrays */
868 info->orig_sim_steps = ir->nsteps;
869 info->pme_order[0] = ir->pme_order;
870 info->rcoulomb[0] = ir->rcoulomb;
871 info->rvdw[0] = ir->rvdw;
872 info->nkx[0] = ir->nkx;
873 info->nky[0] = ir->nky;
874 info->nkz[0] = ir->nkz;
875 info->ewald_rtol[0] = ir->ewald_rtol;
876 info->fracself = fracself;
879 info->ewald_beta[0] = user_beta;
883 info->ewald_beta[0] = calc_ewaldcoeff_q(info->rcoulomb[0], info->ewald_rtol[0]);
886 /* Check if PME was chosen */
887 if (EEL_PME(ir->coulombtype) == FALSE)
889 gmx_fatal(FARGS, "Can only do optimizations for simulations with PME");
892 /* Check if rcoulomb == rlist, which is necessary for PME */
893 if (!(ir->rcoulomb == ir->rlist))
895 gmx_fatal(FARGS, "PME requires rcoulomb (%f) to be equal to rlist (%f).", ir->rcoulomb, ir->rlist);
900 /* Transfer what we need for parallelizing the reciprocal error estimate */
901 static void bcast_info(t_inputinfo *info, t_commrec *cr)
903 nblock_bc(cr, info->n_entries, info->nkx);
904 nblock_bc(cr, info->n_entries, info->nky);
905 nblock_bc(cr, info->n_entries, info->nkz);
906 nblock_bc(cr, info->n_entries, info->ewald_beta);
907 nblock_bc(cr, info->n_entries, info->pme_order);
908 nblock_bc(cr, info->n_entries, info->e_dir);
909 nblock_bc(cr, info->n_entries, info->e_rec);
910 block_bc(cr, info->volume);
911 block_bc(cr, info->recipbox);
912 block_bc(cr, info->natoms);
913 block_bc(cr, info->fracself);
914 block_bc(cr, info->bTUNE);
915 block_bc(cr, info->q2all);
916 block_bc(cr, info->q2allnr);
920 /* Estimate the error of the SPME Ewald sum. This estimate is based upon
921 * a) a homogeneous distribution of the charges
922 * b) a total charge of zero.
924 static void estimate_PME_error(t_inputinfo *info, t_state *state,
925 gmx_mtop_t *mtop, FILE *fp_out, gmx_bool bVerbose, unsigned int seed,
928 rvec *x = NULL; /* The coordinates */
929 real *q = NULL; /* The charges */
930 real edir = 0.0; /* real space error */
931 real erec = 0.0; /* reciprocal space error */
932 real derr = 0.0; /* difference of real and reciprocal space error */
933 real derr0 = 0.0; /* difference of real and reciprocal space error */
934 real beta = 0.0; /* splitting parameter beta */
935 real beta0 = 0.0; /* splitting parameter beta */
936 int ncharges; /* The number of atoms with charges */
937 int nsamples; /* The number of samples used for the calculation of the
938 * self-energy error term */
943 fprintf(fp_out, "\n--- PME ERROR ESTIMATE ---\n");
946 /* Prepare an x and q array with only the charged atoms */
947 ncharges = prepare_x_q(&q, &x, mtop, state->x, cr);
950 calc_q2all(mtop, &(info->q2all), &(info->q2allnr));
951 info->ewald_rtol[0] = gmx_erfc(info->rcoulomb[0]*info->ewald_beta[0]);
952 /* Write some info to log file */
953 fprintf(fp_out, "Box volume : %g nm^3\n", info->volume);
954 fprintf(fp_out, "Number of charged atoms : %d (total atoms %d)\n", ncharges, info->natoms);
955 fprintf(fp_out, "Coulomb radius : %g nm\n", info->rcoulomb[0]);
956 fprintf(fp_out, "Ewald_rtol : %g\n", info->ewald_rtol[0]);
957 fprintf(fp_out, "Ewald parameter beta : %g\n", info->ewald_beta[0]);
958 fprintf(fp_out, "Interpolation order : %d\n", info->pme_order[0]);
959 fprintf(fp_out, "Fourier grid (nx,ny,nz) : %d x %d x %d\n",
960 info->nkx[0], info->nky[0], info->nkz[0]);
967 bcast_info(info, cr);
971 /* Calculate direct space error */
972 info->e_dir[0] = estimate_direct(info);
974 /* Calculate reciprocal space error */
975 info->e_rec[0] = estimate_reciprocal(info, x, q, ncharges, fp_out, bVerbose,
976 seed, &nsamples, cr);
980 bcast_info(info, cr);
985 fprintf(fp_out, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
986 fprintf(fp_out, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
987 fprintf(fp_out, "Self-energy error term was estimated using %d samples\n", nsamples);
989 fprintf(stderr, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
990 fprintf(stderr, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
999 fprintf(stderr, "Starting tuning ...\n");
1001 edir = info->e_dir[0];
1002 erec = info->e_rec[0];
1004 beta0 = info->ewald_beta[0];
1007 info->ewald_beta[0] += 0.1;
1011 info->ewald_beta[0] -= 0.1;
1013 info->e_dir[0] = estimate_direct(info);
1014 info->e_rec[0] = estimate_reciprocal(info, x, q, ncharges, fp_out, bVerbose,
1015 seed, &nsamples, cr);
1019 bcast_info(info, cr);
1023 edir = info->e_dir[0];
1024 erec = info->e_rec[0];
1026 while (fabs(derr/std::min(erec, edir)) > 1e-4)
1029 beta = info->ewald_beta[0];
1030 beta -= derr*(info->ewald_beta[0]-beta0)/(derr-derr0);
1031 beta0 = info->ewald_beta[0];
1032 info->ewald_beta[0] = beta;
1035 info->e_dir[0] = estimate_direct(info);
1036 info->e_rec[0] = estimate_reciprocal(info, x, q, ncharges, fp_out, bVerbose,
1037 seed, &nsamples, cr);
1041 bcast_info(info, cr);
1044 edir = info->e_dir[0];
1045 erec = info->e_rec[0];
1051 fprintf(stderr, "difference between real and rec. space error (step %d): %g\n", i, fabs(derr));
1052 fprintf(stderr, "old beta: %f\n", beta0);
1053 fprintf(stderr, "new beta: %f\n", beta);
1057 info->ewald_rtol[0] = gmx_erfc(info->rcoulomb[0]*info->ewald_beta[0]);
1061 /* Write some info to log file */
1063 fprintf(fp_out, "========= After tuning ========\n");
1064 fprintf(fp_out, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
1065 fprintf(fp_out, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
1066 fprintf(stderr, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
1067 fprintf(stderr, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
1068 fprintf(fp_out, "Ewald_rtol : %g\n", info->ewald_rtol[0]);
1069 fprintf(fp_out, "Ewald parameter beta : %g\n", info->ewald_beta[0]);
1079 int gmx_pme_error(int argc, char *argv[])
1081 const char *desc[] = {
1082 "[THISMODULE] estimates the error of the electrostatic forces",
1083 "if using the sPME algorithm. The flag [TT]-tune[tt] will determine",
1084 "the splitting parameter such that the error is equally",
1085 "distributed over the real and reciprocal space part.",
1086 "The part of the error that stems from self interaction of the particles "
1087 "is computationally demanding. However, a good a approximation is to",
1088 "just use a fraction of the particles for this term which can be",
1089 "indicated by the flag [TT]-self[tt].[PAR]",
1092 real fs = 0.0; /* 0 indicates: not set by the user */
1093 real user_beta = -1.0;
1094 real fracself = 1.0;
1096 t_state state; /* The state from the tpr input file */
1097 gmx_mtop_t mtop; /* The topology from the tpr input file */
1098 t_inputrec *ir = NULL; /* The inputrec from the tpr file */
1101 unsigned long PCA_Flags;
1102 gmx_bool bTUNE = FALSE;
1103 gmx_bool bVerbose = FALSE;
1107 static t_filenm fnm[] = {
1108 { efTPR, "-s", NULL, ffREAD },
1109 { efOUT, "-o", "error", ffWRITE },
1110 { efTPR, "-so", "tuned", ffOPTWR }
1113 output_env_t oenv = NULL;
1116 { "-beta", FALSE, etREAL, {&user_beta},
1117 "If positive, overwrite ewald_beta from [TT].tpr[tt] file with this value" },
1118 { "-tune", FALSE, etBOOL, {&bTUNE},
1119 "Tune the splitting parameter such that the error is equally distributed between real and reciprocal space" },
1120 { "-self", FALSE, etREAL, {&fracself},
1121 "If between 0.0 and 1.0, determine self interaction error from just this fraction of the charged particles" },
1122 { "-seed", FALSE, etINT, {&seed},
1123 "Random number seed used for Monte Carlo algorithm when [TT]-self[tt] is set to a value between 0.0 and 1.0" },
1124 { "-v", FALSE, etBOOL, {&bVerbose},
1125 "Be loud and noisy" }
1129 #define NFILE asize(fnm)
1131 cr = init_commrec();
1133 PCA_Flags = PCA_NOEXIT_ON_ARGS;
1135 if (!parse_common_args(&argc, argv, PCA_Flags,
1136 NFILE, fnm, asize(pa), pa, asize(desc), desc,
1144 bTUNE = opt2bSet("-so", NFILE, fnm);
1149 /* Allocate memory for the inputinfo struct: */
1151 info.fourier_sp[0] = fs;
1153 /* Read in the tpr file and open logfile for reading */
1157 read_tpr_file(opt2fn("-s", NFILE, fnm), &info, &state, &mtop, ir, user_beta, fracself);
1159 fp = fopen(opt2fn("-o", NFILE, fnm), "w");
1162 /* Check consistency if the user provided fourierspacing */
1163 if (fs > 0 && MASTER(cr))
1165 /* Recalculate the grid dimensions using fourierspacing from user input */
1169 calc_grid(stdout, state.box, info.fourier_sp[0], &(info.nkx[0]), &(info.nky[0]), &(info.nkz[0]));
1170 if ( (ir->nkx != info.nkx[0]) || (ir->nky != info.nky[0]) || (ir->nkz != info.nkz[0]) )
1172 gmx_fatal(FARGS, "Wrong fourierspacing %f nm, input file grid = %d x %d x %d, computed grid = %d x %d x %d",
1173 fs, ir->nkx, ir->nky, ir->nkz, info.nkx[0], info.nky[0], info.nkz[0]);
1177 /* Estimate (S)PME force error */
1179 /* Determine the volume of the simulation box */
1182 info.volume = det(state.box);
1183 calc_recipbox(state.box, info.recipbox);
1184 info.natoms = mtop.natoms;
1190 bcast_info(&info, cr);
1193 /* Get an error estimate of the input tpr file and do some tuning if requested */
1194 estimate_PME_error(&info, &state, &mtop, fp, bVerbose, seed, cr);
1198 /* Write out optimized tpr file if requested */
1199 if (opt2bSet("-so", NFILE, fnm) || bTUNE)
1201 ir->ewald_rtol = info.ewald_rtol[0];
1202 write_tpx_state(opt2fn("-so", NFILE, fnm), ir, &state, &mtop);
1204 please_cite(fp, "Wang2010");