src/gromacs/gmxana/gmx_current.c

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