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 <algorithm>
  36
  37 #include "gromacs/commandline/pargs.h"
  38 #include "typedefs.h"
  39 #include "types/commrec.h"
  40 #include "gromacs/utility/smalloc.h"
  41 #include "vec.h"
  42 #include "copyrite.h"
  43 #include "gromacs/fileio/tpxio.h"
  44 #include "readinp.h"
  45 #include "calcgrid.h"
  46 #include "checkpoint.h"
  47 #include "gmx_ana.h"
  48 #include "gromacs/random/random.h"
  49 #include "physics.h"
  50 #include "mdatoms.h"
  51 #include "coulomb.h"
  52 #include "mtop_util.h"
  53 #include "network.h"
  54 #include "main.h"
  55 #include "macros.h"
  56
  57 #include "gromacs/legacyheaders/gmx_fatal.h"
  58
  59 /* We use the same defines as in mvdata.c here */
  60 #define  block_bc(cr,   d) gmx_bcast(     sizeof(d),     &(d), (cr))
  61 #define nblock_bc(cr, nr, d) gmx_bcast((nr)*sizeof((d)[0]), (d), (cr))
  62 #define   snew_bc(cr, d, nr) { if (!MASTER(cr)) {snew((d), (nr)); }}
  63 /* #define TAKETIME */
  64 /* #define DEBUG  */
  65 enum {
  66     ddnoSEL, ddnoINTERLEAVE, ddnoPP_PME, ddnoCARTESIAN, ddnoNR
  67 };
  68
  69 /* Enum for situations that can occur during log file parsing */
  70 enum {
  71     eParselogOK,
  72     eParselogNotFound,
  73     eParselogNoPerfData,
  74     eParselogTerm,
  75     eParselogResetProblem,
  76     eParselogNr
  77 };
  78
  79
  80 typedef struct
  81 {
  82     gmx_int64_t     orig_sim_steps;  /* Number of steps to be done in the real simulation  */
  83     int             n_entries;       /* Number of entries in arrays                        */
  84     real            volume;          /* The volume of the box                              */
  85     matrix          recipbox;        /* The reciprocal box                                 */
  86     int             natoms;          /* The number of atoms in the MD system               */
  87     real           *fac;             /* The scaling factor                                 */
  88     real           *rcoulomb;        /* The coulomb radii [0...nr_inputfiles]              */
  89     real           *rvdw;            /* The vdW radii                                      */
  90     int            *nkx, *nky, *nkz; /* Number of k vectors in each spatial dimension      */
  91     real           *fourier_sp;      /* Fourierspacing                                     */
  92     real           *ewald_rtol;      /* Real space tolerance for Ewald, determines         */
  93                                      /* the real/reciprocal space relative weight          */
  94     real           *ewald_beta;      /* Splitting parameter [1/nm]                         */
  95     real            fracself;        /* fraction of particles for SI error                 */
  96     real            q2all;           /* sum ( q ^2 )                                       */
  97     real            q2allnr;         /* nr of charges                                      */
  98     int            *pme_order;       /* Interpolation order for PME (bsplines)             */
  99     char          **fn_out;          /* Name of the output tpr file                        */
 100     real           *e_dir;           /* Direct space part of PME error with these settings */
 101     real           *e_rec;           /* Reciprocal space part of PME error                 */
 102     gmx_bool        bTUNE;           /* flag for tuning */
 103 } t_inputinfo;
 104
 105
 106 /* Returns TRUE when atom is charged */
 107 static gmx_bool is_charge(real charge)
 108 {
 109     if (charge*charge > GMX_REAL_EPS)
 110     {
 111         return TRUE;
 112     }
 113     else
 114     {
 115         return FALSE;
 116     }
 117 }
 118
 119
 120 /* calculate charge density */
 121 static void calc_q2all(
 122         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 += 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 += pow( tmp, -n );
 219     }
 220
 221     tmp   = m / K;
 222     tmp  *= 2.0*M_PI;
 223     denom = 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 += 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 += 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 += 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 * 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 * 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 += 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 * 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 * 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 += 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    = -sin(2.0 * M_PI * i * K * rcoord);
 376         tmp1   = 2.0 * M_PI * m / K + 2.0 * M_PI * i;
 377         tmp2   = pow(tmp1, -n);
 378         nom   += tmp * tmp2 * i;
 379         denom += tmp2;
 380     }
 381
 382     for (i = SUMORDER; i > 0; i--)
 383     {
 384         tmp    = -sin(2.0 * M_PI * i * K * rcoord);
 385         tmp1   = 2.0 * M_PI * m / K + 2.0 * M_PI * i;
 386         tmp2   = 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         unsigned 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     /* Random number generator */
 458     gmx_rng_t rng     = NULL;
 459     int      *numbers = NULL;
 460
 461     /* Index variables for parallel work distribution */
 462     int startglobal, stopglobal;
 463     int startlocal, stoplocal;
 464     int x_per_core;
 465     int xtot;
 466
 467 #ifdef TAKETIME
 468     double t0 = 0.0;
 469     double t1 = 0.0;
 470 #endif
 471
 472     rng = gmx_rng_init(seed);
 473
 474     clear_rvec(gridpx);
 475     clear_rvec(gridpxy);
 476     clear_rvec(gridp);
 477     clear_rvec(tmpvec);
 478     clear_rvec(tmpvec2);
 479
 480     for (i = 0; i < nr; i++)
 481     {
 482         q2_all += q[i]*q[i];
 483     }
 484
 485     /* Calculate indices for work distribution */
 486     startglobal = -info->nkx[0]/2;
 487     stopglobal  = info->nkx[0]/2;
 488     xtot        = stopglobal*2+1;
 489     if (PAR(cr))
 490     {
 491         x_per_core = static_cast<int>(ceil(static_cast<real>(xtot) / cr->nnodes));
 492         startlocal = startglobal + x_per_core*cr->nodeid;
 493         stoplocal  = startlocal + x_per_core -1;
 494         if (stoplocal > stopglobal)
 495         {
 496             stoplocal = stopglobal;
 497         }
 498     }
 499     else
 500     {
 501         startlocal = startglobal;
 502         stoplocal  = stopglobal;
 503         x_per_core = xtot;
 504     }
 505 /*
 506    #ifdef GMX_LIB_MPI
 507     MPI_Barrier(MPI_COMM_WORLD);
 508    #endif
 509  */
 510
 511 #ifdef 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  = 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         }
 621
 622     }
 623
 624     if (MASTER(cr))
 625     {
 626         fprintf(stderr, "\n");
 627     }
 628
 629     /* Use just a fraction of all charges to estimate the self energy error term? */
 630     bFraction =  (info->fracself > 0.0) && (info->fracself < 1.0);
 631
 632     if (bFraction)
 633     {
 634         /* Here xtot is the number of samples taken for the Monte Carlo calculation
 635          * of the average of term IV of equation 35 in Wang2010. Round up to a
 636          * number of samples that is divisible by the number of nodes */
 637         x_per_core  = static_cast<int>(ceil(info->fracself * nr / cr->nnodes));
 638         xtot        = x_per_core * cr->nnodes;
 639     }
 640     else
 641     {
 642         /* In this case we use all nr particle positions */
 643         xtot       = nr;
 644         x_per_core = static_cast<int>(ceil(static_cast<real>(xtot) / cr->nnodes));
 645     }
 646
 647     startlocal = x_per_core *  cr->nodeid;
 648     stoplocal  = std::min(startlocal + x_per_core, xtot);  /* min needed if xtot == nr */
 649
 650     if (bFraction)
 651     {
 652         /* Make shure we get identical results in serial and parallel. Therefore,
 653          * take the sample indices from a single, global random number array that
 654          * is constructed on the master node and that only depends on the seed */
 655         snew(numbers, xtot);
 656         if (MASTER(cr))
 657         {
 658             for (i = 0; i < xtot; i++)
 659             {
 660                 numbers[i] = static_cast<int>(floor(gmx_rng_uniform_real(rng) * nr));
 661             }
 662         }
 663         /* Broadcast the random number array to the other nodes */
 664         if (PAR(cr))
 665         {
 666             nblock_bc(cr, xtot, numbers);
 667         }
 668
 669         if (bVerbose && MASTER(cr))
 670         {
 671             fprintf(stdout, "Using %d sample%s to approximate the self interaction error term",
 672                     xtot, xtot == 1 ? "" : "s");
 673             if (PAR(cr))
 674             {
 675                 fprintf(stdout, " (%d sample%s per rank)", x_per_core, x_per_core == 1 ? "" : "s");
 676             }
 677             fprintf(stdout, ".\n");
 678         }
 679     }
 680
 681     /* Return the number of positions used for the Monte Carlo algorithm */
 682     *nsamples = xtot;
 683
 684     for (i = startlocal; i < stoplocal; i++)
 685     {
 686         e_rec3x = 0;
 687         e_rec3y = 0;
 688         e_rec3z = 0;
 689
 690         if (bFraction)
 691         {
 692             /* Randomly pick a charge */
 693             ci = numbers[i];
 694         }
 695         else
 696         {
 697             /* Use all charges */
 698             ci = i;
 699         }
 700
 701         /* for(nx=startlocal; nx<=stoplocal; nx++)*/
 702         for (nx = -info->nkx[0]/2; nx < info->nkx[0]/2+1; nx++)
 703         {
 704             svmul(nx, info->recipbox[XX], gridpx);
 705             for (ny = -info->nky[0]/2; ny < info->nky[0]/2+1; ny++)
 706             {
 707                 svmul(ny, info->recipbox[YY], tmpvec);
 708                 rvec_add(gridpx, tmpvec, gridpxy);
 709                 for (nz = -info->nkz[0]/2; nz < info->nkz[0]/2+1; nz++)
 710                 {
 711
 712                     if (0 == nx && 0 == ny && 0 == nz)
 713                     {
 714                         continue;
 715                     }
 716
 717                     svmul(nz, info->recipbox[ZZ], tmpvec);
 718                     rvec_add(gridpxy, tmpvec, gridp);
 719                     tmp      = norm2(gridp);
 720                     coeff    = exp(-1.0 * M_PI * M_PI * tmp / info->ewald_beta[0] / info->ewald_beta[0] );
 721                     coeff   /= tmp;
 722                     e_rec3x += coeff*eps_self(nx, info->nkx[0], info->recipbox[XX], info->pme_order[0], x[ci]);
 723                     e_rec3y += coeff*eps_self(ny, info->nky[0], info->recipbox[YY], info->pme_order[0], x[ci]);
 724                     e_rec3z += coeff*eps_self(nz, info->nkz[0], info->recipbox[ZZ], info->pme_order[0], x[ci]);
 725
 726                 }
 727             }
 728         }
 729
 730         clear_rvec(tmpvec2);
 731
 732         svmul(e_rec3x, info->recipbox[XX], tmpvec);
 733         rvec_inc(tmpvec2, tmpvec);
 734         svmul(e_rec3y, info->recipbox[YY], tmpvec);
 735         rvec_inc(tmpvec2, tmpvec);
 736         svmul(e_rec3z, info->recipbox[ZZ], tmpvec);
 737         rvec_inc(tmpvec2, tmpvec);
 738
 739         e_rec3 += q[ci]*q[ci]*q[ci]*q[ci]*norm2(tmpvec2) / ( xtot * M_PI * info->volume * M_PI * info->volume);
 740         if (MASTER(cr))
 741         {
 742             fprintf(stderr, "\rCalculating reciprocal error part 2 ... %3.0f%%",
 743                     100.0*(i+1)/stoplocal);
 744
 745         }
 746     }
 747
 748     if (MASTER(cr))
 749     {
 750         fprintf(stderr, "\n");
 751     }
 752
 753 #ifdef GMX_LIB_MPI
 754 #ifdef TAKETIME
 755     if (MASTER(cr))
 756     {
 757         t1 = MPI_Wtime() - t0;
 758         fprintf(fp_out, "Recip. err. est. took   : %lf s\n", t1);
 759     }
 760 #endif
 761 #endif
 762
 763 #ifdef DEBUG
 764     if (PAR(cr))
 765     {
 766         fprintf(stderr, "Rank %3d: nx=[%3d...%3d]  e_rec3=%e\n",
 767                 cr->nodeid, startlocal, stoplocal, e_rec3);
 768     }
 769 #endif
 770
 771     if (PAR(cr))
 772     {
 773         gmx_sum(1, &e_rec1, cr);
 774         gmx_sum(1, &e_rec2, cr);
 775         gmx_sum(1, &e_rec3, cr);
 776     }
 777
 778     /* e_rec1*=8.0 * q2_all / info->volume / info->volume / nr ;
 779        e_rec2*=  q2_all / M_PI / M_PI / info->volume / info->volume / nr ;
 780        e_rec3/= M_PI * M_PI * info->volume * info->volume * nr ;
 781      */
 782     e_rec = sqrt(e_rec1+e_rec2+e_rec3);
 783
 784
 785     return ONE_4PI_EPS0 * e_rec;
 786 }
 787
 788
 789 /* Allocate memory for the inputinfo struct: */
 790 static void create_info(t_inputinfo *info)
 791 {
 792     snew(info->fac, info->n_entries);
 793     snew(info->rcoulomb, info->n_entries);
 794     snew(info->rvdw, info->n_entries);
 795     snew(info->nkx, info->n_entries);
 796     snew(info->nky, info->n_entries);
 797     snew(info->nkz, info->n_entries);
 798     snew(info->fourier_sp, info->n_entries);
 799     snew(info->ewald_rtol, info->n_entries);
 800     snew(info->ewald_beta, info->n_entries);
 801     snew(info->pme_order, info->n_entries);
 802     snew(info->fn_out, info->n_entries);
 803     snew(info->e_dir, info->n_entries);
 804     snew(info->e_rec, info->n_entries);
 805 }
 806
 807
 808 /* Allocate and fill an array with coordinates and charges,
 809  * returns the number of charges found
 810  */
 811 static int prepare_x_q(real *q[], rvec *x[], gmx_mtop_t *mtop, rvec x_orig[], t_commrec *cr)
 812 {
 813     int                     i;
 814     int                     nq; /* number of charged particles */
 815     gmx_mtop_atomloop_all_t aloop;
 816     t_atom                 *atom;
 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
 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, NULL, 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, t_state *state,
 922                                gmx_mtop_t *mtop, FILE *fp_out, gmx_bool bVerbose, unsigned int seed,
 923                                t_commrec *cr)
 924 {
 925     rvec *x     = NULL; /* The coordinates */
 926     real *q     = NULL; /* 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, state->x, cr);
 945     if (MASTER(cr))
 946     {
 947         calc_q2all(mtop, &(info->q2all), &(info->q2allnr));
 948         info->ewald_rtol[0] = gmx_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 (fabs(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, fabs(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] = gmx_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 = NULL; /* The inputrec from the tpr file */
1096     FILE           *fp = NULL;
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         { efTPX, "-s",     NULL,    ffREAD },
1106         { efOUT, "-o",    "error",  ffWRITE },
1107         { efTPX, "-so",   "tuned",  ffOPTWR }
1108     };
1109
1110     output_env_t    oenv = NULL;
1111
1112     t_pargs         pa[] = {
1113         { "-beta",     FALSE, etREAL, {&user_beta},
1114           "If positive, overwrite ewald_beta from [TT].tpr[tt] 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 #ifdef GMX_LIB_MPI
1130     MPI_Barrier(MPI_COMM_WORLD);
1131 #endif
1132
1133     PCA_Flags  = PCA_NOEXIT_ON_ARGS;
1134     PCA_Flags |= (MASTER(cr) ? 0 : PCA_QUIET);
1135
1136     if (!parse_common_args(&argc, argv, PCA_Flags,
1137                            NFILE, fnm, asize(pa), pa, asize(desc), desc,
1138                            0, NULL, &oenv))
1139     {
1140         return 0;
1141     }
1142
1143     if (!bTUNE)
1144     {
1145         bTUNE = opt2bSet("-so", NFILE, fnm);
1146     }
1147
1148     info.n_entries = 1;
1149
1150     /* Allocate memory for the inputinfo struct: */
1151     create_info(&info);
1152     info.fourier_sp[0] = fs;
1153
1154     /* Read in the tpr file and open logfile for reading */
1155     if (MASTER(cr))
1156     {
1157         snew(ir, 1);
1158         read_tpr_file(opt2fn("-s", NFILE, fnm), &info, &state, &mtop, ir, user_beta, fracself);
1159
1160         fp = fopen(opt2fn("-o", NFILE, fnm), "w");
1161     }
1162
1163     /* Check consistency if the user provided fourierspacing */
1164     if (fs > 0 && MASTER(cr))
1165     {
1166         /* Recalculate the grid dimensions using fourierspacing from user input */
1167         info.nkx[0] = 0;
1168         info.nky[0] = 0;
1169         info.nkz[0] = 0;
1170         calc_grid(stdout, state.box, info.fourier_sp[0], &(info.nkx[0]), &(info.nky[0]), &(info.nkz[0]));
1171         if ( (ir->nkx != info.nkx[0]) || (ir->nky != info.nky[0]) || (ir->nkz != info.nkz[0]) )
1172         {
1173             gmx_fatal(FARGS, "Wrong fourierspacing %f nm, input file grid = %d x %d x %d, computed grid = %d x %d x %d",
1174                       fs, ir->nkx, ir->nky, ir->nkz, info.nkx[0], info.nky[0], info.nkz[0]);
1175         }
1176     }
1177
1178     /* Estimate (S)PME force error */
1179
1180     /* Determine the volume of the simulation box */
1181     if (MASTER(cr))
1182     {
1183         info.volume = det(state.box);
1184         calc_recipbox(state.box, info.recipbox);
1185         info.natoms = mtop.natoms;
1186         info.bTUNE  = bTUNE;
1187     }
1188
1189     if (PAR(cr))
1190     {
1191         bcast_info(&info, cr);
1192     }
1193
1194     /* Get an error estimate of the input tpr file and do some tuning if requested */
1195     estimate_PME_error(&info, &state, &mtop, fp, bVerbose, seed, cr);
1196
1197     if (MASTER(cr))
1198     {
1199         /* Write out optimized tpr file if requested */
1200         if (opt2bSet("-so", NFILE, fnm) || bTUNE)
1201         {
1202             ir->ewald_rtol = info.ewald_rtol[0];
1203             write_tpx_state(opt2fn("-so", NFILE, fnm), ir, &state, &mtop);
1204         }
1205         please_cite(fp, "Wang2010");
1206         fclose(fp);
1207     }
1208
1209     return 0;
1210 }