*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#include <string.h>
#include <errno.h>
#include <limits.h>
+#include <assert.h>
#include "sysstuff.h"
-#include "smalloc.h"
+#include "gromacs/utility/smalloc.h"
#include "macros.h"
-#include "string2.h"
#include "readir.h"
#include "toputil.h"
#include "topio.h"
#include "symtab.h"
#include "names.h"
#include "grompp-impl.h"
-#include "random.h"
+#include "gromacs/random/random.h"
+#include "gromacs/gmxpreprocess/gen_maxwell_velocities.h"
#include "vec.h"
#include "gromacs/fileio/futil.h"
#include "gromacs/commandline/pargs.h"
#include "splitter.h"
-#include "sortwater.h"
+#include "gromacs/gmxpreprocess/sortwater.h"
#include "convparm.h"
#include "gmx_fatal.h"
#include "warninp.h"
#include "mtop_util.h"
#include "genborn.h"
#include "calc_verletbuf.h"
-
#include "tomorse.h"
+#include "gromacs/imd/imd.h"
+#include "gromacs/utility/cstringutil.h"
+
static int rm_interactions(int ifunc, int nrmols, t_molinfo mols[])
{
}
}
+static void check_shells_inputrec(gmx_mtop_t *mtop,
+ t_inputrec *ir,
+ warninp_t wi)
+{
+ gmx_mtop_atomloop_all_t aloop;
+ t_atom *atom;
+ int a, nshells = 0;
+ char warn_buf[STRLEN];
+
+ aloop = gmx_mtop_atomloop_all_init(mtop);
+ while (gmx_mtop_atomloop_all_next(aloop, &a, &atom))
+ {
+ if (atom->ptype == eptShell ||
+ atom->ptype == eptBond)
+ {
+ nshells++;
+ }
+ }
+ if (IR_TWINRANGE(*ir) && (nshells > 0))
+ {
+ snprintf(warn_buf, STRLEN,
+ "The combination of using shells and a twin-range cut-off is not supported");
+ warning_error(wi, warn_buf);
+ }
+ if ((nshells > 0) && (ir->nstcalcenergy != 1))
+ {
+ set_warning_line(wi, "unknown", -1);
+ snprintf(warn_buf, STRLEN,
+ "There are %d shells, changing nstcalcenergy from %d to 1",
+ nshells, ir->nstcalcenergy);
+ ir->nstcalcenergy = 1;
+ warning(wi, warn_buf);
+ }
+}
+
static gmx_bool nint_ftype(gmx_mtop_t *mtop, t_molinfo *mi, int ftype)
{
int nint, mb;
real *mass;
gmx_mtop_atomloop_all_t aloop;
t_atom *atom;
+ unsigned int useed;
snew(mass, state->natoms);
aloop = gmx_mtop_atomloop_all_init(sys);
mass[i] = atom->m;
}
+ useed = opts->seed;
if (opts->seed == -1)
{
- opts->seed = make_seed();
- fprintf(stderr, "Setting gen_seed to %d\n", opts->seed);
+ useed = (int)gmx_rng_make_seed();
+ fprintf(stderr, "Setting gen_seed to %u\n", useed);
}
- maxwell_speed(opts->tempi, opts->seed, sys, state->v);
+ maxwell_speed(opts->tempi, useed, sys, state->v);
stop_cm(stdout, state->natoms, mass, state->x, state->v);
sfree(mass);
if (rc_scaling != erscNO)
{
+ assert(npbcdim <= DIM);
+
for (mb = 0; mb < mtop->nmolblock; mb++)
{
molb = &mtop->molblock[mb];
int i, j;
read_posres (mtop, mi, FALSE, fnA, rc_scaling, ePBC, com, wi);
- if (strcmp(fnA, fnB) != 0)
- {
- read_posres(mtop, mi, TRUE, fnB, rc_scaling, ePBC, comB, wi);
- }
+ /* It is safer to simply read the b-state posres rather than trying
+ * to be smart and copy the positions.
+ */
+ read_posres(mtop, mi, TRUE, fnB, rc_scaling, ePBC, comB, wi);
}
static void set_wall_atomtype(gpp_atomtype_t at, t_gromppopts *opts,
}
-static void set_verlet_buffer(const gmx_mtop_t *mtop,
- t_inputrec *ir,
- matrix box,
- warninp_t wi)
+static real calc_temp(const gmx_mtop_t *mtop,
+ const t_inputrec *ir,
+ rvec *v)
{
- real ref_T;
- int i;
- verletbuf_list_setup_t ls;
- real rlist_1x1;
- int n_nonlin_vsite;
- char warn_buf[STRLEN];
+ double sum_mv2;
+ gmx_mtop_atomloop_all_t aloop;
+ t_atom *atom;
+ int a;
+ int nrdf, g;
+
+ sum_mv2 = 0;
- ref_T = 0;
+ aloop = gmx_mtop_atomloop_all_init(mtop);
+ while (gmx_mtop_atomloop_all_next(aloop, &a, &atom))
+ {
+ sum_mv2 += atom->m*norm2(v[a]);
+ }
+
+ nrdf = 0;
+ for (g = 0; g < ir->opts.ngtc; g++)
+ {
+ nrdf += ir->opts.nrdf[g];
+ }
+
+ return sum_mv2/(nrdf*BOLTZ);
+}
+
+static real get_max_reference_temp(const t_inputrec *ir,
+ warninp_t wi)
+{
+ real ref_t;
+ int i;
+ gmx_bool bNoCoupl;
+
+ ref_t = 0;
+ bNoCoupl = FALSE;
for (i = 0; i < ir->opts.ngtc; i++)
{
- if (ir->opts.ref_t[i] < 0)
+ if (ir->opts.tau_t[i] < 0)
{
- warning(wi, "Some atom groups do not use temperature coupling. This cannot be accounted for in the energy error estimation for the Verlet buffer size. The energy error and the Verlet buffer might be underestimated.");
+ bNoCoupl = TRUE;
}
else
{
- ref_T = max(ref_T, ir->opts.ref_t[i]);
+ ref_t = max(ref_t, ir->opts.ref_t[i]);
}
}
- printf("Determining Verlet buffer for a tolerance of %g kJ/mol/ps at %g K\n", ir->verletbuf_tol, ref_T);
-
- for (i = 0; i < ir->opts.ngtc; i++)
+ if (bNoCoupl)
{
- if (ir->opts.ref_t[i] >= 0 && ir->opts.ref_t[i] != ref_T)
- {
- sprintf(warn_buf, "ref_T for group of %.1f DOFs is %g K, which is smaller than the maximum of %g K used for the buffer size calculation. The buffer size might be on the conservative (large) side.",
- ir->opts.nrdf[i], ir->opts.ref_t[i], ref_T);
- warning_note(wi, warn_buf);
- }
+ char buf[STRLEN];
+
+ sprintf(buf, "Some temperature coupling groups do not use temperature coupling. We will assume their temperature is not more than %.3f K. If their temperature is higher, the energy error and the Verlet buffer might be underestimated.",
+ ref_t);
+ warning(wi, buf);
}
+ return ref_t;
+}
+
+static void set_verlet_buffer(const gmx_mtop_t *mtop,
+ t_inputrec *ir,
+ real buffer_temp,
+ matrix box,
+ warninp_t wi)
+{
+ int i;
+ verletbuf_list_setup_t ls;
+ real rlist_1x1;
+ int n_nonlin_vsite;
+ char warn_buf[STRLEN];
+
+ printf("Determining Verlet buffer for a tolerance of %g kJ/mol/ps at %g K\n", ir->verletbuf_tol, buffer_temp);
+
/* Calculate the buffer size for simple atom vs atoms list */
ls.cluster_size_i = 1;
ls.cluster_size_j = 1;
- calc_verlet_buffer_size(mtop, det(box), ir,
+ calc_verlet_buffer_size(mtop, det(box), ir, buffer_temp,
&ls, &n_nonlin_vsite, &rlist_1x1);
/* Set the pair-list buffer size in ir */
- verletbuf_get_list_setup(FALSE, &ls);
- calc_verlet_buffer_size(mtop, det(box), ir,
+ verletbuf_get_list_setup(FALSE, FALSE, &ls);
+ calc_verlet_buffer_size(mtop, det(box), ir, buffer_temp,
&ls, &n_nonlin_vsite, &ir->rlist);
if (n_nonlin_vsite > 0)
ls.cluster_size_i, ls.cluster_size_j,
ir->rlist, ir->rlist-max(ir->rvdw, ir->rcoulomb));
+ printf("Note that mdrun will redetermine rlist based on the actual pair-list setup\n");
+
if (sqr(ir->rlistlong) >= max_cutoff2(ir->ePBC, box))
{
gmx_fatal(FARGS, "The pair-list cut-off (%g nm) is longer than half the shortest box vector or longer than the smallest box diagonal element (%g nm). Increase the box size or decrease nstlist or increase verlet-buffer-tolerance.", ir->rlistlong, sqrt(max_cutoff2(ir->ePBC, box)));
"with [TT]-cpi[tt]. If you wish to change the ensemble or things",
"like output frequency, then supplying the checkpoint file to",
"[THISMODULE] with [TT]-t[tt] along with a new [TT].mdp[tt] file",
- "with [TT]-f[tt] is the recommended procedure.[PAR]",
+ "with [TT]-f[tt] is the recommended procedure. Actually preserving",
+ "the ensemble (if possible) still requires passing the checkpoint",
+ "file to [gmx-mdrun] [TT]-cpi[tt].[PAR]",
"By default, all bonded interactions which have constant energy due to",
"virtual site constructions will be removed. If this constant energy is",
t_inputrec *ir;
int natoms, nvsite, comb, mt;
t_params *plist;
- t_state state;
+ t_state *state;
matrix box;
real max_spacing, fudgeQQ;
double reppow;
gmx_bool bVerbose = FALSE;
warninp_t wi;
char warn_buf[STRLEN];
+ unsigned int useed;
+ t_atoms IMDatoms; /* Atoms to be operated on interactively (IMD) */
t_filenm fnm[] = {
{ efMDP, NULL, NULL, ffREAD },
{ efTPX, "-o", NULL, ffWRITE },
{ efTRN, "-t", NULL, ffOPTRD },
{ efEDR, "-e", NULL, ffOPTRD },
- { efTRN, "-ref", "rotref", ffOPTRW }
+ /* This group is needed by the VMD viewer as the start configuration for IMD sessions: */
+ { efGRO, "-imd", "imdgroup", ffOPTWR },
+ { efTRN, "-ref", "rotref", ffOPTRW }
};
#define NFILE asize(fnm)
if (ir->ld_seed == -1)
{
- ir->ld_seed = make_seed();
- fprintf(stderr, "Setting the LD random seed to %d\n", ir->ld_seed);
+ ir->ld_seed = (gmx_int64_t)gmx_rng_make_seed();
+ fprintf(stderr, "Setting the LD random seed to %"GMX_PRId64 "\n", ir->ld_seed);
}
if (ir->expandedvals->lmc_seed == -1)
{
- ir->expandedvals->lmc_seed = make_seed();
+ ir->expandedvals->lmc_seed = (int)gmx_rng_make_seed();
fprintf(stderr, "Setting the lambda MC random seed to %d\n", ir->expandedvals->lmc_seed);
}
{
gmx_fatal(FARGS, "%s does not exist", fn);
}
+ snew(state, 1);
new_status(fn, opt2fn_null("-pp", NFILE, fnm), opt2fn("-c", NFILE, fnm),
- opts, ir, bZero, bGenVel, bVerbose, &state,
+ opts, ir, bZero, bGenVel, bVerbose, state,
atype, sys, &nmi, &mi, plist, &comb, &reppow, &fudgeQQ,
opts->bMorse,
wi);
}
}
+ if (nvsite && ir->eI == eiNM)
+ {
+ 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 (ir->cutoff_scheme == ecutsVERLET)
{
fprintf(stderr, "Removing all charge groups because cutoff-scheme=%s\n",
/* Remove all charge groups */
gmx_mtop_remove_chargegroups(sys);
-
- if (EVDW_PME(ir->vdwtype))
- {
- gmx_fatal(FARGS, "LJ-PME not implemented together with verlet-scheme!");
- }
}
if (count_constraints(sys, mi, wi) && (ir->eConstrAlg == econtSHAKE))
{
warning_error(wi, "AdResS contant weighting function should be between 0 and 1\n\n");
}
- /** \TODO check size of ex+hy width against box size */
+ /** TODO check size of ex+hy width against box size */
}
/* Check for errors in the input now, since they might cause problems
}
/* If we are using CMAP, setup the pre-interpolation grid */
- if (plist->ncmap > 0)
+ if (plist[F_CMAP].ncmap > 0)
{
- init_cmap_grid(&sys->ffparams.cmap_grid, plist->nc, plist->grid_spacing);
- setup_cmap(plist->grid_spacing, plist->nc, plist->cmap, &sys->ffparams.cmap_grid);
+ init_cmap_grid(&sys->ffparams.cmap_grid, plist[F_CMAP].nc, plist[F_CMAP].grid_spacing);
+ setup_cmap(plist[F_CMAP].grid_spacing, plist[F_CMAP].nc, plist[F_CMAP].cmap, &sys->ffparams.cmap_grid);
}
set_wall_atomtype(atype, opts, ir, wi);
/* Check velocity for virtual sites and shells */
if (bGenVel)
{
- check_vel(sys, state.v);
+ check_vel(sys, state->v);
}
+ /* check for shells and inpurecs */
+ check_shells_inputrec(sys, ir, wi);
+
/* check masses */
check_mol(sys, wi);
}
do_index(mdparin, ftp2fn_null(efNDX, NFILE, fnm),
sys, bVerbose, ir,
- bGenVel ? state.v : NULL,
+ bGenVel ? state->v : NULL,
wi);
if (ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0 &&
ir->nstlist > 1)
{
- if (EI_DYNAMICS(ir->eI) &&
- !(EI_MD(ir->eI) && ir->etc == etcNO) &&
- inputrec2nboundeddim(ir) == 3)
+ if (EI_DYNAMICS(ir->eI) && inputrec2nboundeddim(ir) == 3)
{
- set_verlet_buffer(sys, ir, state.box, wi);
+ real buffer_temp;
+
+ if (EI_MD(ir->eI) && ir->etc == etcNO)
+ {
+ if (bGenVel)
+ {
+ buffer_temp = opts->tempi;
+ }
+ else
+ {
+ buffer_temp = calc_temp(sys, ir, state->v);
+ }
+ if (buffer_temp > 0)
+ {
+ sprintf(warn_buf, "NVE simulation: will use the initial temperature of %.3f K for determining the Verlet buffer size", buffer_temp);
+ warning_note(wi, warn_buf);
+ }
+ else
+ {
+ sprintf(warn_buf, "NVE simulation with an initial temperature of zero: will use a Verlet buffer of %d%%. Check your energy drift!",
+ (int)(verlet_buffer_ratio_NVE_T0*100 + 0.5));
+ warning_note(wi, warn_buf);
+ }
+ }
+ else
+ {
+ buffer_temp = get_max_reference_temp(ir, wi);
+ }
+
+ if (EI_MD(ir->eI) && ir->etc == etcNO && buffer_temp == 0)
+ {
+ /* NVE with initial T=0: we add a fixed ratio to rlist.
+ * Since we don't actually use verletbuf_tol, we set it to -1
+ * so it can't be misused later.
+ */
+ ir->rlist *= 1.0 + verlet_buffer_ratio_NVE_T0;
+ ir->verletbuf_tol = -1;
+ }
+ else
+ {
+ /* We warn for NVE simulations with >1(.1)% drift tolerance */
+ const real drift_tol = 0.01;
+ real ener_runtime;
+
+ /* We use 2 DOF per atom = 2kT pot+kin energy, to be on
+ * the safe side with constraints (without constraints: 3 DOF).
+ */
+ ener_runtime = 2*BOLTZ*buffer_temp/(ir->nsteps*ir->delta_t);
+
+ if (EI_MD(ir->eI) && ir->etc == etcNO && ir->nstlist > 1 &&
+ ir->nsteps > 0 &&
+ ir->verletbuf_tol > 1.1*drift_tol*ener_runtime)
+ {
+ sprintf(warn_buf, "You are using a Verlet buffer tolerance of %g kJ/mol/ps for an NVE simulation of length %g ps, which can give a final drift of %d%%. For conserving energy to %d%%, you might need to set verlet-buffer-tolerance to %.1e.",
+ ir->verletbuf_tol, ir->nsteps*ir->delta_t,
+ (int)(ir->verletbuf_tol/ener_runtime*100 + 0.5),
+ (int)(100*drift_tol + 0.5),
+ drift_tol*ener_runtime);
+ warning_note(wi, warn_buf);
+ }
+
+ set_verlet_buffer(sys, ir, buffer_temp, state->box, wi);
+ }
}
}
/* Init the temperature coupling state */
- init_gtc_state(&state, ir->opts.ngtc, 0, ir->opts.nhchainlength); /* need to add nnhpres here? */
+ init_gtc_state(state, ir->opts.ngtc, 0, ir->opts.nhchainlength); /* need to add nnhpres here? */
if (bVerbose)
{
fprintf(stderr, "getting data from old trajectory ...\n");
}
cont_status(ftp2fn(efTRN, NFILE, fnm), ftp2fn_null(efEDR, NFILE, fnm),
- bNeedVel, bGenVel, fr_time, ir, &state, sys, oenv);
+ bNeedVel, bGenVel, fr_time, ir, state, sys, oenv);
}
if (ir->ePBC == epbcXY && ir->nwall != 2)
{
- clear_rvec(state.box[ZZ]);
+ clear_rvec(state->box[ZZ]);
}
if (ir->cutoff_scheme != ecutsVERLET && ir->rlist > 0)
{
set_warning_line(wi, mdparin, -1);
- check_chargegroup_radii(sys, ir, state.x, wi);
+ check_chargegroup_radii(sys, ir, state->x, wi);
}
if (EEL_FULL(ir->coulombtype) || EVDW_PME(ir->vdwtype))
{
/* Calculate the optimal grid dimensions */
- copy_mat(state.box, box);
+ copy_mat(state->box, box);
if (ir->ePBC == epbcXY && ir->nwall == 2)
{
svmul(ir->wall_ewald_zfac, box[ZZ], box[ZZ]);
potentially conflict if not handled correctly. */
if (ir->efep != efepNO)
{
- if (EVDW_PME(ir->vdwtype))
- {
- gmx_fatal(FARGS, "LJ-PME not implemented together with free energy calculations!");
- }
-
- state.fep_state = ir->fepvals->init_fep_state;
+ state->fep_state = ir->fepvals->init_fep_state;
for (i = 0; i < efptNR; i++)
{
/* init_lambda trumps state definitions*/
if (ir->fepvals->init_lambda >= 0)
{
- state.lambda[i] = ir->fepvals->init_lambda;
+ state->lambda[i] = ir->fepvals->init_lambda;
}
else
{
}
else
{
- state.lambda[i] = ir->fepvals->all_lambda[i][state.fep_state];
+ state->lambda[i] = ir->fepvals->all_lambda[i][state->fep_state];
}
}
}
if (ir->ePull != epullNO)
{
- set_pull_init(ir, sys, state.x, state.box, state.lambda[efptMASS], oenv, opts->pull_start);
+ set_pull_init(ir, sys, state->x, state->box, state->lambda[efptMASS], oenv, opts->pull_start);
}
if (ir->bRot)
{
- set_reference_positions(ir->rot, state.x, state.box,
+ set_reference_positions(ir->rot, state->x, state->box,
opt2fn("-ref", NFILE, fnm), opt2bSet("-ref", NFILE, fnm),
wi);
}
if (EEL_PME(ir->coulombtype))
{
- float ratio = pme_load_estimate(sys, ir, state.box);
+ float ratio = pme_load_estimate(sys, ir, state->box);
fprintf(stderr, "Estimate for the relative computational load of the PME mesh part: %.2f\n", ratio);
/* With free energy we might need to do PME both for the A and B state
* charges. This will double the cost, but the optimal performance will
}
done_warning(wi, FARGS);
- write_tpx_state(ftp2fn(efTPX, NFILE, fnm), ir, &state, sys);
+ write_tpx_state(ftp2fn(efTPX, NFILE, fnm), ir, state, sys);
+
+ /* Output IMD group, if bIMD is TRUE */
+ write_IMDgroup_to_file(ir->bIMD, ir, state, sys, NFILE, fnm);
+ done_state(state);
+ sfree(state);
done_atomtype(atype);
done_mtop(sys, TRUE);
done_inputrec_strings();
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff