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