/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2008,2009,2010,2011,2012,2013,2014,2015,2017,2018, by the GROMACS development team, led by
+ * Copyright (c) 2008,2009,2010,2011,2012,2013,2014,2015,2017,2018,2019, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#include "gromacs/utility/gmxassert.h"
#include "gromacs/utility/smalloc.h"
-#define EPSI0 (EPSILON0*E_CHARGE*E_CHARGE*AVOGADRO/(KILO*NANO)) /* EPSILON0 in SI units */
+#define EPSI0 (EPSILON0 * E_CHARGE * E_CHARGE * AVOGADRO / (KILO * NANO)) /* EPSILON0 in SI units */
-static void index_atom2mol(int *n, int *index, t_block *mols)
+static void index_atom2mol(int* n, int* index, t_block* mols)
{
int nat, i, nmol, mol, j;
mol++;
if (mol >= mols->nr)
{
- gmx_fatal(FARGS, "Atom index out of range: %d", index[i]+1);
+ gmx_fatal(FARGS, "Atom index out of range: %d", index[i] + 1);
}
}
- for (j = mols->index[mol]; j < mols->index[mol+1]; j++)
+ for (j = mols->index[mol]; j < mols->index[mol + 1]; j++)
{
if (i >= nat || index[i] != j)
{
qall = 0.0;
-
for (i = 0; i < top.mols.nr; i++)
{
k = top.mols.index[i];
- l = top.mols.index[i+1];
+ l = top.mols.index[i + 1];
mtot = 0.0;
qtot = 0.0;
for (j = k; j < l; j++)
{
- top.atoms.atom[j].q -= top.atoms.atom[j].m*qtot/mtot;
- mass2[j] = top.atoms.atom[j].m/mtot;
- qmol[j] = qtot;
+ top.atoms.atom[j].q -= top.atoms.atom[j].m * qtot / mtot;
+ mass2[j] = top.atoms.atom[j].m / mtot;
+ qmol[j] = qtot;
}
if (std::abs(qall) > 0.01)
{
- printf("\n\nSystem not neutral (q=%f) will not calculate translational part of the dipole moment.\n", qall);
+ printf("\n\nSystem not neutral (q=%f) will not calculate translational part of the dipole "
+ "moment.\n",
+ qall);
bNEU = FALSE;
}
else
}
return bNEU;
-
}
static void remove_jump(matrix box, int natoms, rvec xp[], rvec x[])
for (d = 0; d < DIM; d++)
{
- hbox[d] = 0.5*box[d][d];
+ hbox[d] = 0.5 * box[d][d];
}
for (i = 0; i < natoms; i++)
{
- for (m = DIM-1; m >= 0; m--)
+ for (m = DIM - 1; m >= 0; m--)
{
- while (x[i][m]-xp[i][m] <= -hbox[m])
+ while (x[i][m] - xp[i][m] <= -hbox[m])
{
for (d = 0; d <= m; d++)
{
x[i][d] += box[m][d];
}
}
- while (x[i][m]-xp[i][m] > hbox[m])
+ while (x[i][m] - xp[i][m] > hbox[m])
{
for (d = 0; d <= m; d++)
{
}
}
-static void calc_mj(t_topology top, int ePBC, matrix box, gmx_bool bNoJump, int isize, const int index0[], \
- rvec fr[], rvec mj, real mass2[], real qmol[])
+static void calc_mj(t_topology top,
+ int ePBC,
+ matrix box,
+ gmx_bool bNoJump,
+ int isize,
+ const int index0[],
+ rvec fr[],
+ rvec mj,
+ real mass2[],
+ real qmol[])
{
int i, j, k, l;
clear_rvec(mt1);
clear_rvec(mt2);
k = top.mols.index[index0[i]];
- l = top.mols.index[index0[i+1]];
+ l = top.mols.index[index0[i + 1]];
for (j = k; j < l; j++)
{
svmul(mass2[j], fr[j], tmp);
}
rvec_inc(mj, mt1);
-
}
-
}
if (bCOR)
{
- epsilon = md2-2.0*cor+mj2;
+ epsilon = md2 - 2.0 * cor + mj2;
}
else
{
- epsilon = md2+mj2+2.0*cor;
+ epsilon = md2 + mj2 + 2.0 * cor;
}
if (eps_rf == 0.0)
{
- epsilon = 1.0+prefactor*epsilon;
-
+ epsilon = 1.0 + prefactor * epsilon;
}
else
{
- epsilon = 2.0*eps_rf+1.0+2.0*eps_rf*prefactor*epsilon;
- epsilon /= 2.0*eps_rf+1.0-prefactor*epsilon;
+ epsilon = 2.0 * eps_rf + 1.0 + 2.0 * eps_rf * prefactor * epsilon;
+ epsilon /= 2.0 * eps_rf + 1.0 - prefactor * epsilon;
}
return epsilon;
-
}
-static real calc_cacf(FILE *fcacf, real prefactor, real cacf[], real time[], int nfr, const int vfr[], int ei, int nshift)
+static real calc_cacf(FILE* fcacf, real prefactor, real cacf[], real time[], int nfr, const int vfr[], int ei, int nshift)
{
int i;
{
real corint;
- rfr = static_cast<real>(nfr+i)/nshift;
+ rfr = static_cast<real>(nfr + i) / nshift;
cacf[i] /= rfr;
if (time[vfr[i]] <= time[vfr[ei]])
fprintf(fcacf, "%.3f\t%10.6g\t%10.6g\n", time[vfr[i]], cacf[i], sigma);
- if ((i+1) < nfr)
+ if ((i + 1) < nfr)
{
- deltat = (time[vfr[i+1]]-time[vfr[i]]);
+ deltat = (time[vfr[i + 1]] - time[vfr[i]]);
}
- corint = 2.0*deltat*cacf[i]*prefactor;
- if (i == 0 || (i+1) == nfr)
+ corint = 2.0 * deltat * cacf[i] * prefactor;
+ if (i == 0 || (i + 1) == nfr)
{
corint *= 0.5;
}
i++;
}
-
}
else
{
}
return sigma_ret;
-
}
-static void calc_mjdsp(FILE *fmjdsp, real prefactor, real dsp2[], real time[], int nfr, const real refr[])
+static void calc_mjdsp(FILE* fmjdsp, real prefactor, real dsp2[], real time[], int nfr, const real refr[])
{
- int i;
+ int i;
fprintf(fmjdsp, "#Prefactor fit E-H: 1 / 6.0*V*k_B*T: %g\n", prefactor);
if (refr[i] != 0.0)
{
- dsp2[i] *= prefactor/refr[i];
+ dsp2[i] *= prefactor / refr[i];
fprintf(fmjdsp, "%.3f\t%10.6g\n", time[i], dsp2[i]);
}
-
-
}
-
}
-static void dielectric(FILE *fmj, FILE *fmd, FILE *outf, FILE *fcur, FILE *mcor,
- FILE *fmjdsp, gmx_bool bNoJump, gmx_bool bACF, gmx_bool bINT,
- int ePBC, t_topology top, t_trxframe fr, real temp,
- real bfit, real efit, real bvit, real evit,
- t_trxstatus *status, int isize, int nmols, int nshift,
- const int *index0, int indexm[], real mass2[],
- real qmol[], real eps_rf, const gmx_output_env_t *oenv)
+static void dielectric(FILE* fmj,
+ FILE* fmd,
+ FILE* outf,
+ FILE* fcur,
+ FILE* mcor,
+ FILE* fmjdsp,
+ gmx_bool bNoJump,
+ gmx_bool bACF,
+ gmx_bool bINT,
+ int ePBC,
+ t_topology top,
+ t_trxframe fr,
+ real temp,
+ real bfit,
+ real efit,
+ real bvit,
+ real evit,
+ t_trxstatus* status,
+ int isize,
+ int nmols,
+ int nshift,
+ const int* index0,
+ int indexm[],
+ real mass2[],
+ real qmol[],
+ real eps_rf,
+ const gmx_output_env_t* oenv)
{
- int i, j;
- int valloc, nalloc, nfr, nvfr;
- int vshfr;
- real *xshfr = nullptr;
- int *vfr = nullptr;
- real refr = 0.0;
- real *cacf = nullptr;
- real *time = nullptr;
- real *djc = nullptr;
- real corint = 0.0;
- real prefactorav = 0.0;
- real prefactor = 0.0;
- real volume;
- real volume_av = 0.0;
- real dk_s, dk_d, dk_f;
- real mj = 0.0;
- real mj2 = 0.0;
- real mjd = 0.0;
- real mjdav = 0.0;
- real md2 = 0.0;
- real mdav2 = 0.0;
- real sgk;
- rvec mja_tmp;
- rvec mjd_tmp;
- rvec mdvec;
- rvec *mu = nullptr;
- rvec *xp = nullptr;
- rvec *v0 = nullptr;
- rvec *mjdsp = nullptr;
- real *dsp2 = nullptr;
- real t0;
- real rtmp;
-
- rvec tmp;
- rvec *mtrans = nullptr;
+ int i, j;
+ int valloc, nalloc, nfr, nvfr;
+ int vshfr;
+ real* xshfr = nullptr;
+ int* vfr = nullptr;
+ real refr = 0.0;
+ real* cacf = nullptr;
+ real* time = nullptr;
+ real* djc = nullptr;
+ real corint = 0.0;
+ real prefactorav = 0.0;
+ real prefactor = 0.0;
+ real volume;
+ real volume_av = 0.0;
+ real dk_s, dk_d, dk_f;
+ real mj = 0.0;
+ real mj2 = 0.0;
+ real mjd = 0.0;
+ real mjdav = 0.0;
+ real md2 = 0.0;
+ real mdav2 = 0.0;
+ real sgk;
+ rvec mja_tmp;
+ rvec mjd_tmp;
+ rvec mdvec;
+ rvec* mu = nullptr;
+ rvec* xp = nullptr;
+ rvec* v0 = nullptr;
+ rvec* mjdsp = nullptr;
+ real* dsp2 = nullptr;
+ real t0;
+ real rtmp;
+
+ rvec tmp;
+ rvec* mtrans = nullptr;
/*
* Variables for the least-squares fit for Einstein-Helfand and Green-Kubo
*/
- int bi, ei, ie, ii;
- real rest = 0.0;
- real sigma = 0.0;
- real malt = 0.0;
- real err = 0.0;
- real *xfit;
- real *yfit;
- gmx_rmpbc_t gpbc = nullptr;
+ int bi, ei, ie, ii;
+ real rest = 0.0;
+ real sigma = 0.0;
+ real malt = 0.0;
+ real err = 0.0;
+ real* xfit;
+ real* yfit;
+ gmx_rmpbc_t gpbc = nullptr;
/*
* indices for EH
do
{
- refr = (nfr+1);
+ refr = (nfr + 1);
if (nfr >= nalloc)
{
if (nfr == 0)
{
t0 = fr.time;
-
}
- time[nfr] = fr.time-t0;
+ time[nfr] = fr.time - t0;
if (time[nfr] <= bfit)
{
{
copy_rvec(fr.x[i], xp[i]);
}
-
}
gmx_rmpbc_trxfr(gpbc, &fr);
}
/*if(mod(nfr,nshift)==0){*/
- if (nfr%nshift == 0)
+ if (nfr % nshift == 0)
{
for (j = nfr; j >= 0; j--)
{
rvec_sub(mtrans[nfr], mtrans[j], tmp);
- dsp2[nfr-j] += norm2(tmp);
- xshfr[nfr-j] += 1.0;
+ dsp2[nfr - j] += norm2(tmp);
+ xshfr[nfr - j] += 1.0;
}
}
if (bACF || bINT)
{
/*if(mod(nvfr,nshift)==0){*/
- if (nvfr%nshift == 0)
+ if (nvfr % nshift == 0)
{
for (j = nvfr; j >= 0; j--)
{
if (bACF)
{
- cacf[nvfr-j] += iprod(v0[nvfr], v0[j]);
+ cacf[nvfr - j] += iprod(v0[nvfr], v0[j]);
}
if (bINT)
{
- djc[nvfr-j] += iprod(mu[vfr[j]], v0[nvfr]);
+ djc[nvfr - j] += iprod(mu[vfr[j]], v0[nvfr]);
}
}
vshfr++;
nvfr++;
}
- volume = det(fr.box);
+ volume = det(fr.box);
volume_av += volume;
rvec_inc(mja_tmp, mtrans[nfr]);
mj2 += iprod(mtrans[nfr], mtrans[nfr]);
md2 += iprod(mu[nfr], mu[nfr]);
- 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);
- fprintf(fmd, "%.3f\t%8.5f\t%8.5f\t%8.5f\t%8.5f\t%8.5f\n", \
- time[nfr], mu[nfr][XX], mu[nfr][YY], mu[nfr][ZZ], md2/refr, norm(mdvec)/refr);
+ 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);
+ fprintf(fmd, "%.3f\t%8.5f\t%8.5f\t%8.5f\t%8.5f\t%8.5f\n", time[nfr], mu[nfr][XX],
+ mu[nfr][YY], mu[nfr][ZZ], md2 / refr, norm(mdvec) / refr);
nfr++;
- }
- while (read_next_frame(oenv, status, &fr));
+ } while (read_next_frame(oenv, status, &fr));
gmx_rmpbc_done(gpbc);
volume_av /= refr;
- prefactor = 1.0;
- prefactor /= 3.0*EPSILON0*volume_av*BOLTZ*temp;
+ prefactor = 1.0;
+ prefactor /= 3.0 * EPSILON0 * volume_av * BOLTZ * temp;
- prefactorav = E_CHARGE*E_CHARGE;
- prefactorav /= volume_av*BOLTZMANN*temp*NANO*6.0;
+ prefactorav = E_CHARGE * E_CHARGE;
+ prefactorav /= volume_av * BOLTZMANN * temp * NANO * 6.0;
fprintf(stderr, "Prefactor fit E-H: 1 / 6.0*V*k_B*T: %g\n", prefactorav);
mjd /= refr;
md2 /= refr;
- svmul(1.0/refr, mdvec, mdvec);
- svmul(1.0/refr, mja_tmp, mja_tmp);
+ svmul(1.0 / refr, mdvec, mdvec);
+ svmul(1.0 / refr, mja_tmp, mja_tmp);
mdav2 = norm2(mdvec);
mj = norm2(mja_tmp);
mjdav = iprod(mdvec, mja_tmp);
- 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);
- 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);
+ 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);
+ 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);
if (v0 != nullptr)
{
{
printf("\nCalculating M_D - current correlation integral ... \n");
- corint = calc_cacf(mcor, prefactorav/EPSI0, djc, time, nvfr, vfr, ie, nshift);
-
+ corint = calc_cacf(mcor, prefactorav / EPSI0, djc, time, nvfr, vfr, ie, nshift);
}
if (bACF)
{
printf("\nCalculating current autocorrelation ... \n");
- sgk = calc_cacf(outf, prefactorav/PICO, cacf, time, nvfr, vfr, ie, nshift);
+ sgk = calc_cacf(outf, prefactorav / PICO, cacf, time, nvfr, vfr, ie, nshift);
if (ie > ii)
{
- snew(xfit, ie-ii+1);
- snew(yfit, ie-ii+1);
+ snew(xfit, ie - ii + 1);
+ snew(yfit, ie - ii + 1);
for (i = ii; i <= ie; i++)
{
- xfit[i-ii] = std::log(time[vfr[i]]);
- rtmp = std::abs(cacf[i]);
- yfit[i-ii] = std::log(rtmp);
-
+ xfit[i - ii] = std::log(time[vfr[i]]);
+ rtmp = std::abs(cacf[i]);
+ yfit[i - ii] = std::log(rtmp);
}
- lsq_y_ax_b(ie-ii, xfit, yfit, &sigma, &malt, &err, &rest);
+ lsq_y_ax_b(ie - ii, xfit, yfit, &sigma, &malt, &err, &rest);
malt = std::exp(malt);
sigma += 1.0;
- malt *= prefactorav*2.0e12/sigma;
+ malt *= prefactorav * 2.0e12 / sigma;
sfree(xfit);
sfree(yfit);
-
}
}
}
fprintf(stderr, "********************************************\n");
- dk_f = calceps(prefactor, md2-mdav2, mj2-mj, mjd-mjdav, eps_rf, FALSE);
+ dk_f = calceps(prefactor, md2 - mdav2, mj2 - mj, mjd - mjdav, eps_rf, FALSE);
fprintf(stderr, "\n\nFluctuations:\n epsilon=%f\n\n", dk_f);
- fprintf(stderr, "\n deltaM_D , deltaM_J, deltaM_JD: (%f, %f, %f)\n", md2-mdav2, mj2-mj, mjd-mjdav);
+ fprintf(stderr, "\n deltaM_D , deltaM_J, deltaM_JD: (%f, %f, %f)\n", md2 - mdav2, mj2 - mj,
+ mjd - mjdav);
fprintf(stderr, "\n********************************************\n");
if (bINT)
{
- dk_d = calceps(prefactor, md2-mdav2, mj2-mj, corint, eps_rf, TRUE);
+ dk_d = calceps(prefactor, md2 - mdav2, mj2 - mj, corint, eps_rf, TRUE);
fprintf(stderr, "\nStatic dielectric constant using integral and fluctuations: %f\n", dk_d);
- fprintf(stderr, "\n < M_JM_D > via integral: %.3f\n", -1.0*corint);
+ fprintf(stderr, "\n < M_JM_D > via integral: %.3f\n", -1.0 * corint);
}
fprintf(stderr, "\n***************************************************");
fprintf(stderr, "and corresponding refactor 1.0 / 3.0*V*k_B*T*EPSILON_0: %f \n", prefactor);
-
if (bACF && (ii < nvfr))
{
fprintf(stderr, "Integral and integrated fit to the current acf yields at t=%f:\n", time[vfr[ii]]);
- fprintf(stderr, "sigma=%8.3f (pure integral: %.3f)\n", sgk-malt*std::pow(time[vfr[ii]], sigma), sgk);
+ fprintf(stderr, "sigma=%8.3f (pure integral: %.3f)\n",
+ sgk - malt * std::pow(time[vfr[ii]], sigma), sgk);
}
if (ei > bi)
fprintf(stderr, "\nStart fit at %f ps (%f).\n", time[bi], bfit);
fprintf(stderr, "End fit at %f ps (%f).\n\n", time[ei], efit);
- snew(xfit, ei-bi+1);
- snew(yfit, ei-bi+1);
+ snew(xfit, ei - bi + 1);
+ snew(yfit, ei - bi + 1);
for (i = bi; i <= ei; i++)
{
- xfit[i-bi] = time[i];
- yfit[i-bi] = dsp2[i];
+ xfit[i - bi] = time[i];
+ yfit[i - bi] = dsp2[i];
}
- lsq_y_ax_b(ei-bi, xfit, yfit, &sigma, &malt, &err, &rest);
+ lsq_y_ax_b(ei - bi, xfit, yfit, &sigma, &malt, &err, &rest);
sigma *= 1e12;
- dk_d = calceps(prefactor, md2, 0.5*malt/prefactorav, corint, eps_rf, TRUE);
+ dk_d = calceps(prefactor, md2, 0.5 * malt / prefactorav, corint, eps_rf, TRUE);
- fprintf(stderr, "Einstein-Helfand fit to the MSD of the translational dipole moment yields:\n\n");
+ fprintf(stderr,
+ "Einstein-Helfand fit to the MSD of the translational dipole moment yields:\n\n");
fprintf(stderr, "sigma=%.4f\n", sigma);
- fprintf(stderr, "translational fraction of M^2: %.4f\n", 0.5*malt/prefactorav);
+ fprintf(stderr, "translational fraction of M^2: %.4f\n", 0.5 * malt / prefactorav);
fprintf(stderr, "Dielectric constant using EH: %.4f\n", dk_d);
sfree(xfit);
}
-
-int gmx_current(int argc, char *argv[])
+int gmx_current(int argc, char* argv[])
{
- static int nshift = 1000;
- static real temp = 300.0;
- static real eps_rf = 0.0;
- static gmx_bool bNoJump = TRUE;
- static real bfit = 100.0;
- static real bvit = 0.5;
- static real efit = 400.0;
- static real evit = 5.0;
- t_pargs pa[] = {
- { "-sh", FALSE, etINT, {&nshift},
- "Shift of the frames for averaging the correlation functions and the mean-square displacement."},
- { "-nojump", FALSE, etBOOL, {&bNoJump},
- "Removes jumps of atoms across the box."},
- { "-eps", FALSE, etREAL, {&eps_rf},
- "Dielectric constant of the surrounding medium. The value zero corresponds to infinity (tin-foil boundary conditions)."},
- { "-bfit", FALSE, etREAL, {&bfit},
- "Begin of the fit of the straight line to the MSD of the translational fraction of the dipole moment."},
- { "-efit", FALSE, etREAL, {&efit},
- "End of the fit of the straight line to the MSD of the translational fraction of the dipole moment."},
- { "-bvit", FALSE, etREAL, {&bvit},
- "Begin of the fit of the current autocorrelation function to a*t^b."},
- { "-evit", FALSE, etREAL, {&evit},
- "End of the fit of the current autocorrelation function to a*t^b."},
- { "-temp", FALSE, etREAL, {&temp},
- "Temperature for calculating epsilon."}
+ static int nshift = 1000;
+ static real temp = 300.0;
+ static real eps_rf = 0.0;
+ static gmx_bool bNoJump = TRUE;
+ static real bfit = 100.0;
+ static real bvit = 0.5;
+ static real efit = 400.0;
+ static real evit = 5.0;
+ t_pargs pa[] = {
+ { "-sh",
+ FALSE,
+ etINT,
+ { &nshift },
+ "Shift of the frames for averaging the correlation functions and the mean-square "
+ "displacement." },
+ { "-nojump", FALSE, etBOOL, { &bNoJump }, "Removes jumps of atoms across the box." },
+ { "-eps",
+ FALSE,
+ etREAL,
+ { &eps_rf },
+ "Dielectric constant of the surrounding medium. The value zero corresponds to "
+ "infinity (tin-foil boundary conditions)." },
+ { "-bfit",
+ FALSE,
+ etREAL,
+ { &bfit },
+ "Begin of the fit of the straight line to the MSD of the translational fraction "
+ "of the dipole moment." },
+ { "-efit",
+ FALSE,
+ etREAL,
+ { &efit },
+ "End of the fit of the straight line to the MSD of the translational fraction of "
+ "the dipole moment." },
+ { "-bvit",
+ FALSE,
+ etREAL,
+ { &bvit },
+ "Begin of the fit of the current autocorrelation function to a*t^b." },
+ { "-evit",
+ FALSE,
+ etREAL,
+ { &evit },
+ "End of the fit of the current autocorrelation function to a*t^b." },
+ { "-temp", FALSE, etREAL, { &temp }, "Temperature for calculating epsilon." }
};
- gmx_output_env_t *oenv;
- t_topology top;
- char **grpname = nullptr;
- const char *indexfn;
- t_trxframe fr;
- real *mass2 = nullptr;
- matrix box;
- int *index0;
- int *indexm = nullptr;
- int isize;
- t_trxstatus *status;
- int flags = 0;
- gmx_bool bACF;
- gmx_bool bINT;
- int ePBC = -1;
- int nmols;
- int i;
- real *qmol;
- FILE *outf = nullptr;
- FILE *mcor = nullptr;
- FILE *fmj = nullptr;
- FILE *fmd = nullptr;
- FILE *fmjdsp = nullptr;
- FILE *fcur = nullptr;
- t_filenm fnm[] = {
- { efTPS, nullptr, nullptr, ffREAD }, /* this is for the topology */
- { efNDX, nullptr, nullptr, ffOPTRD },
- { efTRX, "-f", nullptr, ffREAD }, /* and this for the trajectory */
- { efXVG, "-o", "current", ffWRITE },
- { efXVG, "-caf", "caf", ffOPTWR },
- { efXVG, "-dsp", "dsp", ffWRITE },
- { efXVG, "-md", "md", ffWRITE },
- { efXVG, "-mj", "mj", ffWRITE },
- { efXVG, "-mc", "mc", ffOPTWR }
- };
+ gmx_output_env_t* oenv;
+ t_topology top;
+ char** grpname = nullptr;
+ const char* indexfn;
+ t_trxframe fr;
+ real* mass2 = nullptr;
+ matrix box;
+ int* index0;
+ int* indexm = nullptr;
+ int isize;
+ t_trxstatus* status;
+ int flags = 0;
+ gmx_bool bACF;
+ gmx_bool bINT;
+ int ePBC = -1;
+ int nmols;
+ int i;
+ real* qmol;
+ FILE* outf = nullptr;
+ FILE* mcor = nullptr;
+ FILE* fmj = nullptr;
+ FILE* fmd = nullptr;
+ FILE* fmjdsp = nullptr;
+ FILE* fcur = nullptr;
+ t_filenm fnm[] = { { efTPS, nullptr, nullptr, ffREAD }, /* this is for the topology */
+ { efNDX, nullptr, nullptr, ffOPTRD },
+ { efTRX, "-f", nullptr, ffREAD }, /* and this for the trajectory */
+ { efXVG, "-o", "current", ffWRITE },
+ { efXVG, "-caf", "caf", ffOPTWR },
+ { efXVG, "-dsp", "dsp", ffWRITE },
+ { efXVG, "-md", "md", ffWRITE },
+ { efXVG, "-mj", "mj", ffWRITE },
+ { efXVG, "-mc", "mc", ffOPTWR } };
#define NFILE asize(fnm)
- const char *desc[] = {
- "[THISMODULE] is a tool for calculating the current autocorrelation function, the correlation",
+ const char* desc[] = {
+ "[THISMODULE] is a tool for calculating the current autocorrelation function, the "
+ "correlation",
"of the rotational and translational dipole moment of the system, and the resulting static",
"dielectric constant. To obtain a reasonable result, the index group has to be neutral.",
- "Furthermore, the routine is capable of extracting the static conductivity from the current ",
+ "Furthermore, the routine is capable of extracting the static conductivity from the "
+ "current ",
"autocorrelation function, if velocities are given. Additionally, an Einstein-Helfand fit ",
"can be used to obtain the static conductivity.",
"[PAR]",
- "The flag [TT]-caf[tt] is for the output of the current autocorrelation function and [TT]-mc[tt] writes the",
- "correlation of the rotational and translational part of the dipole moment in the corresponding",
+ "The flag [TT]-caf[tt] is for the output of the current autocorrelation function and "
+ "[TT]-mc[tt] writes the",
+ "correlation of the rotational and translational part of the dipole moment in the "
+ "corresponding",
"file. However, this option is only available for trajectories containing velocities.",
- "Options [TT]-sh[tt] and [TT]-tr[tt] are responsible for the averaging and integration of the",
+ "Options [TT]-sh[tt] and [TT]-tr[tt] are responsible for the averaging and integration of "
+ "the",
"autocorrelation functions. Since averaging proceeds by shifting the starting point",
- "through the trajectory, the shift can be modified with [TT]-sh[tt] to enable the choice of uncorrelated",
+ "through the trajectory, the shift can be modified with [TT]-sh[tt] to enable the choice "
+ "of uncorrelated",
"starting points. Towards the end, statistical inaccuracy grows and integrating the",
"correlation function only yields reliable values until a certain point, depending on",
- "the number of frames. The option [TT]-tr[tt] controls the region of the integral taken into account",
+ "the number of frames. The option [TT]-tr[tt] controls the region of the integral taken "
+ "into account",
"for calculating the static dielectric constant.",
"[PAR]",
- "Option [TT]-temp[tt] sets the temperature required for the computation of the static dielectric constant.",
+ "Option [TT]-temp[tt] sets the temperature required for the computation of the static "
+ "dielectric constant.",
"[PAR]",
- "Option [TT]-eps[tt] controls the dielectric constant of the surrounding medium for simulations using",
- "a Reaction Field or dipole corrections of the Ewald summation ([TT]-eps[tt]\\=0 corresponds to",
+ "Option [TT]-eps[tt] controls the dielectric constant of the surrounding medium for "
+ "simulations using",
+ "a Reaction Field or dipole corrections of the Ewald summation ([TT]-eps[tt]\\=0 "
+ "corresponds to",
"tin-foil boundary conditions).",
"[PAR]",
- "[TT]-[no]nojump[tt] unfolds the coordinates to allow free diffusion. This is required to get a continuous",
- "translational dipole moment, required for the Einstein-Helfand fit. The results from the fit allow",
- "the determination of the dielectric constant for system of charged molecules. However, it is also possible to extract",
- "the dielectric constant from the fluctuations of the total dipole moment in folded coordinates. But this",
- "option has to be used with care, since only very short time spans fulfill the approximation that the density",
- "of the molecules is approximately constant and the averages are already converged. To be on the safe side,",
- "the dielectric constant should be calculated with the help of the Einstein-Helfand method for",
+ "[TT]-[no]nojump[tt] unfolds the coordinates to allow free diffusion. This is required to "
+ "get a continuous",
+ "translational dipole moment, required for the Einstein-Helfand fit. The results from the "
+ "fit allow",
+ "the determination of the dielectric constant for system of charged molecules. However, it "
+ "is also possible to extract",
+ "the dielectric constant from the fluctuations of the total dipole moment in folded "
+ "coordinates. But this",
+ "option has to be used with care, since only very short time spans fulfill the "
+ "approximation that the density",
+ "of the molecules is approximately constant and the averages are already converged. To be "
+ "on the safe side,",
+ "the dielectric constant should be calculated with the help of the Einstein-Helfand method "
+ "for",
"the translational part of the dielectric constant."
};
/* At first the arguments will be parsed and the system information processed */
- if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW,
- NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
+ if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW, NFILE, fnm, asize(pa), pa,
+ asize(desc), desc, 0, nullptr, &oenv))
{
return 0;
}
{
if (bACF)
{
- outf = xvgropen(opt2fn("-caf", NFILE, fnm),
- "Current autocorrelation function", output_env_get_xvgr_tlabel(oenv),
- "ACF (e nm/ps)\\S2", oenv);
+ outf = xvgropen(opt2fn("-caf", NFILE, fnm), "Current autocorrelation function",
+ output_env_get_xvgr_tlabel(oenv), "ACF (e nm/ps)\\S2", oenv);
fprintf(outf, "# time\t acf\t average \t std.dev\n");
}
- fcur = xvgropen(opt2fn("-o", NFILE, fnm),
- "Current", output_env_get_xvgr_tlabel(oenv), "J(t) (e nm/ps)", oenv);
+ fcur = xvgropen(opt2fn("-o", NFILE, fnm), "Current", output_env_get_xvgr_tlabel(oenv),
+ "J(t) (e nm/ps)", oenv);
fprintf(fcur, "# time\t Jx\t Jy \t J_z \n");
if (bINT)
{
}
}
- fmj = xvgropen(opt2fn("-mj", NFILE, fnm),
- "Averaged translational part of M", output_env_get_xvgr_tlabel(oenv),
- "< M\\sJ\\N > (enm)", oenv);
+ fmj = xvgropen(opt2fn("-mj", NFILE, fnm), "Averaged translational part of M",
+ output_env_get_xvgr_tlabel(oenv), "< M\\sJ\\N > (enm)", oenv);
fprintf(fmj, "# time\t x\t y \t z \t average of M_J^2 \t std.dev\n");
- fmd = xvgropen(opt2fn("-md", NFILE, fnm),
- "Averaged rotational part of M", output_env_get_xvgr_tlabel(oenv),
- "< M\\sD\\N > (enm)", oenv);
+ fmd = xvgropen(opt2fn("-md", NFILE, fnm), "Averaged rotational part of M",
+ output_env_get_xvgr_tlabel(oenv), "< M\\sD\\N > (enm)", oenv);
fprintf(fmd, "# time\t x\t y \t z \t average of M_D^2 \t std.dev\n");
- fmjdsp = xvgropen(opt2fn("-dsp", NFILE, fnm),
- "MSD of the squared translational dipole moment M",
- output_env_get_xvgr_tlabel(oenv),
- "<|M\\sJ\\N(t)-M\\sJ\\N(0)|\\S2\\N > / 6.0*V*k\\sB\\N*T / Sm\\S-1\\Nps\\S-1\\N",
- oenv);
+ fmjdsp = xvgropen(
+ opt2fn("-dsp", NFILE, fnm), "MSD of the squared translational dipole moment M",
+ output_env_get_xvgr_tlabel(oenv),
+ "<|M\\sJ\\N(t)-M\\sJ\\N(0)|\\S2\\N > / 6.0*V*k\\sB\\N*T / Sm\\S-1\\Nps\\S-1\\N", oenv);
/* System information is read and prepared, dielectric() processes the frames
* and calculates the requested quantities */
- dielectric(fmj, fmd, outf, fcur, mcor, fmjdsp, bNoJump, bACF, bINT, ePBC, top, fr,
- temp, bfit, efit, bvit, evit, status, isize, nmols, nshift,
- index0, indexm, mass2, qmol, eps_rf, oenv);
+ dielectric(fmj, fmd, outf, fcur, mcor, fmjdsp, bNoJump, bACF, bINT, ePBC, top, fr, temp, bfit, efit,
+ bvit, evit, status, isize, nmols, nshift, index0, indexm, mass2, qmol, eps_rf, oenv);
xvgrclose(fmj);
xvgrclose(fmd);
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff