From: Mark Abraham Date: Wed, 24 Jun 2015 21:33:24 +0000 (+0200) Subject: Merge release-4-6 into release-5-0 X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=368bdf6154423ed6e9ece5bf9866d28727b3d0b8;p=alexxy%2Fgromacs.git Merge release-4-6 into release-5-0 Change-Id: I6cc179f767cc4e0fc1dc3fc0b85428df986c16df --- 368bdf6154423ed6e9ece5bf9866d28727b3d0b8 diff --cc src/gromacs/gmxana/gmx_energy.c index de70d3f402,0000000000..68df53d576 mode 100644,000000..100644 --- a/src/gromacs/gmxana/gmx_energy.c +++ b/src/gromacs/gmxana/gmx_energy.c @@@ -1,2845 -1,0 +1,2849 @@@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright (c) 1991-2000, University of Groningen, The Netherlands. + * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2013,2014, by the GROMACS development team, led by ++ * Copyright (c) 2013,2014,2015, 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. + * + * 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. + * + * 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. + * + * 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. + * + * 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 +#endif + +#include +#include + +#include "typedefs.h" +#include "gmx_fatal.h" +#include "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 "gstat.h" +#include "physics.h" +#include "gromacs/fileio/tpxio.h" +#include "gromacs/fileio/trxio.h" +#include "viewit.h" +#include "mtop_util.h" +#include "gmx_ana.h" +#include "mdebin.h" + +static real minthird = -1.0/3.0, minsixth = -1.0/6.0; + +typedef struct { + real sum; + real sum2; +} exactsum_t; + +typedef struct { + real *ener; + exactsum_t *es; + gmx_bool bExactStat; + double av; + double rmsd; + double ee; + double slope; +} enerdat_t; + +typedef struct { + gmx_int64_t nsteps; + gmx_int64_t npoints; + int nframes; + int *step; + int *steps; + int *points; + enerdat_t *s; ++ gmx_bool bHaveSums; +} enerdata_t; + +static double mypow(double x, double y) +{ + if (x > 0) + { + return pow(x, y); + } + else + { + return 0.0; + } +} + +static int *select_it(int nre, char *nm[], int *nset) +{ + gmx_bool *bE; + int n, k, j, i; + int *set; + gmx_bool bVerbose = TRUE; + + if ((getenv("GMX_ENER_VERBOSE")) != NULL) + { + bVerbose = FALSE; + } + + fprintf(stderr, "Select the terms you want from the following list\n"); + fprintf(stderr, "End your selection with 0\n"); + + if (bVerbose) + { + for (k = 0; (k < nre); ) + { + for (j = 0; (j < 4) && (k < nre); j++, k++) + { + fprintf(stderr, " %3d=%14s", k+1, nm[k]); + } + fprintf(stderr, "\n"); + } + } + + snew(bE, nre); + do + { + if (1 != scanf("%d", &n)) + { + gmx_fatal(FARGS, "Error reading user input"); + } + if ((n > 0) && (n <= nre)) + { + bE[n-1] = TRUE; + } + } + while (n != 0); + + snew(set, nre); + for (i = (*nset) = 0; (i < nre); i++) + { + if (bE[i]) + { + set[(*nset)++] = i; + } + } + + sfree(bE); + + return set; +} + +static void chomp(char *buf) +{ + int len = strlen(buf); + + while ((len > 0) && (buf[len-1] == '\n')) + { + buf[len-1] = '\0'; + len--; + } +} + +static int *select_by_name(int nre, gmx_enxnm_t *nm, int *nset) +{ + gmx_bool *bE; + int n, k, kk, j, i, nmatch, nind, nss; + int *set; + gmx_bool bEOF, bVerbose = TRUE, bLong = FALSE; + char *ptr, buf[STRLEN]; + const char *fm4 = "%3d %-14s"; + const char *fm2 = "%3d %-34s"; + char **newnm = NULL; + + if ((getenv("GMX_ENER_VERBOSE")) != NULL) + { + bVerbose = FALSE; + } + + fprintf(stderr, "\n"); + fprintf(stderr, "Select the terms you want from the following list by\n"); + fprintf(stderr, "selecting either (part of) the name or the number or a combination.\n"); + fprintf(stderr, "End your selection with an empty line or a zero.\n"); + fprintf(stderr, "-------------------------------------------------------------------\n"); + + snew(newnm, nre); + j = 0; + for (k = 0; k < nre; k++) + { + newnm[k] = strdup(nm[k].name); + /* Insert dashes in all the names */ + while ((ptr = strchr(newnm[k], ' ')) != NULL) + { + *ptr = '-'; + } + if (bVerbose) + { + if (j == 0) + { + if (k > 0) + { + fprintf(stderr, "\n"); + } + bLong = FALSE; + for (kk = k; kk < k+4; kk++) + { + if (kk < nre && strlen(nm[kk].name) > 14) + { + bLong = TRUE; + } + } + } + else + { + fprintf(stderr, " "); + } + if (!bLong) + { + fprintf(stderr, fm4, k+1, newnm[k]); + j++; + if (j == 4) + { + j = 0; + } + } + else + { + fprintf(stderr, fm2, k+1, newnm[k]); + j++; + if (j == 2) + { + j = 0; + } + } + } + } + if (bVerbose) + { + fprintf(stderr, "\n\n"); + } + + snew(bE, nre); + + bEOF = FALSE; + while (!bEOF && (fgets2(buf, STRLEN-1, stdin))) + { + /* Remove newlines */ + chomp(buf); + + /* Remove spaces */ + trim(buf); + + /* Empty line means end of input */ + bEOF = (strlen(buf) == 0); + if (!bEOF) + { + ptr = buf; + do + { + if (!bEOF) + { + /* First try to read an integer */ + nss = sscanf(ptr, "%d", &nind); + if (nss == 1) + { + /* Zero means end of input */ + if (nind == 0) + { + bEOF = TRUE; + } + else if ((1 <= nind) && (nind <= nre)) + { + bE[nind-1] = TRUE; + } + else + { + fprintf(stderr, "number %d is out of range\n", nind); + } + } + else + { + /* Now try to read a string */ + i = strlen(ptr); + nmatch = 0; + for (nind = 0; nind < nre; nind++) + { + if (gmx_strcasecmp(newnm[nind], ptr) == 0) + { + bE[nind] = TRUE; + nmatch++; + } + } + if (nmatch == 0) + { + i = strlen(ptr); + nmatch = 0; + for (nind = 0; nind < nre; nind++) + { + if (gmx_strncasecmp(newnm[nind], ptr, i) == 0) + { + bE[nind] = TRUE; + nmatch++; + } + } + if (nmatch == 0) + { + fprintf(stderr, "String '%s' does not match anything\n", ptr); + } + } + } + } + /* Look for the first space, and remove spaces from there */ + if ((ptr = strchr(ptr, ' ')) != NULL) + { + trim(ptr); + } + } + while (!bEOF && (ptr && (strlen(ptr) > 0))); + } + } + + snew(set, nre); + for (i = (*nset) = 0; (i < nre); i++) + { + if (bE[i]) + { + set[(*nset)++] = i; + } + } + + sfree(bE); + + if (*nset == 0) + { + gmx_fatal(FARGS, "No energy terms selected"); + } + + for (i = 0; (i < nre); i++) + { + sfree(newnm[i]); + } + sfree(newnm); + + return set; +} + +static void get_dhdl_parms(const char *topnm, t_inputrec *ir) +{ + gmx_mtop_t mtop; + int natoms; + t_iatom *iatom; + matrix box; + + /* all we need is the ir to be able to write the label */ + read_tpx(topnm, ir, box, &natoms, NULL, NULL, NULL, &mtop); +} + +static void get_orires_parms(const char *topnm, + int *nor, int *nex, int **label, real **obs) +{ + gmx_mtop_t mtop; + gmx_localtop_t *top; + t_inputrec ir; + t_iparams *ip; + int natoms, i; + t_iatom *iatom; + int nb; + matrix box; + + read_tpx(topnm, &ir, box, &natoms, NULL, NULL, NULL, &mtop); + top = gmx_mtop_generate_local_top(&mtop, &ir); + + ip = top->idef.iparams; + iatom = top->idef.il[F_ORIRES].iatoms; + + /* Count how many distance restraint there are... */ + nb = top->idef.il[F_ORIRES].nr; + if (nb == 0) + { + gmx_fatal(FARGS, "No orientation restraints in topology!\n"); + } + + *nor = nb/3; + *nex = 0; + snew(*label, *nor); + snew(*obs, *nor); + for (i = 0; i < nb; i += 3) + { + (*label)[i/3] = ip[iatom[i]].orires.label; + (*obs)[i/3] = ip[iatom[i]].orires.obs; + if (ip[iatom[i]].orires.ex >= *nex) + { + *nex = ip[iatom[i]].orires.ex+1; + } + } + fprintf(stderr, "Found %d orientation restraints with %d experiments", + *nor, *nex); +} + +static int get_bounds(const char *topnm, + real **bounds, int **index, int **dr_pair, int *npairs, + gmx_mtop_t *mtop, gmx_localtop_t **ltop, t_inputrec *ir) +{ + gmx_localtop_t *top; + t_functype *functype; + t_iparams *ip; + int natoms, i, j, k, type, ftype, natom; + t_ilist *disres; + t_iatom *iatom; + real *b; + int *ind, *pair; + int nb, label1; + matrix box; + + read_tpx(topnm, ir, box, &natoms, NULL, NULL, NULL, mtop); + snew(*ltop, 1); + top = gmx_mtop_generate_local_top(mtop, ir); + *ltop = top; + + functype = top->idef.functype; + ip = top->idef.iparams; + + /* Count how many distance restraint there are... */ + nb = top->idef.il[F_DISRES].nr; + if (nb == 0) + { + gmx_fatal(FARGS, "No distance restraints in topology!\n"); + } + + /* Allocate memory */ + snew(b, nb); + snew(ind, nb); + snew(pair, nb+1); + + /* Fill the bound array */ + nb = 0; + for (i = 0; (i < top->idef.ntypes); i++) + { + ftype = functype[i]; + if (ftype == F_DISRES) + { + + label1 = ip[i].disres.label; + b[nb] = ip[i].disres.up1; + ind[nb] = label1; + nb++; + } + } + *bounds = b; + + /* Fill the index array */ + label1 = -1; + disres = &(top->idef.il[F_DISRES]); + iatom = disres->iatoms; + for (i = j = k = 0; (i < disres->nr); ) + { + type = iatom[i]; + ftype = top->idef.functype[type]; + natom = interaction_function[ftype].nratoms+1; + if (label1 != top->idef.iparams[type].disres.label) + { + pair[j] = k; + label1 = top->idef.iparams[type].disres.label; + j++; + } + k++; + i += natom; + } + pair[j] = k; + *npairs = k; + if (j != nb) + { + gmx_incons("get_bounds for distance restraints"); + } + + *index = ind; + *dr_pair = pair; + + return nb; +} + +static void calc_violations(real rt[], real rav3[], int nb, int index[], + real bounds[], real *viol, double *st, double *sa) +{ + const real sixth = 1.0/6.0; + int i, j; + double rsum, rav, sumaver, sumt; + + sumaver = 0; + sumt = 0; + for (i = 0; (i < nb); i++) + { + rsum = 0.0; + rav = 0.0; + for (j = index[i]; (j < index[i+1]); j++) + { + if (viol) + { + viol[j] += mypow(rt[j], -3.0); + } + rav += sqr(rav3[j]); + rsum += mypow(rt[j], -6); + } + rsum = max(0.0, mypow(rsum, -sixth)-bounds[i]); + rav = max(0.0, mypow(rav, -sixth)-bounds[i]); + + sumt += rsum; + sumaver += rav; + } + *st = sumt; + *sa = sumaver; +} + +static void analyse_disre(const char *voutfn, int nframes, + real violaver[], real bounds[], int index[], + int pair[], int nbounds, + const output_env_t oenv) +{ + FILE *vout; + double sum, sumt, sumaver; + int i, j; + + /* Subtract bounds from distances, to calculate violations */ + calc_violations(violaver, violaver, + nbounds, pair, bounds, NULL, &sumt, &sumaver); + +#ifdef DEBUG + fprintf(stdout, "\nSum of violations averaged over simulation: %g nm\n", + sumaver); + fprintf(stdout, "Largest violation averaged over simulation: %g nm\n\n", + sumt); +#endif + vout = xvgropen(voutfn, "r\\S-3\\N average violations", "DR Index", "nm", + oenv); + sum = 0.0; + sumt = 0.0; + for (i = 0; (i < nbounds); i++) + { + /* Do ensemble averaging */ + sumaver = 0; + for (j = pair[i]; (j < pair[i+1]); j++) + { + sumaver += sqr(violaver[j]/nframes); + } + sumaver = max(0.0, mypow(sumaver, minsixth)-bounds[i]); + + sumt += sumaver; + sum = max(sum, sumaver); + fprintf(vout, "%10d %10.5e\n", index[i], sumaver); + } +#ifdef DEBUG + for (j = 0; (j < dr.ndr); j++) + { + fprintf(vout, "%10d %10.5e\n", j, mypow(violaver[j]/nframes, minthird)); + } +#endif + gmx_ffclose(vout); + + fprintf(stdout, "\nSum of violations averaged over simulation: %g nm\n", + sumt); + fprintf(stdout, "Largest violation averaged over simulation: %g nm\n\n", sum); + + do_view(oenv, voutfn, "-graphtype bar"); +} + +static void einstein_visco(const char *fn, const char *fni, int nsets, + 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, di; + int i, j, m, nf4; + + nf4 = nint/4 + 1; + + for (i = 0; i <= nsets; i++) + { + avold[i] = 0; + } + fp0 = xvgropen(fni, "Shear viscosity integral", + "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 = 0; i < nf4; i++) + { + for (m = 0; m <= nsets; m++) + { + av[m] = 0; + } + for (j = 0; j < nint - i; j++) + { + for (m = 0; m < nsets; m++) + { + 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)/(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", (i + 0.5)*dt); + for (m = 0; (m <= nsets); m++) + { + fprintf(fp1, " %10g", (av[m]-avold[m])/dt); + avold[m] = av[m]; + } + fprintf(fp1, "\n"); + } + gmx_ffclose(fp0); + gmx_ffclose(fp1); +} + +typedef struct { + gmx_int64_t np; + double sum; + double sav; + double sav2; +} ee_sum_t; + +typedef struct { + int b; + ee_sum_t sum; + gmx_int64_t nst; + gmx_int64_t nst_min; +} ener_ee_t; + +static void clear_ee_sum(ee_sum_t *ees) +{ + ees->sav = 0; + ees->sav2 = 0; + ees->np = 0; + ees->sum = 0; +} + +static void add_ee_sum(ee_sum_t *ees, double sum, int np) +{ + ees->np += np; + ees->sum += sum; +} + +static void add_ee_av(ee_sum_t *ees) +{ + double av; + + av = ees->sum/ees->np; + ees->sav += av; + ees->sav2 += av*av; + ees->np = 0; + ees->sum = 0; +} + +static double calc_ee2(int nb, ee_sum_t *ees) +{ + return (ees->sav2/nb - dsqr(ees->sav/nb))/(nb - 1); +} + +static void set_ee_av(ener_ee_t *eee) +{ + if (debug) + { + char buf[STEPSTRSIZE]; + fprintf(debug, "Storing average for err.est.: %s steps\n", + gmx_step_str(eee->nst, buf)); + } + add_ee_av(&eee->sum); + eee->b++; + if (eee->b == 1 || eee->nst < eee->nst_min) + { + eee->nst_min = eee->nst; + } + eee->nst = 0; +} + +static void calc_averages(int nset, enerdata_t *edat, int nbmin, int nbmax) +{ + int nb, i, f, nee; + double sum, sum2, sump, see2; + gmx_int64_t steps, np, p, bound_nb; + enerdat_t *ed; + exactsum_t *es; + gmx_bool bAllZero; + double x, sx, sy, sxx, sxy; + ener_ee_t *eee; + + /* Check if we have exact statistics over all points */ + for (i = 0; i < nset; i++) + { + ed = &edat->s[i]; + ed->bExactStat = FALSE; - if (edat->npoints > 0) ++ if (edat->bHaveSums) + { + /* All energy file sum entries 0 signals no exact sums. + * But if all energy values are 0, we still have exact sums. + */ + bAllZero = TRUE; + for (f = 0; f < edat->nframes && !ed->bExactStat; f++) + { + if (ed->ener[i] != 0) + { + bAllZero = FALSE; + } + ed->bExactStat = (ed->es[f].sum != 0); + } + if (bAllZero) + { + ed->bExactStat = TRUE; + } + } + } + + snew(eee, nbmax+1); + for (i = 0; i < nset; i++) + { + ed = &edat->s[i]; + + sum = 0; + sum2 = 0; + np = 0; + sx = 0; + sy = 0; + sxx = 0; + sxy = 0; + for (nb = nbmin; nb <= nbmax; nb++) + { + eee[nb].b = 0; + clear_ee_sum(&eee[nb].sum); + eee[nb].nst = 0; + eee[nb].nst_min = 0; + } + for (f = 0; f < edat->nframes; f++) + { + es = &ed->es[f]; + + if (ed->bExactStat) + { + /* Add the sum and the sum of variances to the totals. */ + p = edat->points[f]; + sump = es->sum; + sum2 += es->sum2; + if (np > 0) + { + sum2 += dsqr(sum/np - (sum + es->sum)/(np + p)) + *np*(np + p)/p; + } + } + else + { + /* Add a single value to the sum and sum of squares. */ + p = 1; + sump = ed->ener[f]; + sum2 += dsqr(sump); + } + + /* sum has to be increased after sum2 */ + np += p; + sum += sump; + + /* For the linear regression use variance 1/p. + * Note that sump is the sum, not the average, so we don't need p*. + */ + x = edat->step[f] - 0.5*(edat->steps[f] - 1); + sx += p*x; + sy += sump; + sxx += p*x*x; + sxy += x*sump; + + for (nb = nbmin; nb <= nbmax; nb++) + { + /* Check if the current end step is closer to the desired + * block boundary than the next end step. + */ + bound_nb = (edat->step[0]-1)*nb + edat->nsteps*(eee[nb].b+1); + if (eee[nb].nst > 0 && + bound_nb - edat->step[f-1]*nb < edat->step[f]*nb - bound_nb) + { + set_ee_av(&eee[nb]); + } + if (f == 0) + { + eee[nb].nst = 1; + } + else + { + eee[nb].nst += edat->step[f] - edat->step[f-1]; + } + if (ed->bExactStat) + { + add_ee_sum(&eee[nb].sum, es->sum, edat->points[f]); + } + else + { + add_ee_sum(&eee[nb].sum, edat->s[i].ener[f], 1); + } + bound_nb = (edat->step[0]-1)*nb + edat->nsteps*(eee[nb].b+1); + if (edat->step[f]*nb >= bound_nb) + { + set_ee_av(&eee[nb]); + } + } + } + + edat->s[i].av = sum/np; + if (ed->bExactStat) + { + edat->s[i].rmsd = sqrt(sum2/np); + } + else + { + edat->s[i].rmsd = sqrt(sum2/np - dsqr(edat->s[i].av)); + } + + if (edat->nframes > 1) + { + edat->s[i].slope = (np*sxy - sx*sy)/(np*sxx - sx*sx); + } + else + { + edat->s[i].slope = 0; + } + + nee = 0; + see2 = 0; + for (nb = nbmin; nb <= nbmax; nb++) + { + /* Check if we actually got nb blocks and if the smallest + * block is not shorter than 80% of the average. + */ + if (debug) + { + char buf1[STEPSTRSIZE], buf2[STEPSTRSIZE]; + fprintf(debug, "Requested %d blocks, we have %d blocks, min %s nsteps %s\n", + nb, eee[nb].b, + gmx_step_str(eee[nb].nst_min, buf1), + gmx_step_str(edat->nsteps, buf2)); + } + if (eee[nb].b == nb && 5*nb*eee[nb].nst_min >= 4*edat->nsteps) + { + see2 += calc_ee2(nb, &eee[nb].sum); + nee++; + } + } + if (nee > 0) + { + edat->s[i].ee = sqrt(see2/nee); + } + else + { + edat->s[i].ee = -1; + } + } + sfree(eee); +} + +static enerdata_t *calc_sum(int nset, enerdata_t *edat, int nbmin, int nbmax) +{ + enerdata_t *esum; + enerdat_t *s; + int f, i; + double sum; + + snew(esum, 1); + *esum = *edat; + snew(esum->s, 1); + s = &esum->s[0]; + snew(s->ener, esum->nframes); + snew(s->es, esum->nframes); + + s->bExactStat = TRUE; + s->slope = 0; + for (i = 0; i < nset; i++) + { + if (!edat->s[i].bExactStat) + { + s->bExactStat = FALSE; + } + s->slope += edat->s[i].slope; + } + + for (f = 0; f < edat->nframes; f++) + { + sum = 0; + for (i = 0; i < nset; i++) + { + sum += edat->s[i].ener[f]; + } + s->ener[f] = sum; + sum = 0; + for (i = 0; i < nset; i++) + { + sum += edat->s[i].es[f].sum; + } + s->es[f].sum = sum; + s->es[f].sum2 = 0; + } + + calc_averages(1, esum, nbmin, nbmax); + + return esum; +} + +static char *ee_pr(double ee, char *buf) +{ + char tmp[100]; + double rnd; + + if (ee < 0) + { + sprintf(buf, "%s", "--"); + } + else + { + /* Round to two decimals by printing. */ + sprintf(tmp, "%.1e", ee); + sscanf(tmp, "%lf", &rnd); + sprintf(buf, "%g", rnd); + } + + return buf; +} + +static void remove_drift(int nset, int nbmin, int nbmax, real dt, enerdata_t *edat) +{ +/* Remove the drift by performing a fit to y = ax+b. + Uses 5 iterations. */ + int i, j, k; + double delta, d, sd, sd2; + + edat->npoints = edat->nframes; + edat->nsteps = edat->nframes; + + for (k = 0; (k < 5); k++) + { + for (i = 0; (i < nset); i++) + { + delta = edat->s[i].slope*dt; + + if (NULL != debug) + { + fprintf(debug, "slope for set %d is %g\n", i, edat->s[i].slope); + } + + for (j = 0; (j < edat->nframes); j++) + { + edat->s[i].ener[j] -= j*delta; + edat->s[i].es[j].sum = 0; + edat->s[i].es[j].sum2 = 0; + } + } + calc_averages(nset, edat, nbmin, nbmax); + } +} + +static void calc_fluctuation_props(FILE *fp, + gmx_bool bDriftCorr, real dt, + int nset, int nmol, + char **leg, enerdata_t *edat, + int nbmin, int nbmax) +{ + int i, j; + double vv, v, h, varv, hh, varh, tt, cv, cp, alpha, kappa, dcp, et, varet; + double NANO3; + enum { + eVol, eEnth, eTemp, eEtot, eNR + }; + const char *my_ener[] = { "Volume", "Enthalpy", "Temperature", "Total Energy" }; + int ii[eNR]; + + NANO3 = NANO*NANO*NANO; + if (!bDriftCorr) + { + fprintf(fp, "\nYou may want to use the -driftcorr flag in order to correct\n" + "for spurious drift in the graphs. Note that this is not\n" + "a substitute for proper equilibration and sampling!\n"); + } + else + { + remove_drift(nset, nbmin, nbmax, dt, edat); + } + for (i = 0; (i < eNR); i++) + { + for (ii[i] = 0; (ii[i] < nset && + (gmx_strcasecmp(leg[ii[i]], my_ener[i]) != 0)); ii[i]++) + { + ; + } +/* if (ii[i] < nset) + fprintf(fp,"Found %s data.\n",my_ener[i]); + */ } + /* Compute it all! */ + alpha = kappa = cp = dcp = cv = NOTSET; + + /* Temperature */ + tt = NOTSET; + if (ii[eTemp] < nset) + { + tt = edat->s[ii[eTemp]].av; + } + /* Volume */ + vv = varv = NOTSET; + if ((ii[eVol] < nset) && (ii[eTemp] < nset)) + { + vv = edat->s[ii[eVol]].av*NANO3; + varv = dsqr(edat->s[ii[eVol]].rmsd*NANO3); + kappa = (varv/vv)/(BOLTZMANN*tt); + } + /* Enthalpy */ + hh = varh = NOTSET; + if ((ii[eEnth] < nset) && (ii[eTemp] < nset)) + { + hh = KILO*edat->s[ii[eEnth]].av/AVOGADRO; + varh = dsqr(KILO*edat->s[ii[eEnth]].rmsd/AVOGADRO); + cp = AVOGADRO*((varh/nmol)/(BOLTZMANN*tt*tt)); + } + /* Total energy */ + et = varet = NOTSET; + if ((ii[eEtot] < nset) && (hh == NOTSET) && (tt != NOTSET)) + { + /* Only compute cv in constant volume runs, which we can test + by checking whether the enthalpy was computed. + */ + et = edat->s[ii[eEtot]].av; + varet = sqr(edat->s[ii[eEtot]].rmsd); + cv = KILO*((varet/nmol)/(BOLTZ*tt*tt)); + } + /* Alpha, dcp */ + if ((ii[eVol] < nset) && (ii[eEnth] < nset) && (ii[eTemp] < nset)) + { + 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; + v_sum += v; + h_sum += h; + vh_sum += (v*h); + } + 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) + { + if (nmol < 2) + { + fprintf(fp, "\nWARNING: nmol = %d, this may not be what you want.\n", + nmol); + } + fprintf(fp, "\nTemperature dependent fluctuation properties at T = %g.\n", tt); + fprintf(fp, "\nHeat capacities obtained from fluctuations do *not* include\n"); + fprintf(fp, "quantum corrections. If you want to get a more accurate estimate\n"); + fprintf(fp, "please use the g_dos program.\n\n"); + fprintf(fp, "WARNING: Please verify that your simulations are converged and perform\n" + "a block-averaging error analysis (not implemented in g_energy yet)\n"); + + if (debug != NULL) + { + if (varv != NOTSET) + { + fprintf(fp, "varv = %10g (m^6)\n", varv*AVOGADRO/nmol); + } + } + if (vv != NOTSET) + { + fprintf(fp, "Volume = %10g m^3/mol\n", + vv*AVOGADRO/nmol); + } + if (varh != NOTSET) + { + fprintf(fp, "Enthalpy = %10g kJ/mol\n", + hh*AVOGADRO/(KILO*nmol)); + } + if (alpha != NOTSET) + { + fprintf(fp, "Coefficient of Thermal Expansion Alpha_P = %10g (1/K)\n", + alpha); + } + if (kappa != NOTSET) + { + fprintf(fp, "Isothermal Compressibility Kappa = %10g (J/m^3)\n", + kappa); + fprintf(fp, "Adiabatic bulk modulus = %10g (m^3/J)\n", + 1.0/kappa); + } + if (cp != NOTSET) + { + fprintf(fp, "Heat capacity at constant pressure Cp = %10g J/mol K\n", + cp); + } + if (cv != NOTSET) + { + fprintf(fp, "Heat capacity at constant volume Cv = %10g J/mol K\n", + cv); + } + if (dcp != NOTSET) + { + fprintf(fp, "Cp-Cv = %10g J/mol K\n", + dcp); + } + please_cite(fp, "Allen1987a"); + } + else + { + fprintf(fp, "You should select the temperature in order to obtain fluctuation properties.\n"); + } +} + +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_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 Vaver, real ezero, + int nbmin, int nbmax, + const output_env_t oenv) +{ + FILE *fp; + /* Check out the printed manual for equations! */ + double Dt, aver, stddev, errest, delta_t, totaldrift; + enerdata_t *esum = NULL; + 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_int64_t nsteps; + int nexact, nnotexact; + double x1m, x1mk; + int i, j, k, nout; + real chi2; + char buf[256], eebuf[100]; + + nsteps = step - start_step + 1; + if (nsteps < 1) + { + fprintf(stdout, "Not enough steps (%s) for statistics\n", + gmx_step_str(nsteps, buf)); + } + else + { + /* Calculate the time difference */ + delta_t = t - start_t; + + fprintf(stdout, "\nStatistics over %s steps [ %.4f through %.4f ps ], %d data sets\n", + gmx_step_str(nsteps, buf), start_t, t, nset); + + calc_averages(nset, edat, nbmin, nbmax); + + if (bSum) + { + esum = calc_sum(nset, edat, nbmin, nbmax); + } + - if (edat->npoints == 0) ++ if (!edat->bHaveSums) + { + nexact = 0; + nnotexact = nset; + } + else + { + nexact = 0; + nnotexact = 0; + for (i = 0; (i < nset); i++) + { + if (edat->s[i].bExactStat) + { + nexact++; + } + else + { + nnotexact++; + } + } + } + + if (nnotexact == 0) + { + fprintf(stdout, "All statistics are over %s points\n", + gmx_step_str(edat->npoints, buf)); + } + else if (nexact == 0 || edat->npoints == edat->nframes) + { + fprintf(stdout, "All statistics are over %d points (frames)\n", + edat->nframes); + } + else + { + fprintf(stdout, "The term%s", nnotexact == 1 ? "" : "s"); + for (i = 0; (i < nset); i++) + { + if (!edat->s[i].bExactStat) + { + fprintf(stdout, " '%s'", leg[i]); + } + } + fprintf(stdout, " %s has statistics over %d points (frames)\n", + nnotexact == 1 ? "is" : "are", edat->nframes); + fprintf(stdout, "All other statistics are over %s points\n", + gmx_step_str(edat->npoints, buf)); + } + fprintf(stdout, "\n"); + + fprintf(stdout, "%-24s %10s %10s %10s %10s", + "Energy", "Average", "Err.Est.", "RMSD", "Tot-Drift"); + if (bFee) + { + fprintf(stdout, " %10s\n", "-kT ln"); + } + else + { + fprintf(stdout, "\n"); + } + fprintf(stdout, "-------------------------------------------------------------------------------\n"); + + /* Initiate locals, only used with -sum */ + expEtot = 0; + if (bFee) + { + beta = 1.0/(BOLTZ*reftemp); + snew(fee, nset); + } + for (i = 0; (i < nset); i++) + { + aver = edat->s[i].av; + stddev = edat->s[i].rmsd; + errest = edat->s[i].ee; + + if (bFee) + { + expE = 0; + for (j = 0; (j < edat->nframes); j++) + { + expE += exp(beta*(edat->s[i].ener[j] - aver)/nmol); + } + if (bSum) + { + expEtot += expE/edat->nframes; + } + + fee[i] = log(expE/edat->nframes)/beta + aver/nmol; + } + if (strstr(leg[i], "empera") != NULL) + { + Temp = aver; + } + else if (strstr(leg[i], "olum") != NULL) + { + Vaver = aver; + } + else if (strstr(leg[i], "essure") != NULL) + { + Pres = aver; + } + if (bIsEner[i]) + { + pr_aver = aver/nmol-ezero; + pr_stddev = stddev/nmol; + pr_errest = errest/nmol; + } + else + { + pr_aver = aver; + pr_stddev = stddev; + pr_errest = errest; + } + + /* Multiply the slope in steps with the number of steps taken */ + totaldrift = (edat->nsteps - 1)*edat->s[i].slope; + if (bIsEner[i]) + { + totaldrift /= nmol; + } + + fprintf(stdout, "%-24s %10g %10s %10g %10g", + leg[i], pr_aver, ee_pr(pr_errest, eebuf), pr_stddev, totaldrift); + if (bFee) + { + fprintf(stdout, " %10g", fee[i]); + } + + fprintf(stdout, " (%s)\n", enm[set[i]].unit); + + if (bFluct) + { + for (j = 0; (j < edat->nframes); j++) + { + edat->s[i].ener[j] -= aver; + } + } + } + if (bSum) + { + totaldrift = (edat->nsteps - 1)*esum->s[0].slope; + fprintf(stdout, "%-24s %10g %10s %10s %10g (%s)", + "Total", esum->s[0].av/nmol, ee_pr(esum->s[0].ee/nmol, eebuf), + "--", totaldrift/nmol, enm[set[0]].unit); + /* pr_aver,pr_stddev,a,totaldrift */ + if (bFee) + { + fprintf(stdout, " %10g %10g\n", + log(expEtot)/beta + esum->s[0].av/nmol, log(expEtot)/beta); + } + else + { + fprintf(stdout, "\n"); + } + } + + /* Do correlation function */ + if (edat->nframes > 1) + { + Dt = delta_t/(edat->nframes - 1); + } + else + { + Dt = 0; + } + if (bVisco) + { + const char* leg[] = { "Shear", "Bulk" }; + real factor; + real **eneset; + real **eneint; + + /* Assume pressure tensor is in Pxx Pxy Pxz Pyx Pyy Pyz Pzx Pzy Pzz */ + + /* Symmetrise tensor! (and store in first three elements) + * And subtract average pressure! + */ + snew(eneset, 12); + for (i = 0; i < 12; i++) + { + snew(eneset[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[1][i] = 0.5*(edat->s[2].ener[i]+edat->s[6].ener[i]); + eneset[2][i] = 0.5*(edat->s[5].ener[i]+edat->s[7].ener[i]); + for (j = 3; j <= 11; j++) + { + eneset[j][i] = edat->s[j].ener[i]; + } + eneset[11][i] -= Pres; + } + + /* 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+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 */ + strcpy(buf, "Shear Viscosity"); + low_do_autocorr(corrfn, oenv, buf, edat->nframes, 3, + (edat->nframes+1)/2, eneset, Dt, + eacNormal, 1, TRUE, FALSE, FALSE, 0.0, 0.0, 0); + + /* Now for bulk viscosity */ + strcpy(buf, "Bulk Viscosity"); + low_do_autocorr(corrfn, oenv, buf, edat->nframes, 1, + (edat->nframes+1)/2, &(eneset[11]), Dt, + eacNormal, 1, TRUE, FALSE, FALSE, 0.0, 0.0, 0); + + factor = (Vaver*1e-26/(BOLTZMANN*Temp))*Dt; + fp = xvgropen(visfn, buf, "Time (ps)", "\\8h\\4 (cp)", oenv); + xvgr_legend(fp, asize(leg), leg, oenv); + + /* Use trapezium rule for integration */ + integral = 0; + intBulk = 0; + nout = get_acfnout(); + if ((nout < 2) || (nout >= edat->nframes/2)) + { + nout = edat->nframes/2; + } + for (i = 1; (i < nout); i++) + { + integral += 0.5*(eneset[0][i-1] + eneset[0][i])*factor; + intBulk += 0.5*(eneset[11][i-1] + eneset[11][i])*factor; + 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) + { + if (bFluct) + { + strcpy(buf, "Autocorrelation of Energy Fluctuations"); + } + else + { + strcpy(buf, "Energy Autocorrelation"); + } +#if 0 + do_autocorr(corrfn, oenv, buf, edat->nframes, + bSum ? 1 : nset, + bSum ? &edat->s[nset-1].ener : eneset, + (delta_t/edat->nframes), eacNormal, FALSE); +#endif + } + } +} + +static void print_time(FILE *fp, double t) +{ + fprintf(fp, "%12.6f", t); +} + +static void print1(FILE *fp, gmx_bool bDp, real e) +{ + if (bDp) + { + fprintf(fp, " %16.12f", e); + } + else + { + fprintf(fp, " %10.6f", e); + } +} + +static void fec(const char *ene2fn, const char *runavgfn, + real reftemp, int nset, int set[], char *leg[], + enerdata_t *edat, double time[], + const output_env_t oenv) +{ + const char * ravgleg[] = { + "\\8D\\4E = E\\sB\\N-E\\sA\\N", + "\\s0..t\\N" + }; + FILE *fp; + ener_file_t enx; + int nre, timecheck, step, nenergy, nenergy2, maxenergy; + int i, j; + gmx_bool bCont; + real aver, beta; + real **eneset2; + double dE, sum; + gmx_enxnm_t *enm = NULL; + t_enxframe *fr; + char buf[22]; + + /* read second energy file */ + snew(fr, 1); + enm = NULL; + enx = open_enx(ene2fn, "r"); + do_enxnms(enx, &(fr->nre), &enm); + + snew(eneset2, nset+1); + nenergy2 = 0; + maxenergy = 0; + timecheck = 0; + do + { + /* This loop searches for the first frame (when -b option is given), + * or when this has been found it reads just one energy frame + */ + do + { + bCont = do_enx(enx, fr); + + if (bCont) + { + timecheck = check_times(fr->t); + } + + } + while (bCont && (timecheck < 0)); + + /* Store energies for analysis afterwards... */ + if ((timecheck == 0) && bCont) + { + if (fr->nre > 0) + { + if (nenergy2 >= maxenergy) + { + maxenergy += 1000; + for (i = 0; i <= nset; i++) + { + srenew(eneset2[i], maxenergy); + } + } + if (fr->t != time[nenergy2]) + { + fprintf(stderr, "\nWARNING time mismatch %g!=%g at frame %s\n", + fr->t, time[nenergy2], gmx_step_str(fr->step, buf)); + } + for (i = 0; i < nset; i++) + { + eneset2[i][nenergy2] = fr->ener[set[i]].e; + } + nenergy2++; + } + } + } + while (bCont && (timecheck == 0)); + + /* check */ + if (edat->nframes != nenergy2) + { + fprintf(stderr, "\nWARNING file length mismatch %d!=%d\n", + edat->nframes, nenergy2); + } + nenergy = min(edat->nframes, nenergy2); + + /* calculate fe difference dF = -kT ln < exp(-(E_B-E_A)/kT) >_A */ + fp = NULL; + if (runavgfn) + { + fp = xvgropen(runavgfn, "Running average free energy difference", + "Time (" unit_time ")", "\\8D\\4E (" unit_energy ")", oenv); + xvgr_legend(fp, asize(ravgleg), ravgleg, oenv); + } + fprintf(stdout, "\n%-24s %10s\n", + "Energy", "dF = -kT ln < exp(-(EB-EA)/kT) >A"); + sum = 0; + beta = 1.0/(BOLTZ*reftemp); + for (i = 0; i < nset; i++) + { + if (gmx_strcasecmp(leg[i], enm[set[i]].name) != 0) + { + fprintf(stderr, "\nWARNING energy set name mismatch %s!=%s\n", + leg[i], enm[set[i]].name); + } + for (j = 0; j < nenergy; j++) + { + dE = eneset2[i][j] - edat->s[i].ener[j]; + sum += exp(-dE*beta); + if (fp) + { + fprintf(fp, "%10g %10g %10g\n", + time[j], dE, -BOLTZ*reftemp*log(sum/(j+1)) ); + } + } + aver = -BOLTZ*reftemp*log(sum/nenergy); + fprintf(stdout, "%-24s %10g\n", leg[i], aver); + } + if (fp) + { + gmx_ffclose(fp); + } + sfree(fr); +} + + +static void do_dhdl(t_enxframe *fr, t_inputrec *ir, FILE **fp_dhdl, + const char *filename, gmx_bool bDp, + int *blocks, int *hists, int *samples, int *nlambdas, + const output_env_t oenv) +{ + const char *dhdl = "dH/d\\lambda", *deltag = "\\DeltaH", *lambda = "\\lambda"; + char title[STRLEN], label_x[STRLEN], label_y[STRLEN], legend[STRLEN]; + char buf[STRLEN]; + gmx_bool first = FALSE; + int nblock_hist = 0, nblock_dh = 0, nblock_dhcoll = 0; + int i, j, k; + /* coll data */ + double temp = 0, start_time = 0, delta_time = 0, start_lambda = 0, delta_lambda = 0; + static int setnr = 0; + double *native_lambda_vec = NULL; + const char **lambda_components = NULL; + int n_lambda_vec = 0; + gmx_bool changing_lambda = FALSE; + int lambda_fep_state; + + /* now count the blocks & handle the global dh data */ + for (i = 0; i < fr->nblock; i++) + { + if (fr->block[i].id == enxDHHIST) + { + nblock_hist++; + } + else if (fr->block[i].id == enxDH) + { + nblock_dh++; + } + else if (fr->block[i].id == enxDHCOLL) + { + nblock_dhcoll++; + if ( (fr->block[i].nsub < 1) || + (fr->block[i].sub[0].type != xdr_datatype_double) || + (fr->block[i].sub[0].nr < 5)) + { + gmx_fatal(FARGS, "Unexpected block data"); + } + + /* read the data from the DHCOLL block */ + temp = fr->block[i].sub[0].dval[0]; + start_time = fr->block[i].sub[0].dval[1]; + delta_time = fr->block[i].sub[0].dval[2]; + start_lambda = fr->block[i].sub[0].dval[3]; + delta_lambda = fr->block[i].sub[0].dval[4]; + changing_lambda = (delta_lambda != 0); + if (fr->block[i].nsub > 1) + { + lambda_fep_state = fr->block[i].sub[1].ival[0]; + if (n_lambda_vec == 0) + { + n_lambda_vec = fr->block[i].sub[1].ival[1]; + } + else + { + if (n_lambda_vec != fr->block[i].sub[1].ival[1]) + { + gmx_fatal(FARGS, + "Unexpected change of basis set in lambda"); + } + } + if (lambda_components == NULL) + { + snew(lambda_components, n_lambda_vec); + } + if (native_lambda_vec == NULL) + { + snew(native_lambda_vec, n_lambda_vec); + } + for (j = 0; j < n_lambda_vec; j++) + { + native_lambda_vec[j] = fr->block[i].sub[0].dval[5+j]; + lambda_components[j] = + efpt_singular_names[fr->block[i].sub[1].ival[2+j]]; + } + } + } + } + + if (nblock_hist == 0 && nblock_dh == 0) + { + /* don't do anything */ + return; + } + if (nblock_hist > 0 && nblock_dh > 0) + { + gmx_fatal(FARGS, "This energy file contains both histogram dhdl data and non-histogram dhdl data. Don't know what to do."); + } + if (!*fp_dhdl) + { + if (nblock_dh > 0) + { + /* we have standard, non-histogram data -- + call open_dhdl to open the file */ + /* TODO this is an ugly hack that needs to be fixed: this will only + work if the order of data is always the same and if we're + only using the g_energy compiled with the mdrun that produced + the ener.edr. */ + *fp_dhdl = open_dhdl(filename, ir, oenv); + } + else + { + sprintf(title, "N(%s)", deltag); + sprintf(label_x, "%s (%s)", deltag, unit_energy); + sprintf(label_y, "Samples"); + *fp_dhdl = xvgropen_type(filename, title, label_x, label_y, exvggtXNY, oenv); + sprintf(buf, "T = %g (K), %s = %g", temp, lambda, start_lambda); + xvgr_subtitle(*fp_dhdl, buf, oenv); + } + } + + (*hists) += nblock_hist; + (*blocks) += nblock_dh; + (*nlambdas) = nblock_hist+nblock_dh; + + /* write the data */ + if (nblock_hist > 0) + { + gmx_int64_t sum = 0; + /* histograms */ + for (i = 0; i < fr->nblock; i++) + { + t_enxblock *blk = &(fr->block[i]); + if (blk->id == enxDHHIST) + { + double foreign_lambda, dx; + 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_int64) || + (blk->sub[0].nr < 2) || + (blk->sub[1].nr < 2) ) + { + gmx_fatal(FARGS, "Unexpected block data in file"); + } + foreign_lambda = blk->sub[0].dval[0]; + dx = blk->sub[0].dval[1]; + nhist = blk->sub[1].lval[0]; + derivative = blk->sub[1].lval[1]; + for (j = 0; j < nhist; j++) + { + const char *lg[1]; + x0 = blk->sub[1].lval[2+j]; + + if (!derivative) + { + sprintf(legend, "N(%s(%s=%g) | %s=%g)", + deltag, lambda, foreign_lambda, + lambda, start_lambda); + } + else + { + sprintf(legend, "N(%s | %s=%g)", + dhdl, lambda, start_lambda); + } + + lg[0] = legend; + xvgr_new_dataset(*fp_dhdl, setnr, 1, lg, oenv); + setnr++; + for (k = 0; k < blk->sub[j+2].nr; k++) + { + int hist; + double xmin, xmax; + + hist = blk->sub[j+2].ival[k]; + xmin = (x0+k)*dx; + xmax = (x0+k+1)*dx; + fprintf(*fp_dhdl, "%g %d\n%g %d\n", xmin, hist, + xmax, hist); + sum += hist; + } + /* multiple histogram data blocks in one histogram + mean that the second one is the reverse of the first one: + for dhdl derivatives, it's important to know both the + maximum and minimum values */ + dx = -dx; + } + } + } + (*samples) += (int)(sum/nblock_hist); + } + else + { + /* raw dh */ + int len = 0; + char **setnames = NULL; + int nnames = nblock_dh; + + for (i = 0; i < fr->nblock; i++) + { + t_enxblock *blk = &(fr->block[i]); + if (blk->id == enxDH) + { + if (len == 0) + { + len = blk->sub[2].nr; + } + else + { + if (len != blk->sub[2].nr) + { + gmx_fatal(FARGS, "Length inconsistency in dhdl data"); + } + } + } + } + (*samples) += len; + + for (i = 0; i < len; i++) + { + double time = start_time + delta_time*i; + + fprintf(*fp_dhdl, "%.4f ", time); + + for (j = 0; j < fr->nblock; j++) + { + t_enxblock *blk = &(fr->block[j]); + if (blk->id == enxDH) + { + double value; + if (blk->sub[2].type == xdr_datatype_float) + { + value = blk->sub[2].fval[i]; + } + else + { + value = blk->sub[2].dval[i]; + } + /* we need to decide which data type it is based on the count*/ + + if (j == 1 && ir->bExpanded) + { + fprintf(*fp_dhdl, "%4d", (int)value); /* if expanded ensembles and zero, this is a state value, it's an integer. We need a cleaner conditional than if j==1! */ + } + else + { + if (bDp) + { + fprintf(*fp_dhdl, " %#.12g", value); /* print normal precision */ + } + else + { + fprintf(*fp_dhdl, " %#.8g", value); /* print normal precision */ + } + } + } + } + fprintf(*fp_dhdl, "\n"); + } + } +} + + +int gmx_energy(int argc, char *argv[]) +{ + const char *desc[] = { + "[THISMODULE] extracts energy components or distance restraint", + "data from an energy file. The user is prompted to interactively", + "select the desired energy terms.[PAR]", + + "Average, RMSD, and drift are calculated with full precision from the", + "simulation (see printed manual). Drift is calculated by performing", + "a least-squares fit of the data to a straight line. The reported total drift", + "is the difference of the fit at the first and last point.", + "An error estimate of the average is given based on a block averages", + "over 5 blocks using the full-precision averages. The error estimate", + "can be performed over multiple block lengths with the options", + "[TT]-nbmin[tt] and [TT]-nbmax[tt].", + "[BB]Note[bb] that in most cases the energy files contains averages over all", + "MD steps, or over many more points than the number of frames in", + "energy file. This makes the [THISMODULE] statistics output more accurate", + "than the [TT].xvg[tt] output. When exact averages are not present in the energy", + "file, the statistics mentioned above are simply over the single, per-frame", + "energy values.[PAR]", + + "The term fluctuation gives the RMSD around the least-squares fit.[PAR]", + + "Some fluctuation-dependent properties can be calculated provided", + "the correct energy terms are selected, and that the command line option", + "[TT]-fluct_props[tt] is given. The following properties", + "will be computed:[BR]", + "Property Energy terms needed[BR]", + "---------------------------------------------------[BR]", + "Heat capacity C[SUB]p[sub] (NPT sims): Enthalpy, Temp [BR]", + "Heat capacity C[SUB]v[sub] (NVT sims): Etot, Temp [BR]", + "Thermal expansion coeff. (NPT): Enthalpy, Vol, Temp[BR]", + "Isothermal compressibility: Vol, Temp [BR]", + "Adiabatic bulk modulus: Vol, Temp [BR]", + "---------------------------------------------------[BR]", + "You always need to set the number of molecules [TT]-nmol[tt].", + "The C[SUB]p[sub]/C[SUB]v[sub] computations do [BB]not[bb] include any corrections", + "for quantum effects. Use the [gmx-dos] program if you need that (and you do).[PAR]" + "When the [TT]-viol[tt] option is set, the time averaged", + "violations are plotted and the running time-averaged and", + "instantaneous sum of violations are recalculated. Additionally", + "running time-averaged and instantaneous distances between", + "selected pairs can be plotted with the [TT]-pairs[tt] option.[PAR]", + + "Options [TT]-ora[tt], [TT]-ort[tt], [TT]-oda[tt], [TT]-odr[tt] and", + "[TT]-odt[tt] are used for analyzing orientation restraint data.", + "The first two options plot the orientation, the last three the", + "deviations of the orientations from the experimental values.", + "The options that end on an 'a' plot the average over time", + "as a function of restraint. The options that end on a 't'", + "prompt the user for restraint label numbers and plot the data", + "as a function of time. Option [TT]-odr[tt] plots the RMS", + "deviation as a function of restraint.", + "When the run used time or ensemble averaged orientation restraints,", + "option [TT]-orinst[tt] can be used to analyse the instantaneous,", + "not ensemble-averaged orientations and deviations instead of", + "the time and ensemble averages.[PAR]", + + "Option [TT]-oten[tt] plots the eigenvalues of the molecular order", + "tensor for each orientation restraint experiment. With option", + "[TT]-ovec[tt] also the eigenvectors are plotted.[PAR]", + + "Option [TT]-odh[tt] extracts and plots the free energy data", + "(Hamiltoian differences and/or the Hamiltonian derivative dhdl)", + "from the [TT]ener.edr[tt] file.[PAR]", + + "With [TT]-fee[tt] an estimate is calculated for the free-energy", + "difference with an ideal gas state: [BR]", + " [GRK]Delta[grk] A = A(N,V,T) - A[SUB]idealgas[sub](N,V,T) = kT [LN][CHEVRON][EXP]U[SUB]pot[sub]/kT[exp][chevron][ln][BR]", + " [GRK]Delta[grk] G = G(N,p,T) - G[SUB]idealgas[sub](N,p,T) = kT [LN][CHEVRON][EXP]U[SUB]pot[sub]/kT[exp][chevron][ln][BR]", + "where k is Boltzmann's constant, T is set by [TT]-fetemp[tt] and", + "the average is over the ensemble (or time in a trajectory).", + "Note that this is in principle", + "only correct when averaging over the whole (Boltzmann) ensemble", + "and using the potential energy. This also allows for an entropy", + "estimate using:[BR]", + " [GRK]Delta[grk] S(N,V,T) = S(N,V,T) - S[SUB]idealgas[sub](N,V,T) = ([CHEVRON]U[SUB]pot[sub][chevron] - [GRK]Delta[grk] A)/T[BR]", + " [GRK]Delta[grk] S(N,p,T) = S(N,p,T) - S[SUB]idealgas[sub](N,p,T) = ([CHEVRON]U[SUB]pot[sub][chevron] + pV - [GRK]Delta[grk] G)/T", + "[PAR]", + + "When a second energy file is specified ([TT]-f2[tt]), a free energy", + "difference is calculated [BR] dF = -kT [LN][CHEVRON][EXP]-(E[SUB]B[sub]-E[SUB]A[sub])/kT[exp][chevron][SUB]A[sub][ln] ,", + "where E[SUB]A[sub] and E[SUB]B[sub] are the energies from the first and second energy", + "files, and the average is over the ensemble A. The running average", + "of the free energy difference is printed to a file specified by [TT]-ravg[tt].", + "[BB]Note[bb] that the energies must both be calculated from the same trajectory." + + }; + static gmx_bool bSum = FALSE, bFee = FALSE, bPrAll = FALSE, bFluct = FALSE, bDriftCorr = FALSE; + static gmx_bool bDp = FALSE, bMutot = FALSE, bOrinst = FALSE, bOvec = FALSE, bFluctProps = FALSE; + static int skip = 0, nmol = 1, nbmin = 5, nbmax = 5; + static real reftemp = 300.0, ezero = 0; + t_pargs pa[] = { + { "-fee", FALSE, etBOOL, {&bFee}, + "Do a free energy estimate" }, + { "-fetemp", FALSE, etREAL, {&reftemp}, + "Reference temperature for free energy calculation" }, + { "-zero", FALSE, etREAL, {&ezero}, + "Subtract a zero-point energy" }, + { "-sum", FALSE, etBOOL, {&bSum}, + "Sum the energy terms selected rather than display them all" }, + { "-dp", FALSE, etBOOL, {&bDp}, + "Print energies in high precision" }, + { "-nbmin", FALSE, etINT, {&nbmin}, + "Minimum number of blocks for error estimate" }, + { "-nbmax", FALSE, etINT, {&nbmax}, + "Maximum number of blocks for error estimate" }, + { "-mutot", FALSE, etBOOL, {&bMutot}, + "Compute the total dipole moment from the components" }, + { "-skip", FALSE, etINT, {&skip}, + "Skip number of frames between data points" }, + { "-aver", FALSE, etBOOL, {&bPrAll}, + "Also print the exact average and rmsd stored in the energy frames (only when 1 term is requested)" }, + { "-nmol", FALSE, etINT, {&nmol}, + "Number of molecules in your sample: the energies are divided by this number" }, + { "-fluct_props", FALSE, etBOOL, {&bFluctProps}, + "Compute properties based on energy fluctuations, like heat capacity" }, + { "-driftcorr", FALSE, etBOOL, {&bDriftCorr}, + "Useful only for calculations of fluctuation properties. The drift in the observables will be subtracted before computing the fluctuation properties."}, + { "-fluc", FALSE, etBOOL, {&bFluct}, + "Calculate autocorrelation of energy fluctuations rather than energy itself" }, + { "-orinst", FALSE, etBOOL, {&bOrinst}, + "Analyse instantaneous orientation data" }, + { "-ovec", FALSE, etBOOL, {&bOvec}, + "Also plot the eigenvectors with [TT]-oten[tt]" } + }; + const char * drleg[] = { + "Running average", + "Instantaneous" + }; + static const char *setnm[] = { + "Pres-XX", "Pres-XY", "Pres-XZ", "Pres-YX", "Pres-YY", + "Pres-YZ", "Pres-ZX", "Pres-ZY", "Pres-ZZ", "Temperature", + "Volume", "Pressure" + }; + + FILE *out = NULL, *fp_pairs = NULL, *fort = NULL, *fodt = NULL, *foten = NULL; + FILE *fp_dhdl = NULL; + FILE **drout; + ener_file_t fp; + int timecheck = 0; + gmx_mtop_t mtop; + gmx_localtop_t *top = NULL; + t_inputrec ir; + t_energy **ee; + enerdata_t edat; + gmx_enxnm_t *enm = NULL; + t_enxframe *frame, *fr = NULL; + int cur = 0; +#define NEXT (1-cur) + int nre, teller, teller_disre, nfr; + 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; + int *index = NULL, *pair = NULL, norsel = 0, *orsel = NULL, *or_label = NULL; + int nbounds = 0, npairs; + gmx_bool bDisRe, bDRAll, bORA, bORT, bODA, bODR, bODT, bORIRE, bOTEN, bDHDL; + gmx_bool bFoundStart, bCont, bEDR, bVisco; + double sum, sumaver, sumt, ener, dbl; + double *time = NULL; + real Vaver; + int *set = NULL, i, j, k, nset, sss; + gmx_bool *bIsEner = NULL; + char **pairleg, **odtleg, **otenleg; + char **leg = NULL; + char **nms; + char *anm_j, *anm_k, *resnm_j, *resnm_k; + int resnr_j, resnr_k; + const char *orinst_sub = "@ subtitle \"instantaneous\"\n"; + char buf[256]; + output_env_t oenv; + t_enxblock *blk = NULL; + t_enxblock *blk_disre = NULL; + int ndisre = 0; + int dh_blocks = 0, dh_hists = 0, dh_samples = 0, dh_lambdas = 0; + + t_filenm fnm[] = { + { efEDR, "-f", NULL, ffREAD }, + { efEDR, "-f2", NULL, ffOPTRD }, + { efTPX, "-s", NULL, ffOPTRD }, + { efXVG, "-o", "energy", ffWRITE }, + { efXVG, "-viol", "violaver", ffOPTWR }, + { efXVG, "-pairs", "pairs", ffOPTWR }, + { efXVG, "-ora", "orienta", ffOPTWR }, + { efXVG, "-ort", "orientt", ffOPTWR }, + { efXVG, "-oda", "orideva", ffOPTWR }, + { efXVG, "-odr", "oridevr", ffOPTWR }, + { efXVG, "-odt", "oridevt", ffOPTWR }, + { efXVG, "-oten", "oriten", ffOPTWR }, + { efXVG, "-corr", "enecorr", ffOPTWR }, + { efXVG, "-vis", "visco", ffOPTWR }, + { efXVG, "-ravg", "runavgdf", ffOPTWR }, + { efXVG, "-odh", "dhdl", ffOPTWR } + }; +#define NFILE asize(fnm) + int npargs; + t_pargs *ppa; + + 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, + NFILE, fnm, npargs, ppa, asize(desc), desc, 0, NULL, &oenv)) + { + return 0; + } + + bDRAll = opt2bSet("-pairs", NFILE, fnm); + bDisRe = opt2bSet("-viol", NFILE, fnm) || bDRAll; + bORA = opt2bSet("-ora", NFILE, fnm); + bORT = opt2bSet("-ort", NFILE, fnm); + bODA = opt2bSet("-oda", NFILE, fnm); + bODR = opt2bSet("-odr", NFILE, fnm); + bODT = opt2bSet("-odt", NFILE, fnm); + bORIRE = bORA || bORT || bODA || bODR || bODT; + bOTEN = opt2bSet("-oten", NFILE, fnm); + bDHDL = opt2bSet("-odh", NFILE, fnm); + + nset = 0; + + snew(frame, 2); + fp = open_enx(ftp2fn(efEDR, NFILE, fnm), "r"); + do_enxnms(fp, &nre, &enm); + + Vaver = -1; + + bVisco = opt2bSet("-vis", NFILE, fnm); + + if ((!bDisRe) && (!bDHDL)) + { + if (bVisco) + { + nset = asize(setnm); + snew(set, nset); + /* This is nasty code... To extract Pres tensor, Volume and Temperature */ + for (j = 0; j < nset; j++) + { + for (i = 0; i < nre; i++) + { + if (strstr(enm[i].name, setnm[j])) + { + set[j] = i; + break; + } + } + if (i == nre) + { + if (gmx_strcasecmp(setnm[j], "Volume") == 0) + { + printf("Enter the box volume (" unit_volume "): "); + if (1 != scanf("%lf", &dbl)) + { + gmx_fatal(FARGS, "Error reading user input"); + } + Vaver = dbl; + } + else + { + gmx_fatal(FARGS, "Could not find term %s for viscosity calculation", + setnm[j]); + } + } + } + } + else + { + set = select_by_name(nre, enm, &nset); + } + /* Print all the different units once */ + sprintf(buf, "(%s)", enm[set[0]].unit); + for (i = 1; i < nset; i++) + { + for (j = 0; j < i; j++) + { + if (strcmp(enm[set[i]].unit, enm[set[j]].unit) == 0) + { + break; + } + } + if (j == i) + { + strcat(buf, ", ("); + strcat(buf, enm[set[i]].unit); + strcat(buf, ")"); + } + } + out = xvgropen(opt2fn("-o", NFILE, fnm), "Gromacs Energies", "Time (ps)", buf, + oenv); + + snew(leg, nset+1); + for (i = 0; (i < nset); i++) + { + leg[i] = enm[set[i]].name; + } + if (bSum) + { + leg[nset] = strdup("Sum"); + xvgr_legend(out, nset+1, (const char**)leg, oenv); + } + else + { + xvgr_legend(out, nset, (const char**)leg, oenv); + } + + snew(bIsEner, nset); + for (i = 0; (i < nset); i++) + { + bIsEner[i] = FALSE; + for (j = 0; (j <= F_ETOT); j++) + { + bIsEner[i] = bIsEner[i] || + (gmx_strcasecmp(interaction_function[j].longname, leg[i]) == 0); + } + } + + if (bPrAll && nset > 1) + { + gmx_fatal(FARGS, "Printing averages can only be done when a single set is selected"); + } + + time = NULL; + + if (bORIRE || bOTEN) + { + get_orires_parms(ftp2fn(efTPX, NFILE, fnm), &nor, &nex, &or_label, &oobs); + } + + if (bORIRE) + { + if (bOrinst) + { + enx_i = enxORI; + } + else + { + enx_i = enxOR; + } + + if (bORA || bODA) + { + snew(orient, nor); + } + if (bODR) + { + snew(odrms, nor); + } + if (bORT || bODT) + { + fprintf(stderr, "Select the orientation restraint labels you want (-1 is all)\n"); + fprintf(stderr, "End your selection with 0\n"); + j = -1; + orsel = NULL; + do + { + j++; + srenew(orsel, j+1); + if (1 != scanf("%d", &(orsel[j]))) + { + gmx_fatal(FARGS, "Error reading user input"); + } + } + while (orsel[j] > 0); + if (orsel[0] == -1) + { + fprintf(stderr, "Selecting all %d orientation restraints\n", nor); + norsel = nor; + srenew(orsel, nor); + for (i = 0; i < nor; i++) + { + orsel[i] = i; + } + } + else + { + /* Build the selection */ + norsel = 0; + for (i = 0; i < j; i++) + { + for (k = 0; k < nor; k++) + { + if (or_label[k] == orsel[i]) + { + orsel[norsel] = k; + norsel++; + break; + } + } + if (k == nor) + { + fprintf(stderr, "Orientation restraint label %d not found\n", + orsel[i]); + } + } + } + snew(odtleg, norsel); + for (i = 0; i < norsel; i++) + { + snew(odtleg[i], 256); + sprintf(odtleg[i], "%d", or_label[orsel[i]]); + } + if (bORT) + { + fort = xvgropen(opt2fn("-ort", NFILE, fnm), "Calculated orientations", + "Time (ps)", "", oenv); + if (bOrinst && output_env_get_print_xvgr_codes(oenv)) + { + fprintf(fort, "%s", orinst_sub); + } + xvgr_legend(fort, norsel, (const char**)odtleg, oenv); + } + if (bODT) + { + fodt = xvgropen(opt2fn("-odt", NFILE, fnm), + "Orientation restraint deviation", + "Time (ps)", "", oenv); + if (bOrinst && output_env_get_print_xvgr_codes(oenv)) + { + fprintf(fodt, "%s", orinst_sub); + } + xvgr_legend(fodt, norsel, (const char**)odtleg, oenv); + } + } + } + if (bOTEN) + { + foten = xvgropen(opt2fn("-oten", NFILE, fnm), + "Order tensor", "Time (ps)", "", oenv); + snew(otenleg, bOvec ? nex*12 : nex*3); + for (i = 0; i < nex; i++) + { + for (j = 0; j < 3; j++) + { + sprintf(buf, "eig%d", j+1); + otenleg[(bOvec ? 12 : 3)*i+j] = 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); + } + } + } + xvgr_legend(foten, bOvec ? nex*12 : nex*3, (const char**)otenleg, oenv); + } + } + else if (bDisRe) + { + nbounds = get_bounds(ftp2fn(efTPX, NFILE, fnm), &bounds, &index, &pair, &npairs, + &mtop, &top, &ir); + snew(violaver, npairs); + out = xvgropen(opt2fn("-o", NFILE, fnm), "Sum of Violations", + "Time (ps)", "nm", oenv); + xvgr_legend(out, 2, drleg, oenv); + if (bDRAll) + { + fp_pairs = xvgropen(opt2fn("-pairs", NFILE, fnm), "Pair Distances", + "Time (ps)", "Distance (nm)", oenv); + if (output_env_get_print_xvgr_codes(oenv)) + { + fprintf(fp_pairs, "@ subtitle \"averaged (tau=%g) and instantaneous\"\n", + ir.dr_tau); + } + } + } + else if (bDHDL) + { + get_dhdl_parms(ftp2fn(efTPX, NFILE, fnm), &ir); + } + + /* Initiate energies and set them to zero */ - edat.nsteps = 0; - edat.npoints = 0; - edat.nframes = 0; - edat.step = NULL; - edat.steps = NULL; - edat.points = NULL; ++ edat.nsteps = 0; ++ edat.npoints = 0; ++ edat.nframes = 0; ++ edat.step = NULL; ++ edat.steps = NULL; ++ edat.points = NULL; ++ edat.bHaveSums = TRUE; + snew(edat.s, nset); + + /* Initiate counters */ + teller = 0; + teller_disre = 0; + bFoundStart = FALSE; + start_step = 0; + start_t = 0; + do + { + /* This loop searches for the first frame (when -b option is given), + * or when this has been found it reads just one energy frame + */ + do + { + bCont = do_enx(fp, &(frame[NEXT])); + if (bCont) + { + timecheck = check_times(frame[NEXT].t); + } + } + while (bCont && (timecheck < 0)); + + if ((timecheck == 0) && bCont) + { + /* We read a valid frame, so we can use it */ + fr = &(frame[NEXT]); + + if (fr->nre > 0) + { + /* The frame contains energies, so update cur */ + cur = NEXT; + + if (edat.nframes % 1000 == 0) + { + srenew(edat.step, edat.nframes+1000); + memset(&(edat.step[edat.nframes]), 0, 1000*sizeof(edat.step[0])); + srenew(edat.steps, edat.nframes+1000); + memset(&(edat.steps[edat.nframes]), 0, 1000*sizeof(edat.steps[0])); + srenew(edat.points, edat.nframes+1000); + memset(&(edat.points[edat.nframes]), 0, 1000*sizeof(edat.points[0])); + + for (i = 0; i < nset; i++) + { + srenew(edat.s[i].ener, edat.nframes+1000); + memset(&(edat.s[i].ener[edat.nframes]), 0, + 1000*sizeof(edat.s[i].ener[0])); + srenew(edat.s[i].es, edat.nframes+1000); + memset(&(edat.s[i].es[edat.nframes]), 0, + 1000*sizeof(edat.s[i].es[0])); + } + } + + nfr = edat.nframes; + edat.step[nfr] = fr->step; + + if (!bFoundStart) + { + bFoundStart = TRUE; + /* Initiate the previous step data */ + start_step = fr->step; + start_t = fr->t; + /* Initiate the energy sums */ + edat.steps[nfr] = 1; + edat.points[nfr] = 1; + for (i = 0; i < nset; i++) + { + sss = set[i]; + edat.s[i].es[nfr].sum = fr->ener[sss].e; + edat.s[i].es[nfr].sum2 = 0; + } + edat.nsteps = 1; + edat.npoints = 1; + } + else + { + edat.steps[nfr] = fr->nsteps; ++ ++ if (fr->nsum <= 1) + { - if (fr->step - start_step + 1 == edat.nsteps + fr->nsteps) ++ /* mdrun only calculated the energy at energy output ++ * steps. We don't need to check step intervals. ++ */ ++ edat.points[nfr] = 1; ++ for (i = 0; i < nset; i++) + { - if (fr->nsum <= 1) - { - edat.points[nfr] = 1; - for (i = 0; i < nset; i++) - { - sss = set[i]; - edat.s[i].es[nfr].sum = fr->ener[sss].e; - edat.s[i].es[nfr].sum2 = 0; - } - edat.npoints += 1; - } - else - { - edat.points[nfr] = fr->nsum; - for (i = 0; i < nset; i++) - { - sss = set[i]; - edat.s[i].es[nfr].sum = fr->ener[sss].esum; - edat.s[i].es[nfr].sum2 = fr->ener[sss].eav; - } - edat.npoints += fr->nsum; - } ++ sss = set[i]; ++ edat.s[i].es[nfr].sum = fr->ener[sss].e; ++ edat.s[i].es[nfr].sum2 = 0; + } - else ++ edat.npoints += 1; ++ edat.bHaveSums = FALSE; ++ } ++ else if (fr->step - start_step + 1 == edat.nsteps + fr->nsteps) ++ { ++ /* We have statistics to the previous frame */ ++ edat.points[nfr] = fr->nsum; ++ for (i = 0; i < nset; i++) + { - /* The interval does not match fr->nsteps: - * can not do exact averages. - */ - edat.npoints = 0; ++ sss = set[i]; ++ edat.s[i].es[nfr].sum = fr->ener[sss].esum; ++ edat.s[i].es[nfr].sum2 = fr->ener[sss].eav; + } - edat.nsteps = fr->step - start_step + 1; ++ edat.npoints += fr->nsum; + } ++ else ++ { ++ /* The interval does not match fr->nsteps: ++ * can not do exact averages. ++ */ ++ edat.bHaveSums = FALSE; ++ } ++ ++ edat.nsteps = fr->step - start_step + 1; + } + for (i = 0; i < nset; i++) + { + edat.s[i].ener[nfr] = fr->ener[set[i]].e; + } + } + /* + * Define distance restraint legends. Can only be done after + * the first frame has been read... (Then we know how many there are) + */ + blk_disre = find_block_id_enxframe(fr, enxDISRE, NULL); + if (bDisRe && bDRAll && !leg && blk_disre) + { + t_iatom *fa; + t_iparams *ip; + + fa = top->idef.il[F_DISRES].iatoms; + ip = top->idef.iparams; + if (blk_disre->nsub != 2 || + (blk_disre->sub[0].nr != blk_disre->sub[1].nr) ) + { + gmx_incons("Number of disre sub-blocks not equal to 2"); + } + + ndisre = blk_disre->sub[0].nr; + if (ndisre != top->idef.il[F_DISRES].nr/3) + { + gmx_fatal(FARGS, "Number of disre pairs in the energy file (%d) does not match the number in the run input file (%d)\n", + ndisre, top->idef.il[F_DISRES].nr/3); + } + snew(pairleg, ndisre); + for (i = 0; i < ndisre; i++) + { + snew(pairleg[i], 30); + j = fa[3*i+1]; + k = fa[3*i+2]; + gmx_mtop_atominfo_global(&mtop, j, &anm_j, &resnr_j, &resnm_j); + gmx_mtop_atominfo_global(&mtop, k, &anm_k, &resnr_k, &resnm_k); + sprintf(pairleg[i], "%d %s %d %s (%d)", + resnr_j, anm_j, resnr_k, anm_k, + ip[fa[3*i]].disres.label); + } + set = select_it(ndisre, pairleg, &nset); + snew(leg, 2*nset); + for (i = 0; (i < nset); i++) + { + snew(leg[2*i], 32); + sprintf(leg[2*i], "a %s", pairleg[set[i]]); + snew(leg[2*i+1], 32); + sprintf(leg[2*i+1], "i %s", pairleg[set[i]]); + } + xvgr_legend(fp_pairs, 2*nset, (const char**)leg, oenv); + } + + /* + * Store energies for analysis afterwards... + */ + if (!bDisRe && !bDHDL && (fr->nre > 0)) + { + if (edat.nframes % 1000 == 0) + { + srenew(time, edat.nframes+1000); + } + time[edat.nframes] = fr->t; + edat.nframes++; + } + /* + * Printing time, only when we do not want to skip frames + */ + if (!skip || teller % skip == 0) + { + if (bDisRe) + { + /******************************************* + * D I S T A N C E R E S T R A I N T S + *******************************************/ + if (ndisre > 0) + { + #ifndef GMX_DOUBLE + float *disre_rt = blk_disre->sub[0].fval; + float *disre_rm3tav = blk_disre->sub[1].fval; + #else + double *disre_rt = blk_disre->sub[0].dval; + double *disre_rm3tav = blk_disre->sub[1].dval; + #endif + + print_time(out, fr->t); + if (violaver == NULL) + { + snew(violaver, ndisre); + } + + /* Subtract bounds from distances, to calculate violations */ + calc_violations(disre_rt, disre_rm3tav, + nbounds, pair, bounds, violaver, &sumt, &sumaver); + + fprintf(out, " %8.4f %8.4f\n", sumaver, sumt); + if (bDRAll) + { + print_time(fp_pairs, fr->t); + for (i = 0; (i < nset); i++) + { + sss = set[i]; + fprintf(fp_pairs, " %8.4f", mypow(disre_rm3tav[sss], minthird)); + fprintf(fp_pairs, " %8.4f", disre_rt[sss]); + } + fprintf(fp_pairs, "\n"); + } + teller_disre++; + } + } + else if (bDHDL) + { + do_dhdl(fr, &ir, &fp_dhdl, opt2fn("-odh", NFILE, fnm), bDp, &dh_blocks, &dh_hists, &dh_samples, &dh_lambdas, oenv); + } + + /******************************************* + * E N E R G I E S + *******************************************/ + else + { + if (fr->nre > 0) + { + if (bPrAll) + { + /* We skip frames with single points (usually only the first frame), + * since they would result in an average plot with outliers. + */ + if (fr->nsum > 1) + { + print_time(out, fr->t); + print1(out, bDp, fr->ener[set[0]].e); + print1(out, bDp, fr->ener[set[0]].esum/fr->nsum); + print1(out, bDp, sqrt(fr->ener[set[0]].eav/fr->nsum)); + fprintf(out, "\n"); + } + } + else + { + print_time(out, fr->t); + if (bSum) + { + sum = 0; + for (i = 0; i < nset; i++) + { + sum += fr->ener[set[i]].e; + } + print1(out, bDp, sum/nmol-ezero); + } + else + { + for (i = 0; (i < nset); i++) + { + if (bIsEner[i]) + { + print1(out, bDp, (fr->ener[set[i]].e)/nmol-ezero); + } + else + { + print1(out, bDp, fr->ener[set[i]].e); + } + } + } + fprintf(out, "\n"); + } + } + blk = find_block_id_enxframe(fr, enx_i, NULL); + if (bORIRE && blk) + { +#ifndef GMX_DOUBLE + xdr_datatype dt = xdr_datatype_float; +#else + xdr_datatype dt = xdr_datatype_double; +#endif + real *vals; + + if ( (blk->nsub != 1) || (blk->sub[0].type != dt) ) + { + gmx_fatal(FARGS, "Orientational restraints read in incorrectly"); + } +#ifndef GMX_DOUBLE + vals = blk->sub[0].fval; +#else + vals = blk->sub[0].dval; +#endif + + if (blk->sub[0].nr != (size_t)nor) + { + gmx_fatal(FARGS, "Number of orientation restraints in energy file (%d) does not match with the topology (%d)", blk->sub[0].nr); + } + if (bORA || bODA) + { + for (i = 0; i < nor; i++) + { + orient[i] += vals[i]; + } + } + if (bODR) + { + for (i = 0; i < nor; i++) + { + odrms[i] += sqr(vals[i]-oobs[i]); + } + } + if (bORT) + { + fprintf(fort, " %10f", fr->t); + for (i = 0; i < norsel; i++) + { + fprintf(fort, " %g", vals[orsel[i]]); + } + fprintf(fort, "\n"); + } + if (bODT) + { + fprintf(fodt, " %10f", fr->t); + for (i = 0; i < norsel; i++) + { + fprintf(fodt, " %g", vals[orsel[i]]-oobs[orsel[i]]); + } + fprintf(fodt, "\n"); + } + norfr++; + } + blk = find_block_id_enxframe(fr, enxORT, NULL); + if (bOTEN && blk) + { +#ifndef GMX_DOUBLE + xdr_datatype dt = xdr_datatype_float; +#else + xdr_datatype dt = xdr_datatype_double; +#endif + real *vals; + + if ( (blk->nsub != 1) || (blk->sub[0].type != dt) ) + { + gmx_fatal(FARGS, "Orientational restraints read in incorrectly"); + } +#ifndef GMX_DOUBLE + vals = blk->sub[0].fval; +#else + vals = blk->sub[0].dval; +#endif + + if (blk->sub[0].nr != (size_t)(nex*12)) + { + gmx_fatal(FARGS, "Number of orientation experiments in energy file (%g) does not match with the topology (%d)", + blk->sub[0].nr/12, nex); + } + fprintf(foten, " %10f", fr->t); + for (i = 0; i < nex; i++) + { + for (j = 0; j < (bOvec ? 12 : 3); j++) + { + fprintf(foten, " %g", vals[i*12+j]); + } + } + fprintf(foten, "\n"); + } + } + } + teller++; + } + } + while (bCont && (timecheck == 0)); + + fprintf(stderr, "\n"); + close_enx(fp); + if (out) + { + gmx_ffclose(out); + } + + if (bDRAll) + { + gmx_ffclose(fp_pairs); + } + + if (bORT) + { + gmx_ffclose(fort); + } + if (bODT) + { + gmx_ffclose(fodt); + } + if (bORA) + { + out = xvgropen(opt2fn("-ora", NFILE, fnm), + "Average calculated orientations", + "Restraint label", "", oenv); + if (bOrinst && output_env_get_print_xvgr_codes(oenv)) + { + fprintf(out, "%s", orinst_sub); + } + for (i = 0; i < nor; i++) + { + fprintf(out, "%5d %g\n", or_label[i], orient[i]/norfr); + } + gmx_ffclose(out); + } + if (bODA) + { + out = xvgropen(opt2fn("-oda", NFILE, fnm), + "Average restraint deviation", + "Restraint label", "", oenv); + if (bOrinst && output_env_get_print_xvgr_codes(oenv)) + { + fprintf(out, "%s", orinst_sub); + } + for (i = 0; i < nor; i++) + { + fprintf(out, "%5d %g\n", or_label[i], orient[i]/norfr-oobs[i]); + } + gmx_ffclose(out); + } + if (bODR) + { + out = xvgropen(opt2fn("-odr", NFILE, fnm), + "RMS orientation restraint deviations", + "Restraint label", "", oenv); + if (bOrinst && output_env_get_print_xvgr_codes(oenv)) + { + fprintf(out, "%s", orinst_sub); + } + for (i = 0; i < nor; i++) + { + fprintf(out, "%5d %g\n", or_label[i], sqrt(odrms[i]/norfr)); + } + gmx_ffclose(out); + } + if (bOTEN) + { + gmx_ffclose(foten); + } + + if (bDisRe) + { + analyse_disre(opt2fn("-viol", NFILE, fnm), + teller_disre, violaver, bounds, index, pair, nbounds, oenv); + } + else if (bDHDL) + { + if (fp_dhdl) + { + gmx_ffclose(fp_dhdl); + printf("\n\nWrote %d lambda values with %d samples as ", + dh_lambdas, dh_samples); + if (dh_hists > 0) + { + printf("%d dH histograms ", dh_hists); + } + if (dh_blocks > 0) + { + printf("%d dH data blocks ", dh_blocks); + } + printf("to %s\n", opt2fn("-odh", NFILE, fnm)); + + } + else + { + gmx_fatal(FARGS, "No dH data in %s\n", opt2fn("-f", NFILE, fnm)); + } + + } + else + { + double dt = (frame[cur].t-start_t)/(edat.nframes-1); + analyse_ener(opt2bSet("-corr", NFILE, fnm), opt2fn("-corr", NFILE, fnm), + bFee, bSum, opt2parg_bSet("-nmol", npargs, ppa), + bVisco, opt2fn("-vis", NFILE, fnm), + nmol, + start_step, start_t, frame[cur].step, frame[cur].t, + reftemp, &edat, + nset, set, bIsEner, leg, enm, Vaver, ezero, nbmin, nbmax, + oenv); + if (bFluctProps) + { + calc_fluctuation_props(stdout, bDriftCorr, dt, nset, nmol, leg, &edat, + nbmin, nbmax); + } + } + if (opt2bSet("-f2", NFILE, fnm)) + { + fec(opt2fn("-f2", NFILE, fnm), opt2fn("-ravg", NFILE, fnm), + reftemp, nset, set, leg, &edat, time, oenv); + } + + { + const char *nxy = "-nxy"; + + do_view(oenv, opt2fn("-o", NFILE, fnm), nxy); + do_view(oenv, opt2fn_null("-ravg", NFILE, fnm), nxy); + do_view(oenv, opt2fn_null("-ora", NFILE, fnm), nxy); + do_view(oenv, opt2fn_null("-ort", NFILE, fnm), nxy); + do_view(oenv, opt2fn_null("-oda", NFILE, fnm), nxy); + do_view(oenv, opt2fn_null("-odr", NFILE, fnm), nxy); + do_view(oenv, opt2fn_null("-odt", NFILE, fnm), nxy); + do_view(oenv, opt2fn_null("-oten", NFILE, fnm), nxy); + do_view(oenv, opt2fn_null("-odh", NFILE, fnm), nxy); + } + + return 0; +}