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