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