--- /dev/null
- int nframes, real **sum,
- real V, real T, int nsteps, double time[],
+/*
+ * 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
+ * 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 <config.h>
+#endif
+
+#include <math.h>
+#include <string.h>
+
+#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;
+} 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("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("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,
- double fac, dt, di;
++ int nint, real **eneint,
++ real V, real T, double dt,
+ const output_env_t oenv)
+{
+ FILE *fp0, *fp1;
+ double av[4], avold[4];
- if (nframes < 1)
- {
- return;
- }
-
- dt = (time[1]-time[0]);
- nf4 = nframes/4+1;
++ double fac, di;
+ int i, j, m, nf4;
+
- for (i = 1; i < nf4; i++)
++ 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);
- fac = dt*nframes/nsteps;
++ for (i = 0; i < nf4; i++)
+ {
- for (j = 0; j < nframes-i; j++)
+ for (m = 0; m <= nsets; m++)
+ {
+ av[m] = 0;
+ }
- di = sqr(fac*(sum[m][j+i]-sum[m][j]));
++ for (j = 0; j < nint - i; j++)
+ {
+ for (m = 0; m < nsets; m++)
+ {
- fac = (V*NANO*NANO*NANO*PICO*1e10)/(2*BOLTZMANN*T)/(nframes-i);
- fprintf(fp0, "%10g", time[i]-time[0]);
++ di = sqr(eneint[m][j+i] - eneint[m][j]);
+
+ av[m] += di;
+ av[nsets] += di/nsets;
+ }
+ }
+ /* Convert to SI for the viscosity */
- fprintf(fp1, "%10g", 0.5*(time[i]+time[i-1])-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");
- double time[], real reftemp,
++ 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)
+ {
+ /* 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 **enesum;
++ 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)
+ {
+ 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<e^(E/kT)>");
+ }
+ 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;
- snew(enesum, 3);
- for (i = 0; i < 3; i++)
- {
- snew(enesum[i], edat->nframes);
- }
++ 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);
+ }
- 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);
+ 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;
- 3, edat->nframes, enesum, Vaver, Temp, nsteps, time, oenv);
++ }
++
++ /* 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",
- time, reftemp, &edat,
++ 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",
+ "<e\\S-\\8D\\4E/kT\\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)
+ {
+ 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)
+ {
+ 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;
+ 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->step - start_step + 1 == edat.nsteps + fr->nsteps)
+ {
+ 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;
+ }
+ }
+ else
+ {
+ /* The interval does not match fr->nsteps:
+ * can not do exact averages.
+ */
+ edat.npoints = 0;
+ }
+ 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)
+ {
+ 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)
+ {
+ 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)
+ {
+ 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;
+}
--- /dev/null
+/*
+ * 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
+ * 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.
+ */
+
+#include "ns.h"
+#include "genborn.h"
+#include "qmmmrec.h"
+#include "idef.h"
+#include "nb_verlet.h"
+#include "interaction_const.h"
+#include "hw_info.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+#if 0
+} /* fixes auto-indentation problems */
+#endif
+
+/* Abstract type for PME that is defined only in the routine that use them. */
+typedef struct gmx_pme *gmx_pme_t;
+
+
+
+/* Structure describing the data in a single table */
+typedef struct
+{
+ enum gmx_table_interaction interaction; /* Types of interactions stored in this table */
+ enum gmx_table_format format; /* Interpolation type and data format */
+
+ real r; /* range of the table */
+ int n; /* n+1 is the number of table points */
+ real scale; /* distance (nm) between two table points */
+ real scale_exp; /* distance for exponential part of VdW table, not always used */
+ real * data; /* the actual table data */
+
+ /* Some information about the table layout. This can also be derived from the interpolation
+ * type and the table interactions, but it is convenient to have here for sanity checks, and it makes it
+ * much easier to access the tables in the nonbonded kernels when we can set the data from variables.
+ * It is always true that stride = formatsize*ninteractions
+ */
+ int formatsize; /* Number of fp variables for each table point (1 for F, 2 for VF, 4 for YFGH, etc.) */
+ int ninteractions; /* Number of interactions in table, 1 for coul-only, 3 for coul+rep+disp. */
+ int stride; /* Distance to next table point (number of fp variables per table point in total) */
+} t_forcetable;
+
+typedef struct
+{
+ t_forcetable table_elec;
+ t_forcetable table_vdw;
+ t_forcetable table_elec_vdw;
+
+ /* The actual neighbor lists, short and long range, see enum above
+ * for definition of neighborlist indices.
+ */
+ t_nblist nlist_sr[eNL_NR];
+ t_nblist nlist_lr[eNL_NR];
+} t_nblists;
+
+/* macros for the cginfo data in forcerec
+ *
+ * Since the tpx format support max 256 energy groups, we do the same here.
+ * Note that we thus have bits 8-14 still unused.
+ *
+ * The maximum cg size in cginfo is 63
+ * because we only have space for 6 bits in cginfo,
+ * this cg size entry is actually only read with domain decomposition.
+ * But there is a smaller limit due to the t_excl data structure
+ * which is defined in nblist.h.
+ */
+#define SET_CGINFO_GID(cgi, gid) (cgi) = (((cgi) & ~255) | (gid))
+#define GET_CGINFO_GID(cgi) ( (cgi) & 255)
+#define SET_CGINFO_FEP(cgi) (cgi) = ((cgi) | (1<<15))
+#define GET_CGINFO_FEP(cgi) ( (cgi) & (1<<15))
+#define SET_CGINFO_EXCL_INTRA(cgi) (cgi) = ((cgi) | (1<<16))
+#define GET_CGINFO_EXCL_INTRA(cgi) ( (cgi) & (1<<16))
+#define SET_CGINFO_EXCL_INTER(cgi) (cgi) = ((cgi) | (1<<17))
+#define GET_CGINFO_EXCL_INTER(cgi) ( (cgi) & (1<<17))
+#define SET_CGINFO_SOLOPT(cgi, opt) (cgi) = (((cgi) & ~(3<<18)) | ((opt)<<18))
+#define GET_CGINFO_SOLOPT(cgi) (((cgi)>>18) & 3)
+#define SET_CGINFO_CONSTR(cgi) (cgi) = ((cgi) | (1<<20))
+#define GET_CGINFO_CONSTR(cgi) ( (cgi) & (1<<20))
+#define SET_CGINFO_SETTLE(cgi) (cgi) = ((cgi) | (1<<21))
+#define GET_CGINFO_SETTLE(cgi) ( (cgi) & (1<<21))
+/* This bit is only used with bBondComm in the domain decomposition */
+#define SET_CGINFO_BOND_INTER(cgi) (cgi) = ((cgi) | (1<<22))
+#define GET_CGINFO_BOND_INTER(cgi) ( (cgi) & (1<<22))
+#define SET_CGINFO_HAS_VDW(cgi) (cgi) = ((cgi) | (1<<23))
+#define GET_CGINFO_HAS_VDW(cgi) ( (cgi) & (1<<23))
+#define SET_CGINFO_HAS_Q(cgi) (cgi) = ((cgi) | (1<<24))
+#define GET_CGINFO_HAS_Q(cgi) ( (cgi) & (1<<24))
+#define SET_CGINFO_NATOMS(cgi, opt) (cgi) = (((cgi) & ~(63<<25)) | ((opt)<<25))
+#define GET_CGINFO_NATOMS(cgi) (((cgi)>>25) & 63)
+
+
+/* Value to be used in mdrun for an infinite cut-off.
+ * Since we need to compare with the cut-off squared,
+ * this value should be slighlty smaller than sqrt(GMX_FLOAT_MAX).
+ */
+#define GMX_CUTOFF_INF 1E+18
+
+/* enums for the neighborlist type */
+enum {
+ enbvdwNONE, enbvdwLJ, enbvdwBHAM, enbvdwTAB, enbvdwNR
+};
+/* OOR is "one over r" -- standard coul */
+enum {
+ enbcoulNONE, enbcoulOOR, enbcoulRF, enbcoulTAB, enbcoulGB, enbcoulFEWALD, enbcoulNR
+};
+
+enum {
+ egCOULSR, egLJSR, egBHAMSR, egCOULLR, egLJLR, egBHAMLR,
+ egCOUL14, egLJ14, egGB, egNR
+};
+
+typedef struct {
+ int nener; /* The number of energy group pairs */
+ real *ener[egNR]; /* Energy terms for each pair of groups */
+} gmx_grppairener_t;
+
+typedef struct {
+ real term[F_NRE]; /* The energies for all different interaction types */
+ gmx_grppairener_t grpp;
+ double dvdl_lin[efptNR]; /* Contributions to dvdl with linear lam-dependence */
+ double dvdl_nonlin[efptNR]; /* Idem, but non-linear dependence */
+ int n_lambda;
+ int fep_state; /*current fep state -- just for printing */
+ double *enerpart_lambda; /* Partial energy for lambda and flambda[] */
+ real foreign_term[F_NRE]; /* alternate array for storing foreign lambda energies */
+ gmx_grppairener_t foreign_grpp; /* alternate array for storing foreign lambda energies */
+} gmx_enerdata_t;
+/* The idea is that dvdl terms with linear lambda dependence will be added
+ * automatically to enerpart_lambda. Terms with non-linear lambda dependence
+ * should explicitly determine the energies at foreign lambda points
+ * when n_lambda > 0.
+ */
+
+typedef struct {
+ int cg_start;
+ int cg_end;
+ int cg_mod;
+ int *cginfo;
+} cginfo_mb_t;
+
+
+/* ewald table type */
+typedef struct ewald_tab *ewald_tab_t;
+
+typedef struct {
+ rvec *f;
+ int f_nalloc;
+ unsigned red_mask; /* Mask for marking which parts of f are filled */
+ rvec *fshift;
+ real ener[F_NRE];
+ gmx_grppairener_t grpp;
+ real Vcorr_q;
+ real Vcorr_lj;
+ real dvdl[efptNR];
+ tensor vir_q;
+ tensor vir_lj;
+} f_thread_t;
+
+typedef struct {
+ interaction_const_t *ic;
+
+ /* Domain Decomposition */
+ gmx_bool bDomDec;
+
+ /* PBC stuff */
+ int ePBC;
+ gmx_bool bMolPBC;
+ int rc_scaling;
+ rvec posres_com;
+ rvec posres_comB;
+
+ const gmx_hw_info_t *hwinfo;
+ const gmx_gpu_opt_t *gpu_opt;
+ gmx_bool use_simd_kernels;
+
+ /* Interaction for calculated in kernels. In many cases this is similar to
+ * the electrostatics settings in the inputrecord, but the difference is that
+ * these variables always specify the actual interaction in the kernel - if
+ * we are tabulating reaction-field the inputrec will say reaction-field, but
+ * the kernel interaction will say cubic-spline-table. To be safe we also
+ * have a kernel-specific setting for the modifiers - if the interaction is
+ * tabulated we already included the inputrec modification there, so the kernel
+ * modification setting will say 'none' in that case.
+ */
+ int nbkernel_elec_interaction;
+ int nbkernel_vdw_interaction;
+ int nbkernel_elec_modifier;
+ int nbkernel_vdw_modifier;
+
+ /* Use special N*N kernels? */
+ gmx_bool bAllvsAll;
+ /* Private work data */
+ void *AllvsAll_work;
+ void *AllvsAll_workgb;
+
+ /* Cut-Off stuff.
+ * Infinite cut-off's will be GMX_CUTOFF_INF (unlike in t_inputrec: 0).
+ */
+ real rlist, rlistlong;
+
+ /* Dielectric constant resp. multiplication factor for charges */
+ real zsquare, temp;
+ real epsilon_r, epsilon_rf, epsfac;
+
+ /* Constants for reaction fields */
+ real kappa, k_rf, c_rf;
+
+ /* Charge sum and dipole for topology A/B ([0]/[1]) for Ewald corrections */
+ double qsum[2];
+ double q2sum[2];
+ double c6sum[2];
+ rvec mu_tot[2];
+
+ /* Dispersion correction stuff */
+ int eDispCorr;
+
+ /* The shift of the shift or user potentials */
+ real enershiftsix;
+ real enershifttwelve;
+ /* Integrated differces for energy and virial with cut-off functions */
+ real enerdiffsix;
+ real enerdifftwelve;
+ real virdiffsix;
+ real virdifftwelve;
+ /* Constant for long range dispersion correction (average dispersion)
+ * for topology A/B ([0]/[1]) */
+ real avcsix[2];
+ /* Constant for long range repulsion term. Relative difference of about
+ * 0.1 percent with 0.8 nm cutoffs. But hey, it's cheap anyway...
+ */
+ real avctwelve[2];
+
+ /* Fudge factors */
+ real fudgeQQ;
+
+ /* Table stuff */
+ gmx_bool bcoultab;
+ gmx_bool bvdwtab;
+ /* The normal tables are in the nblists struct(s) below */
+ t_forcetable tab14; /* for 1-4 interactions only */
+
+ /* PPPM & Shifting stuff */
+ int coulomb_modifier;
+ real rcoulomb_switch, rcoulomb;
+ real *phi;
+
+ /* VdW stuff */
+ int vdw_modifier;
+ double reppow;
+ real rvdw_switch, rvdw;
+ real bham_b_max;
+
+ /* Free energy */
+ int efep;
+ real sc_alphavdw;
+ real sc_alphacoul;
+ int sc_power;
+ real sc_r_power;
+ real sc_sigma6_def;
+ real sc_sigma6_min;
+ gmx_bool bSepDVDL;
+
+ /* NS Stuff */
+ int eeltype;
+ int vdwtype;
+ int cg0, hcg;
+ /* solvent_opt contains the enum for the most common solvent
+ * in the system, which will be optimized.
+ * It can be set to esolNO to disable all water optimization */
+ int solvent_opt;
+ int nWatMol;
+ gmx_bool bGrid;
+ gmx_bool bExcl_IntraCGAll_InterCGNone;
+ cginfo_mb_t *cginfo_mb;
+ int *cginfo;
+ rvec *cg_cm;
+ int cg_nalloc;
+ rvec *shift_vec;
+
+ /* The neighborlists including tables */
+ int nnblists;
+ int *gid2nblists;
+ t_nblists *nblists;
+
+ int cutoff_scheme; /* group- or Verlet-style cutoff */
+ gmx_bool bNonbonded; /* true if nonbonded calculations are *not* turned off */
+ nonbonded_verlet_t *nbv;
+
+ /* The wall tables (if used) */
+ int nwall;
+ t_forcetable **wall_tab;
+
+ /* The number of charge groups participating in do_force_lowlevel */
+ int ncg_force;
+ /* The number of atoms participating in do_force_lowlevel */
+ int natoms_force;
+ /* The number of atoms participating in force and constraints */
+ int natoms_force_constr;
+ /* The allocation size of vectors of size natoms_force */
+ int nalloc_force;
+
+ /* Twin Range stuff, f_twin has size natoms_force */
+ gmx_bool bTwinRange;
+ int nlr;
+ rvec *f_twin;
++ /* Constraint virial correction for multiple time stepping */
++ tensor vir_twin_constr;
+
+ /* Forces that should not enter into the virial summation:
+ * PPPM/PME/Ewald/posres
+ */
+ gmx_bool bF_NoVirSum;
+ int f_novirsum_n;
+ int f_novirsum_nalloc;
+ rvec *f_novirsum_alloc;
+ /* Pointer that points to f_novirsum_alloc when pressure is calcaluted,
+ * points to the normal force vectors wen pressure is not requested.
+ */
+ rvec *f_novirsum;
+
+ /* Long-range forces and virial for PPPM/PME/Ewald */
+ gmx_pme_t pmedata;
+ int ljpme_combination_rule;
+ tensor vir_el_recip;
+ tensor vir_lj_recip;
+
+ /* PME/Ewald stuff */
+ gmx_bool bEwald;
+ real ewaldcoeff_q;
+ real ewaldcoeff_lj;
+ ewald_tab_t ewald_table;
+
+ /* Virial Stuff */
+ rvec *fshift;
+ rvec vir_diag_posres;
+ dvec vir_wall_z;
+
+ /* Non bonded Parameter lists */
+ int ntype; /* Number of atom types */
+ gmx_bool bBHAM;
+ real *nbfp;
+ real *ljpme_c6grid; /* C6-values used on grid in LJPME */
+
+ /* Energy group pair flags */
+ int *egp_flags;
+
+ /* Shell molecular dynamics flexible constraints */
+ real fc_stepsize;
+
+ /* Generalized born implicit solvent */
+ gmx_bool bGB;
+ /* Generalized born stuff */
+ real gb_epsilon_solvent;
+ /* Table data for GB */
+ t_forcetable gbtab;
+ /* VdW radius for each atomtype (dim is thus ntype) */
+ real *atype_radius;
+ /* Effective radius (derived from effective volume) for each type */
+ real *atype_vol;
+ /* Implicit solvent - surface tension for each atomtype */
+ real *atype_surftens;
+ /* Implicit solvent - radius for GB calculation */
+ real *atype_gb_radius;
+ /* Implicit solvent - overlap for HCT model */
+ real *atype_S_hct;
+ /* Generalized born interaction data */
+ gmx_genborn_t *born;
+
+ /* Table scale for GB */
+ real gbtabscale;
+ /* Table range for GB */
+ real gbtabr;
+ /* GB neighborlists (the sr list will contain for each atom all other atoms
+ * (for use in the SA calculation) and the lr list will contain
+ * for each atom all atoms 1-4 or greater (for use in the GB calculation)
+ */
+ t_nblist gblist_sr;
+ t_nblist gblist_lr;
+ t_nblist gblist;
+
+ /* Inverse square root of the Born radii for implicit solvent */
+ real *invsqrta;
+ /* Derivatives of the potential with respect to the Born radii */
+ real *dvda;
+ /* Derivatives of the Born radii with respect to coordinates */
+ real *dadx;
+ real *dadx_rawptr;
+ int nalloc_dadx; /* Allocated size of dadx */
+
+ /* If > 0 signals Test Particle Insertion,
+ * the value is the number of atoms of the molecule to insert
+ * Only the energy difference due to the addition of the last molecule
+ * should be calculated.
+ */
+ gmx_bool n_tpi;
+
+ /* Neighbor searching stuff */
+ gmx_ns_t ns;
+
+ /* QMMM stuff */
+ gmx_bool bQMMM;
+ t_QMMMrec *qr;
+
+ /* QM-MM neighborlists */
+ t_nblist QMMMlist;
+
+ /* Limit for printing large forces, negative is don't print */
+ real print_force;
+
+ /* coarse load balancing time measurement */
+ double t_fnbf;
+ double t_wait;
+ int timesteps;
+
+ /* parameter needed for AdResS simulation */
+ int adress_type;
+ gmx_bool badress_tf_full_box;
+ real adress_const_wf;
+ real adress_ex_width;
+ real adress_hy_width;
+ int adress_icor;
+ int adress_site;
+ rvec adress_refs;
+ int n_adress_tf_grps;
+ int * adress_tf_table_index;
+ int *adress_group_explicit;
+ t_forcetable * atf_tabs;
+ real adress_ex_forcecap;
+ gmx_bool adress_do_hybridpairs;
+
+ /* User determined parameters, copied from the inputrec */
+ int userint1;
+ int userint2;
+ int userint3;
+ int userint4;
+ real userreal1;
+ real userreal2;
+ real userreal3;
+ real userreal4;
+
+ /* Thread local force and energy data */
+ /* FIXME move to bonded_thread_data_t */
+ int nthreads;
+ int red_ashift;
+ int red_nblock;
+ f_thread_t *f_t;
+
+ /* Exclusion load distribution over the threads */
+ int *excl_load;
+} t_forcerec;
+
+/* Important: Starting with Gromacs-4.6, the values of c6 and c12 in the nbfp array have
+ * been scaled by 6.0 or 12.0 to save flops in the kernels. We have corrected this everywhere
+ * in the code, but beware if you are using these macros externally.
+ */
+#define C6(nbfp, ntp, ai, aj) (nbfp)[2*((ntp)*(ai)+(aj))]
+#define C12(nbfp, ntp, ai, aj) (nbfp)[2*((ntp)*(ai)+(aj))+1]
+#define BHAMC(nbfp, ntp, ai, aj) (nbfp)[3*((ntp)*(ai)+(aj))]
+#define BHAMA(nbfp, ntp, ai, aj) (nbfp)[3*((ntp)*(ai)+(aj))+1]
+#define BHAMB(nbfp, ntp, ai, aj) (nbfp)[3*((ntp)*(ai)+(aj))+2]
+
+#ifdef __cplusplus
+}
+#endif
--- /dev/null
+/*
+ * 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
+ * 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.
+ */
+
+#ifndef _update_h
+#define _update_h
+
+#include "typedefs.h"
+#include "mshift.h"
+#include "tgroup.h"
+#include "network.h"
+#include "vec.h"
+
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/* Abstract type for stochastic dynamics */
+typedef struct gmx_update *gmx_update_t;
+
+/* Initialize the stochastic dynamics struct */
+gmx_update_t init_update(t_inputrec *ir);
+
+/* Store the random state from sd in state */
+void get_stochd_state(gmx_update_t sd, t_state *state);
+
+/* Set the random in sd from state */
+void set_stochd_state(gmx_update_t sd, t_state *state);
+
+/* Store the box at step step
+ * as a reference state for simulations with box deformation.
+ */
+void set_deform_reference_box(gmx_update_t upd,
+ gmx_int64_t step, matrix box);
+
+void update_tcouple(gmx_int64_t step,
+ t_inputrec *inputrec,
+ t_state *state,
+ gmx_ekindata_t *ekind,
+ t_extmass *MassQ,
+ t_mdatoms *md
+ );
+
+void update_pcouple(FILE *fplog,
+ gmx_int64_t step,
+ t_inputrec *inputrec,
+ t_state *state,
+ matrix pcoupl_mu,
+ matrix M,
+ gmx_bool bInitStep);
+
+void update_coords(FILE *fplog,
+ gmx_int64_t step,
+ t_inputrec *inputrec, /* input record and box stuff */
+ t_mdatoms *md,
+ t_state *state,
+ gmx_bool bMolPBC,
+ rvec *f, /* forces on home particles */
+ gmx_bool bDoLR,
+ rvec *f_lr,
++ tensor *vir_lr_constr,
+ t_fcdata *fcd,
+ gmx_ekindata_t *ekind,
+ matrix M,
+ gmx_update_t upd,
+ gmx_bool bInitStep,
+ int bUpdatePart,
+ t_commrec *cr, /* these shouldn't be here -- need to think about it */
+ t_nrnb *nrnb,
+ gmx_constr_t constr,
+ t_idef *idef);
+
+/* Return TRUE if OK, FALSE in case of Shake Error */
+
+extern gmx_bool update_randomize_velocities(t_inputrec *ir, gmx_int64_t step, const t_commrec *cr, t_mdatoms *md, t_state *state, gmx_update_t upd, gmx_constr_t constr);
+
+void update_constraints(FILE *fplog,
+ gmx_int64_t step,
+ real *dvdlambda, /* FEP stuff */
+ t_inputrec *inputrec, /* input record and box stuff */
+ gmx_ekindata_t *ekind,
+ t_mdatoms *md,
+ t_state *state,
+ gmx_bool bMolPBC,
+ t_graph *graph,
+ rvec force[], /* forces on home particles */
+ t_idef *idef,
+ tensor vir_part,
+ t_commrec *cr,
+ t_nrnb *nrnb,
+ gmx_wallcycle_t wcycle,
+ gmx_update_t upd,
+ gmx_constr_t constr,
+ gmx_bool bFirstHalf,
+ gmx_bool bCalcVir,
+ real vetanew);
+
+/* Return TRUE if OK, FALSE in case of Shake Error */
+
+void update_box(FILE *fplog,
+ gmx_int64_t step,
+ t_inputrec *inputrec, /* input record and box stuff */
+ t_mdatoms *md,
+ t_state *state,
+ rvec force[], /* forces on home particles */
+ matrix *scale_tot,
+ matrix pcoupl_mu,
+ t_nrnb *nrnb,
+ gmx_update_t upd);
+/* Return TRUE if OK, FALSE in case of Shake Error */
+
+void calc_ke_part(t_state *state, t_grpopts *opts, t_mdatoms *md,
+ gmx_ekindata_t *ekind, t_nrnb *nrnb, gmx_bool bEkinAveVel, gmx_bool bSaveOld);
+/*
+ * Compute the partial kinetic energy for home particles;
+ * will be accumulated in the calling routine.
+ * The tensor is
+ *
+ * Ekin = SUM(i) 0.5 m[i] v[i] (x) v[i]
+ *
+ * use v[i] = v[i] - u[i] when calculating temperature
+ *
+ * u must be accumulated already.
+ *
+ * Now also computes the contribution of the kinetic energy to the
+ * free energy
+ *
+ */
+
+
+void
+init_ekinstate(ekinstate_t *ekinstate, const t_inputrec *ir);
+
+void
+update_ekinstate(ekinstate_t *ekinstate, gmx_ekindata_t *ekind);
+
+void
+restore_ekinstate_from_state(t_commrec *cr,
+ gmx_ekindata_t *ekind, ekinstate_t *ekinstate);
+
+void berendsen_tcoupl(t_inputrec *ir, gmx_ekindata_t *ekind, real dt);
+
+void andersen_tcoupl(t_inputrec *ir, gmx_int64_t step,
+ const t_commrec *cr, const t_mdatoms *md, t_state *state, real rate, const gmx_bool *randomize, const real *boltzfac);
+
+void nosehoover_tcoupl(t_grpopts *opts, gmx_ekindata_t *ekind, real dt,
+ double xi[], double vxi[], t_extmass *MassQ);
+
+t_state *init_bufstate(const t_state *template_state);
+
+void destroy_bufstate(t_state *state);
+
+void trotter_update(t_inputrec *ir, gmx_int64_t step, gmx_ekindata_t *ekind,
+ gmx_enerdata_t *enerd, t_state *state, tensor vir, t_mdatoms *md,
+ t_extmass *MassQ, int **trotter_seqlist, int trotter_seqno);
+
+int **init_npt_vars(t_inputrec *ir, t_state *state, t_extmass *Mass, gmx_bool bTrotter);
+
+real NPT_energy(t_inputrec *ir, t_state *state, t_extmass *MassQ);
+/* computes all the pressure/tempertature control energy terms to get a conserved energy */
+
+void NBaroT_trotter(t_grpopts *opts, real dt,
+ double xi[], double vxi[], real *veta, t_extmass *MassQ);
+
+void vrescale_tcoupl(t_inputrec *ir, gmx_int64_t step,
+ gmx_ekindata_t *ekind, real dt,
+ double therm_integral[]);
+/* Compute temperature scaling. For V-rescale it is done in update. */
+
+real vrescale_energy(t_grpopts *opts, double therm_integral[]);
+/* Returns the V-rescale contribution to the conserved energy */
+
+void rescale_velocities(gmx_ekindata_t *ekind, t_mdatoms *mdatoms,
+ int start, int end, rvec v[]);
+/* Rescale the velocities with the scaling factor in ekind */
+
+void update_annealing_target_temp(t_grpopts *opts, real t);
+/* Set reference temp for simulated annealing at time t*/
+
+real calc_temp(real ekin, real nrdf);
+/* Calculate the temperature */
+
+real calc_pres(int ePBC, int nwall, matrix box, tensor ekin, tensor vir,
+ tensor pres);
+/* Calculate the pressure tensor, returns the scalar pressure.
+ * The unit of pressure is bar.
+ */
+
+void parrinellorahman_pcoupl(FILE *fplog, gmx_int64_t step,
+ t_inputrec *ir, real dt, tensor pres,
+ tensor box, tensor box_rel, tensor boxv,
+ tensor M, matrix mu,
+ gmx_bool bFirstStep);
+
+void berendsen_pcoupl(FILE *fplog, gmx_int64_t step,
+ t_inputrec *ir, real dt, tensor pres, matrix box,
+ matrix mu);
+
+
+void berendsen_pscale(t_inputrec *ir, matrix mu,
+ matrix box, matrix box_rel,
+ int start, int nr_atoms,
+ rvec x[], unsigned short cFREEZE[],
+ t_nrnb *nrnb);
+
+void correct_ekin(FILE *log, int start, int end, rvec v[],
+ rvec vcm, real mass[], real tmass, tensor ekin);
+/* Correct ekin for vcm */
+
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif /* _update_h */
--- /dev/null
- static void combine_forces(int nstcalclr,
- gmx_constr_t constr,
- t_inputrec *ir, t_mdatoms *md, t_idef *idef,
- t_commrec *cr,
- gmx_int64_t step,
- t_state *state, gmx_bool bMolPBC,
- int start, int nrend,
- rvec f[], rvec f_lr[],
- t_nrnb *nrnb)
- {
- int i, d, nm1;
-
- /* f contains the short-range forces + the long range forces
- * which are stored separately in f_lr.
- */
-
- if (constr != NULL && !(ir->eConstrAlg == econtSHAKE && ir->epc == epcNO))
- {
- /* We need to constrain the LR forces separately,
- * because due to the different pre-factor for the SR and LR
- * forces in the update algorithm, we can not determine
- * the constraint force for the coordinate constraining.
- * Constrain only the additional LR part of the force.
- */
- /* MRS -- need to make sure this works with trotter integration -- the constraint calls may not be right.*/
- constrain(NULL, FALSE, FALSE, constr, idef, ir, NULL, cr, step, 0, md,
- state->x, f_lr, f_lr, bMolPBC, state->box, state->lambda[efptBONDED], NULL,
- NULL, NULL, nrnb, econqForce, ir->epc == epcMTTK, state->veta, state->veta);
- }
-
- /* Add nstcalclr-1 times the LR force to the sum of both forces
- * and store the result in forces_lr.
- */
- nm1 = nstcalclr - 1;
- for (i = start; i < nrend; i++)
- {
- for (d = 0; d < DIM; d++)
- {
- f_lr[i][d] = f[i][d] + nm1*f_lr[i][d];
- }
- }
- }
-
+/*
+ * 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
+ * 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 <config.h>
+#endif
+
+
+#include <stdio.h>
+#include <math.h>
+
+#include "types/commrec.h"
+#include "sysstuff.h"
+#include "gromacs/utility/smalloc.h"
+#include "typedefs.h"
+#include "nrnb.h"
+#include "physics.h"
+#include "macros.h"
+#include "vec.h"
+#include "main.h"
+#include "update.h"
+#include "gromacs/random/random.h"
+#include "mshift.h"
+#include "tgroup.h"
+#include "force.h"
+#include "names.h"
+#include "txtdump.h"
+#include "mdrun.h"
+#include "constr.h"
+#include "disre.h"
+#include "orires.h"
+#include "gmx_omp_nthreads.h"
+
+#include "gromacs/fileio/confio.h"
+#include "gromacs/fileio/futil.h"
+#include "gromacs/timing/wallcycle.h"
+#include "gromacs/utility/gmxomp.h"
+#include "gromacs/pulling/pull.h"
+
+/*For debugging, start at v(-dt/2) for velolcity verlet -- uncomment next line */
+/*#define STARTFROMDT2*/
+
+typedef struct {
+ double gdt;
+ double eph;
+ double emh;
+ double em;
+ double b;
+ double c;
+ double d;
+} gmx_sd_const_t;
+
+typedef struct {
+ real V;
+ real X;
+ real Yv;
+ real Yx;
+} gmx_sd_sigma_t;
+
+typedef struct {
+ /* BD stuff */
+ real *bd_rf;
+ /* SD stuff */
+ gmx_sd_const_t *sdc;
+ gmx_sd_sigma_t *sdsig;
+ rvec *sd_V;
+ int sd_V_nalloc;
+ /* andersen temperature control stuff */
+ gmx_bool *randomize_group;
+ real *boltzfac;
+} gmx_stochd_t;
+
+typedef struct gmx_update
+{
+ gmx_stochd_t *sd;
+ /* xprime for constraint algorithms */
+ rvec *xp;
+ int xp_nalloc;
+
+ /* Variables for the deform algorithm */
+ gmx_int64_t deformref_step;
+ matrix deformref_box;
+} t_gmx_update;
+
+
+static void do_update_md(int start, int nrend, double dt,
+ t_grp_tcstat *tcstat,
+ double nh_vxi[],
+ gmx_bool bNEMD, t_grp_acc *gstat, rvec accel[],
+ ivec nFreeze[],
+ real invmass[],
+ unsigned short ptype[], unsigned short cFREEZE[],
+ unsigned short cACC[], unsigned short cTC[],
+ rvec x[], rvec xprime[], rvec v[],
+ rvec f[], matrix M,
+ gmx_bool bNH, gmx_bool bPR)
+{
+ double imass, w_dt;
+ int gf = 0, ga = 0, gt = 0;
+ rvec vrel;
+ real vn, vv, va, vb, vnrel;
+ real lg, vxi = 0, u;
+ int n, d;
+
+ if (bNH || bPR)
+ {
+ /* Update with coupling to extended ensembles, used for
+ * Nose-Hoover and Parrinello-Rahman coupling
+ * Nose-Hoover uses the reversible leap-frog integrator from
+ * Holian et al. Phys Rev E 52(3) : 2338, 1995
+ */
+ for (n = start; n < nrend; n++)
+ {
+ imass = invmass[n];
+ if (cFREEZE)
+ {
+ gf = cFREEZE[n];
+ }
+ if (cACC)
+ {
+ ga = cACC[n];
+ }
+ if (cTC)
+ {
+ gt = cTC[n];
+ }
+ lg = tcstat[gt].lambda;
+ if (bNH)
+ {
+ vxi = nh_vxi[gt];
+ }
+ rvec_sub(v[n], gstat[ga].u, vrel);
+
+ for (d = 0; d < DIM; d++)
+ {
+ if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
+ {
+ vnrel = (lg*vrel[d] + dt*(imass*f[n][d] - 0.5*vxi*vrel[d]
+ - iprod(M[d], vrel)))/(1 + 0.5*vxi*dt);
+ /* do not scale the mean velocities u */
+ vn = gstat[ga].u[d] + accel[ga][d]*dt + vnrel;
+ v[n][d] = vn;
+ xprime[n][d] = x[n][d]+vn*dt;
+ }
+ else
+ {
+ v[n][d] = 0.0;
+ xprime[n][d] = x[n][d];
+ }
+ }
+ }
+ }
+ else if (cFREEZE != NULL ||
+ nFreeze[0][XX] || nFreeze[0][YY] || nFreeze[0][ZZ] ||
+ bNEMD)
+ {
+ /* Update with Berendsen/v-rescale coupling and freeze or NEMD */
+ for (n = start; n < nrend; n++)
+ {
+ w_dt = invmass[n]*dt;
+ if (cFREEZE)
+ {
+ gf = cFREEZE[n];
+ }
+ if (cACC)
+ {
+ ga = cACC[n];
+ }
+ if (cTC)
+ {
+ gt = cTC[n];
+ }
+ lg = tcstat[gt].lambda;
+
+ for (d = 0; d < DIM; d++)
+ {
+ vn = v[n][d];
+ if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
+ {
+ vv = lg*vn + f[n][d]*w_dt;
+
+ /* do not scale the mean velocities u */
+ u = gstat[ga].u[d];
+ va = vv + accel[ga][d]*dt;
+ vb = va + (1.0-lg)*u;
+ v[n][d] = vb;
+ xprime[n][d] = x[n][d]+vb*dt;
+ }
+ else
+ {
+ v[n][d] = 0.0;
+ xprime[n][d] = x[n][d];
+ }
+ }
+ }
+ }
+ else
+ {
+ /* Plain update with Berendsen/v-rescale coupling */
+ for (n = start; n < nrend; n++)
+ {
+ if ((ptype[n] != eptVSite) && (ptype[n] != eptShell))
+ {
+ w_dt = invmass[n]*dt;
+ if (cTC)
+ {
+ gt = cTC[n];
+ }
+ lg = tcstat[gt].lambda;
+
+ for (d = 0; d < DIM; d++)
+ {
+ vn = lg*v[n][d] + f[n][d]*w_dt;
+ v[n][d] = vn;
+ xprime[n][d] = x[n][d] + vn*dt;
+ }
+ }
+ else
+ {
+ for (d = 0; d < DIM; d++)
+ {
+ v[n][d] = 0.0;
+ xprime[n][d] = x[n][d];
+ }
+ }
+ }
+ }
+}
+
+static void do_update_vv_vel(int start, int nrend, double dt,
+ rvec accel[], ivec nFreeze[], real invmass[],
+ unsigned short ptype[], unsigned short cFREEZE[],
+ unsigned short cACC[], rvec v[], rvec f[],
+ gmx_bool bExtended, real veta, real alpha)
+{
+ double imass, w_dt;
+ int gf = 0, ga = 0;
+ rvec vrel;
+ real u, vn, vv, va, vb, vnrel;
+ int n, d;
+ double g, mv1, mv2;
+
+ if (bExtended)
+ {
+ g = 0.25*dt*veta*alpha;
+ mv1 = exp(-g);
+ mv2 = series_sinhx(g);
+ }
+ else
+ {
+ mv1 = 1.0;
+ mv2 = 1.0;
+ }
+ for (n = start; n < nrend; n++)
+ {
+ w_dt = invmass[n]*dt;
+ if (cFREEZE)
+ {
+ gf = cFREEZE[n];
+ }
+ if (cACC)
+ {
+ ga = cACC[n];
+ }
+
+ for (d = 0; d < DIM; d++)
+ {
+ if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
+ {
+ v[n][d] = mv1*(mv1*v[n][d] + 0.5*(w_dt*mv2*f[n][d]))+0.5*accel[ga][d]*dt;
+ }
+ else
+ {
+ v[n][d] = 0.0;
+ }
+ }
+ }
+} /* do_update_vv_vel */
+
+static void do_update_vv_pos(int start, int nrend, double dt,
+ ivec nFreeze[],
+ unsigned short ptype[], unsigned short cFREEZE[],
+ rvec x[], rvec xprime[], rvec v[],
+ gmx_bool bExtended, real veta)
+{
+ double imass, w_dt;
+ int gf = 0;
+ int n, d;
+ double g, mr1, mr2;
+
+ /* Would it make more sense if Parrinello-Rahman was put here? */
+ if (bExtended)
+ {
+ g = 0.5*dt*veta;
+ mr1 = exp(g);
+ mr2 = series_sinhx(g);
+ }
+ else
+ {
+ mr1 = 1.0;
+ mr2 = 1.0;
+ }
+
+ for (n = start; n < nrend; n++)
+ {
+
+ if (cFREEZE)
+ {
+ gf = cFREEZE[n];
+ }
+
+ for (d = 0; d < DIM; d++)
+ {
+ if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
+ {
+ xprime[n][d] = mr1*(mr1*x[n][d]+mr2*dt*v[n][d]);
+ }
+ else
+ {
+ xprime[n][d] = x[n][d];
+ }
+ }
+ }
+} /* do_update_vv_pos */
+
+static void do_update_visc(int start, int nrend, double dt,
+ t_grp_tcstat *tcstat,
+ double nh_vxi[],
+ real invmass[],
+ unsigned short ptype[], unsigned short cTC[],
+ rvec x[], rvec xprime[], rvec v[],
+ rvec f[], matrix M, matrix box, real
+ cos_accel, real vcos,
+ gmx_bool bNH, gmx_bool bPR)
+{
+ double imass, w_dt;
+ int gt = 0;
+ real vn, vc;
+ real lg, vxi = 0, vv;
+ real fac, cosz;
+ rvec vrel;
+ int n, d;
+
+ fac = 2*M_PI/(box[ZZ][ZZ]);
+
+ if (bNH || bPR)
+ {
+ /* Update with coupling to extended ensembles, used for
+ * Nose-Hoover and Parrinello-Rahman coupling
+ */
+ for (n = start; n < nrend; n++)
+ {
+ imass = invmass[n];
+ if (cTC)
+ {
+ gt = cTC[n];
+ }
+ lg = tcstat[gt].lambda;
+ cosz = cos(fac*x[n][ZZ]);
+
+ copy_rvec(v[n], vrel);
+
+ vc = cosz*vcos;
+ vrel[XX] -= vc;
+ if (bNH)
+ {
+ vxi = nh_vxi[gt];
+ }
+ for (d = 0; d < DIM; d++)
+ {
+ vn = v[n][d];
+
+ if ((ptype[n] != eptVSite) && (ptype[n] != eptShell))
+ {
+ vn = (lg*vrel[d] + dt*(imass*f[n][d] - 0.5*vxi*vrel[d]
+ - iprod(M[d], vrel)))/(1 + 0.5*vxi*dt);
+ if (d == XX)
+ {
+ vn += vc + dt*cosz*cos_accel;
+ }
+ v[n][d] = vn;
+ xprime[n][d] = x[n][d]+vn*dt;
+ }
+ else
+ {
+ xprime[n][d] = x[n][d];
+ }
+ }
+ }
+ }
+ else
+ {
+ /* Classic version of update, used with berendsen coupling */
+ for (n = start; n < nrend; n++)
+ {
+ w_dt = invmass[n]*dt;
+ if (cTC)
+ {
+ gt = cTC[n];
+ }
+ lg = tcstat[gt].lambda;
+ cosz = cos(fac*x[n][ZZ]);
+
+ for (d = 0; d < DIM; d++)
+ {
+ vn = v[n][d];
+
+ if ((ptype[n] != eptVSite) && (ptype[n] != eptShell))
+ {
+ if (d == XX)
+ {
+ vc = cosz*vcos;
+ /* Do not scale the cosine velocity profile */
+ vv = vc + lg*(vn - vc + f[n][d]*w_dt);
+ /* Add the cosine accelaration profile */
+ vv += dt*cosz*cos_accel;
+ }
+ else
+ {
+ vv = lg*(vn + f[n][d]*w_dt);
+ }
+ v[n][d] = vv;
+ xprime[n][d] = x[n][d]+vv*dt;
+ }
+ else
+ {
+ v[n][d] = 0.0;
+ xprime[n][d] = x[n][d];
+ }
+ }
+ }
+ }
+}
+
+static gmx_stochd_t *init_stochd(t_inputrec *ir)
+{
+ gmx_stochd_t *sd;
+ gmx_sd_const_t *sdc;
+ int ngtc, n, th;
+ real y;
+
+ snew(sd, 1);
+
+ ngtc = ir->opts.ngtc;
+
+ if (ir->eI == eiBD)
+ {
+ snew(sd->bd_rf, ngtc);
+ }
+ else if (EI_SD(ir->eI))
+ {
+ snew(sd->sdc, ngtc);
+ snew(sd->sdsig, ngtc);
+
+ sdc = sd->sdc;
+ for (n = 0; n < ngtc; n++)
+ {
+ if (ir->opts.tau_t[n] > 0)
+ {
+ sdc[n].gdt = ir->delta_t/ir->opts.tau_t[n];
+ sdc[n].eph = exp(sdc[n].gdt/2);
+ sdc[n].emh = exp(-sdc[n].gdt/2);
+ sdc[n].em = exp(-sdc[n].gdt);
+ }
+ else
+ {
+ /* No friction and noise on this group */
+ sdc[n].gdt = 0;
+ sdc[n].eph = 1;
+ sdc[n].emh = 1;
+ sdc[n].em = 1;
+ }
+ if (sdc[n].gdt >= 0.05)
+ {
+ sdc[n].b = sdc[n].gdt*(sdc[n].eph*sdc[n].eph - 1)
+ - 4*(sdc[n].eph - 1)*(sdc[n].eph - 1);
+ sdc[n].c = sdc[n].gdt - 3 + 4*sdc[n].emh - sdc[n].em;
+ sdc[n].d = 2 - sdc[n].eph - sdc[n].emh;
+ }
+ else
+ {
+ y = sdc[n].gdt/2;
+ /* Seventh order expansions for small y */
+ sdc[n].b = y*y*y*y*(1/3.0+y*(1/3.0+y*(17/90.0+y*7/9.0)));
+ sdc[n].c = y*y*y*(2/3.0+y*(-1/2.0+y*(7/30.0+y*(-1/12.0+y*31/1260.0))));
+ sdc[n].d = y*y*(-1+y*y*(-1/12.0-y*y/360.0));
+ }
+ if (debug)
+ {
+ fprintf(debug, "SD const tc-grp %d: b %g c %g d %g\n",
+ n, sdc[n].b, sdc[n].c, sdc[n].d);
+ }
+ }
+ }
+ else if (ETC_ANDERSEN(ir->etc))
+ {
+ int ngtc;
+ t_grpopts *opts;
+ real reft;
+
+ opts = &ir->opts;
+ ngtc = opts->ngtc;
+
+ snew(sd->randomize_group, ngtc);
+ snew(sd->boltzfac, ngtc);
+
+ /* for now, assume that all groups, if randomized, are randomized at the same rate, i.e. tau_t is the same. */
+ /* since constraint groups don't necessarily match up with temperature groups! This is checked in readir.c */
+
+ for (n = 0; n < ngtc; n++)
+ {
+ reft = max(0.0, opts->ref_t[n]);
+ if ((opts->tau_t[n] > 0) && (reft > 0)) /* tau_t or ref_t = 0 means that no randomization is done */
+ {
+ sd->randomize_group[n] = TRUE;
+ sd->boltzfac[n] = BOLTZ*opts->ref_t[n];
+ }
+ else
+ {
+ sd->randomize_group[n] = FALSE;
+ }
+ }
+ }
+ return sd;
+}
+
+gmx_update_t init_update(t_inputrec *ir)
+{
+ t_gmx_update *upd;
+
+ snew(upd, 1);
+
+ if (ir->eI == eiBD || EI_SD(ir->eI) || ir->etc == etcVRESCALE || ETC_ANDERSEN(ir->etc))
+ {
+ upd->sd = init_stochd(ir);
+ }
+
+ upd->xp = NULL;
+ upd->xp_nalloc = 0;
+
+ return upd;
+}
+
+static void do_update_sd1(gmx_stochd_t *sd,
+ int start, int nrend, double dt,
+ rvec accel[], ivec nFreeze[],
+ real invmass[], unsigned short ptype[],
+ unsigned short cFREEZE[], unsigned short cACC[],
+ unsigned short cTC[],
+ rvec x[], rvec xprime[], rvec v[], rvec f[],
+ int ngtc, real tau_t[], real ref_t[],
+ gmx_int64_t step, int seed, int* gatindex)
+{
+ gmx_sd_const_t *sdc;
+ gmx_sd_sigma_t *sig;
+ real kT;
+ int gf = 0, ga = 0, gt = 0;
+ real ism, sd_V;
+ int n, d;
+
+ sdc = sd->sdc;
+ sig = sd->sdsig;
+
+ for (n = 0; n < ngtc; n++)
+ {
+ kT = BOLTZ*ref_t[n];
+ /* The mass is encounted for later, since this differs per atom */
+ sig[n].V = sqrt(kT*(1 - sdc[n].em*sdc[n].em));
+ }
+
+ for (n = start; n < nrend; n++)
+ {
+ real rnd[3];
+ int ng = gatindex ? gatindex[n] : n;
+ ism = sqrt(invmass[n]);
+ if (cFREEZE)
+ {
+ gf = cFREEZE[n];
+ }
+ if (cACC)
+ {
+ ga = cACC[n];
+ }
+ if (cTC)
+ {
+ gt = cTC[n];
+ }
+
+ gmx_rng_cycle_3gaussian_table(step, ng, seed, RND_SEED_UPDATE, rnd);
+ for (d = 0; d < DIM; d++)
+ {
+ if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
+ {
+ sd_V = ism*sig[gt].V*rnd[d];
+
+ v[n][d] = v[n][d]*sdc[gt].em
+ + (invmass[n]*f[n][d] + accel[ga][d])*tau_t[gt]*(1 - sdc[gt].em)
+ + sd_V;
+
+ xprime[n][d] = x[n][d] + v[n][d]*dt;
+ }
+ else
+ {
+ v[n][d] = 0.0;
+ xprime[n][d] = x[n][d];
+ }
+ }
+ }
+}
+
+static void check_sd2_work_data_allocation(gmx_stochd_t *sd, int nrend)
+{
+ if (nrend > sd->sd_V_nalloc)
+ {
+ sd->sd_V_nalloc = over_alloc_dd(nrend);
+ srenew(sd->sd_V, sd->sd_V_nalloc);
+ }
+}
+
+static void do_update_sd2_Tconsts(gmx_stochd_t *sd,
+ int ngtc,
+ const real tau_t[],
+ const real ref_t[])
+{
+ /* This is separated from the update below, because it is single threaded */
+ gmx_sd_const_t *sdc;
+ gmx_sd_sigma_t *sig;
+ int gt;
+ real kT;
+
+ sdc = sd->sdc;
+ sig = sd->sdsig;
+
+ for (gt = 0; gt < ngtc; gt++)
+ {
+ kT = BOLTZ*ref_t[gt];
+ /* The mass is encounted for later, since this differs per atom */
+ sig[gt].V = sqrt(kT*(1-sdc[gt].em));
+ sig[gt].X = sqrt(kT*sqr(tau_t[gt])*sdc[gt].c);
+ sig[gt].Yv = sqrt(kT*sdc[gt].b/sdc[gt].c);
+ sig[gt].Yx = sqrt(kT*sqr(tau_t[gt])*sdc[gt].b/(1-sdc[gt].em));
+ }
+}
+
+static void do_update_sd2(gmx_stochd_t *sd,
+ gmx_bool bInitStep,
+ int start, int nrend,
+ rvec accel[], ivec nFreeze[],
+ real invmass[], unsigned short ptype[],
+ unsigned short cFREEZE[], unsigned short cACC[],
+ unsigned short cTC[],
+ rvec x[], rvec xprime[], rvec v[], rvec f[],
+ rvec sd_X[],
+ const real tau_t[],
+ gmx_bool bFirstHalf, gmx_int64_t step, int seed,
+ int* gatindex)
+{
+ gmx_sd_const_t *sdc;
+ gmx_sd_sigma_t *sig;
+ /* The random part of the velocity update, generated in the first
+ * half of the update, needs to be remembered for the second half.
+ */
+ rvec *sd_V;
+ real kT;
+ int gf = 0, ga = 0, gt = 0;
+ real vn = 0, Vmh, Xmh;
+ real ism;
+ int n, d, ng;
+
+ sdc = sd->sdc;
+ sig = sd->sdsig;
+ sd_V = sd->sd_V;
+
+ for (n = start; n < nrend; n++)
+ {
+ real rnd[6], rndi[3];
+ ng = gatindex ? gatindex[n] : n;
+ ism = sqrt(invmass[n]);
+ if (cFREEZE)
+ {
+ gf = cFREEZE[n];
+ }
+ if (cACC)
+ {
+ ga = cACC[n];
+ }
+ if (cTC)
+ {
+ gt = cTC[n];
+ }
+
+ gmx_rng_cycle_6gaussian_table(step*2+(bFirstHalf ? 1 : 2), ng, seed, RND_SEED_UPDATE, rnd);
+ if (bInitStep)
+ {
+ gmx_rng_cycle_3gaussian_table(step*2, ng, seed, RND_SEED_UPDATE, rndi);
+ }
+ for (d = 0; d < DIM; d++)
+ {
+ if (bFirstHalf)
+ {
+ vn = v[n][d];
+ }
+ if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
+ {
+ if (bFirstHalf)
+ {
+ if (bInitStep)
+ {
+ sd_X[n][d] = ism*sig[gt].X*rndi[d];
+ }
+ Vmh = sd_X[n][d]*sdc[gt].d/(tau_t[gt]*sdc[gt].c)
+ + ism*sig[gt].Yv*rnd[d*2];
+ sd_V[n][d] = ism*sig[gt].V*rnd[d*2+1];
+
+ v[n][d] = vn*sdc[gt].em
+ + (invmass[n]*f[n][d] + accel[ga][d])*tau_t[gt]*(1 - sdc[gt].em)
+ + sd_V[n][d] - sdc[gt].em*Vmh;
+
+ xprime[n][d] = x[n][d] + v[n][d]*tau_t[gt]*(sdc[gt].eph - sdc[gt].emh);
+ }
+ else
+ {
+ /* Correct the velocities for the constraints.
+ * This operation introduces some inaccuracy,
+ * since the velocity is determined from differences in coordinates.
+ */
+ v[n][d] =
+ (xprime[n][d] - x[n][d])/(tau_t[gt]*(sdc[gt].eph - sdc[gt].emh));
+
+ Xmh = sd_V[n][d]*tau_t[gt]*sdc[gt].d/(sdc[gt].em-1)
+ + ism*sig[gt].Yx*rnd[d*2];
+ sd_X[n][d] = ism*sig[gt].X*rnd[d*2+1];
+
+ xprime[n][d] += sd_X[n][d] - Xmh;
+
+ }
+ }
+ else
+ {
+ if (bFirstHalf)
+ {
+ v[n][d] = 0.0;
+ xprime[n][d] = x[n][d];
+ }
+ }
+ }
+ }
+}
+
+static void do_update_bd_Tconsts(double dt, real friction_coefficient,
+ int ngtc, const real ref_t[],
+ real *rf)
+{
+ /* This is separated from the update below, because it is single threaded */
+ int gt;
+
+ if (friction_coefficient != 0)
+ {
+ for (gt = 0; gt < ngtc; gt++)
+ {
+ rf[gt] = sqrt(2.0*BOLTZ*ref_t[gt]/(friction_coefficient*dt));
+ }
+ }
+ else
+ {
+ for (gt = 0; gt < ngtc; gt++)
+ {
+ rf[gt] = sqrt(2.0*BOLTZ*ref_t[gt]);
+ }
+ }
+}
+
+static void do_update_bd(int start, int nrend, double dt,
+ ivec nFreeze[],
+ real invmass[], unsigned short ptype[],
+ unsigned short cFREEZE[], unsigned short cTC[],
+ rvec x[], rvec xprime[], rvec v[],
+ rvec f[], real friction_coefficient,
+ real *rf, gmx_int64_t step, int seed,
+ int* gatindex)
+{
+ /* note -- these appear to be full step velocities . . . */
+ int gf = 0, gt = 0;
+ real vn;
+ real invfr = 0;
+ int n, d;
+
+ if (friction_coefficient != 0)
+ {
+ invfr = 1.0/friction_coefficient;
+ }
+
+ for (n = start; (n < nrend); n++)
+ {
+ real rnd[3];
+ int ng = gatindex ? gatindex[n] : n;
+
+ if (cFREEZE)
+ {
+ gf = cFREEZE[n];
+ }
+ if (cTC)
+ {
+ gt = cTC[n];
+ }
+ gmx_rng_cycle_3gaussian_table(step, ng, seed, RND_SEED_UPDATE, rnd);
+ for (d = 0; (d < DIM); d++)
+ {
+ if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
+ {
+ if (friction_coefficient != 0)
+ {
+ vn = invfr*f[n][d] + rf[gt]*rnd[d];
+ }
+ else
+ {
+ /* NOTE: invmass = 2/(mass*friction_constant*dt) */
+ vn = 0.5*invmass[n]*f[n][d]*dt
+ + sqrt(0.5*invmass[n])*rf[gt]*rnd[d];
+ }
+
+ v[n][d] = vn;
+ xprime[n][d] = x[n][d]+vn*dt;
+ }
+ else
+ {
+ v[n][d] = 0.0;
+ xprime[n][d] = x[n][d];
+ }
+ }
+ }
+}
+
+static void dump_it_all(FILE gmx_unused *fp, const char gmx_unused *title,
+ int gmx_unused natoms, rvec gmx_unused x[], rvec gmx_unused xp[],
+ rvec gmx_unused v[], rvec gmx_unused f[])
+{
+#ifdef DEBUG
+ if (fp)
+ {
+ fprintf(fp, "%s\n", title);
+ pr_rvecs(fp, 0, "x", x, natoms);
+ pr_rvecs(fp, 0, "xp", xp, natoms);
+ pr_rvecs(fp, 0, "v", v, natoms);
+ pr_rvecs(fp, 0, "f", f, natoms);
+ }
+#endif
+}
+
+static void calc_ke_part_normal(rvec v[], t_grpopts *opts, t_mdatoms *md,
+ gmx_ekindata_t *ekind, t_nrnb *nrnb, gmx_bool bEkinAveVel,
+ gmx_bool bSaveEkinOld)
+{
+ int g;
+ t_grp_tcstat *tcstat = ekind->tcstat;
+ t_grp_acc *grpstat = ekind->grpstat;
+ int nthread, thread;
+
+ /* three main: VV with AveVel, vv with AveEkin, leap with AveEkin. Leap with AveVel is also
+ an option, but not supported now. Additionally, if we are doing iterations.
+ bEkinAveVel: If TRUE, we sum into ekin, if FALSE, into ekinh.
+ bSavEkinOld: If TRUE (in the case of iteration = bIterate is TRUE), we don't copy over the ekinh_old.
+ If FALSE, we overrwrite it.
+ */
+
+ /* group velocities are calculated in update_ekindata and
+ * accumulated in acumulate_groups.
+ * Now the partial global and groups ekin.
+ */
+ for (g = 0; (g < opts->ngtc); g++)
+ {
+
+ if (!bSaveEkinOld)
+ {
+ copy_mat(tcstat[g].ekinh, tcstat[g].ekinh_old);
+ }
+ if (bEkinAveVel)
+ {
+ clear_mat(tcstat[g].ekinf);
+ }
+ else
+ {
+ clear_mat(tcstat[g].ekinh);
+ }
+ if (bEkinAveVel)
+ {
+ tcstat[g].ekinscalef_nhc = 1.0; /* need to clear this -- logic is complicated! */
+ }
+ }
+ ekind->dekindl_old = ekind->dekindl;
+
+ nthread = gmx_omp_nthreads_get(emntUpdate);
+
+#pragma omp parallel for num_threads(nthread) schedule(static)
+ for (thread = 0; thread < nthread; thread++)
+ {
+ int start_t, end_t, n;
+ int ga, gt;
+ rvec v_corrt;
+ real hm;
+ int d, m;
+ matrix *ekin_sum;
+ real *dekindl_sum;
+
+ start_t = ((thread+0)*md->homenr)/nthread;
+ end_t = ((thread+1)*md->homenr)/nthread;
+
+ ekin_sum = ekind->ekin_work[thread];
+ dekindl_sum = ekind->dekindl_work[thread];
+
+ for (gt = 0; gt < opts->ngtc; gt++)
+ {
+ clear_mat(ekin_sum[gt]);
+ }
+ *dekindl_sum = 0.0;
+
+ ga = 0;
+ gt = 0;
+ for (n = start_t; n < end_t; n++)
+ {
+ if (md->cACC)
+ {
+ ga = md->cACC[n];
+ }
+ if (md->cTC)
+ {
+ gt = md->cTC[n];
+ }
+ hm = 0.5*md->massT[n];
+
+ for (d = 0; (d < DIM); d++)
+ {
+ v_corrt[d] = v[n][d] - grpstat[ga].u[d];
+ }
+ for (d = 0; (d < DIM); d++)
+ {
+ for (m = 0; (m < DIM); m++)
+ {
+ /* if we're computing a full step velocity, v_corrt[d] has v(t). Otherwise, v(t+dt/2) */
+ ekin_sum[gt][m][d] += hm*v_corrt[m]*v_corrt[d];
+ }
+ }
+ if (md->nMassPerturbed && md->bPerturbed[n])
+ {
+ *dekindl_sum +=
+ 0.5*(md->massB[n] - md->massA[n])*iprod(v_corrt, v_corrt);
+ }
+ }
+ }
+
+ ekind->dekindl = 0;
+ for (thread = 0; thread < nthread; thread++)
+ {
+ for (g = 0; g < opts->ngtc; g++)
+ {
+ if (bEkinAveVel)
+ {
+ m_add(tcstat[g].ekinf, ekind->ekin_work[thread][g],
+ tcstat[g].ekinf);
+ }
+ else
+ {
+ m_add(tcstat[g].ekinh, ekind->ekin_work[thread][g],
+ tcstat[g].ekinh);
+ }
+ }
+
+ ekind->dekindl += *ekind->dekindl_work[thread];
+ }
+
+ inc_nrnb(nrnb, eNR_EKIN, md->homenr);
+}
+
+static void calc_ke_part_visc(matrix box, rvec x[], rvec v[],
+ t_grpopts *opts, t_mdatoms *md,
+ gmx_ekindata_t *ekind,
+ t_nrnb *nrnb, gmx_bool bEkinAveVel)
+{
+ int start = 0, homenr = md->homenr;
+ int g, d, n, m, gt = 0;
+ rvec v_corrt;
+ real hm;
+ t_grp_tcstat *tcstat = ekind->tcstat;
+ t_cos_acc *cosacc = &(ekind->cosacc);
+ real dekindl;
+ real fac, cosz;
+ double mvcos;
+
+ for (g = 0; g < opts->ngtc; g++)
+ {
+ copy_mat(ekind->tcstat[g].ekinh, ekind->tcstat[g].ekinh_old);
+ clear_mat(ekind->tcstat[g].ekinh);
+ }
+ ekind->dekindl_old = ekind->dekindl;
+
+ fac = 2*M_PI/box[ZZ][ZZ];
+ mvcos = 0;
+ dekindl = 0;
+ for (n = start; n < start+homenr; n++)
+ {
+ if (md->cTC)
+ {
+ gt = md->cTC[n];
+ }
+ hm = 0.5*md->massT[n];
+
+ /* Note that the times of x and v differ by half a step */
+ /* MRS -- would have to be changed for VV */
+ cosz = cos(fac*x[n][ZZ]);
+ /* Calculate the amplitude of the new velocity profile */
+ mvcos += 2*cosz*md->massT[n]*v[n][XX];
+
+ copy_rvec(v[n], v_corrt);
+ /* Subtract the profile for the kinetic energy */
+ v_corrt[XX] -= cosz*cosacc->vcos;
+ for (d = 0; (d < DIM); d++)
+ {
+ for (m = 0; (m < DIM); m++)
+ {
+ /* if we're computing a full step velocity, v_corrt[d] has v(t). Otherwise, v(t+dt/2) */
+ if (bEkinAveVel)
+ {
+ tcstat[gt].ekinf[m][d] += hm*v_corrt[m]*v_corrt[d];
+ }
+ else
+ {
+ tcstat[gt].ekinh[m][d] += hm*v_corrt[m]*v_corrt[d];
+ }
+ }
+ }
+ if (md->nPerturbed && md->bPerturbed[n])
+ {
+ dekindl += 0.5*(md->massB[n] - md->massA[n])*iprod(v_corrt, v_corrt);
+ }
+ }
+ ekind->dekindl = dekindl;
+ cosacc->mvcos = mvcos;
+
+ inc_nrnb(nrnb, eNR_EKIN, homenr);
+}
+
+void calc_ke_part(t_state *state, t_grpopts *opts, t_mdatoms *md,
+ gmx_ekindata_t *ekind, t_nrnb *nrnb, gmx_bool bEkinAveVel, gmx_bool bSaveEkinOld)
+{
+ if (ekind->cosacc.cos_accel == 0)
+ {
+ calc_ke_part_normal(state->v, opts, md, ekind, nrnb, bEkinAveVel, bSaveEkinOld);
+ }
+ else
+ {
+ calc_ke_part_visc(state->box, state->x, state->v, opts, md, ekind, nrnb, bEkinAveVel);
+ }
+}
+
+extern void init_ekinstate(ekinstate_t *ekinstate, const t_inputrec *ir)
+{
+ ekinstate->ekin_n = ir->opts.ngtc;
+ snew(ekinstate->ekinh, ekinstate->ekin_n);
+ snew(ekinstate->ekinf, ekinstate->ekin_n);
+ snew(ekinstate->ekinh_old, ekinstate->ekin_n);
+ snew(ekinstate->ekinscalef_nhc, ekinstate->ekin_n);
+ snew(ekinstate->ekinscaleh_nhc, ekinstate->ekin_n);
+ snew(ekinstate->vscale_nhc, ekinstate->ekin_n);
+ ekinstate->dekindl = 0;
+ ekinstate->mvcos = 0;
+}
+
+void update_ekinstate(ekinstate_t *ekinstate, gmx_ekindata_t *ekind)
+{
+ int i;
+
+ for (i = 0; i < ekinstate->ekin_n; i++)
+ {
+ copy_mat(ekind->tcstat[i].ekinh, ekinstate->ekinh[i]);
+ copy_mat(ekind->tcstat[i].ekinf, ekinstate->ekinf[i]);
+ copy_mat(ekind->tcstat[i].ekinh_old, ekinstate->ekinh_old[i]);
+ ekinstate->ekinscalef_nhc[i] = ekind->tcstat[i].ekinscalef_nhc;
+ ekinstate->ekinscaleh_nhc[i] = ekind->tcstat[i].ekinscaleh_nhc;
+ ekinstate->vscale_nhc[i] = ekind->tcstat[i].vscale_nhc;
+ }
+
+ copy_mat(ekind->ekin, ekinstate->ekin_total);
+ ekinstate->dekindl = ekind->dekindl;
+ ekinstate->mvcos = ekind->cosacc.mvcos;
+
+}
+
+void restore_ekinstate_from_state(t_commrec *cr,
+ gmx_ekindata_t *ekind, ekinstate_t *ekinstate)
+{
+ int i, n;
+
+ if (MASTER(cr))
+ {
+ for (i = 0; i < ekinstate->ekin_n; i++)
+ {
+ copy_mat(ekinstate->ekinh[i], ekind->tcstat[i].ekinh);
+ copy_mat(ekinstate->ekinf[i], ekind->tcstat[i].ekinf);
+ copy_mat(ekinstate->ekinh_old[i], ekind->tcstat[i].ekinh_old);
+ ekind->tcstat[i].ekinscalef_nhc = ekinstate->ekinscalef_nhc[i];
+ ekind->tcstat[i].ekinscaleh_nhc = ekinstate->ekinscaleh_nhc[i];
+ ekind->tcstat[i].vscale_nhc = ekinstate->vscale_nhc[i];
+ }
+
+ copy_mat(ekinstate->ekin_total, ekind->ekin);
+
+ ekind->dekindl = ekinstate->dekindl;
+ ekind->cosacc.mvcos = ekinstate->mvcos;
+ n = ekinstate->ekin_n;
+ }
+
+ if (PAR(cr))
+ {
+ gmx_bcast(sizeof(n), &n, cr);
+ for (i = 0; i < n; i++)
+ {
+ gmx_bcast(DIM*DIM*sizeof(ekind->tcstat[i].ekinh[0][0]),
+ ekind->tcstat[i].ekinh[0], cr);
+ gmx_bcast(DIM*DIM*sizeof(ekind->tcstat[i].ekinf[0][0]),
+ ekind->tcstat[i].ekinf[0], cr);
+ gmx_bcast(DIM*DIM*sizeof(ekind->tcstat[i].ekinh_old[0][0]),
+ ekind->tcstat[i].ekinh_old[0], cr);
+
+ gmx_bcast(sizeof(ekind->tcstat[i].ekinscalef_nhc),
+ &(ekind->tcstat[i].ekinscalef_nhc), cr);
+ gmx_bcast(sizeof(ekind->tcstat[i].ekinscaleh_nhc),
+ &(ekind->tcstat[i].ekinscaleh_nhc), cr);
+ gmx_bcast(sizeof(ekind->tcstat[i].vscale_nhc),
+ &(ekind->tcstat[i].vscale_nhc), cr);
+ }
+ gmx_bcast(DIM*DIM*sizeof(ekind->ekin[0][0]),
+ ekind->ekin[0], cr);
+
+ gmx_bcast(sizeof(ekind->dekindl), &ekind->dekindl, cr);
+ gmx_bcast(sizeof(ekind->cosacc.mvcos), &ekind->cosacc.mvcos, cr);
+ }
+}
+
+void set_deform_reference_box(gmx_update_t upd, gmx_int64_t step, matrix box)
+{
+ upd->deformref_step = step;
+ copy_mat(box, upd->deformref_box);
+}
+
+static void deform(gmx_update_t upd,
+ int start, int homenr, rvec x[], matrix box, matrix *scale_tot,
+ const t_inputrec *ir, gmx_int64_t step)
+{
+ matrix bnew, invbox, mu;
+ real elapsed_time;
+ int i, j;
+
+ elapsed_time = (step + 1 - upd->deformref_step)*ir->delta_t;
+ copy_mat(box, bnew);
+ for (i = 0; i < DIM; i++)
+ {
+ for (j = 0; j < DIM; j++)
+ {
+ if (ir->deform[i][j] != 0)
+ {
+ bnew[i][j] =
+ upd->deformref_box[i][j] + elapsed_time*ir->deform[i][j];
+ }
+ }
+ }
+ /* We correct the off-diagonal elements,
+ * which can grow indefinitely during shearing,
+ * so the shifts do not get messed up.
+ */
+ for (i = 1; i < DIM; i++)
+ {
+ for (j = i-1; j >= 0; j--)
+ {
+ while (bnew[i][j] - box[i][j] > 0.5*bnew[j][j])
+ {
+ rvec_dec(bnew[i], bnew[j]);
+ }
+ while (bnew[i][j] - box[i][j] < -0.5*bnew[j][j])
+ {
+ rvec_inc(bnew[i], bnew[j]);
+ }
+ }
+ }
+ m_inv_ur0(box, invbox);
+ copy_mat(bnew, box);
+ mmul_ur0(box, invbox, mu);
+
+ for (i = start; i < start+homenr; i++)
+ {
+ x[i][XX] = mu[XX][XX]*x[i][XX]+mu[YY][XX]*x[i][YY]+mu[ZZ][XX]*x[i][ZZ];
+ x[i][YY] = mu[YY][YY]*x[i][YY]+mu[ZZ][YY]*x[i][ZZ];
+ x[i][ZZ] = mu[ZZ][ZZ]*x[i][ZZ];
+ }
+ if (*scale_tot)
+ {
+ /* The transposes of the scaling matrices are stored,
+ * so we need to do matrix multiplication in the inverse order.
+ */
+ mmul_ur0(*scale_tot, mu, *scale_tot);
+ }
+}
+
- combine_forces(inputrec->nstcalclr, constr, inputrec, md, idef, cr,
+void update_tcouple(gmx_int64_t step,
+ t_inputrec *inputrec,
+ t_state *state,
+ gmx_ekindata_t *ekind,
+ t_extmass *MassQ,
+ t_mdatoms *md)
+
+{
+ gmx_bool bTCouple = FALSE;
+ real dttc;
+ int i, start, end, homenr, offset;
+
+ /* if using vv with trotter decomposition methods, we do this elsewhere in the code */
+ if (inputrec->etc != etcNO &&
+ !(IR_NVT_TROTTER(inputrec) || IR_NPT_TROTTER(inputrec) || IR_NPH_TROTTER(inputrec)))
+ {
+ /* We should only couple after a step where energies were determined (for leapfrog versions)
+ or the step energies are determined, for velocity verlet versions */
+
+ if (EI_VV(inputrec->eI))
+ {
+ offset = 0;
+ }
+ else
+ {
+ offset = 1;
+ }
+ bTCouple = (inputrec->nsttcouple == 1 ||
+ do_per_step(step+inputrec->nsttcouple-offset,
+ inputrec->nsttcouple));
+ }
+
+ if (bTCouple)
+ {
+ dttc = inputrec->nsttcouple*inputrec->delta_t;
+
+ switch (inputrec->etc)
+ {
+ case etcNO:
+ break;
+ case etcBERENDSEN:
+ berendsen_tcoupl(inputrec, ekind, dttc);
+ break;
+ case etcNOSEHOOVER:
+ nosehoover_tcoupl(&(inputrec->opts), ekind, dttc,
+ state->nosehoover_xi, state->nosehoover_vxi, MassQ);
+ break;
+ case etcVRESCALE:
+ vrescale_tcoupl(inputrec, step, ekind, dttc,
+ state->therm_integral);
+ break;
+ }
+ /* rescale in place here */
+ if (EI_VV(inputrec->eI))
+ {
+ rescale_velocities(ekind, md, 0, md->homenr, state->v);
+ }
+ }
+ else
+ {
+ /* Set the T scaling lambda to 1 to have no scaling */
+ for (i = 0; (i < inputrec->opts.ngtc); i++)
+ {
+ ekind->tcstat[i].lambda = 1.0;
+ }
+ }
+}
+
+void update_pcouple(FILE *fplog,
+ gmx_int64_t step,
+ t_inputrec *inputrec,
+ t_state *state,
+ matrix pcoupl_mu,
+ matrix M,
+ gmx_bool bInitStep)
+{
+ gmx_bool bPCouple = FALSE;
+ real dtpc = 0;
+ int i;
+
+ /* if using Trotter pressure, we do this in coupling.c, so we leave it false. */
+ if (inputrec->epc != epcNO && (!(IR_NPT_TROTTER(inputrec) || IR_NPH_TROTTER(inputrec))))
+ {
+ /* We should only couple after a step where energies were determined */
+ bPCouple = (inputrec->nstpcouple == 1 ||
+ do_per_step(step+inputrec->nstpcouple-1,
+ inputrec->nstpcouple));
+ }
+
+ clear_mat(pcoupl_mu);
+ for (i = 0; i < DIM; i++)
+ {
+ pcoupl_mu[i][i] = 1.0;
+ }
+
+ clear_mat(M);
+
+ if (bPCouple)
+ {
+ dtpc = inputrec->nstpcouple*inputrec->delta_t;
+
+ switch (inputrec->epc)
+ {
+ /* We can always pcoupl, even if we did not sum the energies
+ * the previous step, since state->pres_prev is only updated
+ * when the energies have been summed.
+ */
+ case (epcNO):
+ break;
+ case (epcBERENDSEN):
+ if (!bInitStep)
+ {
+ berendsen_pcoupl(fplog, step, inputrec, dtpc, state->pres_prev, state->box,
+ pcoupl_mu);
+ }
+ break;
+ case (epcPARRINELLORAHMAN):
+ parrinellorahman_pcoupl(fplog, step, inputrec, dtpc, state->pres_prev,
+ state->box, state->box_rel, state->boxv,
+ M, pcoupl_mu, bInitStep);
+ break;
+ default:
+ break;
+ }
+ }
+}
+
+static rvec *get_xprime(const t_state *state, gmx_update_t upd)
+{
+ if (state->nalloc > upd->xp_nalloc)
+ {
+ upd->xp_nalloc = state->nalloc;
+ srenew(upd->xp, upd->xp_nalloc);
+ }
+
+ return upd->xp;
+}
+
++static void combine_forces(gmx_update_t upd,
++ int nstcalclr,
++ gmx_constr_t constr,
++ t_inputrec *ir, t_mdatoms *md, t_idef *idef,
++ t_commrec *cr,
++ gmx_int64_t step,
++ t_state *state, gmx_bool bMolPBC,
++ int start, int nrend,
++ rvec f[], rvec f_lr[],
++ tensor *vir_lr_constr,
++ t_nrnb *nrnb)
++{
++ int i, d;
++
++ /* f contains the short-range forces + the long range forces
++ * which are stored separately in f_lr.
++ */
++
++ if (constr != NULL && vir_lr_constr != NULL &&
++ !(ir->eConstrAlg == econtSHAKE && ir->epc == epcNO))
++ {
++ /* We need to constrain the LR forces separately,
++ * because due to the different pre-factor for the SR and LR
++ * forces in the update algorithm, we have to correct
++ * the constraint virial for the nstcalclr-1 extra f_lr.
++ * Constrain only the additional LR part of the force.
++ */
++ /* MRS -- need to make sure this works with trotter integration -- the constraint calls may not be right.*/
++ rvec *xp;
++ real fac;
++ int gf = 0;
++
++ xp = get_xprime(state, upd);
++
++ fac = (nstcalclr - 1)*ir->delta_t*ir->delta_t;
++
++ for (i = 0; i < md->homenr; i++)
++ {
++ if (md->cFREEZE != NULL)
++ {
++ gf = md->cFREEZE[i];
++ }
++ for (d = 0; d < DIM; d++)
++ {
++ if ((md->ptype[i] != eptVSite) &&
++ (md->ptype[i] != eptShell) &&
++ !ir->opts.nFreeze[gf][d])
++ {
++ xp[i][d] = state->x[i][d] + fac*f_lr[i][d]*md->invmass[i];
++ }
++ }
++ }
++ constrain(NULL, FALSE, FALSE, constr, idef, ir, NULL, cr, step, 0, md,
++ state->x, xp, xp, bMolPBC, state->box, state->lambda[efptBONDED], NULL,
++ NULL, vir_lr_constr, nrnb, econqCoord, ir->epc == epcMTTK, state->veta, state->veta);
++ }
++
++ /* Add nstcalclr-1 times the LR force to the sum of both forces
++ * and store the result in forces_lr.
++ */
++ for (i = start; i < nrend; i++)
++ {
++ for (d = 0; d < DIM; d++)
++ {
++ f_lr[i][d] = f[i][d] + (nstcalclr - 1)*f_lr[i][d];
++ }
++ }
++}
++
+void update_constraints(FILE *fplog,
+ gmx_int64_t step,
+ real *dvdlambda, /* the contribution to be added to the bonded interactions */
+ t_inputrec *inputrec, /* input record and box stuff */
+ gmx_ekindata_t *ekind,
+ t_mdatoms *md,
+ t_state *state,
+ gmx_bool bMolPBC,
+ t_graph *graph,
+ rvec force[], /* forces on home particles */
+ t_idef *idef,
+ tensor vir_part,
+ t_commrec *cr,
+ t_nrnb *nrnb,
+ gmx_wallcycle_t wcycle,
+ gmx_update_t upd,
+ gmx_constr_t constr,
+ gmx_bool bFirstHalf,
+ gmx_bool bCalcVir,
+ real vetanew)
+{
+ gmx_bool bExtended, bLastStep, bLog = FALSE, bEner = FALSE, bDoConstr = FALSE;
+ double dt;
+ real dt_1;
+ int start, homenr, nrend, i, n, m, g, d;
+ tensor vir_con;
+ rvec *vbuf, *xprime = NULL;
+ int nth, th;
+
+ if (constr)
+ {
+ bDoConstr = TRUE;
+ }
+ if (bFirstHalf && !EI_VV(inputrec->eI))
+ {
+ bDoConstr = FALSE;
+ }
+
+ /* for now, SD update is here -- though it really seems like it
+ should be reformulated as a velocity verlet method, since it has two parts */
+
+ start = 0;
+ homenr = md->homenr;
+ nrend = start+homenr;
+
+ dt = inputrec->delta_t;
+ dt_1 = 1.0/dt;
+
+ /*
+ * Steps (7C, 8C)
+ * APPLY CONSTRAINTS:
+ * BLOCK SHAKE
+
+ * When doing PR pressure coupling we have to constrain the
+ * bonds in each iteration. If we are only using Nose-Hoover tcoupling
+ * it is enough to do this once though, since the relative velocities
+ * after this will be normal to the bond vector
+ */
+
+ if (bDoConstr)
+ {
+ /* clear out constraints before applying */
+ clear_mat(vir_part);
+
+ xprime = get_xprime(state, upd);
+
+ bLastStep = (step == inputrec->init_step+inputrec->nsteps);
+ bLog = (do_per_step(step, inputrec->nstlog) || bLastStep || (step < 0));
+ bEner = (do_per_step(step, inputrec->nstenergy) || bLastStep);
+ /* Constrain the coordinates xprime */
+ wallcycle_start(wcycle, ewcCONSTR);
+ if (EI_VV(inputrec->eI) && bFirstHalf)
+ {
+ constrain(NULL, bLog, bEner, constr, idef,
+ inputrec, ekind, cr, step, 1, md,
+ state->x, state->v, state->v,
+ bMolPBC, state->box,
+ state->lambda[efptBONDED], dvdlambda,
+ NULL, bCalcVir ? &vir_con : NULL, nrnb, econqVeloc,
+ inputrec->epc == epcMTTK, state->veta, vetanew);
+ }
+ else
+ {
+ constrain(NULL, bLog, bEner, constr, idef,
+ inputrec, ekind, cr, step, 1, md,
+ state->x, xprime, NULL,
+ bMolPBC, state->box,
+ state->lambda[efptBONDED], dvdlambda,
+ state->v, bCalcVir ? &vir_con : NULL, nrnb, econqCoord,
+ inputrec->epc == epcMTTK, state->veta, state->veta);
+ }
+ wallcycle_stop(wcycle, ewcCONSTR);
+
+ where();
+
+ dump_it_all(fplog, "After Shake",
+ state->natoms, state->x, xprime, state->v, force);
+
+ if (bCalcVir)
+ {
+ if (inputrec->eI == eiSD2)
+ {
+ /* A correction factor eph is needed for the SD constraint force */
+ /* Here we can, unfortunately, not have proper corrections
+ * for different friction constants, so we use the first one.
+ */
+ for (i = 0; i < DIM; i++)
+ {
+ for (m = 0; m < DIM; m++)
+ {
+ vir_part[i][m] += upd->sd->sdc[0].eph*vir_con[i][m];
+ }
+ }
+ }
+ else
+ {
+ m_add(vir_part, vir_con, vir_part);
+ }
+ if (debug)
+ {
+ pr_rvecs(debug, 0, "constraint virial", vir_part, DIM);
+ }
+ }
+ }
+
+ where();
+ if ((inputrec->eI == eiSD2) && !(bFirstHalf))
+ {
+ xprime = get_xprime(state, upd);
+
+ nth = gmx_omp_nthreads_get(emntUpdate);
+
+#pragma omp parallel for num_threads(nth) schedule(static)
+ for (th = 0; th < nth; th++)
+ {
+ int start_th, end_th;
+
+ start_th = start + ((nrend-start)* th )/nth;
+ end_th = start + ((nrend-start)*(th+1))/nth;
+
+ /* The second part of the SD integration */
+ do_update_sd2(upd->sd,
+ FALSE, start_th, end_th,
+ inputrec->opts.acc, inputrec->opts.nFreeze,
+ md->invmass, md->ptype,
+ md->cFREEZE, md->cACC, md->cTC,
+ state->x, xprime, state->v, force, state->sd_X,
+ inputrec->opts.tau_t,
+ FALSE, step, inputrec->ld_seed,
+ DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL);
+ }
+ inc_nrnb(nrnb, eNR_UPDATE, homenr);
+
+ if (bDoConstr)
+ {
+ /* Constrain the coordinates xprime */
+ wallcycle_start(wcycle, ewcCONSTR);
+ constrain(NULL, bLog, bEner, constr, idef,
+ inputrec, NULL, cr, step, 1, md,
+ state->x, xprime, NULL,
+ bMolPBC, state->box,
+ state->lambda[efptBONDED], dvdlambda,
+ NULL, NULL, nrnb, econqCoord, FALSE, 0, 0);
+ wallcycle_stop(wcycle, ewcCONSTR);
+ }
+ }
+
+ /* We must always unshift after updating coordinates; if we did not shake
+ x was shifted in do_force */
+
+ if (!(bFirstHalf)) /* in the first half of vv, no shift. */
+ {
+ if (graph && (graph->nnodes > 0))
+ {
+ unshift_x(graph, state->box, state->x, upd->xp);
+ if (TRICLINIC(state->box))
+ {
+ inc_nrnb(nrnb, eNR_SHIFTX, 2*graph->nnodes);
+ }
+ else
+ {
+ inc_nrnb(nrnb, eNR_SHIFTX, graph->nnodes);
+ }
+ }
+ else
+ {
+#pragma omp parallel for num_threads(gmx_omp_nthreads_get(emntUpdate)) schedule(static)
+ for (i = start; i < nrend; i++)
+ {
+ copy_rvec(upd->xp[i], state->x[i]);
+ }
+ }
+
+ dump_it_all(fplog, "After unshift",
+ state->natoms, state->x, upd->xp, state->v, force);
+ }
+/* ############# END the update of velocities and positions ######### */
+}
+
+void update_box(FILE *fplog,
+ gmx_int64_t step,
+ t_inputrec *inputrec, /* input record and box stuff */
+ t_mdatoms *md,
+ t_state *state,
+ rvec force[], /* forces on home particles */
+ matrix *scale_tot,
+ matrix pcoupl_mu,
+ t_nrnb *nrnb,
+ gmx_update_t upd)
+{
+ gmx_bool bExtended, bLastStep, bLog = FALSE, bEner = FALSE;
+ double dt;
+ real dt_1;
+ int start, homenr, nrend, i, n, m, g;
+ tensor vir_con;
+
+ start = 0;
+ homenr = md->homenr;
+ nrend = start+homenr;
+
+ bExtended =
+ (inputrec->etc == etcNOSEHOOVER) ||
+ (inputrec->epc == epcPARRINELLORAHMAN) ||
+ (inputrec->epc == epcMTTK);
+
+ dt = inputrec->delta_t;
+
+ where();
+
+ /* now update boxes */
+ switch (inputrec->epc)
+ {
+ case (epcNO):
+ break;
+ case (epcBERENDSEN):
+ berendsen_pscale(inputrec, pcoupl_mu, state->box, state->box_rel,
+ start, homenr, state->x, md->cFREEZE, nrnb);
+ break;
+ case (epcPARRINELLORAHMAN):
+ /* The box velocities were updated in do_pr_pcoupl in the update
+ * iteration, but we dont change the box vectors until we get here
+ * since we need to be able to shift/unshift above.
+ */
+ for (i = 0; i < DIM; i++)
+ {
+ for (m = 0; m <= i; m++)
+ {
+ state->box[i][m] += dt*state->boxv[i][m];
+ }
+ }
+ preserve_box_shape(inputrec, state->box_rel, state->box);
+
+ /* Scale the coordinates */
+ for (n = start; (n < start+homenr); n++)
+ {
+ tmvmul_ur0(pcoupl_mu, state->x[n], state->x[n]);
+ }
+ break;
+ case (epcMTTK):
+ switch (inputrec->epct)
+ {
+ case (epctISOTROPIC):
+ /* DIM * eta = ln V. so DIM*eta_new = DIM*eta_old + DIM*dt*veta =>
+ ln V_new = ln V_old + 3*dt*veta => V_new = V_old*exp(3*dt*veta) =>
+ Side length scales as exp(veta*dt) */
+
+ msmul(state->box, exp(state->veta*dt), state->box);
+
+ /* Relate veta to boxv. veta = d(eta)/dT = (1/DIM)*1/V dV/dT.
+ o If we assume isotropic scaling, and box length scaling
+ factor L, then V = L^DIM (det(M)). So dV/dt = DIM
+ L^(DIM-1) dL/dt det(M), and veta = (1/L) dL/dt. The
+ determinant of B is L^DIM det(M), and the determinant
+ of dB/dt is (dL/dT)^DIM det (M). veta will be
+ (det(dB/dT)/det(B))^(1/3). Then since M =
+ B_new*(vol_new)^(1/3), dB/dT_new = (veta_new)*B(new). */
+
+ msmul(state->box, state->veta, state->boxv);
+ break;
+ default:
+ break;
+ }
+ break;
+ default:
+ break;
+ }
+
+ if ((!(IR_NPT_TROTTER(inputrec) || IR_NPH_TROTTER(inputrec))) && scale_tot)
+ {
+ /* The transposes of the scaling matrices are stored,
+ * therefore we need to reverse the order in the multiplication.
+ */
+ mmul_ur0(*scale_tot, pcoupl_mu, *scale_tot);
+ }
+
+ if (DEFORM(*inputrec))
+ {
+ deform(upd, start, homenr, state->x, state->box, scale_tot, inputrec, step);
+ }
+ where();
+ dump_it_all(fplog, "After update",
+ state->natoms, state->x, upd->xp, state->v, force);
+}
+
+void update_coords(FILE *fplog,
+ gmx_int64_t step,
+ t_inputrec *inputrec, /* input record and box stuff */
+ t_mdatoms *md,
+ t_state *state,
+ gmx_bool bMolPBC,
+ rvec *f, /* forces on home particles */
+ gmx_bool bDoLR,
+ rvec *f_lr,
++ tensor *vir_lr_constr,
+ t_fcdata *fcd,
+ gmx_ekindata_t *ekind,
+ matrix M,
+ gmx_update_t upd,
+ gmx_bool bInitStep,
+ int UpdatePart,
+ t_commrec *cr, /* these shouldn't be here -- need to think about it */
+ t_nrnb *nrnb,
+ gmx_constr_t constr,
+ t_idef *idef)
+{
+ gmx_bool bNH, bPR, bLastStep, bLog = FALSE, bEner = FALSE;
+ double dt, alpha;
+ real *imass, *imassin;
+ rvec *force;
+ real dt_1;
+ int start, homenr, nrend, i, j, d, n, m, g;
+ int blen0, blen1, iatom, jatom, nshake, nsettle, nconstr, nexpand;
+ int *icom = NULL;
+ tensor vir_con;
+ rvec *vcom, *xcom, *vall, *xall, *xin, *vin, *forcein, *fall, *xpall, *xprimein, *xprime;
+ int nth, th;
+
+ /* Running the velocity half does nothing except for velocity verlet */
+ if ((UpdatePart == etrtVELOCITY1 || UpdatePart == etrtVELOCITY2) &&
+ !EI_VV(inputrec->eI))
+ {
+ gmx_incons("update_coords called for velocity without VV integrator");
+ }
+
+ start = 0;
+ homenr = md->homenr;
+ nrend = start+homenr;
+
+ xprime = get_xprime(state, upd);
+
+ dt = inputrec->delta_t;
+ dt_1 = 1.0/dt;
+
+ /* We need to update the NMR restraint history when time averaging is used */
+ if (state->flags & (1<<estDISRE_RM3TAV))
+ {
+ update_disres_history(fcd, &state->hist);
+ }
+ if (state->flags & (1<<estORIRE_DTAV))
+ {
+ update_orires_history(fcd, &state->hist);
+ }
+
+
+ bNH = inputrec->etc == etcNOSEHOOVER;
+ bPR = ((inputrec->epc == epcPARRINELLORAHMAN) || (inputrec->epc == epcMTTK));
+
+ if (bDoLR && inputrec->nstcalclr > 1 && !EI_VV(inputrec->eI)) /* get this working with VV? */
+ {
+ /* Store the total force + nstcalclr-1 times the LR force
+ * in forces_lr, so it can be used in a normal update algorithm
+ * to produce twin time stepping.
+ */
+ /* is this correct in the new construction? MRS */
- start, nrend, f, f_lr, nrnb);
++ combine_forces(upd,
++ inputrec->nstcalclr, constr, inputrec, md, idef, cr,
+ step, state, bMolPBC,
++ start, nrend, f, f_lr, vir_lr_constr, nrnb);
+ force = f_lr;
+ }
+ else
+ {
+ force = f;
+ }
+
+ /* ############# START The update of velocities and positions ######### */
+ where();
+ dump_it_all(fplog, "Before update",
+ state->natoms, state->x, xprime, state->v, force);
+
+ if (inputrec->eI == eiSD2)
+ {
+ check_sd2_work_data_allocation(upd->sd, nrend);
+
+ do_update_sd2_Tconsts(upd->sd,
+ inputrec->opts.ngtc,
+ inputrec->opts.tau_t,
+ inputrec->opts.ref_t);
+ }
+ if (inputrec->eI == eiBD)
+ {
+ do_update_bd_Tconsts(dt, inputrec->bd_fric,
+ inputrec->opts.ngtc, inputrec->opts.ref_t,
+ upd->sd->bd_rf);
+ }
+
+ nth = gmx_omp_nthreads_get(emntUpdate);
+
+#pragma omp parallel for num_threads(nth) schedule(static) private(alpha)
+ for (th = 0; th < nth; th++)
+ {
+ int start_th, end_th;
+
+ start_th = start + ((nrend-start)* th )/nth;
+ end_th = start + ((nrend-start)*(th+1))/nth;
+
+ switch (inputrec->eI)
+ {
+ case (eiMD):
+ if (ekind->cosacc.cos_accel == 0)
+ {
+ do_update_md(start_th, end_th, dt,
+ ekind->tcstat, state->nosehoover_vxi,
+ ekind->bNEMD, ekind->grpstat, inputrec->opts.acc,
+ inputrec->opts.nFreeze,
+ md->invmass, md->ptype,
+ md->cFREEZE, md->cACC, md->cTC,
+ state->x, xprime, state->v, force, M,
+ bNH, bPR);
+ }
+ else
+ {
+ do_update_visc(start_th, end_th, dt,
+ ekind->tcstat, state->nosehoover_vxi,
+ md->invmass, md->ptype,
+ md->cTC, state->x, xprime, state->v, force, M,
+ state->box,
+ ekind->cosacc.cos_accel,
+ ekind->cosacc.vcos,
+ bNH, bPR);
+ }
+ break;
+ case (eiSD1):
+ do_update_sd1(upd->sd,
+ start_th, end_th, dt,
+ inputrec->opts.acc, inputrec->opts.nFreeze,
+ md->invmass, md->ptype,
+ md->cFREEZE, md->cACC, md->cTC,
+ state->x, xprime, state->v, force,
+ inputrec->opts.ngtc, inputrec->opts.tau_t, inputrec->opts.ref_t,
+ step, inputrec->ld_seed, DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL);
+ break;
+ case (eiSD2):
+ /* The SD update is done in 2 parts, because an extra constraint step
+ * is needed
+ */
+ do_update_sd2(upd->sd,
+ bInitStep, start_th, end_th,
+ inputrec->opts.acc, inputrec->opts.nFreeze,
+ md->invmass, md->ptype,
+ md->cFREEZE, md->cACC, md->cTC,
+ state->x, xprime, state->v, force, state->sd_X,
+ inputrec->opts.tau_t,
+ TRUE, step, inputrec->ld_seed,
+ DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL);
+ break;
+ case (eiBD):
+ do_update_bd(start_th, end_th, dt,
+ inputrec->opts.nFreeze, md->invmass, md->ptype,
+ md->cFREEZE, md->cTC,
+ state->x, xprime, state->v, force,
+ inputrec->bd_fric,
+ upd->sd->bd_rf,
+ step, inputrec->ld_seed, DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL);
+ break;
+ case (eiVV):
+ case (eiVVAK):
+ alpha = 1.0 + DIM/((double)inputrec->opts.nrdf[0]); /* assuming barostat coupled to group 0. */
+ switch (UpdatePart)
+ {
+ case etrtVELOCITY1:
+ case etrtVELOCITY2:
+ do_update_vv_vel(start_th, end_th, dt,
+ inputrec->opts.acc, inputrec->opts.nFreeze,
+ md->invmass, md->ptype,
+ md->cFREEZE, md->cACC,
+ state->v, force,
+ (bNH || bPR), state->veta, alpha);
+ break;
+ case etrtPOSITION:
+ do_update_vv_pos(start_th, end_th, dt,
+ inputrec->opts.nFreeze,
+ md->ptype, md->cFREEZE,
+ state->x, xprime, state->v,
+ (bNH || bPR), state->veta);
+ break;
+ }
+ break;
+ default:
+ gmx_fatal(FARGS, "Don't know how to update coordinates");
+ break;
+ }
+ }
+
+}
+
+
+void correct_ekin(FILE *log, int start, int end, rvec v[], rvec vcm, real mass[],
+ real tmass, tensor ekin)
+{
+ /*
+ * This is a debugging routine. It should not be called for production code
+ *
+ * The kinetic energy should calculated according to:
+ * Ekin = 1/2 m (v-vcm)^2
+ * However the correction is not always applied, since vcm may not be
+ * known in time and we compute
+ * Ekin' = 1/2 m v^2 instead
+ * This can be corrected afterwards by computing
+ * Ekin = Ekin' + 1/2 m ( -2 v vcm + vcm^2)
+ * or in hsorthand:
+ * Ekin = Ekin' - m v vcm + 1/2 m vcm^2
+ */
+ int i, j, k;
+ real m, tm;
+ rvec hvcm, mv;
+ tensor dekin;
+
+ /* Local particles */
+ clear_rvec(mv);
+
+ /* Processor dependent part. */
+ tm = 0;
+ for (i = start; (i < end); i++)
+ {
+ m = mass[i];
+ tm += m;
+ for (j = 0; (j < DIM); j++)
+ {
+ mv[j] += m*v[i][j];
+ }
+ }
+ /* Shortcut */
+ svmul(1/tmass, vcm, vcm);
+ svmul(0.5, vcm, hvcm);
+ clear_mat(dekin);
+ for (j = 0; (j < DIM); j++)
+ {
+ for (k = 0; (k < DIM); k++)
+ {
+ dekin[j][k] += vcm[k]*(tm*hvcm[j]-mv[j]);
+ }
+ }
+ pr_rvecs(log, 0, "dekin", dekin, DIM);
+ pr_rvecs(log, 0, " ekin", ekin, DIM);
+ fprintf(log, "dekin = %g, ekin = %g vcm = (%8.4f %8.4f %8.4f)\n",
+ trace(dekin), trace(ekin), vcm[XX], vcm[YY], vcm[ZZ]);
+ fprintf(log, "mv = (%8.4f %8.4f %8.4f)\n",
+ mv[XX], mv[YY], mv[ZZ]);
+}
+
+extern gmx_bool update_randomize_velocities(t_inputrec *ir, gmx_int64_t step, const t_commrec *cr,
+ t_mdatoms *md, t_state *state, gmx_update_t upd, gmx_constr_t constr)
+{
+
+ int i;
+ real rate = (ir->delta_t)/ir->opts.tau_t[0];
+
+ if (ir->etc == etcANDERSEN && constr != NULL)
+ {
+ gmx_fatal(FARGS, "Normal Andersen is currently not supported with constraints, use massive Andersen instead");
+ }
+
+ /* proceed with andersen if 1) it's fixed probability per
+ particle andersen or 2) it's massive andersen and it's tau_t/dt */
+ if ((ir->etc == etcANDERSEN) || do_per_step(step, (int)(1.0/rate)))
+ {
+ andersen_tcoupl(ir, step, cr, md, state, rate,
+ upd->sd->randomize_group, upd->sd->boltzfac);
+ return TRUE;
+ }
+ return FALSE;
+}
--- /dev/null
- f, bUpdateDoLR, fr->f_twin, fcd,
+/*
+ * 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) 2011,2012,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.
+ *
+ * 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 <config.h>
+#endif
+
+#include "typedefs.h"
+#include "gromacs/utility/smalloc.h"
+#include "sysstuff.h"
+#include "vec.h"
+#include "vcm.h"
+#include "mdebin.h"
+#include "nrnb.h"
+#include "calcmu.h"
+#include "index.h"
+#include "vsite.h"
+#include "update.h"
+#include "ns.h"
+#include "mdrun.h"
+#include "md_support.h"
+#include "md_logging.h"
+#include "network.h"
+#include "xvgr.h"
+#include "physics.h"
+#include "names.h"
+#include "force.h"
+#include "disre.h"
+#include "orires.h"
+#include "pme.h"
+#include "mdatoms.h"
+#include "repl_ex.h"
+#include "deform.h"
+#include "qmmm.h"
+#include "domdec.h"
+#include "domdec_network.h"
+#include "gromacs/gmxlib/topsort.h"
+#include "coulomb.h"
+#include "constr.h"
+#include "shellfc.h"
+#include "gromacs/gmxpreprocess/compute_io.h"
+#include "checkpoint.h"
+#include "mtop_util.h"
+#include "sighandler.h"
+#include "txtdump.h"
+#include "gromacs/utility/cstringutil.h"
+#include "pme_loadbal.h"
+#include "bondf.h"
+#include "membed.h"
+#include "types/nlistheuristics.h"
+#include "types/iteratedconstraints.h"
+#include "nbnxn_cuda_data_mgmt.h"
+
+#include "gromacs/utility/gmxmpi.h"
+#include "gromacs/fileio/confio.h"
+#include "gromacs/fileio/trajectory_writing.h"
+#include "gromacs/fileio/trnio.h"
+#include "gromacs/fileio/trxio.h"
+#include "gromacs/fileio/xtcio.h"
+#include "gromacs/timing/wallcycle.h"
+#include "gromacs/timing/walltime_accounting.h"
+#include "gromacs/pulling/pull.h"
+#include "gromacs/swap/swapcoords.h"
+#include "gromacs/imd/imd.h"
+
+#ifdef GMX_FAHCORE
+#include "corewrap.h"
+#endif
+
+static void reset_all_counters(FILE *fplog, t_commrec *cr,
+ gmx_int64_t step,
+ gmx_int64_t *step_rel, t_inputrec *ir,
+ gmx_wallcycle_t wcycle, t_nrnb *nrnb,
+ gmx_walltime_accounting_t walltime_accounting,
+ nbnxn_cuda_ptr_t cu_nbv)
+{
+ char sbuf[STEPSTRSIZE];
+
+ /* Reset all the counters related to performance over the run */
+ md_print_warn(cr, fplog, "step %s: resetting all time and cycle counters\n",
+ gmx_step_str(step, sbuf));
+
+ if (cu_nbv)
+ {
+ nbnxn_cuda_reset_timings(cu_nbv);
+ }
+
+ wallcycle_stop(wcycle, ewcRUN);
+ wallcycle_reset_all(wcycle);
+ if (DOMAINDECOMP(cr))
+ {
+ reset_dd_statistics_counters(cr->dd);
+ }
+ init_nrnb(nrnb);
+ ir->init_step += *step_rel;
+ ir->nsteps -= *step_rel;
+ *step_rel = 0;
+ wallcycle_start(wcycle, ewcRUN);
+ walltime_accounting_start(walltime_accounting);
+ print_date_and_time(fplog, cr->nodeid, "Restarted time", gmx_gettime());
+}
+
+double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
+ const output_env_t oenv, gmx_bool bVerbose, gmx_bool bCompact,
+ int nstglobalcomm,
+ gmx_vsite_t *vsite, gmx_constr_t constr,
+ int stepout, t_inputrec *ir,
+ gmx_mtop_t *top_global,
+ t_fcdata *fcd,
+ t_state *state_global,
+ t_mdatoms *mdatoms,
+ t_nrnb *nrnb, gmx_wallcycle_t wcycle,
+ gmx_edsam_t ed, t_forcerec *fr,
+ int repl_ex_nst, int repl_ex_nex, int repl_ex_seed, gmx_membed_t membed,
+ real cpt_period, real max_hours,
+ const char gmx_unused *deviceOptions,
+ int imdport,
+ unsigned long Flags,
+ gmx_walltime_accounting_t walltime_accounting)
+{
+ gmx_mdoutf_t outf = NULL;
+ gmx_int64_t step, step_rel;
+ double elapsed_time;
+ double t, t0, lam0[efptNR];
+ gmx_bool bGStatEveryStep, bGStat, bCalcVir, bCalcEner;
+ gmx_bool bNS, bNStList, bSimAnn, bStopCM, bRerunMD, bNotLastFrame = FALSE,
+ bFirstStep, bStateFromCP, bStateFromTPX, bInitStep, bLastStep,
+ bBornRadii, bStartingFromCpt;
+ gmx_bool bDoDHDL = FALSE, bDoFEP = FALSE, bDoExpanded = FALSE;
+ gmx_bool do_ene, do_log, do_verbose, bRerunWarnNoV = TRUE,
+ bForceUpdate = FALSE, bCPT;
+ gmx_bool bMasterState;
+ int force_flags, cglo_flags;
+ tensor force_vir, shake_vir, total_vir, tmp_vir, pres;
+ int i, m;
+ t_trxstatus *status;
+ rvec mu_tot;
+ t_vcm *vcm;
+ t_state *bufstate = NULL;
+ matrix *scale_tot, pcoupl_mu, M, ebox;
+ gmx_nlheur_t nlh;
+ t_trxframe rerun_fr;
+ gmx_repl_ex_t repl_ex = NULL;
+ int nchkpt = 1;
+ gmx_localtop_t *top;
+ t_mdebin *mdebin = NULL;
+ t_state *state = NULL;
+ rvec *f_global = NULL;
+ gmx_enerdata_t *enerd;
+ rvec *f = NULL;
+ gmx_global_stat_t gstat;
+ gmx_update_t upd = NULL;
+ t_graph *graph = NULL;
+ globsig_t gs;
+ gmx_groups_t *groups;
+ gmx_ekindata_t *ekind, *ekind_save;
+ gmx_shellfc_t shellfc;
+ int count, nconverged = 0;
+ real timestep = 0;
+ double tcount = 0;
+ gmx_bool bConverged = TRUE, bOK, bSumEkinhOld, bExchanged, bNeedRepartition;
+ gmx_bool bAppend;
+ gmx_bool bResetCountersHalfMaxH = FALSE;
+ gmx_bool bVV, bIterativeCase, bFirstIterate, bTemp, bPres, bTrotter;
+ gmx_bool bUpdateDoLR;
+ real dvdl_constr;
+ rvec *cbuf = NULL;
+ matrix lastbox;
+ real veta_save, scalevir, tracevir;
+ real vetanew = 0;
+ int lamnew = 0;
+ /* for FEP */
+ int nstfep;
+ double cycles;
+ real saved_conserved_quantity = 0;
+ real last_ekin = 0;
+ int iter_i;
+ t_extmass MassQ;
+ int **trotter_seq;
+ char sbuf[STEPSTRSIZE], sbuf2[STEPSTRSIZE];
+ int handled_stop_condition = gmx_stop_cond_none; /* compare to get_stop_condition*/
+ gmx_iterate_t iterate;
+ gmx_int64_t multisim_nsteps = -1; /* number of steps to do before first multisim
+ simulation stops. If equal to zero, don't
+ communicate any more between multisims.*/
+ /* PME load balancing data for GPU kernels */
+ pme_load_balancing_t pme_loadbal = NULL;
+ double cycles_pmes;
+ gmx_bool bPMETuneTry = FALSE, bPMETuneRunning = FALSE;
+
+ /* Interactive MD */
+ gmx_bool bIMDstep = FALSE;
+
+#ifdef GMX_FAHCORE
+ /* Temporary addition for FAHCORE checkpointing */
+ int chkpt_ret;
+#endif
+
+ /* Check for special mdrun options */
+ bRerunMD = (Flags & MD_RERUN);
+ bAppend = (Flags & MD_APPENDFILES);
+ if (Flags & MD_RESETCOUNTERSHALFWAY)
+ {
+ if (ir->nsteps > 0)
+ {
+ /* Signal to reset the counters half the simulation steps. */
+ wcycle_set_reset_counters(wcycle, ir->nsteps/2);
+ }
+ /* Signal to reset the counters halfway the simulation time. */
+ bResetCountersHalfMaxH = (max_hours > 0);
+ }
+
+ /* md-vv uses averaged full step velocities for T-control
+ md-vv-avek uses averaged half step velocities for T-control (but full step ekin for P control)
+ md uses averaged half step kinetic energies to determine temperature unless defined otherwise by GMX_EKIN_AVE_VEL; */
+ bVV = EI_VV(ir->eI);
+ if (bVV) /* to store the initial velocities while computing virial */
+ {
+ snew(cbuf, top_global->natoms);
+ }
+ /* all the iteratative cases - only if there are constraints */
+ bIterativeCase = ((IR_NPH_TROTTER(ir) || IR_NPT_TROTTER(ir)) && (constr) && (!bRerunMD));
+ gmx_iterate_init(&iterate, FALSE); /* The default value of iterate->bIterationActive is set to
+ false in this step. The correct value, true or false,
+ is set at each step, as it depends on the frequency of temperature
+ and pressure control.*/
+ bTrotter = (bVV && (IR_NPT_TROTTER(ir) || IR_NPH_TROTTER(ir) || IR_NVT_TROTTER(ir)));
+
+ if (bRerunMD)
+ {
+ /* Since we don't know if the frames read are related in any way,
+ * rebuild the neighborlist at every step.
+ */
+ ir->nstlist = 1;
+ ir->nstcalcenergy = 1;
+ nstglobalcomm = 1;
+ }
+
+ check_ir_old_tpx_versions(cr, fplog, ir, top_global);
+
+ nstglobalcomm = check_nstglobalcomm(fplog, cr, nstglobalcomm, ir);
+ bGStatEveryStep = (nstglobalcomm == 1);
+
+ if (!bGStatEveryStep && ir->nstlist == -1 && fplog != NULL)
+ {
+ fprintf(fplog,
+ "To reduce the energy communication with nstlist = -1\n"
+ "the neighbor list validity should not be checked at every step,\n"
+ "this means that exact integration is not guaranteed.\n"
+ "The neighbor list validity is checked after:\n"
+ " <n.list life time> - 2*std.dev.(n.list life time) steps.\n"
+ "In most cases this will result in exact integration.\n"
+ "This reduces the energy communication by a factor of 2 to 3.\n"
+ "If you want less energy communication, set nstlist > 3.\n\n");
+ }
+
+ if (bRerunMD)
+ {
+ ir->nstxout_compressed = 0;
+ }
+ groups = &top_global->groups;
+
+ /* Initial values */
+ init_md(fplog, cr, ir, oenv, &t, &t0, state_global->lambda,
+ &(state_global->fep_state), lam0,
+ nrnb, top_global, &upd,
+ nfile, fnm, &outf, &mdebin,
+ force_vir, shake_vir, mu_tot, &bSimAnn, &vcm, Flags);
+
+ clear_mat(total_vir);
+ clear_mat(pres);
+ /* Energy terms and groups */
+ snew(enerd, 1);
+ init_enerdata(top_global->groups.grps[egcENER].nr, ir->fepvals->n_lambda,
+ enerd);
+ if (DOMAINDECOMP(cr))
+ {
+ f = NULL;
+ }
+ else
+ {
+ snew(f, top_global->natoms);
+ }
+
+ /* Kinetic energy data */
+ snew(ekind, 1);
+ init_ekindata(fplog, top_global, &(ir->opts), ekind);
+ /* needed for iteration of constraints */
+ snew(ekind_save, 1);
+ init_ekindata(fplog, top_global, &(ir->opts), ekind_save);
+ /* Copy the cos acceleration to the groups struct */
+ ekind->cosacc.cos_accel = ir->cos_accel;
+
+ gstat = global_stat_init(ir);
+ debug_gmx();
+
+ /* Check for polarizable models and flexible constraints */
+ shellfc = init_shell_flexcon(fplog, fr->cutoff_scheme == ecutsVERLET,
+ top_global, n_flexible_constraints(constr),
+ (ir->bContinuation ||
+ (DOMAINDECOMP(cr) && !MASTER(cr))) ?
+ NULL : state_global->x);
+
+ if (shellfc && ir->eI == eiNM)
+ {
+ /* Currently shells don't work with Normal Modes */
+ gmx_fatal(FARGS, "Normal Mode analysis is not supported with shells.\nIf you'd like to help with adding support, we have an open discussion at http://redmine.gromacs.org/issues/879\n");
+ }
+
+ if (vsite && ir->eI == eiNM)
+ {
+ /* Currently virtual sites don't work with Normal Modes */
+ gmx_fatal(FARGS, "Normal Mode analysis is not supported with virtual sites.\nIf you'd like to help with adding support, we have an open discussion at http://redmine.gromacs.org/issues/879\n");
+ }
+
+ if (DEFORM(*ir))
+ {
+ tMPI_Thread_mutex_lock(&deform_init_box_mutex);
+ set_deform_reference_box(upd,
+ deform_init_init_step_tpx,
+ deform_init_box_tpx);
+ tMPI_Thread_mutex_unlock(&deform_init_box_mutex);
+ }
+
+ {
+ double io = compute_io(ir, top_global->natoms, groups, mdebin->ebin->nener, 1);
+ if ((io > 2000) && MASTER(cr))
+ {
+ fprintf(stderr,
+ "\nWARNING: This run will generate roughly %.0f Mb of data\n\n",
+ io);
+ }
+ }
+
+ if (DOMAINDECOMP(cr))
+ {
+ top = dd_init_local_top(top_global);
+
+ snew(state, 1);
+ dd_init_local_state(cr->dd, state_global, state);
+
+ if (DDMASTER(cr->dd) && ir->nstfout)
+ {
+ snew(f_global, state_global->natoms);
+ }
+ }
+ else
+ {
+ top = gmx_mtop_generate_local_top(top_global, ir);
+
+ forcerec_set_excl_load(fr, top);
+
+ state = serial_init_local_state(state_global);
+ f_global = f;
+
+ atoms2md(top_global, ir, 0, NULL, top_global->natoms, mdatoms);
+
+ if (vsite)
+ {
+ set_vsite_top(vsite, top, mdatoms, cr);
+ }
+
+ if (ir->ePBC != epbcNONE && !fr->bMolPBC)
+ {
+ graph = mk_graph(fplog, &(top->idef), 0, top_global->natoms, FALSE, FALSE);
+ }
+
+ if (shellfc)
+ {
+ make_local_shells(cr, mdatoms, shellfc);
+ }
+
+ setup_bonded_threading(fr, &top->idef);
+ }
+
+ /* Set up interactive MD (IMD) */
+ init_IMD(ir, cr, top_global, fplog, ir->nstcalcenergy, state_global->x,
+ nfile, fnm, oenv, imdport, Flags);
+
+ if (DOMAINDECOMP(cr))
+ {
+ /* Distribute the charge groups over the nodes from the master node */
+ dd_partition_system(fplog, ir->init_step, cr, TRUE, 1,
+ state_global, top_global, ir,
+ state, &f, mdatoms, top, fr,
+ vsite, shellfc, constr,
+ nrnb, wcycle, FALSE);
+
+ }
+
+ update_mdatoms(mdatoms, state->lambda[efptMASS]);
+
+ if (opt2bSet("-cpi", nfile, fnm))
+ {
+ bStateFromCP = gmx_fexist_master(opt2fn_master("-cpi", nfile, fnm, cr), cr);
+ }
+ else
+ {
+ bStateFromCP = FALSE;
+ }
+
+ if (ir->bExpanded)
+ {
+ init_expanded_ensemble(bStateFromCP, ir, &state->dfhist);
+ }
+
+ if (MASTER(cr))
+ {
+ if (bStateFromCP)
+ {
+ /* Update mdebin with energy history if appending to output files */
+ if (Flags & MD_APPENDFILES)
+ {
+ restore_energyhistory_from_state(mdebin, &state_global->enerhist);
+ }
+ else
+ {
+ /* We might have read an energy history from checkpoint,
+ * free the allocated memory and reset the counts.
+ */
+ done_energyhistory(&state_global->enerhist);
+ init_energyhistory(&state_global->enerhist);
+ }
+ }
+ /* Set the initial energy history in state by updating once */
+ update_energyhistory(&state_global->enerhist, mdebin);
+ }
+
+ /* Initialize constraints */
+ if (constr && !DOMAINDECOMP(cr))
+ {
+ set_constraints(constr, top, ir, mdatoms, cr);
+ }
+
+ if (repl_ex_nst > 0)
+ {
+ /* We need to be sure replica exchange can only occur
+ * when the energies are current */
+ check_nst_param(fplog, cr, "nstcalcenergy", ir->nstcalcenergy,
+ "repl_ex_nst", &repl_ex_nst);
+ /* This check needs to happen before inter-simulation
+ * signals are initialized, too */
+ }
+ if (repl_ex_nst > 0 && MASTER(cr))
+ {
+ repl_ex = init_replica_exchange(fplog, cr->ms, state_global, ir,
+ repl_ex_nst, repl_ex_nex, repl_ex_seed);
+ }
+
+ /* PME tuning is only supported with GPUs or PME nodes and not with rerun.
+ * PME tuning is not supported with PME only for LJ and not for Coulomb.
+ */
+ if ((Flags & MD_TUNEPME) &&
+ EEL_PME(fr->eeltype) &&
+ ( (fr->cutoff_scheme == ecutsVERLET && fr->nbv->bUseGPU) || !(cr->duty & DUTY_PME)) &&
+ !bRerunMD)
+ {
+ pme_loadbal_init(&pme_loadbal, ir, state->box, fr->ic, fr->pmedata);
+ cycles_pmes = 0;
+ if (cr->duty & DUTY_PME)
+ {
+ /* Start tuning right away, as we can't measure the load */
+ bPMETuneRunning = TRUE;
+ }
+ else
+ {
+ /* Separate PME nodes, we can measure the PP/PME load balance */
+ bPMETuneTry = TRUE;
+ }
+ }
+
+ if (!ir->bContinuation && !bRerunMD)
+ {
+ if (mdatoms->cFREEZE && (state->flags & (1<<estV)))
+ {
+ /* Set the velocities of frozen particles to zero */
+ for (i = 0; i < mdatoms->homenr; i++)
+ {
+ for (m = 0; m < DIM; m++)
+ {
+ if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
+ {
+ state->v[i][m] = 0;
+ }
+ }
+ }
+ }
+
+ if (constr)
+ {
+ /* Constrain the initial coordinates and velocities */
+ do_constrain_first(fplog, constr, ir, mdatoms, state,
+ cr, nrnb, fr, top);
+ }
+ if (vsite)
+ {
+ /* Construct the virtual sites for the initial configuration */
+ construct_vsites(vsite, state->x, ir->delta_t, NULL,
+ top->idef.iparams, top->idef.il,
+ fr->ePBC, fr->bMolPBC, cr, state->box);
+ }
+ }
+
+ debug_gmx();
+
+ /* set free energy calculation frequency as the minimum
+ greatest common denominator of nstdhdl, nstexpanded, and repl_ex_nst*/
+ nstfep = ir->fepvals->nstdhdl;
+ if (ir->bExpanded)
+ {
+ nstfep = gmx_greatest_common_divisor(ir->fepvals->nstdhdl, nstfep);
+ }
+ if (repl_ex_nst > 0)
+ {
+ nstfep = gmx_greatest_common_divisor(repl_ex_nst, nstfep);
+ }
+
+ /* I'm assuming we need global communication the first time! MRS */
+ cglo_flags = (CGLO_TEMPERATURE | CGLO_GSTAT
+ | ((ir->comm_mode != ecmNO) ? CGLO_STOPCM : 0)
+ | (bVV ? CGLO_PRESSURE : 0)
+ | (bVV ? CGLO_CONSTRAINT : 0)
+ | (bRerunMD ? CGLO_RERUNMD : 0)
+ | ((Flags & MD_READ_EKIN) ? CGLO_READEKIN : 0));
+
+ bSumEkinhOld = FALSE;
+ compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm,
+ NULL, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
+ constr, NULL, FALSE, state->box,
+ top_global, &bSumEkinhOld, cglo_flags);
+ if (ir->eI == eiVVAK)
+ {
+ /* a second call to get the half step temperature initialized as well */
+ /* we do the same call as above, but turn the pressure off -- internally to
+ compute_globals, this is recognized as a velocity verlet half-step
+ kinetic energy calculation. This minimized excess variables, but
+ perhaps loses some logic?*/
+
+ compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm,
+ NULL, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
+ constr, NULL, FALSE, state->box,
+ top_global, &bSumEkinhOld,
+ cglo_flags &~(CGLO_STOPCM | CGLO_PRESSURE));
+ }
+
+ /* Calculate the initial half step temperature, and save the ekinh_old */
+ if (!(Flags & MD_STARTFROMCPT))
+ {
+ for (i = 0; (i < ir->opts.ngtc); i++)
+ {
+ copy_mat(ekind->tcstat[i].ekinh, ekind->tcstat[i].ekinh_old);
+ }
+ }
+ if (ir->eI != eiVV)
+ {
+ enerd->term[F_TEMP] *= 2; /* result of averages being done over previous and current step,
+ and there is no previous step */
+ }
+
+ /* if using an iterative algorithm, we need to create a working directory for the state. */
+ if (bIterativeCase)
+ {
+ bufstate = init_bufstate(state);
+ }
+
+ /* need to make an initiation call to get the Trotter variables set, as well as other constants for non-trotter
+ temperature control */
+ trotter_seq = init_npt_vars(ir, state, &MassQ, bTrotter);
+
+ if (MASTER(cr))
+ {
+ if (constr && !ir->bContinuation && ir->eConstrAlg == econtLINCS)
+ {
+ fprintf(fplog,
+ "RMS relative constraint deviation after constraining: %.2e\n",
+ constr_rmsd(constr, FALSE));
+ }
+ if (EI_STATE_VELOCITY(ir->eI))
+ {
+ fprintf(fplog, "Initial temperature: %g K\n", enerd->term[F_TEMP]);
+ }
+ if (bRerunMD)
+ {
+ fprintf(stderr, "starting md rerun '%s', reading coordinates from"
+ " input trajectory '%s'\n\n",
+ *(top_global->name), opt2fn("-rerun", nfile, fnm));
+ if (bVerbose)
+ {
+ fprintf(stderr, "Calculated time to finish depends on nsteps from "
+ "run input file,\nwhich may not correspond to the time "
+ "needed to process input trajectory.\n\n");
+ }
+ }
+ else
+ {
+ char tbuf[20];
+ fprintf(stderr, "starting mdrun '%s'\n",
+ *(top_global->name));
+ if (ir->nsteps >= 0)
+ {
+ sprintf(tbuf, "%8.1f", (ir->init_step+ir->nsteps)*ir->delta_t);
+ }
+ else
+ {
+ sprintf(tbuf, "%s", "infinite");
+ }
+ if (ir->init_step > 0)
+ {
+ fprintf(stderr, "%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
+ gmx_step_str(ir->init_step+ir->nsteps, sbuf), tbuf,
+ gmx_step_str(ir->init_step, sbuf2),
+ ir->init_step*ir->delta_t);
+ }
+ else
+ {
+ fprintf(stderr, "%s steps, %s ps.\n",
+ gmx_step_str(ir->nsteps, sbuf), tbuf);
+ }
+ }
+ fprintf(fplog, "\n");
+ }
+
+ walltime_accounting_start(walltime_accounting);
+ wallcycle_start(wcycle, ewcRUN);
+ print_start(fplog, cr, walltime_accounting, "mdrun");
+
+ /* safest point to do file checkpointing is here. More general point would be immediately before integrator call */
+#ifdef GMX_FAHCORE
+ chkpt_ret = fcCheckPointParallel( cr->nodeid,
+ NULL, 0);
+ if (chkpt_ret == 0)
+ {
+ gmx_fatal( 3, __FILE__, __LINE__, "Checkpoint error on step %d\n", 0 );
+ }
+#endif
+
+ debug_gmx();
+ /***********************************************************
+ *
+ * Loop over MD steps
+ *
+ ************************************************************/
+
+ /* if rerunMD then read coordinates and velocities from input trajectory */
+ if (bRerunMD)
+ {
+ if (getenv("GMX_FORCE_UPDATE"))
+ {
+ bForceUpdate = TRUE;
+ }
+
+ rerun_fr.natoms = 0;
+ if (MASTER(cr))
+ {
+ bNotLastFrame = read_first_frame(oenv, &status,
+ opt2fn("-rerun", nfile, fnm),
+ &rerun_fr, TRX_NEED_X | TRX_READ_V);
+ if (rerun_fr.natoms != top_global->natoms)
+ {
+ gmx_fatal(FARGS,
+ "Number of atoms in trajectory (%d) does not match the "
+ "run input file (%d)\n",
+ rerun_fr.natoms, top_global->natoms);
+ }
+ if (ir->ePBC != epbcNONE)
+ {
+ if (!rerun_fr.bBox)
+ {
+ gmx_fatal(FARGS, "Rerun trajectory frame step %d time %f does not contain a box, while pbc is used", rerun_fr.step, rerun_fr.time);
+ }
+ if (max_cutoff2(ir->ePBC, rerun_fr.box) < sqr(fr->rlistlong))
+ {
+ gmx_fatal(FARGS, "Rerun trajectory frame step %d time %f has too small box dimensions", rerun_fr.step, rerun_fr.time);
+ }
+ }
+ }
+
+ if (PAR(cr))
+ {
+ rerun_parallel_comm(cr, &rerun_fr, &bNotLastFrame);
+ }
+
+ if (ir->ePBC != epbcNONE)
+ {
+ /* Set the shift vectors.
+ * Necessary here when have a static box different from the tpr box.
+ */
+ calc_shifts(rerun_fr.box, fr->shift_vec);
+ }
+ }
+
+ /* loop over MD steps or if rerunMD to end of input trajectory */
+ bFirstStep = TRUE;
+ /* Skip the first Nose-Hoover integration when we get the state from tpx */
+ bStateFromTPX = !bStateFromCP;
+ bInitStep = bFirstStep && (bStateFromTPX || bVV);
+ bStartingFromCpt = (Flags & MD_STARTFROMCPT) && bInitStep;
+ bLastStep = FALSE;
+ bSumEkinhOld = FALSE;
+ bExchanged = FALSE;
+ bNeedRepartition = FALSE;
+
+ init_global_signals(&gs, cr, ir, repl_ex_nst);
+
+ step = ir->init_step;
+ step_rel = 0;
+
+ if (ir->nstlist == -1)
+ {
+ init_nlistheuristics(&nlh, bGStatEveryStep, step);
+ }
+
+ if (MULTISIM(cr) && (repl_ex_nst <= 0 ))
+ {
+ /* check how many steps are left in other sims */
+ multisim_nsteps = get_multisim_nsteps(cr, ir->nsteps);
+ }
+
+
+ /* and stop now if we should */
+ bLastStep = (bRerunMD || (ir->nsteps >= 0 && step_rel > ir->nsteps) ||
+ ((multisim_nsteps >= 0) && (step_rel >= multisim_nsteps )));
+ while (!bLastStep || (bRerunMD && bNotLastFrame))
+ {
+
+ wallcycle_start(wcycle, ewcSTEP);
+
+ if (bRerunMD)
+ {
+ if (rerun_fr.bStep)
+ {
+ step = rerun_fr.step;
+ step_rel = step - ir->init_step;
+ }
+ if (rerun_fr.bTime)
+ {
+ t = rerun_fr.time;
+ }
+ else
+ {
+ t = step;
+ }
+ }
+ else
+ {
+ bLastStep = (step_rel == ir->nsteps);
+ t = t0 + step*ir->delta_t;
+ }
+
+ if (ir->efep != efepNO || ir->bSimTemp)
+ {
+ /* find and set the current lambdas. If rerunning, we either read in a state, or a lambda value,
+ requiring different logic. */
+
+ set_current_lambdas(step, ir->fepvals, bRerunMD, &rerun_fr, state_global, state, lam0);
+ bDoDHDL = do_per_step(step, ir->fepvals->nstdhdl);
+ bDoFEP = (do_per_step(step, nstfep) && (ir->efep != efepNO));
+ bDoExpanded = (do_per_step(step, ir->expandedvals->nstexpanded)
+ && (ir->bExpanded) && (step > 0) && (!bStartingFromCpt));
+ }
+
+ if (bSimAnn)
+ {
+ update_annealing_target_temp(&(ir->opts), t);
+ }
+
+ if (bRerunMD)
+ {
+ if (!DOMAINDECOMP(cr) || MASTER(cr))
+ {
+ for (i = 0; i < state_global->natoms; i++)
+ {
+ copy_rvec(rerun_fr.x[i], state_global->x[i]);
+ }
+ if (rerun_fr.bV)
+ {
+ for (i = 0; i < state_global->natoms; i++)
+ {
+ copy_rvec(rerun_fr.v[i], state_global->v[i]);
+ }
+ }
+ else
+ {
+ for (i = 0; i < state_global->natoms; i++)
+ {
+ clear_rvec(state_global->v[i]);
+ }
+ if (bRerunWarnNoV)
+ {
+ fprintf(stderr, "\nWARNING: Some frames do not contain velocities.\n"
+ " Ekin, temperature and pressure are incorrect,\n"
+ " the virial will be incorrect when constraints are present.\n"
+ "\n");
+ bRerunWarnNoV = FALSE;
+ }
+ }
+ }
+ copy_mat(rerun_fr.box, state_global->box);
+ copy_mat(state_global->box, state->box);
+
+ if (vsite && (Flags & MD_RERUN_VSITE))
+ {
+ if (DOMAINDECOMP(cr))
+ {
+ gmx_fatal(FARGS, "Vsite recalculation with -rerun is not implemented with domain decomposition, use a single rank");
+ }
+ if (graph)
+ {
+ /* Following is necessary because the graph may get out of sync
+ * with the coordinates if we only have every N'th coordinate set
+ */
+ mk_mshift(fplog, graph, fr->ePBC, state->box, state->x);
+ shift_self(graph, state->box, state->x);
+ }
+ construct_vsites(vsite, state->x, ir->delta_t, state->v,
+ top->idef.iparams, top->idef.il,
+ fr->ePBC, fr->bMolPBC, cr, state->box);
+ if (graph)
+ {
+ unshift_self(graph, state->box, state->x);
+ }
+ }
+ }
+
+ /* Stop Center of Mass motion */
+ bStopCM = (ir->comm_mode != ecmNO && do_per_step(step, ir->nstcomm));
+
+ if (bRerunMD)
+ {
+ /* for rerun MD always do Neighbour Searching */
+ bNS = (bFirstStep || ir->nstlist != 0);
+ bNStList = bNS;
+ }
+ else
+ {
+ /* Determine whether or not to do Neighbour Searching and LR */
+ bNStList = (ir->nstlist > 0 && step % ir->nstlist == 0);
+
+ bNS = (bFirstStep || bExchanged || bNeedRepartition || bNStList || bDoFEP ||
+ (ir->nstlist == -1 && nlh.nabnsb > 0));
+
+ if (bNS && ir->nstlist == -1)
+ {
+ set_nlistheuristics(&nlh, bFirstStep || bExchanged || bNeedRepartition || bDoFEP, step);
+ }
+ }
+
+ /* check whether we should stop because another simulation has
+ stopped. */
+ if (MULTISIM(cr))
+ {
+ if ( (multisim_nsteps >= 0) && (step_rel >= multisim_nsteps) &&
+ (multisim_nsteps != ir->nsteps) )
+ {
+ if (bNS)
+ {
+ if (MASTER(cr))
+ {
+ fprintf(stderr,
+ "Stopping simulation %d because another one has finished\n",
+ cr->ms->sim);
+ }
+ bLastStep = TRUE;
+ gs.sig[eglsCHKPT] = 1;
+ }
+ }
+ }
+
+ /* < 0 means stop at next step, > 0 means stop at next NS step */
+ if ( (gs.set[eglsSTOPCOND] < 0) ||
+ ( (gs.set[eglsSTOPCOND] > 0) && (bNStList || ir->nstlist == 0) ) )
+ {
+ bLastStep = TRUE;
+ }
+
+ /* Determine whether or not to update the Born radii if doing GB */
+ bBornRadii = bFirstStep;
+ if (ir->implicit_solvent && (step % ir->nstgbradii == 0))
+ {
+ bBornRadii = TRUE;
+ }
+
+ do_log = do_per_step(step, ir->nstlog) || bFirstStep || bLastStep;
+ do_verbose = bVerbose &&
+ (step % stepout == 0 || bFirstStep || bLastStep);
+
+ if (bNS && !(bFirstStep && ir->bContinuation && !bRerunMD))
+ {
+ if (bRerunMD)
+ {
+ bMasterState = TRUE;
+ }
+ else
+ {
+ bMasterState = FALSE;
+ /* Correct the new box if it is too skewed */
+ if (DYNAMIC_BOX(*ir))
+ {
+ if (correct_box(fplog, step, state->box, graph))
+ {
+ bMasterState = TRUE;
+ }
+ }
+ if (DOMAINDECOMP(cr) && bMasterState)
+ {
+ dd_collect_state(cr->dd, state, state_global);
+ }
+ }
+
+ if (DOMAINDECOMP(cr))
+ {
+ /* Repartition the domain decomposition */
+ wallcycle_start(wcycle, ewcDOMDEC);
+ dd_partition_system(fplog, step, cr,
+ bMasterState, nstglobalcomm,
+ state_global, top_global, ir,
+ state, &f, mdatoms, top, fr,
+ vsite, shellfc, constr,
+ nrnb, wcycle,
+ do_verbose && !bPMETuneRunning);
+ wallcycle_stop(wcycle, ewcDOMDEC);
+ /* If using an iterative integrator, reallocate space to match the decomposition */
+ }
+ }
+
+ if (MASTER(cr) && do_log)
+ {
+ print_ebin_header(fplog, step, t, state->lambda[efptFEP]); /* can we improve the information printed here? */
+ }
+
+ if (ir->efep != efepNO)
+ {
+ update_mdatoms(mdatoms, state->lambda[efptMASS]);
+ }
+
+ if ((bRerunMD && rerun_fr.bV) || bExchanged)
+ {
+
+ /* We need the kinetic energy at minus the half step for determining
+ * the full step kinetic energy and possibly for T-coupling.*/
+ /* This may not be quite working correctly yet . . . . */
+ compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm,
+ wcycle, enerd, NULL, NULL, NULL, NULL, mu_tot,
+ constr, NULL, FALSE, state->box,
+ top_global, &bSumEkinhOld,
+ CGLO_RERUNMD | CGLO_GSTAT | CGLO_TEMPERATURE);
+ }
+ clear_mat(force_vir);
+
+ /* We write a checkpoint at this MD step when:
+ * either at an NS step when we signalled through gs,
+ * or at the last step (but not when we do not want confout),
+ * but never at the first step or with rerun.
+ */
+ bCPT = (((gs.set[eglsCHKPT] && (bNS || ir->nstlist == 0)) ||
+ (bLastStep && (Flags & MD_CONFOUT))) &&
+ step > ir->init_step && !bRerunMD);
+ if (bCPT)
+ {
+ gs.set[eglsCHKPT] = 0;
+ }
+
+ /* Determine the energy and pressure:
+ * at nstcalcenergy steps and at energy output steps (set below).
+ */
+ if (EI_VV(ir->eI) && (!bInitStep))
+ {
+ /* for vv, the first half of the integration actually corresponds
+ to the previous step. bCalcEner is only required to be evaluated on the 'next' step,
+ but the virial needs to be calculated on both the current step and the 'next' step. Future
+ reorganization may be able to get rid of one of the bCalcVir=TRUE steps. */
+
+ bCalcEner = do_per_step(step-1, ir->nstcalcenergy);
+ bCalcVir = bCalcEner ||
+ (ir->epc != epcNO && (do_per_step(step, ir->nstpcouple) || do_per_step(step-1, ir->nstpcouple)));
+ }
+ else
+ {
+ bCalcEner = do_per_step(step, ir->nstcalcenergy);
+ bCalcVir = bCalcEner ||
+ (ir->epc != epcNO && do_per_step(step, ir->nstpcouple));
+ }
+
+ /* Do we need global communication ? */
+ bGStat = (bCalcVir || bCalcEner || bStopCM ||
+ do_per_step(step, nstglobalcomm) || (bVV && IR_NVT_TROTTER(ir) && do_per_step(step-1, nstglobalcomm)) ||
+ (ir->nstlist == -1 && !bRerunMD && step >= nlh.step_nscheck));
+
+ do_ene = (do_per_step(step, ir->nstenergy) || bLastStep);
+
+ if (do_ene || do_log)
+ {
+ bCalcVir = TRUE;
+ bCalcEner = TRUE;
+ bGStat = TRUE;
+ }
+
+ /* these CGLO_ options remain the same throughout the iteration */
+ cglo_flags = ((bRerunMD ? CGLO_RERUNMD : 0) |
+ (bGStat ? CGLO_GSTAT : 0)
+ );
+
+ force_flags = (GMX_FORCE_STATECHANGED |
+ ((DYNAMIC_BOX(*ir) || bRerunMD) ? GMX_FORCE_DYNAMICBOX : 0) |
+ GMX_FORCE_ALLFORCES |
+ GMX_FORCE_SEPLRF |
+ (bCalcVir ? GMX_FORCE_VIRIAL : 0) |
+ (bCalcEner ? GMX_FORCE_ENERGY : 0) |
+ (bDoFEP ? GMX_FORCE_DHDL : 0)
+ );
+
+ if (fr->bTwinRange)
+ {
+ if (do_per_step(step, ir->nstcalclr))
+ {
+ force_flags |= GMX_FORCE_DO_LR;
+ }
+ }
+
+ if (shellfc)
+ {
+ /* Now is the time to relax the shells */
+ count = relax_shell_flexcon(fplog, cr, bVerbose, step,
+ ir, bNS, force_flags,
+ top,
+ constr, enerd, fcd,
+ state, f, force_vir, mdatoms,
+ nrnb, wcycle, graph, groups,
+ shellfc, fr, bBornRadii, t, mu_tot,
+ &bConverged, vsite,
+ mdoutf_get_fp_field(outf));
+ tcount += count;
+
+ if (bConverged)
+ {
+ nconverged++;
+ }
+ }
+ else
+ {
+ /* The coordinates (x) are shifted (to get whole molecules)
+ * in do_force.
+ * This is parallellized as well, and does communication too.
+ * Check comments in sim_util.c
+ */
+ do_force(fplog, cr, ir, step, nrnb, wcycle, top, groups,
+ state->box, state->x, &state->hist,
+ f, force_vir, mdatoms, enerd, fcd,
+ state->lambda, graph,
+ fr, vsite, mu_tot, t, mdoutf_get_fp_field(outf), ed, bBornRadii,
+ (bNS ? GMX_FORCE_NS : 0) | force_flags);
+ }
+
+ if (bVV && !bStartingFromCpt && !bRerunMD)
+ /* ############### START FIRST UPDATE HALF-STEP FOR VV METHODS############### */
+ {
+ if (ir->eI == eiVV && bInitStep)
+ {
+ /* if using velocity verlet with full time step Ekin,
+ * take the first half step only to compute the
+ * virial for the first step. From there,
+ * revert back to the initial coordinates
+ * so that the input is actually the initial step.
+ */
+ copy_rvecn(state->v, cbuf, 0, state->natoms); /* should make this better for parallelizing? */
+ }
+ else
+ {
+ /* this is for NHC in the Ekin(t+dt/2) version of vv */
+ trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ1);
+ }
+
+ /* If we are using twin-range interactions where the long-range component
+ * is only evaluated every nstcalclr>1 steps, we should do a special update
+ * step to combine the long-range forces on these steps.
+ * For nstcalclr=1 this is not done, since the forces would have been added
+ * directly to the short-range forces already.
+ *
+ * TODO Remove various aspects of VV+twin-range in master
+ * branch, because VV integrators did not ever support
+ * twin-range multiple time stepping with constraints.
+ */
+ bUpdateDoLR = (fr->bTwinRange && do_per_step(step, ir->nstcalclr));
+
+ update_coords(fplog, step, ir, mdatoms, state, fr->bMolPBC,
- bUpdateDoLR, fr->f_twin, fcd,
++ f, bUpdateDoLR, fr->f_twin, bCalcVir ? &fr->vir_twin_constr : NULL, fcd,
+ ekind, M, upd, bInitStep, etrtVELOCITY1,
+ cr, nrnb, constr, &top->idef);
+
+ if (bIterativeCase && do_per_step(step-1, ir->nstpcouple) && !bInitStep)
+ {
+ gmx_iterate_init(&iterate, TRUE);
+ }
+ /* for iterations, we save these vectors, as we will be self-consistently iterating
+ the calculations */
+
+ /*#### UPDATE EXTENDED VARIABLES IN TROTTER FORMULATION */
+
+ /* save the state */
+ if (iterate.bIterationActive)
+ {
+ copy_coupling_state(state, bufstate, ekind, ekind_save, &(ir->opts));
+ }
+
+ bFirstIterate = TRUE;
+ while (bFirstIterate || iterate.bIterationActive)
+ {
+ if (iterate.bIterationActive)
+ {
+ copy_coupling_state(bufstate, state, ekind_save, ekind, &(ir->opts));
+ if (bFirstIterate && bTrotter)
+ {
+ /* The first time through, we need a decent first estimate
+ of veta(t+dt) to compute the constraints. Do
+ this by computing the box volume part of the
+ trotter integration at this time. Nothing else
+ should be changed by this routine here. If
+ !(first time), we start with the previous value
+ of veta. */
+
+ veta_save = state->veta;
+ trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ0);
+ vetanew = state->veta;
+ state->veta = veta_save;
+ }
+ }
+
+ bOK = TRUE;
+ if (!bRerunMD || rerun_fr.bV || bForceUpdate) /* Why is rerun_fr.bV here? Unclear. */
+ {
+ update_constraints(fplog, step, NULL, ir, ekind, mdatoms,
+ state, fr->bMolPBC, graph, f,
+ &top->idef, shake_vir,
+ cr, nrnb, wcycle, upd, constr,
+ TRUE, bCalcVir, vetanew);
+
++ if (bCalcVir && bUpdateDoLR && ir->nstcalclr > 1)
++ {
++ /* Correct the virial for multiple time stepping */
++ m_sub(shake_vir, fr->vir_twin_constr, shake_vir);
++ }
++
+ if (!bOK)
+ {
+ gmx_fatal(FARGS, "Constraint error: Shake, Lincs or Settle could not solve the constrains");
+ }
+
+ }
+ else if (graph)
+ {
+ /* Need to unshift here if a do_force has been
+ called in the previous step */
+ unshift_self(graph, state->box, state->x);
+ }
+
+ /* if VV, compute the pressure and constraints */
+ /* For VV2, we strictly only need this if using pressure
+ * control, but we really would like to have accurate pressures
+ * printed out.
+ * Think about ways around this in the future?
+ * For now, keep this choice in comments.
+ */
+ /*bPres = (ir->eI==eiVV || IR_NPT_TROTTER(ir)); */
+ /*bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && IR_NPT_TROTTER(ir)));*/
+ bPres = TRUE;
+ bTemp = ((ir->eI == eiVV && (!bInitStep)) || (ir->eI == eiVVAK));
+ if (bCalcEner && ir->eI == eiVVAK) /*MRS: 7/9/2010 -- this still doesn't fix it?*/
+ {
+ bSumEkinhOld = TRUE;
+ }
+ /* for vv, the first half of the integration actually corresponds to the previous step.
+ So we need information from the last step in the first half of the integration */
+ if (bGStat || do_per_step(step-1, nstglobalcomm))
+ {
+ compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm,
+ wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
+ constr, NULL, FALSE, state->box,
+ top_global, &bSumEkinhOld,
+ cglo_flags
+ | CGLO_ENERGY
+ | (bTemp ? CGLO_TEMPERATURE : 0)
+ | (bPres ? CGLO_PRESSURE : 0)
+ | (bPres ? CGLO_CONSTRAINT : 0)
+ | ((iterate.bIterationActive) ? CGLO_ITERATE : 0)
+ | (bFirstIterate ? CGLO_FIRSTITERATE : 0)
+ | CGLO_SCALEEKIN
+ );
+ /* explanation of above:
+ a) We compute Ekin at the full time step
+ if 1) we are using the AveVel Ekin, and it's not the
+ initial step, or 2) if we are using AveEkin, but need the full
+ time step kinetic energy for the pressure (always true now, since we want accurate statistics).
+ b) If we are using EkinAveEkin for the kinetic energy for the temperature control, we still feed in
+ EkinAveVel because it's needed for the pressure */
+ }
+ /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
+ if (!bInitStep)
+ {
+ if (bTrotter)
+ {
+ m_add(force_vir, shake_vir, total_vir); /* we need the un-dispersion corrected total vir here */
+ trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ2);
+ }
+ else
+ {
+ if (bExchanged)
+ {
+
+ /* We need the kinetic energy at minus the half step for determining
+ * the full step kinetic energy and possibly for T-coupling.*/
+ /* This may not be quite working correctly yet . . . . */
+ compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm,
+ wcycle, enerd, NULL, NULL, NULL, NULL, mu_tot,
+ constr, NULL, FALSE, state->box,
+ top_global, &bSumEkinhOld,
+ CGLO_RERUNMD | CGLO_GSTAT | CGLO_TEMPERATURE);
+ }
+ }
+ }
+
+ if (iterate.bIterationActive &&
+ done_iterating(cr, fplog, step, &iterate, bFirstIterate,
+ state->veta, &vetanew))
+ {
+ break;
+ }
+ bFirstIterate = FALSE;
+ }
+
+ if (bTrotter && !bInitStep)
+ {
+ copy_mat(shake_vir, state->svir_prev);
+ copy_mat(force_vir, state->fvir_prev);
+ if (IR_NVT_TROTTER(ir) && ir->eI == eiVV)
+ {
+ /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
+ enerd->term[F_TEMP] = sum_ekin(&(ir->opts), ekind, NULL, (ir->eI == eiVV), FALSE);
+ enerd->term[F_EKIN] = trace(ekind->ekin);
+ }
+ }
+ /* if it's the initial step, we performed this first step just to get the constraint virial */
+ if (bInitStep && ir->eI == eiVV)
+ {
+ copy_rvecn(cbuf, state->v, 0, state->natoms);
+ }
+ }
+
+ /* MRS -- now done iterating -- compute the conserved quantity */
+ if (bVV)
+ {
+ saved_conserved_quantity = compute_conserved_from_auxiliary(ir, state, &MassQ);
+ if (ir->eI == eiVV)
+ {
+ last_ekin = enerd->term[F_EKIN];
+ }
+ if ((ir->eDispCorr != edispcEnerPres) && (ir->eDispCorr != edispcAllEnerPres))
+ {
+ saved_conserved_quantity -= enerd->term[F_DISPCORR];
+ }
+ /* sum up the foreign energy and dhdl terms for vv. currently done every step so that dhdl is correct in the .edr */
+ if (!bRerunMD)
+ {
+ sum_dhdl(enerd, state->lambda, ir->fepvals);
+ }
+ }
+
+ /* ######## END FIRST UPDATE STEP ############## */
+ /* ######## If doing VV, we now have v(dt) ###### */
+ if (bDoExpanded)
+ {
+ /* perform extended ensemble sampling in lambda - we don't
+ actually move to the new state before outputting
+ statistics, but if performing simulated tempering, we
+ do update the velocities and the tau_t. */
+
+ lamnew = ExpandedEnsembleDynamics(fplog, ir, enerd, state, &MassQ, state->fep_state, &state->dfhist, step, state->v, mdatoms);
+ /* history is maintained in state->dfhist, but state_global is what is sent to trajectory and log output */
+ copy_df_history(&state_global->dfhist, &state->dfhist);
+ }
+
+ /* Now we have the energies and forces corresponding to the
+ * coordinates at time t. We must output all of this before
+ * the update.
+ */
+ do_md_trajectory_writing(fplog, cr, nfile, fnm, step, step_rel, t,
+ ir, state, state_global, top_global, fr,
+ outf, mdebin, ekind, f, f_global,
+ wcycle, &nchkpt,
+ bCPT, bRerunMD, bLastStep, (Flags & MD_CONFOUT),
+ bSumEkinhOld);
+ /* Check if IMD step and do IMD communication, if bIMD is TRUE. */
+ bIMDstep = do_IMD(ir->bIMD, step, cr, bNS, state->box, state->x, ir, t, wcycle);
+
+ /* kludge -- virial is lost with restart for NPT control. Must restart */
+ if (bStartingFromCpt && bVV)
+ {
+ copy_mat(state->svir_prev, shake_vir);
+ copy_mat(state->fvir_prev, force_vir);
+ }
+
+ elapsed_time = walltime_accounting_get_current_elapsed_time(walltime_accounting);
+
+ /* Check whether everything is still allright */
+ if (((int)gmx_get_stop_condition() > handled_stop_condition)
+#ifdef GMX_THREAD_MPI
+ && MASTER(cr)
+#endif
+ )
+ {
+ /* this is just make gs.sig compatible with the hack
+ of sending signals around by MPI_Reduce with together with
+ other floats */
+ if (gmx_get_stop_condition() == gmx_stop_cond_next_ns)
+ {
+ gs.sig[eglsSTOPCOND] = 1;
+ }
+ if (gmx_get_stop_condition() == gmx_stop_cond_next)
+ {
+ gs.sig[eglsSTOPCOND] = -1;
+ }
+ /* < 0 means stop at next step, > 0 means stop at next NS step */
+ if (fplog)
+ {
+ fprintf(fplog,
+ "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
+ gmx_get_signal_name(),
+ gs.sig[eglsSTOPCOND] == 1 ? "NS " : "");
+ fflush(fplog);
+ }
+ fprintf(stderr,
+ "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
+ gmx_get_signal_name(),
+ gs.sig[eglsSTOPCOND] == 1 ? "NS " : "");
+ fflush(stderr);
+ handled_stop_condition = (int)gmx_get_stop_condition();
+ }
+ else if (MASTER(cr) && (bNS || ir->nstlist <= 0) &&
+ (max_hours > 0 && elapsed_time > max_hours*60.0*60.0*0.99) &&
+ gs.sig[eglsSTOPCOND] == 0 && gs.set[eglsSTOPCOND] == 0)
+ {
+ /* Signal to terminate the run */
+ gs.sig[eglsSTOPCOND] = 1;
+ if (fplog)
+ {
+ fprintf(fplog, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n", gmx_step_str(step, sbuf), max_hours*0.99);
+ }
+ fprintf(stderr, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n", gmx_step_str(step, sbuf), max_hours*0.99);
+ }
+
+ if (bResetCountersHalfMaxH && MASTER(cr) &&
+ elapsed_time > max_hours*60.0*60.0*0.495)
+ {
+ gs.sig[eglsRESETCOUNTERS] = 1;
+ }
+
+ if (ir->nstlist == -1 && !bRerunMD)
+ {
+ /* When bGStatEveryStep=FALSE, global_stat is only called
+ * when we check the atom displacements, not at NS steps.
+ * This means that also the bonded interaction count check is not
+ * performed immediately after NS. Therefore a few MD steps could
+ * be performed with missing interactions.
+ * But wrong energies are never written to file,
+ * since energies are only written after global_stat
+ * has been called.
+ */
+ if (step >= nlh.step_nscheck)
+ {
+ nlh.nabnsb = natoms_beyond_ns_buffer(ir, fr, &top->cgs,
+ nlh.scale_tot, state->x);
+ }
+ else
+ {
+ /* This is not necessarily true,
+ * but step_nscheck is determined quite conservatively.
+ */
+ nlh.nabnsb = 0;
+ }
+ }
+
+ /* In parallel we only have to check for checkpointing in steps
+ * where we do global communication,
+ * otherwise the other nodes don't know.
+ */
+ if (MASTER(cr) && ((bGStat || !PAR(cr)) &&
+ cpt_period >= 0 &&
+ (cpt_period == 0 ||
+ elapsed_time >= nchkpt*cpt_period*60.0)) &&
+ gs.set[eglsCHKPT] == 0)
+ {
+ gs.sig[eglsCHKPT] = 1;
+ }
+
+ /* at the start of step, randomize or scale the velocities (trotter done elsewhere) */
+ if (EI_VV(ir->eI))
+ {
+ if (!bInitStep)
+ {
+ update_tcouple(step, ir, state, ekind, &MassQ, mdatoms);
+ }
+ if (ETC_ANDERSEN(ir->etc)) /* keep this outside of update_tcouple because of the extra info required to pass */
+ {
+ gmx_bool bIfRandomize;
+ bIfRandomize = update_randomize_velocities(ir, step, cr, mdatoms, state, upd, constr);
+ /* if we have constraints, we have to remove the kinetic energy parallel to the bonds */
+ if (constr && bIfRandomize)
+ {
+ update_constraints(fplog, step, NULL, ir, ekind, mdatoms,
+ state, fr->bMolPBC, graph, f,
+ &top->idef, tmp_vir,
+ cr, nrnb, wcycle, upd, constr,
+ TRUE, bCalcVir, vetanew);
+ }
+ }
+ }
+
+ if (bIterativeCase && do_per_step(step, ir->nstpcouple))
+ {
+ gmx_iterate_init(&iterate, TRUE);
+ /* for iterations, we save these vectors, as we will be redoing the calculations */
+ copy_coupling_state(state, bufstate, ekind, ekind_save, &(ir->opts));
+ }
+
+ bFirstIterate = TRUE;
+ while (bFirstIterate || iterate.bIterationActive)
+ {
+ /* We now restore these vectors to redo the calculation with improved extended variables */
+ if (iterate.bIterationActive)
+ {
+ copy_coupling_state(bufstate, state, ekind_save, ekind, &(ir->opts));
+ }
+
+ /* We make the decision to break or not -after- the calculation of Ekin and Pressure,
+ so scroll down for that logic */
+
+ /* ######### START SECOND UPDATE STEP ################# */
+ /* Box is changed in update() when we do pressure coupling,
+ * but we should still use the old box for energy corrections and when
+ * writing it to the energy file, so it matches the trajectory files for
+ * the same timestep above. Make a copy in a separate array.
+ */
+ copy_mat(state->box, lastbox);
+
+ bOK = TRUE;
+ dvdl_constr = 0;
+
+ if (!(bRerunMD && !rerun_fr.bV && !bForceUpdate))
+ {
+ wallcycle_start(wcycle, ewcUPDATE);
+ /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */
+ if (bTrotter)
+ {
+ if (iterate.bIterationActive)
+ {
+ if (bFirstIterate)
+ {
+ scalevir = 1;
+ }
+ else
+ {
+ /* we use a new value of scalevir to converge the iterations faster */
+ scalevir = tracevir/trace(shake_vir);
+ }
+ msmul(shake_vir, scalevir, shake_vir);
+ m_add(force_vir, shake_vir, total_vir);
+ clear_mat(shake_vir);
+ }
+ trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ3);
+ /* We can only do Berendsen coupling after we have summed
+ * the kinetic energy or virial. Since the happens
+ * in global_state after update, we should only do it at
+ * step % nstlist = 1 with bGStatEveryStep=FALSE.
+ */
+ }
+ else
+ {
+ update_tcouple(step, ir, state, ekind, &MassQ, mdatoms);
+ update_pcouple(fplog, step, ir, state, pcoupl_mu, M, bInitStep);
+ }
+
+ if (bVV)
+ {
+ bUpdateDoLR = (fr->bTwinRange && do_per_step(step, ir->nstcalclr));
+
+ /* velocity half-step update */
+ update_coords(fplog, step, ir, mdatoms, state, fr->bMolPBC, f,
- bUpdateDoLR, fr->f_twin, fcd,
++ bUpdateDoLR, fr->f_twin, bCalcVir ? &fr->vir_twin_constr : NULL, fcd,
+ ekind, M, upd, FALSE, etrtVELOCITY2,
+ cr, nrnb, constr, &top->idef);
+ }
+
+ /* Above, initialize just copies ekinh into ekin,
+ * it doesn't copy position (for VV),
+ * and entire integrator for MD.
+ */
+
+ if (ir->eI == eiVVAK)
+ {
+ copy_rvecn(state->x, cbuf, 0, state->natoms);
+ }
+ bUpdateDoLR = (fr->bTwinRange && do_per_step(step, ir->nstcalclr));
+
+ update_coords(fplog, step, ir, mdatoms, state, fr->bMolPBC, f,
- bUpdateDoLR, fr->f_twin, fcd,
++ bUpdateDoLR, fr->f_twin, bCalcVir ? &fr->vir_twin_constr : NULL, fcd,
+ ekind, M, upd, bInitStep, etrtPOSITION, cr, nrnb, constr, &top->idef);
+ wallcycle_stop(wcycle, ewcUPDATE);
+
+ update_constraints(fplog, step, &dvdl_constr, ir, ekind, mdatoms, state,
+ fr->bMolPBC, graph, f,
+ &top->idef, shake_vir,
+ cr, nrnb, wcycle, upd, constr,
+ FALSE, bCalcVir, state->veta);
+
++ if (bCalcVir && bUpdateDoLR && ir->nstcalclr > 1)
++ {
++ /* Correct the virial for multiple time stepping */
++ m_sub(shake_vir, fr->vir_twin_constr, shake_vir);
++ }
++
+ if (ir->eI == eiVVAK)
+ {
+ /* erase F_EKIN and F_TEMP here? */
+ /* just compute the kinetic energy at the half step to perform a trotter step */
+ compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm,
+ wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
+ constr, NULL, FALSE, lastbox,
+ top_global, &bSumEkinhOld,
+ cglo_flags | CGLO_TEMPERATURE
+ );
+ wallcycle_start(wcycle, ewcUPDATE);
+ trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ4);
+ /* now we know the scaling, we can compute the positions again again */
+ copy_rvecn(cbuf, state->x, 0, state->natoms);
+
+ bUpdateDoLR = (fr->bTwinRange && do_per_step(step, ir->nstcalclr));
+
+ update_coords(fplog, step, ir, mdatoms, state, fr->bMolPBC, f,
++ bUpdateDoLR, fr->f_twin, bCalcVir ? &fr->vir_twin_constr : NULL, fcd,
+ ekind, M, upd, bInitStep, etrtPOSITION, cr, nrnb, constr, &top->idef);
+ wallcycle_stop(wcycle, ewcUPDATE);
+
+ /* do we need an extra constraint here? just need to copy out of state->v to upd->xp? */
+ /* are the small terms in the shake_vir here due
+ * to numerical errors, or are they important
+ * physically? I'm thinking they are just errors, but not completely sure.
+ * For now, will call without actually constraining, constr=NULL*/
+ update_constraints(fplog, step, NULL, ir, ekind, mdatoms,
+ state, fr->bMolPBC, graph, f,
+ &top->idef, tmp_vir,
+ cr, nrnb, wcycle, upd, NULL,
+ FALSE, bCalcVir,
+ state->veta);
+ }
+ if (!bOK)
+ {
+ gmx_fatal(FARGS, "Constraint error: Shake, Lincs or Settle could not solve the constrains");
+ }
+
+ if (fr->bSepDVDL && fplog && do_log)
+ {
+ gmx_print_sepdvdl(fplog, "Constraint dV/dl", 0.0, dvdl_constr);
+ }
+ if (bVV)
+ {
+ /* this factor or 2 correction is necessary
+ because half of the constraint force is removed
+ in the vv step, so we have to double it. See
+ the Redmine issue #1255. It is not yet clear
+ if the factor of 2 is exact, or just a very
+ good approximation, and this will be
+ investigated. The next step is to see if this
+ can be done adding a dhdl contribution from the
+ rattle step, but this is somewhat more
+ complicated with the current code. Will be
+ investigated, hopefully for 4.6.3. However,
+ this current solution is much better than
+ having it completely wrong.
+ */
+ enerd->term[F_DVDL_CONSTR] += 2*dvdl_constr;
+ }
+ else
+ {
+ enerd->term[F_DVDL_CONSTR] += dvdl_constr;
+ }
+ }
+ else if (graph)
+ {
+ /* Need to unshift here */
+ unshift_self(graph, state->box, state->x);
+ }
+
+ if (vsite != NULL)
+ {
+ wallcycle_start(wcycle, ewcVSITECONSTR);
+ if (graph != NULL)
+ {
+ shift_self(graph, state->box, state->x);
+ }
+ construct_vsites(vsite, state->x, ir->delta_t, state->v,
+ top->idef.iparams, top->idef.il,
+ fr->ePBC, fr->bMolPBC, cr, state->box);
+
+ if (graph != NULL)
+ {
+ unshift_self(graph, state->box, state->x);
+ }
+ wallcycle_stop(wcycle, ewcVSITECONSTR);
+ }
+
+ /* ############## IF NOT VV, Calculate globals HERE, also iterate constraints ############ */
+ /* With Leap-Frog we can skip compute_globals at
+ * non-communication steps, but we need to calculate
+ * the kinetic energy one step before communication.
+ */
+ if (bGStat || (!EI_VV(ir->eI) && do_per_step(step+1, nstglobalcomm)))
+ {
+ if (ir->nstlist == -1 && bFirstIterate)
+ {
+ gs.sig[eglsNABNSB] = nlh.nabnsb;
+ }
+ compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm,
+ wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
+ constr,
+ bFirstIterate ? &gs : NULL,
+ (step_rel % gs.nstms == 0) &&
+ (multisim_nsteps < 0 || (step_rel < multisim_nsteps)),
+ lastbox,
+ top_global, &bSumEkinhOld,
+ cglo_flags
+ | (!EI_VV(ir->eI) || bRerunMD ? CGLO_ENERGY : 0)
+ | (!EI_VV(ir->eI) && bStopCM ? CGLO_STOPCM : 0)
+ | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
+ | (!EI_VV(ir->eI) || bRerunMD ? CGLO_PRESSURE : 0)
+ | (iterate.bIterationActive ? CGLO_ITERATE : 0)
+ | (bFirstIterate ? CGLO_FIRSTITERATE : 0)
+ | CGLO_CONSTRAINT
+ );
+ if (ir->nstlist == -1 && bFirstIterate)
+ {
+ nlh.nabnsb = gs.set[eglsNABNSB];
+ gs.set[eglsNABNSB] = 0;
+ }
+ }
+ /* bIterate is set to keep it from eliminating the old ekin kinetic energy terms */
+ /* ############# END CALC EKIN AND PRESSURE ################# */
+
+ /* Note: this is OK, but there are some numerical precision issues with using the convergence of
+ the virial that should probably be addressed eventually. state->veta has better properies,
+ but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
+ generate the new shake_vir, but test the veta value for convergence. This will take some thought. */
+
+ if (iterate.bIterationActive &&
+ done_iterating(cr, fplog, step, &iterate, bFirstIterate,
+ trace(shake_vir), &tracevir))
+ {
+ break;
+ }
+ bFirstIterate = FALSE;
+ }
+
+ if (!bVV || bRerunMD)
+ {
+ /* sum up the foreign energy and dhdl terms for md and sd. currently done every step so that dhdl is correct in the .edr */
+ sum_dhdl(enerd, state->lambda, ir->fepvals);
+ }
+ update_box(fplog, step, ir, mdatoms, state, f,
+ ir->nstlist == -1 ? &nlh.scale_tot : NULL, pcoupl_mu, nrnb, upd);
+
+ /* ################# END UPDATE STEP 2 ################# */
+ /* #### We now have r(t+dt) and v(t+dt/2) ############# */
+
+ /* The coordinates (x) were unshifted in update */
+ if (!bGStat)
+ {
+ /* We will not sum ekinh_old,
+ * so signal that we still have to do it.
+ */
+ bSumEkinhOld = TRUE;
+ }
+
+ /* ######### BEGIN PREPARING EDR OUTPUT ########### */
+
+ /* use the directly determined last velocity, not actually the averaged half steps */
+ if (bTrotter && ir->eI == eiVV)
+ {
+ enerd->term[F_EKIN] = last_ekin;
+ }
+ enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
+
+ if (bVV)
+ {
+ enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + saved_conserved_quantity;
+ }
+ else
+ {
+ enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + compute_conserved_from_auxiliary(ir, state, &MassQ);
+ }
+ /* ######### END PREPARING EDR OUTPUT ########### */
+
+ /* Output stuff */
+ if (MASTER(cr))
+ {
+ gmx_bool do_dr, do_or;
+
+ if (fplog && do_log && bDoExpanded)
+ {
+ /* only needed if doing expanded ensemble */
+ PrintFreeEnergyInfoToFile(fplog, ir->fepvals, ir->expandedvals, ir->bSimTemp ? ir->simtempvals : NULL,
+ &state_global->dfhist, state->fep_state, ir->nstlog, step);
+ }
+ if (!(bStartingFromCpt && (EI_VV(ir->eI))))
+ {
+ if (bCalcEner)
+ {
+ upd_mdebin(mdebin, bDoDHDL, TRUE,
+ t, mdatoms->tmass, enerd, state,
+ ir->fepvals, ir->expandedvals, lastbox,
+ shake_vir, force_vir, total_vir, pres,
+ ekind, mu_tot, constr);
+ }
+ else
+ {
+ upd_mdebin_step(mdebin);
+ }
+
+ do_dr = do_per_step(step, ir->nstdisreout);
+ do_or = do_per_step(step, ir->nstorireout);
+
+ print_ebin(mdoutf_get_fp_ene(outf), do_ene, do_dr, do_or, do_log ? fplog : NULL,
+ step, t,
+ eprNORMAL, bCompact, mdebin, fcd, groups, &(ir->opts));
+ }
+ if (ir->ePull != epullNO)
+ {
+ pull_print_output(ir->pull, step, t);
+ }
+
+ if (do_per_step(step, ir->nstlog))
+ {
+ if (fflush(fplog) != 0)
+ {
+ gmx_fatal(FARGS, "Cannot flush logfile - maybe you are out of disk space?");
+ }
+ }
+ }
+ if (bDoExpanded)
+ {
+ /* Have to do this part _after_ outputting the logfile and the edr file */
+ /* Gets written into the state at the beginning of next loop*/
+ state->fep_state = lamnew;
+ }
+ /* Print the remaining wall clock time for the run */
+ if (MULTIMASTER(cr) && (do_verbose || gmx_got_usr_signal()) && !bPMETuneRunning)
+ {
+ if (shellfc)
+ {
+ fprintf(stderr, "\n");
+ }
+ print_time(stderr, walltime_accounting, step, ir, cr);
+ }
+
+ /* Ion/water position swapping.
+ * Not done in last step since trajectory writing happens before this call
+ * in the MD loop and exchanges would be lost anyway. */
+ bNeedRepartition = FALSE;
+ if ((ir->eSwapCoords != eswapNO) && (step > 0) && !bLastStep &&
+ do_per_step(step, ir->swap->nstswap))
+ {
+ bNeedRepartition = do_swapcoords(cr, step, t, ir, wcycle,
+ bRerunMD ? rerun_fr.x : state->x,
+ bRerunMD ? rerun_fr.box : state->box,
+ top_global, MASTER(cr) && bVerbose, bRerunMD);
+
+ if (bNeedRepartition && DOMAINDECOMP(cr))
+ {
+ dd_collect_state(cr->dd, state, state_global);
+ }
+ }
+
+ /* Replica exchange */
+ bExchanged = FALSE;
+ if ((repl_ex_nst > 0) && (step > 0) && !bLastStep &&
+ do_per_step(step, repl_ex_nst))
+ {
+ bExchanged = replica_exchange(fplog, cr, repl_ex,
+ state_global, enerd,
+ state, step, t);
+ }
+
+ if ( (bExchanged || bNeedRepartition) && DOMAINDECOMP(cr) )
+ {
+ dd_partition_system(fplog, step, cr, TRUE, 1,
+ state_global, top_global, ir,
+ state, &f, mdatoms, top, fr,
+ vsite, shellfc, constr,
+ nrnb, wcycle, FALSE);
+ }
+
+ bFirstStep = FALSE;
+ bInitStep = FALSE;
+ bStartingFromCpt = FALSE;
+
+ /* ####### SET VARIABLES FOR NEXT ITERATION IF THEY STILL NEED IT ###### */
+ /* With all integrators, except VV, we need to retain the pressure
+ * at the current step for coupling at the next step.
+ */
+ if ((state->flags & (1<<estPRES_PREV)) &&
+ (bGStatEveryStep ||
+ (ir->nstpcouple > 0 && step % ir->nstpcouple == 0)))
+ {
+ /* Store the pressure in t_state for pressure coupling
+ * at the next MD step.
+ */
+ copy_mat(pres, state->pres_prev);
+ }
+
+ /* ####### END SET VARIABLES FOR NEXT ITERATION ###### */
+
+ if ( (membed != NULL) && (!bLastStep) )
+ {
+ rescale_membed(step_rel, membed, state_global->x);
+ }
+
+ if (bRerunMD)
+ {
+ if (MASTER(cr))
+ {
+ /* read next frame from input trajectory */
+ bNotLastFrame = read_next_frame(oenv, status, &rerun_fr);
+ }
+
+ if (PAR(cr))
+ {
+ rerun_parallel_comm(cr, &rerun_fr, &bNotLastFrame);
+ }
+ }
+
+ if (!bRerunMD || !rerun_fr.bStep)
+ {
+ /* increase the MD step number */
+ step++;
+ step_rel++;
+ }
+
+ cycles = wallcycle_stop(wcycle, ewcSTEP);
+ if (DOMAINDECOMP(cr) && wcycle)
+ {
+ dd_cycles_add(cr->dd, cycles, ddCyclStep);
+ }
+
+ if (bPMETuneRunning || bPMETuneTry)
+ {
+ /* PME grid + cut-off optimization with GPUs or PME nodes */
+
+ /* Count the total cycles over the last steps */
+ cycles_pmes += cycles;
+
+ /* We can only switch cut-off at NS steps */
+ if (step % ir->nstlist == 0)
+ {
+ /* PME grid + cut-off optimization with GPUs or PME nodes */
+ if (bPMETuneTry)
+ {
+ if (DDMASTER(cr->dd))
+ {
+ /* PME node load is too high, start tuning */
+ bPMETuneRunning = (dd_pme_f_ratio(cr->dd) >= 1.05);
+ }
+ dd_bcast(cr->dd, sizeof(gmx_bool), &bPMETuneRunning);
+
+ if (bPMETuneRunning || step_rel > ir->nstlist*50)
+ {
+ bPMETuneTry = FALSE;
+ }
+ }
+ if (bPMETuneRunning)
+ {
+ /* init_step might not be a multiple of nstlist,
+ * but the first cycle is always skipped anyhow.
+ */
+ bPMETuneRunning =
+ pme_load_balance(pme_loadbal, cr,
+ (bVerbose && MASTER(cr)) ? stderr : NULL,
+ fplog,
+ ir, state, cycles_pmes,
+ fr->ic, fr->nbv, &fr->pmedata,
+ step);
+
+ /* Update constants in forcerec/inputrec to keep them in sync with fr->ic */
+ fr->ewaldcoeff_q = fr->ic->ewaldcoeff_q;
+ fr->ewaldcoeff_lj = fr->ic->ewaldcoeff_lj;
+ fr->rlist = fr->ic->rlist;
+ fr->rlistlong = fr->ic->rlistlong;
+ fr->rcoulomb = fr->ic->rcoulomb;
+ fr->rvdw = fr->ic->rvdw;
+
+ if (ir->eDispCorr != edispcNO)
+ {
+ calc_enervirdiff(NULL, ir->eDispCorr, fr);
+ }
+ }
+ cycles_pmes = 0;
+ }
+ }
+
+ if (step_rel == wcycle_get_reset_counters(wcycle) ||
+ gs.set[eglsRESETCOUNTERS] != 0)
+ {
+ /* Reset all the counters related to performance over the run */
+ reset_all_counters(fplog, cr, step, &step_rel, ir, wcycle, nrnb, walltime_accounting,
+ fr->nbv != NULL && fr->nbv->bUseGPU ? fr->nbv->cu_nbv : NULL);
+ wcycle_set_reset_counters(wcycle, -1);
+ if (!(cr->duty & DUTY_PME))
+ {
+ /* Tell our PME node to reset its counters */
+ gmx_pme_send_resetcounters(cr, step);
+ }
+ /* Correct max_hours for the elapsed time */
+ max_hours -= elapsed_time/(60.0*60.0);
+ bResetCountersHalfMaxH = FALSE;
+ gs.set[eglsRESETCOUNTERS] = 0;
+ }
+
+ /* If bIMD is TRUE, the master updates the IMD energy record and sends positions to VMD client */
+ IMD_prep_energies_send_positions(ir->bIMD && MASTER(cr), bIMDstep, ir->imd, enerd, step, bCalcEner, wcycle);
+
+ }
+ /* End of main MD loop */
+ debug_gmx();
+
+ /* Stop measuring walltime */
+ walltime_accounting_end(walltime_accounting);
+
+ if (bRerunMD && MASTER(cr))
+ {
+ close_trj(status);
+ }
+
+ if (!(cr->duty & DUTY_PME))
+ {
+ /* Tell the PME only node to finish */
+ gmx_pme_send_finish(cr);
+ }
+
+ if (MASTER(cr))
+ {
+ if (ir->nstcalcenergy > 0 && !bRerunMD)
+ {
+ print_ebin(mdoutf_get_fp_ene(outf), FALSE, FALSE, FALSE, fplog, step, t,
+ eprAVER, FALSE, mdebin, fcd, groups, &(ir->opts));
+ }
+ }
+
+ done_mdoutf(outf);
+ debug_gmx();
+
+ if (ir->nstlist == -1 && nlh.nns > 0 && fplog)
+ {
+ fprintf(fplog, "Average neighborlist lifetime: %.1f steps, std.dev.: %.1f steps\n", nlh.s1/nlh.nns, sqrt(nlh.s2/nlh.nns - sqr(nlh.s1/nlh.nns)));
+ fprintf(fplog, "Average number of atoms that crossed the half buffer length: %.1f\n\n", nlh.ab/nlh.nns);
+ }
+
+ if (pme_loadbal != NULL)
+ {
+ pme_loadbal_done(pme_loadbal, cr, fplog,
+ fr->nbv != NULL && fr->nbv->bUseGPU);
+ }
+
+ if (shellfc && fplog)
+ {
+ fprintf(fplog, "Fraction of iterations that converged: %.2f %%\n",
+ (nconverged*100.0)/step_rel);
+ fprintf(fplog, "Average number of force evaluations per MD step: %.2f\n\n",
+ tcount/step_rel);
+ }
+
+ if (repl_ex_nst > 0 && MASTER(cr))
+ {
+ print_replica_exchange_statistics(fplog, repl_ex);
+ }
+
+ /* IMD cleanup, if bIMD is TRUE. */
+ IMD_finalize(ir->bIMD, ir->imd);
+
+ walltime_accounting_set_nsteps_done(walltime_accounting, step_rel);
+
+ return 0;
+}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff