src/gromacs/tools/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.
   5  * Copyright (c) 2015,2016,2017,2018,2019 by the GROMACS development team.
   6  * Copyright (c) 2019,2020,2021, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include "pme_error.h"
  40
  41 #include "config.h"
  42
  43 #include <cmath>
  44
  45 #include <algorithm>
  46
  47 #include "gromacs/commandline/pargs.h"
  48 #include "gromacs/ewald/ewald_utils.h"
  49 #include "gromacs/ewald/pme.h"
  50 #include "gromacs/fft/calcgrid.h"
  51 #include "gromacs/fileio/checkpoint.h"
  52 #include "gromacs/fileio/tpxio.h"
  53 #include "gromacs/gmxlib/network.h"
  54 #include "gromacs/math/functions.h"
  55 #include "gromacs/math/units.h"
  56 #include "gromacs/math/vec.h"
  57 #include "gromacs/mdlib/broadcaststructs.h"
  58 #include "gromacs/mdtypes/commrec.h"
  59 #include "gromacs/mdtypes/inputrec.h"
  60 #include "gromacs/mdtypes/md_enums.h"
  61 #include "gromacs/mdtypes/state.h"
  62 #include "gromacs/random/threefry.h"
  63 #include "gromacs/random/uniformintdistribution.h"
  64 #include "gromacs/topology/mtop_util.h"
  65 #include "gromacs/topology/topology.h"
  66 #include "gromacs/utility/arraysize.h"
  67 #include "gromacs/utility/fatalerror.h"
  68 #include "gromacs/utility/pleasecite.h"
  69 #include "gromacs/utility/smalloc.h"
  70
  71 /* #define TAKETIME */
  72 /* #define DEBUG  */
  73
  74 /* Enum for situations that can occur during log file parsing */
  75 enum
  76 {
  77     eParselogOK,
  78     eParselogNotFound,
  79     eParselogNoPerfData,
  80     eParselogTerm,
  81     eParselogResetProblem,
  82     eParselogNr
  83 };
  84
  85
  86 typedef struct
  87 {
  88     int64_t orig_sim_steps;  /* Number of steps to be done in the real simulation  */
  89     int     n_entries;       /* Number of entries in arrays                        */
  90     real    volume;          /* The volume of the box                              */
  91     matrix  recipbox;        /* The reciprocal box                                 */
  92     int     natoms;          /* The number of atoms in the MD system               */
  93     real*   fac;             /* The scaling factor                                 */
  94     real*   rcoulomb;        /* The coulomb radii [0...nr_inputfiles]              */
  95     real*   rvdw;            /* The vdW radii                                      */
  96     int *   nkx, *nky, *nkz; /* Number of k vectors in each spatial dimension      */
  97     real*   fourier_sp;      /* Fourierspacing                                     */
  98     real*   ewald_rtol;      /* Real space tolerance for Ewald, determines         */
  99                              /* the real/reciprocal space relative weight          */
 100     real*    ewald_beta;     /* Splitting parameter [1/nm]                         */
 101     real     fracself;       /* fraction of particles for SI error                 */
 102     real     q2all;          /* sum ( q ^2 )                                       */
 103     real     q2allnr;        /* nr of charges                                      */
 104     int*     pme_order;      /* Interpolation order for PME (bsplines)             */
 105     char**   fn_out;         /* Name of the output tpr file                        */
 106     real*    e_dir;          /* Direct space part of PME error with these settings */
 107     real*    e_rec;          /* Reciprocal space part of PME error                 */
 108     gmx_bool bTUNE;          /* flag for tuning */
 109 } t_inputinfo;
 110
 111
 112 /* Returns TRUE when atom is charged */
 113 static gmx_bool is_charge(real charge)
 114 {
 115     return charge * charge > GMX_REAL_EPS;
 116 }
 117
 118
 119 /* calculate charge density */
 120 static void calc_q2all(const gmx_mtop_t* mtop, /* molecular topology */
 121                        real*             q2all,
 122                        real*             q2allnr)
 123 {
 124     real q2_all = 0; /* Sum of squared charges */
 125     int  nrq_mol;    /* Number of charges in a single molecule */
 126     int  nrq_all;    /* Total number of charges in the MD system */
 127     real qi, q2_mol;
 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 (const gmx_molblock_t& molblock : mtop->molblock) /* 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         const gmx_moltype_t& molecule = mtop->moltype[molblock.type];
 139         for (int i = 0; i < molecule.atoms.nr; i++)
 140         {
 141             qi = molecule.atoms.atom[i].q;
 142             /* Is this charge worth to be considered? */
 143             if (is_charge(qi))
 144             {
 145                 q2_mol += qi * qi;
 146                 nrq_mol++;
 147             }
 148         }
 149         /* Multiply with the number of molecules present of this type and add */
 150         q2_all += q2_mol * molblock.nmol;
 151         nrq_all += nrq_mol * molblock.nmol;
 152 #ifdef DEBUG
 153         fprintf(stderr,
 154                 "Molecule %2d (%5d atoms) q2_mol=%10.3e nr.mol.charges=%5d (%6dx)  q2_all=%10.3e  "
 155                 "tot.charges=%d\n",
 156                 imol,
 157                 molecule.atoms.nr,
 158                 q2_mol,
 159                 nrq_mol,
 160                 molblock.nmol,
 161                 q2_all,
 162                 nrq_all);
 163 #endif
 164     }
 165
 166     *q2all   = q2_all;
 167     *q2allnr = nrq_all;
 168 }
 169
 170
 171 /* Estimate the direct space part error of the SPME Ewald sum */
 172 static real estimate_direct(t_inputinfo* info)
 173 {
 174     real e_dir     = 0; /* Error estimate */
 175     real beta      = 0; /* Splitting parameter (1/nm) */
 176     real r_coulomb = 0; /* Cut-off in direct space */
 177
 178
 179     beta      = info->ewald_beta[0];
 180     r_coulomb = info->rcoulomb[0];
 181
 182     e_dir = 2.0 * info->q2all * gmx::invsqrt(info->q2allnr * r_coulomb * info->volume);
 183     e_dir *= std::exp(-beta * beta * r_coulomb * r_coulomb);
 184
 185     return gmx::c_one4PiEps0 * e_dir;
 186 }
 187
 188 #define SUMORDER 6
 189
 190 /* the following 4 functions determine polynomials required for the reciprocal error estimate */
 191
 192 static inline real eps_poly1(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 += std::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 += std::pow(tmp, -n);
 218     }
 219
 220     tmp = m / K;
 221     tmp *= 2.0 * M_PI;
 222     denom = std::pow(tmp, -n) + nom;
 223
 224     return -nom / denom;
 225 }
 226
 227 static inline real eps_poly2(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 static inline real eps_poly3(real m, /* grid coordinate in certain direction */
 266                              real K, /* grid size in corresponding direction */
 267                              real n) /* spline interpolation order of the SPME */
 268 {
 269     int  i;
 270     real nom   = 0; /* nominator */
 271     real denom = 0; /* denominator */
 272     real tmp   = 0;
 273
 274     if (m == 0.0)
 275     {
 276         return 0.0;
 277     }
 278
 279     for (i = -SUMORDER; i < 0; i++)
 280     {
 281         tmp = m / K + i;
 282         tmp *= 2.0 * M_PI;
 283         nom += i * std::pow(tmp, -2 * n);
 284     }
 285
 286     for (i = SUMORDER; i > 0; i--)
 287     {
 288         tmp = m / K + i;
 289         tmp *= 2.0 * M_PI;
 290         nom += i * std::pow(tmp, -2 * n);
 291     }
 292
 293     for (i = -SUMORDER; i < SUMORDER + 1; i++)
 294     {
 295         tmp = m / K + i;
 296         tmp *= 2.0 * M_PI;
 297         denom += std::pow(tmp, -n);
 298     }
 299
 300     return 2.0 * M_PI * nom / denom / denom;
 301 }
 302
 303 static inline real eps_poly4(real m, /* grid coordinate in certain direction */
 304                              real K, /* grid size in corresponding direction */
 305                              real n) /* spline interpolation order of the SPME */
 306 {
 307     int  i;
 308     real nom   = 0; /* nominator */
 309     real denom = 0; /* denominator */
 310     real tmp   = 0;
 311
 312     if (m == 0.0)
 313     {
 314         return 0.0;
 315     }
 316
 317     for (i = -SUMORDER; i < 0; i++)
 318     {
 319         tmp = m / K + i;
 320         tmp *= 2.0 * M_PI;
 321         nom += i * i * std::pow(tmp, -2 * n);
 322     }
 323
 324     for (i = SUMORDER; i > 0; i--)
 325     {
 326         tmp = m / K + i;
 327         tmp *= 2.0 * M_PI;
 328         nom += i * i * std::pow(tmp, -2 * n);
 329     }
 330
 331     for (i = -SUMORDER; i < SUMORDER + 1; i++)
 332     {
 333         tmp = m / K + i;
 334         tmp *= 2.0 * M_PI;
 335         denom += std::pow(tmp, -n);
 336     }
 337
 338     return 4.0 * M_PI * M_PI * nom / denom / denom;
 339 }
 340
 341 static inline real eps_self(real m,     /* grid coordinate in certain direction */
 342                             real K,     /* grid size in corresponding direction */
 343                             rvec rboxv, /* reciprocal box vector */
 344                             real n,     /* spline interpolation order of the SPME */
 345                             rvec x)     /* coordinate of charge */
 346 {
 347     int  i;
 348     real tmp    = 0; /* temporary variables for computations */
 349     real tmp1   = 0; /* temporary variables for computations */
 350     real tmp2   = 0; /* temporary variables for computations */
 351     real rcoord = 0; /* coordinate in certain reciprocal space direction */
 352     real nom    = 0; /* nominator */
 353     real denom  = 0; /* denominator */
 354
 355
 356     if (m == 0.0)
 357     {
 358         return 0.0;
 359     }
 360
 361     rcoord = iprod(rboxv, x);
 362
 363
 364     for (i = -SUMORDER; i < 0; i++)
 365     {
 366         tmp  = -std::sin(2.0 * M_PI * i * K * rcoord);
 367         tmp1 = 2.0 * M_PI * m / K + 2.0 * M_PI * i;
 368         tmp2 = std::pow(tmp1, -n);
 369         nom += tmp * tmp2 * i;
 370         denom += tmp2;
 371     }
 372
 373     for (i = SUMORDER; i > 0; i--)
 374     {
 375         tmp  = -std::sin(2.0 * M_PI * i * K * rcoord);
 376         tmp1 = 2.0 * M_PI * m / K + 2.0 * M_PI * i;
 377         tmp2 = std::pow(tmp1, -n);
 378         nom += tmp * tmp2 * i;
 379         denom += tmp2;
 380     }
 381
 382
 383     tmp  = 2.0 * M_PI * m / K;
 384     tmp1 = pow(tmp, -n);
 385     denom += tmp1;
 386
 387     return 2.0 * M_PI * nom / denom * K;
 388 }
 389
 390 #undef SUMORDER
 391
 392 /* The following routine is just a copy from pme.c */
 393
 394 static void calc_recipbox(matrix box, matrix recipbox)
 395 {
 396     /* Save some time by assuming upper right part is zero */
 397
 398     real tmp = 1.0 / (box[XX][XX] * box[YY][YY] * box[ZZ][ZZ]);
 399
 400     recipbox[XX][XX] = box[YY][YY] * box[ZZ][ZZ] * tmp;
 401     recipbox[XX][YY] = 0;
 402     recipbox[XX][ZZ] = 0;
 403     recipbox[YY][XX] = -box[YY][XX] * box[ZZ][ZZ] * tmp;
 404     recipbox[YY][YY] = box[XX][XX] * box[ZZ][ZZ] * tmp;
 405     recipbox[YY][ZZ] = 0;
 406     recipbox[ZZ][XX] = (box[YY][XX] * box[ZZ][YY] - box[YY][YY] * box[ZZ][XX]) * tmp;
 407     recipbox[ZZ][YY] = -box[ZZ][YY] * box[XX][XX] * tmp;
 408     recipbox[ZZ][ZZ] = box[XX][XX] * box[YY][YY] * tmp;
 409 }
 410
 411
 412 /* Estimate the reciprocal space part error of the SPME Ewald sum. */
 413 static real estimate_reciprocal(t_inputinfo* info,
 414                                 rvec         x[], /* array of particles */
 415                                 const real   q[], /* array of charges */
 416                                 int          nr,  /* number of charges = size of the charge array */
 417                                 FILE gmx_unused* fp_out,
 418                                 gmx_bool         bVerbose,
 419                                 int  seed,     /* The seed for the random number generator */
 420                                 int* nsamples, /* Return the number of samples used if Monte Carlo
 421                                                 * algorithm is used for self energy error estimate */
 422                                 t_commrec* cr)
 423 {
 424     real     e_rec   = 0; /* reciprocal error estimate */
 425     real     e_rec1  = 0; /* Error estimate term 1*/
 426     real     e_rec2  = 0; /* Error estimate term 2*/
 427     real     e_rec3  = 0; /* Error estimate term 3 */
 428     real     e_rec3x = 0; /* part of Error estimate term 3 in x */
 429     real     e_rec3y = 0; /* part of Error estimate term 3 in y */
 430     real     e_rec3z = 0; /* part of Error estimate term 3 in z */
 431     int      i, ci;
 432     int      nx, ny, nz; /* grid coordinates */
 433     real     q2_all = 0; /* sum of squared charges */
 434     rvec     gridpx;     /* reciprocal grid point in x direction*/
 435     rvec     gridpxy;    /* reciprocal grid point in x and y direction*/
 436     rvec     gridp;      /* complete reciprocal grid point in 3 directions*/
 437     rvec     tmpvec;     /* template to create points from basis vectors */
 438     rvec     tmpvec2;    /* template to create points from basis vectors */
 439     real     coeff  = 0; /* variable to compute coefficients of the error estimate */
 440     real     coeff2 = 0; /* variable to compute coefficients of the error estimate */
 441     real     tmp    = 0; /* variables to compute different factors from vectors */
 442     real     tmp1   = 0;
 443     real     tmp2   = 0;
 444     gmx_bool bFraction;
 445
 446     int* numbers = nullptr;
 447
 448     /* Index variables for parallel work distribution */
 449     int startglobal, stopglobal;
 450     int startlocal, stoplocal;
 451     int x_per_core;
 452     int xtot;
 453
 454 #ifdef TAKETIME
 455     double t0 = 0.0;
 456     double t1 = 0.0;
 457 #endif
 458
 459     GMX_RELEASE_ASSERT(q != nullptr, "Must have charges");
 460
 461     if (seed == 0)
 462     {
 463         seed = static_cast<int>(gmx::makeRandomSeed());
 464     }
 465     fprintf(stderr, "Using random seed %d.\n", seed);
 466
 467     gmx::DefaultRandomEngine         rng(seed);
 468     gmx::UniformIntDistribution<int> dist(0, nr - 1);
 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>(std::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 #if GMX_LIB_MPI
 503 #    ifdef TAKETIME
 504     if (MASTER(cr))
 505     {
 506         t0 = MPI_Wtime();
 507     }
 508 #    endif
 509 #endif
 510
 511     if (MASTER(cr))
 512     {
 513
 514         fprintf(stderr, "Calculating reciprocal error part 1 ...");
 515     }
 516
 517     for (nx = startlocal; nx <= stoplocal; nx++)
 518     {
 519         svmul(nx, info->recipbox[XX], gridpx);
 520         for (ny = -info->nky[0] / 2; ny < info->nky[0] / 2 + 1; ny++)
 521         {
 522             svmul(ny, info->recipbox[YY], tmpvec);
 523             rvec_add(gridpx, tmpvec, gridpxy);
 524             for (nz = -info->nkz[0] / 2; nz < info->nkz[0] / 2 + 1; nz++)
 525             {
 526                 if (0 == nx && 0 == ny && 0 == nz)
 527                 {
 528                     continue;
 529                 }
 530                 svmul(nz, info->recipbox[ZZ], tmpvec);
 531                 rvec_add(gridpxy, tmpvec, gridp);
 532                 tmp = norm2(gridp);
 533                 coeff = std::exp(-1.0 * M_PI * M_PI * tmp / info->ewald_beta[0] / info->ewald_beta[0]);
 534                 coeff /= 2.0 * M_PI * info->volume * tmp;
 535                 coeff2 = tmp;
 536
 537
 538                 tmp = eps_poly2(nx, info->nkx[0], info->pme_order[0]);
 539                 tmp += eps_poly2(ny, info->nkx[0], info->pme_order[0]);
 540                 tmp += eps_poly2(nz, info->nkx[0], info->pme_order[0]);
 541
 542                 tmp1 = eps_poly1(nx, info->nkx[0], info->pme_order[0]);
 543                 tmp2 = eps_poly1(ny, info->nky[0], info->pme_order[0]);
 544
 545                 tmp += 2.0 * tmp1 * tmp2;
 546
 547                 tmp1 = eps_poly1(nz, info->nkz[0], info->pme_order[0]);
 548                 tmp2 = eps_poly1(ny, info->nky[0], info->pme_order[0]);
 549
 550                 tmp += 2.0 * tmp1 * tmp2;
 551
 552                 tmp1 = eps_poly1(nz, info->nkz[0], info->pme_order[0]);
 553                 tmp2 = eps_poly1(nx, info->nkx[0], info->pme_order[0]);
 554
 555                 tmp += 2.0 * tmp1 * tmp2;
 556
 557                 tmp1 = eps_poly1(nx, info->nkx[0], info->pme_order[0]);
 558                 tmp1 += eps_poly1(ny, info->nky[0], info->pme_order[0]);
 559                 tmp1 += eps_poly1(nz, info->nkz[0], info->pme_order[0]);
 560
 561                 tmp += tmp1 * tmp1;
 562
 563                 e_rec1 += 32.0 * M_PI * M_PI * coeff * coeff * coeff2 * tmp * q2_all * q2_all / nr;
 564
 565                 tmp1 = eps_poly3(nx, info->nkx[0], info->pme_order[0]);
 566                 tmp1 *= info->nkx[0];
 567                 tmp2 = iprod(gridp, info->recipbox[XX]);
 568
 569                 tmp = tmp1 * tmp2;
 570
 571                 tmp1 = eps_poly3(ny, info->nky[0], info->pme_order[0]);
 572                 tmp1 *= info->nky[0];
 573                 tmp2 = iprod(gridp, info->recipbox[YY]);
 574
 575                 tmp += tmp1 * tmp2;
 576
 577                 tmp1 = eps_poly3(nz, info->nkz[0], info->pme_order[0]);
 578                 tmp1 *= info->nkz[0];
 579                 tmp2 = iprod(gridp, info->recipbox[ZZ]);
 580
 581                 tmp += tmp1 * tmp2;
 582
 583                 tmp *= 4.0 * M_PI;
 584
 585                 tmp1 = eps_poly4(nx, info->nkx[0], info->pme_order[0]);
 586                 tmp1 *= norm2(info->recipbox[XX]);
 587                 tmp1 *= info->nkx[0] * info->nkx[0];
 588
 589                 tmp += tmp1;
 590
 591                 tmp1 = eps_poly4(ny, info->nky[0], info->pme_order[0]);
 592                 tmp1 *= norm2(info->recipbox[YY]);
 593                 tmp1 *= info->nky[0] * info->nky[0];
 594
 595                 tmp += tmp1;
 596
 597                 tmp1 = eps_poly4(nz, info->nkz[0], info->pme_order[0]);
 598                 tmp1 *= norm2(info->recipbox[ZZ]);
 599                 tmp1 *= info->nkz[0] * info->nkz[0];
 600
 601                 tmp += tmp1;
 602
 603                 e_rec2 += 4.0 * coeff * coeff * tmp * q2_all * q2_all / nr;
 604             }
 605         }
 606         if (MASTER(cr))
 607         {
 608             fprintf(stderr,
 609                     "\rCalculating reciprocal error part 1 ... %3.0f%%",
 610                     100.0 * (nx - startlocal + 1) / (x_per_core));
 611             fflush(stderr);
 612         }
 613     }
 614
 615     if (MASTER(cr))
 616     {
 617         fprintf(stderr, "\n");
 618     }
 619
 620     /* Use just a fraction of all charges to estimate the self energy error term? */
 621     bFraction = (info->fracself > 0.0) && (info->fracself < 1.0);
 622
 623     if (bFraction)
 624     {
 625         /* Here xtot is the number of samples taken for the Monte Carlo calculation
 626          * of the average of term IV of equation 35 in Wang2010. Round up to a
 627          * number of samples that is divisible by the number of nodes */
 628         x_per_core = static_cast<int>(std::ceil(info->fracself * nr / cr->nnodes));
 629         xtot       = x_per_core * cr->nnodes;
 630     }
 631     else
 632     {
 633         /* In this case we use all nr particle positions */
 634         xtot       = nr;
 635         x_per_core = static_cast<int>(std::ceil(static_cast<real>(xtot) / cr->nnodes));
 636     }
 637
 638     startlocal = x_per_core * cr->nodeid;
 639     stoplocal  = std::min(startlocal + x_per_core, xtot); /* min needed if xtot == nr */
 640
 641     if (bFraction)
 642     {
 643         /* Make shure we get identical results in serial and parallel. Therefore,
 644          * take the sample indices from a single, global random number array that
 645          * is constructed on the master node and that only depends on the seed */
 646         snew(numbers, xtot);
 647         if (MASTER(cr))
 648         {
 649             for (i = 0; i < xtot; i++)
 650             {
 651                 numbers[i] = dist(rng); // [0,nr-1]
 652             }
 653         }
 654         /* Broadcast the random number array to the other nodes */
 655         if (PAR(cr))
 656         {
 657             nblock_bc(cr->mpi_comm_mygroup, xtot, numbers);
 658         }
 659
 660         if (bVerbose && MASTER(cr))
 661         {
 662             fprintf(stdout,
 663                     "Using %d sample%s to approximate the self interaction error term",
 664                     xtot,
 665                     xtot == 1 ? "" : "s");
 666             if (PAR(cr))
 667             {
 668                 fprintf(stdout, " (%d sample%s per rank)", x_per_core, x_per_core == 1 ? "" : "s");
 669             }
 670             fprintf(stdout, ".\n");
 671         }
 672     }
 673
 674     /* Return the number of positions used for the Monte Carlo algorithm */
 675     *nsamples = xtot;
 676
 677     for (i = startlocal; i < stoplocal; i++)
 678     {
 679         e_rec3x = 0;
 680         e_rec3y = 0;
 681         e_rec3z = 0;
 682
 683         if (bFraction)
 684         {
 685             /* Randomly pick a charge */
 686             ci = numbers[i];
 687         }
 688         else
 689         {
 690             /* Use all charges */
 691             ci = i;
 692         }
 693
 694         /* for(nx=startlocal; nx<=stoplocal; nx++)*/
 695         for (nx = -info->nkx[0] / 2; nx < info->nkx[0] / 2 + 1; nx++)
 696         {
 697             svmul(nx, info->recipbox[XX], gridpx);
 698             for (ny = -info->nky[0] / 2; ny < info->nky[0] / 2 + 1; ny++)
 699             {
 700                 svmul(ny, info->recipbox[YY], tmpvec);
 701                 rvec_add(gridpx, tmpvec, gridpxy);
 702                 for (nz = -info->nkz[0] / 2; nz < info->nkz[0] / 2 + 1; nz++)
 703                 {
 704
 705                     if (0 == nx && 0 == ny && 0 == nz)
 706                     {
 707                         continue;
 708                     }
 709
 710                     svmul(nz, info->recipbox[ZZ], tmpvec);
 711                     rvec_add(gridpxy, tmpvec, gridp);
 712                     tmp = norm2(gridp);
 713                     coeff = std::exp(-1.0 * M_PI * M_PI * tmp / info->ewald_beta[0] / info->ewald_beta[0]);
 714                     coeff /= tmp;
 715                     e_rec3x += coeff
 716                                * eps_self(nx, info->nkx[0], info->recipbox[XX], info->pme_order[0], x[ci]);
 717                     e_rec3y += coeff
 718                                * eps_self(ny, info->nky[0], info->recipbox[YY], info->pme_order[0], x[ci]);
 719                     e_rec3z += coeff
 720                                * eps_self(nz, info->nkz[0], info->recipbox[ZZ], info->pme_order[0], x[ci]);
 721                 }
 722             }
 723         }
 724
 725         clear_rvec(tmpvec2);
 726
 727         svmul(e_rec3x, info->recipbox[XX], tmpvec);
 728         rvec_inc(tmpvec2, tmpvec);
 729         svmul(e_rec3y, info->recipbox[YY], tmpvec);
 730         rvec_inc(tmpvec2, tmpvec);
 731         svmul(e_rec3z, info->recipbox[ZZ], tmpvec);
 732         rvec_inc(tmpvec2, tmpvec);
 733
 734         e_rec3 += q[ci] * q[ci] * q[ci] * q[ci] * norm2(tmpvec2)
 735                   / (xtot * M_PI * info->volume * M_PI * info->volume);
 736         if (MASTER(cr))
 737         {
 738             fprintf(stderr, "\rCalculating reciprocal error part 2 ... %3.0f%%", 100.0 * (i + 1) / stoplocal);
 739             fflush(stderr);
 740         }
 741     }
 742
 743     if (MASTER(cr))
 744     {
 745         fprintf(stderr, "\n");
 746     }
 747
 748 #if GMX_LIB_MPI
 749 #    ifdef TAKETIME
 750     if (MASTER(cr))
 751     {
 752         t1 = MPI_Wtime() - t0;
 753         fprintf(fp_out, "Recip. err. est. took   : %lf s\n", t1);
 754     }
 755 #    endif
 756 #endif
 757
 758 #ifdef DEBUG
 759     if (PAR(cr))
 760     {
 761         fprintf(stderr, "Rank %3d: nx=[%3d...%3d]  e_rec3=%e\n", cr->nodeid, startlocal, stoplocal, e_rec3);
 762     }
 763 #endif
 764
 765     if (PAR(cr))
 766     {
 767         gmx_sum(1, &e_rec1, cr);
 768         gmx_sum(1, &e_rec2, cr);
 769         gmx_sum(1, &e_rec3, cr);
 770     }
 771
 772     /* e_rec1*=8.0 * q2_all / info->volume / info->volume / nr ;
 773        e_rec2*=  q2_all / M_PI / M_PI / info->volume / info->volume / nr ;
 774        e_rec3/= M_PI * M_PI * info->volume * info->volume * nr ;
 775      */
 776     e_rec = std::sqrt(e_rec1 + e_rec2 + e_rec3);
 777
 778
 779     return gmx::c_one4PiEps0 * e_rec;
 780 }
 781
 782
 783 /* Allocate memory for the inputinfo struct: */
 784 static void create_info(t_inputinfo* info)
 785 {
 786     snew(info->fac, info->n_entries);
 787     snew(info->rcoulomb, info->n_entries);
 788     snew(info->rvdw, info->n_entries);
 789     snew(info->nkx, info->n_entries);
 790     snew(info->nky, info->n_entries);
 791     snew(info->nkz, info->n_entries);
 792     snew(info->fourier_sp, info->n_entries);
 793     snew(info->ewald_rtol, info->n_entries);
 794     snew(info->ewald_beta, info->n_entries);
 795     snew(info->pme_order, info->n_entries);
 796     snew(info->fn_out, info->n_entries);
 797     snew(info->e_dir, info->n_entries);
 798     snew(info->e_rec, info->n_entries);
 799 }
 800
 801
 802 /* Allocate and fill an array with coordinates and charges,
 803  * returns the number of charges found
 804  */
 805 static int prepare_x_q(real* q[], rvec* x[], const gmx_mtop_t* mtop, const rvec x_orig[], t_commrec* cr)
 806 {
 807     int nq; /* number of charged particles */
 808
 809
 810     if (MASTER(cr))
 811     {
 812         snew(*q, mtop->natoms);
 813         snew(*x, mtop->natoms);
 814         nq = 0;
 815
 816         for (const AtomProxy atomP : AtomRange(*mtop))
 817         {
 818             const t_atom& local = atomP.atom();
 819             int           i     = atomP.globalAtomNumber();
 820             if (is_charge(local.q))
 821             {
 822                 (*q)[nq]     = local.q;
 823                 (*x)[nq][XX] = x_orig[i][XX];
 824                 (*x)[nq][YY] = x_orig[i][YY];
 825                 (*x)[nq][ZZ] = x_orig[i][ZZ];
 826                 nq++;
 827             }
 828         }
 829         /* Give back some unneeded memory */
 830         srenew(*q, nq);
 831         srenew(*x, nq);
 832     }
 833     /* Broadcast x and q in the parallel case */
 834     if (PAR(cr))
 835     {
 836         /* Transfer the number of charges */
 837         block_bc(cr->mpi_comm_mygroup, nq);
 838         snew_bc(MASTER(cr), *x, nq);
 839         snew_bc(MASTER(cr), *q, nq);
 840         nblock_bc(cr->mpi_comm_mygroup, nq, *x);
 841         nblock_bc(cr->mpi_comm_mygroup, nq, *q);
 842     }
 843
 844     return nq;
 845 }
 846
 847
 848 /* Read in the tpr file and save information we need later in info */
 849 static void read_tpr_file(const char*  fn_sim_tpr,
 850                           t_inputinfo* info,
 851                           t_state*     state,
 852                           gmx_mtop_t*  mtop,
 853                           t_inputrec*  ir,
 854                           real         user_beta,
 855                           real         fracself)
 856 {
 857     read_tpx_state(fn_sim_tpr, ir, state, 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_q(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, const t_commrec* cr)
 895 {
 896     nblock_bc(cr->mpi_comm_mygroup, info->n_entries, info->nkx);
 897     nblock_bc(cr->mpi_comm_mygroup, info->n_entries, info->nky);
 898     nblock_bc(cr->mpi_comm_mygroup, info->n_entries, info->nkz);
 899     nblock_bc(cr->mpi_comm_mygroup, info->n_entries, info->ewald_beta);
 900     nblock_bc(cr->mpi_comm_mygroup, info->n_entries, info->pme_order);
 901     nblock_bc(cr->mpi_comm_mygroup, info->n_entries, info->e_dir);
 902     nblock_bc(cr->mpi_comm_mygroup, info->n_entries, info->e_rec);
 903     block_bc(cr->mpi_comm_mygroup, info->volume);
 904     block_bc(cr->mpi_comm_mygroup, info->recipbox);
 905     block_bc(cr->mpi_comm_mygroup, info->natoms);
 906     block_bc(cr->mpi_comm_mygroup, info->fracself);
 907     block_bc(cr->mpi_comm_mygroup, info->bTUNE);
 908     block_bc(cr->mpi_comm_mygroup, info->q2all);
 909     block_bc(cr->mpi_comm_mygroup, 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,
 918                                const t_state*    state,
 919                                const gmx_mtop_t* mtop,
 920                                FILE*             fp_out,
 921                                gmx_bool          bVerbose,
 922                                unsigned int      seed,
 923                                t_commrec*        cr)
 924 {
 925     rvec* x     = nullptr; /* The coordinates */
 926     real* q     = nullptr; /* The charges     */
 927     real  edir  = 0.0;     /* real space error */
 928     real  erec  = 0.0;     /* reciprocal space error */
 929     real  derr  = 0.0;     /* difference of real and reciprocal space error */
 930     real  derr0 = 0.0;     /* difference of real and reciprocal space error */
 931     real  beta  = 0.0;     /* splitting parameter beta */
 932     real  beta0 = 0.0;     /* splitting parameter beta */
 933     int   ncharges;        /* The number of atoms with charges */
 934     int   nsamples;        /* The number of samples used for the calculation of the
 935                             * self-energy error term */
 936     int i = 0;
 937
 938     if (MASTER(cr))
 939     {
 940         fprintf(fp_out, "\n--- PME ERROR ESTIMATE ---\n");
 941     }
 942
 943     /* Prepare an x and q array with only the charged atoms */
 944     ncharges = prepare_x_q(&q, &x, mtop, state->x.rvec_array(), cr);
 945     if (MASTER(cr))
 946     {
 947         calc_q2all(mtop, &(info->q2all), &(info->q2allnr));
 948         info->ewald_rtol[0] = std::erfc(info->rcoulomb[0] * info->ewald_beta[0]);
 949         /* Write some info to log file */
 950         fprintf(fp_out, "Box volume              : %g nm^3\n", info->volume);
 951         fprintf(fp_out, "Number of charged atoms : %d (total atoms %d)\n", ncharges, info->natoms);
 952         fprintf(fp_out, "Coulomb radius          : %g nm\n", info->rcoulomb[0]);
 953         fprintf(fp_out, "Ewald_rtol              : %g\n", info->ewald_rtol[0]);
 954         fprintf(fp_out, "Ewald parameter beta    : %g\n", info->ewald_beta[0]);
 955         fprintf(fp_out, "Interpolation order     : %d\n", info->pme_order[0]);
 956         fprintf(fp_out, "Fourier grid (nx,ny,nz) : %d x %d x %d\n", info->nkx[0], info->nky[0], info->nkz[0]);
 957         fflush(fp_out);
 958     }
 959
 960     if (PAR(cr))
 961     {
 962         bcast_info(info, cr);
 963     }
 964
 965
 966     /* Calculate direct space error */
 967     info->e_dir[0] = estimate_direct(info);
 968
 969     /* Calculate reciprocal space error */
 970     info->e_rec[0] = estimate_reciprocal(info, x, q, ncharges, fp_out, bVerbose, 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, 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 (std::abs(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] =
1030                     estimate_reciprocal(info, x, q, ncharges, fp_out, bVerbose, 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,
1045                         "difference between real and rec. space error (step %d): %g\n",
1046                         i,
1047                         std::abs(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] = std::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 int gmx_pme_error(int argc, char* argv[])
1073 {
1074     const char* desc[] = {
1075         "[THISMODULE] 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     FILE*         fp = nullptr;
1092     unsigned long PCA_Flags;
1093     gmx_bool      bTUNE    = FALSE;
1094     gmx_bool      bVerbose = FALSE;
1095     int           seed     = 0;
1096
1097
1098     static t_filenm fnm[] = { { efTPR, "-s", nullptr, ffREAD },
1099                               { efOUT, "-o", "error", ffWRITE },
1100                               { efTPR, "-so", "tuned", ffOPTWR } };
1101
1102     gmx_output_env_t* oenv = nullptr;
1103
1104     t_pargs pa[] = {
1105         { "-beta",
1106           FALSE,
1107           etREAL,
1108           { &user_beta },
1109           "If positive, overwrite ewald_beta from [REF].tpr[ref] file with this value" },
1110         { "-tune",
1111           FALSE,
1112           etBOOL,
1113           { &bTUNE },
1114           "Tune the splitting parameter such that the error is equally distributed between "
1115           "real and reciprocal space" },
1116         { "-self",
1117           FALSE,
1118           etREAL,
1119           { &fracself },
1120           "If between 0.0 and 1.0, determine self interaction error from just this "
1121           "fraction of the charged particles" },
1122         { "-seed",
1123           FALSE,
1124           etINT,
1125           { &seed },
1126           "Random number seed used for Monte Carlo algorithm when [TT]-self[tt] is set to "
1127           "a value between 0.0 and 1.0" },
1128         { "-v", FALSE, etBOOL, { &bVerbose }, "Be loud and noisy" }
1129     };
1130
1131
1132 #define NFILE asize(fnm)
1133
1134     CommrecHandle commrecHandle = init_commrec(MPI_COMM_WORLD);
1135     t_commrec*    cr            = commrecHandle.get();
1136     PCA_Flags                   = PCA_NOEXIT_ON_ARGS;
1137
1138     if (!parse_common_args(
1139                 &argc, argv, PCA_Flags, NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
1140     {
1141         return 0;
1142     }
1143
1144     if (!bTUNE)
1145     {
1146         bTUNE = opt2bSet("-so", NFILE, fnm);
1147     }
1148
1149     info.n_entries = 1;
1150
1151     /* Allocate memory for the inputinfo struct: */
1152     create_info(&info);
1153     info.fourier_sp[0] = fs;
1154
1155     t_inputrec ir;
1156     if (MASTER(cr))
1157     {
1158         read_tpr_file(opt2fn("-s", NFILE, fnm), &info, &state, &mtop, &ir, user_beta, fracself);
1159         /* Open logfile for reading */
1160         fp = fopen(opt2fn("-o", NFILE, fnm), "w");
1161
1162         /* Determine the volume of the simulation box */
1163         info.volume = det(state.box);
1164         calc_recipbox(state.box, info.recipbox);
1165         info.natoms = mtop.natoms;
1166         info.bTUNE  = bTUNE;
1167     }
1168
1169     /* Check consistency if the user provided fourierspacing */
1170     if (fs > 0 && MASTER(cr))
1171     {
1172         /* Recalculate the grid dimensions using fourierspacing from user input */
1173         info.nkx[0] = 0;
1174         info.nky[0] = 0;
1175         info.nkz[0] = 0;
1176         calcFftGrid(stdout,
1177                     state.box,
1178                     info.fourier_sp[0],
1179                     minimalPmeGridSize(info.pme_order[0]),
1180                     &(info.nkx[0]),
1181                     &(info.nky[0]),
1182                     &(info.nkz[0]));
1183         if ((ir.nkx != info.nkx[0]) || (ir.nky != info.nky[0]) || (ir.nkz != info.nkz[0]))
1184         {
1185             gmx_fatal(FARGS,
1186                       "Wrong fourierspacing %f nm, input file grid = %d x %d x %d, computed grid = "
1187                       "%d x %d x %d",
1188                       fs,
1189                       ir.nkx,
1190                       ir.nky,
1191                       ir.nkz,
1192                       info.nkx[0],
1193                       info.nky[0],
1194                       info.nkz[0]);
1195         }
1196     }
1197
1198     /* Estimate (S)PME force error */
1199
1200     if (PAR(cr))
1201     {
1202         bcast_info(&info, cr);
1203     }
1204
1205     /* Get an error estimate of the input tpr file and do some tuning if requested */
1206     estimate_PME_error(&info, &state, &mtop, fp, bVerbose, seed, cr);
1207
1208     if (MASTER(cr))
1209     {
1210         /* Write out optimized tpr file if requested */
1211         if (opt2bSet("-so", NFILE, fnm) || bTUNE)
1212         {
1213             ir.ewald_rtol = info.ewald_rtol[0];
1214             write_tpx_state(opt2fn("-so", NFILE, fnm), &ir, &state, mtop);
1215         }
1216         please_cite(fp, "Wang2010");
1217         fclose(fp);
1218     }
1219
1220     return 0;
1221 }