From 9ffb129916e065182fba1bae02fc7cfa72f3b6d6 Mon Sep 17 00:00:00 2001 From: Rossen Apostolov Date: Thu, 26 Aug 2010 22:51:48 +0200 Subject: [PATCH] Added the g_pme_error tool from Florian Dommert. Thanks! --- src/gmxlib/copyrite.c | 8 +- src/tools/CMakeLists.txt | 5 +- src/tools/Makefile.am | 4 +- src/tools/g_pme_error.c | 50 ++ src/tools/gmx_pme_error.c | 1177 +++++++++++++++++++++++++++++++++++++ 5 files changed, 1238 insertions(+), 6 deletions(-) create mode 100644 src/tools/g_pme_error.c create mode 100644 src/tools/gmx_pme_error.c diff --git a/src/gmxlib/copyrite.c b/src/gmxlib/copyrite.c index 8e4abfead0..b3c8e053c3 100644 --- a/src/gmxlib/copyrite.c +++ b/src/gmxlib/copyrite.c @@ -507,7 +507,13 @@ void please_cite(FILE *fp,const char *key) "O. Engin, A. Villa, M. Sayar and B. Hess", "Driving Forces for Adsorption of Amphiphilic Peptides to Air-Water Interface", "J. Phys. Chem. B", - 0, 2010, "???" } + 0, 2010, "???" }, + { "Wang2010", + "H. Wang, F. Dommert, C.Holm", + "Optimizing working parameters of the smooth particle mesh Ewald algorithm in terms of accuracy and efficiency", + "J. Chem. Phys. B", + 133, 2010, "034117" + } }; #define NSTR (int)asize(citedb) diff --git a/src/tools/CMakeLists.txt b/src/tools/CMakeLists.txt index 3b5dceadad..fe631f0711 100644 --- a/src/tools/CMakeLists.txt +++ b/src/tools/CMakeLists.txt @@ -29,7 +29,7 @@ add_library(gmxana gmx_editconf.c gmx_genbox.c gmx_genion.c gmx_genconf.c gmx_genpr.c gmx_eneconv.c gmx_vanhove.c gmx_wheel.c addconf.c calcpot.c edittop.c gmx_bar.c - gmx_membed.c ) + gmx_membed.c gmx_pme_error.c ) target_link_libraries(gmxana md gmx) @@ -51,8 +51,7 @@ set(GMX_TOOLS_PROGRAMS g_rmsf g_rotacf g_saltbr g_sas g_select g_sgangle g_sham g_sorient g_spol g_sdf g_spatial g_tcaf g_traj g_tune_pme g_vanhove g_velacc g_clustsize g_mdmat g_wham g_sigeps g_bar - g_membed - ) + g_membed g_pme_error) diff --git a/src/tools/Makefile.am b/src/tools/Makefile.am index d21cedc68f..eae8cb9cf5 100644 --- a/src/tools/Makefile.am +++ b/src/tools/Makefile.am @@ -37,7 +37,7 @@ libgmxana@LIBSUFFIX@_la_SOURCES = \ gmx_polystat.c gmx_potential.c gmx_rama.c \ gmx_rdf.c gmx_rms.c gmx_rmsdist.c gmx_rmsf.c \ gmx_rotacf.c gmx_rotmat.c gmx_saltbr.c gmx_sas.c \ - gmx_select.c gmx_sdf.c \ + gmx_select.c gmx_sdf.c gmx_pme_error.c \ gmx_sgangle.c gmx_sorient.c gmx_spol.c gmx_tcaf.c \ gmx_traj.c gmx_velacc.c gmx_helixorient.c \ gmx_clustsize.c gmx_mdmat.c gmx_wham.c eigio.h \ @@ -70,7 +70,7 @@ bin_PROGRAMS = \ g_rotacf g_rotmat g_saltbr g_sas \ g_select g_sgangle \ g_sham g_sorient g_spol \ - g_sdf g_spatial \ + g_sdf g_spatial g_pme_error \ g_tcaf g_traj g_tune_pme \ g_vanhove g_velacc g_membed \ g_clustsize g_mdmat g_wham \ diff --git a/src/tools/g_pme_error.c b/src/tools/g_pme_error.c new file mode 100644 index 0000000000..6eb36b1daa --- /dev/null +++ b/src/tools/g_pme_error.c @@ -0,0 +1,50 @@ +/* + * + * This source code is part of + * + * G R O M A C S + * + * GROningen MAchine for Chemical Simulations + * + * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others. + * Copyright (c) 1991-2000, University of Groningen, The Netherlands. + * Copyright (c) 2001-2008, The GROMACS development team, + * check out http://www.gromacs.org for more information. + + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. + * + * If you want to redistribute modifications, please consider that + * scientific software is very special. Version control is crucial - + * bugs must be traceable. We will be happy to consider code for + * inclusion in the official distribution, but derived work must not + * be called official GROMACS. Details are found in the README & COPYING + * files - if they are missing, get the official version at www.gromacs.org. + * + * To help us fund GROMACS development, we humbly ask that you cite + * the papers on the package - you can find them in the top README file. + * + * For more info, check our website at http://www.gromacs.org + * + * And Hey: + * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include + + +/* This is just a wrapper binary. +*/ +int +main(int argc, char *argv[]) +{ + gmx_pme_error(argc,argv); + return 0; +} + diff --git a/src/tools/gmx_pme_error.c b/src/tools/gmx_pme_error.c new file mode 100644 index 0000000000..5d78ea4896 --- /dev/null +++ b/src/tools/gmx_pme_error.c @@ -0,0 +1,1177 @@ +/* $Id: gmx_tune_pme.c 9 2009-08-11 09:43:30Z dommert $ + * + * This source code is part of + * + * G R O M A C S + * + * GROningen MAchine for Chemical Simulations + * + * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others. + * Copyright (c) 1991-2000, University of Groningen, The Netherlands. + * Copyright (c) 2001-2008, The GROMACS development team, + * check out http://www.gromacs.org for more information. + + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. + * + * If you want to redistribute modifications, please consider that + * scientific software is very special. Version control is crucial - + * bugs must be traceable. We will be happy to consider code for + * inclusion in the official distribution, but derived work must not + * be called official GROMACS. Details are found in the README & COPYING + * files - if they are missing, get the official version at www.gromacs.org. + * + * To help us fund GROMACS development, we humbly ask that you cite + * the papers on the package - you can find them in the top README file. + * + * For more info, check our website at http://www.gromacs.org + * + * And Hey: + * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon + */ +#include "statutil.h" +#include "typedefs.h" +#include "smalloc.h" +#include "vec.h" +#include "copyrite.h" +#include "tpxio.h" +#include "string2.h" +#include "readinp.h" +#include "calcgrid.h" +#include "checkpoint.h" +#include "gmx_ana.h" +#include "gmx_random.h" +#include "physics.h" +#include "mdatoms.h" +#include "coulomb.h" +#include "mtop_util.h" +#include "network.h" + +/* We use the same defines as in mvdata.c here */ +#define block_bc(cr, d) gmx_bcast( sizeof(d), &(d),(cr)) +#define nblock_bc(cr,nr,d) gmx_bcast((nr)*sizeof((d)[0]), (d),(cr)) +#define snew_bc(cr,d,nr) { if (!MASTER(cr)) snew((d),(nr)); } +/* #define TAKETIME */ +/* #define DEBUG */ +enum { + ddnoSEL, ddnoINTERLEAVE, ddnoPP_PME, ddnoCARTESIAN, ddnoNR +}; + +/* Enum for situations that can occur during log file parsing */ +enum { + eParselogOK, + eParselogNotFound, + eParselogNoPerfData, + eParselogTerm, + eParselogResetProblem, + eParselogNr +}; + + +typedef struct +{ + int nPMEnodes; /* number of PME only nodes used in this test */ + int nx, ny, nz; /* DD grid */ + int guessPME; /* if nPMEnodes == -1, this is the guessed number of PME nodes */ + float *Gcycles; /* This can contain more than one value if doing multiple tests */ + float Gcycles_Av; + float *ns_per_day; + float ns_per_day_Av; + float *PME_f_load; /* PME mesh/force load average*/ + float PME_f_load_Av; /* Average average ;) ... */ + char *mdrun_cmd_line; /* Mdrun command line used for this test */ +} t_perf; + + +typedef struct +{ + gmx_large_int_t orig_sim_steps; /* Number of steps to be done in the real simulation */ + int n_entries; /* Number of entries in arrays */ + real volume; /* The volume of the box */ + matrix recipbox; /* The reciprocal box */ + int natoms; /* The number of atoms in the MD system */ + real *fac; /* The scaling factor */ + real *rcoulomb; /* The coulomb radii [0...nr_inputfiles] */ + real *rvdw; /* The vdW radii */ + int *nkx, *nky, *nkz; /* Number of k vectors in each spatial dimension */ + real *fourier_sp; /* Fourierspacing */ + real *ewald_rtol; /* Real space tolerance for Ewald, determines */ + /* the real/reciprocal space relative weight */ + real *ewald_beta; /* Splitting parameter [1/nm] */ + real fracself; /* fraction of particles for SI error */ + real q2all; /* sum ( q ^2 ) */ + real q2allnr; /* nr of charges */ + int *pme_order; /* Interpolation order for PME (bsplines) */ + char **fn_out; /* Name of the output tpr file */ + real *e_dir; /* Direct space part of PME error with these settings */ + real *e_rec; /* Reciprocal space part of PME error */ + gmx_bool bTUNE; /* flag for tuning */ +} t_inputinfo; + + +static void sep_line(FILE *fp) +{ + fprintf(fp, "\n------------------------------------------------------------\n"); +} + + + +static gmx_bool is_equal(real a, real b) +{ + real diff, eps=1.0e-6; + + + diff = a - b; + + if (diff < 0.0) diff = -diff; + + if (diff < eps) + return TRUE; + else + return FALSE; +} + + +/* Returns TRUE when atom is charged */ +static gmx_bool is_charge(real charge) +{ + if (charge*charge > GMX_REAL_EPS) + return TRUE; + else + return FALSE; +} + + +/* calculate charge density */ +static void calc_q2all( + gmx_mtop_t *mtop, /* molecular topology */ + real *q2all, real *q2allnr) +{ + int imol,iatom; /* indices for loops */ + real q2_all=0; /* Sum of squared charges */ + int nrq_mol; /* Number of charges in a single molecule */ + int nrq_all; /* Total number of charges in the MD system */ + real nrq_all_r; /* No of charges in real format */ + real qi,q2_mol; + gmx_moltype_t *molecule; + gmx_molblock_t *molblock; + +#ifdef DEBUG + fprintf(stderr, "\nCharge density:\n"); +#endif + q2_all = 0.0; /* total q squared */ + nrq_all = 0; /* total number of charges in the system */ + for (imol=0; imolnmolblock; imol++) /* Loop over molecule types */ + { + q2_mol=0.0; /* q squared value of this molecule */ + nrq_mol=0; /* number of charges this molecule carries */ + molecule = &(mtop->moltype[imol]); + molblock = &(mtop->molblock[imol]); + for (iatom=0; iatomnatoms_mol; iatom++) /* Loop over atoms in this molecule */ + { + qi = molecule->atoms.atom[iatom].q; /* Charge of this atom */ + /* Is this charge worth to be considered? */ + if (is_charge(qi)) + { + q2_mol += qi*qi; + nrq_mol++; + } + } + /* Multiply with the number of molecules present of this type and add */ + q2_all += q2_mol*molblock->nmol; + nrq_all += nrq_mol*molblock->nmol; +#ifdef DEBUG + fprintf(stderr, "Molecule %2d (%5d atoms) q2_mol=%10.3e nr.mol.charges=%5d (%6dx) q2_all=%10.3e tot.charges=%d\n", + imol,molblock->natoms_mol,q2_mol,nrq_mol,molblock->nmol,q2_all,nrq_all); +#endif + } + nrq_all_r = nrq_all; + + *q2all=q2_all; + *q2allnr=nrq_all; + +} + + +/* Estimate the direct space part error of the SPME Ewald sum */ +static real estimate_direct( + t_inputinfo *info + ) +{ + real e_dir=0; /* Error estimate */ + real beta=0; /* Splitting parameter (1/nm) */ + real r_coulomb=0; /* Cut-off in direct space */ + + + beta = info->ewald_beta[0]; + r_coulomb = info->rcoulomb[0]; + + e_dir = 2.0 * info->q2all * gmx_invsqrt( info->q2allnr * r_coulomb * info->volume ); + e_dir *= exp (-beta*beta*r_coulomb*r_coulomb); + + return ONE_4PI_EPS0*e_dir; +} + +#define SUMORDER 6 + +/* the following 4 functions determine polynomials required for the reciprocal error estimate */ + +static inline real eps_poly1( + real m, /* grid coordinate in certain direction */ + real K, /* grid size in corresponding direction */ + real n) /* spline interpolation order of the SPME */ +{ + int i; + real nom=0; /* nominator */ + real denom=0; /* denominator */ + real tmp=0; + + if ( m == 0.0 ) + return 0.0 ; + + for(i=-SUMORDER ; i<0 ; i++) + { + tmp=m / K + i; + tmp*=2.0*M_PI; + nom+=pow( tmp , -n ); + } + + for(i=SUMORDER ; i>0 ; i--) + { + tmp=m / K + i; + tmp*=2.0*M_PI; + nom+=pow( tmp , -n ); + } + + tmp=m / K; + tmp*=2.0*M_PI; + denom=pow( tmp , -n )+nom; + + return -nom/denom; + +} + +static inline real eps_poly2( + real m, /* grid coordinate in certain direction */ + real K, /* grid size in corresponding direction */ + real n) /* spline interpolation order of the SPME */ +{ + int i; + real nom=0; /* nominator */ + real denom=0; /* denominator */ + real tmp=0; + + if ( m == 0.0 ) + return 0.0 ; + + for(i=-SUMORDER ; i<0 ; i++) + { + tmp=m / K + i; + tmp*=2.0*M_PI; + nom+=pow( tmp , -2.0*n ); + } + + for(i=SUMORDER ; i>0 ; i--) + { + tmp=m / K + i; + tmp*=2.0*M_PI; + nom+=pow( tmp , -2.0*n ); + } + + for(i=-SUMORDER ; i0 ; i--) + { + tmp=m / K + i; + tmp*=2.0*M_PI; + nom+= i * pow( tmp , -2.0*n ); + } + + for(i=-SUMORDER ; i0 ; i--) + { + tmp=m / K + i; + tmp*=2.0*M_PI; + nom+= i * i * pow( tmp , -2.0*n ); + } + + for(i=-SUMORDER ; i0;i--) + { + tmp=-sin(2.0 * M_PI * i * K * rcoord); + tmp1=2.0 * M_PI * m / K + 2.0 * M_PI * i; + tmp2=pow(tmp1,-1.0*n); + nom+=tmp * tmp2 * i; + denom+=tmp2; + } + + + tmp=2.0 * M_PI * m / K; + tmp1=pow(tmp,-1.0*n); + denom+=tmp1; + + return 2.0 * M_PI * nom / denom * K ; + +} + +#undef SUMORDER + +/* The following routine is just a copy from pme.c */ + +static void calc_recipbox(matrix box,matrix recipbox) +{ + /* Save some time by assuming upper right part is zero */ + + real tmp=1.0/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]); + + recipbox[XX][XX]=box[YY][YY]*box[ZZ][ZZ]*tmp; + recipbox[XX][YY]=0; + recipbox[XX][ZZ]=0; + recipbox[YY][XX]=-box[YY][XX]*box[ZZ][ZZ]*tmp; + recipbox[YY][YY]=box[XX][XX]*box[ZZ][ZZ]*tmp; + recipbox[YY][ZZ]=0; + recipbox[ZZ][XX]=(box[YY][XX]*box[ZZ][YY]-box[YY][YY]*box[ZZ][XX])*tmp; + recipbox[ZZ][YY]=-box[ZZ][YY]*box[XX][XX]*tmp; + recipbox[ZZ][ZZ]=box[XX][XX]*box[YY][YY]*tmp; +} + + +/* Estimate the reciprocal space part error of the SPME Ewald sum. */ +static real estimate_reciprocal( + t_inputinfo *info, + rvec x[], /* array of particles */ + real q[], /* array of charges */ + int nr, /* number of charges = size of the charge array */ + FILE *fp_out, + t_commrec *cr) +{ + real e_rec=0; /* reciprocal error estimate */ + real e_rec1=0; /* Error estimate term 1*/ + real e_rec2=0; /* Error estimate term 2*/ + real e_rec3=0; /* Error estimate term 3 */ + real e_rec3x=0; /* part of Error estimate term 3 in x */ + real e_rec3y=0; /* part of Error estimate term 3 in y */ + real e_rec3z=0; /* part of Error estimate term 3 in z */ + int i,ci; + int nx,ny,nz; /* grid coordinates */ + real q2_all=0; /* sum of squared charges */ + rvec gridpx; /* reciprocal grid point in x direction*/ + rvec gridpxy; /* reciprocal grid point in x and y direction*/ + rvec gridp; /* complete reciprocal grid point in 3 directions*/ + rvec tmpvec; /* template to create points from basis vectors */ + rvec tmpvec2; /* template to create points from basis vectors */ + real coeff=0; /* variable to compute coefficients of the error estimate */ + real coeff2=0; /* variable to compute coefficients of the error estimate */ + real tmp=0; /* variables to compute different factors from vectors */ + real tmp1=0; + real tmp2=0; + real xtmp=0; + real ytmp=0; + real ztmp=0; + double ewald_error; + + /* Random number generator */ + + gmx_rng_t rng=NULL; + /*rng=gmx_rng_init(gmx_rng_make_seed()); */ + rng=gmx_rng_init(cr->nodeid); + + /* Index variables for parallel work distribution */ + int startglobal,stopglobal; + int startlocal, stoplocal; + int x_per_core; + int nrsamples; + real xtot; + +/* #define TAKETIME */ +#ifdef TAKETIME + double t0=0.0; + double t1=0.0; + double t2=0.0; +#endif + + clear_rvec(gridpx); + clear_rvec(gridpxy); + clear_rvec(gridp); + clear_rvec(tmpvec); + clear_rvec(tmpvec2); + + for(i=0;inkx[0]/2; + stopglobal = info->nkx[0]/2; + xtot = stopglobal*2+1; + if (PAR(cr)) + { + x_per_core = ceil(xtot / cr->nnodes); + startlocal = startglobal + x_per_core*cr->nodeid; + stoplocal = startlocal + x_per_core -1; + if (stoplocal > stopglobal) + stoplocal = stopglobal; + } + else + { + startlocal = startglobal; + stoplocal = stopglobal; + x_per_core = xtot; + } +/* +#ifdef GMX_MPI + MPI_Barrier(MPI_COMM_WORLD); +#endif +*/ + +#ifdef TAKETIME + if (MASTER(cr)) + t0 = MPI_Wtime(); +#endif + + if (MASTER(cr)){ + + fprintf(stderr, "Calculating reciprocal error part 1 ..."); + + } + + for(nx=startlocal; nx<=stoplocal; nx++) + { + svmul(nx,info->recipbox[XX],gridpx); + for(ny=-info->nky[0]/2; nynky[0]/2+1; ny++) + { + svmul(ny,info->recipbox[YY],tmpvec); + rvec_add(gridpx,tmpvec,gridpxy); + for(nz=-info->nkz[0]/2; nznkz[0]/2+1; nz++) + { + if ( 0 == nx && 0 == ny && 0 == nz ) + continue; + svmul(nz,info->recipbox[ZZ],tmpvec); + rvec_add(gridpxy,tmpvec,gridp); + tmp=norm2(gridp); + coeff=exp(-1.0 * M_PI * M_PI * tmp / info->ewald_beta[0] / info->ewald_beta[0] ) ; + coeff/= 2.0 * M_PI * info->volume * tmp; + coeff2=tmp ; + + + tmp=eps_poly2(nx,info->nkx[0],info->pme_order[0]); + tmp+=eps_poly2(ny,info->nkx[0],info->pme_order[0]); + tmp+=eps_poly2(nz,info->nkx[0],info->pme_order[0]); + + tmp1=eps_poly1(nx,info->nkx[0],info->pme_order[0]); + tmp2=eps_poly1(ny,info->nky[0],info->pme_order[0]); + + tmp+=2.0 * tmp1 * tmp2; + + tmp1=eps_poly1(nz,info->nkz[0],info->pme_order[0]); + tmp2=eps_poly1(ny,info->nky[0],info->pme_order[0]); + + tmp+=2.0 * tmp1 * tmp2; + + tmp1=eps_poly1(nz,info->nkz[0],info->pme_order[0]); + tmp2=eps_poly1(nx,info->nkx[0],info->pme_order[0]); + + tmp+=2.0 * tmp1 * tmp2; + + tmp1=eps_poly1(nx,info->nkx[0],info->pme_order[0]); + tmp1+=eps_poly1(ny,info->nky[0],info->pme_order[0]); + tmp1+=eps_poly1(nz,info->nkz[0],info->pme_order[0]); + + tmp+= tmp1 * tmp1; + + e_rec1+= 32.0 * M_PI * M_PI * coeff * coeff * coeff2 * tmp * q2_all * q2_all / nr ; + + tmp1=eps_poly3(nx,info->nkx[0],info->pme_order[0]); + tmp1*=info->nkx[0]; + tmp2=iprod(gridp,info->recipbox[XX]); + + tmp=tmp1*tmp2; + + tmp1=eps_poly3(ny,info->nky[0],info->pme_order[0]); + tmp1*=info->nky[0]; + tmp2=iprod(gridp,info->recipbox[YY]); + + tmp+=tmp1*tmp2; + + tmp1=eps_poly3(nz,info->nkz[0],info->pme_order[0]); + tmp1*=info->nkz[0]; + tmp2=iprod(gridp,info->recipbox[ZZ]); + + tmp+=tmp1*tmp2; + + tmp*=4.0 * M_PI; + + tmp1=eps_poly4(nx,info->nkx[0],info->pme_order[0]); + tmp1*=norm2(info->recipbox[XX]); + tmp1*=info->nkx[0] * info->nkx[0]; + + tmp+=tmp1; + + tmp1=eps_poly4(ny,info->nky[0],info->pme_order[0]); + tmp1*=norm2(info->recipbox[YY]); + tmp1*=info->nky[0] * info->nky[0]; + + tmp+=tmp1; + + tmp1=eps_poly4(nz,info->nkz[0],info->pme_order[0]); + tmp1*=norm2(info->recipbox[ZZ]); + tmp1*=info->nkz[0] * info->nkz[0]; + + tmp+=tmp1; + + e_rec2+= 4.0 * coeff * coeff * tmp * q2_all * q2_all / nr ; + + } + } + if (MASTER(cr)) + fprintf(stderr, "\rCalculating reciprocal error part 1 ... %3.0f%%", 100.0*(nx-startlocal+1)/(x_per_core)); + + } + +/* +#ifdef GMX_MPI + MPI_Barrier(MPI_COMM_WORLD); +#endif +*/ + + if (MASTER(cr)) + fprintf(stderr, "\n"); + if (info->fracself>0) + { + nrsamples=ceil(info->fracself*nr); + } + else + { + nrsamples=nr; + } + + + xtot=nrsamples; + + + startglobal=0; + stopglobal=nr; + + if(PAR(cr)) + { + x_per_core=ceil(xtot/cr->nnodes); + startlocal=startglobal+x_per_core*cr->nodeid; + stoplocal=startglobal+x_per_core*(cr->nodeid+1); + if (stoplocal>stopglobal) + stoplocal=stopglobal; + } + else + { + startlocal=startglobal; + stoplocal=stopglobal; + x_per_core=xtot; + } + + + + for(i=startlocal;ifracself<0) { + ci=i; + }else { + ci=floor(gmx_rng_uniform_real(rng) * nr ); + if (ci==nr) + { + ci=nr-1; + } + } + + /* for(nx=startlocal; nx<=stoplocal; nx++)*/ + for(nx=-info->nkx[0]/2; nxnkx[0]/2+1; nx++) + { + svmul(nx,info->recipbox[XX],gridpx); + for(ny=-info->nky[0]/2; nynky[0]/2+1; ny++) + { + svmul(ny,info->recipbox[YY],tmpvec); + rvec_add(gridpx,tmpvec,gridpxy); + for(nz=-info->nkz[0]/2; nznkz[0]/2+1; nz++) + { + + if ( 0 == nx && 0 == ny && 0 == nz) + continue; + + svmul(nz,info->recipbox[ZZ],tmpvec); + rvec_add(gridpxy,tmpvec,gridp); + tmp=norm2(gridp); + coeff=exp(-1.0 * M_PI * M_PI * tmp / info->ewald_beta[0] / info->ewald_beta[0] ); + coeff/= tmp ; + e_rec3x+=coeff*eps_self(nx,info->nkx[0],info->recipbox[XX],info->pme_order[0],x[ci]); + e_rec3y+=coeff*eps_self(ny,info->nky[0],info->recipbox[YY],info->pme_order[0],x[ci]); + e_rec3z+=coeff*eps_self(nz,info->nkz[0],info->recipbox[ZZ],info->pme_order[0],x[ci]); + + } + } + } + + clear_rvec(tmpvec2); + + svmul(e_rec3x,info->recipbox[XX],tmpvec); + rvec_inc(tmpvec2,tmpvec); + svmul(e_rec3y,info->recipbox[YY],tmpvec); + rvec_inc(tmpvec2,tmpvec); + svmul(e_rec3z,info->recipbox[ZZ],tmpvec); + rvec_inc(tmpvec2,tmpvec); + + e_rec3 += q[ci]*q[ci]*q[ci]*q[ci]*norm2(tmpvec2) / ( nrsamples * M_PI * info->volume * M_PI * info->volume); + if (MASTER(cr)){ + fprintf(stderr, "\rCalculating reciprocal error part 2 ... %3.0f%%", + 100.0*(i+1)/stoplocal); + + } + + } + + if (MASTER(cr)) + fprintf(stderr, "\n"); + + +#ifdef TAKETIME + if (MASTER(cr)) + { + t1= MPI_Wtime() - t0; + fprintf(fp_out, "Recip. err. est. took : %lf s\n", t1); + } +#endif + +#ifdef DEBUG + if (PAR(cr)) + { + fprintf(stderr, "Node %3d: nx=[%3d...%3d] e_rec3=%e\n", + cr->nodeid, startlocal, stoplocal, e_rec3); + } +#endif + + +/* +#ifdef GMX_MPI + MPI_Barrier(MPI_COMM_WORLD); +#endif + */ + +#ifdef TAKETIME + if (MASTER(cr)) + { + t2= MPI_Wtime() - t0; + fprintf(fp_out, "barrier : %lf s\n", t2-t1); + } +#endif + + if (PAR(cr)) + { + gmx_sum(1,&e_rec1,cr); + gmx_sum(1,&e_rec2,cr); + gmx_sum(1,&e_rec3,cr); + } + +#ifdef TAKETIME + if (MASTER(cr)) + fprintf(fp_out, "final reduce : %lf s\n", MPI_Wtime() - t0-t2); +#endif + /* e_rec1*=8.0 * q2_all / info->volume / info->volume / nr ; + e_rec2*= q2_all / M_PI / M_PI / info->volume / info->volume / nr ; + e_rec3/= M_PI * M_PI * info->volume * info->volume * nr ; + */ + e_rec=sqrt(e_rec1+e_rec2+e_rec3); + + + return ONE_4PI_EPS0 * e_rec; +} + + +/* Allocate memory for the inputinfo struct: */ +static void create_info(t_inputinfo *info) +{ + snew(info->fac , info->n_entries); + snew(info->rcoulomb , info->n_entries); + snew(info->rvdw , info->n_entries); + snew(info->nkx , info->n_entries); + snew(info->nky , info->n_entries); + snew(info->nkz , info->n_entries); + snew(info->fourier_sp, info->n_entries); + snew(info->ewald_rtol, info->n_entries); + snew(info->ewald_beta, info->n_entries); + snew(info->pme_order , info->n_entries); + snew(info->fn_out , info->n_entries); + snew(info->e_dir , info->n_entries); + snew(info->e_rec , info->n_entries); +} + + +/* Allocate and fill an array with coordinates and charges, + * returns the number of charges found + */ +static int prepare_x_q(real *q[], rvec *x[], gmx_mtop_t *mtop, rvec x_orig[], t_commrec *cr) +{ + int i,anr_global; + int nq; /* number of charged particles */ + t_atom *atom; + + + if (MASTER(cr)) + { + snew(*q, mtop->natoms); + snew(*x, mtop->natoms); + nq=0; + for (i=0; inatoms; i++) + { + anr_global = i; + gmx_mtop_atomnr_to_atom(mtop,anr_global,&atom); + if (is_charge(atom->q)) + { + (*q)[nq] = atom->q; + (*x)[nq][XX] = x_orig[i][XX]; + (*x)[nq][YY] = x_orig[i][YY]; + (*x)[nq][ZZ] = x_orig[i][ZZ]; + nq++; + } + } + /* Give back some unneeded memory */ + srenew(*q, nq); + srenew(*x, nq); + } + /* Broadcast x and q in the parallel case */ + if (PAR(cr)) + { + /* Transfer the number of charges */ + block_bc(cr,nq); + snew_bc(cr, *x, nq); + snew_bc(cr, *q, nq); + nblock_bc(cr,nq,*x); + nblock_bc(cr,nq,*q); + } + + return nq; +} + + + +/* Read in the tpr file and save information we need later in info */ +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) +{ + read_tpx_state(fn_sim_tpr,ir,state,NULL,mtop); + + /* The values of the original tpr input file are save in the first + * place [0] of the arrays */ + info->orig_sim_steps = ir->nsteps; + info->pme_order[0] = ir->pme_order; + info->rcoulomb[0] = ir->rcoulomb; + info->rvdw[0] = ir->rvdw; + info->nkx[0] = ir->nkx; + info->nky[0] = ir->nky; + info->nkz[0] = ir->nkz; + info->ewald_rtol[0] = ir->ewald_rtol; + info->fracself = fracself; + if (user_beta > 0) + info->ewald_beta[0] = user_beta; + else + info->ewald_beta[0] = calc_ewaldcoeff(info->rcoulomb[0],info->ewald_rtol[0]); + + /* Check if PME was chosen */ + if (EEL_PME(ir->coulombtype) == FALSE) + gmx_fatal(FARGS, "Can only do optimizations for simulations with PME"); + + /* Check if rcoulomb == rlist, which is necessary for PME */ + if (!(ir->rcoulomb == ir->rlist)) + gmx_fatal(FARGS, "PME requires rcoulomb (%f) to be equal to rlist (%f).", ir->rcoulomb, ir->rlist); +} + + +/* Transfer what we need for parallelizing the reciprocal error estimate */ +static void bcast_info(t_inputinfo *info, t_commrec *cr) +{ + nblock_bc(cr, info->n_entries, info->nkx); + nblock_bc(cr, info->n_entries, info->nky); + nblock_bc(cr, info->n_entries, info->nkz); + nblock_bc(cr, info->n_entries, info->ewald_beta); + nblock_bc(cr, info->n_entries, info->pme_order); + nblock_bc(cr, info->n_entries, info->e_dir); + nblock_bc(cr, info->n_entries, info->e_rec); + block_bc(cr, info->volume); + block_bc(cr, info->recipbox); + block_bc(cr, info->natoms); + block_bc(cr, info->fracself); + block_bc(cr, info->bTUNE); + block_bc(cr, info->q2all); + block_bc(cr, info->q2allnr); +} + + +/* Estimate the error of the SPME Ewald sum. This estimate is based upon + * a) a homogeneous distribution of the charges + * b) a total charge of zero. + */ +static void estimate_PME_error(t_inputinfo *info, t_state *state, + gmx_mtop_t *mtop, FILE *fp_out, t_commrec *cr) +{ + rvec *x=NULL; /* The coordinates */ + real *q=NULL; /* The charges */ + real q2all=0.0; + real q2allnr=0.0; + real edir=0.0; /* real space error */ + real erec=0.0; /* reciprocal space error */ + real derr=0.0; /* difference of real and reciprocal space error */ + real derr0=0.0; /* difference of real and reciprocal space error */ + real beta=0.0; /* splitting parameter beta */ + real beta0=0.0; /* splitting parameter beta */ + int ncharges; /* The number of atoms with charges */ + int i=0; + + if (MASTER(cr)) + fprintf(fp_out, "\n--- PME ERROR ESTIMATE ---\n"); + + /* Prepare an x and q array with only the charged atoms */ + ncharges = prepare_x_q(&q, &x, mtop, state->x, cr); + if (MASTER(cr)) + { + calc_q2all(mtop, &(info->q2all), &(info->q2allnr)); + info->ewald_rtol[0]=gmx_erfc(info->rcoulomb[0]*info->ewald_beta[0]); + /* Write some info to log file */ + fprintf(fp_out, "Box volume : %g nm^3\n", info->volume); + fprintf(fp_out, "Number of charged atoms : %d (total atoms %d)\n",ncharges, info->natoms); + fprintf(fp_out, "Coulomb radius : %g nm\n", info->rcoulomb[0]); + fprintf(fp_out, "Ewald_rtol : %g\n", info->ewald_rtol[0]); + fprintf(fp_out, "Ewald parameter beta : %g\n", info->ewald_beta[0]); + fprintf(fp_out, "Interpolation order : %d\n", info->pme_order[0]); + fprintf(fp_out, "Fourier grid (nx,ny,nz) : %d x %d x %d\n", + info->nkx[0],info->nky[0],info->nkz[0]); + fflush(fp_out); + + } + + if (PAR(cr)) + bcast_info(info, cr); + + + /* Calculate direct space error */ + info->e_dir[0] = estimate_direct(info); + + /* Calculate reciprocal space error */ + info->e_rec[0] = estimate_reciprocal(info, x, q, ncharges, fp_out, cr); + + if (PAR(cr)) + bcast_info(info, cr); + + if (MASTER(cr)) + { + fprintf(fp_out, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]); + fprintf(fp_out, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]); + fflush(fp_out); + fprintf(stderr, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]); + fprintf(stderr, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]); + } + + i=0; + + if (info->bTUNE) + { + if(MASTER(cr)) + fprintf(stderr,"Starting tuning ...\n"); + edir=info->e_dir[0]; + erec=info->e_rec[0]; + derr0=edir-erec; + beta0=info->ewald_beta[0]; + if (derr>0.0) + info->ewald_beta[0]+=0.1; + else + info->ewald_beta[0]-=0.1; + info->e_dir[0] = estimate_direct(info); + info->e_rec[0] = estimate_reciprocal(info, x, q, ncharges, fp_out, cr); + + if (PAR(cr)) + bcast_info(info, cr); + + + edir=info->e_dir[0]; + erec=info->e_rec[0]; + derr=edir-erec; + while ( fabs(derr/min(erec,edir)) > 1e-4) + { + + beta=info->ewald_beta[0]; + beta-=derr*(info->ewald_beta[0]-beta0)/(derr-derr0); + beta0=info->ewald_beta[0]; + info->ewald_beta[0]=beta; + derr0=derr; + + info->e_dir[0] = estimate_direct(info); + info->e_rec[0] = estimate_reciprocal(info, x, q, ncharges, fp_out, cr); + + if (PAR(cr)) + bcast_info(info, cr); + + edir=info->e_dir[0]; + erec=info->e_rec[0]; + derr=edir-erec; + + if (MASTER(cr)) + { + i++; + fprintf(stderr,"difference between real and rec. space error (step %d): %g\n",i,fabs(derr)); + fprintf(stderr,"old beta: %f\n",beta0); + fprintf(stderr,"new beta: %f\n",beta); + } + } + + info->ewald_rtol[0]=gmx_erfc(info->rcoulomb[0]*info->ewald_beta[0]); + + if (MASTER(cr)) + { + /* Write some info to log file */ + fflush(fp_out); + fprintf(fp_out, "========= After tuning ========\n"); + fprintf(fp_out, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]); + fprintf(fp_out, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]); + fprintf(stderr, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]); + fprintf(stderr, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]); + fprintf(fp_out, "Ewald_rtol : %g\n", info->ewald_rtol[0]); + fprintf(fp_out, "Ewald parameter beta : %g\n", info->ewald_beta[0]); + fflush(fp_out); + + } + + } + +} + +#define BENCHSTEPS (1000) + +int gmx_pme_error(int argc,char *argv[]) +{ + const char *desc[] = { + "g_pme_error estimates the error of the electrostatic forces", + "if using the SPME algorithm. The flag [TT]-tune[tt] will determine", + "the splitting parameter such that the error is equally", + "distributed over the real and reciprocal space part.", + "As a part of the error stems from self interaction of the particles " + "and is computationally very demanding a good a approximation is possible", + "if just a fraction of the particles is used to calculate the average", + "of this error by using the flag [TT]-self[tt].[PAR]", + }; + + int repeats=2; + real fs=0.0; /* 0 indicates: not set by the user */ + + real user_beta=-1.0; + real fracself=-1.0; + + + t_perf **perfdata; + t_inputinfo info; + t_state state; /* The state from the tpr input file */ + gmx_mtop_t mtop; /* The topology from the tpr input file */ + t_inputrec *ir=NULL; /* The inputrec from the tpr file */ + FILE *fp=NULL; + t_commrec *cr; + unsigned long PCA_Flags; + gmx_bool bTUNE=FALSE; + + + static t_filenm fnm[] = { + /* g_tune_pme */ + { efTPX, "-s", NULL, ffREAD }, + { efOUT, "-o", "error", ffWRITE }, + { efTPX, "-so", "tuned", ffOPTWR } + }; + + + output_env_t oenv=NULL; + + t_pargs pa[] = { + /***********************/ + /* g_tune_pme options: */ + /***********************/ + { "-beta", FALSE, etREAL, {&user_beta}, + "If positive, overwrite ewald_beta from tpr file with this value" }, + { "-tune", FALSE, etBOOL, {&bTUNE}, + "If flag is set the splitting parameter will be tuned to distribute the error equally in real and rec. space" }, + { "-self", FALSE, etREAL, {&fracself}, + "If positive, determine selfinteraction error just over this fraction (default=1.0)" } + }; + + +#define NFILE asize(fnm) + + cr = init_par(&argc,&argv); + + if (MASTER(cr)) + CopyRight(stderr,argv[0]); + + PCA_Flags = PCA_NOEXIT_ON_ARGS; + PCA_Flags |= (MASTER(cr) ? 0 : PCA_QUIET); + + parse_common_args(&argc,argv,PCA_Flags, + NFILE,fnm,asize(pa),pa,asize(desc),desc, + 0,NULL,&oenv); + + if (!bTUNE) + bTUNE = opt2bSet("-so",NFILE,fnm); + + info.n_entries = 1; + + /* Allocate memory for the inputinfo struct: */ + create_info(&info); + info.fourier_sp[0] = fs; + + /* Read in the tpr file and open logfile for reading */ + if (MASTER(cr)) + { + snew(ir,1); + read_tpr_file(opt2fn("-s",NFILE,fnm), &info, &state, &mtop, ir, user_beta,fracself); + + fp=fopen(opt2fn("-o",NFILE,fnm),"w"); + } + + /* Check consistency if the user provided fourierspacing */ + if (fs > 0 && MASTER(cr)) + { + /* Recalculate the grid dimensions using fourierspacing from user input */ + info.nkx[0] = 0; + info.nky[0] = 0; + info.nkz[0] = 0; + calc_grid(stdout,state.box,info.fourier_sp[0],&(info.nkx[0]),&(info.nky[0]),&(info.nkz[0])); + if ( (ir->nkx != info.nkx[0]) || (ir->nky != info.nky[0]) || (ir->nkz != info.nkz[0]) ) + gmx_fatal(FARGS, "Wrong fourierspacing %f nm, input file grid = %d x %d x %d, computed grid = %d x %d x %d", + fs,ir->nkx,ir->nky,ir->nkz,info.nkx[0],info.nky[0],info.nkz[0]); + } + + /* Estimate (S)PME force error */ + + /* Determine the volume of the simulation box */ + if (MASTER(cr)) + { + info.volume = det(state.box); + calc_recipbox(state.box,info.recipbox); + info.natoms = mtop.natoms; + info.bTUNE = bTUNE; + } + + if (PAR(cr)) + bcast_info(&info, cr); + + /* Get an error estimate of the input tpr file */ + estimate_PME_error(&info, &state, &mtop, fp, cr); + + if (MASTER(cr)) + { + ir->ewald_rtol=info.ewald_rtol[0]; + write_tpx_state(opt2fn("-so",NFILE,fnm),ir,&state,&mtop); + please_cite(fp,"Wang2010"); + fclose(fp); + } + + gmx_finalize(); + + return 0; +} -- 2.22.0