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