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.
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"
55 #include "gromacs/legacyheaders/macros.h"
57 #include "gromacs/utility/fatalerror.h"
59 #define SQR(x) (pow(x, 2.0))
60 #define EPSI0 (EPSILON0*E_CHARGE*E_CHARGE*AVOGADRO/(KILO*NANO)) /* EPSILON0 in SI units */
62 static void index_atom2mol(int *n, int *index, t_block *mols)
64 int nat, i, nmol, mol, j;
72 while (index[i] > mols->index[mol])
77 gmx_fatal(FARGS, "Atom index out of range: %d", index[i]+1);
80 for (j = mols->index[mol]; j < mols->index[mol+1]; j++)
82 if (i >= nat || index[i] != j)
84 gmx_fatal(FARGS, "The index group does not consist of whole molecules");
91 fprintf(stderr, "\nSplit group of %d atoms into %d molecules\n", nat, nmol);
96 static gmx_bool precalc(t_topology top, real mass2[], real qmol[])
111 for (i = 0; i < top.mols.nr; i++)
113 k = top.mols.index[i];
114 l = top.mols.index[i+1];
118 for (j = k; j < l; j++)
120 mtot += top.atoms.atom[j].m;
121 qtot += top.atoms.atom[j].q;
124 for (j = k; j < l; j++)
126 top.atoms.atom[j].q -= top.atoms.atom[j].m*qtot/mtot;
127 mass2[j] = top.atoms.atom[j].m/mtot;
144 if (fabs(qall) > 0.01)
146 printf("\n\nSystem not neutral (q=%f) will not calculate translational part of the dipole moment.\n", qall);
158 static void remove_jump(matrix box, int natoms, rvec xp[], rvec x[])
164 for (d = 0; d < DIM; d++)
166 hbox[d] = 0.5*box[d][d];
168 for (i = 0; i < natoms; i++)
170 for (m = DIM-1; m >= 0; m--)
172 while (x[i][m]-xp[i][m] <= -hbox[m])
174 for (d = 0; d <= m; d++)
176 x[i][d] += box[m][d];
179 while (x[i][m]-xp[i][m] > hbox[m])
181 for (d = 0; d <= m; d++)
183 x[i][d] -= box[m][d];
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[])
201 calc_box_center(ecenterRECT, box, center);
205 set_pbc(&pbc, ePBC, box);
211 for (i = 0; i < isize; i++)
215 k = top.mols.index[index0[i]];
216 l = top.mols.index[index0[i+1]];
217 for (j = k; j < l; j++)
219 svmul(mass2[j], fr[j], tmp);
225 svmul(qmol[k], mt1, mt1);
229 pbc_dx(&pbc, mt1, center, mt2);
230 svmul(qmol[k], mt2, mt1);
240 static real calceps(real prefactor, real md2, real mj2, real cor, real eps_rf, gmx_bool bCOR)
243 /* bCOR determines if the correlation is computed via
244 * static properties (FALSE) or the correlation integral (TRUE)
252 epsilon = md2-2.0*cor+mj2;
256 epsilon = md2+mj2+2.0*cor;
261 epsilon = 1.0+prefactor*epsilon;
266 epsilon = 2.0*eps_rf+1.0+2.0*eps_rf*prefactor*epsilon;
267 epsilon /= 2.0*eps_rf+1.0-prefactor*epsilon;
276 static real calc_cacf(FILE *fcacf, real prefactor, real cacf[], real time[], int nfr, int vfr[], int ei, int nshift)
284 real sigma_ret = 0.0;
294 rfr = (real) (nfr/nshift-i/nshift);
297 if (time[vfr[i]] <= time[vfr[ei]])
302 fprintf(fcacf, "%.3f\t%10.6g\t%10.6g\n", time[vfr[i]], cacf[i], sigma);
306 deltat = (time[vfr[i+1]]-time[vfr[i]]);
308 corint = 2.0*deltat*cacf[i]*prefactor;
309 if (i == 0 || (i+1) == nfr)
321 printf("Too less points.\n");
328 static void calc_mjdsp(FILE *fmjdsp, real prefactor, real dsp2[], real time[], int nfr, real refr[])
336 fprintf(fmjdsp, "#Prefactor fit E-H: 1 / 6.0*V*k_B*T: %g\n", prefactor);
340 for (i = 0; i < nfr; i++)
345 dsp2[i] *= prefactor/refr[i];
346 fprintf(fmjdsp, "%.3f\t%10.6g\n", time[i], dsp2[i]);
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)
364 int valloc, nalloc, nfr, nvfr, m, itrust = 0;
374 real prefactorav = 0.0;
375 real prefactor = 0.0;
377 real volume_av = 0.0;
378 real dk_s, dk_d, dk_f;
405 * Variables for the least-squares fit for Einstein-Helfand and Green-Kubo
415 gmx_rmpbc_t gpbc = NULL;
434 /* This is the main loop over frames */
449 gpbc = gmx_rmpbc_init(&top.idef, ePBC, fr.natoms);
454 refr = (real)(nfr+1);
459 srenew(time, nalloc);
461 srenew(mjdsp, nalloc);
462 srenew(dsp2, nalloc);
463 srenew(mtrans, nalloc);
464 srenew(xshfr, nalloc);
466 for (i = nfr; i < nalloc; i++)
468 clear_rvec(mjdsp[i]);
470 clear_rvec(mtrans[i]);
475 assert(time != NULL);
484 time[nfr] = fr.time-t0;
486 if (time[nfr] <= bfit)
490 if (time[nfr] <= efit)
500 remove_jump(fr.box, fr.natoms, xp, fr.x);
507 for (i = 0; i < fr.natoms; i++)
509 copy_rvec(fr.x[i], xp[i]);
514 gmx_rmpbc_trxfr(gpbc, &fr);
516 calc_mj(top, ePBC, fr.box, bNoJump, nmols, indexm, fr.x, mtrans[nfr], mass2, qmol);
518 for (i = 0; i < isize; i++)
521 svmul(top.atoms.atom[j].q, fr.x[j], fr.x[j]);
522 rvec_inc(mu[nfr], fr.x[j]);
525 /*if(mod(nfr,nshift)==0){*/
528 for (j = nfr; j >= 0; j--)
530 rvec_sub(mtrans[nfr], mtrans[j], tmp);
531 dsp2[nfr-j] += norm2(tmp);
549 srenew(cacf, valloc);
552 if (time[nfr] <= bvit)
556 if (time[nfr] <= evit)
561 clear_rvec(v0[nvfr]);
570 for (i = 0; i < isize; 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]);
578 fprintf(fcur, "%.3f\t%.6f\t%.6f\t%.6f\n", time[nfr], v0[nfr][XX], v0[nfr][YY], v0[nfr][ZZ]);
581 /*if(mod(nvfr,nshift)==0){*/
582 if (nvfr%nshift == 0)
584 for (j = nvfr; j >= 0; j--)
588 cacf[nvfr-j] += iprod(v0[nvfr], v0[j]);
592 djc[nvfr-j] += iprod(mu[vfr[j]], v0[nvfr]);
601 volume = det(fr.box);
604 rvec_inc(mja_tmp, mtrans[nfr]);
605 mjd += iprod(mu[nfr], mtrans[nfr]);
606 rvec_inc(mdvec, mu[nfr]);
608 mj2 += iprod(mtrans[nfr], mtrans[nfr]);
609 md2 += iprod(mu[nfr], mu[nfr]);
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);
618 while (read_next_frame(oenv, status, &fr));
620 gmx_rmpbc_done(gpbc);
625 prefactor /= 3.0*EPSILON0*volume_av*BOLTZ*temp;
628 prefactorav = E_CHARGE*E_CHARGE;
629 prefactorav /= volume_av*BOLTZMANN*temp*NANO*6.0;
631 fprintf(stderr, "Prefactor fit E-H: 1 / 6.0*V*k_B*T: %g\n", prefactorav);
633 calc_mjdsp(fmjdsp, prefactorav, dsp2, time, nfr, xshfr);
636 * Now we can average and calculate the correlation functions
644 svmul(1.0/refr, mdvec, mdvec);
645 svmul(1.0/refr, mja_tmp, mja_tmp);
647 mdav2 = norm2(mdvec);
649 mjdav = iprod(mdvec, mja_tmp);
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);
657 trust *= (real) nvfr;
658 itrust = (int) trust;
663 printf("\nCalculating M_D - current correlation integral ... \n");
664 corint = calc_cacf(mcor, prefactorav/EPSI0, djc, time, nvfr, vfr, ie, nshift);
671 printf("\nCalculating current autocorrelation ... \n");
672 sgk = calc_cacf(outf, prefactorav/PICO, cacf, time, nvfr, vfr, ie, nshift);
680 for (i = ii; i <= ie; i++)
683 xfit[i-ii] = log(time[vfr[i]]);
684 rtmp = fabs(cacf[i]);
685 yfit[i-ii] = log(rtmp);
689 lsq_y_ax_b(ie-ii, xfit, yfit, &sigma, &malt, &err, &rest);
695 malt *= prefactorav*2.0e12/sigma;
705 /* Calculation of the dielectric constant */
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");
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");
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);
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);
731 if (bACF && (ii < nvfr))
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);
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);
745 for (i = bi; i <= ei; i++)
747 xfit[i-bi] = time[i];
748 yfit[i-bi] = dsp2[i];
751 lsq_y_ax_b(ei-bi, xfit, yfit, &sigma, &malt, &err, &rest);
754 dk_d = calceps(prefactor, md2, 0.5*malt/prefactorav, corint, eps_rf, TRUE);
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);
766 fprintf(stderr, "Too few points for a fit.\n");
792 int gmx_current(int argc, char *argv[])
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;
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."}
828 char **grpname = NULL;
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 }
871 #define NFILE asize(fnm)
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."
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.",
893 "Option [TT]-temp[tt] sets the temperature required for the computation of the static dielectric constant.",
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).",
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."
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))
917 bACF = opt2bSet("-caf", NFILE, fnm);
918 bINT = opt2bSet("-mc", NFILE, fnm);
920 bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC, &xtop, &vtop, box, TRUE);
926 indexfn = ftp2fn_null(efNDX, NFILE, fnm);
929 get_index(&(top.atoms), indexfn, 1, &isize, &index0, grpname);
931 flags = flags | TRX_READ_X | TRX_READ_V;
933 read_first_frame(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &fr, flags);
935 snew(mass2, top.atoms.nr);
936 snew(qmol, top.atoms.nr);
938 bNEU = precalc(top, mass2, qmol);
943 for (i = 0; i < isize; i++)
945 indexm[i] = index0[i];
951 index_atom2mol(&nmols, indexm, &top.mols);
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");
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");
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");
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");
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",
991 /* System information is read and prepared, dielectric() processes the frames
992 * and calculates the requested quantities */
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);
1000 gmx_ffclose(fmjdsp);