src/tools/gmx_pme_error.c

   1 /*
   2  *                This source code is part of
   3  *
   4  *                 G   R   O   M   A   C   S
   5  *
   6  *          GROningen MAchine for Chemical Simulations
   7  *
   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.
  12
  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.
  17  *
  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.
  24  *
  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.
  27  *
  28  * For more info, check our website at http://www.gromacs.org
  29  *
  30  * And Hey:
  31  * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
  32  */
  33 #include "statutil.h"
  34 #include "typedefs.h"
  35 #include "smalloc.h"
  36 #include "vec.h"
  37 #include "copyrite.h"
  38 #include "tpxio.h"
  39 #include "string2.h"
  40 #include "readinp.h"
  41 #include "calcgrid.h"
  42 #include "checkpoint.h"
  43 #include "gmx_ana.h"
  44 #include "gmx_random.h"
  45 #include "physics.h"
  46 #include "mdatoms.h"
  47 #include "coulomb.h"
  48 #include "mtop_util.h"
  49 #include "network.h"
  50 #include "main.h"
  51
  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 */
  57 /* #define DEBUG  */
  58 enum {
  59   ddnoSEL, ddnoINTERLEAVE, ddnoPP_PME, ddnoCARTESIAN, ddnoNR
  60 };
  61
  62 /* Enum for situations that can occur during log file parsing */
  63 enum {
  64     eParselogOK,
  65     eParselogNotFound,
  66     eParselogNoPerfData,
  67     eParselogTerm,
  68     eParselogResetProblem,
  69     eParselogNr
  70 };
  71
  72
  73 typedef struct
  74 {
  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 */
  79     float Gcycles_Av;
  80     float *ns_per_day;
  81     float ns_per_day_Av;
  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 */
  85 } t_perf;
  86
  87
  88 typedef struct
  89 {
  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 */
 111 } t_inputinfo;
 112
 113
 114 /* Returns TRUE when atom is charged */
 115 static gmx_bool is_charge(real charge)
 116 {
 117     if (charge*charge > GMX_REAL_EPS)
 118         return TRUE;
 119     else
 120         return FALSE;
 121 }
 122
 123
 124 /* calculate charge density */
 125 static void calc_q2all(
 126         gmx_mtop_t *mtop,   /* molecular topology */
 127         real *q2all, real *q2allnr)
 128 {
 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 */
 134     real qi,q2_mol;
 135     gmx_moltype_t *molecule;
 136     gmx_molblock_t *molblock;
 137
 138 #ifdef DEBUG
 139         fprintf(stderr, "\nCharge density:\n");
 140 #endif
 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 */
 144     {
 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 */
 150         {
 151             qi = molecule->atoms.atom[iatom].q; /* Charge of this atom */
 152             /* Is this charge worth to be considered? */
 153             if (is_charge(qi))
 154             {
 155                 q2_mol += qi*qi;
 156                 nrq_mol++;
 157             }
 158         }
 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;
 162 #ifdef DEBUG
 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);
 165 #endif
 166     }
 167     nrq_all_r = nrq_all;
 168
 169     *q2all=q2_all;
 170     *q2allnr=nrq_all;
 171
 172 }
 173
 174
 175 /* Estimate the direct space part error of the SPME Ewald sum */
 176 static real estimate_direct(
 177         t_inputinfo *info
 178         )
 179 {
 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 */
 183
 184
 185     beta      = info->ewald_beta[0];
 186     r_coulomb = info->rcoulomb[0];
 187
 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);
 190
 191     return ONE_4PI_EPS0*e_dir;
 192 }
 193
 194 #define SUMORDER 6
 195
 196 /* the following 4 functions determine polynomials required for the reciprocal error estimate */
 197
 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 */
 202 {
 203     int i;
 204     real nom=0;  /* nominator */
 205     real denom=0; /* denominator */
 206     real tmp=0;
 207
 208     if ( m == 0.0 )
 209         return 0.0 ;
 210
 211     for(i=-SUMORDER ; i<0 ; i++)
 212     {
 213         tmp=m / K + i;
 214         tmp*=2.0*M_PI;
 215         nom+=pow( tmp , -n );
 216     }
 217
 218     for(i=SUMORDER ; i>0 ; i--)
 219     {
 220         tmp=m / K + i;
 221         tmp*=2.0*M_PI;
 222         nom+=pow( tmp , -n );
 223     }
 224
 225     tmp=m / K;
 226     tmp*=2.0*M_PI;
 227     denom=pow( tmp , -n )+nom;
 228
 229     return -nom/denom;
 230
 231 }
 232
 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 */
 237 {
 238     int i;
 239     real nom=0;  /* nominator */
 240     real denom=0; /* denominator */
 241     real tmp=0;
 242
 243     if ( m == 0.0 )
 244         return 0.0 ;
 245
 246     for(i=-SUMORDER ; i<0 ; i++)
 247     {
 248         tmp=m / K + i;
 249         tmp*=2.0*M_PI;
 250         nom+=pow( tmp , -2.0*n );
 251     }
 252
 253     for(i=SUMORDER ; i>0 ; i--)
 254     {
 255         tmp=m / K + i;
 256         tmp*=2.0*M_PI;
 257         nom+=pow( tmp , -2.0*n );
 258     }
 259
 260     for(i=-SUMORDER ; i<SUMORDER+1 ; i++)
 261     {
 262         tmp=m / K + i;
 263         tmp*=2.0*M_PI;
 264         denom+=pow( tmp , -n );
 265     }
 266     tmp=eps_poly1(m,K,n);
 267     return nom / denom / denom + tmp*tmp ;
 268
 269 }
 270
 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 */
 275 {
 276     int i;
 277     real nom=0;  /* nominator */
 278     real denom=0; /* denominator */
 279     real tmp=0;
 280
 281     if ( m == 0.0 )
 282         return 0.0 ;
 283
 284     for(i=-SUMORDER ; i<0 ; i++)
 285     {
 286         tmp=m / K + i;
 287         tmp*=2.0*M_PI;
 288         nom+= i * pow( tmp , -2.0*n );
 289     }
 290
 291     for(i=SUMORDER ; i>0 ; i--)
 292     {
 293         tmp=m / K + i;
 294         tmp*=2.0*M_PI;
 295         nom+= i * pow( tmp , -2.0*n );
 296     }
 297
 298     for(i=-SUMORDER ; i<SUMORDER+1 ; i++)
 299     {
 300         tmp=m / K + i;
 301         tmp*=2.0*M_PI;
 302         denom+=pow( tmp , -n );
 303     }
 304
 305     return 2.0 * M_PI * nom / denom / denom;
 306
 307 }
 308
 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 */
 313 {
 314     int i;
 315     real nom=0;  /* nominator */
 316     real denom=0; /* denominator */
 317     real tmp=0;
 318
 319     if ( m == 0.0 )
 320         return 0.0 ;
 321
 322     for(i=-SUMORDER ; i<0 ; i++)
 323     {
 324         tmp=m / K + i;
 325         tmp*=2.0*M_PI;
 326         nom+= i * i * pow( tmp , -2.0*n );
 327     }
 328
 329     for(i=SUMORDER ; i>0 ; i--)
 330     {
 331         tmp=m / K + i;
 332         tmp*=2.0*M_PI;
 333         nom+= i * i * pow( tmp , -2.0*n );
 334     }
 335
 336     for(i=-SUMORDER ; i<SUMORDER+1 ; i++)
 337     {
 338         tmp=m / K + i;
 339         tmp*=2.0*M_PI;
 340         denom+=pow( tmp , -n );
 341     }
 342
 343     return 4.0 * M_PI * M_PI * nom / denom / denom;
 344
 345 }
 346
 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 */
 353 {
 354     int i;
 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 */
 361
 362
 363     if ( m == 0.0 )
 364         return 0.0 ;
 365
 366     rcoord=iprod(rboxv,x);
 367
 368
 369     for(i=-SUMORDER;i<0;i++)
 370     {
 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);
 374         nom+=tmp * tmp2 * i;
 375         denom+=tmp2;
 376     }
 377
 378     for(i=SUMORDER;i>0;i--)
 379     {
 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);
 383         nom+=tmp * tmp2 * i;
 384         denom+=tmp2;
 385     }
 386
 387
 388     tmp=2.0 * M_PI * m / K;
 389     tmp1=pow(tmp,-1.0*n);
 390     denom+=tmp1;
 391
 392    return 2.0 * M_PI * nom / denom * K  ;
 393
 394 }
 395
 396 #undef SUMORDER
 397
 398 /* The following routine is just a copy from pme.c */
 399
 400 static void calc_recipbox(matrix box,matrix recipbox)
 401 {
 402   /* Save some time by assuming upper right part is zero */
 403
 404   real tmp=1.0/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
 405
 406   recipbox[XX][XX]=box[YY][YY]*box[ZZ][ZZ]*tmp;
 407   recipbox[XX][YY]=0;
 408   recipbox[XX][ZZ]=0;
 409   recipbox[YY][XX]=-box[YY][XX]*box[ZZ][ZZ]*tmp;
 410   recipbox[YY][YY]=box[XX][XX]*box[ZZ][ZZ]*tmp;
 411   recipbox[YY][ZZ]=0;
 412   recipbox[ZZ][XX]=(box[YY][XX]*box[ZZ][YY]-box[YY][YY]*box[ZZ][XX])*tmp;
 413   recipbox[ZZ][YY]=-box[ZZ][YY]*box[XX][XX]*tmp;
 414   recipbox[ZZ][ZZ]=box[XX][XX]*box[YY][YY]*tmp;
 415 }
 416
 417
 418 /* Estimate the reciprocal space part error of the SPME Ewald sum. */
 419 static real estimate_reciprocal(
 420         t_inputinfo *info,
 421         rvec x[],           /* array of particles */
 422         real q[],           /* array of charges */
 423         int nr,             /* number of charges = size of the charge array */
 424         FILE *fp_out,
 425         gmx_bool bVerbose,
 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 */
 429         t_commrec *cr)
 430 {
 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 */
 438     int i,ci;
 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 */
 449     real tmp1=0;
 450     real tmp2=0;
 451     gmx_bool bFraction;
 452
 453     /* Random number generator */
 454     gmx_rng_t rng=NULL;
 455     int *numbers=NULL;
 456
 457     /* Index variables for parallel work distribution */
 458     int startglobal,stopglobal;
 459     int startlocal, stoplocal;
 460     int x_per_core;
 461     int xtot;
 462
 463 #ifdef TAKETIME
 464     double t0=0.0;
 465     double t1=0.0;
 466 #endif
 467
 468     rng=gmx_rng_init(seed);
 469
 470     clear_rvec(gridpx);
 471     clear_rvec(gridpxy);
 472     clear_rvec(gridp);
 473     clear_rvec(tmpvec);
 474     clear_rvec(tmpvec2);
 475
 476     for(i=0;i<nr;i++)
 477     {
 478         q2_all += q[i]*q[i];
 479     }
 480
 481     /* Calculate indices for work distribution */
 482     startglobal=-info->nkx[0]/2;
 483     stopglobal = info->nkx[0]/2;
 484     xtot = stopglobal*2+1;
 485     if (PAR(cr))
 486     {
 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;
 492     }
 493     else
 494     {
 495         startlocal = startglobal;
 496         stoplocal  = stopglobal;
 497         x_per_core = xtot;
 498     }
 499 /*
 500 #ifdef GMX_LIB_MPI
 501     MPI_Barrier(MPI_COMM_WORLD);
 502 #endif
 503 */
 504
 505 #ifdef GMX_LIB_MPI
 506 #ifdef TAKETIME
 507     if (MASTER(cr))
 508         t0 = MPI_Wtime();
 509 #endif
 510 #endif
 511
 512     if (MASTER(cr)){
 513
 514         fprintf(stderr, "Calculating reciprocal error part 1 ...");
 515
 516     }
 517
 518     for(nx=startlocal; nx<=stoplocal; nx++)
 519     {
 520         svmul(nx,info->recipbox[XX],gridpx);
 521         for(ny=-info->nky[0]/2; ny<info->nky[0]/2+1; ny++)
 522         {
 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++)
 526             {
 527                 if (  0 == nx &&  0 == ny &&  0 == nz )
 528                     continue;
 529                 svmul(nz,info->recipbox[ZZ],tmpvec);
 530                 rvec_add(gridpxy,tmpvec,gridp);
 531                 tmp=norm2(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;
 534                 coeff2=tmp ;
 535
 536
 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]);
 540
 541                 tmp1=eps_poly1(nx,info->nkx[0],info->pme_order[0]);
 542                 tmp2=eps_poly1(ny,info->nky[0],info->pme_order[0]);
 543
 544                 tmp+=2.0 * tmp1 * tmp2;
 545
 546                 tmp1=eps_poly1(nz,info->nkz[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(nx,info->nkx[0],info->pme_order[0]);
 553
 554                 tmp+=2.0 * tmp1 * tmp2;
 555
 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]);
 559
 560                 tmp+= tmp1 * tmp1;
 561
 562                 e_rec1+= 32.0 * M_PI * M_PI * coeff * coeff * coeff2 * tmp  * q2_all * q2_all / nr ;
 563
 564                 tmp1=eps_poly3(nx,info->nkx[0],info->pme_order[0]);
 565                 tmp1*=info->nkx[0];
 566                 tmp2=iprod(gridp,info->recipbox[XX]);
 567
 568                 tmp=tmp1*tmp2;
 569
 570                 tmp1=eps_poly3(ny,info->nky[0],info->pme_order[0]);
 571                 tmp1*=info->nky[0];
 572                 tmp2=iprod(gridp,info->recipbox[YY]);
 573
 574                 tmp+=tmp1*tmp2;
 575
 576                 tmp1=eps_poly3(nz,info->nkz[0],info->pme_order[0]);
 577                 tmp1*=info->nkz[0];
 578                 tmp2=iprod(gridp,info->recipbox[ZZ]);
 579
 580                 tmp+=tmp1*tmp2;
 581
 582                 tmp*=4.0 * M_PI;
 583
 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];
 587
 588                 tmp+=tmp1;
 589
 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];
 593
 594                 tmp+=tmp1;
 595
 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];
 599
 600                 tmp+=tmp1;
 601
 602                 e_rec2+= 4.0 * coeff * coeff * tmp * q2_all * q2_all / nr ;
 603
 604             }
 605         }
 606         if (MASTER(cr))
 607             fprintf(stderr, "\rCalculating reciprocal error part 1 ... %3.0f%%", 100.0*(nx-startlocal+1)/(x_per_core));
 608
 609     }
 610
 611     if (MASTER(cr))
 612         fprintf(stderr, "\n");
 613
 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);
 616
 617     if (bFraction)
 618     {
 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;
 624     }
 625     else
 626     {
 627         /* In this case we use all nr particle positions */
 628         xtot = nr;
 629         x_per_core = ceil( (real)xtot / (real)cr->nnodes );
 630     }
 631
 632     startlocal = x_per_core *  cr->nodeid;
 633     stoplocal  = min(startlocal + x_per_core, xtot);  /* min needed if xtot == nr */
 634
 635     if (bFraction)
 636     {
 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 */
 640         snew(numbers, xtot);
 641         if (MASTER(cr))
 642         {
 643             for (i=0; i<xtot; i++)
 644             {
 645                 numbers[i] = floor(gmx_rng_uniform_real(rng) * nr );
 646             }
 647         }
 648         /* Broadcast the random number array to the other nodes */
 649         if (PAR(cr))
 650         {
 651             nblock_bc(cr,xtot,numbers);
 652         }
 653
 654         if (bVerbose && MASTER(cr))
 655         {
 656             fprintf(stdout, "Using %d sample%s to approximate the self interaction error term",
 657                     xtot, xtot==1?"":"s");
 658             if (PAR(cr))
 659                 fprintf(stdout, " (%d sample%s per node)", x_per_core, x_per_core==1?"":"s");
 660             fprintf(stdout, ".\n");
 661         }
 662     }
 663
 664     /* Return the number of positions used for the Monte Carlo algorithm */
 665     *nsamples = xtot;
 666
 667     for(i=startlocal;i<stoplocal;i++)
 668     {
 669         e_rec3x=0;
 670         e_rec3y=0;
 671         e_rec3z=0;
 672
 673         if (bFraction)
 674         {
 675             /* Randomly pick a charge */
 676             ci = numbers[i];
 677         }
 678         else
 679         {
 680             /* Use all charges */
 681             ci = i;
 682         }
 683
 684         /* for(nx=startlocal; nx<=stoplocal; nx++)*/
 685         for(nx=-info->nkx[0]/2; nx<info->nkx[0]/2+1; nx++)
 686         {
 687             svmul(nx,info->recipbox[XX],gridpx);
 688             for(ny=-info->nky[0]/2; ny<info->nky[0]/2+1; ny++)
 689             {
 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++)
 693                 {
 694
 695                     if (  0 == nx && 0 == ny && 0 == nz)
 696                         continue;
 697
 698                     svmul(nz,info->recipbox[ZZ],tmpvec);
 699                     rvec_add(gridpxy,tmpvec,gridp);
 700                     tmp=norm2(gridp);
 701                     coeff=exp(-1.0 * M_PI * M_PI * tmp / info->ewald_beta[0] / info->ewald_beta[0] );
 702                     coeff/= tmp ;
 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]);
 706
 707                 }
 708             }
 709         }
 710
 711         clear_rvec(tmpvec2);
 712
 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);
 719
 720         e_rec3 += q[ci]*q[ci]*q[ci]*q[ci]*norm2(tmpvec2) / ( xtot * M_PI * info->volume * M_PI * info->volume);
 721         if (MASTER(cr)){
 722             fprintf(stderr, "\rCalculating reciprocal error part 2 ... %3.0f%%",
 723                     100.0*(i+1)/stoplocal);
 724
 725         }
 726     }
 727
 728     if (MASTER(cr))
 729         fprintf(stderr, "\n");
 730
 731 #ifdef GMX_LIB_MPI
 732 #ifdef TAKETIME
 733     if (MASTER(cr))
 734     {
 735         t1= MPI_Wtime() - t0;
 736         fprintf(fp_out, "Recip. err. est. took   : %lf s\n", t1);
 737     }
 738 #endif
 739 #endif
 740
 741 #ifdef DEBUG
 742     if (PAR(cr))
 743     {
 744         fprintf(stderr, "Node %3d: nx=[%3d...%3d]  e_rec3=%e\n",
 745                 cr->nodeid, startlocal, stoplocal, e_rec3);
 746     }
 747 #endif
 748
 749     if (PAR(cr))
 750     {
 751         gmx_sum(1,&e_rec1,cr);
 752         gmx_sum(1,&e_rec2,cr);
 753         gmx_sum(1,&e_rec3,cr);
 754     }
 755
 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 ;
 759      */
 760     e_rec=sqrt(e_rec1+e_rec2+e_rec3);
 761
 762
 763     return ONE_4PI_EPS0 * e_rec;
 764 }
 765
 766
 767 /* Allocate memory for the inputinfo struct: */
 768 static void create_info(t_inputinfo *info)
 769 {
 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);
 783 }
 784
 785
 786 /* Allocate and fill an array with coordinates and charges,
 787  * returns the number of charges found
 788  */
 789 static int prepare_x_q(real *q[], rvec *x[], gmx_mtop_t *mtop, rvec x_orig[], t_commrec *cr)
 790 {
 791     int i;
 792     int nq; /* number of charged particles */
 793     gmx_mtop_atomloop_all_t aloop;
 794     t_atom *atom;
 795
 796
 797     if (MASTER(cr))
 798     {
 799         snew(*q, mtop->natoms);
 800         snew(*x, mtop->natoms);
 801         nq=0;
 802
 803         aloop = gmx_mtop_atomloop_all_init(mtop);
 804
 805         while (gmx_mtop_atomloop_all_next(aloop,&i,&atom))
 806         {
 807             if (is_charge(atom->q))
 808             {
 809                 (*q)[nq] = atom->q;
 810                 (*x)[nq][XX] = x_orig[i][XX];
 811                 (*x)[nq][YY] = x_orig[i][YY];
 812                 (*x)[nq][ZZ] = x_orig[i][ZZ];
 813                 nq++;
 814             }
 815         }
 816         /* Give back some unneeded memory */
 817         srenew(*q, nq);
 818         srenew(*x, nq);
 819     }
 820     /* Broadcast x and q in the parallel case */
 821     if (PAR(cr))
 822     {
 823         /* Transfer the number of charges */
 824         block_bc(cr,nq);
 825         snew_bc(cr, *x, nq);
 826         snew_bc(cr, *q, nq);
 827         nblock_bc(cr,nq,*x);
 828         nblock_bc(cr,nq,*q);
 829     }
 830
 831     return nq;
 832 }
 833
 834
 835
 836 /* Read in the tpr file and save information we need later in info */
 837 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)
 838 {
 839     read_tpx_state(fn_sim_tpr,ir,state,NULL,mtop);
 840
 841     /* The values of the original tpr input file are save in the first
 842      * place [0] of the arrays */
 843     info->orig_sim_steps = ir->nsteps;
 844     info->pme_order[0]   = ir->pme_order;
 845     info->rcoulomb[0]    = ir->rcoulomb;
 846     info->rvdw[0]        = ir->rvdw;
 847     info->nkx[0]         = ir->nkx;
 848     info->nky[0]         = ir->nky;
 849     info->nkz[0]         = ir->nkz;
 850     info->ewald_rtol[0]  = ir->ewald_rtol;
 851     info->fracself       = fracself;
 852     if (user_beta > 0)
 853         info->ewald_beta[0] = user_beta;
 854     else
 855         info->ewald_beta[0]  = calc_ewaldcoeff(info->rcoulomb[0],info->ewald_rtol[0]);
 856
 857     /* Check if PME was chosen */
 858     if (EEL_PME(ir->coulombtype) == FALSE)
 859         gmx_fatal(FARGS, "Can only do optimizations for simulations with PME");
 860
 861     /* Check if rcoulomb == rlist, which is necessary for PME */
 862     if (!(ir->rcoulomb == ir->rlist))
 863         gmx_fatal(FARGS, "PME requires rcoulomb (%f) to be equal to rlist (%f).", ir->rcoulomb, ir->rlist);
 864 }
 865
 866
 867 /* Transfer what we need for parallelizing the reciprocal error estimate */
 868 static void bcast_info(t_inputinfo *info, t_commrec *cr)
 869 {
 870     nblock_bc(cr, info->n_entries, info->nkx);
 871     nblock_bc(cr, info->n_entries, info->nky);
 872     nblock_bc(cr, info->n_entries, info->nkz);
 873     nblock_bc(cr, info->n_entries, info->ewald_beta);
 874     nblock_bc(cr, info->n_entries, info->pme_order);
 875     nblock_bc(cr, info->n_entries, info->e_dir);
 876     nblock_bc(cr, info->n_entries, info->e_rec);
 877     block_bc(cr, info->volume);
 878     block_bc(cr, info->recipbox);
 879     block_bc(cr, info->natoms);
 880     block_bc(cr, info->fracself);
 881     block_bc(cr, info->bTUNE);
 882     block_bc(cr, info->q2all);
 883     block_bc(cr, info->q2allnr);
 884 }
 885
 886
 887 /* Estimate the error of the SPME Ewald sum. This estimate is based upon
 888  * a) a homogeneous distribution of the charges
 889  * b) a total charge of zero.
 890  */
 891 static void estimate_PME_error(t_inputinfo *info, t_state *state,
 892         gmx_mtop_t *mtop, FILE *fp_out, gmx_bool bVerbose, unsigned int seed,
 893         t_commrec *cr)
 894 {
 895     rvec *x=NULL; /* The coordinates */
 896     real *q=NULL; /* The charges     */
 897     real  edir=0.0; /* real space error */
 898     real  erec=0.0; /* reciprocal space error */
 899     real  derr=0.0; /* difference of real and reciprocal space error */
 900     real  derr0=0.0; /* difference of real and reciprocal space error */
 901     real  beta=0.0; /* splitting parameter beta */
 902     real  beta0=0.0; /* splitting parameter beta */
 903     int ncharges; /* The number of atoms with charges */
 904     int nsamples; /* The number of samples used for the calculation of the
 905                    * self-energy error term */
 906     int i=0;
 907
 908     if (MASTER(cr))
 909         fprintf(fp_out, "\n--- PME ERROR ESTIMATE ---\n");
 910
 911     /* Prepare an x and q array with only the charged atoms */
 912     ncharges = prepare_x_q(&q, &x, mtop, state->x, cr);
 913     if (MASTER(cr))
 914     {
 915         calc_q2all(mtop, &(info->q2all), &(info->q2allnr));
 916         info->ewald_rtol[0]=gmx_erfc(info->rcoulomb[0]*info->ewald_beta[0]);
 917         /* Write some info to log file */
 918         fprintf(fp_out, "Box volume              : %g nm^3\n", info->volume);
 919         fprintf(fp_out, "Number of charged atoms : %d (total atoms %d)\n",ncharges, info->natoms);
 920         fprintf(fp_out, "Coulomb radius          : %g nm\n", info->rcoulomb[0]);
 921         fprintf(fp_out, "Ewald_rtol              : %g\n", info->ewald_rtol[0]);
 922         fprintf(fp_out, "Ewald parameter beta    : %g\n", info->ewald_beta[0]);
 923         fprintf(fp_out, "Interpolation order     : %d\n", info->pme_order[0]);
 924         fprintf(fp_out, "Fourier grid (nx,ny,nz) : %d x %d x %d\n",
 925                 info->nkx[0],info->nky[0],info->nkz[0]);
 926         fflush(fp_out);
 927
 928     }
 929
 930     if (PAR(cr))
 931         bcast_info(info, cr);
 932
 933
 934     /* Calculate direct space error */
 935     info->e_dir[0] = estimate_direct(info);
 936
 937     /* Calculate reciprocal space error */
 938     info->e_rec[0] = estimate_reciprocal(info, x, q, ncharges, fp_out, bVerbose,
 939                                          seed, &nsamples, cr);
 940
 941     if (PAR(cr))
 942         bcast_info(info, cr);
 943
 944     if (MASTER(cr))
 945     {
 946         fprintf(fp_out, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
 947         fprintf(fp_out, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
 948         fprintf(fp_out, "Self-energy error term was estimated using %d samples\n", nsamples);
 949         fflush(fp_out);
 950         fprintf(stderr, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
 951         fprintf(stderr, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
 952     }
 953
 954     i=0;
 955
 956     if (info->bTUNE)
 957     {
 958         if(MASTER(cr))
 959             fprintf(stderr,"Starting tuning ...\n");
 960         edir=info->e_dir[0];
 961         erec=info->e_rec[0];
 962         derr0=edir-erec;
 963         beta0=info->ewald_beta[0];
 964         if (derr>0.0)
 965             info->ewald_beta[0]+=0.1;
 966         else
 967             info->ewald_beta[0]-=0.1;
 968         info->e_dir[0] = estimate_direct(info);
 969         info->e_rec[0] = estimate_reciprocal(info, x, q, ncharges, fp_out, bVerbose,
 970                                              seed, &nsamples, cr);
 971
 972         if (PAR(cr))
 973             bcast_info(info, cr);
 974
 975
 976         edir=info->e_dir[0];
 977         erec=info->e_rec[0];
 978         derr=edir-erec;
 979         while ( fabs(derr/min(erec,edir)) > 1e-4)
 980         {
 981
 982             beta=info->ewald_beta[0];
 983             beta-=derr*(info->ewald_beta[0]-beta0)/(derr-derr0);
 984             beta0=info->ewald_beta[0];
 985             info->ewald_beta[0]=beta;
 986             derr0=derr;
 987
 988             info->e_dir[0] = estimate_direct(info);
 989             info->e_rec[0] = estimate_reciprocal(info, x, q, ncharges, fp_out, bVerbose,
 990                                                  seed, &nsamples, cr);
 991
 992             if (PAR(cr))
 993                 bcast_info(info, cr);
 994
 995             edir=info->e_dir[0];
 996             erec=info->e_rec[0];
 997             derr=edir-erec;
 998
 999             if (MASTER(cr))
1000             {
1001                 i++;
1002                 fprintf(stderr,"difference between real and rec. space error (step %d): %g\n",i,fabs(derr));
1003                 fprintf(stderr,"old beta: %f\n",beta0);
1004                 fprintf(stderr,"new beta: %f\n",beta);
1005             }
1006         }
1007
1008         info->ewald_rtol[0]=gmx_erfc(info->rcoulomb[0]*info->ewald_beta[0]);
1009
1010         if (MASTER(cr))
1011         {
1012             /* Write some info to log file */
1013             fflush(fp_out);
1014             fprintf(fp_out, "=========  After tuning ========\n");
1015             fprintf(fp_out, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
1016             fprintf(fp_out, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
1017             fprintf(stderr, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
1018             fprintf(stderr, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
1019             fprintf(fp_out, "Ewald_rtol              : %g\n", info->ewald_rtol[0]);
1020             fprintf(fp_out, "Ewald parameter beta    : %g\n", info->ewald_beta[0]);
1021             fflush(fp_out);
1022
1023         }
1024
1025     }
1026
1027 }
1028
1029
1030 int gmx_pme_error(int argc,char *argv[])
1031 {
1032     const char *desc[] = {
1033             "[TT]g_pme_error[tt] estimates the error of the electrostatic forces",
1034             "if using the sPME algorithm. The flag [TT]-tune[tt] will determine",
1035             "the splitting parameter such that the error is equally",
1036             "distributed over the real and reciprocal space part.",
1037             "The part of the error that stems from self interaction of the particles "
1038             "is computationally demanding. However, a good a approximation is to",
1039             "just use a fraction of the particles for this term which can be",
1040             "indicated by the flag [TT]-self[tt].[PAR]",
1041     };
1042
1043     real        fs=0.0;             /* 0 indicates: not set by the user */
1044     real        user_beta=-1.0;
1045     real        fracself=1.0;
1046     t_inputinfo info;
1047     t_state     state;     /* The state from the tpr input file */
1048     gmx_mtop_t  mtop;      /* The topology from the tpr input file */
1049     t_inputrec  *ir=NULL;  /* The inputrec from the tpr file */
1050     FILE        *fp=NULL;
1051     t_commrec   *cr;
1052     unsigned long PCA_Flags;
1053     gmx_bool        bTUNE=FALSE;
1054     gmx_bool    bVerbose=FALSE;
1055     int         seed=0;
1056
1057
1058     static t_filenm fnm[] = {
1059       { efTPX, "-s",     NULL,    ffREAD },
1060       { efOUT, "-o",    "error",  ffWRITE },
1061       { efTPX, "-so",   "tuned",  ffOPTWR }
1062     };
1063
1064     output_env_t oenv=NULL;
1065
1066     t_pargs pa[] = {
1067         { "-beta",     FALSE, etREAL, {&user_beta},
1068             "If positive, overwrite ewald_beta from [TT].tpr[tt] file with this value" },
1069         { "-tune",     FALSE, etBOOL, {&bTUNE},
1070             "Tune the splitting parameter such that the error is equally distributed between real and reciprocal space" },
1071         { "-self",     FALSE, etREAL, {&fracself},
1072             "If between 0.0 and 1.0, determine self interaction error from just this fraction of the charged particles" },
1073         { "-seed",     FALSE, etINT,  {&seed},
1074           "Random number seed used for Monte Carlo algorithm when [TT]-self[tt] is set to a value between 0.0 and 1.0" },
1075         { "-v",        FALSE, etBOOL, {&bVerbose},
1076             "Be loud and noisy" }
1077     };
1078
1079
1080 #define NFILE asize(fnm)
1081
1082     cr = init_par(&argc,&argv);
1083
1084 #ifdef GMX_LIB_MPI
1085     MPI_Barrier(MPI_COMM_WORLD);
1086 #endif
1087
1088     if (MASTER(cr))
1089       CopyRight(stderr,argv[0]);
1090
1091     PCA_Flags = PCA_NOEXIT_ON_ARGS;
1092     PCA_Flags |= (MASTER(cr) ? 0 : PCA_QUIET);
1093
1094     parse_common_args(&argc,argv,PCA_Flags,
1095                       NFILE,fnm,asize(pa),pa,asize(desc),desc,
1096                       0,NULL,&oenv);
1097
1098     if (!bTUNE)
1099         bTUNE = opt2bSet("-so",NFILE,fnm);
1100
1101     info.n_entries = 1;
1102
1103     /* Allocate memory for the inputinfo struct: */
1104     create_info(&info);
1105     info.fourier_sp[0] = fs;
1106
1107     /* Read in the tpr file and open logfile for reading */
1108     if (MASTER(cr))
1109     {
1110         snew(ir,1);
1111         read_tpr_file(opt2fn("-s",NFILE,fnm), &info, &state, &mtop, ir, user_beta,fracself);
1112
1113         fp=fopen(opt2fn("-o",NFILE,fnm),"w");
1114     }
1115
1116     /* Check consistency if the user provided fourierspacing */
1117     if (fs > 0 && MASTER(cr))
1118     {
1119         /* Recalculate the grid dimensions using fourierspacing from user input */
1120         info.nkx[0] = 0;
1121         info.nky[0] = 0;
1122         info.nkz[0] = 0;
1123         calc_grid(stdout,state.box,info.fourier_sp[0],&(info.nkx[0]),&(info.nky[0]),&(info.nkz[0]));
1124         if ( (ir->nkx != info.nkx[0]) || (ir->nky != info.nky[0]) || (ir->nkz != info.nkz[0]) )
1125             gmx_fatal(FARGS, "Wrong fourierspacing %f nm, input file grid = %d x %d x %d, computed grid = %d x %d x %d",
1126                       fs,ir->nkx,ir->nky,ir->nkz,info.nkx[0],info.nky[0],info.nkz[0]);
1127     }
1128
1129     /* Estimate (S)PME force error */
1130
1131     /* Determine the volume of the simulation box */
1132     if (MASTER(cr))
1133     {
1134         info.volume = det(state.box);
1135         calc_recipbox(state.box,info.recipbox);
1136         info.natoms = mtop.natoms;
1137         info.bTUNE  = bTUNE;
1138     }
1139
1140     if (PAR(cr))
1141         bcast_info(&info, cr);
1142
1143     /* Get an error estimate of the input tpr file and do some tuning if requested */
1144     estimate_PME_error(&info, &state, &mtop, fp, bVerbose, seed, cr);
1145
1146     if (MASTER(cr))
1147     {
1148         /* Write out optimized tpr file if requested */
1149         if ( opt2bSet("-so",NFILE,fnm) || bTUNE )
1150         {
1151             ir->ewald_rtol=info.ewald_rtol[0];
1152             write_tpx_state(opt2fn("-so",NFILE,fnm),ir,&state,&mtop);
1153         }
1154         please_cite(fp,"Wang2010");
1155         fclose(fp);
1156     }
1157
1158     gmx_finalize_par();
1159
1160     return 0;
1161 }