src/gromacs/gmxana/gmx_pme_error.cpp

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