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