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