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