src/gromacs/gmxana/gmx_current.c

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