src/gromacs/gmxana/gmx_pme_error.cpp

   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 <algorithm>
  34
  35 #include "statutil.h"
  36 #include "typedefs.h"
  37 #include "smalloc.h"
  38 #include "vec.h"
  39 #include "copyrite.h"
  40 #include "tpxio.h"
  41 #include "string2.h"
  42 #include "readinp.h"
  43 #include "calcgrid.h"
  44 #include "checkpoint.h"
  45 #include "gmx_ana.h"
  46 #include "gmx_random.h"
  47 #include "physics.h"
  48 #include "mdatoms.h"
  49 #include "coulomb.h"
  50 #include "mtop_util.h"
  51 #include "network.h"
  52 #include "main.h"
  53 #include "macros.h"
  54
  55 /* We use the same defines as in mvdata.c here */
  56 #define  block_bc(cr,   d) gmx_bcast(     sizeof(d),     &(d), (cr))
  57 #define nblock_bc(cr, nr, d) gmx_bcast((nr)*sizeof((d)[0]), (d), (cr))
  58 #define   snew_bc(cr, d, nr) { if (!MASTER(cr)) {snew((d), (nr)); }}
  59 /* #define TAKETIME */
  60 /* #define DEBUG  */
  61 enum {
  62     ddnoSEL, ddnoINTERLEAVE, ddnoPP_PME, ddnoCARTESIAN, ddnoNR
  63 };
  64
  65 /* Enum for situations that can occur during log file parsing */
  66 enum {
  67     eParselogOK,
  68     eParselogNotFound,
  69     eParselogNoPerfData,
  70     eParselogTerm,
  71     eParselogResetProblem,
  72     eParselogNr
  73 };
  74
  75
  76 typedef struct
  77 {
  78     gmx_large_int_t orig_sim_steps;  /* Number of steps to be done in the real simulation  */
  79     int             n_entries;       /* Number of entries in arrays                        */
  80     real            volume;          /* The volume of the box                              */
  81     matrix          recipbox;        /* The reciprocal box                                 */
  82     int             natoms;          /* The number of atoms in the MD system               */
  83     real           *fac;             /* The scaling factor                                 */
  84     real           *rcoulomb;        /* The coulomb radii [0...nr_inputfiles]              */
  85     real           *rvdw;            /* The vdW radii                                      */
  86     int            *nkx, *nky, *nkz; /* Number of k vectors in each spatial dimension      */
  87     real           *fourier_sp;      /* Fourierspacing                                     */
  88     real           *ewald_rtol;      /* Real space tolerance for Ewald, determines         */
  89                                      /* the real/reciprocal space relative weight          */
  90     real           *ewald_beta;      /* Splitting parameter [1/nm]                         */
  91     real            fracself;        /* fraction of particles for SI error                 */
  92     real            q2all;           /* sum ( q ^2 )                                       */
  93     real            q2allnr;         /* nr of charges                                      */
  94     int            *pme_order;       /* Interpolation order for PME (bsplines)             */
  95     char          **fn_out;          /* Name of the output tpr file                        */
  96     real           *e_dir;           /* Direct space part of PME error with these settings */
  97     real           *e_rec;           /* Reciprocal space part of PME error                 */
  98     gmx_bool        bTUNE;           /* flag for tuning */
  99 } t_inputinfo;
 100
 101
 102 /* Returns TRUE when atom is charged */
 103 static gmx_bool is_charge(real charge)
 104 {
 105     if (charge*charge > GMX_REAL_EPS)
 106     {
 107         return TRUE;
 108     }
 109     else
 110     {
 111         return FALSE;
 112     }
 113 }
 114
 115
 116 /* calculate charge density */
 117 static void calc_q2all(
 118         gmx_mtop_t *mtop,   /* molecular topology */
 119         real *q2all, real *q2allnr)
 120 {
 121     int             imol, iatom; /* indices for loops */
 122     real            q2_all = 0;  /* Sum of squared charges */
 123     int             nrq_mol;     /* Number of charges in a single molecule */
 124     int             nrq_all;     /* Total number of charges in the MD system */
 125     real            qi, q2_mol;
 126     gmx_moltype_t  *molecule;
 127     gmx_molblock_t *molblock;
 128
 129 #ifdef DEBUG
 130     fprintf(stderr, "\nCharge density:\n");
 131 #endif
 132     q2_all  = 0.0;                                 /* total q squared */
 133     nrq_all = 0;                                   /* total number of charges in the system */
 134     for (imol = 0; imol < mtop->nmolblock; imol++) /* Loop over molecule types */
 135     {
 136         q2_mol   = 0.0;                            /* q squared value of this molecule */
 137         nrq_mol  = 0;                              /* number of charges this molecule carries */
 138         molecule = &(mtop->moltype[imol]);
 139         molblock = &(mtop->molblock[imol]);
 140         for (iatom = 0; iatom < molblock->natoms_mol; iatom++) /* Loop over atoms in this molecule */
 141         {
 142             qi = molecule->atoms.atom[iatom].q;                /* Charge of this atom */
 143             /* Is this charge worth to be considered? */
 144             if (is_charge(qi))
 145             {
 146                 q2_mol += qi*qi;
 147                 nrq_mol++;
 148             }
 149         }
 150         /* Multiply with the number of molecules present of this type and add */
 151         q2_all  += q2_mol*molblock->nmol;
 152         nrq_all += nrq_mol*molblock->nmol;
 153 #ifdef DEBUG
 154         fprintf(stderr, "Molecule %2d (%5d atoms) q2_mol=%10.3e nr.mol.charges=%5d (%6dx)  q2_all=%10.3e  tot.charges=%d\n",
 155                 imol, molblock->natoms_mol, q2_mol, nrq_mol, molblock->nmol, q2_all, nrq_all);
 156 #endif
 157     }
 158
 159     *q2all   = q2_all;
 160     *q2allnr = nrq_all;
 161
 162 }
 163
 164
 165 /* Estimate the direct space part error of the SPME Ewald sum */
 166 static real estimate_direct(
 167         t_inputinfo *info
 168         )
 169 {
 170     real e_dir     = 0; /* Error estimate */
 171     real beta      = 0; /* Splitting parameter (1/nm) */
 172     real r_coulomb = 0; /* Cut-off in direct space */
 173
 174
 175     beta      = info->ewald_beta[0];
 176     r_coulomb = info->rcoulomb[0];
 177
 178     e_dir  = 2.0 * info->q2all * gmx_invsqrt( info->q2allnr  *  r_coulomb * info->volume );
 179     e_dir *= exp (-beta*beta*r_coulomb*r_coulomb);
 180
 181     return ONE_4PI_EPS0*e_dir;
 182 }
 183
 184 #define SUMORDER 6
 185
 186 /* the following 4 functions determine polynomials required for the reciprocal error estimate */
 187
 188 static inline real eps_poly1(
 189         real m,          /* grid coordinate in certain direction */
 190         real K,          /* grid size in corresponding direction */
 191         real n)          /* spline interpolation order of the SPME */
 192 {
 193     int  i;
 194     real nom   = 0; /* nominator */
 195     real denom = 0; /* denominator */
 196     real tmp   = 0;
 197
 198     if (m == 0.0)
 199     {
 200         return 0.0;
 201     }
 202
 203     for (i = -SUMORDER; i < 0; i++)
 204     {
 205         tmp  = m / K + i;
 206         tmp *= 2.0*M_PI;
 207         nom += pow( tmp, -n );
 208     }
 209
 210     for (i = SUMORDER; i > 0; i--)
 211     {
 212         tmp  = m / K + i;
 213         tmp *= 2.0*M_PI;
 214         nom += pow( tmp, -n );
 215     }
 216
 217     tmp   = m / K;
 218     tmp  *= 2.0*M_PI;
 219     denom = pow( tmp, -n )+nom;
 220
 221     return -nom/denom;
 222
 223 }
 224
 225 static inline real eps_poly2(
 226         real m,          /* grid coordinate in certain direction */
 227         real K,          /* grid size in corresponding direction */
 228         real n)          /* spline interpolation order of the SPME */
 229 {
 230     int  i;
 231     real nom   = 0; /* nominator */
 232     real denom = 0; /* denominator */
 233     real tmp   = 0;
 234
 235     if (m == 0.0)
 236     {
 237         return 0.0;
 238     }
 239
 240     for (i = -SUMORDER; i < 0; i++)
 241     {
 242         tmp  = m / K + i;
 243         tmp *= 2.0*M_PI;
 244         nom += pow( tmp, -2*n );
 245     }
 246
 247     for (i = SUMORDER; i > 0; i--)
 248     {
 249         tmp  = m / K + i;
 250         tmp *= 2.0*M_PI;
 251         nom += pow( tmp, -2*n );
 252     }
 253
 254     for (i = -SUMORDER; i < SUMORDER+1; i++)
 255     {
 256         tmp    = m / K + i;
 257         tmp   *= 2.0*M_PI;
 258         denom += pow( tmp, -n );
 259     }
 260     tmp = eps_poly1(m, K, n);
 261     return nom / denom / denom + tmp*tmp;
 262
 263 }
 264
 265 static inline real eps_poly3(
 266         real m,          /* grid coordinate in certain direction */
 267         real K,          /* grid size in corresponding direction */
 268         real n)          /* spline interpolation order of the SPME */
 269 {
 270     int  i;
 271     real nom   = 0; /* nominator */
 272     real denom = 0; /* denominator */
 273     real tmp   = 0;
 274
 275     if (m == 0.0)
 276     {
 277         return 0.0;
 278     }
 279
 280     for (i = -SUMORDER; i < 0; i++)
 281     {
 282         tmp  = m / K + i;
 283         tmp *= 2.0*M_PI;
 284         nom += i * pow( tmp, -2*n );
 285     }
 286
 287     for (i = SUMORDER; i > 0; i--)
 288     {
 289         tmp  = m / K + i;
 290         tmp *= 2.0*M_PI;
 291         nom += i * pow( tmp, -2*n );
 292     }
 293
 294     for (i = -SUMORDER; i < SUMORDER+1; i++)
 295     {
 296         tmp    = m / K + i;
 297         tmp   *= 2.0*M_PI;
 298         denom += pow( tmp, -n );
 299     }
 300
 301     return 2.0 * M_PI * nom / denom / denom;
 302
 303 }
 304
 305 static inline real eps_poly4(
 306         real m,          /* grid coordinate in certain direction */
 307         real K,          /* grid size in corresponding direction */
 308         real n)          /* spline interpolation order of the SPME */
 309 {
 310     int  i;
 311     real nom   = 0; /* nominator */
 312     real denom = 0; /* denominator */
 313     real tmp   = 0;
 314
 315     if (m == 0.0)
 316     {
 317         return 0.0;
 318     }
 319
 320     for (i = -SUMORDER; i < 0; i++)
 321     {
 322         tmp  = m / K + i;
 323         tmp *= 2.0*M_PI;
 324         nom += i * i * pow( tmp, -2*n );
 325     }
 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*n );
 332     }
 333
 334     for (i = -SUMORDER; i < SUMORDER+1; i++)
 335     {
 336         tmp    = m / K + i;
 337         tmp   *= 2.0*M_PI;
 338         denom += pow( tmp, -n );
 339     }
 340
 341     return 4.0 * M_PI * M_PI * nom / denom / denom;
 342
 343 }
 344
 345 static inline real eps_self(
 346         real m,       /* grid coordinate in certain direction */
 347         real K,       /* grid size in corresponding direction */
 348         rvec rboxv,   /* reciprocal box vector */
 349         real n,       /* spline interpolation order of the SPME */
 350         rvec x)       /* coordinate of charge */
 351 {
 352     int  i;
 353     real tmp    = 0; /* temporary variables for computations */
 354     real tmp1   = 0; /* temporary variables for computations */
 355     real tmp2   = 0; /* temporary variables for computations */
 356     real rcoord = 0; /* coordinate in certain reciprocal space direction */
 357     real nom    = 0; /* nominator */
 358     real denom  = 0; /* denominator */
 359
 360
 361     if (m == 0.0)
 362     {
 363         return 0.0;
 364     }
 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, -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, -n);
 383         nom   += tmp * tmp2 * i;
 384         denom += tmp2;
 385     }
 386
 387
 388     tmp    = 2.0 * M_PI * m / K;
 389     tmp1   = pow(tmp, -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  gmx_unused  *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 = static_cast<int>(ceil(static_cast<real>(xtot) / cr->nnodes));
 488         startlocal = startglobal + x_per_core*cr->nodeid;
 489         stoplocal  = startlocal + x_per_core -1;
 490         if (stoplocal > stopglobal)
 491         {
 492             stoplocal = stopglobal;
 493         }
 494     }
 495     else
 496     {
 497         startlocal = startglobal;
 498         stoplocal  = stopglobal;
 499         x_per_core = xtot;
 500     }
 501 /*
 502    #ifdef GMX_LIB_MPI
 503     MPI_Barrier(MPI_COMM_WORLD);
 504    #endif
 505  */
 506
 507 #ifdef GMX_LIB_MPI
 508 #ifdef TAKETIME
 509     if (MASTER(cr))
 510     {
 511         t0 = MPI_Wtime();
 512     }
 513 #endif
 514 #endif
 515
 516     if (MASTER(cr))
 517     {
 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                 {
 534                     continue;
 535                 }
 536                 svmul(nz, info->recipbox[ZZ], tmpvec);
 537                 rvec_add(gridpxy, tmpvec, gridp);
 538                 tmp    = norm2(gridp);
 539                 coeff  = exp(-1.0 * M_PI * M_PI * tmp / info->ewald_beta[0] / info->ewald_beta[0] );
 540                 coeff /= 2.0 * M_PI * info->volume * tmp;
 541                 coeff2 = tmp;
 542
 543
 544                 tmp  = eps_poly2(nx, info->nkx[0], info->pme_order[0]);
 545                 tmp += eps_poly2(ny, info->nkx[0], info->pme_order[0]);
 546                 tmp += eps_poly2(nz, info->nkx[0], info->pme_order[0]);
 547
 548                 tmp1 = eps_poly1(nx, info->nkx[0], info->pme_order[0]);
 549                 tmp2 = eps_poly1(ny, info->nky[0], info->pme_order[0]);
 550
 551                 tmp += 2.0 * tmp1 * tmp2;
 552
 553                 tmp1 = eps_poly1(nz, info->nkz[0], info->pme_order[0]);
 554                 tmp2 = eps_poly1(ny, info->nky[0], info->pme_order[0]);
 555
 556                 tmp += 2.0 * tmp1 * tmp2;
 557
 558                 tmp1 = eps_poly1(nz, info->nkz[0], info->pme_order[0]);
 559                 tmp2 = eps_poly1(nx, info->nkx[0], info->pme_order[0]);
 560
 561                 tmp += 2.0 * tmp1 * tmp2;
 562
 563                 tmp1  = eps_poly1(nx, info->nkx[0], info->pme_order[0]);
 564                 tmp1 += eps_poly1(ny, info->nky[0], info->pme_order[0]);
 565                 tmp1 += eps_poly1(nz, info->nkz[0], info->pme_order[0]);
 566
 567                 tmp += tmp1 * tmp1;
 568
 569                 e_rec1 += 32.0 * M_PI * M_PI * coeff * coeff * coeff2 * tmp  * q2_all * q2_all / nr;
 570
 571                 tmp1  = eps_poly3(nx, info->nkx[0], info->pme_order[0]);
 572                 tmp1 *= info->nkx[0];
 573                 tmp2  = iprod(gridp, info->recipbox[XX]);
 574
 575                 tmp = tmp1*tmp2;
 576
 577                 tmp1  = eps_poly3(ny, info->nky[0], info->pme_order[0]);
 578                 tmp1 *= info->nky[0];
 579                 tmp2  = iprod(gridp, info->recipbox[YY]);
 580
 581                 tmp += tmp1*tmp2;
 582
 583                 tmp1  = eps_poly3(nz, info->nkz[0], info->pme_order[0]);
 584                 tmp1 *= info->nkz[0];
 585                 tmp2  = iprod(gridp, info->recipbox[ZZ]);
 586
 587                 tmp += tmp1*tmp2;
 588
 589                 tmp *= 4.0 * M_PI;
 590
 591                 tmp1  = eps_poly4(nx, info->nkx[0], info->pme_order[0]);
 592                 tmp1 *= norm2(info->recipbox[XX]);
 593                 tmp1 *= info->nkx[0] * info->nkx[0];
 594
 595                 tmp += tmp1;
 596
 597                 tmp1  = eps_poly4(ny, info->nky[0], info->pme_order[0]);
 598                 tmp1 *= norm2(info->recipbox[YY]);
 599                 tmp1 *= info->nky[0] * info->nky[0];
 600
 601                 tmp += tmp1;
 602
 603                 tmp1  = eps_poly4(nz, info->nkz[0], info->pme_order[0]);
 604                 tmp1 *= norm2(info->recipbox[ZZ]);
 605                 tmp1 *= info->nkz[0] * info->nkz[0];
 606
 607                 tmp += tmp1;
 608
 609                 e_rec2 += 4.0 * coeff * coeff * tmp * q2_all * q2_all / nr;
 610
 611             }
 612         }
 613         if (MASTER(cr))
 614         {
 615             fprintf(stderr, "\rCalculating reciprocal error part 1 ... %3.0f%%", 100.0*(nx-startlocal+1)/(x_per_core));
 616         }
 617
 618     }
 619
 620     if (MASTER(cr))
 621     {
 622         fprintf(stderr, "\n");
 623     }
 624
 625     /* Use just a fraction of all charges to estimate the self energy error term? */
 626     bFraction =  (info->fracself > 0.0) && (info->fracself < 1.0);
 627
 628     if (bFraction)
 629     {
 630         /* Here xtot is the number of samples taken for the Monte Carlo calculation
 631          * of the average of term IV of equation 35 in Wang2010. Round up to a
 632          * number of samples that is divisible by the number of nodes */
 633         x_per_core  = static_cast<int>(ceil(info->fracself * nr / cr->nnodes));
 634         xtot        = x_per_core * cr->nnodes;
 635     }
 636     else
 637     {
 638         /* In this case we use all nr particle positions */
 639         xtot       = nr;
 640         x_per_core = static_cast<int>(ceil(static_cast<real>(xtot) / cr->nnodes));
 641     }
 642
 643     startlocal = x_per_core *  cr->nodeid;
 644     stoplocal  = std::min(startlocal + x_per_core, xtot);  /* min needed if xtot == nr */
 645
 646     if (bFraction)
 647     {
 648         /* Make shure we get identical results in serial and parallel. Therefore,
 649          * take the sample indices from a single, global random number array that
 650          * is constructed on the master node and that only depends on the seed */
 651         snew(numbers, xtot);
 652         if (MASTER(cr))
 653         {
 654             for (i = 0; i < xtot; i++)
 655             {
 656                 numbers[i] = static_cast<int>(floor(gmx_rng_uniform_real(rng) * nr));
 657             }
 658         }
 659         /* Broadcast the random number array to the other nodes */
 660         if (PAR(cr))
 661         {
 662             nblock_bc(cr, xtot, numbers);
 663         }
 664
 665         if (bVerbose && MASTER(cr))
 666         {
 667             fprintf(stdout, "Using %d sample%s to approximate the self interaction error term",
 668                     xtot, xtot == 1 ? "" : "s");
 669             if (PAR(cr))
 670             {
 671                 fprintf(stdout, " (%d sample%s per node)", x_per_core, x_per_core == 1 ? "" : "s");
 672             }
 673             fprintf(stdout, ".\n");
 674         }
 675     }
 676
 677     /* Return the number of positions used for the Monte Carlo algorithm */
 678     *nsamples = xtot;
 679
 680     for (i = startlocal; i < stoplocal; i++)
 681     {
 682         e_rec3x = 0;
 683         e_rec3y = 0;
 684         e_rec3z = 0;
 685
 686         if (bFraction)
 687         {
 688             /* Randomly pick a charge */
 689             ci = numbers[i];
 690         }
 691         else
 692         {
 693             /* Use all charges */
 694             ci = i;
 695         }
 696
 697         /* for(nx=startlocal; nx<=stoplocal; nx++)*/
 698         for (nx = -info->nkx[0]/2; nx < info->nkx[0]/2+1; nx++)
 699         {
 700             svmul(nx, info->recipbox[XX], gridpx);
 701             for (ny = -info->nky[0]/2; ny < info->nky[0]/2+1; ny++)
 702             {
 703                 svmul(ny, info->recipbox[YY], tmpvec);
 704                 rvec_add(gridpx, tmpvec, gridpxy);
 705                 for (nz = -info->nkz[0]/2; nz < info->nkz[0]/2+1; nz++)
 706                 {
 707
 708                     if (0 == nx && 0 == ny && 0 == nz)
 709                     {
 710                         continue;
 711                     }
 712
 713                     svmul(nz, info->recipbox[ZZ], tmpvec);
 714                     rvec_add(gridpxy, tmpvec, gridp);
 715                     tmp      = norm2(gridp);
 716                     coeff    = exp(-1.0 * M_PI * M_PI * tmp / info->ewald_beta[0] / info->ewald_beta[0] );
 717                     coeff   /= tmp;
 718                     e_rec3x += coeff*eps_self(nx, info->nkx[0], info->recipbox[XX], info->pme_order[0], x[ci]);
 719                     e_rec3y += coeff*eps_self(ny, info->nky[0], info->recipbox[YY], info->pme_order[0], x[ci]);
 720                     e_rec3z += coeff*eps_self(nz, info->nkz[0], info->recipbox[ZZ], info->pme_order[0], x[ci]);
 721
 722                 }
 723             }
 724         }
 725
 726         clear_rvec(tmpvec2);
 727
 728         svmul(e_rec3x, info->recipbox[XX], tmpvec);
 729         rvec_inc(tmpvec2, tmpvec);
 730         svmul(e_rec3y, info->recipbox[YY], tmpvec);
 731         rvec_inc(tmpvec2, tmpvec);
 732         svmul(e_rec3z, info->recipbox[ZZ], tmpvec);
 733         rvec_inc(tmpvec2, tmpvec);
 734
 735         e_rec3 += q[ci]*q[ci]*q[ci]*q[ci]*norm2(tmpvec2) / ( xtot * M_PI * info->volume * M_PI * info->volume);
 736         if (MASTER(cr))
 737         {
 738             fprintf(stderr, "\rCalculating reciprocal error part 2 ... %3.0f%%",
 739                     100.0*(i+1)/stoplocal);
 740
 741         }
 742     }
 743
 744     if (MASTER(cr))
 745     {
 746         fprintf(stderr, "\n");
 747     }
 748
 749 #ifdef GMX_LIB_MPI
 750 #ifdef TAKETIME
 751     if (MASTER(cr))
 752     {
 753         t1 = MPI_Wtime() - t0;
 754         fprintf(fp_out, "Recip. err. est. took   : %lf s\n", t1);
 755     }
 756 #endif
 757 #endif
 758
 759 #ifdef DEBUG
 760     if (PAR(cr))
 761     {
 762         fprintf(stderr, "Node %3d: nx=[%3d...%3d]  e_rec3=%e\n",
 763                 cr->nodeid, startlocal, stoplocal, e_rec3);
 764     }
 765 #endif
 766
 767     if (PAR(cr))
 768     {
 769         gmx_sum(1, &e_rec1, cr);
 770         gmx_sum(1, &e_rec2, cr);
 771         gmx_sum(1, &e_rec3, cr);
 772     }
 773
 774     /* e_rec1*=8.0 * q2_all / info->volume / info->volume / nr ;
 775        e_rec2*=  q2_all / M_PI / M_PI / info->volume / info->volume / nr ;
 776        e_rec3/= M_PI * M_PI * info->volume * info->volume * nr ;
 777      */
 778     e_rec = sqrt(e_rec1+e_rec2+e_rec3);
 779
 780
 781     return ONE_4PI_EPS0 * e_rec;
 782 }
 783
 784
 785 /* Allocate memory for the inputinfo struct: */
 786 static void create_info(t_inputinfo *info)
 787 {
 788     snew(info->fac, info->n_entries);
 789     snew(info->rcoulomb, info->n_entries);
 790     snew(info->rvdw, info->n_entries);
 791     snew(info->nkx, info->n_entries);
 792     snew(info->nky, info->n_entries);
 793     snew(info->nkz, info->n_entries);
 794     snew(info->fourier_sp, info->n_entries);
 795     snew(info->ewald_rtol, info->n_entries);
 796     snew(info->ewald_beta, info->n_entries);
 797     snew(info->pme_order, info->n_entries);
 798     snew(info->fn_out, info->n_entries);
 799     snew(info->e_dir, info->n_entries);
 800     snew(info->e_rec, info->n_entries);
 801 }
 802
 803
 804 /* Allocate and fill an array with coordinates and charges,
 805  * returns the number of charges found
 806  */
 807 static int prepare_x_q(real *q[], rvec *x[], gmx_mtop_t *mtop, rvec x_orig[], t_commrec *cr)
 808 {
 809     int                     i;
 810     int                     nq; /* number of charged particles */
 811     gmx_mtop_atomloop_all_t aloop;
 812     t_atom                 *atom;
 813
 814
 815     if (MASTER(cr))
 816     {
 817         snew(*q, mtop->natoms);
 818         snew(*x, mtop->natoms);
 819         nq = 0;
 820
 821         aloop = gmx_mtop_atomloop_all_init(mtop);
 822
 823         while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
 824         {
 825             if (is_charge(atom->q))
 826             {
 827                 (*q)[nq]     = atom->q;
 828                 (*x)[nq][XX] = x_orig[i][XX];
 829                 (*x)[nq][YY] = x_orig[i][YY];
 830                 (*x)[nq][ZZ] = x_orig[i][ZZ];
 831                 nq++;
 832             }
 833         }
 834         /* Give back some unneeded memory */
 835         srenew(*q, nq);
 836         srenew(*x, nq);
 837     }
 838     /* Broadcast x and q in the parallel case */
 839     if (PAR(cr))
 840     {
 841         /* Transfer the number of charges */
 842         block_bc(cr, nq);
 843         snew_bc(cr, *x, nq);
 844         snew_bc(cr, *q, nq);
 845         nblock_bc(cr, nq, *x);
 846         nblock_bc(cr, nq, *q);
 847     }
 848
 849     return nq;
 850 }
 851
 852
 853
 854 /* Read in the tpr file and save information we need later in info */
 855 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)
 856 {
 857     read_tpx_state(fn_sim_tpr, ir, state, NULL, mtop);
 858
 859     /* The values of the original tpr input file are save in the first
 860      * place [0] of the arrays */
 861     info->orig_sim_steps = ir->nsteps;
 862     info->pme_order[0]   = ir->pme_order;
 863     info->rcoulomb[0]    = ir->rcoulomb;
 864     info->rvdw[0]        = ir->rvdw;
 865     info->nkx[0]         = ir->nkx;
 866     info->nky[0]         = ir->nky;
 867     info->nkz[0]         = ir->nkz;
 868     info->ewald_rtol[0]  = ir->ewald_rtol;
 869     info->fracself       = fracself;
 870     if (user_beta > 0)
 871     {
 872         info->ewald_beta[0] = user_beta;
 873     }
 874     else
 875     {
 876         info->ewald_beta[0]  = calc_ewaldcoeff(info->rcoulomb[0], info->ewald_rtol[0]);
 877     }
 878
 879     /* Check if PME was chosen */
 880     if (EEL_PME(ir->coulombtype) == FALSE)
 881     {
 882         gmx_fatal(FARGS, "Can only do optimizations for simulations with PME");
 883     }
 884
 885     /* Check if rcoulomb == rlist, which is necessary for PME */
 886     if (!(ir->rcoulomb == ir->rlist))
 887     {
 888         gmx_fatal(FARGS, "PME requires rcoulomb (%f) to be equal to rlist (%f).", ir->rcoulomb, ir->rlist);
 889     }
 890 }
 891
 892
 893 /* Transfer what we need for parallelizing the reciprocal error estimate */
 894 static void bcast_info(t_inputinfo *info, t_commrec *cr)
 895 {
 896     nblock_bc(cr, info->n_entries, info->nkx);
 897     nblock_bc(cr, info->n_entries, info->nky);
 898     nblock_bc(cr, info->n_entries, info->nkz);
 899     nblock_bc(cr, info->n_entries, info->ewald_beta);
 900     nblock_bc(cr, info->n_entries, info->pme_order);
 901     nblock_bc(cr, info->n_entries, info->e_dir);
 902     nblock_bc(cr, info->n_entries, info->e_rec);
 903     block_bc(cr, info->volume);
 904     block_bc(cr, info->recipbox);
 905     block_bc(cr, info->natoms);
 906     block_bc(cr, info->fracself);
 907     block_bc(cr, info->bTUNE);
 908     block_bc(cr, info->q2all);
 909     block_bc(cr, info->q2allnr);
 910 }
 911
 912
 913 /* Estimate the error of the SPME Ewald sum. This estimate is based upon
 914  * a) a homogeneous distribution of the charges
 915  * b) a total charge of zero.
 916  */
 917 static void estimate_PME_error(t_inputinfo *info, t_state *state,
 918                                gmx_mtop_t *mtop, FILE *fp_out, gmx_bool bVerbose, unsigned int seed,
 919                                t_commrec *cr)
 920 {
 921     rvec *x     = NULL; /* The coordinates */
 922     real *q     = NULL; /* The charges     */
 923     real  edir  = 0.0;  /* real space error */
 924     real  erec  = 0.0;  /* reciprocal space error */
 925     real  derr  = 0.0;  /* difference of real and reciprocal space error */
 926     real  derr0 = 0.0;  /* difference of real and reciprocal space error */
 927     real  beta  = 0.0;  /* splitting parameter beta */
 928     real  beta0 = 0.0;  /* splitting parameter beta */
 929     int   ncharges;     /* The number of atoms with charges */
 930     int   nsamples;     /* The number of samples used for the calculation of the
 931                          * self-energy error term */
 932     int   i = 0;
 933
 934     if (MASTER(cr))
 935     {
 936         fprintf(fp_out, "\n--- PME ERROR ESTIMATE ---\n");
 937     }
 938
 939     /* Prepare an x and q array with only the charged atoms */
 940     ncharges = prepare_x_q(&q, &x, mtop, state->x, cr);
 941     if (MASTER(cr))
 942     {
 943         calc_q2all(mtop, &(info->q2all), &(info->q2allnr));
 944         info->ewald_rtol[0] = gmx_erfc(info->rcoulomb[0]*info->ewald_beta[0]);
 945         /* Write some info to log file */
 946         fprintf(fp_out, "Box volume              : %g nm^3\n", info->volume);
 947         fprintf(fp_out, "Number of charged atoms : %d (total atoms %d)\n", ncharges, info->natoms);
 948         fprintf(fp_out, "Coulomb radius          : %g nm\n", info->rcoulomb[0]);
 949         fprintf(fp_out, "Ewald_rtol              : %g\n", info->ewald_rtol[0]);
 950         fprintf(fp_out, "Ewald parameter beta    : %g\n", info->ewald_beta[0]);
 951         fprintf(fp_out, "Interpolation order     : %d\n", info->pme_order[0]);
 952         fprintf(fp_out, "Fourier grid (nx,ny,nz) : %d x %d x %d\n",
 953                 info->nkx[0], info->nky[0], info->nkz[0]);
 954         fflush(fp_out);
 955
 956     }
 957
 958     if (PAR(cr))
 959     {
 960         bcast_info(info, cr);
 961     }
 962
 963
 964     /* Calculate direct space error */
 965     info->e_dir[0] = estimate_direct(info);
 966
 967     /* Calculate reciprocal space error */
 968     info->e_rec[0] = estimate_reciprocal(info, x, q, ncharges, fp_out, bVerbose,
 969                                          seed, &nsamples, cr);
 970
 971     if (PAR(cr))
 972     {
 973         bcast_info(info, cr);
 974     }
 975
 976     if (MASTER(cr))
 977     {
 978         fprintf(fp_out, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
 979         fprintf(fp_out, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
 980         fprintf(fp_out, "Self-energy error term was estimated using %d samples\n", nsamples);
 981         fflush(fp_out);
 982         fprintf(stderr, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
 983         fprintf(stderr, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
 984     }
 985
 986     i = 0;
 987
 988     if (info->bTUNE)
 989     {
 990         if (MASTER(cr))
 991         {
 992             fprintf(stderr, "Starting tuning ...\n");
 993         }
 994         edir  = info->e_dir[0];
 995         erec  = info->e_rec[0];
 996         derr0 = edir-erec;
 997         beta0 = info->ewald_beta[0];
 998         if (derr > 0.0)
 999         {
1000             info->ewald_beta[0] += 0.1;
1001         }
1002         else
1003         {
1004             info->ewald_beta[0] -= 0.1;
1005         }
1006         info->e_dir[0] = estimate_direct(info);
1007         info->e_rec[0] = estimate_reciprocal(info, x, q, ncharges, fp_out, bVerbose,
1008                                              seed, &nsamples, cr);
1009
1010         if (PAR(cr))
1011         {
1012             bcast_info(info, cr);
1013         }
1014
1015
1016         edir = info->e_dir[0];
1017         erec = info->e_rec[0];
1018         derr = edir-erec;
1019         while (fabs(derr/std::min(erec, edir)) > 1e-4)
1020         {
1021
1022             beta                = info->ewald_beta[0];
1023             beta               -= derr*(info->ewald_beta[0]-beta0)/(derr-derr0);
1024             beta0               = info->ewald_beta[0];
1025             info->ewald_beta[0] = beta;
1026             derr0               = derr;
1027
1028             info->e_dir[0] = estimate_direct(info);
1029             info->e_rec[0] = estimate_reciprocal(info, x, q, ncharges, fp_out, bVerbose,
1030                                                  seed, &nsamples, cr);
1031
1032             if (PAR(cr))
1033             {
1034                 bcast_info(info, cr);
1035             }
1036
1037             edir = info->e_dir[0];
1038             erec = info->e_rec[0];
1039             derr = edir-erec;
1040
1041             if (MASTER(cr))
1042             {
1043                 i++;
1044                 fprintf(stderr, "difference between real and rec. space error (step %d): %g\n", i, fabs(derr));
1045                 fprintf(stderr, "old beta: %f\n", beta0);
1046                 fprintf(stderr, "new beta: %f\n", beta);
1047             }
1048         }
1049
1050         info->ewald_rtol[0] = gmx_erfc(info->rcoulomb[0]*info->ewald_beta[0]);
1051
1052         if (MASTER(cr))
1053         {
1054             /* Write some info to log file */
1055             fflush(fp_out);
1056             fprintf(fp_out, "=========  After tuning ========\n");
1057             fprintf(fp_out, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
1058             fprintf(fp_out, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
1059             fprintf(stderr, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
1060             fprintf(stderr, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
1061             fprintf(fp_out, "Ewald_rtol              : %g\n", info->ewald_rtol[0]);
1062             fprintf(fp_out, "Ewald parameter beta    : %g\n", info->ewald_beta[0]);
1063             fflush(fp_out);
1064
1065         }
1066
1067     }
1068
1069 }
1070
1071
1072 int gmx_pme_error(int argc, char *argv[])
1073 {
1074     const char     *desc[] = {
1075         "[TT]g_pme_error[tt] estimates the error of the electrostatic forces",
1076         "if using the sPME algorithm. The flag [TT]-tune[tt] will determine",
1077         "the splitting parameter such that the error is equally",
1078         "distributed over the real and reciprocal space part.",
1079         "The part of the error that stems from self interaction of the particles "
1080         "is computationally demanding. However, a good a approximation is to",
1081         "just use a fraction of the particles for this term which can be",
1082         "indicated by the flag [TT]-self[tt].[PAR]",
1083     };
1084
1085     real            fs        = 0.0; /* 0 indicates: not set by the user */
1086     real            user_beta = -1.0;
1087     real            fracself  = 1.0;
1088     t_inputinfo     info;
1089     t_state         state;     /* The state from the tpr input file */
1090     gmx_mtop_t      mtop;      /* The topology from the tpr input file */
1091     t_inputrec     *ir = NULL; /* The inputrec from the tpr file */
1092     FILE           *fp = NULL;
1093     t_commrec      *cr;
1094     unsigned long   PCA_Flags;
1095     gmx_bool        bTUNE    = FALSE;
1096     gmx_bool        bVerbose = FALSE;
1097     int             seed     = 0;
1098
1099
1100     static t_filenm fnm[] = {
1101         { efTPX, "-s",     NULL,    ffREAD },
1102         { efOUT, "-o",    "error",  ffWRITE },
1103         { efTPX, "-so",   "tuned",  ffOPTWR }
1104     };
1105
1106     output_env_t    oenv = NULL;
1107
1108     t_pargs         pa[] = {
1109         { "-beta",     FALSE, etREAL, {&user_beta},
1110           "If positive, overwrite ewald_beta from [TT].tpr[tt] file with this value" },
1111         { "-tune",     FALSE, etBOOL, {&bTUNE},
1112           "Tune the splitting parameter such that the error is equally distributed between real and reciprocal space" },
1113         { "-self",     FALSE, etREAL, {&fracself},
1114           "If between 0.0 and 1.0, determine self interaction error from just this fraction of the charged particles" },
1115         { "-seed",     FALSE, etINT,  {&seed},
1116           "Random number seed used for Monte Carlo algorithm when [TT]-self[tt] is set to a value between 0.0 and 1.0" },
1117         { "-v",        FALSE, etBOOL, {&bVerbose},
1118           "Be loud and noisy" }
1119     };
1120
1121
1122 #define NFILE asize(fnm)
1123
1124     cr = init_par();
1125 #ifdef GMX_LIB_MPI
1126     MPI_Barrier(MPI_COMM_WORLD);
1127 #endif
1128
1129     PCA_Flags  = PCA_NOEXIT_ON_ARGS;
1130     PCA_Flags |= (MASTER(cr) ? 0 : PCA_QUIET);
1131
1132     if (!parse_common_args(&argc, argv, PCA_Flags,
1133                            NFILE, fnm, asize(pa), pa, asize(desc), desc,
1134                            0, NULL, &oenv))
1135     {
1136         return 0;
1137     }
1138
1139     if (!bTUNE)
1140     {
1141         bTUNE = opt2bSet("-so", NFILE, fnm);
1142     }
1143
1144     info.n_entries = 1;
1145
1146     /* Allocate memory for the inputinfo struct: */
1147     create_info(&info);
1148     info.fourier_sp[0] = fs;
1149
1150     /* Read in the tpr file and open logfile for reading */
1151     if (MASTER(cr))
1152     {
1153         snew(ir, 1);
1154         read_tpr_file(opt2fn("-s", NFILE, fnm), &info, &state, &mtop, ir, user_beta, fracself);
1155
1156         fp = fopen(opt2fn("-o", NFILE, fnm), "w");
1157     }
1158
1159     /* Check consistency if the user provided fourierspacing */
1160     if (fs > 0 && MASTER(cr))
1161     {
1162         /* Recalculate the grid dimensions using fourierspacing from user input */
1163         info.nkx[0] = 0;
1164         info.nky[0] = 0;
1165         info.nkz[0] = 0;
1166         calc_grid(stdout, state.box, info.fourier_sp[0], &(info.nkx[0]), &(info.nky[0]), &(info.nkz[0]));
1167         if ( (ir->nkx != info.nkx[0]) || (ir->nky != info.nky[0]) || (ir->nkz != info.nkz[0]) )
1168         {
1169             gmx_fatal(FARGS, "Wrong fourierspacing %f nm, input file grid = %d x %d x %d, computed grid = %d x %d x %d",
1170                       fs, ir->nkx, ir->nky, ir->nkz, info.nkx[0], info.nky[0], info.nkz[0]);
1171         }
1172     }
1173
1174     /* Estimate (S)PME force error */
1175
1176     /* Determine the volume of the simulation box */
1177     if (MASTER(cr))
1178     {
1179         info.volume = det(state.box);
1180         calc_recipbox(state.box, info.recipbox);
1181         info.natoms = mtop.natoms;
1182         info.bTUNE  = bTUNE;
1183     }
1184
1185     if (PAR(cr))
1186     {
1187         bcast_info(&info, cr);
1188     }
1189
1190     /* Get an error estimate of the input tpr file and do some tuning if requested */
1191     estimate_PME_error(&info, &state, &mtop, fp, bVerbose, seed, cr);
1192
1193     if (MASTER(cr))
1194     {
1195         /* Write out optimized tpr file if requested */
1196         if (opt2bSet("-so", NFILE, fnm) || bTUNE)
1197         {
1198             ir->ewald_rtol = info.ewald_rtol[0];
1199             write_tpx_state(opt2fn("-so", NFILE, fnm), ir, &state, &mtop);
1200         }
1201         please_cite(fp, "Wang2010");
1202         fclose(fp);
1203     }
1204
1205     return 0;
1206 }