src/gromacs/gmxana/gmx_pme_error.c

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