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                     int        ePBC,
 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, ePBC, 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                        int                     ePBC,
 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, ePBC, 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, ePBC, 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, "%.3f\t%8.5f\t%8.5f\t%8.5f\t%8.5f\t%8.5f\n", time[nfr], mtrans[nfr][XX],
 618                 mtrans[nfr][YY], mtrans[nfr][ZZ], mj2 / refr, norm(mja_tmp) / refr);
 619         fprintf(fmd, "%.3f\t%8.5f\t%8.5f\t%8.5f\t%8.5f\t%8.5f\n", time[nfr], mu[nfr][XX],
 620                 mu[nfr][YY], mu[nfr][ZZ], md2 / refr, norm(mdvec) / refr);
 621
 622         nfr++;
 623
 624     } while (read_next_frame(oenv, status, &fr));
 625
 626     gmx_rmpbc_done(gpbc);
 627
 628     volume_av /= refr;
 629
 630     prefactor = 1.0;
 631     prefactor /= 3.0 * EPSILON0 * volume_av * BOLTZ * temp;
 632
 633
 634     prefactorav = E_CHARGE * E_CHARGE;
 635     prefactorav /= volume_av * BOLTZMANN * temp * NANO * 6.0;
 636
 637     fprintf(stderr, "Prefactor fit E-H: 1 / 6.0*V*k_B*T: %g\n", prefactorav);
 638
 639     calc_mjdsp(fmjdsp, prefactorav, dsp2, time, nfr, xshfr);
 640
 641     /*
 642      * Now we can average and calculate the correlation functions
 643      */
 644
 645
 646     mj2 /= refr;
 647     mjd /= refr;
 648     md2 /= refr;
 649
 650     svmul(1.0 / refr, mdvec, mdvec);
 651     svmul(1.0 / refr, mja_tmp, mja_tmp);
 652
 653     mdav2 = norm2(mdvec);
 654     mj    = norm2(mja_tmp);
 655     mjdav = iprod(mdvec, mja_tmp);
 656
 657
 658     printf("\n\nAverage translational dipole moment M_J [enm] after %d frames (|M|^2): %f %f %f "
 659            "(%f)\n",
 660            nfr, mja_tmp[XX], mja_tmp[YY], mja_tmp[ZZ], mj2);
 661     printf("\n\nAverage molecular dipole moment M_D [enm] after %d frames (|M|^2): %f %f %f (%f)\n",
 662            nfr, mdvec[XX], mdvec[YY], mdvec[ZZ], md2);
 663
 664     if (v0 != nullptr)
 665     {
 666         if (bINT)
 667         {
 668
 669             printf("\nCalculating M_D - current correlation integral ... \n");
 670             corint = calc_cacf(mcor, prefactorav / EPSI0, djc, time, nvfr, vfr, ie, nshift);
 671         }
 672
 673         if (bACF)
 674         {
 675
 676             printf("\nCalculating current autocorrelation ... \n");
 677             sgk = calc_cacf(outf, prefactorav / PICO, cacf, time, nvfr, vfr, ie, nshift);
 678
 679             if (ie > ii)
 680             {
 681
 682                 snew(xfit, ie - ii + 1);
 683                 snew(yfit, ie - ii + 1);
 684
 685                 for (i = ii; i <= ie; i++)
 686                 {
 687
 688                     xfit[i - ii] = std::log(time[vfr[i]]);
 689                     rtmp         = std::abs(cacf[i]);
 690                     yfit[i - ii] = std::log(rtmp);
 691                 }
 692
 693                 lsq_y_ax_b(ie - ii, xfit, yfit, &sigma, &malt, &err, &rest);
 694
 695                 malt = std::exp(malt);
 696
 697                 sigma += 1.0;
 698
 699                 malt *= prefactorav * 2.0e12 / sigma;
 700
 701                 sfree(xfit);
 702                 sfree(yfit);
 703             }
 704         }
 705     }
 706
 707
 708     /* Calculation of the dielectric constant */
 709
 710     fprintf(stderr, "\n********************************************\n");
 711     dk_s = calceps(prefactor, md2, mj2, mjd, eps_rf, FALSE);
 712     fprintf(stderr, "\nAbsolute values:\n epsilon=%f\n", dk_s);
 713     fprintf(stderr, " <M_D^2> , <M_J^2>, <(M_J*M_D)^2>:  (%f, %f, %f)\n\n", md2, mj2, mjd);
 714     fprintf(stderr, "********************************************\n");
 715
 716
 717     dk_f = calceps(prefactor, md2 - mdav2, mj2 - mj, mjd - mjdav, eps_rf, FALSE);
 718     fprintf(stderr, "\n\nFluctuations:\n epsilon=%f\n\n", dk_f);
 719     fprintf(stderr, "\n deltaM_D , deltaM_J, deltaM_JD:  (%f, %f, %f)\n", md2 - mdav2, mj2 - mj,
 720             mjd - mjdav);
 721     fprintf(stderr, "\n********************************************\n");
 722     if (bINT)
 723     {
 724         dk_d = calceps(prefactor, md2 - mdav2, mj2 - mj, corint, eps_rf, TRUE);
 725         fprintf(stderr, "\nStatic dielectric constant using integral and fluctuations: %f\n", dk_d);
 726         fprintf(stderr, "\n < M_JM_D > via integral:  %.3f\n", -1.0 * corint);
 727     }
 728
 729     fprintf(stderr, "\n***************************************************");
 730     fprintf(stderr, "\n\nAverage volume V=%f nm^3 at T=%f K\n", volume_av, temp);
 731     fprintf(stderr, "and corresponding refactor 1.0 / 3.0*V*k_B*T*EPSILON_0: %f \n", prefactor);
 732
 733
 734     if (bACF && (ii < nvfr))
 735     {
 736         fprintf(stderr, "Integral and integrated fit to the current acf yields at t=%f:\n", time[vfr[ii]]);
 737         fprintf(stderr, "sigma=%8.3f (pure integral: %.3f)\n",
 738                 sgk - malt * std::pow(time[vfr[ii]], sigma), sgk);
 739     }
 740
 741     if (ei > bi)
 742     {
 743         fprintf(stderr, "\nStart fit at %f ps (%f).\n", time[bi], bfit);
 744         fprintf(stderr, "End fit at %f ps (%f).\n\n", time[ei], efit);
 745
 746         snew(xfit, ei - bi + 1);
 747         snew(yfit, ei - bi + 1);
 748
 749         for (i = bi; i <= ei; i++)
 750         {
 751             xfit[i - bi] = time[i];
 752             yfit[i - bi] = dsp2[i];
 753         }
 754
 755         lsq_y_ax_b(ei - bi, xfit, yfit, &sigma, &malt, &err, &rest);
 756
 757         sigma *= 1e12;
 758         dk_d = calceps(prefactor, md2, 0.5 * malt / prefactorav, corint, eps_rf, TRUE);
 759
 760         fprintf(stderr,
 761                 "Einstein-Helfand fit to the MSD of the translational dipole moment yields:\n\n");
 762         fprintf(stderr, "sigma=%.4f\n", sigma);
 763         fprintf(stderr, "translational fraction of M^2: %.4f\n", 0.5 * malt / prefactorav);
 764         fprintf(stderr, "Dielectric constant using EH: %.4f\n", dk_d);
 765
 766         sfree(xfit);
 767         sfree(yfit);
 768     }
 769     else
 770     {
 771         fprintf(stderr, "Too few points for a fit.\n");
 772     }
 773
 774
 775     if (v0 != nullptr)
 776     {
 777         sfree(v0);
 778     }
 779     if (bACF)
 780     {
 781         sfree(cacf);
 782     }
 783     if (bINT)
 784     {
 785         sfree(djc);
 786     }
 787
 788     sfree(time);
 789
 790
 791     sfree(mjdsp);
 792     sfree(mu);
 793 }
 794
 795
 796 int gmx_current(int argc, char* argv[])
 797 {
 798
 799     static int      nshift  = 1000;
 800     static real     temp    = 300.0;
 801     static real     eps_rf  = 0.0;
 802     static gmx_bool bNoJump = TRUE;
 803     static real     bfit    = 100.0;
 804     static real     bvit    = 0.5;
 805     static real     efit    = 400.0;
 806     static real     evit    = 5.0;
 807     t_pargs         pa[]    = {
 808         { "-sh",
 809           FALSE,
 810           etINT,
 811           { &nshift },
 812           "Shift of the frames for averaging the correlation functions and the mean-square "
 813           "displacement." },
 814         { "-nojump", FALSE, etBOOL, { &bNoJump }, "Removes jumps of atoms across the box." },
 815         { "-eps",
 816           FALSE,
 817           etREAL,
 818           { &eps_rf },
 819           "Dielectric constant of the surrounding medium. The value zero corresponds to "
 820           "infinity (tin-foil boundary conditions)." },
 821         { "-bfit",
 822           FALSE,
 823           etREAL,
 824           { &bfit },
 825           "Begin of the fit of the straight line to the MSD of the translational fraction "
 826           "of the dipole moment." },
 827         { "-efit",
 828           FALSE,
 829           etREAL,
 830           { &efit },
 831           "End of the fit of the straight line to the MSD of the translational fraction of "
 832           "the dipole moment." },
 833         { "-bvit",
 834           FALSE,
 835           etREAL,
 836           { &bvit },
 837           "Begin of the fit of the current autocorrelation function to a*t^b." },
 838         { "-evit",
 839           FALSE,
 840           etREAL,
 841           { &evit },
 842           "End of the fit of the current autocorrelation function to a*t^b." },
 843         { "-temp", FALSE, etREAL, { &temp }, "Temperature for calculating epsilon." }
 844     };
 845
 846     gmx_output_env_t* oenv;
 847     t_topology        top;
 848     char**            grpname = nullptr;
 849     const char*       indexfn;
 850     t_trxframe        fr;
 851     real*             mass2 = nullptr;
 852     matrix            box;
 853     int*              index0;
 854     int*              indexm = nullptr;
 855     int               isize;
 856     t_trxstatus*      status;
 857     int               flags = 0;
 858     gmx_bool          bACF;
 859     gmx_bool          bINT;
 860     int               ePBC = -1;
 861     int               nmols;
 862     int               i;
 863     real*             qmol;
 864     FILE*             outf   = nullptr;
 865     FILE*             mcor   = nullptr;
 866     FILE*             fmj    = nullptr;
 867     FILE*             fmd    = nullptr;
 868     FILE*             fmjdsp = nullptr;
 869     FILE*             fcur   = nullptr;
 870     t_filenm          fnm[]  = { { efTPS, nullptr, nullptr, ffREAD }, /* this is for the topology */
 871                        { efNDX, nullptr, nullptr, ffOPTRD },
 872                        { efTRX, "-f", nullptr, ffREAD }, /* and this for the trajectory */
 873                        { efXVG, "-o", "current", ffWRITE },
 874                        { efXVG, "-caf", "caf", ffOPTWR },
 875                        { efXVG, "-dsp", "dsp", ffWRITE },
 876                        { efXVG, "-md", "md", ffWRITE },
 877                        { efXVG, "-mj", "mj", ffWRITE },
 878                        { efXVG, "-mc", "mc", ffOPTWR } };
 879
 880 #define NFILE asize(fnm)
 881
 882
 883     const char* desc[] = {
 884         "[THISMODULE] is a tool for calculating the current autocorrelation function, the "
 885         "correlation",
 886         "of the rotational and translational dipole moment of the system, and the resulting static",
 887         "dielectric constant. To obtain a reasonable result, the index group has to be neutral.",
 888         "Furthermore, the routine is capable of extracting the static conductivity from the "
 889         "current ",
 890         "autocorrelation function, if velocities are given. Additionally, an Einstein-Helfand fit ",
 891         "can be used to obtain the static conductivity.",
 892         "[PAR]",
 893         "The flag [TT]-caf[tt] is for the output of the current autocorrelation function and "
 894         "[TT]-mc[tt] writes the",
 895         "correlation of the rotational and translational part of the dipole moment in the "
 896         "corresponding",
 897         "file. However, this option is only available for trajectories containing velocities.",
 898         "Options [TT]-sh[tt] and [TT]-tr[tt] are responsible for the averaging and integration of "
 899         "the",
 900         "autocorrelation functions. Since averaging proceeds by shifting the starting point",
 901         "through the trajectory, the shift can be modified with [TT]-sh[tt] to enable the choice "
 902         "of uncorrelated",
 903         "starting points. Towards the end, statistical inaccuracy grows and integrating the",
 904         "correlation function only yields reliable values until a certain point, depending on",
 905         "the number of frames. The option [TT]-tr[tt] controls the region of the integral taken "
 906         "into account",
 907         "for calculating the static dielectric constant.",
 908         "[PAR]",
 909         "Option [TT]-temp[tt] sets the temperature required for the computation of the static "
 910         "dielectric constant.",
 911         "[PAR]",
 912         "Option [TT]-eps[tt] controls the dielectric constant of the surrounding medium for "
 913         "simulations using",
 914         "a Reaction Field or dipole corrections of the Ewald summation ([TT]-eps[tt]\\=0 "
 915         "corresponds to",
 916         "tin-foil boundary conditions).",
 917         "[PAR]",
 918         "[TT]-[no]nojump[tt] unfolds the coordinates to allow free diffusion. This is required to "
 919         "get a continuous",
 920         "translational dipole moment, required for the Einstein-Helfand fit. The results from the "
 921         "fit allow",
 922         "the determination of the dielectric constant for system of charged molecules. However, it "
 923         "is also possible to extract",
 924         "the dielectric constant from the fluctuations of the total dipole moment in folded "
 925         "coordinates. But this",
 926         "option has to be used with care, since only very short time spans fulfill the "
 927         "approximation that the density",
 928         "of the molecules is approximately constant and the averages are already converged. To be "
 929         "on the safe side,",
 930         "the dielectric constant should be calculated with the help of the Einstein-Helfand method "
 931         "for",
 932         "the translational part of the dielectric constant."
 933     };
 934
 935
 936     /* At first the arguments will be parsed and the system information processed */
 937     if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW, NFILE, fnm, asize(pa), pa,
 938                            asize(desc), desc, 0, nullptr, &oenv))
 939     {
 940         return 0;
 941     }
 942
 943     bACF = opt2bSet("-caf", NFILE, fnm);
 944     bINT = opt2bSet("-mc", NFILE, fnm);
 945
 946     read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &ePBC, nullptr, nullptr, box, TRUE);
 947
 948     indexfn = ftp2fn_null(efNDX, NFILE, fnm);
 949     snew(grpname, 1);
 950
 951     get_index(&(top.atoms), indexfn, 1, &isize, &index0, grpname);
 952
 953     flags = flags | TRX_READ_X | TRX_READ_V;
 954
 955     read_first_frame(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &fr, flags);
 956
 957     snew(mass2, top.atoms.nr);
 958     snew(qmol, top.atoms.nr);
 959
 960     precalc(top, mass2, qmol);
 961
 962
 963     snew(indexm, isize);
 964
 965     for (i = 0; i < isize; i++)
 966     {
 967         indexm[i] = index0[i];
 968     }
 969
 970     nmols = isize;
 971
 972
 973     index_atom2mol(&nmols, indexm, &top.mols);
 974
 975     if (fr.bV)
 976     {
 977         if (bACF)
 978         {
 979             outf = xvgropen(opt2fn("-caf", NFILE, fnm), "Current autocorrelation function",
 980                             output_env_get_xvgr_tlabel(oenv), "ACF (e nm/ps)\\S2", oenv);
 981             fprintf(outf, "# time\t acf\t average \t std.dev\n");
 982         }
 983         fcur = xvgropen(opt2fn("-o", NFILE, fnm), "Current", output_env_get_xvgr_tlabel(oenv),
 984                         "J(t) (e nm/ps)", oenv);
 985         fprintf(fcur, "# time\t Jx\t Jy \t J_z \n");
 986         if (bINT)
 987         {
 988             mcor = xvgropen(opt2fn("-mc", NFILE, fnm),
 989                             "M\\sD\\N - current  autocorrelation function",
 990                             output_env_get_xvgr_tlabel(oenv),
 991                             "< M\\sD\\N (0)\\c7\\CJ(t) >  (e nm/ps)\\S2", oenv);
 992             fprintf(mcor, "# time\t M_D(0) J(t) acf \t Integral acf\n");
 993         }
 994     }
 995
 996     fmj = xvgropen(opt2fn("-mj", NFILE, fnm), "Averaged translational part of M",
 997                    output_env_get_xvgr_tlabel(oenv), "< M\\sJ\\N > (enm)", oenv);
 998     fprintf(fmj, "# time\t x\t y \t z \t average of M_J^2 \t std.dev\n");
 999     fmd = xvgropen(opt2fn("-md", NFILE, fnm), "Averaged rotational part of M",
1000                    output_env_get_xvgr_tlabel(oenv), "< M\\sD\\N > (enm)", oenv);
1001     fprintf(fmd, "# time\t x\t y \t z \t average of M_D^2 \t std.dev\n");
1002
1003     fmjdsp = xvgropen(
1004             opt2fn("-dsp", NFILE, fnm), "MSD of the squared translational dipole moment M",
1005             output_env_get_xvgr_tlabel(oenv),
1006             "<|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);
1007
1008
1009     /* System information is read and prepared, dielectric() processes the frames
1010      * and calculates the requested quantities */
1011
1012     dielectric(fmj, fmd, outf, fcur, mcor, fmjdsp, bNoJump, bACF, bINT, ePBC, top, fr, temp, bfit, efit,
1013                bvit, evit, status, isize, nmols, nshift, index0, indexm, mass2, qmol, eps_rf, oenv);
1014
1015     xvgrclose(fmj);
1016     xvgrclose(fmd);
1017     xvgrclose(fmjdsp);
1018     if (fr.bV)
1019     {
1020         if (bACF)
1021         {
1022             xvgrclose(outf);
1023         }
1024         xvgrclose(fcur);
1025         if (bINT)
1026         {
1027             xvgrclose(mcor);
1028         }
1029     }
1030
1031     return 0;
1032 }