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