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"
43 #include "gromacs/fileio/tpxio.h"
44 #include "gromacs/fileio/trxio.h"
50 #include "gromacs/statistics/statistics.h"
54 #define SQR(x) (pow(x, 2.0))
55 #define EPSI0 (EPSILON0*E_CHARGE*E_CHARGE*AVOGADRO/(KILO*NANO)) /* EPSILON0 in SI units */
58 static void index_atom2mol(int *n, int *index, t_block *mols)
60 int nat, i, nmol, mol, j;
68 while (index[i] > mols->index[mol])
73 gmx_fatal(FARGS, "Atom index out of range: %d", index[i]+1);
76 for (j = mols->index[mol]; j < mols->index[mol+1]; j++)
78 if (i >= nat || index[i] != j)
80 gmx_fatal(FARGS, "The index group does not consist of whole molecules");
87 fprintf(stderr, "\nSplit group of %d atoms into %d molecules\n", nat, nmol);
92 static gmx_bool precalc(t_topology top, real mass2[], real qmol[])
107 for (i = 0; i < top.mols.nr; i++)
109 k = top.mols.index[i];
110 l = top.mols.index[i+1];
114 for (j = k; j < l; j++)
116 mtot += top.atoms.atom[j].m;
117 qtot += top.atoms.atom[j].q;
120 for (j = k; j < l; j++)
122 top.atoms.atom[j].q -= top.atoms.atom[j].m*qtot/mtot;
123 mass2[j] = top.atoms.atom[j].m/mtot;
140 if (abs(qall) > 0.01)
142 printf("\n\nSystem not neutral (q=%f) will not calculate translational part of the dipole moment.\n", qall);
154 static void remove_jump(matrix box, int natoms, rvec xp[], rvec x[])
160 for (d = 0; d < DIM; d++)
162 hbox[d] = 0.5*box[d][d];
164 for (i = 0; i < natoms; i++)
166 for (m = DIM-1; m >= 0; m--)
168 while (x[i][m]-xp[i][m] <= -hbox[m])
170 for (d = 0; d <= m; d++)
172 x[i][d] += box[m][d];
175 while (x[i][m]-xp[i][m] > hbox[m])
177 for (d = 0; d <= m; d++)
179 x[i][d] -= box[m][d];
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[])
197 calc_box_center(ecenterRECT, box, center);
201 set_pbc(&pbc, ePBC, box);
207 for (i = 0; i < isize; i++)
211 k = top.mols.index[index0[i]];
212 l = top.mols.index[index0[i+1]];
213 for (j = k; j < l; j++)
215 svmul(mass2[j], fr[j], tmp);
221 svmul(qmol[k], mt1, mt1);
225 pbc_dx(&pbc, mt1, center, mt2);
226 svmul(qmol[k], mt2, mt1);
236 static real calceps(real prefactor, real md2, real mj2, real cor, real eps_rf, gmx_bool bCOR)
239 /* bCOR determines if the correlation is computed via
240 * static properties (FALSE) or the correlation integral (TRUE)
248 epsilon = md2-2.0*cor+mj2;
252 epsilon = md2+mj2+2.0*cor;
257 epsilon = 1.0+prefactor*epsilon;
262 epsilon = 2.0*eps_rf+1.0+2.0*eps_rf*prefactor*epsilon;
263 epsilon /= 2.0*eps_rf+1.0-prefactor*epsilon;
272 static real calc_cacf(FILE *fcacf, real prefactor, real cacf[], real time[], int nfr, int vfr[], int ei, int nshift)
280 real sigma_ret = 0.0;
290 rfr = (real) (nfr/nshift-i/nshift);
293 if (time[vfr[i]] <= time[vfr[ei]])
298 fprintf(fcacf, "%.3f\t%10.6g\t%10.6g\n", time[vfr[i]], cacf[i], sigma);
302 deltat = (time[vfr[i+1]]-time[vfr[i]]);
304 corint = 2.0*deltat*cacf[i]*prefactor;
305 if (i == 0 || (i+1) == nfr)
317 printf("Too less points.\n");
324 static void calc_mjdsp(FILE *fmjdsp, real prefactor, real dsp2[], real time[], int nfr, real refr[])
332 fprintf(fmjdsp, "#Prefactor fit E-H: 1 / 6.0*V*k_B*T: %g\n", prefactor);
336 for (i = 0; i < nfr; i++)
341 dsp2[i] *= prefactor/refr[i];
342 fprintf(fmjdsp, "%.3f\t%10.6g\n", time[i], dsp2[i]);
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)
360 int valloc, nalloc, nfr, nvfr, m, itrust = 0;
370 real prefactorav = 0.0;
371 real prefactor = 0.0;
373 real volume_av = 0.0;
374 real dk_s, dk_d, dk_f;
401 * Variables for the least-squares fit for Einstein-Helfand and Green-Kubo
411 gmx_rmpbc_t gpbc = NULL;
430 /* This is the main loop over frames */
445 gpbc = gmx_rmpbc_init(&top.idef, ePBC, fr.natoms);
450 refr = (real)(nfr+1);
455 srenew(time, nalloc);
457 srenew(mjdsp, nalloc);
458 srenew(dsp2, nalloc);
459 srenew(mtrans, nalloc);
460 srenew(xshfr, nalloc);
462 for (i = nfr; i < nalloc; i++)
464 clear_rvec(mjdsp[i]);
466 clear_rvec(mtrans[i]);
480 time[nfr] = fr.time-t0;
482 if (time[nfr] <= bfit)
486 if (time[nfr] <= efit)
496 remove_jump(fr.box, fr.natoms, xp, fr.x);
503 for (i = 0; i < fr.natoms; i++)
505 copy_rvec(fr.x[i], xp[i]);
510 gmx_rmpbc_trxfr(gpbc, &fr);
512 calc_mj(top, ePBC, fr.box, bNoJump, nmols, indexm, fr.x, mtrans[nfr], mass2, qmol);
514 for (i = 0; i < isize; i++)
517 svmul(top.atoms.atom[j].q, fr.x[j], fr.x[j]);
518 rvec_inc(mu[nfr], fr.x[j]);
521 /*if(mod(nfr,nshift)==0){*/
524 for (j = nfr; j >= 0; j--)
526 rvec_sub(mtrans[nfr], mtrans[j], tmp);
527 dsp2[nfr-j] += norm2(tmp);
545 srenew(cacf, valloc);
548 if (time[nfr] <= bvit)
552 if (time[nfr] <= evit)
557 clear_rvec(v0[nvfr]);
566 for (i = 0; i < isize; 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]);
574 fprintf(fcur, "%.3f\t%.6f\t%.6f\t%.6f\n", time[nfr], v0[nfr][XX], v0[nfr][YY], v0[nfr][ZZ]);
577 /*if(mod(nvfr,nshift)==0){*/
578 if (nvfr%nshift == 0)
580 for (j = nvfr; j >= 0; j--)
584 cacf[nvfr-j] += iprod(v0[nvfr], v0[j]);
588 djc[nvfr-j] += iprod(mu[vfr[j]], v0[nvfr]);
597 volume = det(fr.box);
600 rvec_inc(mja_tmp, mtrans[nfr]);
601 mjd += iprod(mu[nfr], mtrans[nfr]);
602 rvec_inc(mdvec, mu[nfr]);
604 mj2 += iprod(mtrans[nfr], mtrans[nfr]);
605 md2 += iprod(mu[nfr], mu[nfr]);
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);
614 while (read_next_frame(oenv, status, &fr));
616 gmx_rmpbc_done(gpbc);
621 prefactor /= 3.0*EPSILON0*volume_av*BOLTZ*temp;
624 prefactorav = E_CHARGE*E_CHARGE;
625 prefactorav /= volume_av*BOLTZMANN*temp*NANO*6.0;
627 fprintf(stderr, "Prefactor fit E-H: 1 / 6.0*V*k_B*T: %g\n", prefactorav);
629 calc_mjdsp(fmjdsp, prefactorav, dsp2, time, nfr, xshfr);
632 * Now we can average and calculate the correlation functions
640 svmul(1.0/refr, mdvec, mdvec);
641 svmul(1.0/refr, mja_tmp, mja_tmp);
643 mdav2 = norm2(mdvec);
645 mjdav = iprod(mdvec, mja_tmp);
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);
653 trust *= (real) nvfr;
654 itrust = (int) trust;
659 printf("\nCalculating M_D - current correlation integral ... \n");
660 corint = calc_cacf(mcor, prefactorav/EPSI0, djc, time, nvfr, vfr, ie, nshift);
667 printf("\nCalculating current autocorrelation ... \n");
668 sgk = calc_cacf(outf, prefactorav/PICO, cacf, time, nvfr, vfr, ie, nshift);
676 for (i = ii; i <= ie; i++)
679 xfit[i-ii] = log(time[vfr[i]]);
680 rtmp = fabs(cacf[i]);
681 yfit[i-ii] = log(rtmp);
685 lsq_y_ax_b(ie-ii, xfit, yfit, &sigma, &malt, &err, &rest);
691 malt *= prefactorav*2.0e12/sigma;
701 /* Calculation of the dielectric constant */
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");
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");
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);
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);
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);
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);
741 for (i = bi; i <= ei; i++)
743 xfit[i-bi] = time[i];
744 yfit[i-bi] = dsp2[i];
747 lsq_y_ax_b(ei-bi, xfit, yfit, &sigma, &malt, &err, &rest);
750 dk_d = calceps(prefactor, md2, 0.5*malt/prefactorav, corint, eps_rf, TRUE);
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);
762 fprintf(stderr, "Too less points for a fit.\n");
788 int gmx_current(int argc, char *argv[])
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;
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."}
824 char **grpname = NULL;
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 }
867 #define NFILE asize(fnm)
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."
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.",
889 "Option [TT]-temp[tt] sets the temperature required for the computation of the static dielectric constant.",
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).",
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."
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))
913 bACF = opt2bSet("-caf", NFILE, fnm);
914 bINT = opt2bSet("-mc", NFILE, fnm);
916 bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC, &xtop, &vtop, box, TRUE);
922 indexfn = ftp2fn_null(efNDX, NFILE, fnm);
925 get_index(&(top.atoms), indexfn, 1, &isize, &index0, grpname);
927 flags = flags | TRX_READ_X | TRX_READ_V;
929 read_first_frame(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &fr, flags);
931 snew(mass2, top.atoms.nr);
932 snew(qmol, top.atoms.nr);
934 bNEU = precalc(top, mass2, qmol);
939 for (i = 0; i < isize; i++)
941 indexm[i] = index0[i];
947 index_atom2mol(&nmols, indexm, &top.mols);
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");
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");
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");
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");
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",
987 /* System information is read and prepared, dielectric() processes the frames
988 * and calculates the requested quantities */
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);