From: Mark Abraham Date: Thu, 5 Jun 2014 21:34:10 +0000 (+0200) Subject: Merge release-4-6 into release-5-0 X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=15184a1bb5e9ac037b2456628484cd1f26a58ecd;p=alexxy%2Fgromacs.git Merge release-4-6 into release-5-0 Conflicts: src/gromacs/gmxana/gmx_energy.c Eliminated time parameter (removed in release-5-0), used gmx_ffclose (fixed in release-5-0). src/gromacs/mdlib/update.c Eliminated fplog parameter (removed in release-5-0), otherwise added new MTS-fixing functionality from release-4-6. src/programs/mdrun/md.c Added new code path for MTS (from release-4-6), eliminated wcycle and bInitStep parameters (from release-5-0). Change-Id: I303bb649ed2373656aa4420d0afb27f738a9ea30 --- 15184a1bb5e9ac037b2456628484cd1f26a58ecd diff --cc src/gromacs/gmxana/gmx_energy.c index 38555952e0,0000000000..8f199dbe7e mode 100644,000000..100644 --- a/src/gromacs/gmxana/gmx_energy.c +++ b/src/gromacs/gmxana/gmx_energy.c @@@ -1,2832 -1,0 +1,2845 @@@ +/* + * 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 +#endif + +#include +#include + +#include "typedefs.h" +#include "gmx_fatal.h" +#include "vec.h" +#include "gromacs/utility/cstringutil.h" +#include "gromacs/utility/smalloc.h" +#include "gromacs/fileio/enxio.h" +#include "gromacs/commandline/pargs.h" +#include "names.h" +#include "copyrite.h" +#include "macros.h" +#include "xvgr.h" +#include "gstat.h" +#include "physics.h" +#include "gromacs/fileio/tpxio.h" +#include "gromacs/fileio/trxio.h" +#include "viewit.h" +#include "mtop_util.h" +#include "gmx_ana.h" +#include "mdebin.h" + +static real minthird = -1.0/3.0, minsixth = -1.0/6.0; + +typedef struct { + real sum; + real sum2; +} exactsum_t; + +typedef struct { + real *ener; + exactsum_t *es; + gmx_bool bExactStat; + double av; + double rmsd; + double ee; + double slope; +} enerdat_t; + +typedef struct { + gmx_int64_t nsteps; + gmx_int64_t npoints; + int nframes; + int *step; + int *steps; + int *points; + enerdat_t *s; +} 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, - int nframes, real **sum, - real V, real T, int nsteps, double time[], ++ int nint, real **eneint, ++ real V, real T, double dt, + const output_env_t oenv) +{ + FILE *fp0, *fp1; + double av[4], avold[4]; - double fac, dt, di; ++ double fac, di; + int i, j, m, nf4; + - if (nframes < 1) - { - return; - } - - dt = (time[1]-time[0]); - nf4 = nframes/4+1; ++ nf4 = nint/4 + 1; + + for (i = 0; i <= nsets; i++) + { + avold[i] = 0; + } + fp0 = xvgropen(fni, "Shear viscosity integral", + "Time (ps)", "(kg m\\S-1\\N s\\S-1\\N ps)", oenv); + fp1 = xvgropen(fn, "Shear viscosity using Einstein relation", + "Time (ps)", "(kg m\\S-1\\N s\\S-1\\N)", oenv); - for (i = 1; i < nf4; i++) ++ for (i = 0; i < nf4; i++) + { - fac = dt*nframes/nsteps; + for (m = 0; m <= nsets; m++) + { + av[m] = 0; + } - for (j = 0; j < nframes-i; j++) ++ for (j = 0; j < nint - i; j++) + { + for (m = 0; m < nsets; m++) + { - di = sqr(fac*(sum[m][j+i]-sum[m][j])); ++ di = sqr(eneint[m][j+i] - eneint[m][j]); + + av[m] += di; + av[nsets] += di/nsets; + } + } + /* Convert to SI for the viscosity */ - fac = (V*NANO*NANO*NANO*PICO*1e10)/(2*BOLTZMANN*T)/(nframes-i); - fprintf(fp0, "%10g", time[i]-time[0]); ++ fac = (V*NANO*NANO*NANO*PICO*1e10)/(2*BOLTZMANN*T)/(nint - i); ++ fprintf(fp0, "%10g", i*dt); + for (m = 0; (m <= nsets); m++) + { + av[m] = fac*av[m]; + fprintf(fp0, " %10g", av[m]); + } + fprintf(fp0, "\n"); - fprintf(fp1, "%10g", 0.5*(time[i]+time[i-1])-time[0]); ++ fprintf(fp1, "%10g", (i + 0.5)*dt); + for (m = 0; (m <= nsets); m++) + { + fprintf(fp1, " %10g", (av[m]-avold[m])/dt); + 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, - double time[], real reftemp, ++ 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"); + } + else + { + fprintf(stdout, "\n"); + } + fprintf(stdout, "-------------------------------------------------------------------------------\n"); + + /* Initiate locals, only used with -sum */ + expEtot = 0; + if (bFee) + { + beta = 1.0/(BOLTZ*reftemp); + snew(fee, nset); + } + for (i = 0; (i < nset); i++) + { + aver = edat->s[i].av; + stddev = edat->s[i].rmsd; + errest = edat->s[i].ee; + + if (bFee) + { + expE = 0; + for (j = 0; (j < edat->nframes); j++) + { + expE += exp(beta*(edat->s[i].ener[j] - aver)/nmol); + } + if (bSum) + { + expEtot += expE/edat->nframes; + } + + fee[i] = log(expE/edat->nframes)/beta + aver/nmol; + } + if (strstr(leg[i], "empera") != NULL) + { + Temp = aver; + } + else if (strstr(leg[i], "olum") != NULL) + { + Vaver = aver; + } + else if (strstr(leg[i], "essure") != NULL) + { + Pres = aver; + } + if (bIsEner[i]) + { + pr_aver = aver/nmol-ezero; + pr_stddev = stddev/nmol; + pr_errest = errest/nmol; + } + else + { + pr_aver = aver; + pr_stddev = stddev; + pr_errest = errest; + } + + /* Multiply the slope in steps with the number of steps taken */ + totaldrift = (edat->nsteps - 1)*edat->s[i].slope; + if (bIsEner[i]) + { + totaldrift /= nmol; + } + + fprintf(stdout, "%-24s %10g %10s %10g %10g", + leg[i], pr_aver, ee_pr(pr_errest, eebuf), pr_stddev, totaldrift); + if (bFee) + { + fprintf(stdout, " %10g", fee[i]); + } + + fprintf(stdout, " (%s)\n", enm[set[i]].unit); + + if (bFluct) + { + for (j = 0; (j < edat->nframes); j++) + { + edat->s[i].ener[j] -= aver; + } + } + } + if (bSum) + { + totaldrift = (edat->nsteps - 1)*esum->s[0].slope; + fprintf(stdout, "%-24s %10g %10s %10s %10g (%s)", + "Total", esum->s[0].av/nmol, ee_pr(esum->s[0].ee/nmol, eebuf), + "--", totaldrift/nmol, enm[set[0]].unit); + /* pr_aver,pr_stddev,a,totaldrift */ + if (bFee) + { + fprintf(stdout, " %10g %10g\n", + log(expEtot)/beta + esum->s[0].av/nmol, log(expEtot)/beta); + } + else + { + fprintf(stdout, "\n"); + } + } + + /* Do correlation function */ + if (edat->nframes > 1) + { + Dt = delta_t/(edat->nframes - 1); + } + else + { + Dt = 0; + } + if (bVisco) + { + const char* leg[] = { "Shear", "Bulk" }; + real factor; + real **eneset; - real **enesum; ++ 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); + } - snew(enesum, 3); - for (i = 0; i < 3; i++) - { - snew(enesum[i], edat->nframes); - } + for (i = 0; (i < edat->nframes); i++) + { + eneset[0][i] = 0.5*(edat->s[1].ener[i]+edat->s[3].ener[i]); + eneset[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; - enesum[0][i] = 0.5*(edat->s[1].es[i].sum+edat->s[3].es[i].sum); - enesum[1][i] = 0.5*(edat->s[2].es[i].sum+edat->s[6].es[i].sum); - enesum[2][i] = 0.5*(edat->s[5].es[i].sum+edat->s[7].es[i].sum); ++ } ++ ++ /* Determine integrals of the off-diagonal pressure elements */ ++ snew(eneint, 3); ++ for (i = 0; i < 3; i++) ++ { ++ snew(eneint[i], edat->nframes + 1); ++ } ++ eneint[0][0] = 0; ++ eneint[1][0] = 0; ++ eneint[2][0] = 0; ++ for (i = 0; i < edat->nframes; i++) ++ { ++ eneint[0][i+1] = eneint[0][i] + 0.5*(edat->s[1].es[i].sum + edat->s[3].es[i].sum)*Dt/edat->points[i]; ++ eneint[1][i+1] = eneint[1][i] + 0.5*(edat->s[2].es[i].sum + edat->s[6].es[i].sum)*Dt/edat->points[i]; ++ eneint[2][i+1] = eneint[2][i] + 0.5*(edat->s[5].es[i].sum + edat->s[7].es[i].sum)*Dt/edat->points[i]; + } + + einstein_visco("evisco.xvg", "eviscoi.xvg", - 3, edat->nframes, enesum, Vaver, Temp, nsteps, time, oenv); ++ 3, edat->nframes+1, eneint, Vaver, Temp, Dt, oenv); ++ ++ for (i = 0; i < 3; i++) ++ { ++ sfree(eneint[i]); ++ } ++ sfree(eneint); + + /*do_autocorr(corrfn,buf,nenergy,3,eneset,Dt,eacNormal,TRUE);*/ + /* Do it for shear viscosity */ + strcpy(buf, "Shear Viscosity"); + low_do_autocorr(corrfn, oenv, buf, edat->nframes, 3, + (edat->nframes+1)/2, eneset, Dt, + eacNormal, 1, TRUE, FALSE, FALSE, 0.0, 0.0, 0); + + /* Now for bulk viscosity */ + strcpy(buf, "Bulk Viscosity"); + low_do_autocorr(corrfn, oenv, buf, edat->nframes, 1, + (edat->nframes+1)/2, &(eneset[11]), Dt, + eacNormal, 1, TRUE, FALSE, FALSE, 0.0, 0.0, 0); + + factor = (Vaver*1e-26/(BOLTZMANN*Temp))*Dt; + fp = xvgropen(visfn, buf, "Time (ps)", "\\8h\\4 (cp)", oenv); + xvgr_legend(fp, asize(leg), leg, oenv); + + /* Use trapezium rule for integration */ + integral = 0; + intBulk = 0; + nout = get_acfnout(); + if ((nout < 2) || (nout >= edat->nframes/2)) + { + nout = edat->nframes/2; + } + for (i = 1; (i < nout); i++) + { + integral += 0.5*(eneset[0][i-1] + eneset[0][i])*factor; + intBulk += 0.5*(eneset[11][i-1] + eneset[11][i])*factor; + fprintf(fp, "%10g %10g %10g\n", (i*Dt), integral, intBulk); + } + gmx_ffclose(fp); ++ ++ for (i = 0; i < 12; i++) ++ { ++ sfree(eneset[i]); ++ } ++ sfree(eneset); + } + else if (bCorr) + { + if (bFluct) + { + strcpy(buf, "Autocorrelation of Energy Fluctuations"); + } + else + { + strcpy(buf, "Energy Autocorrelation"); + } +#if 0 + do_autocorr(corrfn, oenv, buf, edat->nframes, + bSum ? 1 : nset, + bSum ? &edat->s[nset-1].ener : eneset, + (delta_t/edat->nframes), eacNormal, FALSE); +#endif + } + } +} + +static void print_time(FILE *fp, double t) +{ + fprintf(fp, "%12.6f", t); +} + +static void print1(FILE *fp, gmx_bool bDp, real e) +{ + if (bDp) + { + fprintf(fp, " %16.12f", e); + } + else + { + fprintf(fp, " %10.6f", e); + } +} + +static void fec(const char *ene2fn, const char *runavgfn, + real reftemp, int nset, int set[], char *leg[], + enerdata_t *edat, double time[], + const output_env_t oenv) +{ + const char * ravgleg[] = { + "\\8D\\4E = E\\sB\\N-E\\sA\\N", + "\\s0..t\\N" + }; + FILE *fp; + ener_file_t enx; + int nre, timecheck, step, nenergy, nenergy2, maxenergy; + int i, j; + gmx_bool bCont; + real aver, beta; + real **eneset2; + double dE, sum; + gmx_enxnm_t *enm = NULL; + t_enxframe *fr; + char buf[22]; + + /* read second energy file */ + snew(fr, 1); + enm = NULL; + enx = open_enx(ene2fn, "r"); + do_enxnms(enx, &(fr->nre), &enm); + + snew(eneset2, nset+1); + nenergy2 = 0; + maxenergy = 0; + timecheck = 0; + do + { + /* This loop searches for the first frame (when -b option is given), + * or when this has been found it reads just one energy frame + */ + do + { + bCont = do_enx(enx, fr); + + if (bCont) + { + timecheck = check_times(fr->t); + } + + } + while (bCont && (timecheck < 0)); + + /* Store energies for analysis afterwards... */ + if ((timecheck == 0) && bCont) + { + if (fr->nre > 0) + { + if (nenergy2 >= maxenergy) + { + maxenergy += 1000; + for (i = 0; i <= nset; i++) + { + srenew(eneset2[i], maxenergy); + } + } + if (fr->t != time[nenergy2]) + { + fprintf(stderr, "\nWARNING time mismatch %g!=%g at frame %s\n", + fr->t, time[nenergy2], gmx_step_str(fr->step, buf)); + } + for (i = 0; i < nset; i++) + { + eneset2[i][nenergy2] = fr->ener[set[i]].e; + } + nenergy2++; + } + } + } + while (bCont && (timecheck == 0)); + + /* check */ + if (edat->nframes != nenergy2) + { + fprintf(stderr, "\nWARNING file length mismatch %d!=%d\n", + edat->nframes, nenergy2); + } + nenergy = min(edat->nframes, nenergy2); + + /* calculate fe difference dF = -kT ln < exp(-(E_B-E_A)/kT) >_A */ + fp = NULL; + if (runavgfn) + { + fp = xvgropen(runavgfn, "Running average free energy difference", + "Time (" unit_time ")", "\\8D\\4E (" unit_energy ")", oenv); + xvgr_legend(fp, asize(ravgleg), ravgleg, oenv); + } + fprintf(stdout, "\n%-24s %10s\n", + "Energy", "dF = -kT ln < exp(-(EB-EA)/kT) >A"); + sum = 0; + beta = 1.0/(BOLTZ*reftemp); + for (i = 0; i < nset; i++) + { + if (gmx_strcasecmp(leg[i], enm[set[i]].name) != 0) + { + fprintf(stderr, "\nWARNING energy set name mismatch %s!=%s\n", + leg[i], enm[set[i]].name); + } + for (j = 0; j < nenergy; j++) + { + dE = eneset2[i][j] - edat->s[i].ener[j]; + sum += exp(-dE*beta); + if (fp) + { + fprintf(fp, "%10g %10g %10g\n", + time[j], dE, -BOLTZ*reftemp*log(sum/(j+1)) ); + } + } + aver = -BOLTZ*reftemp*log(sum/nenergy); + fprintf(stdout, "%-24s %10g\n", leg[i], aver); + } + if (fp) + { + gmx_ffclose(fp); + } + sfree(fr); +} + + +static void do_dhdl(t_enxframe *fr, t_inputrec *ir, FILE **fp_dhdl, + const char *filename, gmx_bool bDp, + int *blocks, int *hists, int *samples, int *nlambdas, + const output_env_t oenv) +{ + const char *dhdl = "dH/d\\lambda", *deltag = "\\DeltaH", *lambda = "\\lambda"; + char title[STRLEN], label_x[STRLEN], label_y[STRLEN], legend[STRLEN]; + char buf[STRLEN]; + gmx_bool first = FALSE; + int nblock_hist = 0, nblock_dh = 0, nblock_dhcoll = 0; + int i, j, k; + /* coll data */ + double temp = 0, start_time = 0, delta_time = 0, start_lambda = 0, delta_lambda = 0; + static int setnr = 0; + double *native_lambda_vec = NULL; + const char **lambda_components = NULL; + int n_lambda_vec = 0; + gmx_bool changing_lambda = FALSE; + int lambda_fep_state; + + /* now count the blocks & handle the global dh data */ + for (i = 0; i < fr->nblock; i++) + { + if (fr->block[i].id == enxDHHIST) + { + nblock_hist++; + } + else if (fr->block[i].id == enxDH) + { + nblock_dh++; + } + else if (fr->block[i].id == enxDHCOLL) + { + nblock_dhcoll++; + if ( (fr->block[i].nsub < 1) || + (fr->block[i].sub[0].type != xdr_datatype_double) || + (fr->block[i].sub[0].nr < 5)) + { + gmx_fatal(FARGS, "Unexpected block data"); + } + + /* read the data from the DHCOLL block */ + temp = fr->block[i].sub[0].dval[0]; + start_time = fr->block[i].sub[0].dval[1]; + delta_time = fr->block[i].sub[0].dval[2]; + start_lambda = fr->block[i].sub[0].dval[3]; + delta_lambda = fr->block[i].sub[0].dval[4]; + changing_lambda = (delta_lambda != 0); + if (fr->block[i].nsub > 1) + { + lambda_fep_state = fr->block[i].sub[1].ival[0]; + if (n_lambda_vec == 0) + { + n_lambda_vec = fr->block[i].sub[1].ival[1]; + } + else + { + if (n_lambda_vec != fr->block[i].sub[1].ival[1]) + { + gmx_fatal(FARGS, + "Unexpected change of basis set in lambda"); + } + } + if (lambda_components == NULL) + { + snew(lambda_components, n_lambda_vec); + } + if (native_lambda_vec == NULL) + { + snew(native_lambda_vec, n_lambda_vec); + } + for (j = 0; j < n_lambda_vec; j++) + { + native_lambda_vec[j] = fr->block[i].sub[0].dval[5+j]; + lambda_components[j] = + efpt_singular_names[fr->block[i].sub[1].ival[2+j]]; + } + } + } + } + + if (nblock_hist == 0 && nblock_dh == 0) + { + /* don't do anything */ + return; + } + if (nblock_hist > 0 && nblock_dh > 0) + { + gmx_fatal(FARGS, "This energy file contains both histogram dhdl data and non-histogram dhdl data. Don't know what to do."); + } + if (!*fp_dhdl) + { + if (nblock_dh > 0) + { + /* we have standard, non-histogram data -- + call open_dhdl to open the file */ + /* TODO this is an ugly hack that needs to be fixed: this will only + work if the order of data is always the same and if we're + only using the g_energy compiled with the mdrun that produced + the ener.edr. */ + *fp_dhdl = open_dhdl(filename, ir, oenv); + } + else + { + sprintf(title, "N(%s)", deltag); + sprintf(label_x, "%s (%s)", deltag, unit_energy); + sprintf(label_y, "Samples"); + *fp_dhdl = xvgropen_type(filename, title, label_x, label_y, exvggtXNY, oenv); + sprintf(buf, "T = %g (K), %s = %g", temp, lambda, start_lambda); + xvgr_subtitle(*fp_dhdl, buf, oenv); + } + } + + (*hists) += nblock_hist; + (*blocks) += nblock_dh; + (*nlambdas) = nblock_hist+nblock_dh; + + /* write the data */ + if (nblock_hist > 0) + { + gmx_int64_t sum = 0; + /* histograms */ + for (i = 0; i < fr->nblock; i++) + { + t_enxblock *blk = &(fr->block[i]); + if (blk->id == enxDHHIST) + { + double foreign_lambda, dx; + gmx_int64_t x0; + int nhist, derivative; + + /* check the block types etc. */ + if ( (blk->nsub < 2) || + (blk->sub[0].type != xdr_datatype_double) || + (blk->sub[1].type != xdr_datatype_int64) || + (blk->sub[0].nr < 2) || + (blk->sub[1].nr < 2) ) + { + gmx_fatal(FARGS, "Unexpected block data in file"); + } + foreign_lambda = blk->sub[0].dval[0]; + dx = blk->sub[0].dval[1]; + nhist = blk->sub[1].lval[0]; + derivative = blk->sub[1].lval[1]; + for (j = 0; j < nhist; j++) + { + const char *lg[1]; + x0 = blk->sub[1].lval[2+j]; + + if (!derivative) + { + sprintf(legend, "N(%s(%s=%g) | %s=%g)", + deltag, lambda, foreign_lambda, + lambda, start_lambda); + } + else + { + sprintf(legend, "N(%s | %s=%g)", + dhdl, lambda, start_lambda); + } + + lg[0] = legend; + xvgr_new_dataset(*fp_dhdl, setnr, 1, lg, oenv); + setnr++; + for (k = 0; k < blk->sub[j+2].nr; k++) + { + int hist; + double xmin, xmax; + + hist = blk->sub[j+2].ival[k]; + xmin = (x0+k)*dx; + xmax = (x0+k+1)*dx; + fprintf(*fp_dhdl, "%g %d\n%g %d\n", xmin, hist, + xmax, hist); + sum += hist; + } + /* multiple histogram data blocks in one histogram + mean that the second one is the reverse of the first one: + for dhdl derivatives, it's important to know both the + maximum and minimum values */ + dx = -dx; + } + } + } + (*samples) += (int)(sum/nblock_hist); + } + else + { + /* raw dh */ + int len = 0; + char **setnames = NULL; + int nnames = nblock_dh; + + for (i = 0; i < fr->nblock; i++) + { + t_enxblock *blk = &(fr->block[i]); + if (blk->id == enxDH) + { + if (len == 0) + { + len = blk->sub[2].nr; + } + else + { + if (len != blk->sub[2].nr) + { + gmx_fatal(FARGS, "Length inconsistency in dhdl data"); + } + } + } + } + (*samples) += len; + + for (i = 0; i < len; i++) + { + double time = start_time + delta_time*i; + + fprintf(*fp_dhdl, "%.4f ", time); + + for (j = 0; j < fr->nblock; j++) + { + t_enxblock *blk = &(fr->block[j]); + if (blk->id == enxDH) + { + double value; + if (blk->sub[2].type == xdr_datatype_float) + { + value = blk->sub[2].fval[i]; + } + else + { + value = blk->sub[2].dval[i]; + } + /* we need to decide which data type it is based on the count*/ + + if (j == 1 && ir->bExpanded) + { + fprintf(*fp_dhdl, "%4d", (int)value); /* if expanded ensembles and zero, this is a state value, it's an integer. We need a cleaner conditional than if j==1! */ + } + else + { + if (bDp) + { + fprintf(*fp_dhdl, " %#.12g", value); /* print normal precision */ + } + else + { + fprintf(*fp_dhdl, " %#.8g", value); /* print normal precision */ + } + } + } + } + fprintf(*fp_dhdl, "\n"); + } + } +} + + +int gmx_energy(int argc, char *argv[]) +{ + const char *desc[] = { + "[THISMODULE] extracts energy components or distance restraint", + "data from an energy file. The user is prompted to interactively", + "select the desired energy terms.[PAR]", + + "Average, RMSD, and drift are calculated with full precision from the", + "simulation (see printed manual). Drift is calculated by performing", + "a least-squares fit of the data to a straight line. The reported total drift", + "is the difference of the fit at the first and last point.", + "An error estimate of the average is given based on a block averages", + "over 5 blocks using the full-precision averages. The error estimate", + "can be performed over multiple block lengths with the options", + "[TT]-nbmin[tt] and [TT]-nbmax[tt].", + "[BB]Note[bb] that in most cases the energy files contains averages over all", + "MD steps, or over many more points than the number of frames in", + "energy file. This makes the [THISMODULE] statistics output more accurate", + "than the [TT].xvg[tt] output. When exact averages are not present in the energy", + "file, the statistics mentioned above are simply over the single, per-frame", + "energy values.[PAR]", + + "The term fluctuation gives the RMSD around the least-squares fit.[PAR]", + + "Some fluctuation-dependent properties can be calculated provided", + "the correct energy terms are selected, and that the command line option", + "[TT]-fluct_props[tt] is given. The following properties", + "will be computed:[BR]", + "Property Energy terms needed[BR]", + "---------------------------------------------------[BR]", + "Heat capacity C[SUB]p[sub] (NPT sims): Enthalpy, Temp [BR]", + "Heat capacity C[SUB]v[sub] (NVT sims): Etot, Temp [BR]", + "Thermal expansion coeff. (NPT): Enthalpy, Vol, Temp[BR]", + "Isothermal compressibility: Vol, Temp [BR]", + "Adiabatic bulk modulus: Vol, Temp [BR]", + "---------------------------------------------------[BR]", + "You always need to set the number of molecules [TT]-nmol[tt].", + "The C[SUB]p[sub]/C[SUB]v[sub] computations do [BB]not[bb] include any corrections", + "for quantum effects. Use the [gmx-dos] program if you need that (and you do).[PAR]" + "When the [TT]-viol[tt] option is set, the time averaged", + "violations are plotted and the running time-averaged and", + "instantaneous sum of violations are recalculated. Additionally", + "running time-averaged and instantaneous distances between", + "selected pairs can be plotted with the [TT]-pairs[tt] option.[PAR]", + + "Options [TT]-ora[tt], [TT]-ort[tt], [TT]-oda[tt], [TT]-odr[tt] and", + "[TT]-odt[tt] are used for analyzing orientation restraint data.", + "The first two options plot the orientation, the last three the", + "deviations of the orientations from the experimental values.", + "The options that end on an 'a' plot the average over time", + "as a function of restraint. The options that end on a 't'", + "prompt the user for restraint label numbers and plot the data", + "as a function of time. Option [TT]-odr[tt] plots the RMS", + "deviation as a function of restraint.", + "When the run used time or ensemble averaged orientation restraints,", + "option [TT]-orinst[tt] can be used to analyse the instantaneous,", + "not ensemble-averaged orientations and deviations instead of", + "the time and ensemble averages.[PAR]", + + "Option [TT]-oten[tt] plots the eigenvalues of the molecular order", + "tensor for each orientation restraint experiment. With option", + "[TT]-ovec[tt] also the eigenvectors are plotted.[PAR]", + + "Option [TT]-odh[tt] extracts and plots the free energy data", + "(Hamiltoian differences and/or the Hamiltonian derivative dhdl)", + "from the [TT]ener.edr[tt] file.[PAR]", + + "With [TT]-fee[tt] an estimate is calculated for the free-energy", + "difference with an ideal gas state: [BR]", + " [GRK]Delta[grk] A = A(N,V,T) - A[SUB]idealgas[sub](N,V,T) = kT [LN][CHEVRON][EXP]U[SUB]pot[sub]/kT[exp][chevron][ln][BR]", + " [GRK]Delta[grk] G = G(N,p,T) - G[SUB]idealgas[sub](N,p,T) = kT [LN][CHEVRON][EXP]U[SUB]pot[sub]/kT[exp][chevron][ln][BR]", + "where k is Boltzmann's constant, T is set by [TT]-fetemp[tt] and", + "the average is over the ensemble (or time in a trajectory).", + "Note that this is in principle", + "only correct when averaging over the whole (Boltzmann) ensemble", + "and using the potential energy. This also allows for an entropy", + "estimate using:[BR]", + " [GRK]Delta[grk] S(N,V,T) = S(N,V,T) - S[SUB]idealgas[sub](N,V,T) = ([CHEVRON]U[SUB]pot[sub][chevron] - [GRK]Delta[grk] A)/T[BR]", + " [GRK]Delta[grk] S(N,p,T) = S(N,p,T) - S[SUB]idealgas[sub](N,p,T) = ([CHEVRON]U[SUB]pot[sub][chevron] + pV - [GRK]Delta[grk] G)/T", + "[PAR]", + + "When a second energy file is specified ([TT]-f2[tt]), a free energy", + "difference is calculated [BR] dF = -kT [LN][CHEVRON][EXP]-(E[SUB]B[sub]-E[SUB]A[sub])/kT[exp][chevron][SUB]A[sub][ln] ,", + "where E[SUB]A[sub] and E[SUB]B[sub] are the energies from the first and second energy", + "files, and the average is over the ensemble A. The running average", + "of the free energy difference is printed to a file specified by [TT]-ravg[tt].", + "[BB]Note[bb] that the energies must both be calculated from the same trajectory." + + }; + static gmx_bool bSum = FALSE, bFee = FALSE, bPrAll = FALSE, bFluct = FALSE, bDriftCorr = FALSE; + static gmx_bool bDp = FALSE, bMutot = FALSE, bOrinst = FALSE, bOvec = FALSE, bFluctProps = FALSE; + static int skip = 0, nmol = 1, nbmin = 5, nbmax = 5; + static real reftemp = 300.0, ezero = 0; + t_pargs pa[] = { + { "-fee", FALSE, etBOOL, {&bFee}, + "Do a free energy estimate" }, + { "-fetemp", FALSE, etREAL, {&reftemp}, + "Reference temperature for free energy calculation" }, + { "-zero", FALSE, etREAL, {&ezero}, + "Subtract a zero-point energy" }, + { "-sum", FALSE, etBOOL, {&bSum}, + "Sum the energy terms selected rather than display them all" }, + { "-dp", FALSE, etBOOL, {&bDp}, + "Print energies in high precision" }, + { "-nbmin", FALSE, etINT, {&nbmin}, + "Minimum number of blocks for error estimate" }, + { "-nbmax", FALSE, etINT, {&nbmax}, + "Maximum number of blocks for error estimate" }, + { "-mutot", FALSE, etBOOL, {&bMutot}, + "Compute the total dipole moment from the components" }, + { "-skip", FALSE, etINT, {&skip}, + "Skip number of frames between data points" }, + { "-aver", FALSE, etBOOL, {&bPrAll}, + "Also print the exact average and rmsd stored in the energy frames (only when 1 term is requested)" }, + { "-nmol", FALSE, etINT, {&nmol}, + "Number of molecules in your sample: the energies are divided by this number" }, + { "-fluct_props", FALSE, etBOOL, {&bFluctProps}, + "Compute properties based on energy fluctuations, like heat capacity" }, + { "-driftcorr", FALSE, etBOOL, {&bDriftCorr}, + "Useful only for calculations of fluctuation properties. The drift in the observables will be subtracted before computing the fluctuation properties."}, + { "-fluc", FALSE, etBOOL, {&bFluct}, + "Calculate autocorrelation of energy fluctuations rather than energy itself" }, + { "-orinst", FALSE, etBOOL, {&bOrinst}, + "Analyse instantaneous orientation data" }, + { "-ovec", FALSE, etBOOL, {&bOvec}, + "Also plot the eigenvectors with [TT]-oten[tt]" } + }; + const char * drleg[] = { + "Running average", + "Instantaneous" + }; + static const char *setnm[] = { + "Pres-XX", "Pres-XY", "Pres-XZ", "Pres-YX", "Pres-YY", + "Pres-YZ", "Pres-ZX", "Pres-ZY", "Pres-ZZ", "Temperature", + "Volume", "Pressure" + }; + + FILE *out = NULL, *fp_pairs = NULL, *fort = NULL, *fodt = NULL, *foten = NULL; + FILE *fp_dhdl = NULL; + FILE **drout; + ener_file_t fp; + int timecheck = 0; + gmx_mtop_t mtop; + gmx_localtop_t *top = NULL; + t_inputrec ir; + t_energy **ee; + enerdata_t edat; + gmx_enxnm_t *enm = NULL; + t_enxframe *frame, *fr = NULL; + int cur = 0; +#define NEXT (1-cur) + int nre, teller, teller_disre, nfr; + gmx_int64_t start_step; + int nor = 0, nex = 0, norfr = 0, enx_i = 0; + real start_t; + real *bounds = NULL, *violaver = NULL, *oobs = NULL, *orient = NULL, *odrms = NULL; + int *index = NULL, *pair = NULL, norsel = 0, *orsel = NULL, *or_label = NULL; + int nbounds = 0, npairs; + gmx_bool bDisRe, bDRAll, bORA, bORT, bODA, bODR, bODT, bORIRE, bOTEN, bDHDL; + gmx_bool bFoundStart, bCont, bEDR, bVisco; + double sum, sumaver, sumt, ener, dbl; + double *time = NULL; + real Vaver; + int *set = NULL, i, j, k, nset, sss; + gmx_bool *bIsEner = NULL; + char **pairleg, **odtleg, **otenleg; + char **leg = NULL; + char **nms; + char *anm_j, *anm_k, *resnm_j, *resnm_k; + int resnr_j, resnr_k; + const char *orinst_sub = "@ subtitle \"instantaneous\"\n"; + char buf[256]; + output_env_t oenv; + t_enxblock *blk = NULL; + t_enxblock *blk_disre = NULL; + int ndisre = 0; + int dh_blocks = 0, dh_hists = 0, dh_samples = 0, dh_lambdas = 0; + + t_filenm fnm[] = { + { efEDR, "-f", NULL, ffREAD }, + { efEDR, "-f2", NULL, ffOPTRD }, + { efTPX, "-s", NULL, ffOPTRD }, + { efXVG, "-o", "energy", ffWRITE }, + { efXVG, "-viol", "violaver", ffOPTWR }, + { efXVG, "-pairs", "pairs", ffOPTWR }, + { efXVG, "-ora", "orienta", ffOPTWR }, + { efXVG, "-ort", "orientt", ffOPTWR }, + { efXVG, "-oda", "orideva", ffOPTWR }, + { efXVG, "-odr", "oridevr", ffOPTWR }, + { efXVG, "-odt", "oridevt", ffOPTWR }, + { efXVG, "-oten", "oriten", ffOPTWR }, + { efXVG, "-corr", "enecorr", ffOPTWR }, + { efXVG, "-vis", "visco", ffOPTWR }, + { efXVG, "-ravg", "runavgdf", ffOPTWR }, + { efXVG, "-odh", "dhdl", ffOPTWR } + }; +#define NFILE asize(fnm) + int npargs; + t_pargs *ppa; + + npargs = asize(pa); + ppa = add_acf_pargs(&npargs, pa); + if (!parse_common_args(&argc, argv, + PCA_CAN_VIEW | PCA_CAN_BEGIN | PCA_CAN_END | PCA_BE_NICE, + NFILE, fnm, npargs, ppa, asize(desc), desc, 0, NULL, &oenv)) + { + return 0; + } + + bDRAll = opt2bSet("-pairs", NFILE, fnm); + bDisRe = opt2bSet("-viol", NFILE, fnm) || bDRAll; + bORA = opt2bSet("-ora", NFILE, fnm); + bORT = opt2bSet("-ort", NFILE, fnm); + bODA = opt2bSet("-oda", NFILE, fnm); + bODR = opt2bSet("-odr", NFILE, fnm); + bODT = opt2bSet("-odt", NFILE, fnm); + bORIRE = bORA || bORT || bODA || bODR || bODT; + bOTEN = opt2bSet("-oten", NFILE, fnm); + bDHDL = opt2bSet("-odh", NFILE, fnm); + + nset = 0; + + snew(frame, 2); + fp = open_enx(ftp2fn(efEDR, NFILE, fnm), "r"); + do_enxnms(fp, &nre, &enm); + + Vaver = -1; + + bVisco = opt2bSet("-vis", NFILE, fnm); + + if ((!bDisRe) && (!bDHDL)) + { + if (bVisco) + { + nset = asize(setnm); + snew(set, nset); + /* This is nasty code... To extract Pres tensor, Volume and Temperature */ + for (j = 0; j < nset; j++) + { + for (i = 0; i < nre; i++) + { + if (strstr(enm[i].name, setnm[j])) + { + set[j] = i; + break; + } + } + if (i == nre) + { + if (gmx_strcasecmp(setnm[j], "Volume") == 0) + { + printf("Enter the box volume (" unit_volume "): "); + if (1 != scanf("%lf", &dbl)) + { + gmx_fatal(FARGS, "Error reading user input"); + } + Vaver = dbl; + } + else + { + gmx_fatal(FARGS, "Could not find term %s for viscosity calculation", + setnm[j]); + } + } + } + } + else + { + set = select_by_name(nre, enm, &nset); + } + /* Print all the different units once */ + sprintf(buf, "(%s)", enm[set[0]].unit); + for (i = 1; i < nset; i++) + { + for (j = 0; j < i; j++) + { + if (strcmp(enm[set[i]].unit, enm[set[j]].unit) == 0) + { + break; + } + } + if (j == i) + { + strcat(buf, ", ("); + strcat(buf, enm[set[i]].unit); + strcat(buf, ")"); + } + } + out = xvgropen(opt2fn("-o", NFILE, fnm), "Gromacs Energies", "Time (ps)", buf, + oenv); + + snew(leg, nset+1); + for (i = 0; (i < nset); i++) + { + leg[i] = enm[set[i]].name; + } + if (bSum) + { + leg[nset] = strdup("Sum"); + xvgr_legend(out, nset+1, (const char**)leg, oenv); + } + else + { + xvgr_legend(out, nset, (const char**)leg, oenv); + } + + snew(bIsEner, nset); + for (i = 0; (i < nset); i++) + { + bIsEner[i] = FALSE; + for (j = 0; (j <= F_ETOT); j++) + { + bIsEner[i] = bIsEner[i] || + (gmx_strcasecmp(interaction_function[j].longname, leg[i]) == 0); + } + } + + if (bPrAll && nset > 1) + { + gmx_fatal(FARGS, "Printing averages can only be done when a single set is selected"); + } + + time = NULL; + + if (bORIRE || bOTEN) + { + get_orires_parms(ftp2fn(efTPX, NFILE, fnm), &nor, &nex, &or_label, &oobs); + } + + if (bORIRE) + { + if (bOrinst) + { + enx_i = enxORI; + } + else + { + enx_i = enxOR; + } + + if (bORA || bODA) + { + snew(orient, nor); + } + if (bODR) + { + snew(odrms, nor); + } + if (bORT || bODT) + { + fprintf(stderr, "Select the orientation restraint labels you want (-1 is all)\n"); + fprintf(stderr, "End your selection with 0\n"); + j = -1; + orsel = NULL; + do + { + j++; + srenew(orsel, j+1); + if (1 != scanf("%d", &(orsel[j]))) + { + gmx_fatal(FARGS, "Error reading user input"); + } + } + while (orsel[j] > 0); + if (orsel[0] == -1) + { + fprintf(stderr, "Selecting all %d orientation restraints\n", nor); + norsel = nor; + srenew(orsel, nor); + for (i = 0; i < nor; i++) + { + orsel[i] = i; + } + } + else + { + /* Build the selection */ + norsel = 0; + for (i = 0; i < j; i++) + { + for (k = 0; k < nor; k++) + { + if (or_label[k] == orsel[i]) + { + orsel[norsel] = k; + norsel++; + break; + } + } + if (k == nor) + { + fprintf(stderr, "Orientation restraint label %d not found\n", + orsel[i]); + } + } + } + snew(odtleg, norsel); + for (i = 0; i < norsel; i++) + { + snew(odtleg[i], 256); + sprintf(odtleg[i], "%d", or_label[orsel[i]]); + } + if (bORT) + { + fort = xvgropen(opt2fn("-ort", NFILE, fnm), "Calculated orientations", + "Time (ps)", "", oenv); + if (bOrinst) + { + 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, - time, reftemp, &edat, ++ 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; +} diff --cc src/gromacs/legacyheaders/types/forcerec.h index 4adcbf0d2a,0000000000..71ffd993e1 mode 100644,000000..100644 --- a/src/gromacs/legacyheaders/types/forcerec.h +++ b/src/gromacs/legacyheaders/types/forcerec.h @@@ -1,499 -1,0 +1,501 @@@ +/* + * 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 diff --cc src/gromacs/legacyheaders/update.h index 398e505af5,0000000000..e589e725d1 mode 100644,000000..100644 --- a/src/gromacs/legacyheaders/update.h +++ b/src/gromacs/legacyheaders/update.h @@@ -1,248 -1,0 +1,249 @@@ +/* + * 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 */ diff --cc src/gromacs/mdlib/update.c index 0f96a35777,0000000000..3d6d040bd4 mode 100644,000000..100644 --- a/src/gromacs/mdlib/update.c +++ b/src/gromacs/mdlib/update.c @@@ -1,1997 -1,0 +1,2025 @@@ +/* + * 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 +#endif + + +#include +#include + +#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); + } +} + - 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]; - } - } - } - +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<hist); + } + if (state->flags & (1<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 */ - combine_forces(inputrec->nstcalclr, constr, inputrec, md, idef, cr, ++ combine_forces(upd, ++ inputrec->nstcalclr, constr, inputrec, md, idef, cr, + step, state, bMolPBC, - start, nrend, f, f_lr, nrnb); ++ 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; +} diff --cc src/programs/mdrun/md.c index 018daf5d0e,0000000000..60b01edd93 mode 100644,000000..100644 --- a/src/programs/mdrun/md.c +++ b/src/programs/mdrun/md.c @@@ -1,2009 -1,0 +1,2021 @@@ +/* + * 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 +#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" + " - 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<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, - f, 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, 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); + + /* 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<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; +}