src/gromacs/gmxana/gmx_pme_error.cpp

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