*
* 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 <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 "gromacs/fileio/futil.h"
#include "gromacs/commandline/pargs.h"
#include "splitter.h"
-#include "gromacs/gmxana/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;
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,
&ls, &n_nonlin_vsite, &rlist_1x1);
/* Set the pair-list buffer size in ir */
- verletbuf_get_list_setup(FALSE, &ls);
+ verletbuf_get_list_setup(FALSE, FALSE, &ls);
calc_verlet_buffer_size(mtop, det(box), ir, buffer_temp,
&ls, &n_nonlin_vsite, &ir->rlist);
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;
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)
{
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",
{
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 &&
}
else
{
- buffer_temp = calc_temp(sys, ir, state.v);
+ buffer_temp = calc_temp(sys, ir, state->v);
}
if (buffer_temp > 0)
{
warning_note(wi, warn_buf);
}
- set_verlet_buffer(sys, ir, buffer_temp, state.box, wi);
+ 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