* To help us fund GROMACS development, we humbly ask that you cite
* the research papers on the package. Check out http://www.gromacs.org.
*/
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
+#include "gmxpre.h"
#include <math.h>
#include <stdlib.h>
#include <string.h>
-#include "typedefs.h"
+#include "gromacs/legacyheaders/typedefs.h"
#include "gromacs/utility/fatalerror.h"
#include "gromacs/math/vec.h"
#include "gromacs/utility/cstringutil.h"
#include "gromacs/utility/smalloc.h"
#include "gromacs/fileio/enxio.h"
#include "gromacs/commandline/pargs.h"
-#include "names.h"
-#include "copyrite.h"
-#include "macros.h"
+#include "gromacs/legacyheaders/names.h"
+#include "gromacs/legacyheaders/copyrite.h"
+#include "gromacs/legacyheaders/macros.h"
#include "gromacs/fileio/xvgr.h"
#include "gstat.h"
#include "gromacs/math/units.h"
#include "gromacs/fileio/tpxio.h"
#include "gromacs/fileio/trxio.h"
-#include "viewit.h"
+#include "gromacs/legacyheaders/viewit.h"
#include "gromacs/topology/mtop_util.h"
#include "gmx_ana.h"
-#include "mdebin.h"
+#include "gromacs/legacyheaders/mdebin.h"
static real minthird = -1.0/3.0, minsixth = -1.0/6.0;
int *set;
gmx_bool bVerbose = TRUE;
- if ((getenv("VERBOSE")) != NULL)
+ if ((getenv("GMX_ENER_VERBOSE")) != NULL)
{
bVerbose = FALSE;
}
const char *fm2 = "%3d %-34s";
char **newnm = NULL;
- if ((getenv("VERBOSE")) != NULL)
+ if ((getenv("GMX_ENER_VERBOSE")) != NULL)
{
bVerbose = FALSE;
}
j = 0;
for (k = 0; k < nre; k++)
{
- newnm[k] = strdup(nm[k].name);
+ newnm[k] = gmx_strdup(nm[k].name);
/* Insert dashes in all the names */
while ((ptr = strchr(newnm[k], ' ')) != NULL)
{
}
static void einstein_visco(const char *fn, const char *fni, int nsets,
- int nframes, real **sum,
- real V, real T, int nsteps, double time[],
+ int nint, real **eneint,
+ real V, real T, double dt,
const output_env_t oenv)
{
FILE *fp0, *fp1;
double av[4], avold[4];
- double fac, dt, di;
+ double fac, di;
int i, j, m, nf4;
- if (nframes < 1)
- {
- return;
- }
-
- dt = (time[1]-time[0]);
- nf4 = nframes/4+1;
+ nf4 = nint/4 + 1;
for (i = 0; i <= nsets; i++)
{
"Time (ps)", "(kg m\\S-1\\N s\\S-1\\N ps)", oenv);
fp1 = xvgropen(fn, "Shear viscosity using Einstein relation",
"Time (ps)", "(kg m\\S-1\\N s\\S-1\\N)", oenv);
- for (i = 1; i < nf4; i++)
+ for (i = 0; i < nf4; i++)
{
- fac = dt*nframes/nsteps;
for (m = 0; m <= nsets; m++)
{
av[m] = 0;
}
- for (j = 0; j < nframes-i; j++)
+ for (j = 0; j < nint - i; j++)
{
for (m = 0; m < nsets; m++)
{
- di = sqr(fac*(sum[m][j+i]-sum[m][j]));
+ di = sqr(eneint[m][j+i] - eneint[m][j]);
av[m] += di;
av[nsets] += di/nsets;
}
}
/* Convert to SI for the viscosity */
- fac = (V*NANO*NANO*NANO*PICO*1e10)/(2*BOLTZMANN*T)/(nframes-i);
- fprintf(fp0, "%10g", time[i]-time[0]);
+ fac = (V*NANO*NANO*NANO*PICO*1e10)/(2*BOLTZMANN*T)/(nint - i);
+ fprintf(fp0, "%10g", i*dt);
for (m = 0; (m <= nsets); m++)
{
av[m] = fac*av[m];
fprintf(fp0, " %10g", av[m]);
}
fprintf(fp0, "\n");
- fprintf(fp1, "%10g", 0.5*(time[i]+time[i-1])-time[0]);
+ fprintf(fp1, "%10g", (i + 0.5)*dt);
for (m = 0; (m <= nsets); m++)
{
fprintf(fp1, " %10g", (av[m]-avold[m])/dt);
gmx_bool bVisco, const char *visfn, int nmol,
gmx_int64_t start_step, double start_t,
gmx_int64_t step, double t,
- double time[], real reftemp,
+ real reftemp,
enerdata_t *edat,
int nset, int set[], gmx_bool *bIsEner,
char **leg, gmx_enxnm_t *enm,
const char* leg[] = { "Shear", "Bulk" };
real factor;
real **eneset;
- real **enesum;
+ real **eneint;
/* Assume pressure tensor is in Pxx Pxy Pxz Pyx Pyy Pyz Pzx Pzy Pzz */
{
snew(eneset[i], edat->nframes);
}
- snew(enesum, 3);
- for (i = 0; i < 3; i++)
- {
- snew(enesum[i], edat->nframes);
- }
for (i = 0; (i < edat->nframes); i++)
{
eneset[0][i] = 0.5*(edat->s[1].ener[i]+edat->s[3].ener[i]);
eneset[j][i] = edat->s[j].ener[i];
}
eneset[11][i] -= Pres;
- enesum[0][i] = 0.5*(edat->s[1].es[i].sum+edat->s[3].es[i].sum);
- enesum[1][i] = 0.5*(edat->s[2].es[i].sum+edat->s[6].es[i].sum);
- enesum[2][i] = 0.5*(edat->s[5].es[i].sum+edat->s[7].es[i].sum);
+ }
+
+ /* Determine integrals of the off-diagonal pressure elements */
+ snew(eneint, 3);
+ for (i = 0; i < 3; i++)
+ {
+ snew(eneint[i], edat->nframes + 1);
+ }
+ eneint[0][0] = 0;
+ eneint[1][0] = 0;
+ eneint[2][0] = 0;
+ for (i = 0; i < edat->nframes; i++)
+ {
+ eneint[0][i+1] = eneint[0][i] + 0.5*(edat->s[1].es[i].sum + edat->s[3].es[i].sum)*Dt/edat->points[i];
+ eneint[1][i+1] = eneint[1][i] + 0.5*(edat->s[2].es[i].sum + edat->s[6].es[i].sum)*Dt/edat->points[i];
+ eneint[2][i+1] = eneint[2][i] + 0.5*(edat->s[5].es[i].sum + edat->s[7].es[i].sum)*Dt/edat->points[i];
}
einstein_visco("evisco.xvg", "eviscoi.xvg",
- 3, edat->nframes, enesum, Vaver, Temp, nsteps, time, oenv);
+ 3, edat->nframes+1, eneint, Vaver, Temp, Dt, oenv);
+
+ for (i = 0; i < 3; i++)
+ {
+ sfree(eneint[i]);
+ }
+ sfree(eneint);
/*do_autocorr(corrfn,buf,nenergy,3,eneset,Dt,eacNormal,TRUE);*/
/* Do it for shear viscosity */
fprintf(fp, "%10g %10g %10g\n", (i*Dt), integral, intBulk);
}
gmx_ffclose(fp);
+
+ for (i = 0; i < 12; i++)
+ {
+ sfree(eneset[i]);
+ }
+ sfree(eneset);
}
else if (bCorr)
{
t_filenm fnm[] = {
{ efEDR, "-f", NULL, ffREAD },
{ efEDR, "-f2", NULL, ffOPTRD },
- { efTPX, "-s", NULL, ffOPTRD },
+ { efTPR, "-s", NULL, ffOPTRD },
{ efXVG, "-o", "energy", ffWRITE },
{ efXVG, "-viol", "violaver", ffOPTWR },
{ efXVG, "-pairs", "pairs", ffOPTWR },
npargs = asize(pa);
ppa = add_acf_pargs(&npargs, pa);
if (!parse_common_args(&argc, argv,
- PCA_CAN_VIEW | PCA_CAN_BEGIN | PCA_CAN_END | PCA_BE_NICE,
+ PCA_CAN_VIEW | PCA_CAN_BEGIN | PCA_CAN_END,
NFILE, fnm, npargs, ppa, asize(desc), desc, 0, NULL, &oenv))
{
return 0;
}
if (bSum)
{
- leg[nset] = strdup("Sum");
+ leg[nset] = gmx_strdup("Sum");
xvgr_legend(out, nset+1, (const char**)leg, oenv);
}
else
if (bORIRE || bOTEN)
{
- get_orires_parms(ftp2fn(efTPX, NFILE, fnm), &nor, &nex, &or_label, &oobs);
+ get_orires_parms(ftp2fn(efTPR, NFILE, fnm), &nor, &nex, &or_label, &oobs);
}
if (bORIRE)
{
fort = xvgropen(opt2fn("-ort", NFILE, fnm), "Calculated orientations",
"Time (ps)", "", oenv);
- if (bOrinst)
+ if (bOrinst && output_env_get_print_xvgr_codes(oenv))
{
fprintf(fort, "%s", orinst_sub);
}
fodt = xvgropen(opt2fn("-odt", NFILE, fnm),
"Orientation restraint deviation",
"Time (ps)", "", oenv);
- if (bOrinst)
+ if (bOrinst && output_env_get_print_xvgr_codes(oenv))
{
fprintf(fodt, "%s", orinst_sub);
}
for (j = 0; j < 3; j++)
{
sprintf(buf, "eig%d", j+1);
- otenleg[(bOvec ? 12 : 3)*i+j] = strdup(buf);
+ otenleg[(bOvec ? 12 : 3)*i+j] = gmx_strdup(buf);
}
if (bOvec)
{
for (j = 0; j < 9; j++)
{
sprintf(buf, "vec%d%s", j/3+1, j%3 == 0 ? "x" : (j%3 == 1 ? "y" : "z"));
- otenleg[12*i+3+j] = strdup(buf);
+ otenleg[12*i+3+j] = gmx_strdup(buf);
}
}
}
}
else if (bDisRe)
{
- nbounds = get_bounds(ftp2fn(efTPX, NFILE, fnm), &bounds, &index, &pair, &npairs,
+ nbounds = get_bounds(ftp2fn(efTPR, NFILE, fnm), &bounds, &index, &pair, &npairs,
&mtop, &top, &ir);
snew(violaver, npairs);
out = xvgropen(opt2fn("-o", NFILE, fnm), "Sum of Violations",
}
else if (bDHDL)
{
- get_dhdl_parms(ftp2fn(efTPX, NFILE, fnm), &ir);
+ get_dhdl_parms(ftp2fn(efTPR, NFILE, fnm), &ir);
}
/* Initiate energies and set them to zero */
out = xvgropen(opt2fn("-ora", NFILE, fnm),
"Average calculated orientations",
"Restraint label", "", oenv);
- if (bOrinst)
+ if (bOrinst && output_env_get_print_xvgr_codes(oenv))
{
fprintf(out, "%s", orinst_sub);
}
out = xvgropen(opt2fn("-oda", NFILE, fnm),
"Average restraint deviation",
"Restraint label", "", oenv);
- if (bOrinst)
+ if (bOrinst && output_env_get_print_xvgr_codes(oenv))
{
fprintf(out, "%s", orinst_sub);
}
out = xvgropen(opt2fn("-odr", NFILE, fnm),
"RMS orientation restraint deviations",
"Restraint label", "", oenv);
- if (bOrinst)
+ if (bOrinst && output_env_get_print_xvgr_codes(oenv))
{
fprintf(out, "%s", orinst_sub);
}
bVisco, opt2fn("-vis", NFILE, fnm),
nmol,
start_step, start_t, frame[cur].step, frame[cur].t,
- time, reftemp, &edat,
+ reftemp, &edat,
nset, set, bIsEner, leg, enm, Vaver, ezero, nbmin, nbmax,
oenv);
if (bFluctProps)
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff