src/gromacs/gmxana/gmx_pme_error.cpp

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