src/gromacs/gmxana/gmx_current.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2008,2009,2010,2011,2012 by the GROMACS development team.
   5  * Copyright (c) 2013,2014,2015,2017,2018 by the GROMACS development team.
   6  * Copyright (c) 2019,2020,2021, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include <cmath>
  40 #include <cstdlib>
  41
  42 #include "gromacs/commandline/pargs.h"
  43 #include "gromacs/fileio/confio.h"
  44 #include "gromacs/fileio/trxio.h"
  45 #include "gromacs/fileio/xvgr.h"
  46 #include "gromacs/gmxana/gmx_ana.h"
  47 #include "gromacs/math/units.h"
  48 #include "gromacs/math/vec.h"
  49 #include "gromacs/pbcutil/pbc.h"
  50 #include "gromacs/pbcutil/rmpbc.h"
  51 #include "gromacs/statistics/statistics.h"
  52 #include "gromacs/topology/index.h"
  53 #include "gromacs/topology/topology.h"
  54 #include "gromacs/trajectory/trajectoryframe.h"
  55 #include "gromacs/utility/arraysize.h"
  56 #include "gromacs/utility/cstringutil.h"
  57 #include "gromacs/utility/fatalerror.h"
  58 #include "gromacs/utility/gmxassert.h"
  59 #include "gromacs/utility/smalloc.h"
  60
  61 constexpr double EPSI0 = (gmx::c_epsilon0 * gmx::c_electronCharge * gmx::c_electronCharge * gmx::c_avogadro
  62                           / (gmx::c_kilo * gmx::c_nano)); /* c_epsilon0 in SI units */
  63
  64 static void index_atom2mol(int* n, int* index, t_block* mols)
  65 {
  66     int nat, i, nmol, mol, j;
  67
  68     nat  = *n;
  69     i    = 0;
  70     nmol = 0;
  71     mol  = 0;
  72     while (i < nat)
  73     {
  74         while (index[i] > mols->index[mol])
  75         {
  76             mol++;
  77             if (mol >= mols->nr)
  78             {
  79                 gmx_fatal(FARGS, "Atom index out of range: %d", index[i] + 1);
  80             }
  81         }
  82         for (j = mols->index[mol]; j < mols->index[mol + 1]; j++)
  83         {
  84             if (i >= nat || index[i] != j)
  85             {
  86                 gmx_fatal(FARGS, "The index group does not consist of whole molecules");
  87             }
  88             i++;
  89         }
  90         index[nmol++] = mol;
  91     }
  92
  93     fprintf(stderr, "\nSplit group of %d atoms into %d molecules\n", nat, nmol);
  94
  95     *n = nmol;
  96 }
  97
  98 static gmx_bool precalc(t_topology top, real mass2[], real qmol[])
  99 {
 100
 101     real     mtot;
 102     real     qtot;
 103     real     qall;
 104     int      i, j, k, l;
 105     int      ai, ci;
 106     gmx_bool bNEU;
 107     ai   = 0;
 108     ci   = 0;
 109     qall = 0.0;
 110
 111
 112     for (i = 0; i < top.mols.nr; i++)
 113     {
 114         k    = top.mols.index[i];
 115         l    = top.mols.index[i + 1];
 116         mtot = 0.0;
 117         qtot = 0.0;
 118
 119         for (j = k; j < l; j++)
 120         {
 121             mtot += top.atoms.atom[j].m;
 122             qtot += top.atoms.atom[j].q;
 123         }
 124
 125         for (j = k; j < l; j++)
 126         {
 127             top.atoms.atom[j].q -= top.atoms.atom[j].m * qtot / mtot;
 128             mass2[j] = top.atoms.atom[j].m / mtot;
 129             qmol[j]  = qtot;
 130         }
 131
 132
 133         qall += qtot;
 134
 135         if (qtot < 0.0)
 136         {
 137             ai++;
 138         }
 139         if (qtot > 0.0)
 140         {
 141             ci++;
 142         }
 143     }
 144
 145     if (std::abs(qall) > 0.01)
 146     {
 147         printf("\n\nSystem not neutral (q=%f) will not calculate translational part of the dipole "
 148                "moment.\n",
 149                qall);
 150         bNEU = FALSE;
 151     }
 152     else
 153     {
 154         bNEU = TRUE;
 155     }
 156
 157     return bNEU;
 158 }
 159
 160 static void remove_jump(matrix box, int natoms, rvec xp[], rvec x[])
 161 {
 162
 163     rvec hbox;
 164     int  d, i, m;
 165
 166     for (d = 0; d < DIM; d++)
 167     {
 168         hbox[d] = 0.5 * box[d][d];
 169     }
 170     for (i = 0; i < natoms; i++)
 171     {
 172         for (m = DIM - 1; m >= 0; m--)
 173         {
 174             while (x[i][m] - xp[i][m] <= -hbox[m])
 175             {
 176                 for (d = 0; d <= m; d++)
 177                 {
 178                     x[i][d] += box[m][d];
 179                 }
 180             }
 181             while (x[i][m] - xp[i][m] > hbox[m])
 182             {
 183                 for (d = 0; d <= m; d++)
 184                 {
 185                     x[i][d] -= box[m][d];
 186                 }
 187             }
 188         }
 189     }
 190 }
 191
 192 static void calc_mj(t_topology top,
 193                     PbcType    pbcType,
 194                     matrix     box,
 195                     gmx_bool   bNoJump,
 196                     int        isize,
 197                     const int  index0[],
 198                     rvec       fr[],
 199                     rvec       mj,
 200                     real       mass2[],
 201                     real       qmol[])
 202 {
 203
 204     int   i, j, k, l;
 205     rvec  tmp;
 206     rvec  center;
 207     rvec  mt1, mt2;
 208     t_pbc pbc;
 209
 210
 211     calc_box_center(ecenterRECT, box, center);
 212
 213     if (!bNoJump)
 214     {
 215         set_pbc(&pbc, pbcType, box);
 216     }
 217
 218     clear_rvec(tmp);
 219
 220
 221     for (i = 0; i < isize; i++)
 222     {
 223         clear_rvec(mt1);
 224         clear_rvec(mt2);
 225         k = top.mols.index[index0[i]];
 226         l = top.mols.index[index0[i + 1]];
 227         for (j = k; j < l; j++)
 228         {
 229             svmul(mass2[j], fr[j], tmp);
 230             rvec_inc(mt1, tmp);
 231         }
 232
 233         if (bNoJump)
 234         {
 235             svmul(qmol[k], mt1, mt1);
 236         }
 237         else
 238         {
 239             pbc_dx(&pbc, mt1, center, mt2);
 240             svmul(qmol[k], mt2, mt1);
 241         }
 242
 243         rvec_inc(mj, mt1);
 244     }
 245 }
 246
 247
 248 static real calceps(real prefactor, real md2, real mj2, real cor, real eps_rf, gmx_bool bCOR)
 249 {
 250
 251     /* bCOR determines if the correlation is computed via
 252      * static properties (FALSE) or the correlation integral (TRUE)
 253      */
 254
 255     real epsilon = 0.0;
 256
 257
 258     if (bCOR)
 259     {
 260         epsilon = md2 - 2.0 * cor + mj2;
 261     }
 262     else
 263     {
 264         epsilon = md2 + mj2 + 2.0 * cor;
 265     }
 266
 267     if (eps_rf == 0.0)
 268     {
 269         epsilon = 1.0 + prefactor * epsilon;
 270     }
 271     else
 272     {
 273         epsilon = 2.0 * eps_rf + 1.0 + 2.0 * eps_rf * prefactor * epsilon;
 274         epsilon /= 2.0 * eps_rf + 1.0 - prefactor * epsilon;
 275     }
 276
 277
 278     return epsilon;
 279 }
 280
 281
 282 static real calc_cacf(FILE* fcacf, real prefactor, real cacf[], real time[], int nfr, const int vfr[], int ei, int nshift)
 283 {
 284
 285     int  i;
 286     real deltat = 0.0;
 287     real rfr;
 288     real sigma     = 0.0;
 289     real sigma_ret = 0.0;
 290
 291     if (nfr > 1)
 292     {
 293         i = 0;
 294
 295         while (i < nfr)
 296         {
 297             real corint;
 298
 299             rfr = static_cast<real>(nfr + i) / nshift;
 300             cacf[i] /= rfr;
 301
 302             if (time[vfr[i]] <= time[vfr[ei]])
 303             {
 304                 sigma_ret = sigma;
 305             }
 306
 307             fprintf(fcacf, "%.3f\t%10.6g\t%10.6g\n", time[vfr[i]], cacf[i], sigma);
 308
 309             if ((i + 1) < nfr)
 310             {
 311                 deltat = (time[vfr[i + 1]] - time[vfr[i]]);
 312             }
 313             corint = 2.0 * deltat * cacf[i] * prefactor;
 314             if (i == 0 || (i + 1) == nfr)
 315             {
 316                 corint *= 0.5;
 317             }
 318             sigma += corint;
 319
 320             i++;
 321         }
 322     }
 323     else
 324     {
 325         printf("Too less points.\n");
 326     }
 327
 328     return sigma_ret;
 329 }
 330
 331 static void calc_mjdsp(FILE* fmjdsp, real prefactor, real dsp2[], real time[], int nfr, const real refr[])
 332 {
 333
 334     int i;
 335
 336     fprintf(fmjdsp, "#Prefactor fit E-H: 1 / 6.0*V*k_B*T: %g\n", prefactor);
 337
 338     for (i = 0; i < nfr; i++)
 339     {
 340
 341         if (refr[i] != 0.0)
 342         {
 343             dsp2[i] *= prefactor / refr[i];
 344             fprintf(fmjdsp, "%.3f\t%10.6g\n", time[i], dsp2[i]);
 345         }
 346     }
 347 }
 348
 349
 350 static void dielectric(FILE*                   fmj,
 351                        FILE*                   fmd,
 352                        FILE*                   outf,
 353                        FILE*                   fcur,
 354                        FILE*                   mcor,
 355                        FILE*                   fmjdsp,
 356                        gmx_bool                bNoJump,
 357                        gmx_bool                bACF,
 358                        gmx_bool                bINT,
 359                        PbcType                 pbcType,
 360                        t_topology              top,
 361                        t_trxframe              fr,
 362                        real                    temp,
 363                        real                    bfit,
 364                        real                    efit,
 365                        real                    bvit,
 366                        real                    evit,
 367                        t_trxstatus*            status,
 368                        int                     isize,
 369                        int                     nmols,
 370                        int                     nshift,
 371                        const int*              index0,
 372                        int                     indexm[],
 373                        real                    mass2[],
 374                        real                    qmol[],
 375                        real                    eps_rf,
 376                        const gmx_output_env_t* oenv)
 377 {
 378     int   i, j;
 379     int   valloc, nalloc, nfr, nvfr;
 380     int   vshfr;
 381     real* xshfr       = nullptr;
 382     int*  vfr         = nullptr;
 383     real  refr        = 0.0;
 384     real* cacf        = nullptr;
 385     real* time        = nullptr;
 386     real* djc         = nullptr;
 387     real  corint      = 0.0;
 388     real  prefactorav = 0.0;
 389     real  prefactor   = 0.0;
 390     real  volume;
 391     real  volume_av = 0.0;
 392     real  dk_s, dk_d, dk_f;
 393     real  mj    = 0.0;
 394     real  mj2   = 0.0;
 395     real  mjd   = 0.0;
 396     real  mjdav = 0.0;
 397     real  md2   = 0.0;
 398     real  mdav2 = 0.0;
 399     real  sgk;
 400     rvec  mja_tmp;
 401     rvec  mjd_tmp;
 402     rvec  mdvec;
 403     rvec* mu    = nullptr;
 404     rvec* xp    = nullptr;
 405     rvec* v0    = nullptr;
 406     rvec* mjdsp = nullptr;
 407     real* dsp2  = nullptr;
 408     real  t0;
 409     real  rtmp;
 410
 411     rvec  tmp;
 412     rvec* mtrans = nullptr;
 413
 414     /*
 415      * Variables for the least-squares fit for Einstein-Helfand and Green-Kubo
 416      */
 417
 418     int         bi, ei, ie, ii;
 419     real        rest  = 0.0;
 420     real        sigma = 0.0;
 421     real        malt  = 0.0;
 422     real        err   = 0.0;
 423     real*       xfit;
 424     real*       yfit;
 425     gmx_rmpbc_t gpbc = nullptr;
 426
 427     /*
 428      * indices for EH
 429      */
 430
 431     ei = 0;
 432     bi = 0;
 433
 434     /*
 435      * indices for GK
 436      */
 437
 438     ii  = 0;
 439     ie  = 0;
 440     t0  = 0;
 441     sgk = 0.0;
 442
 443
 444     /* This is the main loop over frames */
 445
 446
 447     nfr = 0;
 448
 449
 450     nvfr   = 0;
 451     vshfr  = 0;
 452     nalloc = 0;
 453     valloc = 0;
 454
 455     clear_rvec(mja_tmp);
 456     clear_rvec(mjd_tmp);
 457     clear_rvec(mdvec);
 458     clear_rvec(tmp);
 459     gpbc = gmx_rmpbc_init(&top.idef, pbcType, fr.natoms);
 460
 461     do
 462     {
 463
 464         refr = (nfr + 1);
 465
 466         if (nfr >= nalloc)
 467         {
 468             nalloc += 100;
 469             srenew(time, nalloc);
 470             srenew(mu, nalloc);
 471             srenew(mjdsp, nalloc);
 472             srenew(dsp2, nalloc);
 473             srenew(mtrans, nalloc);
 474             srenew(xshfr, nalloc);
 475
 476             for (i = nfr; i < nalloc; i++)
 477             {
 478                 clear_rvec(mjdsp[i]);
 479                 clear_rvec(mu[i]);
 480                 clear_rvec(mtrans[i]);
 481                 dsp2[i]  = 0.0;
 482                 xshfr[i] = 0.0;
 483             }
 484         }
 485         GMX_RELEASE_ASSERT(time != nullptr, "Memory not allocated correctly - time array is NULL");
 486
 487         if (nfr == 0)
 488         {
 489             t0 = fr.time;
 490         }
 491
 492         time[nfr] = fr.time - t0;
 493
 494         if (time[nfr] <= bfit)
 495         {
 496             bi = nfr;
 497         }
 498         if (time[nfr] <= efit)
 499         {
 500             ei = nfr;
 501         }
 502
 503         if (bNoJump)
 504         {
 505
 506             if (xp)
 507             {
 508                 remove_jump(fr.box, fr.natoms, xp, fr.x);
 509             }
 510             else
 511             {
 512                 snew(xp, fr.natoms);
 513             }
 514
 515             for (i = 0; i < fr.natoms; i++)
 516             {
 517                 copy_rvec(fr.x[i], xp[i]);
 518             }
 519         }
 520
 521         gmx_rmpbc_trxfr(gpbc, &fr);
 522
 523         calc_mj(top, pbcType, fr.box, bNoJump, nmols, indexm, fr.x, mtrans[nfr], mass2, qmol);
 524
 525         for (i = 0; i < isize; i++)
 526         {
 527             j = index0[i];
 528             svmul(top.atoms.atom[j].q, fr.x[j], fr.x[j]);
 529             rvec_inc(mu[nfr], fr.x[j]);
 530         }
 531
 532         /*if(mod(nfr,nshift)==0){*/
 533         if (nfr % nshift == 0)
 534         {
 535             for (j = nfr; j >= 0; j--)
 536             {
 537                 rvec_sub(mtrans[nfr], mtrans[j], tmp);
 538                 dsp2[nfr - j] += norm2(tmp);
 539                 xshfr[nfr - j] += 1.0;
 540             }
 541         }
 542
 543         if (fr.bV)
 544         {
 545             if (nvfr >= valloc)
 546             {
 547                 valloc += 100;
 548                 srenew(vfr, valloc);
 549                 if (bINT)
 550                 {
 551                     srenew(djc, valloc);
 552                 }
 553                 srenew(v0, valloc);
 554                 if (bACF)
 555                 {
 556                     srenew(cacf, valloc);
 557                 }
 558             }
 559             if (time[nfr] <= bvit)
 560             {
 561                 ii = nvfr;
 562             }
 563             if (time[nfr] <= evit)
 564             {
 565                 ie = nvfr;
 566             }
 567             vfr[nvfr] = nfr;
 568             clear_rvec(v0[nvfr]);
 569             if (bACF)
 570             {
 571                 cacf[nvfr] = 0.0;
 572             }
 573             if (bINT)
 574             {
 575                 djc[nvfr] = 0.0;
 576             }
 577             for (i = 0; i < isize; i++)
 578             {
 579                 j = index0[i];
 580                 svmul(mass2[j], fr.v[j], fr.v[j]);
 581                 svmul(qmol[j], fr.v[j], fr.v[j]);
 582                 rvec_inc(v0[nvfr], fr.v[j]);
 583             }
 584
 585             fprintf(fcur, "%.3f\t%.6f\t%.6f\t%.6f\n", time[nfr], v0[nfr][XX], v0[nfr][YY], v0[nfr][ZZ]);
 586             if (bACF || bINT)
 587             {
 588                 /*if(mod(nvfr,nshift)==0){*/
 589                 if (nvfr % nshift == 0)
 590                 {
 591                     for (j = nvfr; j >= 0; j--)
 592                     {
 593                         if (bACF)
 594                         {
 595                             cacf[nvfr - j] += iprod(v0[nvfr], v0[j]);
 596                         }
 597                         if (bINT)
 598                         {
 599                             djc[nvfr - j] += iprod(mu[vfr[j]], v0[nvfr]);
 600                         }
 601                     }
 602                     vshfr++;
 603                 }
 604             }
 605             nvfr++;
 606         }
 607
 608         volume = det(fr.box);
 609         volume_av += volume;
 610
 611         rvec_inc(mja_tmp, mtrans[nfr]);
 612         mjd += iprod(mu[nfr], mtrans[nfr]);
 613         rvec_inc(mdvec, mu[nfr]);
 614
 615         mj2 += iprod(mtrans[nfr], mtrans[nfr]);
 616         md2 += iprod(mu[nfr], mu[nfr]);
 617
 618         fprintf(fmj,
 619                 "%.3f\t%8.5f\t%8.5f\t%8.5f\t%8.5f\t%8.5f\n",
 620                 time[nfr],
 621                 mtrans[nfr][XX],
 622                 mtrans[nfr][YY],
 623                 mtrans[nfr][ZZ],
 624                 mj2 / refr,
 625                 norm(mja_tmp) / refr);
 626         fprintf(fmd,
 627                 "%.3f\t%8.5f\t%8.5f\t%8.5f\t%8.5f\t%8.5f\n",
 628                 time[nfr],
 629                 mu[nfr][XX],
 630                 mu[nfr][YY],
 631                 mu[nfr][ZZ],
 632                 md2 / refr,
 633                 norm(mdvec) / refr);
 634
 635         nfr++;
 636
 637     } while (read_next_frame(oenv, status, &fr));
 638
 639     gmx_rmpbc_done(gpbc);
 640
 641     volume_av /= refr;
 642
 643     prefactor = 1.0;
 644     prefactor /= 3.0 * gmx::c_epsilon0 * volume_av * gmx::c_boltz * temp;
 645
 646
 647     prefactorav = gmx::c_electronCharge * gmx::c_electronCharge;
 648     prefactorav /= volume_av * gmx::c_boltzmann * temp * gmx::c_nano * 6.0;
 649
 650     fprintf(stderr, "Prefactor fit E-H: 1 / 6.0*V*k_B*T: %g\n", prefactorav);
 651
 652     calc_mjdsp(fmjdsp, prefactorav, dsp2, time, nfr, xshfr);
 653
 654     /*
 655      * Now we can average and calculate the correlation functions
 656      */
 657
 658
 659     mj2 /= refr;
 660     mjd /= refr;
 661     md2 /= refr;
 662
 663     svmul(1.0 / refr, mdvec, mdvec);
 664     svmul(1.0 / refr, mja_tmp, mja_tmp);
 665
 666     mdav2 = norm2(mdvec);
 667     mj    = norm2(mja_tmp);
 668     mjdav = iprod(mdvec, mja_tmp);
 669
 670
 671     printf("\n\nAverage translational dipole moment M_J [enm] after %d frames (|M|^2): %f %f %f "
 672            "(%f)\n",
 673            nfr,
 674            mja_tmp[XX],
 675            mja_tmp[YY],
 676            mja_tmp[ZZ],
 677            mj2);
 678     printf("\n\nAverage molecular dipole moment M_D [enm] after %d frames (|M|^2): %f %f %f (%f)\n",
 679            nfr,
 680            mdvec[XX],
 681            mdvec[YY],
 682            mdvec[ZZ],
 683            md2);
 684
 685     if (v0 != nullptr)
 686     {
 687         if (bINT)
 688         {
 689
 690             printf("\nCalculating M_D - current correlation integral ... \n");
 691             corint = calc_cacf(mcor, prefactorav / EPSI0, djc, time, nvfr, vfr, ie, nshift);
 692         }
 693
 694         if (bACF)
 695         {
 696
 697             printf("\nCalculating current autocorrelation ... \n");
 698             sgk = calc_cacf(outf, prefactorav / gmx::c_pico, cacf, time, nvfr, vfr, ie, nshift);
 699
 700             if (ie > ii)
 701             {
 702
 703                 snew(xfit, ie - ii + 1);
 704                 snew(yfit, ie - ii + 1);
 705
 706                 for (i = ii; i <= ie; i++)
 707                 {
 708
 709                     xfit[i - ii] = std::log(time[vfr[i]]);
 710                     rtmp         = std::abs(cacf[i]);
 711                     yfit[i - ii] = std::log(rtmp);
 712                 }
 713
 714                 lsq_y_ax_b(ie - ii, xfit, yfit, &sigma, &malt, &err, &rest);
 715
 716                 malt = std::exp(malt);
 717
 718                 sigma += 1.0;
 719
 720                 malt *= prefactorav * 2.0e12 / sigma;
 721
 722                 sfree(xfit);
 723                 sfree(yfit);
 724             }
 725         }
 726     }
 727
 728
 729     /* Calculation of the dielectric constant */
 730
 731     fprintf(stderr, "\n********************************************\n");
 732     dk_s = calceps(prefactor, md2, mj2, mjd, eps_rf, FALSE);
 733     fprintf(stderr, "\nAbsolute values:\n epsilon=%f\n", dk_s);
 734     fprintf(stderr, " <M_D^2> , <M_J^2>, <(M_J*M_D)^2>:  (%f, %f, %f)\n\n", md2, mj2, mjd);
 735     fprintf(stderr, "********************************************\n");
 736
 737
 738     dk_f = calceps(prefactor, md2 - mdav2, mj2 - mj, mjd - mjdav, eps_rf, FALSE);
 739     fprintf(stderr, "\n\nFluctuations:\n epsilon=%f\n\n", dk_f);
 740     fprintf(stderr, "\n deltaM_D , deltaM_J, deltaM_JD:  (%f, %f, %f)\n", md2 - mdav2, mj2 - mj, mjd - mjdav);
 741     fprintf(stderr, "\n********************************************\n");
 742     if (bINT)
 743     {
 744         dk_d = calceps(prefactor, md2 - mdav2, mj2 - mj, corint, eps_rf, TRUE);
 745         fprintf(stderr, "\nStatic dielectric constant using integral and fluctuations: %f\n", dk_d);
 746         fprintf(stderr, "\n < M_JM_D > via integral:  %.3f\n", -1.0 * corint);
 747     }
 748
 749     fprintf(stderr, "\n***************************************************");
 750     fprintf(stderr, "\n\nAverage volume V=%f nm^3 at T=%f K\n", volume_av, temp);
 751     fprintf(stderr, "and corresponding refactor 1.0 / 3.0*V*k_B*T*EPSILON_0: %f \n", prefactor);
 752
 753
 754     if (bACF && (ii < nvfr))
 755     {
 756         fprintf(stderr, "Integral and integrated fit to the current acf yields at t=%f:\n", time[vfr[ii]]);
 757         fprintf(stderr, "sigma=%8.3f (pure integral: %.3f)\n", sgk - malt * std::pow(time[vfr[ii]], sigma), sgk);
 758     }
 759
 760     if (ei > bi)
 761     {
 762         fprintf(stderr, "\nStart fit at %f ps (%f).\n", time[bi], bfit);
 763         fprintf(stderr, "End fit at %f ps (%f).\n\n", time[ei], efit);
 764
 765         snew(xfit, ei - bi + 1);
 766         snew(yfit, ei - bi + 1);
 767
 768         for (i = bi; i <= ei; i++)
 769         {
 770             xfit[i - bi] = time[i];
 771             yfit[i - bi] = dsp2[i];
 772         }
 773
 774         lsq_y_ax_b(ei - bi, xfit, yfit, &sigma, &malt, &err, &rest);
 775
 776         sigma *= 1e12;
 777         dk_d = calceps(prefactor, md2, 0.5 * malt / prefactorav, corint, eps_rf, TRUE);
 778
 779         fprintf(stderr,
 780                 "Einstein-Helfand fit to the MSD of the translational dipole moment yields:\n\n");
 781         fprintf(stderr, "sigma=%.4f\n", sigma);
 782         fprintf(stderr, "translational fraction of M^2: %.4f\n", 0.5 * malt / prefactorav);
 783         fprintf(stderr, "Dielectric constant using EH: %.4f\n", dk_d);
 784
 785         sfree(xfit);
 786         sfree(yfit);
 787     }
 788     else
 789     {
 790         fprintf(stderr, "Too few points for a fit.\n");
 791     }
 792
 793
 794     if (v0 != nullptr)
 795     {
 796         sfree(v0);
 797     }
 798     if (bACF)
 799     {
 800         sfree(cacf);
 801     }
 802     if (bINT)
 803     {
 804         sfree(djc);
 805     }
 806
 807     sfree(time);
 808
 809
 810     sfree(mjdsp);
 811     sfree(mu);
 812 }
 813
 814
 815 int gmx_current(int argc, char* argv[])
 816 {
 817
 818     static int      nshift  = 1000;
 819     static real     temp    = 300.0;
 820     static real     eps_rf  = 0.0;
 821     static gmx_bool bNoJump = TRUE;
 822     static real     bfit    = 100.0;
 823     static real     bvit    = 0.5;
 824     static real     efit    = 400.0;
 825     static real     evit    = 5.0;
 826     t_pargs         pa[]    = {
 827         { "-sh",
 828           FALSE,
 829           etINT,
 830           { &nshift },
 831           "Shift of the frames for averaging the correlation functions and the mean-square "
 832           "displacement." },
 833         { "-nojump", FALSE, etBOOL, { &bNoJump }, "Removes jumps of atoms across the box." },
 834         { "-eps",
 835           FALSE,
 836           etREAL,
 837           { &eps_rf },
 838           "Dielectric constant of the surrounding medium. The value zero corresponds to "
 839           "infinity (tin-foil boundary conditions)." },
 840         { "-bfit",
 841           FALSE,
 842           etREAL,
 843           { &bfit },
 844           "Begin of the fit of the straight line to the MSD of the translational fraction "
 845           "of the dipole moment." },
 846         { "-efit",
 847           FALSE,
 848           etREAL,
 849           { &efit },
 850           "End of the fit of the straight line to the MSD of the translational fraction of "
 851           "the dipole moment." },
 852         { "-bvit",
 853           FALSE,
 854           etREAL,
 855           { &bvit },
 856           "Begin of the fit of the current autocorrelation function to a*t^b." },
 857         { "-evit",
 858           FALSE,
 859           etREAL,
 860           { &evit },
 861           "End of the fit of the current autocorrelation function to a*t^b." },
 862         { "-temp", FALSE, etREAL, { &temp }, "Temperature for calculating epsilon." }
 863     };
 864
 865     gmx_output_env_t* oenv;
 866     t_topology        top;
 867     char**            grpname = nullptr;
 868     const char*       indexfn;
 869     t_trxframe        fr;
 870     real*             mass2 = nullptr;
 871     matrix            box;
 872     int*              index0;
 873     int*              indexm = nullptr;
 874     int               isize;
 875     t_trxstatus*      status;
 876     int               flags = 0;
 877     gmx_bool          bACF;
 878     gmx_bool          bINT;
 879     PbcType           pbcType = PbcType::Unset;
 880     int               nmols;
 881     int               i;
 882     real*             qmol;
 883     FILE*             outf   = nullptr;
 884     FILE*             mcor   = nullptr;
 885     FILE*             fmj    = nullptr;
 886     FILE*             fmd    = nullptr;
 887     FILE*             fmjdsp = nullptr;
 888     FILE*             fcur   = nullptr;
 889     t_filenm          fnm[]  = { { efTPS, nullptr, nullptr, ffREAD }, /* this is for the topology */
 890                        { efNDX, nullptr, nullptr, ffOPTRD },
 891                        { efTRX, "-f", nullptr, ffREAD }, /* and this for the trajectory */
 892                        { efXVG, "-o", "current", ffWRITE },
 893                        { efXVG, "-caf", "caf", ffOPTWR },
 894                        { efXVG, "-dsp", "dsp", ffWRITE },
 895                        { efXVG, "-md", "md", ffWRITE },
 896                        { efXVG, "-mj", "mj", ffWRITE },
 897                        { efXVG, "-mc", "mc", ffOPTWR } };
 898
 899 #define NFILE asize(fnm)
 900
 901
 902     const char* desc[] = {
 903         "[THISMODULE] is a tool for calculating the current autocorrelation function, the "
 904         "correlation",
 905         "of the rotational and translational dipole moment of the system, and the resulting static",
 906         "dielectric constant. To obtain a reasonable result, the index group has to be neutral.",
 907         "Furthermore, the routine is capable of extracting the static conductivity from the "
 908         "current ",
 909         "autocorrelation function, if velocities are given. Additionally, an Einstein-Helfand fit ",
 910         "can be used to obtain the static conductivity.",
 911         "[PAR]",
 912         "The flag [TT]-caf[tt] is for the output of the current autocorrelation function and "
 913         "[TT]-mc[tt] writes the",
 914         "correlation of the rotational and translational part of the dipole moment in the "
 915         "corresponding",
 916         "file. However, this option is only available for trajectories containing velocities.",
 917         "Options [TT]-sh[tt] and [TT]-tr[tt] are responsible for the averaging and integration of "
 918         "the",
 919         "autocorrelation functions. Since averaging proceeds by shifting the starting point",
 920         "through the trajectory, the shift can be modified with [TT]-sh[tt] to enable the choice "
 921         "of uncorrelated",
 922         "starting points. Towards the end, statistical inaccuracy grows and integrating the",
 923         "correlation function only yields reliable values until a certain point, depending on",
 924         "the number of frames. The option [TT]-tr[tt] controls the region of the integral taken "
 925         "into account",
 926         "for calculating the static dielectric constant.",
 927         "[PAR]",
 928         "Option [TT]-temp[tt] sets the temperature required for the computation of the static "
 929         "dielectric constant.",
 930         "[PAR]",
 931         "Option [TT]-eps[tt] controls the dielectric constant of the surrounding medium for "
 932         "simulations using",
 933         "a Reaction Field or dipole corrections of the Ewald summation ([TT]-eps[tt]\\=0 "
 934         "corresponds to",
 935         "tin-foil boundary conditions).",
 936         "[PAR]",
 937         "[TT]-[no]nojump[tt] unfolds the coordinates to allow free diffusion. This is required to "
 938         "get a continuous",
 939         "translational dipole moment, required for the Einstein-Helfand fit. The results from the "
 940         "fit allow",
 941         "the determination of the dielectric constant for system of charged molecules. However, it "
 942         "is also possible to extract",
 943         "the dielectric constant from the fluctuations of the total dipole moment in folded "
 944         "coordinates. But this",
 945         "option has to be used with care, since only very short time spans fulfill the "
 946         "approximation that the density",
 947         "of the molecules is approximately constant and the averages are already converged. To be "
 948         "on the safe side,",
 949         "the dielectric constant should be calculated with the help of the Einstein-Helfand method "
 950         "for",
 951         "the translational part of the dielectric constant."
 952     };
 953
 954
 955     /* At first the arguments will be parsed and the system information processed */
 956     if (!parse_common_args(
 957                 &argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW, NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
 958     {
 959         return 0;
 960     }
 961
 962     bACF = opt2bSet("-caf", NFILE, fnm);
 963     bINT = opt2bSet("-mc", NFILE, fnm);
 964
 965     read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &pbcType, nullptr, nullptr, box, TRUE);
 966
 967     indexfn = ftp2fn_null(efNDX, NFILE, fnm);
 968     snew(grpname, 1);
 969
 970     get_index(&(top.atoms), indexfn, 1, &isize, &index0, grpname);
 971
 972     flags = flags | TRX_READ_X | TRX_READ_V;
 973
 974     read_first_frame(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &fr, flags);
 975
 976     snew(mass2, top.atoms.nr);
 977     snew(qmol, top.atoms.nr);
 978
 979     precalc(top, mass2, qmol);
 980
 981
 982     snew(indexm, isize);
 983
 984     for (i = 0; i < isize; i++)
 985     {
 986         indexm[i] = index0[i];
 987     }
 988
 989     nmols = isize;
 990
 991
 992     index_atom2mol(&nmols, indexm, &top.mols);
 993
 994     if (fr.bV)
 995     {
 996         if (bACF)
 997         {
 998             outf = xvgropen(opt2fn("-caf", NFILE, fnm),
 999                             "Current autocorrelation function",
1000                             output_env_get_xvgr_tlabel(oenv),
1001                             "ACF (e nm/ps)\\S2",
1002                             oenv);
1003             fprintf(outf, "# time\t acf\t average \t std.dev\n");
1004         }
1005         fcur = xvgropen(opt2fn("-o", NFILE, fnm),
1006                         "Current",
1007                         output_env_get_xvgr_tlabel(oenv),
1008                         "J(t) (e nm/ps)",
1009                         oenv);
1010         fprintf(fcur, "# time\t Jx\t Jy \t J_z \n");
1011         if (bINT)
1012         {
1013             mcor = xvgropen(opt2fn("-mc", NFILE, fnm),
1014                             "M\\sD\\N - current  autocorrelation function",
1015                             output_env_get_xvgr_tlabel(oenv),
1016                             "< M\\sD\\N (0)\\c7\\CJ(t) >  (e nm/ps)\\S2",
1017                             oenv);
1018             fprintf(mcor, "# time\t M_D(0) J(t) acf \t Integral acf\n");
1019         }
1020     }
1021
1022     fmj = xvgropen(opt2fn("-mj", NFILE, fnm),
1023                    "Averaged translational part of M",
1024                    output_env_get_xvgr_tlabel(oenv),
1025                    "< M\\sJ\\N > (enm)",
1026                    oenv);
1027     fprintf(fmj, "# time\t x\t y \t z \t average of M_J^2 \t std.dev\n");
1028     fmd = xvgropen(opt2fn("-md", NFILE, fnm),
1029                    "Averaged rotational part of M",
1030                    output_env_get_xvgr_tlabel(oenv),
1031                    "< M\\sD\\N > (enm)",
1032                    oenv);
1033     fprintf(fmd, "# time\t x\t y \t z \t average of M_D^2 \t std.dev\n");
1034
1035     fmjdsp = xvgropen(
1036             opt2fn("-dsp", NFILE, fnm),
1037             "MSD of the squared translational dipole moment M",
1038             output_env_get_xvgr_tlabel(oenv),
1039             "<|M\\sJ\\N(t)-M\\sJ\\N(0)|\\S2\\N > / 6.0*V*k\\sB\\N*T / Sm\\S-1\\Nps\\S-1\\N",
1040             oenv);
1041
1042
1043     /* System information is read and prepared, dielectric() processes the frames
1044      * and calculates the requested quantities */
1045
1046     dielectric(fmj,
1047                fmd,
1048                outf,
1049                fcur,
1050                mcor,
1051                fmjdsp,
1052                bNoJump,
1053                bACF,
1054                bINT,
1055                pbcType,
1056                top,
1057                fr,
1058                temp,
1059                bfit,
1060                efit,
1061                bvit,
1062                evit,
1063                status,
1064                isize,
1065                nmols,
1066                nshift,
1067                index0,
1068                indexm,
1069                mass2,
1070                qmol,
1071                eps_rf,
1072                oenv);
1073
1074     xvgrclose(fmj);
1075     xvgrclose(fmd);
1076     xvgrclose(fmjdsp);
1077     if (fr.bV)
1078     {
1079         if (bACF)
1080         {
1081             xvgrclose(outf);
1082         }
1083         xvgrclose(fcur);
1084         if (bINT)
1085         {
1086             xvgrclose(mcor);
1087         }
1088     }
1089
1090     return 0;
1091 }