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