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