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