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.
40 #include "gromacs/commandline/pargs.h"
41 #include "gromacs/legacyheaders/typedefs.h"
42 #include "gromacs/utility/smalloc.h"
43 #include "gromacs/math/vec.h"
44 #include "gromacs/fileio/tpxio.h"
45 #include "gromacs/fileio/trxio.h"
46 #include "gromacs/fileio/xvgr.h"
47 #include "gromacs/pbcutil/rmpbc.h"
48 #include "gromacs/pbcutil/pbc.h"
49 #include "gromacs/math/units.h"
50 #include "gromacs/topology/index.h"
51 #include "gromacs/statistics/statistics.h"
53 #include "gromacs/legacyheaders/macros.h"
55 #include "gromacs/utility/fatalerror.h"
57 #define SQR(x) (pow(x, 2.0))
58 #define EPSI0 (EPSILON0*E_CHARGE*E_CHARGE*AVOGADRO/(KILO*NANO)) /* EPSILON0 in SI units */
60 static void index_atom2mol(int *n, int *index, t_block *mols)
62 int nat, i, nmol, mol, j;
70 while (index[i] > mols->index[mol])
75 gmx_fatal(FARGS, "Atom index out of range: %d", index[i]+1);
78 for (j = mols->index[mol]; j < mols->index[mol+1]; j++)
80 if (i >= nat || index[i] != j)
82 gmx_fatal(FARGS, "The index group does not consist of whole molecules");
89 fprintf(stderr, "\nSplit group of %d atoms into %d molecules\n", nat, nmol);
94 static gmx_bool precalc(t_topology top, real mass2[], real qmol[])
109 for (i = 0; i < top.mols.nr; i++)
111 k = top.mols.index[i];
112 l = top.mols.index[i+1];
116 for (j = k; j < l; j++)
118 mtot += top.atoms.atom[j].m;
119 qtot += top.atoms.atom[j].q;
122 for (j = k; j < l; j++)
124 top.atoms.atom[j].q -= top.atoms.atom[j].m*qtot/mtot;
125 mass2[j] = top.atoms.atom[j].m/mtot;
142 if (fabs(qall) > 0.01)
144 printf("\n\nSystem not neutral (q=%f) will not calculate translational part of the dipole moment.\n", qall);
156 static void remove_jump(matrix box, int natoms, rvec xp[], rvec x[])
162 for (d = 0; d < DIM; d++)
164 hbox[d] = 0.5*box[d][d];
166 for (i = 0; i < natoms; i++)
168 for (m = DIM-1; m >= 0; m--)
170 while (x[i][m]-xp[i][m] <= -hbox[m])
172 for (d = 0; d <= m; d++)
174 x[i][d] += box[m][d];
177 while (x[i][m]-xp[i][m] > hbox[m])
179 for (d = 0; d <= m; d++)
181 x[i][d] -= box[m][d];
188 static void calc_mj(t_topology top, int ePBC, matrix box, gmx_bool bNoJump, int isize, int index0[], \
189 rvec fr[], rvec mj, real mass2[], real qmol[])
199 calc_box_center(ecenterRECT, box, center);
203 set_pbc(&pbc, ePBC, box);
209 for (i = 0; i < isize; i++)
213 k = top.mols.index[index0[i]];
214 l = top.mols.index[index0[i+1]];
215 for (j = k; j < l; j++)
217 svmul(mass2[j], fr[j], tmp);
223 svmul(qmol[k], mt1, mt1);
227 pbc_dx(&pbc, mt1, center, mt2);
228 svmul(qmol[k], mt2, mt1);
238 static real calceps(real prefactor, real md2, real mj2, real cor, real eps_rf, gmx_bool bCOR)
241 /* bCOR determines if the correlation is computed via
242 * static properties (FALSE) or the correlation integral (TRUE)
250 epsilon = md2-2.0*cor+mj2;
254 epsilon = md2+mj2+2.0*cor;
259 epsilon = 1.0+prefactor*epsilon;
264 epsilon = 2.0*eps_rf+1.0+2.0*eps_rf*prefactor*epsilon;
265 epsilon /= 2.0*eps_rf+1.0-prefactor*epsilon;
274 static real calc_cacf(FILE *fcacf, real prefactor, real cacf[], real time[], int nfr, int vfr[], int ei, int nshift)
282 real sigma_ret = 0.0;
292 rfr = (real) (nfr/nshift-i/nshift);
295 if (time[vfr[i]] <= time[vfr[ei]])
300 fprintf(fcacf, "%.3f\t%10.6g\t%10.6g\n", time[vfr[i]], cacf[i], sigma);
304 deltat = (time[vfr[i+1]]-time[vfr[i]]);
306 corint = 2.0*deltat*cacf[i]*prefactor;
307 if (i == 0 || (i+1) == nfr)
319 printf("Too less points.\n");
326 static void calc_mjdsp(FILE *fmjdsp, real prefactor, real dsp2[], real time[], int nfr, real refr[])
334 fprintf(fmjdsp, "#Prefactor fit E-H: 1 / 6.0*V*k_B*T: %g\n", prefactor);
338 for (i = 0; i < nfr; i++)
343 dsp2[i] *= prefactor/refr[i];
344 fprintf(fmjdsp, "%.3f\t%10.6g\n", time[i], dsp2[i]);
353 static void dielectric(FILE *fmj, FILE *fmd, FILE *outf, FILE *fcur, FILE *mcor,
354 FILE *fmjdsp, gmx_bool bNoJump, gmx_bool bACF, gmx_bool bINT,
355 int ePBC, t_topology top, t_trxframe fr, real temp,
356 real trust, real bfit, real efit, real bvit, real evit,
357 t_trxstatus *status, int isize, int nmols, int nshift,
358 atom_id *index0, int indexm[], real mass2[],
359 real qmol[], real eps_rf, const output_env_t oenv)
362 int valloc, nalloc, nfr, nvfr, m, itrust = 0;
372 real prefactorav = 0.0;
373 real prefactor = 0.0;
375 real volume_av = 0.0;
376 real dk_s, dk_d, dk_f;
403 * Variables for the least-squares fit for Einstein-Helfand and Green-Kubo
413 gmx_rmpbc_t gpbc = NULL;
432 /* This is the main loop over frames */
447 gpbc = gmx_rmpbc_init(&top.idef, ePBC, fr.natoms);
452 refr = (real)(nfr+1);
457 srenew(time, nalloc);
459 srenew(mjdsp, nalloc);
460 srenew(dsp2, nalloc);
461 srenew(mtrans, nalloc);
462 srenew(xshfr, nalloc);
464 for (i = nfr; i < nalloc; i++)
466 clear_rvec(mjdsp[i]);
468 clear_rvec(mtrans[i]);
473 assert(time != NULL);
482 time[nfr] = fr.time-t0;
484 if (time[nfr] <= bfit)
488 if (time[nfr] <= efit)
498 remove_jump(fr.box, fr.natoms, xp, fr.x);
505 for (i = 0; i < fr.natoms; i++)
507 copy_rvec(fr.x[i], xp[i]);
512 gmx_rmpbc_trxfr(gpbc, &fr);
514 calc_mj(top, ePBC, fr.box, bNoJump, nmols, indexm, fr.x, mtrans[nfr], mass2, qmol);
516 for (i = 0; i < isize; i++)
519 svmul(top.atoms.atom[j].q, fr.x[j], fr.x[j]);
520 rvec_inc(mu[nfr], fr.x[j]);
523 /*if(mod(nfr,nshift)==0){*/
526 for (j = nfr; j >= 0; j--)
528 rvec_sub(mtrans[nfr], mtrans[j], tmp);
529 dsp2[nfr-j] += norm2(tmp);
547 srenew(cacf, valloc);
550 if (time[nfr] <= bvit)
554 if (time[nfr] <= evit)
559 clear_rvec(v0[nvfr]);
568 for (i = 0; i < isize; i++)
571 svmul(mass2[j], fr.v[j], fr.v[j]);
572 svmul(qmol[j], fr.v[j], fr.v[j]);
573 rvec_inc(v0[nvfr], fr.v[j]);
576 fprintf(fcur, "%.3f\t%.6f\t%.6f\t%.6f\n", time[nfr], v0[nfr][XX], v0[nfr][YY], v0[nfr][ZZ]);
579 /*if(mod(nvfr,nshift)==0){*/
580 if (nvfr%nshift == 0)
582 for (j = nvfr; j >= 0; j--)
586 cacf[nvfr-j] += iprod(v0[nvfr], v0[j]);
590 djc[nvfr-j] += iprod(mu[vfr[j]], v0[nvfr]);
599 volume = det(fr.box);
602 rvec_inc(mja_tmp, mtrans[nfr]);
603 mjd += iprod(mu[nfr], mtrans[nfr]);
604 rvec_inc(mdvec, mu[nfr]);
606 mj2 += iprod(mtrans[nfr], mtrans[nfr]);
607 md2 += iprod(mu[nfr], mu[nfr]);
609 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);
610 fprintf(fmd, "%.3f\t%8.5f\t%8.5f\t%8.5f\t%8.5f\t%8.5f\n", \
611 time[nfr], mu[nfr][XX], mu[nfr][YY], mu[nfr][ZZ], md2/refr, norm(mdvec)/refr);
616 while (read_next_frame(oenv, status, &fr));
618 gmx_rmpbc_done(gpbc);
623 prefactor /= 3.0*EPSILON0*volume_av*BOLTZ*temp;
626 prefactorav = E_CHARGE*E_CHARGE;
627 prefactorav /= volume_av*BOLTZMANN*temp*NANO*6.0;
629 fprintf(stderr, "Prefactor fit E-H: 1 / 6.0*V*k_B*T: %g\n", prefactorav);
631 calc_mjdsp(fmjdsp, prefactorav, dsp2, time, nfr, xshfr);
634 * Now we can average and calculate the correlation functions
642 svmul(1.0/refr, mdvec, mdvec);
643 svmul(1.0/refr, mja_tmp, mja_tmp);
645 mdav2 = norm2(mdvec);
647 mjdav = iprod(mdvec, mja_tmp);
650 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);
651 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);
655 trust *= (real) nvfr;
656 itrust = (int) trust;
661 printf("\nCalculating M_D - current correlation integral ... \n");
662 corint = calc_cacf(mcor, prefactorav/EPSI0, djc, time, nvfr, vfr, ie, nshift);
669 printf("\nCalculating current autocorrelation ... \n");
670 sgk = calc_cacf(outf, prefactorav/PICO, cacf, time, nvfr, vfr, ie, nshift);
678 for (i = ii; i <= ie; i++)
681 xfit[i-ii] = log(time[vfr[i]]);
682 rtmp = fabs(cacf[i]);
683 yfit[i-ii] = log(rtmp);
687 lsq_y_ax_b(ie-ii, xfit, yfit, &sigma, &malt, &err, &rest);
693 malt *= prefactorav*2.0e12/sigma;
703 /* Calculation of the dielectric constant */
705 fprintf(stderr, "\n********************************************\n");
706 dk_s = calceps(prefactor, md2, mj2, mjd, eps_rf, FALSE);
707 fprintf(stderr, "\nAbsolute values:\n epsilon=%f\n", dk_s);
708 fprintf(stderr, " <M_D^2> , <M_J^2>, <(M_J*M_D)^2>: (%f, %f, %f)\n\n", md2, mj2, mjd);
709 fprintf(stderr, "********************************************\n");
712 dk_f = calceps(prefactor, md2-mdav2, mj2-mj, mjd-mjdav, eps_rf, FALSE);
713 fprintf(stderr, "\n\nFluctuations:\n epsilon=%f\n\n", dk_f);
714 fprintf(stderr, "\n deltaM_D , deltaM_J, deltaM_JD: (%f, %f, %f)\n", md2-mdav2, mj2-mj, mjd-mjdav);
715 fprintf(stderr, "\n********************************************\n");
718 dk_d = calceps(prefactor, md2-mdav2, mj2-mj, corint, eps_rf, TRUE);
719 fprintf(stderr, "\nStatic dielectric constant using integral and fluctuations: %f\n", dk_d);
720 fprintf(stderr, "\n < M_JM_D > via integral: %.3f\n", -1.0*corint);
723 fprintf(stderr, "\n***************************************************");
724 fprintf(stderr, "\n\nAverage volume V=%f nm^3 at T=%f K\n", volume_av, temp);
725 fprintf(stderr, "and corresponding refactor 1.0 / 3.0*V*k_B*T*EPSILON_0: %f \n", prefactor);
729 if (bACF && (ii < nvfr))
731 fprintf(stderr, "Integral and integrated fit to the current acf yields at t=%f:\n", time[vfr[ii]]);
732 fprintf(stderr, "sigma=%8.3f (pure integral: %.3f)\n", sgk-malt*pow(time[vfr[ii]], sigma), sgk);
737 fprintf(stderr, "\nStart fit at %f ps (%f).\n", time[bi], bfit);
738 fprintf(stderr, "End fit at %f ps (%f).\n\n", time[ei], efit);
743 for (i = bi; i <= ei; i++)
745 xfit[i-bi] = time[i];
746 yfit[i-bi] = dsp2[i];
749 lsq_y_ax_b(ei-bi, xfit, yfit, &sigma, &malt, &err, &rest);
752 dk_d = calceps(prefactor, md2, 0.5*malt/prefactorav, corint, eps_rf, TRUE);
754 fprintf(stderr, "Einstein-Helfand fit to the MSD of the translational dipole moment yields:\n\n");
755 fprintf(stderr, "sigma=%.4f\n", sigma);
756 fprintf(stderr, "translational fraction of M^2: %.4f\n", 0.5*malt/prefactorav);
757 fprintf(stderr, "Dielectric constant using EH: %.4f\n", dk_d);
764 fprintf(stderr, "Too few points for a fit.\n");
790 int gmx_current(int argc, char *argv[])
793 static int nshift = 1000;
794 static real temp = 300.0;
795 static real trust = 0.25;
796 static real eps_rf = 0.0;
797 static gmx_bool bNoJump = TRUE;
798 static real bfit = 100.0;
799 static real bvit = 0.5;
800 static real efit = 400.0;
801 static real evit = 5.0;
803 { "-sh", FALSE, etINT, {&nshift},
804 "Shift of the frames for averaging the correlation functions and the mean-square displacement."},
805 { "-nojump", FALSE, etBOOL, {&bNoJump},
806 "Removes jumps of atoms across the box."},
807 { "-eps", FALSE, etREAL, {&eps_rf},
808 "Dielectric constant of the surrounding medium. The value zero corresponds to infinity (tin-foil boundary conditions)."},
809 { "-bfit", FALSE, etREAL, {&bfit},
810 "Begin of the fit of the straight line to the MSD of the translational fraction of the dipole moment."},
811 { "-efit", FALSE, etREAL, {&efit},
812 "End of the fit of the straight line to the MSD of the translational fraction of the dipole moment."},
813 { "-bvit", FALSE, etREAL, {&bvit},
814 "Begin of the fit of the current autocorrelation function to a*t^b."},
815 { "-evit", FALSE, etREAL, {&evit},
816 "End of the fit of the current autocorrelation function to a*t^b."},
817 { "-tr", FALSE, etREAL, {&trust},
818 "Fraction of the trajectory taken into account for the integral."},
819 { "-temp", FALSE, etREAL, {&temp},
820 "Temperature for calculating epsilon."}
826 char **grpname = NULL;
858 { efTPS, NULL, NULL, ffREAD }, /* this is for the topology */
859 { efNDX, NULL, NULL, ffOPTRD },
860 { efTRX, "-f", NULL, ffREAD }, /* and this for the trajectory */
861 { efXVG, "-o", "current", ffWRITE },
862 { efXVG, "-caf", "caf", ffOPTWR },
863 { efXVG, "-dsp", "dsp", ffWRITE },
864 { efXVG, "-md", "md", ffWRITE },
865 { efXVG, "-mj", "mj", ffWRITE },
866 { efXVG, "-mc", "mc", ffOPTWR }
869 #define NFILE asize(fnm)
872 const char *desc[] = {
873 "[THISMODULE] is a tool for calculating the current autocorrelation function, the correlation",
874 "of the rotational and translational dipole moment of the system, and the resulting static",
875 "dielectric constant. To obtain a reasonable result, the index group has to be neutral.",
876 "Furthermore, the routine is capable of extracting the static conductivity from the current ",
877 "autocorrelation function, if velocities are given. Additionally, an Einstein-Helfand fit ",
878 "can be used to obtain the static conductivity."
880 "The flag [TT]-caf[tt] is for the output of the current autocorrelation function and [TT]-mc[tt] writes the",
881 "correlation of the rotational and translational part of the dipole moment in the corresponding",
882 "file. However, this option is only available for trajectories containing velocities.",
883 "Options [TT]-sh[tt] and [TT]-tr[tt] are responsible for the averaging and integration of the",
884 "autocorrelation functions. Since averaging proceeds by shifting the starting point",
885 "through the trajectory, the shift can be modified with [TT]-sh[tt] to enable the choice of uncorrelated",
886 "starting points. Towards the end, statistical inaccuracy grows and integrating the",
887 "correlation function only yields reliable values until a certain point, depending on",
888 "the number of frames. The option [TT]-tr[tt] controls the region of the integral taken into account",
889 "for calculating the static dielectric constant.",
891 "Option [TT]-temp[tt] sets the temperature required for the computation of the static dielectric constant.",
893 "Option [TT]-eps[tt] controls the dielectric constant of the surrounding medium for simulations using",
894 "a Reaction Field or dipole corrections of the Ewald summation ([TT]-eps[tt]=0 corresponds to",
895 "tin-foil boundary conditions).",
897 "[TT]-[no]nojump[tt] unfolds the coordinates to allow free diffusion. This is required to get a continuous",
898 "translational dipole moment, required for the Einstein-Helfand fit. The results from the fit allow",
899 "the determination of the dielectric constant for system of charged molecules. However, it is also possible to extract",
900 "the dielectric constant from the fluctuations of the total dipole moment in folded coordinates. But this",
901 "option has to be used with care, since only very short time spans fulfill the approximation that the density",
902 "of the molecules is approximately constant and the averages are already converged. To be on the safe side,",
903 "the dielectric constant should be calculated with the help of the Einstein-Helfand method for",
904 "the translational part of the dielectric constant."
908 /* At first the arguments will be parsed and the system information processed */
909 if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW,
910 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
915 bACF = opt2bSet("-caf", NFILE, fnm);
916 bINT = opt2bSet("-mc", NFILE, fnm);
918 bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC, &xtop, &vtop, box, TRUE);
924 indexfn = ftp2fn_null(efNDX, NFILE, fnm);
927 get_index(&(top.atoms), indexfn, 1, &isize, &index0, grpname);
929 flags = flags | TRX_READ_X | TRX_READ_V;
931 read_first_frame(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &fr, flags);
933 snew(mass2, top.atoms.nr);
934 snew(qmol, top.atoms.nr);
936 bNEU = precalc(top, mass2, qmol);
941 for (i = 0; i < isize; i++)
943 indexm[i] = index0[i];
949 index_atom2mol(&nmols, indexm, &top.mols);
955 outf = xvgropen(opt2fn("-caf", NFILE, fnm),
956 "Current autocorrelation function", output_env_get_xvgr_tlabel(oenv),
957 "ACF (e nm/ps)\\S2", oenv);
958 fprintf(outf, "# time\t acf\t average \t std.dev\n");
960 fcur = xvgropen(opt2fn("-o", NFILE, fnm),
961 "Current", output_env_get_xvgr_tlabel(oenv), "J(t) (e nm/ps)", oenv);
962 fprintf(fcur, "# time\t Jx\t Jy \t J_z \n");
965 mcor = xvgropen(opt2fn("-mc", NFILE, fnm),
966 "M\\sD\\N - current autocorrelation function",
967 output_env_get_xvgr_tlabel(oenv),
968 "< M\\sD\\N (0)\\c7\\CJ(t) > (e nm/ps)\\S2", oenv);
969 fprintf(mcor, "# time\t M_D(0) J(t) acf \t Integral acf\n");
973 fmj = xvgropen(opt2fn("-mj", NFILE, fnm),
974 "Averaged translational part of M", output_env_get_xvgr_tlabel(oenv),
975 "< M\\sJ\\N > (enm)", oenv);
976 fprintf(fmj, "# time\t x\t y \t z \t average of M_J^2 \t std.dev\n");
977 fmd = xvgropen(opt2fn("-md", NFILE, fnm),
978 "Averaged rotational part of M", output_env_get_xvgr_tlabel(oenv),
979 "< M\\sD\\N > (enm)", oenv);
980 fprintf(fmd, "# time\t x\t y \t z \t average of M_D^2 \t std.dev\n");
982 fmjdsp = xvgropen(opt2fn("-dsp", NFILE, fnm),
983 "MSD of the squared translational dipole moment M",
984 output_env_get_xvgr_tlabel(oenv),
985 "<|M\\sJ\\N(t)-M\\sJ\\N(0)|\\S2\\N > / 6.0*V*k\\sB\\N*T / Sm\\S-1\\Nps\\S-1\\N",
989 /* System information is read and prepared, dielectric() processes the frames
990 * and calculates the requested quantities */
992 dielectric(fmj, fmd, outf, fcur, mcor, fmjdsp, bNoJump, bACF, bINT, ePBC, top, fr,
993 temp, trust, bfit, efit, bvit, evit, status, isize, nmols, nshift,
994 index0, indexm, mass2, qmol, eps_rf, oenv);