3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Copyright (c) 1991-2001
11 * BIOSON Research Institute, Dept. of Biophysical Chemistry
12 * University of Groningen, The Netherlands
14 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * Do check out http://www.gromacs.org , or mail us at gromacs@gromacs.org .
33 * Gyas ROwers Mature At Cryogenic Speed
35 * finished FD 09/07/08
48 #include "gromacs/fileio/tpxio.h"
49 #include "gromacs/fileio/trxio.h"
55 #include "gmx_statistics.h"
59 #define SQR(x) (pow(x, 2.0))
60 #define EPSI0 (EPSILON0*E_CHARGE*E_CHARGE*AVOGADRO/(KILO*NANO)) /* EPSILON0 in SI units */
63 static void index_atom2mol(int *n, int *index, t_block *mols)
65 int nat, i, nmol, mol, j;
73 while (index[i] > mols->index[mol])
78 gmx_fatal(FARGS, "Atom index out of range: %d", index[i]+1);
81 for (j = mols->index[mol]; j < mols->index[mol+1]; j++)
83 if (i >= nat || index[i] != j)
85 gmx_fatal(FARGS, "The index group does not consist of whole molecules");
92 fprintf(stderr, "\nSplit group of %d atoms into %d molecules\n", nat, nmol);
97 static gmx_bool precalc(t_topology top, real mass2[], real qmol[])
112 for (i = 0; i < top.mols.nr; i++)
114 k = top.mols.index[i];
115 l = top.mols.index[i+1];
119 for (j = k; j < l; j++)
121 mtot += top.atoms.atom[j].m;
122 qtot += top.atoms.atom[j].q;
125 for (j = k; j < l; j++)
127 top.atoms.atom[j].q -= top.atoms.atom[j].m*qtot/mtot;
128 mass2[j] = top.atoms.atom[j].m/mtot;
145 if (abs(qall) > 0.01)
147 printf("\n\nSystem not neutral (q=%f) will not calculate translational part of the dipole moment.\n", qall);
159 static void remove_jump(matrix box, int natoms, rvec xp[], rvec x[])
165 for (d = 0; d < DIM; d++)
167 hbox[d] = 0.5*box[d][d];
169 for (i = 0; i < natoms; i++)
171 for (m = DIM-1; m >= 0; m--)
173 while (x[i][m]-xp[i][m] <= -hbox[m])
175 for (d = 0; d <= m; d++)
177 x[i][d] += box[m][d];
180 while (x[i][m]-xp[i][m] > hbox[m])
182 for (d = 0; d <= m; d++)
184 x[i][d] -= box[m][d];
191 static void calc_mj(t_topology top, int ePBC, matrix box, gmx_bool bNoJump, int isize, int index0[], \
192 rvec fr[], rvec mj, real mass2[], real qmol[])
202 calc_box_center(ecenterRECT, box, center);
206 set_pbc(&pbc, ePBC, box);
212 for (i = 0; i < isize; i++)
216 k = top.mols.index[index0[i]];
217 l = top.mols.index[index0[i+1]];
218 for (j = k; j < l; j++)
220 svmul(mass2[j], fr[j], tmp);
226 svmul(qmol[k], mt1, mt1);
230 pbc_dx(&pbc, mt1, center, mt2);
231 svmul(qmol[k], mt2, mt1);
241 static real calceps(real prefactor, real md2, real mj2, real cor, real eps_rf, gmx_bool bCOR)
244 /* bCOR determines if the correlation is computed via
245 * static properties (FALSE) or the correlation integral (TRUE)
253 epsilon = md2-2.0*cor+mj2;
257 epsilon = md2+mj2+2.0*cor;
262 epsilon = 1.0+prefactor*epsilon;
267 epsilon = 2.0*eps_rf+1.0+2.0*eps_rf*prefactor*epsilon;
268 epsilon /= 2.0*eps_rf+1.0-prefactor*epsilon;
277 static real calc_cacf(FILE *fcacf, real prefactor, real cacf[], real time[], int nfr, int vfr[], int ei, int nshift)
285 real sigma_ret = 0.0;
295 rfr = (real) (nfr/nshift-i/nshift);
298 if (time[vfr[i]] <= time[vfr[ei]])
303 fprintf(fcacf, "%.3f\t%10.6g\t%10.6g\n", time[vfr[i]], cacf[i], sigma);
307 deltat = (time[vfr[i+1]]-time[vfr[i]]);
309 corint = 2.0*deltat*cacf[i]*prefactor;
310 if (i == 0 || (i+1) == nfr)
322 printf("Too less points.\n");
329 static void calc_mjdsp(FILE *fmjdsp, real prefactor, real dsp2[], real time[], int nfr, real refr[])
337 fprintf(fmjdsp, "#Prefactor fit E-H: 1 / 6.0*V*k_B*T: %g\n", prefactor);
341 for (i = 0; i < nfr; i++)
346 dsp2[i] *= prefactor/refr[i];
347 fprintf(fmjdsp, "%.3f\t%10.6g\n", time[i], dsp2[i]);
356 static void dielectric(FILE *fmj, FILE *fmd, FILE *outf, FILE *fcur, FILE *mcor,
357 FILE *fmjdsp, gmx_bool bNoJump, gmx_bool bACF, gmx_bool bINT,
358 int ePBC, t_topology top, t_trxframe fr, real temp,
359 real trust, real bfit, real efit, real bvit, real evit,
360 t_trxstatus *status, int isize, int nmols, int nshift,
361 atom_id *index0, int indexm[], real mass2[],
362 real qmol[], real eps_rf, const output_env_t oenv)
365 int valloc, nalloc, nfr, nvfr, m, itrust = 0;
375 real prefactorav = 0.0;
376 real prefactor = 0.0;
378 real volume_av = 0.0;
379 real dk_s, dk_d, dk_f;
406 * Variables for the least-squares fit for Einstein-Helfand and Green-Kubo
416 gmx_rmpbc_t gpbc = NULL;
435 /* This is the main loop over frames */
450 gpbc = gmx_rmpbc_init(&top.idef, ePBC, fr.natoms);
455 refr = (real)(nfr+1);
460 srenew(time, nalloc);
462 srenew(mjdsp, nalloc);
463 srenew(dsp2, nalloc);
464 srenew(mtrans, nalloc);
465 srenew(xshfr, nalloc);
467 for (i = nfr; i < nalloc; i++)
469 clear_rvec(mjdsp[i]);
471 clear_rvec(mtrans[i]);
485 time[nfr] = fr.time-t0;
487 if (time[nfr] <= bfit)
491 if (time[nfr] <= efit)
501 remove_jump(fr.box, fr.natoms, xp, fr.x);
508 for (i = 0; i < fr.natoms; i++)
510 copy_rvec(fr.x[i], xp[i]);
515 gmx_rmpbc_trxfr(gpbc, &fr);
517 calc_mj(top, ePBC, fr.box, bNoJump, nmols, indexm, fr.x, mtrans[nfr], mass2, qmol);
519 for (i = 0; i < isize; i++)
522 svmul(top.atoms.atom[j].q, fr.x[j], fr.x[j]);
523 rvec_inc(mu[nfr], fr.x[j]);
526 /*if(mod(nfr,nshift)==0){*/
529 for (j = nfr; j >= 0; j--)
531 rvec_sub(mtrans[nfr], mtrans[j], tmp);
532 dsp2[nfr-j] += norm2(tmp);
550 srenew(cacf, valloc);
553 if (time[nfr] <= bvit)
557 if (time[nfr] <= evit)
562 clear_rvec(v0[nvfr]);
571 for (i = 0; i < isize; i++)
574 svmul(mass2[j], fr.v[j], fr.v[j]);
575 svmul(qmol[j], fr.v[j], fr.v[j]);
576 rvec_inc(v0[nvfr], fr.v[j]);
579 fprintf(fcur, "%.3f\t%.6f\t%.6f\t%.6f\n", time[nfr], v0[nfr][XX], v0[nfr][YY], v0[nfr][ZZ]);
582 /*if(mod(nvfr,nshift)==0){*/
583 if (nvfr%nshift == 0)
585 for (j = nvfr; j >= 0; j--)
589 cacf[nvfr-j] += iprod(v0[nvfr], v0[j]);
593 djc[nvfr-j] += iprod(mu[vfr[j]], v0[nvfr]);
602 volume = det(fr.box);
605 rvec_inc(mja_tmp, mtrans[nfr]);
606 mjd += iprod(mu[nfr], mtrans[nfr]);
607 rvec_inc(mdvec, mu[nfr]);
609 mj2 += iprod(mtrans[nfr], mtrans[nfr]);
610 md2 += iprod(mu[nfr], mu[nfr]);
612 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);
613 fprintf(fmd, "%.3f\t%8.5f\t%8.5f\t%8.5f\t%8.5f\t%8.5f\n", \
614 time[nfr], mu[nfr][XX], mu[nfr][YY], mu[nfr][ZZ], md2/refr, norm(mdvec)/refr);
619 while (read_next_frame(oenv, status, &fr));
621 gmx_rmpbc_done(gpbc);
626 prefactor /= 3.0*EPSILON0*volume_av*BOLTZ*temp;
629 prefactorav = E_CHARGE*E_CHARGE;
630 prefactorav /= volume_av*BOLTZMANN*temp*NANO*6.0;
632 fprintf(stderr, "Prefactor fit E-H: 1 / 6.0*V*k_B*T: %g\n", prefactorav);
634 calc_mjdsp(fmjdsp, prefactorav, dsp2, time, nfr, xshfr);
637 * Now we can average and calculate the correlation functions
645 svmul(1.0/refr, mdvec, mdvec);
646 svmul(1.0/refr, mja_tmp, mja_tmp);
648 mdav2 = norm2(mdvec);
650 mjdav = iprod(mdvec, mja_tmp);
653 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);
654 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);
658 trust *= (real) nvfr;
659 itrust = (int) trust;
664 printf("\nCalculating M_D - current correlation integral ... \n");
665 corint = calc_cacf(mcor, prefactorav/EPSI0, djc, time, nvfr, vfr, ie, nshift);
672 printf("\nCalculating current autocorrelation ... \n");
673 sgk = calc_cacf(outf, prefactorav/PICO, cacf, time, nvfr, vfr, ie, nshift);
681 for (i = ii; i <= ie; i++)
684 xfit[i-ii] = log(time[vfr[i]]);
685 rtmp = fabs(cacf[i]);
686 yfit[i-ii] = log(rtmp);
690 lsq_y_ax_b(ie-ii, xfit, yfit, &sigma, &malt, &err, &rest);
696 malt *= prefactorav*2.0e12/sigma;
706 /* Calculation of the dielectric constant */
708 fprintf(stderr, "\n********************************************\n");
709 dk_s = calceps(prefactor, md2, mj2, mjd, eps_rf, FALSE);
710 fprintf(stderr, "\nAbsolute values:\n epsilon=%f\n", dk_s);
711 fprintf(stderr, " <M_D^2> , <M_J^2>, <(M_J*M_D)^2>: (%f, %f, %f)\n\n", md2, mj2, mjd);
712 fprintf(stderr, "********************************************\n");
715 dk_f = calceps(prefactor, md2-mdav2, mj2-mj, mjd-mjdav, eps_rf, FALSE);
716 fprintf(stderr, "\n\nFluctuations:\n epsilon=%f\n\n", dk_f);
717 fprintf(stderr, "\n deltaM_D , deltaM_J, deltaM_JD: (%f, %f, %f)\n", md2-mdav2, mj2-mj, mjd-mjdav);
718 fprintf(stderr, "\n********************************************\n");
721 dk_d = calceps(prefactor, md2-mdav2, mj2-mj, corint, eps_rf, TRUE);
722 fprintf(stderr, "\nStatic dielectric constant using integral and fluctuations: %f\n", dk_d);
723 fprintf(stderr, "\n < M_JM_D > via integral: %.3f\n", -1.0*corint);
726 fprintf(stderr, "\n***************************************************");
727 fprintf(stderr, "\n\nAverage volume V=%f nm^3 at T=%f K\n", volume_av, temp);
728 fprintf(stderr, "and corresponding refactor 1.0 / 3.0*V*k_B*T*EPSILON_0: %f \n", prefactor);
734 fprintf(stderr, "Integral and integrated fit to the current acf yields at t=%f:\n", time[vfr[ii]]);
735 fprintf(stderr, "sigma=%8.3f (pure integral: %.3f)\n", sgk-malt*pow(time[vfr[ii]], sigma), sgk);
740 fprintf(stderr, "\nStart fit at %f ps (%f).\n", time[bi], bfit);
741 fprintf(stderr, "End fit at %f ps (%f).\n\n", time[ei], efit);
746 for (i = bi; i <= ei; i++)
748 xfit[i-bi] = time[i];
749 yfit[i-bi] = dsp2[i];
752 lsq_y_ax_b(ei-bi, xfit, yfit, &sigma, &malt, &err, &rest);
755 dk_d = calceps(prefactor, md2, 0.5*malt/prefactorav, corint, eps_rf, TRUE);
757 fprintf(stderr, "Einstein-Helfand fit to the MSD of the translational dipole moment yields:\n\n");
758 fprintf(stderr, "sigma=%.4f\n", sigma);
759 fprintf(stderr, "translational fraction of M^2: %.4f\n", 0.5*malt/prefactorav);
760 fprintf(stderr, "Dielectric constant using EH: %.4f\n", dk_d);
767 fprintf(stderr, "Too less points for a fit.\n");
793 int gmx_current(int argc, char *argv[])
796 static int nshift = 1000;
797 static real temp = 300.0;
798 static real trust = 0.25;
799 static real eps_rf = 0.0;
800 static gmx_bool bNoJump = TRUE;
801 static real bfit = 100.0;
802 static real bvit = 0.5;
803 static real efit = 400.0;
804 static real evit = 5.0;
806 { "-sh", FALSE, etINT, {&nshift},
807 "Shift of the frames for averaging the correlation functions and the mean-square displacement."},
808 { "-nojump", FALSE, etBOOL, {&bNoJump},
809 "Removes jumps of atoms across the box."},
810 { "-eps", FALSE, etREAL, {&eps_rf},
811 "Dielectric constant of the surrounding medium. The value zero corresponds to infinity (tin-foil boundary conditions)."},
812 { "-bfit", FALSE, etREAL, {&bfit},
813 "Begin of the fit of the straight line to the MSD of the translational fraction of the dipole moment."},
814 { "-efit", FALSE, etREAL, {&efit},
815 "End of the fit of the straight line to the MSD of the translational fraction of the dipole moment."},
816 { "-bvit", FALSE, etREAL, {&bvit},
817 "Begin of the fit of the current autocorrelation function to a*t^b."},
818 { "-evit", FALSE, etREAL, {&evit},
819 "End of the fit of the current autocorrelation function to a*t^b."},
820 { "-tr", FALSE, etREAL, {&trust},
821 "Fraction of the trajectory taken into account for the integral."},
822 { "-temp", FALSE, etREAL, {&temp},
823 "Temperature for calculating epsilon."}
829 char **grpname = NULL;
861 { efTPS, NULL, NULL, ffREAD }, /* this is for the topology */
862 { efNDX, NULL, NULL, ffOPTRD },
863 { efTRX, "-f", NULL, ffREAD }, /* and this for the trajectory */
864 { efXVG, "-o", "current.xvg", ffWRITE },
865 { efXVG, "-caf", "caf.xvg", ffOPTWR },
866 { efXVG, "-dsp", "dsp.xvg", ffWRITE },
867 { efXVG, "-md", "md.xvg", ffWRITE },
868 { efXVG, "-mj", "mj.xvg", ffWRITE},
869 { efXVG, "-mc", "mc.xvg", ffOPTWR }
872 #define NFILE asize(fnm)
875 const char *desc[] = {
876 "[THISMODULE] is a tool for calculating the current autocorrelation function, the correlation",
877 "of the rotational and translational dipole moment of the system, and the resulting static",
878 "dielectric constant. To obtain a reasonable result, the index group has to be neutral.",
879 "Furthermore, the routine is capable of extracting the static conductivity from the current ",
880 "autocorrelation function, if velocities are given. Additionally, an Einstein-Helfand fit ",
881 "can be used to obtain the static conductivity."
883 "The flag [TT]-caf[tt] is for the output of the current autocorrelation function and [TT]-mc[tt] writes the",
884 "correlation of the rotational and translational part of the dipole moment in the corresponding",
885 "file. However, this option is only available for trajectories containing velocities.",
886 "Options [TT]-sh[tt] and [TT]-tr[tt] are responsible for the averaging and integration of the",
887 "autocorrelation functions. Since averaging proceeds by shifting the starting point",
888 "through the trajectory, the shift can be modified with [TT]-sh[tt] to enable the choice of uncorrelated",
889 "starting points. Towards the end, statistical inaccuracy grows and integrating the",
890 "correlation function only yields reliable values until a certain point, depending on",
891 "the number of frames. The option [TT]-tr[tt] controls the region of the integral taken into account",
892 "for calculating the static dielectric constant.",
894 "Option [TT]-temp[tt] sets the temperature required for the computation of the static dielectric constant.",
896 "Option [TT]-eps[tt] controls the dielectric constant of the surrounding medium for simulations using",
897 "a Reaction Field or dipole corrections of the Ewald summation ([TT]-eps[tt]=0 corresponds to",
898 "tin-foil boundary conditions).",
900 "[TT]-[no]nojump[tt] unfolds the coordinates to allow free diffusion. This is required to get a continuous",
901 "translational dipole moment, required for the Einstein-Helfand fit. The results from the fit allow",
902 "the determination of the dielectric constant for system of charged molecules. However, it is also possible to extract",
903 "the dielectric constant from the fluctuations of the total dipole moment in folded coordinates. But this",
904 "option has to be used with care, since only very short time spans fulfill the approximation that the density",
905 "of the molecules is approximately constant and the averages are already converged. To be on the safe side,",
906 "the dielectric constant should be calculated with the help of the Einstein-Helfand method for",
907 "the translational part of the dielectric constant."
911 /* At first the arguments will be parsed and the system information processed */
912 if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW,
913 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
918 bACF = opt2bSet("-caf", NFILE, fnm);
919 bINT = opt2bSet("-mc", NFILE, fnm);
921 bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC, &xtop, &vtop, box, TRUE);
927 indexfn = ftp2fn_null(efNDX, NFILE, fnm);
930 get_index(&(top.atoms), indexfn, 1, &isize, &index0, grpname);
932 flags = flags | TRX_READ_X | TRX_READ_V;
934 read_first_frame(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &fr, flags);
936 snew(mass2, top.atoms.nr);
937 snew(qmol, top.atoms.nr);
939 bNEU = precalc(top, mass2, qmol);
944 for (i = 0; i < isize; i++)
946 indexm[i] = index0[i];
952 index_atom2mol(&nmols, indexm, &top.mols);
958 outf = xvgropen(opt2fn("-caf", NFILE, fnm),
959 "Current autocorrelation function", output_env_get_xvgr_tlabel(oenv),
960 "ACF (e nm/ps)\\S2", oenv);
961 fprintf(outf, "# time\t acf\t average \t std.dev\n");
963 fcur = xvgropen(opt2fn("-o", NFILE, fnm),
964 "Current", output_env_get_xvgr_tlabel(oenv), "J(t) (e nm/ps)", oenv);
965 fprintf(fcur, "# time\t Jx\t Jy \t J_z \n");
968 mcor = xvgropen(opt2fn("-mc", NFILE, fnm),
969 "M\\sD\\N - current autocorrelation function",
970 output_env_get_xvgr_tlabel(oenv),
971 "< M\\sD\\N (0)\\c7\\CJ(t) > (e nm/ps)\\S2", oenv);
972 fprintf(mcor, "# time\t M_D(0) J(t) acf \t Integral acf\n");
976 fmj = xvgropen(opt2fn("-mj", NFILE, fnm),
977 "Averaged translational part of M", output_env_get_xvgr_tlabel(oenv),
978 "< M\\sJ\\N > (enm)", oenv);
979 fprintf(fmj, "# time\t x\t y \t z \t average of M_J^2 \t std.dev\n");
980 fmd = xvgropen(opt2fn("-md", NFILE, fnm),
981 "Averaged rotational part of M", output_env_get_xvgr_tlabel(oenv),
982 "< M\\sD\\N > (enm)", oenv);
983 fprintf(fmd, "# time\t x\t y \t z \t average of M_D^2 \t std.dev\n");
985 fmjdsp = xvgropen(opt2fn("-dsp", NFILE, fnm),
986 "MSD of the squared translational dipole moment M",
987 output_env_get_xvgr_tlabel(oenv),
988 "<|M\\sJ\\N(t)-M\\sJ\\N(0)|\\S2\\N > / 6.0*V*k\\sB\\N*T / Sm\\S-1\\Nps\\S-1\\N",
992 /* System information is read and prepared, dielectric() processes the frames
993 * and calculates the requested quantities */
995 dielectric(fmj, fmd, outf, fcur, mcor, fmjdsp, bNoJump, bACF, bINT, ePBC, top, fr,
996 temp, trust, bfit, efit, bvit, evit, status, isize, nmols, nshift,
997 index0, indexm, mass2, qmol, eps_rf, oenv);