2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
41 #include "gromacs/commandline/pargs.h"
43 #include "gromacs/utility/smalloc.h"
44 #include "gromacs/math/vec.h"
45 #include "gromacs/fileio/tpxio.h"
46 #include "gromacs/fileio/trxio.h"
47 #include "gromacs/fileio/xvgr.h"
52 #include "gromacs/statistics/statistics.h"
56 #include "gromacs/utility/fatalerror.h"
58 #define SQR(x) (pow(x, 2.0))
59 #define EPSI0 (EPSILON0*E_CHARGE*E_CHARGE*AVOGADRO/(KILO*NANO)) /* EPSILON0 in SI units */
61 static void index_atom2mol(int *n, int *index, t_block *mols)
63 int nat, i, nmol, mol, j;
71 while (index[i] > mols->index[mol])
76 gmx_fatal(FARGS, "Atom index out of range: %d", index[i]+1);
79 for (j = mols->index[mol]; j < mols->index[mol+1]; j++)
81 if (i >= nat || index[i] != j)
83 gmx_fatal(FARGS, "The index group does not consist of whole molecules");
90 fprintf(stderr, "\nSplit group of %d atoms into %d molecules\n", nat, nmol);
95 static gmx_bool precalc(t_topology top, real mass2[], real qmol[])
110 for (i = 0; i < top.mols.nr; i++)
112 k = top.mols.index[i];
113 l = top.mols.index[i+1];
117 for (j = k; j < l; j++)
119 mtot += top.atoms.atom[j].m;
120 qtot += top.atoms.atom[j].q;
123 for (j = k; j < l; j++)
125 top.atoms.atom[j].q -= top.atoms.atom[j].m*qtot/mtot;
126 mass2[j] = top.atoms.atom[j].m/mtot;
143 if (abs(qall) > 0.01)
145 printf("\n\nSystem not neutral (q=%f) will not calculate translational part of the dipole moment.\n", qall);
157 static void remove_jump(matrix box, int natoms, rvec xp[], rvec x[])
163 for (d = 0; d < DIM; d++)
165 hbox[d] = 0.5*box[d][d];
167 for (i = 0; i < natoms; i++)
169 for (m = DIM-1; m >= 0; m--)
171 while (x[i][m]-xp[i][m] <= -hbox[m])
173 for (d = 0; d <= m; d++)
175 x[i][d] += box[m][d];
178 while (x[i][m]-xp[i][m] > hbox[m])
180 for (d = 0; d <= m; d++)
182 x[i][d] -= box[m][d];
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[])
200 calc_box_center(ecenterRECT, box, center);
204 set_pbc(&pbc, ePBC, box);
210 for (i = 0; i < isize; i++)
214 k = top.mols.index[index0[i]];
215 l = top.mols.index[index0[i+1]];
216 for (j = k; j < l; j++)
218 svmul(mass2[j], fr[j], tmp);
224 svmul(qmol[k], mt1, mt1);
228 pbc_dx(&pbc, mt1, center, mt2);
229 svmul(qmol[k], mt2, mt1);
239 static real calceps(real prefactor, real md2, real mj2, real cor, real eps_rf, gmx_bool bCOR)
242 /* bCOR determines if the correlation is computed via
243 * static properties (FALSE) or the correlation integral (TRUE)
251 epsilon = md2-2.0*cor+mj2;
255 epsilon = md2+mj2+2.0*cor;
260 epsilon = 1.0+prefactor*epsilon;
265 epsilon = 2.0*eps_rf+1.0+2.0*eps_rf*prefactor*epsilon;
266 epsilon /= 2.0*eps_rf+1.0-prefactor*epsilon;
275 static real calc_cacf(FILE *fcacf, real prefactor, real cacf[], real time[], int nfr, int vfr[], int ei, int nshift)
283 real sigma_ret = 0.0;
293 rfr = (real) (nfr/nshift-i/nshift);
296 if (time[vfr[i]] <= time[vfr[ei]])
301 fprintf(fcacf, "%.3f\t%10.6g\t%10.6g\n", time[vfr[i]], cacf[i], sigma);
305 deltat = (time[vfr[i+1]]-time[vfr[i]]);
307 corint = 2.0*deltat*cacf[i]*prefactor;
308 if (i == 0 || (i+1) == nfr)
320 printf("Too less points.\n");
327 static void calc_mjdsp(FILE *fmjdsp, real prefactor, real dsp2[], real time[], int nfr, real refr[])
335 fprintf(fmjdsp, "#Prefactor fit E-H: 1 / 6.0*V*k_B*T: %g\n", prefactor);
339 for (i = 0; i < nfr; i++)
344 dsp2[i] *= prefactor/refr[i];
345 fprintf(fmjdsp, "%.3f\t%10.6g\n", time[i], dsp2[i]);
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)
363 int valloc, nalloc, nfr, nvfr, m, itrust = 0;
373 real prefactorav = 0.0;
374 real prefactor = 0.0;
376 real volume_av = 0.0;
377 real dk_s, dk_d, dk_f;
404 * Variables for the least-squares fit for Einstein-Helfand and Green-Kubo
414 gmx_rmpbc_t gpbc = NULL;
433 /* This is the main loop over frames */
448 gpbc = gmx_rmpbc_init(&top.idef, ePBC, fr.natoms);
453 refr = (real)(nfr+1);
458 srenew(time, nalloc);
460 srenew(mjdsp, nalloc);
461 srenew(dsp2, nalloc);
462 srenew(mtrans, nalloc);
463 srenew(xshfr, nalloc);
465 for (i = nfr; i < nalloc; i++)
467 clear_rvec(mjdsp[i]);
469 clear_rvec(mtrans[i]);
483 time[nfr] = fr.time-t0;
485 if (time[nfr] <= bfit)
489 if (time[nfr] <= efit)
499 remove_jump(fr.box, fr.natoms, xp, fr.x);
506 for (i = 0; i < fr.natoms; i++)
508 copy_rvec(fr.x[i], xp[i]);
513 gmx_rmpbc_trxfr(gpbc, &fr);
515 calc_mj(top, ePBC, fr.box, bNoJump, nmols, indexm, fr.x, mtrans[nfr], mass2, qmol);
517 for (i = 0; i < isize; i++)
520 svmul(top.atoms.atom[j].q, fr.x[j], fr.x[j]);
521 rvec_inc(mu[nfr], fr.x[j]);
524 /*if(mod(nfr,nshift)==0){*/
527 for (j = nfr; j >= 0; j--)
529 rvec_sub(mtrans[nfr], mtrans[j], tmp);
530 dsp2[nfr-j] += norm2(tmp);
548 srenew(cacf, valloc);
551 if (time[nfr] <= bvit)
555 if (time[nfr] <= evit)
560 clear_rvec(v0[nvfr]);
569 for (i = 0; i < isize; 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]);
577 fprintf(fcur, "%.3f\t%.6f\t%.6f\t%.6f\n", time[nfr], v0[nfr][XX], v0[nfr][YY], v0[nfr][ZZ]);
580 /*if(mod(nvfr,nshift)==0){*/
581 if (nvfr%nshift == 0)
583 for (j = nvfr; j >= 0; j--)
587 cacf[nvfr-j] += iprod(v0[nvfr], v0[j]);
591 djc[nvfr-j] += iprod(mu[vfr[j]], v0[nvfr]);
600 volume = det(fr.box);
603 rvec_inc(mja_tmp, mtrans[nfr]);
604 mjd += iprod(mu[nfr], mtrans[nfr]);
605 rvec_inc(mdvec, mu[nfr]);
607 mj2 += iprod(mtrans[nfr], mtrans[nfr]);
608 md2 += iprod(mu[nfr], mu[nfr]);
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);
617 while (read_next_frame(oenv, status, &fr));
619 gmx_rmpbc_done(gpbc);
624 prefactor /= 3.0*EPSILON0*volume_av*BOLTZ*temp;
627 prefactorav = E_CHARGE*E_CHARGE;
628 prefactorav /= volume_av*BOLTZMANN*temp*NANO*6.0;
630 fprintf(stderr, "Prefactor fit E-H: 1 / 6.0*V*k_B*T: %g\n", prefactorav);
632 calc_mjdsp(fmjdsp, prefactorav, dsp2, time, nfr, xshfr);
635 * Now we can average and calculate the correlation functions
643 svmul(1.0/refr, mdvec, mdvec);
644 svmul(1.0/refr, mja_tmp, mja_tmp);
646 mdav2 = norm2(mdvec);
648 mjdav = iprod(mdvec, mja_tmp);
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);
656 trust *= (real) nvfr;
657 itrust = (int) trust;
662 printf("\nCalculating M_D - current correlation integral ... \n");
663 corint = calc_cacf(mcor, prefactorav/EPSI0, djc, time, nvfr, vfr, ie, nshift);
670 printf("\nCalculating current autocorrelation ... \n");
671 sgk = calc_cacf(outf, prefactorav/PICO, cacf, time, nvfr, vfr, ie, nshift);
679 for (i = ii; i <= ie; i++)
682 xfit[i-ii] = log(time[vfr[i]]);
683 rtmp = fabs(cacf[i]);
684 yfit[i-ii] = log(rtmp);
688 lsq_y_ax_b(ie-ii, xfit, yfit, &sigma, &malt, &err, &rest);
694 malt *= prefactorav*2.0e12/sigma;
704 /* Calculation of the dielectric constant */
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");
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");
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);
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);
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);
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);
744 for (i = bi; i <= ei; i++)
746 xfit[i-bi] = time[i];
747 yfit[i-bi] = dsp2[i];
750 lsq_y_ax_b(ei-bi, xfit, yfit, &sigma, &malt, &err, &rest);
753 dk_d = calceps(prefactor, md2, 0.5*malt/prefactorav, corint, eps_rf, TRUE);
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);
765 fprintf(stderr, "Too less points for a fit.\n");
791 int gmx_current(int argc, char *argv[])
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;
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."}
827 char **grpname = NULL;
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 }
870 #define NFILE asize(fnm)
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."
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.",
892 "Option [TT]-temp[tt] sets the temperature required for the computation of the static dielectric constant.",
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).",
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."
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))
916 bACF = opt2bSet("-caf", NFILE, fnm);
917 bINT = opt2bSet("-mc", NFILE, fnm);
919 bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC, &xtop, &vtop, box, TRUE);
925 indexfn = ftp2fn_null(efNDX, NFILE, fnm);
928 get_index(&(top.atoms), indexfn, 1, &isize, &index0, grpname);
930 flags = flags | TRX_READ_X | TRX_READ_V;
932 read_first_frame(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &fr, flags);
934 snew(mass2, top.atoms.nr);
935 snew(qmol, top.atoms.nr);
937 bNEU = precalc(top, mass2, qmol);
942 for (i = 0; i < isize; i++)
944 indexm[i] = index0[i];
950 index_atom2mol(&nmols, indexm, &top.mols);
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");
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");
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");
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");
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",
990 /* System information is read and prepared, dielectric() processes the frames
991 * and calculates the requested quantities */
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);