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
42 #include "checkpoint.h"
44 #include "gmx_random.h"
48 #include "mtop_util.h"
52 /* We use the same defines as in mvdata.c here */
53 #define block_bc(cr, d) gmx_bcast( sizeof(d), &(d),(cr))
54 #define nblock_bc(cr,nr,d) gmx_bcast((nr)*sizeof((d)[0]), (d),(cr))
55 #define snew_bc(cr,d,nr) { if (!MASTER(cr)) snew((d),(nr)); }
56 /* #define TAKETIME */
59 ddnoSEL, ddnoINTERLEAVE, ddnoPP_PME, ddnoCARTESIAN, ddnoNR
62 /* Enum for situations that can occur during log file parsing */
68 eParselogResetProblem,
75 int nPMEnodes; /* number of PME only nodes used in this test */
76 int nx, ny, nz; /* DD grid */
77 int guessPME; /* if nPMEnodes == -1, this is the guessed number of PME nodes */
78 float *Gcycles; /* This can contain more than one value if doing multiple tests */
82 float *PME_f_load; /* PME mesh/force load average*/
83 float PME_f_load_Av; /* Average average ;) ... */
84 char *mdrun_cmd_line; /* Mdrun command line used for this test */
90 gmx_large_int_t orig_sim_steps; /* Number of steps to be done in the real simulation */
91 int n_entries; /* Number of entries in arrays */
92 real volume; /* The volume of the box */
93 matrix recipbox; /* The reciprocal box */
94 int natoms; /* The number of atoms in the MD system */
95 real *fac; /* The scaling factor */
96 real *rcoulomb; /* The coulomb radii [0...nr_inputfiles] */
97 real *rvdw; /* The vdW radii */
98 int *nkx, *nky, *nkz; /* Number of k vectors in each spatial dimension */
99 real *fourier_sp; /* Fourierspacing */
100 real *ewald_rtol; /* Real space tolerance for Ewald, determines */
101 /* the real/reciprocal space relative weight */
102 real *ewald_beta; /* Splitting parameter [1/nm] */
103 real fracself; /* fraction of particles for SI error */
104 real q2all; /* sum ( q ^2 ) */
105 real q2allnr; /* nr of charges */
106 int *pme_order; /* Interpolation order for PME (bsplines) */
107 char **fn_out; /* Name of the output tpr file */
108 real *e_dir; /* Direct space part of PME error with these settings */
109 real *e_rec; /* Reciprocal space part of PME error */
110 gmx_bool bTUNE; /* flag for tuning */
114 /* Returns TRUE when atom is charged */
115 static gmx_bool is_charge(real charge)
117 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 */
133 real nrq_all_r; /* No of charges in real format */
135 gmx_moltype_t *molecule;
136 gmx_molblock_t *molblock;
139 fprintf(stderr, "\nCharge density:\n");
141 q2_all = 0.0; /* total q squared */
142 nrq_all = 0; /* total number of charges in the system */
143 for (imol=0; imol<mtop->nmolblock; imol++) /* Loop over molecule types */
145 q2_mol=0.0; /* q squared value of this molecule */
146 nrq_mol=0; /* number of charges this molecule carries */
147 molecule = &(mtop->moltype[imol]);
148 molblock = &(mtop->molblock[imol]);
149 for (iatom=0; iatom<molblock->natoms_mol; iatom++) /* Loop over atoms in this molecule */
151 qi = molecule->atoms.atom[iatom].q; /* Charge of this atom */
152 /* Is this charge worth to be considered? */
159 /* Multiply with the number of molecules present of this type and add */
160 q2_all += q2_mol*molblock->nmol;
161 nrq_all += nrq_mol*molblock->nmol;
163 fprintf(stderr, "Molecule %2d (%5d atoms) q2_mol=%10.3e nr.mol.charges=%5d (%6dx) q2_all=%10.3e tot.charges=%d\n",
164 imol,molblock->natoms_mol,q2_mol,nrq_mol,molblock->nmol,q2_all,nrq_all);
175 /* Estimate the direct space part error of the SPME Ewald sum */
176 static real estimate_direct(
180 real e_dir=0; /* Error estimate */
181 real beta=0; /* Splitting parameter (1/nm) */
182 real r_coulomb=0; /* Cut-off in direct space */
185 beta = info->ewald_beta[0];
186 r_coulomb = info->rcoulomb[0];
188 e_dir = 2.0 * info->q2all * gmx_invsqrt( info->q2allnr * r_coulomb * info->volume );
189 e_dir *= exp (-beta*beta*r_coulomb*r_coulomb);
191 return ONE_4PI_EPS0*e_dir;
196 /* the following 4 functions determine polynomials required for the reciprocal error estimate */
198 static inline real eps_poly1(
199 real m, /* grid coordinate in certain direction */
200 real K, /* grid size in corresponding direction */
201 real n) /* spline interpolation order of the SPME */
204 real nom=0; /* nominator */
205 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 */
246 for(i=-SUMORDER ; i<0 ; i++)
250 nom+=pow( tmp , -2.0*n );
253 for(i=SUMORDER ; i>0 ; i--)
257 nom+=pow( tmp , -2.0*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 */
284 for(i=-SUMORDER ; i<0 ; i++)
288 nom+= i * pow( tmp , -2.0*n );
291 for(i=SUMORDER ; i>0 ; i--)
295 nom+= i * pow( tmp , -2.0*n );
298 for(i=-SUMORDER ; i<SUMORDER+1 ; i++)
302 denom+=pow( tmp , -n );
305 return 2.0 * M_PI * nom / denom / denom;
309 static inline real eps_poly4(
310 real m, /* grid coordinate in certain direction */
311 real K, /* grid size in corresponding direction */
312 real n) /* spline interpolation order of the SPME */
315 real nom=0; /* nominator */
316 real denom=0; /* denominator */
322 for(i=-SUMORDER ; i<0 ; i++)
326 nom+= i * i * pow( tmp , -2.0*n );
329 for(i=SUMORDER ; i>0 ; i--)
333 nom+= i * i * pow( tmp , -2.0*n );
336 for(i=-SUMORDER ; i<SUMORDER+1 ; i++)
340 denom+=pow( tmp , -n );
343 return 4.0 * M_PI * M_PI * nom / denom / denom;
347 static inline real eps_self(
348 real m, /* grid coordinate in certain direction */
349 real K, /* grid size in corresponding direction */
350 rvec rboxv, /* reciprocal box vector */
351 real n, /* spline interpolation order of the SPME */
352 rvec x) /* coordinate of charge */
355 real tmp=0; /* temporary variables for computations */
356 real tmp1=0; /* temporary variables for computations */
357 real tmp2=0; /* temporary variables for computations */
358 real rcoord=0; /* coordinate in certain reciprocal space direction */
359 real nom=0; /* nominator */
360 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,-1.0*n);
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,-1.0*n);
388 tmp=2.0 * M_PI * m / K;
389 tmp1=pow(tmp,-1.0*n);
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;
409 recipbox[YY][XX]=-box[YY][XX]*box[ZZ][ZZ]*tmp;
410 recipbox[YY][YY]=box[XX][XX]*box[ZZ][ZZ]*tmp;
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 */
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 */
457 /* Index variables for parallel work distribution */
458 int startglobal,stopglobal;
459 int startlocal, stoplocal;
468 rng=gmx_rng_init(seed);
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 = ceil((real)xtot / (real)cr->nnodes);
488 startlocal = startglobal + x_per_core*cr->nodeid;
489 stoplocal = startlocal + x_per_core -1;
490 if (stoplocal > stopglobal)
491 stoplocal = stopglobal;
495 startlocal = startglobal;
496 stoplocal = stopglobal;
501 MPI_Barrier(MPI_COMM_WORLD);
514 fprintf(stderr, "Calculating reciprocal error part 1 ...");
518 for(nx=startlocal; nx<=stoplocal; nx++)
520 svmul(nx,info->recipbox[XX],gridpx);
521 for(ny=-info->nky[0]/2; ny<info->nky[0]/2+1; ny++)
523 svmul(ny,info->recipbox[YY],tmpvec);
524 rvec_add(gridpx,tmpvec,gridpxy);
525 for(nz=-info->nkz[0]/2; nz<info->nkz[0]/2+1; nz++)
527 if ( 0 == nx && 0 == ny && 0 == nz )
529 svmul(nz,info->recipbox[ZZ],tmpvec);
530 rvec_add(gridpxy,tmpvec,gridp);
532 coeff=exp(-1.0 * M_PI * M_PI * tmp / info->ewald_beta[0] / info->ewald_beta[0] ) ;
533 coeff/= 2.0 * M_PI * info->volume * tmp;
537 tmp=eps_poly2(nx,info->nkx[0],info->pme_order[0]);
538 tmp+=eps_poly2(ny,info->nkx[0],info->pme_order[0]);
539 tmp+=eps_poly2(nz,info->nkx[0],info->pme_order[0]);
541 tmp1=eps_poly1(nx,info->nkx[0],info->pme_order[0]);
542 tmp2=eps_poly1(ny,info->nky[0],info->pme_order[0]);
544 tmp+=2.0 * tmp1 * tmp2;
546 tmp1=eps_poly1(nz,info->nkz[0],info->pme_order[0]);
547 tmp2=eps_poly1(ny,info->nky[0],info->pme_order[0]);
549 tmp+=2.0 * tmp1 * tmp2;
551 tmp1=eps_poly1(nz,info->nkz[0],info->pme_order[0]);
552 tmp2=eps_poly1(nx,info->nkx[0],info->pme_order[0]);
554 tmp+=2.0 * tmp1 * tmp2;
556 tmp1=eps_poly1(nx,info->nkx[0],info->pme_order[0]);
557 tmp1+=eps_poly1(ny,info->nky[0],info->pme_order[0]);
558 tmp1+=eps_poly1(nz,info->nkz[0],info->pme_order[0]);
562 e_rec1+= 32.0 * M_PI * M_PI * coeff * coeff * coeff2 * tmp * q2_all * q2_all / nr ;
564 tmp1=eps_poly3(nx,info->nkx[0],info->pme_order[0]);
566 tmp2=iprod(gridp,info->recipbox[XX]);
570 tmp1=eps_poly3(ny,info->nky[0],info->pme_order[0]);
572 tmp2=iprod(gridp,info->recipbox[YY]);
576 tmp1=eps_poly3(nz,info->nkz[0],info->pme_order[0]);
578 tmp2=iprod(gridp,info->recipbox[ZZ]);
584 tmp1=eps_poly4(nx,info->nkx[0],info->pme_order[0]);
585 tmp1*=norm2(info->recipbox[XX]);
586 tmp1*=info->nkx[0] * info->nkx[0];
590 tmp1=eps_poly4(ny,info->nky[0],info->pme_order[0]);
591 tmp1*=norm2(info->recipbox[YY]);
592 tmp1*=info->nky[0] * info->nky[0];
596 tmp1=eps_poly4(nz,info->nkz[0],info->pme_order[0]);
597 tmp1*=norm2(info->recipbox[ZZ]);
598 tmp1*=info->nkz[0] * info->nkz[0];
602 e_rec2+= 4.0 * coeff * coeff * tmp * q2_all * q2_all / nr ;
607 fprintf(stderr, "\rCalculating reciprocal error part 1 ... %3.0f%%", 100.0*(nx-startlocal+1)/(x_per_core));
612 fprintf(stderr, "\n");
614 /* Use just a fraction of all charges to estimate the self energy error term? */
615 bFraction = (info->fracself > 0.0) && (info->fracself < 1.0);
619 /* Here xtot is the number of samples taken for the Monte Carlo calculation
620 * of the average of term IV of equation 35 in Wang2010. Round up to a
621 * number of samples that is divisible by the number of nodes */
622 x_per_core = ceil(info->fracself * nr / (real)cr->nnodes);
623 xtot = x_per_core * cr->nnodes;
627 /* In this case we use all nr particle positions */
629 x_per_core = ceil( (real)xtot / (real)cr->nnodes );
632 startlocal = x_per_core * cr->nodeid;
633 stoplocal = min(startlocal + x_per_core, xtot); /* min needed if xtot == nr */
637 /* Make shure we get identical results in serial and parallel. Therefore,
638 * take the sample indices from a single, global random number array that
639 * is constructed on the master node and that only depends on the seed */
643 for (i=0; i<xtot; i++)
645 numbers[i] = floor(gmx_rng_uniform_real(rng) * nr );
648 /* Broadcast the random number array to the other nodes */
651 nblock_bc(cr,xtot,numbers);
654 if (bVerbose && MASTER(cr))
656 fprintf(stdout, "Using %d sample%s to approximate the self interaction error term",
657 xtot, xtot==1?"":"s");
659 fprintf(stdout, " (%d sample%s per node)", x_per_core, x_per_core==1?"":"s");
660 fprintf(stdout, ".\n");
664 /* Return the number of positions used for the Monte Carlo algorithm */
667 for(i=startlocal;i<stoplocal;i++)
675 /* Randomly pick a charge */
680 /* Use all charges */
684 /* for(nx=startlocal; nx<=stoplocal; nx++)*/
685 for(nx=-info->nkx[0]/2; nx<info->nkx[0]/2+1; nx++)
687 svmul(nx,info->recipbox[XX],gridpx);
688 for(ny=-info->nky[0]/2; ny<info->nky[0]/2+1; ny++)
690 svmul(ny,info->recipbox[YY],tmpvec);
691 rvec_add(gridpx,tmpvec,gridpxy);
692 for(nz=-info->nkz[0]/2; nz<info->nkz[0]/2+1; nz++)
695 if ( 0 == nx && 0 == ny && 0 == nz)
698 svmul(nz,info->recipbox[ZZ],tmpvec);
699 rvec_add(gridpxy,tmpvec,gridp);
701 coeff=exp(-1.0 * M_PI * M_PI * tmp / info->ewald_beta[0] / info->ewald_beta[0] );
703 e_rec3x+=coeff*eps_self(nx,info->nkx[0],info->recipbox[XX],info->pme_order[0],x[ci]);
704 e_rec3y+=coeff*eps_self(ny,info->nky[0],info->recipbox[YY],info->pme_order[0],x[ci]);
705 e_rec3z+=coeff*eps_self(nz,info->nkz[0],info->recipbox[ZZ],info->pme_order[0],x[ci]);
713 svmul(e_rec3x,info->recipbox[XX],tmpvec);
714 rvec_inc(tmpvec2,tmpvec);
715 svmul(e_rec3y,info->recipbox[YY],tmpvec);
716 rvec_inc(tmpvec2,tmpvec);
717 svmul(e_rec3z,info->recipbox[ZZ],tmpvec);
718 rvec_inc(tmpvec2,tmpvec);
720 e_rec3 += q[ci]*q[ci]*q[ci]*q[ci]*norm2(tmpvec2) / ( xtot * M_PI * info->volume * M_PI * info->volume);
722 fprintf(stderr, "\rCalculating reciprocal error part 2 ... %3.0f%%",
723 100.0*(i+1)/stoplocal);
729 fprintf(stderr, "\n");
735 t1= MPI_Wtime() - t0;
736 fprintf(fp_out, "Recip. err. est. took : %lf s\n", t1);
744 fprintf(stderr, "Node %3d: nx=[%3d...%3d] e_rec3=%e\n",
745 cr->nodeid, startlocal, stoplocal, e_rec3);
751 gmx_sum(1,&e_rec1,cr);
752 gmx_sum(1,&e_rec2,cr);
753 gmx_sum(1,&e_rec3,cr);
756 /* e_rec1*=8.0 * q2_all / info->volume / info->volume / nr ;
757 e_rec2*= q2_all / M_PI / M_PI / info->volume / info->volume / nr ;
758 e_rec3/= M_PI * M_PI * info->volume * info->volume * nr ;
760 e_rec=sqrt(e_rec1+e_rec2+e_rec3);
763 return ONE_4PI_EPS0 * e_rec;
767 /* Allocate memory for the inputinfo struct: */
768 static void create_info(t_inputinfo *info)
770 snew(info->fac , info->n_entries);
771 snew(info->rcoulomb , info->n_entries);
772 snew(info->rvdw , info->n_entries);
773 snew(info->nkx , info->n_entries);
774 snew(info->nky , info->n_entries);
775 snew(info->nkz , info->n_entries);
776 snew(info->fourier_sp, info->n_entries);
777 snew(info->ewald_rtol, info->n_entries);
778 snew(info->ewald_beta, info->n_entries);
779 snew(info->pme_order , info->n_entries);
780 snew(info->fn_out , info->n_entries);
781 snew(info->e_dir , info->n_entries);
782 snew(info->e_rec , info->n_entries);
786 /* Allocate and fill an array with coordinates and charges,
787 * returns the number of charges found
789 static int prepare_x_q(real *q[], rvec *x[], gmx_mtop_t *mtop, rvec x_orig[], t_commrec *cr)
792 int nq; /* number of charged particles */
798 snew(*q, mtop->natoms);
799 snew(*x, mtop->natoms);
801 for (i=0; i<mtop->natoms; i++)
804 gmx_mtop_atomnr_to_atom(mtop,anr_global,&atom);
805 if (is_charge(atom->q))
808 (*x)[nq][XX] = x_orig[i][XX];
809 (*x)[nq][YY] = x_orig[i][YY];
810 (*x)[nq][ZZ] = x_orig[i][ZZ];
814 /* Give back some unneeded memory */
818 /* Broadcast x and q in the parallel case */
821 /* Transfer the number of charges */
834 /* Read in the tpr file and save information we need later in info */
835 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)
837 read_tpx_state(fn_sim_tpr,ir,state,NULL,mtop);
839 /* The values of the original tpr input file are save in the first
840 * place [0] of the arrays */
841 info->orig_sim_steps = ir->nsteps;
842 info->pme_order[0] = ir->pme_order;
843 info->rcoulomb[0] = ir->rcoulomb;
844 info->rvdw[0] = ir->rvdw;
845 info->nkx[0] = ir->nkx;
846 info->nky[0] = ir->nky;
847 info->nkz[0] = ir->nkz;
848 info->ewald_rtol[0] = ir->ewald_rtol;
849 info->fracself = fracself;
851 info->ewald_beta[0] = user_beta;
853 info->ewald_beta[0] = calc_ewaldcoeff(info->rcoulomb[0],info->ewald_rtol[0]);
855 /* Check if PME was chosen */
856 if (EEL_PME(ir->coulombtype) == FALSE)
857 gmx_fatal(FARGS, "Can only do optimizations for simulations with PME");
859 /* Check if rcoulomb == rlist, which is necessary for PME */
860 if (!(ir->rcoulomb == ir->rlist))
861 gmx_fatal(FARGS, "PME requires rcoulomb (%f) to be equal to rlist (%f).", ir->rcoulomb, ir->rlist);
865 /* Transfer what we need for parallelizing the reciprocal error estimate */
866 static void bcast_info(t_inputinfo *info, t_commrec *cr)
868 nblock_bc(cr, info->n_entries, info->nkx);
869 nblock_bc(cr, info->n_entries, info->nky);
870 nblock_bc(cr, info->n_entries, info->nkz);
871 nblock_bc(cr, info->n_entries, info->ewald_beta);
872 nblock_bc(cr, info->n_entries, info->pme_order);
873 nblock_bc(cr, info->n_entries, info->e_dir);
874 nblock_bc(cr, info->n_entries, info->e_rec);
875 block_bc(cr, info->volume);
876 block_bc(cr, info->recipbox);
877 block_bc(cr, info->natoms);
878 block_bc(cr, info->fracself);
879 block_bc(cr, info->bTUNE);
880 block_bc(cr, info->q2all);
881 block_bc(cr, info->q2allnr);
885 /* Estimate the error of the SPME Ewald sum. This estimate is based upon
886 * a) a homogeneous distribution of the charges
887 * b) a total charge of zero.
889 static void estimate_PME_error(t_inputinfo *info, t_state *state,
890 gmx_mtop_t *mtop, FILE *fp_out, gmx_bool bVerbose, unsigned int seed,
893 rvec *x=NULL; /* The coordinates */
894 real *q=NULL; /* The charges */
895 real edir=0.0; /* real space error */
896 real erec=0.0; /* reciprocal space error */
897 real derr=0.0; /* difference of real and reciprocal space error */
898 real derr0=0.0; /* difference of real and reciprocal space error */
899 real beta=0.0; /* splitting parameter beta */
900 real beta0=0.0; /* splitting parameter beta */
901 int ncharges; /* The number of atoms with charges */
902 int nsamples; /* The number of samples used for the calculation of the
903 * self-energy error term */
907 fprintf(fp_out, "\n--- PME ERROR ESTIMATE ---\n");
909 /* Prepare an x and q array with only the charged atoms */
910 ncharges = prepare_x_q(&q, &x, mtop, state->x, cr);
913 calc_q2all(mtop, &(info->q2all), &(info->q2allnr));
914 info->ewald_rtol[0]=gmx_erfc(info->rcoulomb[0]*info->ewald_beta[0]);
915 /* Write some info to log file */
916 fprintf(fp_out, "Box volume : %g nm^3\n", info->volume);
917 fprintf(fp_out, "Number of charged atoms : %d (total atoms %d)\n",ncharges, info->natoms);
918 fprintf(fp_out, "Coulomb radius : %g nm\n", info->rcoulomb[0]);
919 fprintf(fp_out, "Ewald_rtol : %g\n", info->ewald_rtol[0]);
920 fprintf(fp_out, "Ewald parameter beta : %g\n", info->ewald_beta[0]);
921 fprintf(fp_out, "Interpolation order : %d\n", info->pme_order[0]);
922 fprintf(fp_out, "Fourier grid (nx,ny,nz) : %d x %d x %d\n",
923 info->nkx[0],info->nky[0],info->nkz[0]);
929 bcast_info(info, cr);
932 /* Calculate direct space error */
933 info->e_dir[0] = estimate_direct(info);
935 /* Calculate reciprocal space error */
936 info->e_rec[0] = estimate_reciprocal(info, x, q, ncharges, fp_out, bVerbose,
937 seed, &nsamples, cr);
940 bcast_info(info, cr);
944 fprintf(fp_out, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
945 fprintf(fp_out, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
946 fprintf(fp_out, "Self-energy error term was estimated using %d samples\n", nsamples);
948 fprintf(stderr, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
949 fprintf(stderr, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
957 fprintf(stderr,"Starting tuning ...\n");
961 beta0=info->ewald_beta[0];
963 info->ewald_beta[0]+=0.1;
965 info->ewald_beta[0]-=0.1;
966 info->e_dir[0] = estimate_direct(info);
967 info->e_rec[0] = estimate_reciprocal(info, x, q, ncharges, fp_out, bVerbose,
968 seed, &nsamples, cr);
971 bcast_info(info, cr);
977 while ( fabs(derr/min(erec,edir)) > 1e-4)
980 beta=info->ewald_beta[0];
981 beta-=derr*(info->ewald_beta[0]-beta0)/(derr-derr0);
982 beta0=info->ewald_beta[0];
983 info->ewald_beta[0]=beta;
986 info->e_dir[0] = estimate_direct(info);
987 info->e_rec[0] = estimate_reciprocal(info, x, q, ncharges, fp_out, bVerbose,
988 seed, &nsamples, cr);
991 bcast_info(info, cr);
1000 fprintf(stderr,"difference between real and rec. space error (step %d): %g\n",i,fabs(derr));
1001 fprintf(stderr,"old beta: %f\n",beta0);
1002 fprintf(stderr,"new beta: %f\n",beta);
1006 info->ewald_rtol[0]=gmx_erfc(info->rcoulomb[0]*info->ewald_beta[0]);
1010 /* Write some info to log file */
1012 fprintf(fp_out, "========= After tuning ========\n");
1013 fprintf(fp_out, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
1014 fprintf(fp_out, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
1015 fprintf(stderr, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
1016 fprintf(stderr, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
1017 fprintf(fp_out, "Ewald_rtol : %g\n", info->ewald_rtol[0]);
1018 fprintf(fp_out, "Ewald parameter beta : %g\n", info->ewald_beta[0]);
1028 int gmx_pme_error(int argc,char *argv[])
1030 const char *desc[] = {
1031 "[TT]g_pme_error[tt] estimates the error of the electrostatic forces",
1032 "if using the sPME algorithm. The flag [TT]-tune[tt] will determine",
1033 "the splitting parameter such that the error is equally",
1034 "distributed over the real and reciprocal space part.",
1035 "The part of the error that stems from self interaction of the particles "
1036 "is computationally demanding. However, a good a approximation is to",
1037 "just use a fraction of the particles for this term which can be",
1038 "indicated by the flag [TT]-self[tt].[PAR]",
1041 real fs=0.0; /* 0 indicates: not set by the user */
1042 real user_beta=-1.0;
1045 t_state state; /* The state from the tpr input file */
1046 gmx_mtop_t mtop; /* The topology from the tpr input file */
1047 t_inputrec *ir=NULL; /* The inputrec from the tpr file */
1050 unsigned long PCA_Flags;
1051 gmx_bool bTUNE=FALSE;
1052 gmx_bool bVerbose=FALSE;
1056 static t_filenm fnm[] = {
1057 { efTPX, "-s", NULL, ffREAD },
1058 { efOUT, "-o", "error", ffWRITE },
1059 { efTPX, "-so", "tuned", ffOPTWR }
1062 output_env_t oenv=NULL;
1065 { "-beta", FALSE, etREAL, {&user_beta},
1066 "If positive, overwrite ewald_beta from [TT].tpr[tt] file with this value" },
1067 { "-tune", FALSE, etBOOL, {&bTUNE},
1068 "Tune the splitting parameter such that the error is equally distributed between real and reciprocal space" },
1069 { "-self", FALSE, etREAL, {&fracself},
1070 "If between 0.0 and 1.0, determine self interaction error from just this fraction of the charged particles" },
1071 { "-seed", FALSE, etINT, {&seed},
1072 "Random number seed used for Monte Carlo algorithm when [TT]-self[tt] is set to a value between 0.0 and 1.0" },
1073 { "-v", FALSE, etBOOL, {&bVerbose},
1074 "Be loud and noisy" }
1078 #define NFILE asize(fnm)
1080 cr = init_par(&argc,&argv);
1083 MPI_Barrier(MPI_COMM_WORLD);
1087 CopyRight(stderr,argv[0]);
1089 PCA_Flags = PCA_NOEXIT_ON_ARGS;
1090 PCA_Flags |= (MASTER(cr) ? 0 : PCA_QUIET);
1092 parse_common_args(&argc,argv,PCA_Flags,
1093 NFILE,fnm,asize(pa),pa,asize(desc),desc,
1097 bTUNE = opt2bSet("-so",NFILE,fnm);
1101 /* Allocate memory for the inputinfo struct: */
1103 info.fourier_sp[0] = fs;
1105 /* Read in the tpr file and open logfile for reading */
1109 read_tpr_file(opt2fn("-s",NFILE,fnm), &info, &state, &mtop, ir, user_beta,fracself);
1111 fp=fopen(opt2fn("-o",NFILE,fnm),"w");
1114 /* Check consistency if the user provided fourierspacing */
1115 if (fs > 0 && MASTER(cr))
1117 /* Recalculate the grid dimensions using fourierspacing from user input */
1121 calc_grid(stdout,state.box,info.fourier_sp[0],&(info.nkx[0]),&(info.nky[0]),&(info.nkz[0]));
1122 if ( (ir->nkx != info.nkx[0]) || (ir->nky != info.nky[0]) || (ir->nkz != info.nkz[0]) )
1123 gmx_fatal(FARGS, "Wrong fourierspacing %f nm, input file grid = %d x %d x %d, computed grid = %d x %d x %d",
1124 fs,ir->nkx,ir->nky,ir->nkz,info.nkx[0],info.nky[0],info.nkz[0]);
1127 /* Estimate (S)PME force error */
1129 /* Determine the volume of the simulation box */
1132 info.volume = det(state.box);
1133 calc_recipbox(state.box,info.recipbox);
1134 info.natoms = mtop.natoms;
1139 bcast_info(&info, cr);
1141 /* Get an error estimate of the input tpr file and do some tuning if requested */
1142 estimate_PME_error(&info, &state, &mtop, fp, bVerbose, seed, cr);
1146 /* Write out optimized tpr file if requested */
1147 if ( opt2bSet("-so",NFILE,fnm) || bTUNE )
1149 ir->ewald_rtol=info.ewald_rtol[0];
1150 write_tpx_state(opt2fn("-so",NFILE,fnm),ir,&state,&mtop);
1152 please_cite(fp,"Wang2010");