src/gromacs/gmxana/gmx_current.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *                        VERSION 3.0
   9  *
  10  * Copyright (c) 1991-2001
  11  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  12  * University of Groningen, The Netherlands
  13  *
  14  * This program is free software; you can redistribute it and/or
  15  *
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * Do check out http://www.gromacs.org , or mail us at gromacs@gromacs.org .
  31  *
  32  * And Hey:
  33  * Gyas ROwers Mature At Cryogenic Speed
  34  *
  35  * finished FD 09/07/08
  36  *
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42
  43 #include "statutil.h"
  44 #include "typedefs.h"
  45 #include "smalloc.h"
  46 #include "vec.h"
  47 #include "copyrite.h"
  48 #include "statutil.h"
  49 #include "tpxio.h"
  50 #include "xvgr.h"
  51 #include "rmpbc.h"
  52 #include "pbc.h"
  53 #include "physics.h"
  54 #include "index.h"
  55 #include "gmx_statistics.h"
  56 #include "gmx_ana.h"
  57 #include "macros.h"
  58
  59 #define SQR(x) (pow(x, 2.0))
  60 #define EPSI0 (EPSILON0*E_CHARGE*E_CHARGE*AVOGADRO/(KILO*NANO)) /* EPSILON0 in SI units */
  61
  62
  63 static void index_atom2mol(int *n, int *index, t_block *mols)
  64 {
  65     int nat, i, nmol, mol, j;
  66
  67     nat  = *n;
  68     i    = 0;
  69     nmol = 0;
  70     mol  = 0;
  71     while (i < nat)
  72     {
  73         while (index[i] > mols->index[mol])
  74         {
  75             mol++;
  76             if (mol >= mols->nr)
  77             {
  78                 gmx_fatal(FARGS, "Atom index out of range: %d", index[i]+1);
  79             }
  80         }
  81         for (j = mols->index[mol]; j < mols->index[mol+1]; j++)
  82         {
  83             if (i >= nat || index[i] != j)
  84             {
  85                 gmx_fatal(FARGS, "The index group does not consist of whole molecules");
  86             }
  87             i++;
  88         }
  89         index[nmol++] = mol;
  90     }
  91
  92     fprintf(stderr, "\nSplit group of %d atoms into %d molecules\n", nat, nmol);
  93
  94     *n = nmol;
  95 }
  96
  97 static gmx_bool precalc(t_topology top, real mass2[], real qmol[])
  98 {
  99
 100     real     mtot;
 101     real     qtot;
 102     real     qall;
 103     int      i, j, k, l;
 104     int      ai, ci;
 105     gmx_bool bNEU;
 106     ai   = 0;
 107     ci   = 0;
 108     qall = 0.0;
 109
 110
 111
 112     for (i = 0; i < top.mols.nr; i++)
 113     {
 114         k    = top.mols.index[i];
 115         l    = top.mols.index[i+1];
 116         mtot = 0.0;
 117         qtot = 0.0;
 118
 119         for (j = k; j < l; j++)
 120         {
 121             mtot += top.atoms.atom[j].m;
 122             qtot += top.atoms.atom[j].q;
 123         }
 124
 125         for (j = k; j < l; j++)
 126         {
 127             top.atoms.atom[j].q -= top.atoms.atom[j].m*qtot/mtot;
 128             mass2[j]             = top.atoms.atom[j].m/mtot;
 129             qmol[j]              = qtot;
 130         }
 131
 132
 133         qall += qtot;
 134
 135         if (qtot < 0.0)
 136         {
 137             ai++;
 138         }
 139         if (qtot > 0.0)
 140         {
 141             ci++;
 142         }
 143     }
 144
 145     if (abs(qall) > 0.01)
 146     {
 147         printf("\n\nSystem not neutral (q=%f) will not calculate translational part of the dipole moment.\n", qall);
 148         bNEU = FALSE;
 149     }
 150     else
 151     {
 152         bNEU = TRUE;
 153     }
 154
 155     return bNEU;
 156
 157 }
 158
 159 static void remove_jump(matrix box, int natoms, rvec xp[], rvec x[])
 160 {
 161
 162     rvec hbox;
 163     int  d, i, m;
 164
 165     for (d = 0; d < DIM; d++)
 166     {
 167         hbox[d] = 0.5*box[d][d];
 168     }
 169     for (i = 0; i < natoms; i++)
 170     {
 171         for (m = DIM-1; m >= 0; m--)
 172         {
 173             while (x[i][m]-xp[i][m] <= -hbox[m])
 174             {
 175                 for (d = 0; d <= m; d++)
 176                 {
 177                     x[i][d] += box[m][d];
 178                 }
 179             }
 180             while (x[i][m]-xp[i][m] > hbox[m])
 181             {
 182                 for (d = 0; d <= m; d++)
 183                 {
 184                     x[i][d] -= box[m][d];
 185                 }
 186             }
 187         }
 188     }
 189 }
 190
 191 static void calc_mj(t_topology top, int ePBC, matrix box, gmx_bool bNoJump, int isize, int index0[], \
 192                     rvec fr[], rvec mj, real mass2[], real qmol[])
 193 {
 194
 195     int   i, j, k, l;
 196     rvec  tmp;
 197     rvec  center;
 198     rvec  mt1, mt2;
 199     t_pbc pbc;
 200
 201
 202     calc_box_center(ecenterRECT, box, center);
 203
 204     if (!bNoJump)
 205     {
 206         set_pbc(&pbc, ePBC, box);
 207     }
 208
 209     clear_rvec(tmp);
 210
 211
 212     for (i = 0; i < isize; i++)
 213     {
 214         clear_rvec(mt1);
 215         clear_rvec(mt2);
 216         k = top.mols.index[index0[i]];
 217         l = top.mols.index[index0[i+1]];
 218         for (j = k; j < l; j++)
 219         {
 220             svmul(mass2[j], fr[j], tmp);
 221             rvec_inc(mt1, tmp);
 222         }
 223
 224         if (bNoJump)
 225         {
 226             svmul(qmol[k], mt1, mt1);
 227         }
 228         else
 229         {
 230             pbc_dx(&pbc, mt1, center, mt2);
 231             svmul(qmol[k], mt2, mt1);
 232         }
 233
 234         rvec_inc(mj, mt1);
 235
 236     }
 237
 238 }
 239
 240
 241 static real calceps(real prefactor, real md2, real mj2, real cor, real eps_rf, gmx_bool bCOR)
 242 {
 243
 244     /* bCOR determines if the correlation is computed via
 245      * static properties (FALSE) or the correlation integral (TRUE)
 246      */
 247
 248     real epsilon = 0.0;
 249
 250
 251     if (bCOR)
 252     {
 253         epsilon = md2-2.0*cor+mj2;
 254     }
 255     else
 256     {
 257         epsilon = md2+mj2+2.0*cor;
 258     }
 259
 260     if (eps_rf == 0.0)
 261     {
 262         epsilon = 1.0+prefactor*epsilon;
 263
 264     }
 265     else
 266     {
 267         epsilon  = 2.0*eps_rf+1.0+2.0*eps_rf*prefactor*epsilon;
 268         epsilon /= 2.0*eps_rf+1.0-prefactor*epsilon;
 269     }
 270
 271
 272     return epsilon;
 273
 274 }
 275
 276
 277 static real calc_cacf(FILE *fcacf, real prefactor, real cacf[], real time[], int nfr, int vfr[], int ei, int nshift)
 278 {
 279
 280     int  i;
 281     real corint;
 282     real deltat = 0.0;
 283     real rfr;
 284     real sigma     = 0.0;
 285     real sigma_ret = 0.0;
 286     corint = 0.0;
 287
 288     if (nfr > 1)
 289     {
 290         i = 0;
 291
 292         while (i < nfr)
 293         {
 294
 295             rfr      = (real) (nfr/nshift-i/nshift);
 296             cacf[i] /= rfr;
 297
 298             if (time[vfr[i]] <= time[vfr[ei]])
 299             {
 300                 sigma_ret = sigma;
 301             }
 302
 303             fprintf(fcacf, "%.3f\t%10.6g\t%10.6g\n", time[vfr[i]], cacf[i], sigma);
 304
 305             if ((i+1) < nfr)
 306             {
 307                 deltat = (time[vfr[i+1]]-time[vfr[i]]);
 308             }
 309             corint = 2.0*deltat*cacf[i]*prefactor;
 310             if (i == 0 || (i+1) == nfr)
 311             {
 312                 corint *= 0.5;
 313             }
 314             sigma += corint;
 315
 316             i++;
 317         }
 318
 319     }
 320     else
 321     {
 322         printf("Too less points.\n");
 323     }
 324
 325     return sigma_ret;
 326
 327 }
 328
 329 static void calc_mjdsp(FILE *fmjdsp, real vol, real temp, real prefactor, rvec mjdsp[], real dsp2[], real time[], int nfr, real refr[])
 330 {
 331
 332     int     i;
 333     real    rtmp;
 334     real    rfr;
 335
 336
 337     fprintf(fmjdsp, "#Prefactor fit E-H: 1 / 6.0*V*k_B*T: %g\n", prefactor);
 338
 339
 340
 341     for (i = 0; i < nfr; i++)
 342     {
 343
 344         if (refr[i] != 0.0)
 345         {
 346             dsp2[i] *= prefactor/refr[i];
 347             fprintf(fmjdsp, "%.3f\t%10.6g\n", time[i], dsp2[i]);
 348         }
 349
 350
 351     }
 352
 353 }
 354
 355
 356 static void dielectric(FILE *fmj, FILE *fmd, FILE *outf, FILE *fcur, FILE *mcor,
 357                        FILE *fmjdsp, gmx_bool bNoJump, gmx_bool bACF, gmx_bool bINT,
 358                        int ePBC, t_topology top, t_trxframe fr, real temp,
 359                        real trust, real bfit, real efit, real bvit, real evit,
 360                        t_trxstatus *status, int isize, int nmols, int nshift,
 361                        atom_id *index0, int indexm[], real mass2[],
 362                        real qmol[], real eps_rf, const output_env_t oenv)
 363 {
 364     int       i, j, k, l, f;
 365     int       valloc, nalloc, nfr, nvfr, m, itrust = 0;
 366     int       vshfr;
 367     real     *xshfr       = NULL;
 368     int      *vfr         = NULL;
 369     real      refr        = 0.0;
 370     real      rfr         = 0.0;
 371     real     *cacf        = NULL;
 372     real     *time        = NULL;
 373     real     *djc         = NULL;
 374     real      corint      = 0.0;
 375     real      prefactorav = 0.0;
 376     real      prefactor   = 0.0;
 377     real      volume;
 378     real      volume_av = 0.0;
 379     real      dk_s, dk_d, dk_f;
 380     real      dm_s, dm_d;
 381     real      mj    = 0.0;
 382     real      mj2   = 0.0;
 383     real      mjd   = 0.0;
 384     real      mjdav = 0.0;
 385     real      md2   = 0.0;
 386     real      mdav2 = 0.0;
 387     real      sgk;
 388     rvec      mja_tmp;
 389     rvec      mjd_tmp;
 390     rvec      mdvec;
 391     rvec     *mu    = NULL;
 392     rvec     *xp    = NULL;
 393     rvec     *v0    = NULL;
 394     rvec     *mjdsp = NULL;
 395     real     *dsp2  = NULL;
 396     real      t0;
 397     real      rtmp;
 398     real      qtmp;
 399
 400
 401
 402     rvec   tmp;
 403     rvec  *mtrans = NULL;
 404
 405     /*
 406      * Variables for the least-squares fit for Einstein-Helfand and Green-Kubo
 407      */
 408
 409     int          bi, ei, ie, ii;
 410     real         rest  = 0.0;
 411     real         sigma = 0.0;
 412     real         malt  = 0.0;
 413     real         err   = 0.0;
 414     real        *xfit;
 415     real        *yfit;
 416     gmx_rmpbc_t  gpbc = NULL;
 417
 418     /*
 419      * indices for EH
 420      */
 421
 422     ei = 0;
 423     bi = 0;
 424
 425     /*
 426      * indices for GK
 427      */
 428
 429     ii  = 0;
 430     ie  = 0;
 431     t0  = 0;
 432     sgk = 0.0;
 433
 434
 435     /* This is the main loop over frames */
 436
 437
 438     nfr = 0;
 439
 440
 441     nvfr   = 0;
 442     vshfr  = 0;
 443     nalloc = 0;
 444     valloc = 0;
 445
 446     clear_rvec(mja_tmp);
 447     clear_rvec(mjd_tmp);
 448     clear_rvec(mdvec);
 449     clear_rvec(tmp);
 450     gpbc = gmx_rmpbc_init(&top.idef, ePBC, fr.natoms, fr.box);
 451
 452     do
 453     {
 454
 455         refr = (real)(nfr+1);
 456
 457         if (nfr >= nalloc)
 458         {
 459             nalloc += 100;
 460             srenew(time, nalloc);
 461             srenew(mu, nalloc);
 462             srenew(mjdsp, nalloc);
 463             srenew(dsp2, nalloc);
 464             srenew(mtrans, nalloc);
 465             srenew(xshfr, nalloc);
 466
 467             for (i = nfr; i < nalloc; i++)
 468             {
 469                 clear_rvec(mjdsp[i]);
 470                 clear_rvec(mu[i]);
 471                 clear_rvec(mtrans[i]);
 472                 dsp2[i]  = 0.0;
 473                 xshfr[i] = 0.0;
 474             }
 475         }
 476
 477
 478
 479         if (nfr == 0)
 480         {
 481             t0 = fr.time;
 482
 483         }
 484
 485         time[nfr] = fr.time-t0;
 486
 487         if (time[nfr] <= bfit)
 488         {
 489             bi = nfr;
 490         }
 491         if (time[nfr] <= efit)
 492         {
 493             ei = nfr;
 494         }
 495
 496         if (bNoJump)
 497         {
 498
 499             if (xp)
 500             {
 501                 remove_jump(fr.box, fr.natoms, xp, fr.x);
 502             }
 503             else
 504             {
 505                 snew(xp, fr.natoms);
 506             }
 507
 508             for (i = 0; i < fr.natoms; i++)
 509             {
 510                 copy_rvec(fr.x[i], xp[i]);
 511             }
 512
 513         }
 514
 515         gmx_rmpbc_trxfr(gpbc, &fr);
 516
 517         calc_mj(top, ePBC, fr.box, bNoJump, nmols, indexm, fr.x, mtrans[nfr], mass2, qmol);
 518
 519         for (i = 0; i < isize; i++)
 520         {
 521             j = index0[i];
 522             svmul(top.atoms.atom[j].q, fr.x[j], fr.x[j]);
 523             rvec_inc(mu[nfr], fr.x[j]);
 524         }
 525
 526         /*if(mod(nfr,nshift)==0){*/
 527         if (nfr%nshift == 0)
 528         {
 529             for (j = nfr; j >= 0; j--)
 530             {
 531                 rvec_sub(mtrans[nfr], mtrans[j], tmp);
 532                 dsp2[nfr-j]  += norm2(tmp);
 533                 xshfr[nfr-j] += 1.0;
 534             }
 535         }
 536
 537         if (fr.bV)
 538         {
 539             if (nvfr >= valloc)
 540             {
 541                 valloc += 100;
 542                 srenew(vfr, valloc);
 543                 if (bINT)
 544                 {
 545                     srenew(djc, valloc);
 546                 }
 547                 srenew(v0, valloc);
 548                 if (bACF)
 549                 {
 550                     srenew(cacf, valloc);
 551                 }
 552             }
 553             if (time[nfr] <= bvit)
 554             {
 555                 ii = nvfr;
 556             }
 557             if (time[nfr] <= evit)
 558             {
 559                 ie = nvfr;
 560             }
 561             vfr[nvfr] = nfr;
 562             clear_rvec(v0[nvfr]);
 563             if (bACF)
 564             {
 565                 cacf[nvfr] = 0.0;
 566             }
 567             if (bINT)
 568             {
 569                 djc[nvfr] = 0.0;
 570             }
 571             for (i = 0; i < isize; i++)
 572             {
 573                 j = index0[i];
 574                 svmul(mass2[j], fr.v[j], fr.v[j]);
 575                 svmul(qmol[j], fr.v[j], fr.v[j]);
 576                 rvec_inc(v0[nvfr], fr.v[j]);
 577             }
 578
 579             fprintf(fcur, "%.3f\t%.6f\t%.6f\t%.6f\n", time[nfr], v0[nfr][XX], v0[nfr][YY], v0[nfr][ZZ]);
 580             if (bACF || bINT)
 581             {
 582                 /*if(mod(nvfr,nshift)==0){*/
 583                 if (nvfr%nshift == 0)
 584                 {
 585                     for (j = nvfr; j >= 0; j--)
 586                     {
 587                         if (bACF)
 588                         {
 589                             cacf[nvfr-j] += iprod(v0[nvfr], v0[j]);
 590                         }
 591                         if (bINT)
 592                         {
 593                             djc[nvfr-j] += iprod(mu[vfr[j]], v0[nvfr]);
 594                         }
 595                     }
 596                     vshfr++;
 597                 }
 598             }
 599             nvfr++;
 600         }
 601
 602         volume     = det(fr.box);
 603         volume_av += volume;
 604
 605         rvec_inc(mja_tmp, mtrans[nfr]);
 606         mjd += iprod(mu[nfr], mtrans[nfr]);
 607         rvec_inc(mdvec, mu[nfr]);
 608
 609         mj2 += iprod(mtrans[nfr], mtrans[nfr]);
 610         md2 += iprod(mu[nfr], mu[nfr]);
 611
 612         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);
 613         fprintf(fmd, "%.3f\t%8.5f\t%8.5f\t%8.5f\t%8.5f\t%8.5f\n",    \
 614                 time[nfr], mu[nfr][XX], mu[nfr][YY], mu[nfr][ZZ], md2/refr, norm(mdvec)/refr);
 615
 616         nfr++;
 617
 618     }
 619     while (read_next_frame(oenv, status, &fr));
 620
 621     gmx_rmpbc_done(gpbc);
 622
 623     volume_av /= refr;
 624
 625     prefactor  = 1.0;
 626     prefactor /= 3.0*EPSILON0*volume_av*BOLTZ*temp;
 627
 628
 629     prefactorav  = E_CHARGE*E_CHARGE;
 630     prefactorav /= volume_av*BOLTZMANN*temp*NANO*6.0;
 631
 632     fprintf(stderr, "Prefactor fit E-H: 1 / 6.0*V*k_B*T: %g\n", prefactorav);
 633
 634     calc_mjdsp(fmjdsp, volume_av, temp, prefactorav, mjdsp, dsp2, time, nfr, xshfr);
 635
 636     /*
 637      * Now we can average and calculate the correlation functions
 638      */
 639
 640
 641     mj2 /= refr;
 642     mjd /= refr;
 643     md2 /= refr;
 644
 645     svmul(1.0/refr, mdvec, mdvec);
 646     svmul(1.0/refr, mja_tmp, mja_tmp);
 647
 648     mdav2 = norm2(mdvec);
 649     mj    = norm2(mja_tmp);
 650     mjdav = iprod(mdvec, mja_tmp);
 651
 652
 653     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);
 654     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);
 655
 656     if (v0 != NULL)
 657     {
 658         trust *= (real) nvfr;
 659         itrust = (int) trust;
 660
 661         if (bINT)
 662         {
 663
 664             printf("\nCalculating M_D - current correlation integral ... \n");
 665             corint = calc_cacf(mcor, prefactorav/EPSI0, djc, time, nvfr, vfr, ie, nshift);
 666
 667         }
 668
 669         if (bACF)
 670         {
 671
 672             printf("\nCalculating current autocorrelation ... \n");
 673             sgk = calc_cacf(outf, prefactorav/PICO, cacf, time, nvfr, vfr, ie, nshift);
 674
 675             if (ie > ii)
 676             {
 677
 678                 snew(xfit, ie-ii+1);
 679                 snew(yfit, ie-ii+1);
 680
 681                 for (i = ii; i <= ie; i++)
 682                 {
 683
 684                     xfit[i-ii] = log(time[vfr[i]]);
 685                     rtmp       = fabs(cacf[i]);
 686                     yfit[i-ii] = log(rtmp);
 687
 688                 }
 689
 690                 lsq_y_ax_b(ie-ii, xfit, yfit, &sigma, &malt, &err, &rest);
 691
 692                 malt = exp(malt);
 693
 694                 sigma += 1.0;
 695
 696                 malt *= prefactorav*2.0e12/sigma;
 697
 698                 sfree(xfit);
 699                 sfree(yfit);
 700
 701             }
 702         }
 703     }
 704
 705
 706     /* Calculation of the dielectric constant */
 707
 708     fprintf(stderr, "\n********************************************\n");
 709     dk_s = calceps(prefactor, md2, mj2, mjd, eps_rf, FALSE);
 710     fprintf(stderr, "\nAbsolute values:\n epsilon=%f\n", dk_s);
 711     fprintf(stderr, " <M_D^2> , <M_J^2>, <(M_J*M_D)^2>:  (%f, %f, %f)\n\n", md2, mj2, mjd);
 712     fprintf(stderr, "********************************************\n");
 713
 714
 715     dk_f = calceps(prefactor, md2-mdav2, mj2-mj, mjd-mjdav, eps_rf, FALSE);
 716     fprintf(stderr, "\n\nFluctuations:\n epsilon=%f\n\n", dk_f);
 717     fprintf(stderr, "\n deltaM_D , deltaM_J, deltaM_JD:  (%f, %f, %f)\n", md2-mdav2, mj2-mj, mjd-mjdav);
 718     fprintf(stderr, "\n********************************************\n");
 719     if (bINT)
 720     {
 721         dk_d = calceps(prefactor, md2-mdav2, mj2-mj, corint, eps_rf, TRUE);
 722         fprintf(stderr, "\nStatic dielectric constant using integral and fluctuations: %f\n", dk_d);
 723         fprintf(stderr, "\n < M_JM_D > via integral:  %.3f\n", -1.0*corint);
 724     }
 725
 726     fprintf(stderr, "\n***************************************************");
 727     fprintf(stderr, "\n\nAverage volume V=%f nm^3 at T=%f K\n", volume_av, temp);
 728     fprintf(stderr, "and corresponding refactor 1.0 / 3.0*V*k_B*T*EPSILON_0: %f \n", prefactor);
 729
 730
 731
 732     if (bACF)
 733     {
 734         fprintf(stderr, "Integral and integrated fit to the current acf yields at t=%f:\n", time[vfr[ii]]);
 735         fprintf(stderr, "sigma=%8.3f (pure integral: %.3f)\n", sgk-malt*pow(time[vfr[ii]], sigma), sgk);
 736     }
 737
 738     if (ei > bi)
 739     {
 740         fprintf(stderr, "\nStart fit at %f ps (%f).\n", time[bi], bfit);
 741         fprintf(stderr, "End fit at %f ps (%f).\n\n", time[ei], efit);
 742
 743         snew(xfit, ei-bi+1);
 744         snew(yfit, ei-bi+1);
 745
 746         for (i = bi; i <= ei; i++)
 747         {
 748             xfit[i-bi] = time[i];
 749             yfit[i-bi] = dsp2[i];
 750         }
 751
 752         lsq_y_ax_b(ei-bi, xfit, yfit, &sigma, &malt, &err, &rest);
 753
 754         sigma *= 1e12;
 755         dk_d   = calceps(prefactor, md2, 0.5*malt/prefactorav, corint, eps_rf, TRUE);
 756
 757         fprintf(stderr, "Einstein-Helfand fit to the MSD of the translational dipole moment yields:\n\n");
 758         fprintf(stderr, "sigma=%.4f\n", sigma);
 759         fprintf(stderr, "translational fraction of M^2: %.4f\n", 0.5*malt/prefactorav);
 760         fprintf(stderr, "Dielectric constant using EH: %.4f\n", dk_d);
 761
 762         sfree(xfit);
 763         sfree(yfit);
 764     }
 765     else
 766     {
 767         fprintf(stderr, "Too less points for a fit.\n");
 768     }
 769
 770
 771     if (v0 != NULL)
 772     {
 773         sfree(v0);
 774     }
 775     if (bACF)
 776     {
 777         sfree(cacf);
 778     }
 779     if (bINT)
 780     {
 781         sfree(djc);
 782     }
 783
 784     sfree(time);
 785
 786
 787     sfree(mjdsp);
 788     sfree(mu);
 789 }
 790
 791
 792
 793 int gmx_current(int argc, char *argv[])
 794 {
 795
 796     static int             nshift  = 1000;
 797     static real            temp    = 300.0;
 798     static real            trust   = 0.25;
 799     static real            eps_rf  = 0.0;
 800     static gmx_bool        bNoJump = TRUE;
 801     static real            bfit    = 100.0;
 802     static real            bvit    = 0.5;
 803     static real            efit    = 400.0;
 804     static real            evit    = 5.0;
 805     t_pargs                pa[]    = {
 806         { "-sh", FALSE, etINT, {&nshift},
 807           "Shift of the frames for averaging the correlation functions and the mean-square displacement."},
 808         { "-nojump", FALSE, etBOOL, {&bNoJump},
 809           "Removes jumps of atoms across the box."},
 810         { "-eps", FALSE, etREAL, {&eps_rf},
 811           "Dielectric constant of the surrounding medium. The value zero corresponds to infinity (tin-foil boundary conditions)."},
 812         { "-bfit", FALSE, etREAL, {&bfit},
 813           "Begin of the fit of the straight line to the MSD of the translational fraction of the dipole moment."},
 814         { "-efit", FALSE, etREAL, {&efit},
 815           "End of the fit of the straight line to the MSD of the translational fraction of the dipole moment."},
 816         { "-bvit", FALSE, etREAL, {&bvit},
 817           "Begin of the fit of the current autocorrelation function to a*t^b."},
 818         { "-evit", FALSE, etREAL, {&evit},
 819           "End of the fit of the current autocorrelation function to a*t^b."},
 820         { "-tr", FALSE, etREAL, {&trust},
 821           "Fraction of the trajectory taken into account for the integral."},
 822         { "-temp", FALSE, etREAL, {&temp},
 823           "Temperature for calculating epsilon."}
 824     };
 825
 826     output_env_t           oenv;
 827     t_topology             top;
 828     char                   title[STRLEN];
 829     char                 **grpname = NULL;
 830     const char            *indexfn;
 831     t_trxframe             fr;
 832     real                  *mass2 = NULL;
 833     rvec                  *xtop, *vtop;
 834     matrix                 box;
 835     atom_id               *index0;
 836     int                   *indexm = NULL;
 837     int                    isize;
 838     t_trxstatus           *status;
 839     int                    flags = 0;
 840     gmx_bool               bTop;
 841     gmx_bool               bNEU;
 842     gmx_bool               bACF;
 843     gmx_bool               bINT;
 844     int                    ePBC = -1;
 845     int                    natoms;
 846     int                    nmols;
 847     int                    i, j, k = 0, l;
 848     int                    step;
 849     real                   t;
 850     real                   lambda;
 851     real                  *qmol;
 852     FILE                  *outf   = NULL;
 853     FILE                  *outi   = NULL;
 854     FILE                  *tfile  = NULL;
 855     FILE                  *mcor   = NULL;
 856     FILE                  *fmj    = NULL;
 857     FILE                  *fmd    = NULL;
 858     FILE                  *fmjdsp = NULL;
 859     FILE                  *fcur   = NULL;
 860     t_filenm               fnm[]  = {
 861         { efTPS,  NULL,  NULL, ffREAD }, /* this is for the topology */
 862         { efNDX, NULL, NULL, ffOPTRD },
 863         { efTRX, "-f", NULL, ffREAD },   /* and this for the trajectory */
 864         { efXVG, "-o", "current.xvg", ffWRITE },
 865         { efXVG, "-caf", "caf.xvg", ffOPTWR },
 866         { efXVG, "-dsp", "dsp.xvg", ffWRITE },
 867         { efXVG, "-md", "md.xvg", ffWRITE },
 868         { efXVG, "-mj", "mj.xvg", ffWRITE},
 869         { efXVG, "-mc", "mc.xvg", ffOPTWR }
 870     };
 871
 872 #define NFILE asize(fnm)
 873
 874
 875     const char *desc[] = {
 876         "[TT]g_current[tt] is a tool for calculating the current autocorrelation function, the correlation",
 877         "of the rotational and translational dipole moment of the system, and the resulting static",
 878         "dielectric constant. To obtain a reasonable result, the index group has to be neutral.",
 879         "Furthermore, the routine is capable of extracting the static conductivity from the current ",
 880         "autocorrelation function, if velocities are given. Additionally, an Einstein-Helfand fit ",
 881         "can be used to obtain the static conductivity."
 882         "[PAR]",
 883         "The flag [TT]-caf[tt] is for the output of the current autocorrelation function and [TT]-mc[tt] writes the",
 884         "correlation of the rotational and translational part of the dipole moment in the corresponding",
 885         "file. However, this option is only available for trajectories containing velocities.",
 886         "Options [TT]-sh[tt] and [TT]-tr[tt] are responsible for the averaging and integration of the",
 887         "autocorrelation functions. Since averaging proceeds by shifting the starting point",
 888         "through the trajectory, the shift can be modified with [TT]-sh[tt] to enable the choice of uncorrelated",
 889         "starting points. Towards the end, statistical inaccuracy grows and integrating the",
 890         "correlation function only yields reliable values until a certain point, depending on",
 891         "the number of frames. The option [TT]-tr[tt] controls the region of the integral taken into account",
 892         "for calculating the static dielectric constant.",
 893         "[PAR]",
 894         "Option [TT]-temp[tt] sets the temperature required for the computation of the static dielectric constant.",
 895         "[PAR]",
 896         "Option [TT]-eps[tt] controls the dielectric constant of the surrounding medium for simulations using",
 897         "a Reaction Field or dipole corrections of the Ewald summation ([TT]-eps[tt]=0 corresponds to",
 898         "tin-foil boundary conditions).",
 899         "[PAR]",
 900         "[TT]-[no]nojump[tt] unfolds the coordinates to allow free diffusion. This is required to get a continuous",
 901         "translational dipole moment, required for the Einstein-Helfand fit. The results from the fit allow",
 902         "the determination of the dielectric constant for system of charged molecules. However, it is also possible to extract",
 903         "the dielectric constant from the fluctuations of the total dipole moment in folded coordinates. But this",
 904         "option has to be used with care, since only very short time spans fulfill the approximation that the density",
 905         "of the molecules is approximately constant and the averages are already converged. To be on the safe side,",
 906         "the dielectric constant should be calculated with the help of the Einstein-Helfand method for",
 907         "the translational part of the dielectric constant."
 908     };
 909
 910
 911     /* At first the arguments will be parsed and the system information processed */
 912     parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW,
 913                       NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv);
 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     ffclose(fmj);
 997     ffclose(fmd);
 998     ffclose(fmjdsp);
 999     if (bACF)
1000     {
1001         ffclose(outf);
1002     }
1003     if (bINT)
1004     {
1005         ffclose(mcor);
1006     }
1007
1008     thanx(stderr);
1009
1010     return 0;
1011 }