#include "gromacs/utility/pleasecite.h"
#include "gromacs/utility/smalloc.h"
-static void dump_dih_trr(int nframes, int nangles, real **dih, const char *fn,
- real *time)
+static void dump_dih_trr(int nframes, int nangles, real** dih, const char* fn, real* time)
{
int i, j, k, l, m, na;
- struct t_fileio *fio;
- rvec *x;
- matrix box = {{2, 0, 0}, {0, 2, 0}, {0, 0, 2}};
+ struct t_fileio* fio;
+ rvec* x;
+ matrix box = { { 2, 0, 0 }, { 0, 2, 0 }, { 0, 0, 2 } };
- na = (nangles*2);
+ na = (nangles * 2);
if ((na % 3) != 0)
{
- na = 1+na/3;
+ na = 1 + na / 3;
}
else
{
- na = na/3;
+ na = na / 3;
}
- printf("There are %d dihedrals. Will fill %d atom positions with cos/sin\n",
- nangles, na);
+ printf("There are %d dihedrals. Will fill %d atom positions with cos/sin\n", nangles, na);
snew(x, na);
fio = gmx_trr_open(fn, "w");
for (i = 0; (i < nframes); i++)
sfree(x);
}
-int gmx_g_angle(int argc, char *argv[])
+int gmx_g_angle(int argc, char* argv[])
{
- static const char *desc[] = {
+ static const char* desc[] = {
"[THISMODULE] computes the angle distribution for a number of angles",
"or dihedrals.[PAR]",
"With option [TT]-ov[tt], you can plot the average angle of",
"records a histogram of the times between such transitions,",
"assuming the input trajectory frames are equally spaced in time."
};
- static const char *opt[] = { nullptr, "angle", "dihedral", "improper", "ryckaert-bellemans", nullptr };
- static gmx_bool bALL = FALSE, bChandler = FALSE, bAverCorr = FALSE, bPBC = TRUE;
+ static const char* opt[] = { nullptr, "angle", "dihedral", "improper", "ryckaert-bellemans",
+ nullptr };
+ static gmx_bool bALL = FALSE, bChandler = FALSE, bAverCorr = FALSE, bPBC = TRUE;
static real binwidth = 1;
t_pargs pa[] = {
- { "-type", FALSE, etENUM, {opt},
- "Type of angle to analyse" },
- { "-all", FALSE, etBOOL, {&bALL},
- "Plot all angles separately in the averages file, in the order of appearance in the index file." },
- { "-binwidth", FALSE, etREAL, {&binwidth},
+ { "-type", FALSE, etENUM, { opt }, "Type of angle to analyse" },
+ { "-all",
+ FALSE,
+ etBOOL,
+ { &bALL },
+ "Plot all angles separately in the averages file, in the order of appearance in the "
+ "index file." },
+ { "-binwidth",
+ FALSE,
+ etREAL,
+ { &binwidth },
"binwidth (degrees) for calculating the distribution" },
- { "-periodic", FALSE, etBOOL, {&bPBC},
- "Print dihedral angles modulo 360 degrees" },
- { "-chandler", FALSE, etBOOL, {&bChandler},
- "Use Chandler correlation function (N[trans] = 1, N[gauche] = 0) rather than cosine correlation function. Trans is defined as phi < -60 or phi > 60." },
- { "-avercorr", FALSE, etBOOL, {&bAverCorr},
+ { "-periodic", FALSE, etBOOL, { &bPBC }, "Print dihedral angles modulo 360 degrees" },
+ { "-chandler",
+ FALSE,
+ etBOOL,
+ { &bChandler },
+ "Use Chandler correlation function (N[trans] = 1, N[gauche] = 0) rather than cosine "
+ "correlation function. Trans is defined as phi < -60 or phi > 60." },
+ { "-avercorr",
+ FALSE,
+ etBOOL,
+ { &bAverCorr },
"Average the correlation functions for the individual angles/dihedrals" }
};
- static const char *bugs[] = {
+ static const char* bugs[] = {
"Counting transitions only works for dihedrals with multiplicity 3"
};
- FILE *out;
- real dt;
- int isize;
- int *index;
- char *grpname;
- real maxang, S2, norm_fac, maxstat;
- unsigned long mode;
- int nframes, maxangstat, mult, *angstat;
- int i, j, nangles, first, last;
- gmx_bool bAver, bRb, bPeriodic,
- bFrac, /* calculate fraction too? */
- bTrans, /* worry about transtions too? */
- bCorr; /* correlation function ? */
- double tfrac = 0;
- char title[256];
- real **dih = nullptr; /* mega array with all dih. angles at all times*/
- real *time, *trans_frac, *aver_angle;
- t_filenm fnm[] = {
- { efTRX, "-f", nullptr, ffREAD },
- { efNDX, nullptr, "angle", ffREAD },
- { efXVG, "-od", "angdist", ffWRITE },
- { efXVG, "-ov", "angaver", ffOPTWR },
- { efXVG, "-of", "dihfrac", ffOPTWR },
- { efXVG, "-ot", "dihtrans", ffOPTWR },
- { efXVG, "-oh", "trhisto", ffOPTWR },
- { efXVG, "-oc", "dihcorr", ffOPTWR },
- { efTRR, "-or", nullptr, ffOPTWR }
- };
+ FILE* out;
+ real dt;
+ int isize;
+ int* index;
+ char* grpname;
+ real maxang, S2, norm_fac, maxstat;
+ unsigned long mode;
+ int nframes, maxangstat, mult, *angstat;
+ int i, j, nangles, first, last;
+ gmx_bool bAver, bRb, bPeriodic, bFrac, /* calculate fraction too? */
+ bTrans, /* worry about transtions too? */
+ bCorr; /* correlation function ? */
+ double tfrac = 0;
+ char title[256];
+ real** dih = nullptr; /* mega array with all dih. angles at all times*/
+ real * time, *trans_frac, *aver_angle;
+ t_filenm fnm[] = { { efTRX, "-f", nullptr, ffREAD }, { efNDX, nullptr, "angle", ffREAD },
+ { efXVG, "-od", "angdist", ffWRITE }, { efXVG, "-ov", "angaver", ffOPTWR },
+ { efXVG, "-of", "dihfrac", ffOPTWR }, { efXVG, "-ot", "dihtrans", ffOPTWR },
+ { efXVG, "-oh", "trhisto", ffOPTWR }, { efXVG, "-oc", "dihcorr", ffOPTWR },
+ { efTRR, "-or", nullptr, ffOPTWR } };
#define NFILE asize(fnm)
- int npargs;
- t_pargs *ppa;
- gmx_output_env_t *oenv;
+ int npargs;
+ t_pargs* ppa;
+ gmx_output_env_t* oenv;
npargs = asize(pa);
ppa = add_acf_pargs(&npargs, pa);
- if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME,
- NFILE, fnm, npargs, ppa, asize(desc), desc, asize(bugs), bugs,
- &oenv))
+ if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME, NFILE, fnm, npargs, ppa,
+ asize(desc), desc, asize(bugs), bugs, &oenv))
{
sfree(ppa);
return 0;
maxang = 360.0;
bRb = FALSE;
- GMX_RELEASE_ASSERT(opt[0] != nullptr, "Internal option inconsistency; opt[0]==NULL after processing");
+ GMX_RELEASE_ASSERT(opt[0] != nullptr,
+ "Internal option inconsistency; opt[0]==NULL after processing");
switch (opt[0][0])
{
mult = 3;
maxang = 180.0;
break;
- case 'd':
- break;
- case 'i':
- break;
- case 'r':
- bRb = TRUE;
- break;
+ case 'd': break;
+ case 'i': break;
+ case 'r': bRb = TRUE; break;
}
if (opt2bSet("-or", NFILE, fnm))
}
/* Calculate bin size */
- maxangstat = gmx::roundToInt(maxang/binwidth);
- binwidth = maxang/maxangstat;
+ maxangstat = gmx::roundToInt(maxang / binwidth);
+ binwidth = maxang / maxangstat;
rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
- nangles = isize/mult;
+ nangles = isize / mult;
if ((isize % mult) != 0)
{
- gmx_fatal(FARGS, "number of index elements not multiple of %d, "
+ gmx_fatal(FARGS,
+ "number of index elements not multiple of %d, "
"these can not be %s\n",
mult, (mult == 3) ? "angle triplets" : "dihedral quadruplets");
}
if (bFrac && !bRb)
{
- fprintf(stderr, "Warning:"
+ fprintf(stderr,
+ "Warning:"
" calculating fractions as defined in this program\n"
"makes sense for Ryckaert Bellemans dihs. only. Ignoring -of\n\n");
bFrac = FALSE;
}
- if ( (bTrans || bFrac || bCorr) && mult == 3)
+ if ((bTrans || bFrac || bCorr) && mult == 3)
{
- gmx_fatal(FARGS, "Can only do transition, fraction or correlation\n"
+ gmx_fatal(FARGS,
+ "Can only do transition, fraction or correlation\n"
"on dihedrals. Select -d\n");
}
* We need to know the nr of frames so we can allocate memory for an array
* with all dihedral angles at all timesteps. Works for me.
*/
- if (bTrans || bCorr || bALL || opt2bSet("-or", NFILE, fnm))
+ if (bTrans || bCorr || bALL || opt2bSet("-or", NFILE, fnm))
{
snew(dih, nangles);
}
snew(angstat, maxangstat);
read_ang_dih(ftp2fn(efTRX, NFILE, fnm), (mult == 3),
- bALL || bCorr || bTrans || opt2bSet("-or", NFILE, fnm),
- bRb, bPBC, maxangstat, angstat,
- &nframes, &time, isize, index, &trans_frac, &aver_angle, dih,
- oenv);
+ bALL || bCorr || bTrans || opt2bSet("-or", NFILE, fnm), bRb, bPBC, maxangstat,
+ angstat, &nframes, &time, isize, index, &trans_frac, &aver_angle, dih, oenv);
- dt = (time[nframes-1]-time[0])/(nframes-1);
+ dt = (time[nframes - 1] - time[0]) / (nframes - 1);
if (bAver)
{
sprintf(title, "Average Angle: %s", grpname);
- out = xvgropen(opt2fn("-ov", NFILE, fnm),
- title, "Time (ps)", "Angle (degrees)", oenv);
+ out = xvgropen(opt2fn("-ov", NFILE, fnm), title, "Time (ps)", "Angle (degrees)", oenv);
for (i = 0; (i < nframes); i++)
{
- fprintf(out, "%10.5f %8.3f", time[i], aver_angle[i]*RAD2DEG);
+ fprintf(out, "%10.5f %8.3f", time[i], aver_angle[i] * RAD2DEG);
if (bALL)
{
for (j = 0; (j < nangles); j++)
if (bPBC)
{
real dd = dih[j][i];
- fprintf(out, " %8.3f", std::atan2(std::sin(dd), std::cos(dd))*RAD2DEG);
+ fprintf(out, " %8.3f", std::atan2(std::sin(dd), std::cos(dd)) * RAD2DEG);
}
else
{
- fprintf(out, " %8.3f", dih[j][i]*RAD2DEG);
+ fprintf(out, " %8.3f", dih[j][i] * RAD2DEG);
}
}
}
if (bFrac)
{
sprintf(title, "Trans fraction: %s", grpname);
- out = xvgropen(opt2fn("-of", NFILE, fnm),
- title, "Time (ps)", "Fraction", oenv);
+ out = xvgropen(opt2fn("-of", NFILE, fnm), title, "Time (ps)", "Fraction", oenv);
tfrac = 0.0;
for (i = 0; (i < nframes); i++)
{
if (bTrans)
{
- ana_dih_trans(opt2fn("-ot", NFILE, fnm), opt2fn("-oh", NFILE, fnm),
- dih, nframes, nangles, grpname, time, bRb, oenv);
+ ana_dih_trans(opt2fn("-ot", NFILE, fnm), opt2fn("-oh", NFILE, fnm), dih, nframes, nangles,
+ grpname, time, bRb, oenv);
}
if (bCorr)
if (bChandler)
{
- real dval, sixty = DEG2RAD*60;
+ real dval, sixty = DEG2RAD * 60;
gmx_bool bTest;
for (i = 0; (i < nangles); i++)
}
if (bTest)
{
- dih[i][j] = dval-tfrac;
+ dih[i][j] = dval - tfrac;
}
else
{
{
mode = eacCos;
}
- do_autocorr(opt2fn("-oc", NFILE, fnm), oenv,
- "Dihedral Autocorrelation Function",
+ do_autocorr(opt2fn("-oc", NFILE, fnm), oenv, "Dihedral Autocorrelation Function",
nframes, nangles, dih, dt, mode, bAverCorr);
}
}
/* Determine the non-zero part of the distribution */
- for (first = 0; (first < maxangstat-1) && (angstat[first+1] == 0); first++)
- {
- ;
- }
- for (last = maxangstat-1; (last > 0) && (angstat[last-1] == 0); last--)
- {
- ;
- }
+ for (first = 0; (first < maxangstat - 1) && (angstat[first + 1] == 0); first++) {}
+ for (last = maxangstat - 1; (last > 0) && (angstat[last - 1] == 0); last--) {}
- double aver = 0;
- printf("Found points in the range from %d to %d (max %d)\n",
- first, last, maxangstat);
- if (bTrans || bCorr || bALL || opt2bSet("-or", NFILE, fnm))
- { /* It's better to re-calculate Std. Dev per sample */
+ double aver = 0;
+ printf("Found points in the range from %d to %d (max %d)\n", first, last, maxangstat);
+ if (bTrans || bCorr || bALL || opt2bSet("-or", NFILE, fnm))
+ { /* It's better to re-calculate Std. Dev per sample */
real b_aver = aver_angle[0];
real b = dih[0][0];
real delta;
for (int i = 0; (i < nframes); i++)
{
- delta = correctRadianAngleRange(aver_angle[i] - b_aver);
+ delta = correctRadianAngleRange(aver_angle[i] - b_aver);
b_aver += delta;
- aver += b_aver;
+ aver += b_aver;
for (int j = 0; (j < nangles); j++)
{
- delta = correctRadianAngleRange(dih[j][i] - b);
- b += delta;
+ delta = correctRadianAngleRange(dih[j][i] - b);
+ b += delta;
}
}
}
else
- { /* Incorrect for Std. Dev. */
+ { /* Incorrect for Std. Dev. */
real delta, b_aver = aver_angle[0];
for (i = 0; (i < nframes); i++)
{
- delta = correctRadianAngleRange(aver_angle[i] - b_aver);
+ delta = correctRadianAngleRange(aver_angle[i] - b_aver);
b_aver += delta;
- aver += b_aver;
+ aver += b_aver;
}
}
- aver /= nframes;
+ aver /= nframes;
double aversig = correctRadianAngleRange(aver);
aversig *= RAD2DEG;
- aver *= RAD2DEG;
+ aver *= RAD2DEG;
printf(" < angle > = %g\n", aversig);
if (mult == 3)
fprintf(stderr, "Order parameter S^2 = %g\n", S2);
}
- bPeriodic = (mult == 4) && (first == 0) && (last == maxangstat-1);
+ bPeriodic = (mult == 4) && (first == 0) && (last == maxangstat - 1);
out = xvgropen(opt2fn("-od", NFILE, fnm), title, "Degrees", "", oenv);
if (output_env_get_print_xvgr_codes(oenv))
{
fprintf(out, "@ subtitle \"average angle: %g\\So\\N\"\n", aver);
}
- norm_fac = 1.0/(nangles*nframes*binwidth);
+ norm_fac = 1.0 / (nangles * nframes * binwidth);
if (bPeriodic)
{
maxstat = 0;
for (i = first; (i <= last); i++)
{
- maxstat = std::max(maxstat, angstat[i]*norm_fac);
+ maxstat = std::max(maxstat, angstat[i] * norm_fac);
}
if (output_env_get_print_xvgr_codes(oenv))
{
fprintf(out, "@ world xmin -180\n");
fprintf(out, "@ world xmax 180\n");
fprintf(out, "@ world ymin 0\n");
- fprintf(out, "@ world ymax %g\n", maxstat*1.05);
+ fprintf(out, "@ world ymax %g\n", maxstat * 1.05);
fprintf(out, "@ xaxis tick major 60\n");
fprintf(out, "@ xaxis tick minor 30\n");
fprintf(out, "@ yaxis tick major 0.005\n");
}
for (i = first; (i <= last); i++)
{
- fprintf(out, "%10g %10f\n", i*binwidth+180.0-maxang, angstat[i]*norm_fac);
+ fprintf(out, "%10g %10f\n", i * binwidth + 180.0 - maxang, angstat[i] * norm_fac);
}
if (bPeriodic)
{
/* print first bin again as last one */
- fprintf(out, "%10g %10f\n", 180.0, angstat[0]*norm_fac);
+ fprintf(out, "%10g %10f\n", 180.0, angstat[0] * norm_fac);
}
xvgrclose(out);
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff