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