-/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
+/*
+ * This file is part of the GROMACS molecular simulation package.
*
- *
- * This source code is part of
- *
- * G R O M A C S
- *
- * GROningen MAchine for Chemical Simulations
- *
- * VERSION 3.2.0
- * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2004, The GROMACS development team,
- * check out http://www.gromacs.org for more information.
-
- * This program is free software; you can redistribute it and/or
- * modify it under the terms of the GNU General Public License
- * as published by the Free Software Foundation; either version 2
+ * Copyright (c) 2001-2004, The GROMACS development team.
+ * Copyright (c) 2013,2014, 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.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
* of the License, or (at your option) any later version.
*
- * If you want to redistribute modifications, please consider that
- * scientific software is very special. Version control is crucial -
- * bugs must be traceable. We will be happy to consider code for
- * inclusion in the official distribution, but derived work must not
- * be called official GROMACS. Details are found in the README & COPYING
- * files - if they are missing, get the official version at www.gromacs.org.
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
*
- * To help us fund GROMACS development, we humbly ask that you cite
- * the papers on the package - you can find them in the top README file.
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
*
- * For more info, check our website at http://www.gromacs.org
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
*
- * And Hey:
- * Green Red Orange Magenta Azure Cyan Skyblue
+ * 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 <string.h>
#include <math.h>
-#include <ctype.h>
+#include <stdlib.h>
+#include <string.h>
-#include "typedefs.h"
-#include "gmx_fatal.h"
-#include "vec.h"
-#include "string2.h"
-#include "smalloc.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 "xvgr.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 "physics.h"
+#include "gromacs/math/units.h"
#include "gromacs/fileio/tpxio.h"
#include "gromacs/fileio/trxio.h"
-#include "viewit.h"
-#include "mtop_util.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;
} enerdat_t;
typedef struct {
- gmx_large_int_t nsteps;
- gmx_large_int_t npoints;
+ gmx_int64_t nsteps;
+ gmx_int64_t npoints;
int nframes;
int *step;
int *steps;
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)
{
fprintf(vout, "%10d %10.5e\n", j, mypow(violaver[j]/nframes, minthird));
}
#endif
- ffclose(vout);
+ gmx_ffclose(vout);
fprintf(stdout, "\nSum of violations averaged over simulation: %g nm\n",
sumt);
}
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);
}
fprintf(fp1, "\n");
}
- ffclose(fp0);
- ffclose(fp1);
+ gmx_ffclose(fp0);
+ gmx_ffclose(fp1);
}
typedef struct {
- gmx_large_int_t np;
+ gmx_int64_t np;
double sum;
double sav;
double sav2;
typedef struct {
int b;
ee_sum_t sum;
- gmx_large_int_t nst;
- gmx_large_int_t nst_min;
+ gmx_int64_t nst;
+ gmx_int64_t nst_min;
} ener_ee_t;
static void clear_ee_sum(ee_sum_t *ees)
{
int nb, i, f, nee;
double sum, sum2, sump, see2;
- gmx_large_int_t steps, np, p, bound_nb;
+ gmx_int64_t steps, np, p, bound_nb;
enerdat_t *ed;
exactsum_t *es;
gmx_bool bAllZero;
int nbmin, int nbmax)
{
int i, j;
- double vvhh, vv, v, h, hh2, vv2, varv, hh, varh, tt, cv, cp, alpha, kappa, dcp, et, varet;
+ double vv, v, h, varv, hh, varh, tt, cv, cp, alpha, kappa, dcp, et, varet;
double NANO3;
enum {
eVol, eEnth, eTemp, eEtot, eNR
fprintf(fp,"Found %s data.\n",my_ener[i]);
*/ }
/* Compute it all! */
- vvhh = alpha = kappa = cp = dcp = cv = NOTSET;
+ alpha = kappa = cp = dcp = cv = NOTSET;
/* Temperature */
tt = NOTSET;
/* Alpha, dcp */
if ((ii[eVol] < nset) && (ii[eEnth] < nset) && (ii[eTemp] < nset))
{
- vvhh = 0;
+ double v_sum, h_sum, vh_sum, v_aver, h_aver, vh_aver;
+ vh_sum = v_sum = h_sum = 0;
for (j = 0; (j < edat->nframes); j++)
{
- v = edat->s[ii[eVol]].ener[j]*NANO3;
- h = KILO*edat->s[ii[eEnth]].ener[j]/AVOGADRO;
- vvhh += (v*h);
+ v = edat->s[ii[eVol]].ener[j]*NANO3;
+ h = KILO*edat->s[ii[eEnth]].ener[j]/AVOGADRO;
+ v_sum += v;
+ h_sum += h;
+ vh_sum += (v*h);
}
- vvhh /= edat->nframes;
- alpha = (vvhh-vv*hh)/(vv*BOLTZMANN*tt*tt);
- dcp = (vv*AVOGADRO/nmol)*tt*sqr(alpha)/(kappa);
+ vh_aver = vh_sum / edat->nframes;
+ v_aver = v_sum / edat->nframes;
+ h_aver = h_sum / edat->nframes;
+ alpha = (vh_aver-v_aver*h_aver)/(v_aver*BOLTZMANN*tt*tt);
+ dcp = (v_aver*AVOGADRO/nmol)*tt*sqr(alpha)/(kappa);
}
if (tt != NOTSET)
{
fprintf(fp, "varv = %10g (m^6)\n", varv*AVOGADRO/nmol);
}
- if (vvhh != NOTSET)
- {
- fprintf(fp, "vvhh = %10g (m^3 J)\n", vvhh);
- }
}
if (vv != NOTSET)
{
static void analyse_ener(gmx_bool bCorr, const char *corrfn,
gmx_bool bFee, gmx_bool bSum, gmx_bool bFluct,
gmx_bool bVisco, const char *visfn, int nmol,
- gmx_large_int_t start_step, double start_t,
- gmx_large_int_t step, double t,
- double time[], real reftemp,
+ gmx_int64_t start_step, double start_t,
+ gmx_int64_t step, double t,
+ real reftemp,
enerdata_t *edat,
int nset, int set[], gmx_bool *bIsEner,
char **leg, gmx_enxnm_t *enm,
real xxx, integral, intBulk, Temp = 0, Pres = 0;
real sfrac, oldfrac, diffsum, diffav, fstep, pr_aver, pr_stddev, pr_errest;
double beta = 0, expE, expEtot, *fee = NULL;
- gmx_large_int_t nsteps;
+ gmx_int64_t nsteps;
int nexact, nnotexact;
double x1m, x1mk;
int i, j, k, nout;
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 */
intBulk += 0.5*(eneset[11][i-1] + eneset[11][i])*factor;
fprintf(fp, "%10g %10g %10g\n", (i*Dt), integral, intBulk);
}
- ffclose(fp);
+ gmx_ffclose(fp);
+
+ for (i = 0; i < 12; i++)
+ {
+ sfree(eneset[i]);
+ }
+ sfree(eneset);
}
else if (bCorr)
{
}
if (fp)
{
- ffclose(fp);
+ gmx_ffclose(fp);
}
sfree(fr);
}
/* write the data */
if (nblock_hist > 0)
{
- gmx_large_int_t sum = 0;
+ gmx_int64_t sum = 0;
/* histograms */
for (i = 0; i < fr->nblock; i++)
{
if (blk->id == enxDHHIST)
{
double foreign_lambda, dx;
- gmx_large_int_t x0;
+ gmx_int64_t x0;
int nhist, derivative;
/* check the block types etc. */
if ( (blk->nsub < 2) ||
(blk->sub[0].type != xdr_datatype_double) ||
- (blk->sub[1].type != xdr_datatype_large_int) ||
+ (blk->sub[1].type != xdr_datatype_int64) ||
(blk->sub[0].nr < 2) ||
(blk->sub[1].nr < 2) )
{
int cur = 0;
#define NEXT (1-cur)
int nre, teller, teller_disre, nfr;
- gmx_large_int_t start_step;
+ gmx_int64_t start_step;
int nor = 0, nex = 0, norfr = 0, enx_i = 0;
real start_t;
real *bounds = NULL, *violaver = NULL, *oobs = NULL, *orient = NULL, *odrms = NULL;
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
{
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);
}
}
}
close_enx(fp);
if (out)
{
- ffclose(out);
+ gmx_ffclose(out);
}
if (bDRAll)
{
- ffclose(fp_pairs);
+ gmx_ffclose(fp_pairs);
}
if (bORT)
{
- ffclose(fort);
+ gmx_ffclose(fort);
}
if (bODT)
{
- ffclose(fodt);
+ gmx_ffclose(fodt);
}
if (bORA)
{
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);
}
{
fprintf(out, "%5d %g\n", or_label[i], orient[i]/norfr);
}
- ffclose(out);
+ gmx_ffclose(out);
}
if (bODA)
{
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);
}
{
fprintf(out, "%5d %g\n", or_label[i], orient[i]/norfr-oobs[i]);
}
- ffclose(out);
+ gmx_ffclose(out);
}
if (bODR)
{
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);
}
{
fprintf(out, "%5d %g\n", or_label[i], sqrt(odrms[i]/norfr));
}
- ffclose(out);
+ gmx_ffclose(out);
}
if (bOTEN)
{
- ffclose(foten);
+ gmx_ffclose(foten);
}
if (bDisRe)
{
if (fp_dhdl)
{
- ffclose(fp_dhdl);
+ gmx_ffclose(fp_dhdl);
printf("\n\nWrote %d lambda values with %d samples as ",
dh_lambdas, dh_samples);
if (dh_hists > 0)
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