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.
39 #include "gromacs/commandline/pargs.h"
41 #include "gromacs/utility/smalloc.h"
43 #include "gromacs/fileio/tpxio.h"
44 #include "gromacs/fileio/trxio.h"
50 #include "gromacs/statistics/statistics.h"
54 #include "gromacs/legacyheaders/gmx_fatal.h"
56 #define SQR(x) (pow(x, 2.0))
57 #define EPSI0 (EPSILON0*E_CHARGE*E_CHARGE*AVOGADRO/(KILO*NANO)) /* EPSILON0 in SI units */
59 static void index_atom2mol(int *n, int *index, t_block *mols)
61 int nat, i, nmol, mol, j;
69 while (index[i] > mols->index[mol])
74 gmx_fatal(FARGS, "Atom index out of range: %d", index[i]+1);
77 for (j = mols->index[mol]; j < mols->index[mol+1]; j++)
79 if (i >= nat || index[i] != j)
81 gmx_fatal(FARGS, "The index group does not consist of whole molecules");
88 fprintf(stderr, "\nSplit group of %d atoms into %d molecules\n", nat, nmol);
93 static gmx_bool precalc(t_topology top, real mass2[], real qmol[])
108 for (i = 0; i < top.mols.nr; i++)
110 k = top.mols.index[i];
111 l = top.mols.index[i+1];
115 for (j = k; j < l; j++)
117 mtot += top.atoms.atom[j].m;
118 qtot += top.atoms.atom[j].q;
121 for (j = k; j < l; j++)
123 top.atoms.atom[j].q -= top.atoms.atom[j].m*qtot/mtot;
124 mass2[j] = top.atoms.atom[j].m/mtot;
141 if (fabs(qall) > 0.01)
143 printf("\n\nSystem not neutral (q=%f) will not calculate translational part of the dipole moment.\n", qall);
155 static void remove_jump(matrix box, int natoms, rvec xp[], rvec x[])
161 for (d = 0; d < DIM; d++)
163 hbox[d] = 0.5*box[d][d];
165 for (i = 0; i < natoms; i++)
167 for (m = DIM-1; m >= 0; m--)
169 while (x[i][m]-xp[i][m] <= -hbox[m])
171 for (d = 0; d <= m; d++)
173 x[i][d] += box[m][d];
176 while (x[i][m]-xp[i][m] > hbox[m])
178 for (d = 0; d <= m; d++)
180 x[i][d] -= box[m][d];
187 static void calc_mj(t_topology top, int ePBC, matrix box, gmx_bool bNoJump, int isize, int index0[], \
188 rvec fr[], rvec mj, real mass2[], real qmol[])
198 calc_box_center(ecenterRECT, box, center);
202 set_pbc(&pbc, ePBC, box);
208 for (i = 0; i < isize; i++)
212 k = top.mols.index[index0[i]];
213 l = top.mols.index[index0[i+1]];
214 for (j = k; j < l; j++)
216 svmul(mass2[j], fr[j], tmp);
222 svmul(qmol[k], mt1, mt1);
226 pbc_dx(&pbc, mt1, center, mt2);
227 svmul(qmol[k], mt2, mt1);
237 static real calceps(real prefactor, real md2, real mj2, real cor, real eps_rf, gmx_bool bCOR)
240 /* bCOR determines if the correlation is computed via
241 * static properties (FALSE) or the correlation integral (TRUE)
249 epsilon = md2-2.0*cor+mj2;
253 epsilon = md2+mj2+2.0*cor;
258 epsilon = 1.0+prefactor*epsilon;
263 epsilon = 2.0*eps_rf+1.0+2.0*eps_rf*prefactor*epsilon;
264 epsilon /= 2.0*eps_rf+1.0-prefactor*epsilon;
273 static real calc_cacf(FILE *fcacf, real prefactor, real cacf[], real time[], int nfr, int vfr[], int ei, int nshift)
281 real sigma_ret = 0.0;
291 rfr = (real) (nfr/nshift-i/nshift);
294 if (time[vfr[i]] <= time[vfr[ei]])
299 fprintf(fcacf, "%.3f\t%10.6g\t%10.6g\n", time[vfr[i]], cacf[i], sigma);
303 deltat = (time[vfr[i+1]]-time[vfr[i]]);
305 corint = 2.0*deltat*cacf[i]*prefactor;
306 if (i == 0 || (i+1) == nfr)
318 printf("Too less points.\n");
325 static void calc_mjdsp(FILE *fmjdsp, real prefactor, real dsp2[], real time[], int nfr, real refr[])
333 fprintf(fmjdsp, "#Prefactor fit E-H: 1 / 6.0*V*k_B*T: %g\n", prefactor);
337 for (i = 0; i < nfr; i++)
342 dsp2[i] *= prefactor/refr[i];
343 fprintf(fmjdsp, "%.3f\t%10.6g\n", time[i], dsp2[i]);
352 static void dielectric(FILE *fmj, FILE *fmd, FILE *outf, FILE *fcur, FILE *mcor,
353 FILE *fmjdsp, gmx_bool bNoJump, gmx_bool bACF, gmx_bool bINT,
354 int ePBC, t_topology top, t_trxframe fr, real temp,
355 real trust, real bfit, real efit, real bvit, real evit,
356 t_trxstatus *status, int isize, int nmols, int nshift,
357 atom_id *index0, int indexm[], real mass2[],
358 real qmol[], real eps_rf, const output_env_t oenv)
361 int valloc, nalloc, nfr, nvfr, m, itrust = 0;
371 real prefactorav = 0.0;
372 real prefactor = 0.0;
374 real volume_av = 0.0;
375 real dk_s, dk_d, dk_f;
402 * Variables for the least-squares fit for Einstein-Helfand and Green-Kubo
412 gmx_rmpbc_t gpbc = NULL;
431 /* This is the main loop over frames */
446 gpbc = gmx_rmpbc_init(&top.idef, ePBC, fr.natoms);
451 refr = (real)(nfr+1);
456 srenew(time, nalloc);
458 srenew(mjdsp, nalloc);
459 srenew(dsp2, nalloc);
460 srenew(mtrans, nalloc);
461 srenew(xshfr, nalloc);
463 for (i = nfr; i < nalloc; i++)
465 clear_rvec(mjdsp[i]);
467 clear_rvec(mtrans[i]);
481 time[nfr] = fr.time-t0;
483 if (time[nfr] <= bfit)
487 if (time[nfr] <= efit)
497 remove_jump(fr.box, fr.natoms, xp, fr.x);
504 for (i = 0; i < fr.natoms; i++)
506 copy_rvec(fr.x[i], xp[i]);
511 gmx_rmpbc_trxfr(gpbc, &fr);
513 calc_mj(top, ePBC, fr.box, bNoJump, nmols, indexm, fr.x, mtrans[nfr], mass2, qmol);
515 for (i = 0; i < isize; i++)
518 svmul(top.atoms.atom[j].q, fr.x[j], fr.x[j]);
519 rvec_inc(mu[nfr], fr.x[j]);
522 /*if(mod(nfr,nshift)==0){*/
525 for (j = nfr; j >= 0; j--)
527 rvec_sub(mtrans[nfr], mtrans[j], tmp);
528 dsp2[nfr-j] += norm2(tmp);
546 srenew(cacf, valloc);
549 if (time[nfr] <= bvit)
553 if (time[nfr] <= evit)
558 clear_rvec(v0[nvfr]);
567 for (i = 0; i < isize; i++)
570 svmul(mass2[j], fr.v[j], fr.v[j]);
571 svmul(qmol[j], fr.v[j], fr.v[j]);
572 rvec_inc(v0[nvfr], fr.v[j]);
575 fprintf(fcur, "%.3f\t%.6f\t%.6f\t%.6f\n", time[nfr], v0[nfr][XX], v0[nfr][YY], v0[nfr][ZZ]);
578 /*if(mod(nvfr,nshift)==0){*/
579 if (nvfr%nshift == 0)
581 for (j = nvfr; j >= 0; j--)
585 cacf[nvfr-j] += iprod(v0[nvfr], v0[j]);
589 djc[nvfr-j] += iprod(mu[vfr[j]], v0[nvfr]);
598 volume = det(fr.box);
601 rvec_inc(mja_tmp, mtrans[nfr]);
602 mjd += iprod(mu[nfr], mtrans[nfr]);
603 rvec_inc(mdvec, mu[nfr]);
605 mj2 += iprod(mtrans[nfr], mtrans[nfr]);
606 md2 += iprod(mu[nfr], mu[nfr]);
608 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);
609 fprintf(fmd, "%.3f\t%8.5f\t%8.5f\t%8.5f\t%8.5f\t%8.5f\n", \
610 time[nfr], mu[nfr][XX], mu[nfr][YY], mu[nfr][ZZ], md2/refr, norm(mdvec)/refr);
615 while (read_next_frame(oenv, status, &fr));
617 gmx_rmpbc_done(gpbc);
622 prefactor /= 3.0*EPSILON0*volume_av*BOLTZ*temp;
625 prefactorav = E_CHARGE*E_CHARGE;
626 prefactorav /= volume_av*BOLTZMANN*temp*NANO*6.0;
628 fprintf(stderr, "Prefactor fit E-H: 1 / 6.0*V*k_B*T: %g\n", prefactorav);
630 calc_mjdsp(fmjdsp, prefactorav, dsp2, time, nfr, xshfr);
633 * Now we can average and calculate the correlation functions
641 svmul(1.0/refr, mdvec, mdvec);
642 svmul(1.0/refr, mja_tmp, mja_tmp);
644 mdav2 = norm2(mdvec);
646 mjdav = iprod(mdvec, mja_tmp);
649 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);
650 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 trust *= (real) nvfr;
655 itrust = (int) trust;
660 printf("\nCalculating M_D - current correlation integral ... \n");
661 corint = calc_cacf(mcor, prefactorav/EPSI0, djc, time, nvfr, vfr, ie, nshift);
668 printf("\nCalculating current autocorrelation ... \n");
669 sgk = calc_cacf(outf, prefactorav/PICO, cacf, time, nvfr, vfr, ie, nshift);
677 for (i = ii; i <= ie; i++)
680 xfit[i-ii] = log(time[vfr[i]]);
681 rtmp = fabs(cacf[i]);
682 yfit[i-ii] = log(rtmp);
686 lsq_y_ax_b(ie-ii, xfit, yfit, &sigma, &malt, &err, &rest);
692 malt *= prefactorav*2.0e12/sigma;
702 /* Calculation of the dielectric constant */
704 fprintf(stderr, "\n********************************************\n");
705 dk_s = calceps(prefactor, md2, mj2, mjd, eps_rf, FALSE);
706 fprintf(stderr, "\nAbsolute values:\n epsilon=%f\n", dk_s);
707 fprintf(stderr, " <M_D^2> , <M_J^2>, <(M_J*M_D)^2>: (%f, %f, %f)\n\n", md2, mj2, mjd);
708 fprintf(stderr, "********************************************\n");
711 dk_f = calceps(prefactor, md2-mdav2, mj2-mj, mjd-mjdav, eps_rf, FALSE);
712 fprintf(stderr, "\n\nFluctuations:\n epsilon=%f\n\n", dk_f);
713 fprintf(stderr, "\n deltaM_D , deltaM_J, deltaM_JD: (%f, %f, %f)\n", md2-mdav2, mj2-mj, mjd-mjdav);
714 fprintf(stderr, "\n********************************************\n");
717 dk_d = calceps(prefactor, md2-mdav2, mj2-mj, corint, eps_rf, TRUE);
718 fprintf(stderr, "\nStatic dielectric constant using integral and fluctuations: %f\n", dk_d);
719 fprintf(stderr, "\n < M_JM_D > via integral: %.3f\n", -1.0*corint);
722 fprintf(stderr, "\n***************************************************");
723 fprintf(stderr, "\n\nAverage volume V=%f nm^3 at T=%f K\n", volume_av, temp);
724 fprintf(stderr, "and corresponding refactor 1.0 / 3.0*V*k_B*T*EPSILON_0: %f \n", prefactor);
730 fprintf(stderr, "Integral and integrated fit to the current acf yields at t=%f:\n", time[vfr[ii]]);
731 fprintf(stderr, "sigma=%8.3f (pure integral: %.3f)\n", sgk-malt*pow(time[vfr[ii]], sigma), sgk);
736 fprintf(stderr, "\nStart fit at %f ps (%f).\n", time[bi], bfit);
737 fprintf(stderr, "End fit at %f ps (%f).\n\n", time[ei], efit);
742 for (i = bi; i <= ei; i++)
744 xfit[i-bi] = time[i];
745 yfit[i-bi] = dsp2[i];
748 lsq_y_ax_b(ei-bi, xfit, yfit, &sigma, &malt, &err, &rest);
751 dk_d = calceps(prefactor, md2, 0.5*malt/prefactorav, corint, eps_rf, TRUE);
753 fprintf(stderr, "Einstein-Helfand fit to the MSD of the translational dipole moment yields:\n\n");
754 fprintf(stderr, "sigma=%.4f\n", sigma);
755 fprintf(stderr, "translational fraction of M^2: %.4f\n", 0.5*malt/prefactorav);
756 fprintf(stderr, "Dielectric constant using EH: %.4f\n", dk_d);
763 fprintf(stderr, "Too less points for a fit.\n");
789 int gmx_current(int argc, char *argv[])
792 static int nshift = 1000;
793 static real temp = 300.0;
794 static real trust = 0.25;
795 static real eps_rf = 0.0;
796 static gmx_bool bNoJump = TRUE;
797 static real bfit = 100.0;
798 static real bvit = 0.5;
799 static real efit = 400.0;
800 static real evit = 5.0;
802 { "-sh", FALSE, etINT, {&nshift},
803 "Shift of the frames for averaging the correlation functions and the mean-square displacement."},
804 { "-nojump", FALSE, etBOOL, {&bNoJump},
805 "Removes jumps of atoms across the box."},
806 { "-eps", FALSE, etREAL, {&eps_rf},
807 "Dielectric constant of the surrounding medium. The value zero corresponds to infinity (tin-foil boundary conditions)."},
808 { "-bfit", FALSE, etREAL, {&bfit},
809 "Begin of the fit of the straight line to the MSD of the translational fraction of the dipole moment."},
810 { "-efit", FALSE, etREAL, {&efit},
811 "End of the fit of the straight line to the MSD of the translational fraction of the dipole moment."},
812 { "-bvit", FALSE, etREAL, {&bvit},
813 "Begin of the fit of the current autocorrelation function to a*t^b."},
814 { "-evit", FALSE, etREAL, {&evit},
815 "End of the fit of the current autocorrelation function to a*t^b."},
816 { "-tr", FALSE, etREAL, {&trust},
817 "Fraction of the trajectory taken into account for the integral."},
818 { "-temp", FALSE, etREAL, {&temp},
819 "Temperature for calculating epsilon."}
825 char **grpname = NULL;
857 { efTPS, NULL, NULL, ffREAD }, /* this is for the topology */
858 { efNDX, NULL, NULL, ffOPTRD },
859 { efTRX, "-f", NULL, ffREAD }, /* and this for the trajectory */
860 { efXVG, "-o", "current.xvg", ffWRITE },
861 { efXVG, "-caf", "caf.xvg", ffOPTWR },
862 { efXVG, "-dsp", "dsp.xvg", ffWRITE },
863 { efXVG, "-md", "md.xvg", ffWRITE },
864 { efXVG, "-mj", "mj.xvg", ffWRITE},
865 { efXVG, "-mc", "mc.xvg", ffOPTWR }
868 #define NFILE asize(fnm)
871 const char *desc[] = {
872 "[THISMODULE] is a tool for calculating the current autocorrelation function, the correlation",
873 "of the rotational and translational dipole moment of the system, and the resulting static",
874 "dielectric constant. To obtain a reasonable result, the index group has to be neutral.",
875 "Furthermore, the routine is capable of extracting the static conductivity from the current ",
876 "autocorrelation function, if velocities are given. Additionally, an Einstein-Helfand fit ",
877 "can be used to obtain the static conductivity."
879 "The flag [TT]-caf[tt] is for the output of the current autocorrelation function and [TT]-mc[tt] writes the",
880 "correlation of the rotational and translational part of the dipole moment in the corresponding",
881 "file. However, this option is only available for trajectories containing velocities.",
882 "Options [TT]-sh[tt] and [TT]-tr[tt] are responsible for the averaging and integration of the",
883 "autocorrelation functions. Since averaging proceeds by shifting the starting point",
884 "through the trajectory, the shift can be modified with [TT]-sh[tt] to enable the choice of uncorrelated",
885 "starting points. Towards the end, statistical inaccuracy grows and integrating the",
886 "correlation function only yields reliable values until a certain point, depending on",
887 "the number of frames. The option [TT]-tr[tt] controls the region of the integral taken into account",
888 "for calculating the static dielectric constant.",
890 "Option [TT]-temp[tt] sets the temperature required for the computation of the static dielectric constant.",
892 "Option [TT]-eps[tt] controls the dielectric constant of the surrounding medium for simulations using",
893 "a Reaction Field or dipole corrections of the Ewald summation ([TT]-eps[tt]=0 corresponds to",
894 "tin-foil boundary conditions).",
896 "[TT]-[no]nojump[tt] unfolds the coordinates to allow free diffusion. This is required to get a continuous",
897 "translational dipole moment, required for the Einstein-Helfand fit. The results from the fit allow",
898 "the determination of the dielectric constant for system of charged molecules. However, it is also possible to extract",
899 "the dielectric constant from the fluctuations of the total dipole moment in folded coordinates. But this",
900 "option has to be used with care, since only very short time spans fulfill the approximation that the density",
901 "of the molecules is approximately constant and the averages are already converged. To be on the safe side,",
902 "the dielectric constant should be calculated with the help of the Einstein-Helfand method for",
903 "the translational part of the dielectric constant."
907 /* At first the arguments will be parsed and the system information processed */
908 if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW,
909 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
914 bACF = opt2bSet("-caf", NFILE, fnm);
915 bINT = opt2bSet("-mc", NFILE, fnm);
917 bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC, &xtop, &vtop, box, TRUE);
923 indexfn = ftp2fn_null(efNDX, NFILE, fnm);
926 get_index(&(top.atoms), indexfn, 1, &isize, &index0, grpname);
928 flags = flags | TRX_READ_X | TRX_READ_V;
930 read_first_frame(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &fr, flags);
932 snew(mass2, top.atoms.nr);
933 snew(qmol, top.atoms.nr);
935 bNEU = precalc(top, mass2, qmol);
940 for (i = 0; i < isize; i++)
942 indexm[i] = index0[i];
948 index_atom2mol(&nmols, indexm, &top.mols);
954 outf = xvgropen(opt2fn("-caf", NFILE, fnm),
955 "Current autocorrelation function", output_env_get_xvgr_tlabel(oenv),
956 "ACF (e nm/ps)\\S2", oenv);
957 fprintf(outf, "# time\t acf\t average \t std.dev\n");
959 fcur = xvgropen(opt2fn("-o", NFILE, fnm),
960 "Current", output_env_get_xvgr_tlabel(oenv), "J(t) (e nm/ps)", oenv);
961 fprintf(fcur, "# time\t Jx\t Jy \t J_z \n");
964 mcor = xvgropen(opt2fn("-mc", NFILE, fnm),
965 "M\\sD\\N - current autocorrelation function",
966 output_env_get_xvgr_tlabel(oenv),
967 "< M\\sD\\N (0)\\c7\\CJ(t) > (e nm/ps)\\S2", oenv);
968 fprintf(mcor, "# time\t M_D(0) J(t) acf \t Integral acf\n");
972 fmj = xvgropen(opt2fn("-mj", NFILE, fnm),
973 "Averaged translational part of M", output_env_get_xvgr_tlabel(oenv),
974 "< M\\sJ\\N > (enm)", oenv);
975 fprintf(fmj, "# time\t x\t y \t z \t average of M_J^2 \t std.dev\n");
976 fmd = xvgropen(opt2fn("-md", NFILE, fnm),
977 "Averaged rotational part of M", output_env_get_xvgr_tlabel(oenv),
978 "< M\\sD\\N > (enm)", oenv);
979 fprintf(fmd, "# time\t x\t y \t z \t average of M_D^2 \t std.dev\n");
981 fmjdsp = xvgropen(opt2fn("-dsp", NFILE, fnm),
982 "MSD of the squared translational dipole moment M",
983 output_env_get_xvgr_tlabel(oenv),
984 "<|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 /* System information is read and prepared, dielectric() processes the frames
989 * and calculates the requested quantities */
991 dielectric(fmj, fmd, outf, fcur, mcor, fmjdsp, bNoJump, bACF, bINT, ePBC, top, fr,
992 temp, trust, bfit, efit, bvit, evit, status, isize, nmols, nshift,
993 index0, indexm, mass2, qmol, eps_rf, oenv);