src/tools/gmx_pme_error.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2008, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012,2013, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #include "statutil.h"
  39 #include "typedefs.h"
  40 #include "smalloc.h"
  41 #include "vec.h"
  42 #include "copyrite.h"
  43 #include "tpxio.h"
  44 #include "string2.h"
  45 #include "readinp.h"
  46 #include "calcgrid.h"
  47 #include "checkpoint.h"
  48 #include "gmx_ana.h"
  49 #include "gmx_random.h"
  50 #include "physics.h"
  51 #include "mdatoms.h"
  52 #include "coulomb.h"
  53 #include "mtop_util.h"
  54 #include "network.h"
  55 #include "main.h"
  56
  57 /* We use the same defines as in mvdata.c here */
  58 #define  block_bc(cr,   d) gmx_bcast(     sizeof(d),     &(d),(cr))
  59 #define nblock_bc(cr,nr,d) gmx_bcast((nr)*sizeof((d)[0]), (d),(cr))
  60 #define   snew_bc(cr,d,nr) { if (!MASTER(cr)) snew((d),(nr)); }
  61 /* #define TAKETIME */
  62 /* #define DEBUG  */
  63 enum {
  64   ddnoSEL, ddnoINTERLEAVE, ddnoPP_PME, ddnoCARTESIAN, ddnoNR
  65 };
  66
  67 /* Enum for situations that can occur during log file parsing */
  68 enum {
  69     eParselogOK,
  70     eParselogNotFound,
  71     eParselogNoPerfData,
  72     eParselogTerm,
  73     eParselogResetProblem,
  74     eParselogNr
  75 };
  76
  77
  78 typedef struct
  79 {
  80     int  nPMEnodes;       /* number of PME only nodes used in this test */
  81     int  nx, ny, nz;      /* DD grid */
  82     int  guessPME;        /* if nPMEnodes == -1, this is the guessed number of PME nodes */
  83     float *Gcycles;        /* This can contain more than one value if doing multiple tests */
  84     float Gcycles_Av;
  85     float *ns_per_day;
  86     float ns_per_day_Av;
  87     float *PME_f_load;     /* PME mesh/force load average*/
  88     float PME_f_load_Av;   /* Average average ;) ... */
  89     char *mdrun_cmd_line; /* Mdrun command line used for this test */
  90 } t_perf;
  91
  92
  93 typedef struct
  94 {
  95     gmx_large_int_t orig_sim_steps;  /* Number of steps to be done in the real simulation  */
  96     int  n_entries;             /* Number of entries in arrays                        */
  97     real volume;                /* The volume of the box                              */
  98     matrix recipbox;            /* The reciprocal box                                 */
  99     int  natoms;                /* The number of atoms in the MD system               */
 100     real *fac;                  /* The scaling factor                                 */
 101     real *rcoulomb;             /* The coulomb radii [0...nr_inputfiles]              */
 102     real *rvdw;                 /* The vdW radii                                      */
 103     int  *nkx, *nky, *nkz;      /* Number of k vectors in each spatial dimension      */
 104     real *fourier_sp;           /* Fourierspacing                                     */
 105     real *ewald_rtol;           /* Real space tolerance for Ewald, determines         */
 106                                 /* the real/reciprocal space relative weight          */
 107     real *ewald_beta;           /* Splitting parameter [1/nm]                         */
 108     real fracself;              /* fraction of particles for SI error                 */
 109     real q2all;                 /* sum ( q ^2 )                                       */
 110     real q2allnr;               /* nr of charges                                      */
 111     int  *pme_order;            /* Interpolation order for PME (bsplines)             */
 112     char **fn_out;              /* Name of the output tpr file                        */
 113     real *e_dir;                /* Direct space part of PME error with these settings */
 114     real *e_rec;                /* Reciprocal space part of PME error                 */
 115     gmx_bool bTUNE;                 /* flag for tuning */
 116 } t_inputinfo;
 117
 118
 119 /* Returns TRUE when atom is charged */
 120 static gmx_bool is_charge(real charge)
 121 {
 122     if (charge*charge > GMX_REAL_EPS)
 123         return TRUE;
 124     else
 125         return FALSE;
 126 }
 127
 128
 129 /* calculate charge density */
 130 static void calc_q2all(
 131         gmx_mtop_t *mtop,   /* molecular topology */
 132         real *q2all, real *q2allnr)
 133 {
 134     int imol,iatom;  /* indices for loops */
 135     real q2_all=0;   /* Sum of squared charges */
 136     int  nrq_mol;    /* Number of charges in a single molecule */
 137     int  nrq_all;    /* Total number of charges in the MD system */
 138     real nrq_all_r;  /* No of charges in real format */
 139     real qi,q2_mol;
 140     gmx_moltype_t *molecule;
 141     gmx_molblock_t *molblock;
 142
 143 #ifdef DEBUG
 144         fprintf(stderr, "\nCharge density:\n");
 145 #endif
 146     q2_all = 0.0;  /* total q squared */
 147     nrq_all = 0;   /* total number of charges in the system */
 148     for (imol=0; imol<mtop->nmolblock; imol++) /* Loop over molecule types */
 149     {
 150         q2_mol=0.0; /* q squared value of this molecule */
 151         nrq_mol=0;  /* number of charges this molecule carries */
 152         molecule = &(mtop->moltype[imol]);
 153         molblock = &(mtop->molblock[imol]);
 154         for (iatom=0; iatom<molblock->natoms_mol; iatom++) /* Loop over atoms in this molecule */
 155         {
 156             qi = molecule->atoms.atom[iatom].q; /* Charge of this atom */
 157             /* Is this charge worth to be considered? */
 158             if (is_charge(qi))
 159             {
 160                 q2_mol += qi*qi;
 161                 nrq_mol++;
 162             }
 163         }
 164         /* Multiply with the number of molecules present of this type and add */
 165         q2_all  += q2_mol*molblock->nmol;
 166         nrq_all += nrq_mol*molblock->nmol;
 167 #ifdef DEBUG
 168         fprintf(stderr, "Molecule %2d (%5d atoms) q2_mol=%10.3e nr.mol.charges=%5d (%6dx)  q2_all=%10.3e  tot.charges=%d\n",
 169                 imol,molblock->natoms_mol,q2_mol,nrq_mol,molblock->nmol,q2_all,nrq_all);
 170 #endif
 171     }
 172     nrq_all_r = nrq_all;
 173
 174     *q2all=q2_all;
 175     *q2allnr=nrq_all;
 176
 177 }
 178
 179
 180 /* Estimate the direct space part error of the SPME Ewald sum */
 181 static real estimate_direct(
 182         t_inputinfo *info
 183         )
 184 {
 185     real e_dir=0;    /* Error estimate */
 186     real beta=0;     /* Splitting parameter (1/nm) */
 187     real r_coulomb=0;  /* Cut-off in direct space */
 188
 189
 190     beta      = info->ewald_beta[0];
 191     r_coulomb = info->rcoulomb[0];
 192
 193     e_dir  = 2.0 * info->q2all * gmx_invsqrt( info->q2allnr  *  r_coulomb * info->volume );
 194     e_dir *= exp (-beta*beta*r_coulomb*r_coulomb);
 195
 196     return ONE_4PI_EPS0*e_dir;
 197 }
 198
 199 #define SUMORDER 6
 200
 201 /* the following 4 functions determine polynomials required for the reciprocal error estimate */
 202
 203 static inline real eps_poly1(
 204         real m,          /* grid coordinate in certain direction */
 205         real K,          /* grid size in corresponding direction */
 206         real n)          /* spline interpolation order of the SPME */
 207 {
 208     int i;
 209     real nom=0;  /* nominator */
 210     real denom=0; /* denominator */
 211     real tmp=0;
 212
 213     if ( m == 0.0 )
 214         return 0.0 ;
 215
 216     for(i=-SUMORDER ; i<0 ; i++)
 217     {
 218         tmp=m / K + i;
 219         tmp*=2.0*M_PI;
 220         nom+=pow( tmp , -n );
 221     }
 222
 223     for(i=SUMORDER ; i>0 ; i--)
 224     {
 225         tmp=m / K + i;
 226         tmp*=2.0*M_PI;
 227         nom+=pow( tmp , -n );
 228     }
 229
 230     tmp=m / K;
 231     tmp*=2.0*M_PI;
 232     denom=pow( tmp , -n )+nom;
 233
 234     return -nom/denom;
 235
 236 }
 237
 238 static inline real eps_poly2(
 239         real m,          /* grid coordinate in certain direction */
 240         real K,          /* grid size in corresponding direction */
 241         real n)          /* spline interpolation order of the SPME */
 242 {
 243     int i;
 244     real nom=0;  /* nominator */
 245     real denom=0; /* denominator */
 246     real tmp=0;
 247
 248     if ( m == 0.0 )
 249         return 0.0 ;
 250
 251     for(i=-SUMORDER ; i<0 ; i++)
 252     {
 253         tmp=m / K + i;
 254         tmp*=2.0*M_PI;
 255         nom+=pow( tmp , -2.0*n );
 256     }
 257
 258     for(i=SUMORDER ; i>0 ; i--)
 259     {
 260         tmp=m / K + i;
 261         tmp*=2.0*M_PI;
 262         nom+=pow( tmp , -2.0*n );
 263     }
 264
 265     for(i=-SUMORDER ; i<SUMORDER+1 ; i++)
 266     {
 267         tmp=m / K + i;
 268         tmp*=2.0*M_PI;
 269         denom+=pow( tmp , -n );
 270     }
 271     tmp=eps_poly1(m,K,n);
 272     return nom / denom / denom + tmp*tmp ;
 273
 274 }
 275
 276 static inline real eps_poly3(
 277         real m,          /* grid coordinate in certain direction */
 278         real K,          /* grid size in corresponding direction */
 279         real n)          /* spline interpolation order of the SPME */
 280 {
 281     int i;
 282     real nom=0;  /* nominator */
 283     real denom=0; /* denominator */
 284     real tmp=0;
 285
 286     if ( m == 0.0 )
 287         return 0.0 ;
 288
 289     for(i=-SUMORDER ; i<0 ; i++)
 290     {
 291         tmp=m / K + i;
 292         tmp*=2.0*M_PI;
 293         nom+= i * pow( tmp , -2.0*n );
 294     }
 295
 296     for(i=SUMORDER ; i>0 ; i--)
 297     {
 298         tmp=m / K + i;
 299         tmp*=2.0*M_PI;
 300         nom+= i * pow( tmp , -2.0*n );
 301     }
 302
 303     for(i=-SUMORDER ; i<SUMORDER+1 ; i++)
 304     {
 305         tmp=m / K + i;
 306         tmp*=2.0*M_PI;
 307         denom+=pow( tmp , -n );
 308     }
 309
 310     return 2.0 * M_PI * nom / denom / denom;
 311
 312 }
 313
 314 static inline real eps_poly4(
 315         real m,          /* grid coordinate in certain direction */
 316         real K,          /* grid size in corresponding direction */
 317         real n)          /* spline interpolation order of the SPME */
 318 {
 319     int i;
 320     real nom=0;  /* nominator */
 321     real denom=0; /* denominator */
 322     real tmp=0;
 323
 324     if ( m == 0.0 )
 325         return 0.0 ;
 326
 327     for(i=-SUMORDER ; i<0 ; i++)
 328     {
 329         tmp=m / K + i;
 330         tmp*=2.0*M_PI;
 331         nom+= i * i * pow( tmp , -2.0*n );
 332     }
 333
 334     for(i=SUMORDER ; i>0 ; i--)
 335     {
 336         tmp=m / K + i;
 337         tmp*=2.0*M_PI;
 338         nom+= i * i * pow( tmp , -2.0*n );
 339     }
 340
 341     for(i=-SUMORDER ; i<SUMORDER+1 ; i++)
 342     {
 343         tmp=m / K + i;
 344         tmp*=2.0*M_PI;
 345         denom+=pow( tmp , -n );
 346     }
 347
 348     return 4.0 * M_PI * M_PI * nom / denom / denom;
 349
 350 }
 351
 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 */
 358 {
 359     int i;
 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 */
 366
 367
 368     if ( m == 0.0 )
 369         return 0.0 ;
 370
 371     rcoord=iprod(rboxv,x);
 372
 373
 374     for(i=-SUMORDER;i<0;i++)
 375     {
 376         tmp=-sin(2.0 * M_PI * i * K * rcoord);
 377         tmp1=2.0 * M_PI * m / K + 2.0 * M_PI * i;
 378         tmp2=pow(tmp1,-1.0*n);
 379         nom+=tmp * tmp2 * i;
 380         denom+=tmp2;
 381     }
 382
 383     for(i=SUMORDER;i>0;i--)
 384     {
 385         tmp=-sin(2.0 * M_PI * i * K * rcoord);
 386         tmp1=2.0 * M_PI * m / K + 2.0 * M_PI * i;
 387         tmp2=pow(tmp1,-1.0*n);
 388         nom+=tmp * tmp2 * i;
 389         denom+=tmp2;
 390     }
 391
 392
 393     tmp=2.0 * M_PI * m / K;
 394     tmp1=pow(tmp,-1.0*n);
 395     denom+=tmp1;
 396
 397    return 2.0 * M_PI * nom / denom * K  ;
 398
 399 }
 400
 401 #undef SUMORDER
 402
 403 /* The following routine is just a copy from pme.c */
 404
 405 static void calc_recipbox(matrix box,matrix recipbox)
 406 {
 407   /* Save some time by assuming upper right part is zero */
 408
 409   real tmp=1.0/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
 410
 411   recipbox[XX][XX]=box[YY][YY]*box[ZZ][ZZ]*tmp;
 412   recipbox[XX][YY]=0;
 413   recipbox[XX][ZZ]=0;
 414   recipbox[YY][XX]=-box[YY][XX]*box[ZZ][ZZ]*tmp;
 415   recipbox[YY][YY]=box[XX][XX]*box[ZZ][ZZ]*tmp;
 416   recipbox[YY][ZZ]=0;
 417   recipbox[ZZ][XX]=(box[YY][XX]*box[ZZ][YY]-box[YY][YY]*box[ZZ][XX])*tmp;
 418   recipbox[ZZ][YY]=-box[ZZ][YY]*box[XX][XX]*tmp;
 419   recipbox[ZZ][ZZ]=box[XX][XX]*box[YY][YY]*tmp;
 420 }
 421
 422
 423 /* Estimate the reciprocal space part error of the SPME Ewald sum. */
 424 static real estimate_reciprocal(
 425         t_inputinfo *info,
 426         rvec x[],           /* array of particles */
 427         real q[],           /* array of charges */
 428         int nr,             /* number of charges = size of the charge array */
 429         FILE *fp_out,
 430         gmx_bool bVerbose,
 431         unsigned int seed,  /* The seed for the random number generator */
 432         int *nsamples,      /* Return the number of samples used if Monte Carlo
 433                              * algorithm is used for self energy error estimate */
 434         t_commrec *cr)
 435 {
 436     real e_rec=0;   /* reciprocal error estimate */
 437     real e_rec1=0;  /* Error estimate term 1*/
 438     real e_rec2=0;  /* Error estimate term 2*/
 439     real e_rec3=0;  /* Error estimate term 3 */
 440     real e_rec3x=0; /* part of Error estimate term 3 in x */
 441     real e_rec3y=0; /* part of Error estimate term 3 in y */
 442     real e_rec3z=0; /* part of Error estimate term 3 in z */
 443     int i,ci;
 444     int nx,ny,nz;   /* grid coordinates */
 445     real q2_all=0;  /* sum of squared charges */
 446     rvec gridpx;    /* reciprocal grid point in x direction*/
 447     rvec gridpxy;   /* reciprocal grid point in x and y direction*/
 448     rvec gridp;     /* complete reciprocal grid point in 3 directions*/
 449     rvec tmpvec;    /* template to create points from basis vectors */
 450     rvec tmpvec2;   /* template to create points from basis vectors */
 451     real coeff=0;   /* variable to compute coefficients of the error estimate */
 452     real coeff2=0;   /* variable to compute coefficients of the error estimate */
 453     real tmp=0;     /* variables to compute different factors from vectors */
 454     real tmp1=0;
 455     real tmp2=0;
 456     gmx_bool bFraction;
 457
 458     /* Random number generator */
 459     gmx_rng_t rng=NULL;
 460     int *numbers=NULL;
 461
 462     /* Index variables for parallel work distribution */
 463     int startglobal,stopglobal;
 464     int startlocal, stoplocal;
 465     int x_per_core;
 466     int xtot;
 467
 468 #ifdef TAKETIME
 469     double t0=0.0;
 470     double t1=0.0;
 471 #endif
 472
 473     rng=gmx_rng_init(seed);
 474
 475     clear_rvec(gridpx);
 476     clear_rvec(gridpxy);
 477     clear_rvec(gridp);
 478     clear_rvec(tmpvec);
 479     clear_rvec(tmpvec2);
 480
 481     for(i=0;i<nr;i++)
 482     {
 483         q2_all += q[i]*q[i];
 484     }
 485
 486     /* Calculate indices for work distribution */
 487     startglobal=-info->nkx[0]/2;
 488     stopglobal = info->nkx[0]/2;
 489     xtot = stopglobal*2+1;
 490     if (PAR(cr))
 491     {
 492         x_per_core = ceil((real)xtot / (real)cr->nnodes);
 493         startlocal = startglobal + x_per_core*cr->nodeid;
 494         stoplocal = startlocal + x_per_core -1;
 495         if (stoplocal > stopglobal)
 496              stoplocal = stopglobal;
 497     }
 498     else
 499     {
 500         startlocal = startglobal;
 501         stoplocal  = stopglobal;
 502         x_per_core = xtot;
 503     }
 504 /*
 505 #ifdef GMX_LIB_MPI
 506     MPI_Barrier(MPI_COMM_WORLD);
 507 #endif
 508 */
 509
 510 #ifdef GMX_LIB_MPI
 511 #ifdef TAKETIME
 512     if (MASTER(cr))
 513         t0 = MPI_Wtime();
 514 #endif
 515 #endif
 516
 517     if (MASTER(cr)){
 518
 519         fprintf(stderr, "Calculating reciprocal error part 1 ...");
 520
 521     }
 522
 523     for(nx=startlocal; nx<=stoplocal; nx++)
 524     {
 525         svmul(nx,info->recipbox[XX],gridpx);
 526         for(ny=-info->nky[0]/2; ny<info->nky[0]/2+1; ny++)
 527         {
 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++)
 531             {
 532                 if (  0 == nx &&  0 == ny &&  0 == nz )
 533                     continue;
 534                 svmul(nz,info->recipbox[ZZ],tmpvec);
 535                 rvec_add(gridpxy,tmpvec,gridp);
 536                 tmp=norm2(gridp);
 537                 coeff=exp(-1.0 * M_PI * M_PI * tmp / info->ewald_beta[0] / info->ewald_beta[0] ) ;
 538                 coeff/= 2.0 * M_PI * info->volume * tmp;
 539                 coeff2=tmp ;
 540
 541
 542                 tmp=eps_poly2(nx,info->nkx[0],info->pme_order[0]);
 543                 tmp+=eps_poly2(ny,info->nkx[0],info->pme_order[0]);
 544                 tmp+=eps_poly2(nz,info->nkx[0],info->pme_order[0]);
 545
 546                 tmp1=eps_poly1(nx,info->nkx[0],info->pme_order[0]);
 547                 tmp2=eps_poly1(ny,info->nky[0],info->pme_order[0]);
 548
 549                 tmp+=2.0 * tmp1 * tmp2;
 550
 551                 tmp1=eps_poly1(nz,info->nkz[0],info->pme_order[0]);
 552                 tmp2=eps_poly1(ny,info->nky[0],info->pme_order[0]);
 553
 554                 tmp+=2.0 * tmp1 * tmp2;
 555
 556                 tmp1=eps_poly1(nz,info->nkz[0],info->pme_order[0]);
 557                 tmp2=eps_poly1(nx,info->nkx[0],info->pme_order[0]);
 558
 559                 tmp+=2.0 * tmp1 * tmp2;
 560
 561                 tmp1=eps_poly1(nx,info->nkx[0],info->pme_order[0]);
 562                 tmp1+=eps_poly1(ny,info->nky[0],info->pme_order[0]);
 563                 tmp1+=eps_poly1(nz,info->nkz[0],info->pme_order[0]);
 564
 565                 tmp+= tmp1 * tmp1;
 566
 567                 e_rec1+= 32.0 * M_PI * M_PI * coeff * coeff * coeff2 * tmp  * q2_all * q2_all / nr ;
 568
 569                 tmp1=eps_poly3(nx,info->nkx[0],info->pme_order[0]);
 570                 tmp1*=info->nkx[0];
 571                 tmp2=iprod(gridp,info->recipbox[XX]);
 572
 573                 tmp=tmp1*tmp2;
 574
 575                 tmp1=eps_poly3(ny,info->nky[0],info->pme_order[0]);
 576                 tmp1*=info->nky[0];
 577                 tmp2=iprod(gridp,info->recipbox[YY]);
 578
 579                 tmp+=tmp1*tmp2;
 580
 581                 tmp1=eps_poly3(nz,info->nkz[0],info->pme_order[0]);
 582                 tmp1*=info->nkz[0];
 583                 tmp2=iprod(gridp,info->recipbox[ZZ]);
 584
 585                 tmp+=tmp1*tmp2;
 586
 587                 tmp*=4.0 * M_PI;
 588
 589                 tmp1=eps_poly4(nx,info->nkx[0],info->pme_order[0]);
 590                 tmp1*=norm2(info->recipbox[XX]);
 591                 tmp1*=info->nkx[0] * info->nkx[0];
 592
 593                 tmp+=tmp1;
 594
 595                 tmp1=eps_poly4(ny,info->nky[0],info->pme_order[0]);
 596                 tmp1*=norm2(info->recipbox[YY]);
 597                 tmp1*=info->nky[0] * info->nky[0];
 598
 599                 tmp+=tmp1;
 600
 601                 tmp1=eps_poly4(nz,info->nkz[0],info->pme_order[0]);
 602                 tmp1*=norm2(info->recipbox[ZZ]);
 603                 tmp1*=info->nkz[0] * info->nkz[0];
 604
 605                 tmp+=tmp1;
 606
 607                 e_rec2+= 4.0 * coeff * coeff * tmp * q2_all * q2_all / nr ;
 608
 609             }
 610         }
 611         if (MASTER(cr))
 612             fprintf(stderr, "\rCalculating reciprocal error part 1 ... %3.0f%%", 100.0*(nx-startlocal+1)/(x_per_core));
 613
 614     }
 615
 616     if (MASTER(cr))
 617         fprintf(stderr, "\n");
 618
 619     /* Use just a fraction of all charges to estimate the self energy error term? */
 620     bFraction =  (info->fracself > 0.0) && (info->fracself < 1.0);
 621
 622     if (bFraction)
 623     {
 624         /* Here xtot is the number of samples taken for the Monte Carlo calculation
 625          * of the average of term IV of equation 35 in Wang2010. Round up to a
 626          * number of samples that is divisible by the number of nodes */
 627         x_per_core  = ceil(info->fracself * nr / (real)cr->nnodes);
 628         xtot = x_per_core * cr->nnodes;
 629     }
 630     else
 631     {
 632         /* In this case we use all nr particle positions */
 633         xtot = nr;
 634         x_per_core = ceil( (real)xtot / (real)cr->nnodes );
 635     }
 636
 637     startlocal = x_per_core *  cr->nodeid;
 638     stoplocal  = min(startlocal + x_per_core, xtot);  /* min needed if xtot == nr */
 639
 640     if (bFraction)
 641     {
 642         /* Make shure we get identical results in serial and parallel. Therefore,
 643          * take the sample indices from a single, global random number array that
 644          * is constructed on the master node and that only depends on the seed */
 645         snew(numbers, xtot);
 646         if (MASTER(cr))
 647         {
 648             for (i=0; i<xtot; i++)
 649             {
 650                 numbers[i] = floor(gmx_rng_uniform_real(rng) * nr );
 651             }
 652         }
 653         /* Broadcast the random number array to the other nodes */
 654         if (PAR(cr))
 655         {
 656             nblock_bc(cr,xtot,numbers);
 657         }
 658
 659         if (bVerbose && MASTER(cr))
 660         {
 661             fprintf(stdout, "Using %d sample%s to approximate the self interaction error term",
 662                     xtot, xtot==1?"":"s");
 663             if (PAR(cr))
 664                 fprintf(stdout, " (%d sample%s per node)", x_per_core, x_per_core==1?"":"s");
 665             fprintf(stdout, ".\n");
 666         }
 667     }
 668
 669     /* Return the number of positions used for the Monte Carlo algorithm */
 670     *nsamples = xtot;
 671
 672     for(i=startlocal;i<stoplocal;i++)
 673     {
 674         e_rec3x=0;
 675         e_rec3y=0;
 676         e_rec3z=0;
 677
 678         if (bFraction)
 679         {
 680             /* Randomly pick a charge */
 681             ci = numbers[i];
 682         }
 683         else
 684         {
 685             /* Use all charges */
 686             ci = i;
 687         }
 688
 689         /* for(nx=startlocal; nx<=stoplocal; nx++)*/
 690         for(nx=-info->nkx[0]/2; nx<info->nkx[0]/2+1; nx++)
 691         {
 692             svmul(nx,info->recipbox[XX],gridpx);
 693             for(ny=-info->nky[0]/2; ny<info->nky[0]/2+1; ny++)
 694             {
 695                 svmul(ny,info->recipbox[YY],tmpvec);
 696                 rvec_add(gridpx,tmpvec,gridpxy);
 697                 for(nz=-info->nkz[0]/2; nz<info->nkz[0]/2+1; nz++)
 698                 {
 699
 700                     if (  0 == nx && 0 == ny && 0 == nz)
 701                         continue;
 702
 703                     svmul(nz,info->recipbox[ZZ],tmpvec);
 704                     rvec_add(gridpxy,tmpvec,gridp);
 705                     tmp=norm2(gridp);
 706                     coeff=exp(-1.0 * M_PI * M_PI * tmp / info->ewald_beta[0] / info->ewald_beta[0] );
 707                     coeff/= tmp ;
 708                     e_rec3x+=coeff*eps_self(nx,info->nkx[0],info->recipbox[XX],info->pme_order[0],x[ci]);
 709                     e_rec3y+=coeff*eps_self(ny,info->nky[0],info->recipbox[YY],info->pme_order[0],x[ci]);
 710                     e_rec3z+=coeff*eps_self(nz,info->nkz[0],info->recipbox[ZZ],info->pme_order[0],x[ci]);
 711
 712                 }
 713             }
 714         }
 715
 716         clear_rvec(tmpvec2);
 717
 718         svmul(e_rec3x,info->recipbox[XX],tmpvec);
 719         rvec_inc(tmpvec2,tmpvec);
 720         svmul(e_rec3y,info->recipbox[YY],tmpvec);
 721         rvec_inc(tmpvec2,tmpvec);
 722         svmul(e_rec3z,info->recipbox[ZZ],tmpvec);
 723         rvec_inc(tmpvec2,tmpvec);
 724
 725         e_rec3 += q[ci]*q[ci]*q[ci]*q[ci]*norm2(tmpvec2) / ( xtot * M_PI * info->volume * M_PI * info->volume);
 726         if (MASTER(cr)){
 727             fprintf(stderr, "\rCalculating reciprocal error part 2 ... %3.0f%%",
 728                     100.0*(i+1)/stoplocal);
 729
 730         }
 731     }
 732
 733     if (MASTER(cr))
 734         fprintf(stderr, "\n");
 735
 736 #ifdef GMX_LIB_MPI
 737 #ifdef TAKETIME
 738     if (MASTER(cr))
 739     {
 740         t1= MPI_Wtime() - t0;
 741         fprintf(fp_out, "Recip. err. est. took   : %lf s\n", t1);
 742     }
 743 #endif
 744 #endif
 745
 746 #ifdef DEBUG
 747     if (PAR(cr))
 748     {
 749         fprintf(stderr, "Node %3d: nx=[%3d...%3d]  e_rec3=%e\n",
 750                 cr->nodeid, startlocal, stoplocal, e_rec3);
 751     }
 752 #endif
 753
 754     if (PAR(cr))
 755     {
 756         gmx_sum(1,&e_rec1,cr);
 757         gmx_sum(1,&e_rec2,cr);
 758         gmx_sum(1,&e_rec3,cr);
 759     }
 760
 761     /* e_rec1*=8.0 * q2_all / info->volume / info->volume / nr ;
 762        e_rec2*=  q2_all / M_PI / M_PI / info->volume / info->volume / nr ;
 763        e_rec3/= M_PI * M_PI * info->volume * info->volume * nr ;
 764      */
 765     e_rec=sqrt(e_rec1+e_rec2+e_rec3);
 766
 767
 768     return ONE_4PI_EPS0 * e_rec;
 769 }
 770
 771
 772 /* Allocate memory for the inputinfo struct: */
 773 static void create_info(t_inputinfo *info)
 774 {
 775     snew(info->fac       , info->n_entries);
 776     snew(info->rcoulomb  , info->n_entries);
 777     snew(info->rvdw      , info->n_entries);
 778     snew(info->nkx       , info->n_entries);
 779     snew(info->nky       , info->n_entries);
 780     snew(info->nkz       , info->n_entries);
 781     snew(info->fourier_sp, info->n_entries);
 782     snew(info->ewald_rtol, info->n_entries);
 783     snew(info->ewald_beta, info->n_entries);
 784     snew(info->pme_order , info->n_entries);
 785     snew(info->fn_out    , info->n_entries);
 786     snew(info->e_dir     , info->n_entries);
 787     snew(info->e_rec     , info->n_entries);
 788 }
 789
 790
 791 /* Allocate and fill an array with coordinates and charges,
 792  * returns the number of charges found
 793  */
 794 static int prepare_x_q(real *q[], rvec *x[], gmx_mtop_t *mtop, rvec x_orig[], t_commrec *cr)
 795 {
 796     int i;
 797     int nq; /* number of charged particles */
 798     gmx_mtop_atomloop_all_t aloop;
 799     t_atom *atom;
 800
 801
 802     if (MASTER(cr))
 803     {
 804         snew(*q, mtop->natoms);
 805         snew(*x, mtop->natoms);
 806         nq=0;
 807
 808         aloop = gmx_mtop_atomloop_all_init(mtop);
 809
 810         while (gmx_mtop_atomloop_all_next(aloop,&i,&atom))
 811         {
 812             if (is_charge(atom->q))
 813             {
 814                 (*q)[nq] = atom->q;
 815                 (*x)[nq][XX] = x_orig[i][XX];
 816                 (*x)[nq][YY] = x_orig[i][YY];
 817                 (*x)[nq][ZZ] = x_orig[i][ZZ];
 818                 nq++;
 819             }
 820         }
 821         /* Give back some unneeded memory */
 822         srenew(*q, nq);
 823         srenew(*x, nq);
 824     }
 825     /* Broadcast x and q in the parallel case */
 826     if (PAR(cr))
 827     {
 828         /* Transfer the number of charges */
 829         block_bc(cr,nq);
 830         snew_bc(cr, *x, nq);
 831         snew_bc(cr, *q, nq);
 832         nblock_bc(cr,nq,*x);
 833         nblock_bc(cr,nq,*q);
 834     }
 835
 836     return nq;
 837 }
 838
 839
 840
 841 /* Read in the tpr file and save information we need later in info */
 842 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)
 843 {
 844     read_tpx_state(fn_sim_tpr,ir,state,NULL,mtop);
 845
 846     /* The values of the original tpr input file are save in the first
 847      * place [0] of the arrays */
 848     info->orig_sim_steps = ir->nsteps;
 849     info->pme_order[0]   = ir->pme_order;
 850     info->rcoulomb[0]    = ir->rcoulomb;
 851     info->rvdw[0]        = ir->rvdw;
 852     info->nkx[0]         = ir->nkx;
 853     info->nky[0]         = ir->nky;
 854     info->nkz[0]         = ir->nkz;
 855     info->ewald_rtol[0]  = ir->ewald_rtol;
 856     info->fracself       = fracself;
 857     if (user_beta > 0)
 858         info->ewald_beta[0] = user_beta;
 859     else
 860         info->ewald_beta[0]  = calc_ewaldcoeff(info->rcoulomb[0],info->ewald_rtol[0]);
 861
 862     /* Check if PME was chosen */
 863     if (EEL_PME(ir->coulombtype) == FALSE)
 864         gmx_fatal(FARGS, "Can only do optimizations for simulations with PME");
 865
 866     /* Check if rcoulomb == rlist, which is necessary for PME */
 867     if (!(ir->rcoulomb == ir->rlist))
 868         gmx_fatal(FARGS, "PME requires rcoulomb (%f) to be equal to rlist (%f).", ir->rcoulomb, ir->rlist);
 869 }
 870
 871
 872 /* Transfer what we need for parallelizing the reciprocal error estimate */
 873 static void bcast_info(t_inputinfo *info, t_commrec *cr)
 874 {
 875     nblock_bc(cr, info->n_entries, info->nkx);
 876     nblock_bc(cr, info->n_entries, info->nky);
 877     nblock_bc(cr, info->n_entries, info->nkz);
 878     nblock_bc(cr, info->n_entries, info->ewald_beta);
 879     nblock_bc(cr, info->n_entries, info->pme_order);
 880     nblock_bc(cr, info->n_entries, info->e_dir);
 881     nblock_bc(cr, info->n_entries, info->e_rec);
 882     block_bc(cr, info->volume);
 883     block_bc(cr, info->recipbox);
 884     block_bc(cr, info->natoms);
 885     block_bc(cr, info->fracself);
 886     block_bc(cr, info->bTUNE);
 887     block_bc(cr, info->q2all);
 888     block_bc(cr, info->q2allnr);
 889 }
 890
 891
 892 /* Estimate the error of the SPME Ewald sum. This estimate is based upon
 893  * a) a homogeneous distribution of the charges
 894  * b) a total charge of zero.
 895  */
 896 static void estimate_PME_error(t_inputinfo *info, t_state *state,
 897         gmx_mtop_t *mtop, FILE *fp_out, gmx_bool bVerbose, unsigned int seed,
 898         t_commrec *cr)
 899 {
 900     rvec *x=NULL; /* The coordinates */
 901     real *q=NULL; /* The charges     */
 902     real  edir=0.0; /* real space error */
 903     real  erec=0.0; /* reciprocal space error */
 904     real  derr=0.0; /* difference of real and reciprocal space error */
 905     real  derr0=0.0; /* difference of real and reciprocal space error */
 906     real  beta=0.0; /* splitting parameter beta */
 907     real  beta0=0.0; /* splitting parameter beta */
 908     int ncharges; /* The number of atoms with charges */
 909     int nsamples; /* The number of samples used for the calculation of the
 910                    * self-energy error term */
 911     int i=0;
 912
 913     if (MASTER(cr))
 914         fprintf(fp_out, "\n--- PME ERROR ESTIMATE ---\n");
 915
 916     /* Prepare an x and q array with only the charged atoms */
 917     ncharges = prepare_x_q(&q, &x, mtop, state->x, cr);
 918     if (MASTER(cr))
 919     {
 920         calc_q2all(mtop, &(info->q2all), &(info->q2allnr));
 921         info->ewald_rtol[0]=gmx_erfc(info->rcoulomb[0]*info->ewald_beta[0]);
 922         /* Write some info to log file */
 923         fprintf(fp_out, "Box volume              : %g nm^3\n", info->volume);
 924         fprintf(fp_out, "Number of charged atoms : %d (total atoms %d)\n",ncharges, info->natoms);
 925         fprintf(fp_out, "Coulomb radius          : %g nm\n", info->rcoulomb[0]);
 926         fprintf(fp_out, "Ewald_rtol              : %g\n", info->ewald_rtol[0]);
 927         fprintf(fp_out, "Ewald parameter beta    : %g\n", info->ewald_beta[0]);
 928         fprintf(fp_out, "Interpolation order     : %d\n", info->pme_order[0]);
 929         fprintf(fp_out, "Fourier grid (nx,ny,nz) : %d x %d x %d\n",
 930                 info->nkx[0],info->nky[0],info->nkz[0]);
 931         fflush(fp_out);
 932
 933     }
 934
 935     if (PAR(cr))
 936         bcast_info(info, cr);
 937
 938
 939     /* Calculate direct space error */
 940     info->e_dir[0] = estimate_direct(info);
 941
 942     /* Calculate reciprocal space error */
 943     info->e_rec[0] = estimate_reciprocal(info, x, q, ncharges, fp_out, bVerbose,
 944                                          seed, &nsamples, cr);
 945
 946     if (PAR(cr))
 947         bcast_info(info, cr);
 948
 949     if (MASTER(cr))
 950     {
 951         fprintf(fp_out, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
 952         fprintf(fp_out, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
 953         fprintf(fp_out, "Self-energy error term was estimated using %d samples\n", nsamples);
 954         fflush(fp_out);
 955         fprintf(stderr, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
 956         fprintf(stderr, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
 957     }
 958
 959     i=0;
 960
 961     if (info->bTUNE)
 962     {
 963         if(MASTER(cr))
 964             fprintf(stderr,"Starting tuning ...\n");
 965         edir=info->e_dir[0];
 966         erec=info->e_rec[0];
 967         derr0=edir-erec;
 968         beta0=info->ewald_beta[0];
 969         if (derr>0.0)
 970             info->ewald_beta[0]+=0.1;
 971         else
 972             info->ewald_beta[0]-=0.1;
 973         info->e_dir[0] = estimate_direct(info);
 974         info->e_rec[0] = estimate_reciprocal(info, x, q, ncharges, fp_out, bVerbose,
 975                                              seed, &nsamples, cr);
 976
 977         if (PAR(cr))
 978             bcast_info(info, cr);
 979
 980
 981         edir=info->e_dir[0];
 982         erec=info->e_rec[0];
 983         derr=edir-erec;
 984         while ( fabs(derr/min(erec,edir)) > 1e-4)
 985         {
 986
 987             beta=info->ewald_beta[0];
 988             beta-=derr*(info->ewald_beta[0]-beta0)/(derr-derr0);
 989             beta0=info->ewald_beta[0];
 990             info->ewald_beta[0]=beta;
 991             derr0=derr;
 992
 993             info->e_dir[0] = estimate_direct(info);
 994             info->e_rec[0] = estimate_reciprocal(info, x, q, ncharges, fp_out, bVerbose,
 995                                                  seed, &nsamples, cr);
 996
 997             if (PAR(cr))
 998                 bcast_info(info, cr);
 999
1000             edir=info->e_dir[0];
1001             erec=info->e_rec[0];
1002             derr=edir-erec;
1003
1004             if (MASTER(cr))
1005             {
1006                 i++;
1007                 fprintf(stderr,"difference between real and rec. space error (step %d): %g\n",i,fabs(derr));
1008                 fprintf(stderr,"old beta: %f\n",beta0);
1009                 fprintf(stderr,"new beta: %f\n",beta);
1010             }
1011         }
1012
1013         info->ewald_rtol[0]=gmx_erfc(info->rcoulomb[0]*info->ewald_beta[0]);
1014
1015         if (MASTER(cr))
1016         {
1017             /* Write some info to log file */
1018             fflush(fp_out);
1019             fprintf(fp_out, "=========  After tuning ========\n");
1020             fprintf(fp_out, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
1021             fprintf(fp_out, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
1022             fprintf(stderr, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
1023             fprintf(stderr, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
1024             fprintf(fp_out, "Ewald_rtol              : %g\n", info->ewald_rtol[0]);
1025             fprintf(fp_out, "Ewald parameter beta    : %g\n", info->ewald_beta[0]);
1026             fflush(fp_out);
1027
1028         }
1029
1030     }
1031
1032 }
1033
1034
1035 int gmx_pme_error(int argc,char *argv[])
1036 {
1037     const char *desc[] = {
1038             "[TT]g_pme_error[tt] estimates the error of the electrostatic forces",
1039             "if using the sPME algorithm. The flag [TT]-tune[tt] will determine",
1040             "the splitting parameter such that the error is equally",
1041             "distributed over the real and reciprocal space part.",
1042             "The part of the error that stems from self interaction of the particles "
1043             "is computationally demanding. However, a good a approximation is to",
1044             "just use a fraction of the particles for this term which can be",
1045             "indicated by the flag [TT]-self[tt].[PAR]",
1046     };
1047
1048     real        fs=0.0;             /* 0 indicates: not set by the user */
1049     real        user_beta=-1.0;
1050     real        fracself=1.0;
1051     t_inputinfo info;
1052     t_state     state;     /* The state from the tpr input file */
1053     gmx_mtop_t  mtop;      /* The topology from the tpr input file */
1054     t_inputrec  *ir=NULL;  /* The inputrec from the tpr file */
1055     FILE        *fp=NULL;
1056     t_commrec   *cr;
1057     unsigned long PCA_Flags;
1058     gmx_bool        bTUNE=FALSE;
1059     gmx_bool    bVerbose=FALSE;
1060     int         seed=0;
1061
1062
1063     static t_filenm fnm[] = {
1064       { efTPX, "-s",     NULL,    ffREAD },
1065       { efOUT, "-o",    "error",  ffWRITE },
1066       { efTPX, "-so",   "tuned",  ffOPTWR }
1067     };
1068
1069     output_env_t oenv=NULL;
1070
1071     t_pargs pa[] = {
1072         { "-beta",     FALSE, etREAL, {&user_beta},
1073             "If positive, overwrite ewald_beta from [TT].tpr[tt] file with this value" },
1074         { "-tune",     FALSE, etBOOL, {&bTUNE},
1075             "Tune the splitting parameter such that the error is equally distributed between real and reciprocal space" },
1076         { "-self",     FALSE, etREAL, {&fracself},
1077             "If between 0.0 and 1.0, determine self interaction error from just this fraction of the charged particles" },
1078         { "-seed",     FALSE, etINT,  {&seed},
1079           "Random number seed used for Monte Carlo algorithm when [TT]-self[tt] is set to a value between 0.0 and 1.0" },
1080         { "-v",        FALSE, etBOOL, {&bVerbose},
1081             "Be loud and noisy" }
1082     };
1083
1084
1085 #define NFILE asize(fnm)
1086
1087     cr = init_par(&argc,&argv);
1088
1089 #ifdef GMX_LIB_MPI
1090     MPI_Barrier(MPI_COMM_WORLD);
1091 #endif
1092
1093     if (MASTER(cr))
1094       CopyRight(stderr,argv[0]);
1095
1096     PCA_Flags = PCA_NOEXIT_ON_ARGS;
1097     PCA_Flags |= (MASTER(cr) ? 0 : PCA_QUIET);
1098
1099     parse_common_args(&argc,argv,PCA_Flags,
1100                       NFILE,fnm,asize(pa),pa,asize(desc),desc,
1101                       0,NULL,&oenv);
1102
1103     if (!bTUNE)
1104         bTUNE = opt2bSet("-so",NFILE,fnm);
1105
1106     info.n_entries = 1;
1107
1108     /* Allocate memory for the inputinfo struct: */
1109     create_info(&info);
1110     info.fourier_sp[0] = fs;
1111
1112     /* Read in the tpr file and open logfile for reading */
1113     if (MASTER(cr))
1114     {
1115         snew(ir,1);
1116         read_tpr_file(opt2fn("-s",NFILE,fnm), &info, &state, &mtop, ir, user_beta,fracself);
1117
1118         fp=fopen(opt2fn("-o",NFILE,fnm),"w");
1119     }
1120
1121     /* Check consistency if the user provided fourierspacing */
1122     if (fs > 0 && MASTER(cr))
1123     {
1124         /* Recalculate the grid dimensions using fourierspacing from user input */
1125         info.nkx[0] = 0;
1126         info.nky[0] = 0;
1127         info.nkz[0] = 0;
1128         calc_grid(stdout,state.box,info.fourier_sp[0],&(info.nkx[0]),&(info.nky[0]),&(info.nkz[0]));
1129         if ( (ir->nkx != info.nkx[0]) || (ir->nky != info.nky[0]) || (ir->nkz != info.nkz[0]) )
1130             gmx_fatal(FARGS, "Wrong fourierspacing %f nm, input file grid = %d x %d x %d, computed grid = %d x %d x %d",
1131                       fs,ir->nkx,ir->nky,ir->nkz,info.nkx[0],info.nky[0],info.nkz[0]);
1132     }
1133
1134     /* Estimate (S)PME force error */
1135
1136     /* Determine the volume of the simulation box */
1137     if (MASTER(cr))
1138     {
1139         info.volume = det(state.box);
1140         calc_recipbox(state.box,info.recipbox);
1141         info.natoms = mtop.natoms;
1142         info.bTUNE  = bTUNE;
1143     }
1144
1145     if (PAR(cr))
1146         bcast_info(&info, cr);
1147
1148     /* Get an error estimate of the input tpr file and do some tuning if requested */
1149     estimate_PME_error(&info, &state, &mtop, fp, bVerbose, seed, cr);
1150
1151     if (MASTER(cr))
1152     {
1153         /* Write out optimized tpr file if requested */
1154         if ( opt2bSet("-so",NFILE,fnm) || bTUNE )
1155         {
1156             ir->ewald_rtol=info.ewald_rtol[0];
1157             write_tpx_state(opt2fn("-so",NFILE,fnm),ir,&state,&mtop);
1158         }
1159         please_cite(fp,"Wang2010");
1160         fclose(fp);
1161     }
1162
1163     gmx_finalize_par();
1164
1165     return 0;
1166 }