src/gromacs/gmxana/gmx_pme_error.cpp

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