* 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 "grompp.h"
+#include "gmxpre.h"
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
+#include "grompp.h"
-#include <sys/types.h>
-#include <math.h>
-#include <string.h>
+#include <assert.h>
#include <errno.h>
#include <limits.h>
+#include <math.h>
+#include <string.h>
+
+#include <sys/types.h>
-#include "sysstuff.h"
-#include "smalloc.h"
-#include "macros.h"
-#include "string2.h"
-#include "readir.h"
-#include "toputil.h"
-#include "topio.h"
-#include "gromacs/fileio/confio.h"
-#include "readir.h"
-#include "symtab.h"
-#include "names.h"
-#include "grompp-impl.h"
-#include "random.h"
-#include "vec.h"
-#include "gromacs/fileio/futil.h"
#include "gromacs/commandline/pargs.h"
-#include "splitter.h"
-#include "sortwater.h"
-#include "convparm.h"
-#include "gmx_fatal.h"
-#include "warninp.h"
-#include "index.h"
+#include "gromacs/fileio/confio.h"
+#include "gromacs/fileio/enxio.h"
#include "gromacs/fileio/gmxfio.h"
-#include "gromacs/fileio/trnio.h"
#include "gromacs/fileio/tpxio.h"
+#include "gromacs/fileio/trnio.h"
#include "gromacs/fileio/trxio.h"
-#include "vsite_parm.h"
-#include "txtdump.h"
-#include "calcgrid.h"
-#include "add_par.h"
-#include "gromacs/fileio/enxio.h"
-#include "perf_est.h"
-#include "compute_io.h"
-#include "gpp_atomtype.h"
-#include "mtop_util.h"
-#include "genborn.h"
-#include "calc_verletbuf.h"
-
-#include "tomorse.h"
+#include "gromacs/gmxpreprocess/add_par.h"
+#include "gromacs/gmxpreprocess/calc_verletbuf.h"
+#include "gromacs/gmxpreprocess/compute_io.h"
+#include "gromacs/gmxpreprocess/convparm.h"
+#include "gromacs/gmxpreprocess/gen_maxwell_velocities.h"
+#include "gromacs/gmxpreprocess/gpp_atomtype.h"
+#include "gromacs/gmxpreprocess/grompp-impl.h"
+#include "gromacs/gmxpreprocess/readir.h"
+#include "gromacs/gmxpreprocess/sortwater.h"
+#include "gromacs/gmxpreprocess/tomorse.h"
+#include "gromacs/gmxpreprocess/topio.h"
+#include "gromacs/gmxpreprocess/toputil.h"
+#include "gromacs/gmxpreprocess/vsite_parm.h"
+#include "gromacs/imd/imd.h"
+#include "gromacs/legacyheaders/calcgrid.h"
+#include "gromacs/legacyheaders/genborn.h"
+#include "gromacs/legacyheaders/macros.h"
+#include "gromacs/legacyheaders/names.h"
+#include "gromacs/legacyheaders/perf_est.h"
+#include "gromacs/legacyheaders/splitter.h"
+#include "gromacs/legacyheaders/txtdump.h"
+#include "gromacs/legacyheaders/warninp.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/random/random.h"
+#include "gromacs/topology/mtop_util.h"
+#include "gromacs/topology/symtab.h"
+#include "gromacs/topology/topology.h"
+#include "gromacs/utility/cstringutil.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/futil.h"
+#include "gromacs/utility/smalloc.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;
+
+ 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];
+ }
- ref_T = 0;
+ 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,
+ calc_verlet_buffer_size(mtop, det(box), ir, buffer_temp,
&ls, &n_nonlin_vsite, &ir->rlist);
if (n_nonlin_vsite > 0)
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 },
{ efNDX, NULL, NULL, ffOPTRD },
{ efTOP, NULL, NULL, ffREAD },
{ efTOP, "-pp", "processed", ffOPTWR },
- { efTPX, "-o", NULL, ffWRITE },
+ { efTPR, "-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);
}
}
}
+ 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_vel(sys, state.v);
}
+ /* check for shells and inpurecs */
+ check_shells_inputrec(sys, ir, wi);
+
/* check masses */
check_mol(sys, 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);
+ }
}
}
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;
for (i = 0; i < efptNR; i++)
{
}
done_warning(wi, FARGS);
- write_tpx_state(ftp2fn(efTPX, NFILE, fnm), ir, &state, sys);
+ write_tpx_state(ftp2fn(efTPR, NFILE, fnm), ir, &state, sys);
+
+ /* Output IMD group, if bIMD is TRUE */
+ write_IMDgroup_to_file(ir->bIMD, ir, &state, sys, NFILE, fnm);
done_atomtype(atype);
done_mtop(sys, TRUE);
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff