-/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
+/*
+ * This file is part of the GROMACS molecular simulation package.
*
- *
- * This source code is part of
- *
- * G R O M A C S
- *
- * GROningen MAchine for Chemical Simulations
- *
- * VERSION 3.2.0
- * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2004, The GROMACS development team,
- * check out http://www.gromacs.org for more information.
-
- * This program is free software; you can redistribute it and/or
- * modify it under the terms of the GNU General Public License
- * as published by the Free Software Foundation; either version 2
+ * 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.
*
- * If you want to redistribute modifications, 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 www.gromacs.org.
+ * 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.
*
- * To help us fund GROMACS development, we humbly ask that you cite
- * the papers on the package - you can find them in the top README file.
+ * 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.
*
- * For more info, check our website at http://www.gromacs.org
+ * 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.
*
- * And Hey:
- * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
*/
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
+#include "gmxpre.h"
#include <ctype.h>
-#include <stdlib.h>
#include <limits.h>
-#include "sysstuff.h"
-#include "smalloc.h"
-#include "typedefs.h"
-#include "physics.h"
-#include "names.h"
-#include "gmx_fatal.h"
-#include "macros.h"
-#include "index.h"
-#include "symtab.h"
-#include "string2.h"
-#include "readinp.h"
-#include "warninp.h"
+#include <stdlib.h>
+
+#include "gromacs/legacyheaders/typedefs.h"
+#include "gromacs/math/units.h"
+#include "gromacs/legacyheaders/names.h"
+#include "gromacs/legacyheaders/macros.h"
+#include "gromacs/topology/index.h"
+#include "gromacs/utility/cstringutil.h"
+#include "gromacs/legacyheaders/readinp.h"
+#include "gromacs/legacyheaders/warninp.h"
#include "readir.h"
#include "toputil.h"
-#include "index.h"
-#include "network.h"
-#include "vec.h"
-#include "pbc.h"
-#include "mtop_util.h"
-#include "chargegroup.h"
-#include "inputrec.h"
+#include "gromacs/legacyheaders/network.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/topology/mtop_util.h"
+#include "gromacs/legacyheaders/chargegroup.h"
+#include "gromacs/legacyheaders/inputrec.h"
+#include "calc_verletbuf.h"
+
+#include "gromacs/topology/block.h"
+#include "gromacs/topology/symtab.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/smalloc.h"
#define MAXPTR 254
#define NOGID 255
-#define MAXLAMBDAS 1024
/* Resource parameters
* Do not change any of these until you read the instruction
* message.
*/
-static char tcgrps[STRLEN], tau_t[STRLEN], ref_t[STRLEN],
- acc[STRLEN], accgrps[STRLEN], freeze[STRLEN], frdim[STRLEN],
- energy[STRLEN], user1[STRLEN], user2[STRLEN], vcm[STRLEN], xtc_grps[STRLEN],
- couple_moltype[STRLEN], orirefitgrp[STRLEN], egptable[STRLEN], egpexcl[STRLEN],
- wall_atomtype[STRLEN], wall_density[STRLEN], deform[STRLEN], QMMM[STRLEN];
-static char fep_lambda[efptNR][STRLEN];
-static char lambda_weights[STRLEN];
-static char **pull_grp;
-static char **rot_grp;
-static char anneal[STRLEN], anneal_npoints[STRLEN],
- anneal_time[STRLEN], anneal_temp[STRLEN];
-static char QMmethod[STRLEN], QMbasis[STRLEN], QMcharge[STRLEN], QMmult[STRLEN],
- bSH[STRLEN], CASorbitals[STRLEN], CASelectrons[STRLEN], SAon[STRLEN],
- SAoff[STRLEN], SAsteps[STRLEN], bTS[STRLEN], bOPT[STRLEN];
-static char efield_x[STRLEN], efield_xt[STRLEN], efield_y[STRLEN],
- efield_yt[STRLEN], efield_z[STRLEN], efield_zt[STRLEN];
+typedef struct t_inputrec_strings
+{
+ char tcgrps[STRLEN], tau_t[STRLEN], ref_t[STRLEN],
+ acc[STRLEN], accgrps[STRLEN], freeze[STRLEN], frdim[STRLEN],
+ energy[STRLEN], user1[STRLEN], user2[STRLEN], vcm[STRLEN], x_compressed_groups[STRLEN],
+ couple_moltype[STRLEN], orirefitgrp[STRLEN], egptable[STRLEN], egpexcl[STRLEN],
+ wall_atomtype[STRLEN], wall_density[STRLEN], deform[STRLEN], QMMM[STRLEN],
+ imd_grp[STRLEN];
+ char fep_lambda[efptNR][STRLEN];
+ char lambda_weights[STRLEN];
+ char **pull_grp;
+ char **rot_grp;
+ char anneal[STRLEN], anneal_npoints[STRLEN],
+ anneal_time[STRLEN], anneal_temp[STRLEN];
+ char QMmethod[STRLEN], QMbasis[STRLEN], QMcharge[STRLEN], QMmult[STRLEN],
+ bSH[STRLEN], CASorbitals[STRLEN], CASelectrons[STRLEN], SAon[STRLEN],
+ SAoff[STRLEN], SAsteps[STRLEN], bTS[STRLEN], bOPT[STRLEN];
+ char efield_x[STRLEN], efield_xt[STRLEN], efield_y[STRLEN],
+ efield_yt[STRLEN], efield_z[STRLEN], efield_zt[STRLEN];
+
+} gmx_inputrec_strings;
+
+static gmx_inputrec_strings *is = NULL;
+
+void init_inputrec_strings()
+{
+ if (is)
+ {
+ gmx_incons("Attempted to call init_inputrec_strings before calling done_inputrec_strings. Only one inputrec (i.e. .mdp file) can be parsed at a time.");
+ }
+ snew(is, 1);
+}
+
+void done_inputrec_strings()
+{
+ sfree(is);
+ is = NULL;
+}
+
+static char swapgrp[STRLEN], splitgrp0[STRLEN], splitgrp1[STRLEN], solgrp[STRLEN];
enum {
egrptpALL, /* All particles have to be a member of a group. */
void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
warninp_t wi)
-/* Check internal consistency */
+/* Check internal consistency.
+ * NOTE: index groups are not set here yet, don't check things
+ * like temperature coupling group options here, but in triple_check
+ */
{
/* Strange macro: first one fills the err_buf, and then one can check
* the condition, which will print the message and increase the error
warning_error(wi, "rvdw should be >= 0");
}
if (ir->rlist < 0 &&
- !(ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_drift > 0))
+ !(ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0))
{
warning_error(wi, "rlist should be >= 0");
}
if (ir->cutoff_scheme == ecutsGROUP)
{
+ warning_note(wi,
+ "The group cutoff scheme is deprecated in Gromacs 5.0 and will be removed in a future "
+ "release when all interaction forms are supported for the verlet scheme. The verlet "
+ "scheme already scales better, and it is compatible with GPUs and other accelerators.");
+
/* BASIC CUT-OFF STUFF */
if (ir->rlist == 0 ||
- !((EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype) && ir->rcoulomb > ir->rlist) ||
- (EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype) && ir->rvdw > ir->rlist)))
+ !((ir_coulomb_might_be_zero_at_cutoff(ir) && ir->rcoulomb > ir->rlist) ||
+ (ir_vdw_might_be_zero_at_cutoff(ir) && ir->rvdw > ir->rlist)))
{
/* No switched potential and/or no twin-range:
* we can set the long-range cut-off to the maximum of the other cut-offs.
{
warning_error(wi, "With Verlet lists rcoulomb!=rvdw is not supported");
}
- if (ir->vdwtype != evdwCUT)
+ if (ir->vdwtype == evdwSHIFT || ir->vdwtype == evdwSWITCH)
+ {
+ if (ir->vdw_modifier == eintmodNONE ||
+ ir->vdw_modifier == eintmodPOTSHIFT)
+ {
+ ir->vdw_modifier = (ir->vdwtype == evdwSHIFT ? eintmodFORCESWITCH : eintmodPOTSWITCH);
+
+ sprintf(warn_buf, "Replacing vdwtype=%s by the equivalent combination of vdwtype=%s and vdw_modifier=%s", evdw_names[ir->vdwtype], evdw_names[evdwCUT], eintmod_names[ir->vdw_modifier]);
+ warning_note(wi, warn_buf);
+
+ ir->vdwtype = evdwCUT;
+ }
+ else
+ {
+ sprintf(warn_buf, "Unsupported combination of vdwtype=%s and vdw_modifier=%s", evdw_names[ir->vdwtype], eintmod_names[ir->vdw_modifier]);
+ warning_error(wi, warn_buf);
+ }
+ }
+
+ if (!(ir->vdwtype == evdwCUT || ir->vdwtype == evdwPME))
{
- warning_error(wi, "With Verlet lists only cut-off LJ interactions are supported");
+ warning_error(wi, "With Verlet lists only cut-off and PME LJ interactions are supported");
}
if (!(ir->coulombtype == eelCUT ||
(EEL_RF(ir->coulombtype) && ir->coulombtype != eelRF_NEC) ||
{
warning_error(wi, "With Verlet lists only cut-off, reaction-field, PME and Ewald electrostatics are supported");
}
+ if (!(ir->coulomb_modifier == eintmodNONE ||
+ ir->coulomb_modifier == eintmodPOTSHIFT))
+ {
+ sprintf(warn_buf, "coulomb_modifier=%s is not supported with the Verlet cut-off scheme", eintmod_names[ir->coulomb_modifier]);
+ warning_error(wi, warn_buf);
+ }
if (ir->nstlist <= 0)
{
rc_max = max(ir->rvdw, ir->rcoulomb);
- if (ir->verletbuf_drift <= 0)
+ if (ir->verletbuf_tol <= 0)
{
- if (ir->verletbuf_drift == 0)
+ if (ir->verletbuf_tol == 0)
{
- warning_error(wi, "Can not have an energy drift of exactly 0");
+ warning_error(wi, "Can not have Verlet buffer tolerance of exactly 0");
}
if (ir->rlist < rc_max)
{
if (ir->rlist > rc_max)
{
- warning_note(wi, "You have set rlist larger than the interaction cut-off, but you also have verlet-buffer-drift > 0. Will set rlist using verlet-buffer-drift.");
+ warning_note(wi, "You have set rlist larger than the interaction cut-off, but you also have verlet-buffer-tolerance > 0. Will set rlist using verlet-buffer-tolerance.");
}
if (ir->nstlist == 1)
{
if (EI_DYNAMICS(ir->eI))
{
- if (EI_MD(ir->eI) && ir->etc == etcNO)
- {
- warning_error(wi, "Temperature coupling is required for calculating rlist using the energy drift with verlet-buffer-drift > 0. Either use temperature coupling or set rlist yourself together with verlet-buffer-drift = -1.");
- }
-
if (inputrec2nboundeddim(ir) < 3)
{
- warning_error(wi, "The box volume is required for calculating rlist from the energy drift with verlet-buffer-drift > 0. You are using at least one unbounded dimension, so no volume can be computed. Either use a finite box, or set rlist yourself together with verlet-buffer-drift = -1.");
+ warning_error(wi, "The box volume is required for calculating rlist from the energy drift with verlet-buffer-tolerance > 0. You are using at least one unbounded dimension, so no volume can be computed. Either use a finite box, or set rlist yourself together with verlet-buffer-tolerance = -1.");
}
/* Set rlist temporarily so we can continue processing */
ir->rlist = rc_max;
else
{
/* Set the buffer to 5% of the cut-off */
- ir->rlist = 1.05*rc_max;
+ ir->rlist = (1.0 + verlet_buffer_ratio_nodynamics)*rc_max;
}
}
}
/* find the smallest of ( nstenergy, nstdhdl ) */
if (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
- (ir->fepvals->nstdhdl < ir->nstenergy) )
+ (ir->nstenergy == 0 || ir->fepvals->nstdhdl < ir->nstenergy))
{
min_nst = ir->fepvals->nstdhdl;
min_name = nstdh;
}
}
+ if (ir->nsteps == 0 && !ir->bContinuation)
+ {
+ warning_note(wi, "For a correct single-point energy evaluation with nsteps = 0, use continuation = yes to avoid constraining the input coordinates.");
+ }
+
/* LD STUFF */
if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
ir->bContinuation && ir->ld_seed != -1)
CHECK(ir->nstlist <= 0);
sprintf(err_buf, "TPI does not work with full electrostatics other than PME");
CHECK(EEL_FULL(ir->coulombtype) && !EEL_PME(ir->coulombtype));
+ sprintf(err_buf, "TPI does not work (yet) with the Verlet cut-off scheme");
+ CHECK(ir->cutoff_scheme == ecutsVERLET);
}
/* SHAKE / LINCS */
sprintf(err_buf, "Can't use postive delta-lambda (%g) with expanded ensemble simulations", fep->delta_lambda);
CHECK(fep->delta_lambda > 0 && (ir->efep == efepEXPANDED));
+ sprintf(err_buf, "Can only use expanded ensemble with md-vv for now; should be supported for other integrators in 5.0");
+ CHECK(!(EI_VV(ir->eI)) && (ir->efep == efepEXPANDED));
+
sprintf(err_buf, "Free-energy not implemented for Ewald");
CHECK(ir->coulombtype == eelEWALD);
}
if (ir->nstlist > 0)
{
- warning_note(wi, "Simulating without cut-offs is usually (slightly) faster with nstlist=0, nstype=simple and particle decomposition");
+ warning_note(wi, "Simulating without cut-offs can be (slightly) faster with nstlist=0, nstype=simple and only one MPI rank");
}
}
sprintf(err_buf, "%s temperature control not supported for integrator %s.", etcoupl_names[ir->etc], ei_names[ir->eI]);
CHECK(!(EI_VV(ir->eI)));
- for (i = 0; i < ir->opts.ngtc; i++)
- {
- sprintf(err_buf, "all tau_t must currently be equal using Andersen temperature control, violated for group %d", i);
- CHECK(ir->opts.tau_t[0] != ir->opts.tau_t[i]);
- sprintf(err_buf, "all tau_t must be postive using Andersen temperature control, tau_t[%d]=%10.6f",
- i, ir->opts.tau_t[i]);
- CHECK(ir->opts.tau_t[i] < 0);
- }
if (ir->nstcomm > 0 && (ir->etc == etcANDERSEN))
{
sprintf(warn_buf, "Center of mass removal not necessary for %s. All velocities of coupled groups are rerandomized periodically, so flying ice cube errors will not occur.", etcoupl_names[ir->etc]);
sprintf(err_buf, "nstcomm must be 1, not %d for %s, as velocities of atoms in coupled groups are randomized every time step", ir->nstcomm, etcoupl_names[ir->etc]);
CHECK(ir->nstcomm > 1 && (ir->etc == etcANDERSEN));
-
- for (i = 0; i < ir->opts.ngtc; i++)
- {
- int nsteps = (int)(ir->opts.tau_t[i]/ir->delta_t);
- sprintf(err_buf, "tau_t/delta_t for group %d for temperature control method %s must be a multiple of nstcomm (%d), as velocities of atoms in coupled groups are randomized every time step. The input tau_t (%8.3f) leads to %d steps per randomization", i, etcoupl_names[ir->etc], ir->nstcomm, ir->opts.tau_t[i], nsteps);
- CHECK((nsteps % ir->nstcomm) && (ir->etc == etcANDERSENMASSIVE));
- }
}
+
if (ir->etc == etcBERENDSEN)
{
sprintf(warn_buf, "The %s thermostat does not generate the correct kinetic energy distribution. You might want to consider using the %s thermostat.",
sprintf(err_buf, "tau-p must be > 0 instead of %g\n", ir->tau_p);
CHECK(ir->tau_p <= 0);
- if (ir->tau_p/dt_pcoupl < pcouple_min_integration_steps(ir->epc))
+ if (ir->tau_p/dt_pcoupl < pcouple_min_integration_steps(ir->epc) - 10*GMX_REAL_EPS)
{
sprintf(warn_buf, "For proper integration of the %s barostat, tau-p (%g) should be at least %d times larger than nstpcouple*dt (%g)",
EPCOUPLTYPE(ir->epc), ir->tau_p, pcouple_min_integration_steps(ir->epc), dt_pcoupl);
}
}
}
+ else
+ {
+ if (ir->epc == epcMTTK)
+ {
+ warning_error(wi, "MTTK pressure coupling requires a Velocity-verlet integrator");
+ }
+ }
/* ELECTROSTATICS */
/* More checks are in triple check (grompp.c) */
ir->epsilon_r = 1.0;
}
- if (getenv("GALACTIC_DYNAMICS") == NULL)
+ if (getenv("GMX_DO_GALACTIC_DYNAMICS") == NULL)
{
sprintf(err_buf, "epsilon-r must be >= 0 instead of %g\n", ir->epsilon_r);
CHECK(ir->epsilon_r < 0);
* means the interaction is zero outside rcoulomb, but it helps to
* provide accurate energy conservation.
*/
- if (EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype))
+ if (ir_coulomb_might_be_zero_at_cutoff(ir))
{
- if (EEL_SWITCHED(ir->coulombtype))
+ if (ir_coulomb_switched(ir))
{
sprintf(err_buf,
"With coulombtype = %s rcoulomb_switch must be < rcoulomb. Or, better: Use the potential modifier options!",
}
}
+ if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT)
+ {
+ sprintf(err_buf,
+ "Explicit switch/shift coulomb interactions cannot be used in combination with a secondary coulomb-modifier.");
+ CHECK( ir->coulomb_modifier != eintmodNONE);
+ }
+ if (ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
+ {
+ sprintf(err_buf,
+ "Explicit switch/shift vdw interactions cannot be used in combination with a secondary vdw-modifier.");
+ CHECK( ir->vdw_modifier != eintmodNONE);
+ }
+
if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT ||
ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
{
warning_note(wi, warn_buf);
}
- if (ir->coulombtype == eelPMESWITCH)
+ if (ir->coulombtype == eelPMESWITCH || ir->coulomb_modifier == eintmodPOTSWITCH)
{
if (ir->rcoulomb_switch/ir->rcoulomb < 0.9499)
{
- sprintf(warn_buf, "The switching range for %s should be 5%% or less, energy conservation will be good anyhow, since ewald_rtol = %g",
- eel_names[ir->coulombtype],
- ir->ewald_rtol);
+ real percentage = 100*(ir->rcoulomb-ir->rcoulomb_switch)/ir->rcoulomb;
+ sprintf(warn_buf, "The switching range should be 5%% or less (currently %.2f%% using a switching range of %4f-%4f) for accurate electrostatic energies, energy conservation will be good regardless, since ewald_rtol = %g.",
+ percentage, ir->rcoulomb_switch, ir->rcoulomb, ir->ewald_rtol);
+ warning(wi, warn_buf);
+ }
+ }
+
+ if (ir->vdwtype == evdwSWITCH || ir->vdw_modifier == eintmodPOTSWITCH)
+ {
+ if (ir->rvdw_switch == 0)
+ {
+ sprintf(warn_buf, "rvdw-switch is equal 0 even though you are using a switched Lennard-Jones potential. This suggests it was not set in the mdp, which can lead to large energy errors. In GROMACS, 0.05 to 0.1 nm is often a reasonable vdw switching range.");
warning(wi, warn_buf);
}
}
}
}
- if (EEL_PME(ir->coulombtype))
+ if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
{
if (ir->pme_order < 3)
{
CHECK(ir->wall_ewald_zfac < 2);
}
- if (EVDW_SWITCHED(ir->vdwtype))
+ if (ir_vdw_switched(ir))
{
- sprintf(err_buf, "With vdwtype = %s rvdw-switch must be < rvdw. Or, better - use a potential modifier.",
- evdw_names[ir->vdwtype]);
+ sprintf(err_buf, "With switched vdw forces or potentials, rvdw-switch must be < rvdw");
CHECK(ir->rvdw_switch >= ir->rvdw);
+
+ if (ir->rvdw_switch < 0.5*ir->rvdw)
+ {
+ sprintf(warn_buf, "You are applying a switch function to vdw forces or potentials from %g to %g nm, which is more than half the interaction range, whereas switch functions are intended to act only close to the cut-off.",
+ ir->rvdw_switch, ir->rvdw);
+ warning_note(wi, warn_buf);
+ }
}
- else if (ir->vdwtype == evdwCUT)
+ else if (ir->vdwtype == evdwCUT || ir->vdwtype == evdwPME)
{
if (ir->cutoff_scheme == ecutsGROUP && ir->vdw_modifier == eintmodNONE)
{
CHECK(ir->rlist > ir->rvdw);
}
}
+
+ if (ir->vdwtype == evdwPME)
+ {
+ if (!(ir->vdw_modifier == eintmodNONE || ir->vdw_modifier == eintmodPOTSHIFT))
+ {
+ sprintf(err_buf, "With vdwtype = %s, the only supported modifiers are %s a\
+nd %s",
+ evdw_names[ir->vdwtype],
+ eintmod_names[eintmodPOTSHIFT],
+ eintmod_names[eintmodNONE]);
+ }
+ }
+
if (ir->cutoff_scheme == ecutsGROUP)
{
if (((ir->coulomb_modifier != eintmodNONE && ir->rcoulomb == ir->rlist) ||
"between neighborsearch steps");
}
- if (EEL_IS_ZERO_AT_CUTOFF(ir->coulombtype)
- && (ir->rlistlong <= ir->rcoulomb))
+ if (ir_coulomb_is_zero_at_cutoff(ir) && ir->rlistlong <= ir->rcoulomb)
{
sprintf(warn_buf, "For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rcoulomb.",
IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
warning_note(wi, warn_buf);
}
- if (EVDW_SWITCHED(ir->vdwtype) && (ir->rlistlong <= ir->rvdw))
+ if (ir_vdw_switched(ir) && (ir->rlistlong <= ir->rvdw))
{
sprintf(warn_buf, "For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rvdw.",
IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
/* ENERGY CONSERVATION */
if (ir_NVE(ir) && ir->cutoff_scheme == ecutsGROUP)
{
- if (!EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype) && ir->rvdw > 0 && ir->vdw_modifier == eintmodNONE)
+ if (!ir_vdw_might_be_zero_at_cutoff(ir) && ir->rvdw > 0 && ir->vdw_modifier == eintmodNONE)
{
sprintf(warn_buf, "You are using a cut-off for VdW interactions with NVE, for good energy conservation use vdwtype = %s (possibly with DispCorr)",
evdw_names[evdwSHIFT]);
warning_note(wi, warn_buf);
}
- if (!EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype) && ir->rcoulomb > 0 && ir->coulomb_modifier == eintmodNONE)
+ if (!ir_coulomb_might_be_zero_at_cutoff(ir) && ir->rcoulomb > 0)
{
sprintf(warn_buf, "You are using a cut-off for electrostatics with NVE, for good energy conservation use coulombtype = %s or %s",
eel_names[eelPMESWITCH], eel_names[eelRF_ZERO]);
}
}
+ if (EI_VV(ir->eI) && IR_TWINRANGE(*ir) && ir->nstlist > 1)
+ {
+ sprintf(warn_buf, "Twin-range multiple time stepping does not work with integrator %s.", ei_names[ir->eI]);
+ warning_error(wi, warn_buf);
+ }
+
/* IMPLICIT SOLVENT */
if (ir->coulombtype == eelGB_NOTUSED)
{
int np = 0;
char *copy0, *copy;
- copy0 = strdup(str);
+ copy0 = gmx_strdup(str);
copy = copy0;
ltrim(copy);
while (*copy != '\0')
parse_n_real(weights, &nweights, &(expand->init_lambda_weights));
if (nweights == 0)
{
- expand->bInit_weights = FALSE;
snew(expand->init_lambda_weights, fep->n_lambda); /* initialize to zero */
}
else if (nweights != fep->n_lambda)
gmx_fatal(FARGS, "Number of weights (%d) is not equal to number of lambda values (%d)",
nweights, fep->n_lambda);
}
- else
- {
- expand->bInit_weights = TRUE;
- }
if ((expand->nstexpanded < 0) && (ir->efep != efepNO))
{
expand->nstexpanded = fep->nstdhdl;
}
for (i = 0; i < ir->nwall; i++)
{
- opts->wall_atomtype[i] = strdup(names[i]);
+ opts->wall_atomtype[i] = gmx_strdup(names[i]);
}
if (ir->wall_type == ewt93 || ir->wall_type == ewt104)
t_lambda *fep = ir->fepvals;
t_expanded *expand = ir->expandedvals;
+ init_inputrec_strings();
inp = read_inpfile(mdparin, &ninp, wi);
snew(dumstr[0], STRLEN);
snew(dumstr[1], STRLEN);
- /* remove the following deprecated commands */
+ if (-1 == search_einp(ninp, inp, "cutoff-scheme"))
+ {
+ sprintf(warn_buf,
+ "%s did not specify a value for the .mdp option "
+ "\"cutoff-scheme\". Probably it was first intended for use "
+ "with GROMACS before 4.6. In 4.6, the Verlet scheme was "
+ "introduced, but the group scheme was still the default. "
+ "The default is now the Verlet scheme, so you will observe "
+ "different behaviour.", mdparin);
+ warning_note(wi, warn_buf);
+ }
+
+ /* ignore the following deprecated commands */
REM_TYPE("title");
REM_TYPE("cpp");
REM_TYPE("domain-decomposition");
REM_TYPE("dihre-tau");
REM_TYPE("nstdihreout");
REM_TYPE("nstcheckpoint");
+ REM_TYPE("optimize-fft");
/* replace the following commands with the clearer new versions*/
REPL_TYPE("unconstrained-start", "continuation");
REPL_TYPE("foreign-lambda", "fep-lambdas");
+ REPL_TYPE("verlet-buffer-drift", "verlet-buffer-tolerance");
+ REPL_TYPE("nstxtcout", "nstxout-compressed");
+ REPL_TYPE("xtc-grps", "compressed-x-grps");
+ REPL_TYPE("xtc-precision", "compressed-x-precision");
CCTYPE ("VARIOUS PREPROCESSING OPTIONS");
CTYPE ("Preprocessor information: use cpp syntax.");
CTYPE ("number of steps for center of mass motion removal");
ITYPE ("nstcomm", ir->nstcomm, 100);
CTYPE ("group(s) for center of mass motion removal");
- STYPE ("comm-grps", vcm, NULL);
+ STYPE ("comm-grps", is->vcm, NULL);
CCTYPE ("LANGEVIN DYNAMICS OPTIONS");
CTYPE ("Friction coefficient (amu/ps) and random seed");
RTYPE ("bd-fric", ir->bd_fric, 0.0);
- ITYPE ("ld-seed", ir->ld_seed, 1993);
+ STEPTYPE ("ld-seed", ir->ld_seed, -1);
/* Em stuff */
CCTYPE ("ENERGY MINIMIZATION OPTIONS");
ITYPE ("nstxout", ir->nstxout, 0);
ITYPE ("nstvout", ir->nstvout, 0);
ITYPE ("nstfout", ir->nstfout, 0);
- ir->nstcheckpoint = 1000;
CTYPE ("Output frequency for energies to log file and energy file");
ITYPE ("nstlog", ir->nstlog, 1000);
ITYPE ("nstcalcenergy", ir->nstcalcenergy, 100);
ITYPE ("nstenergy", ir->nstenergy, 1000);
CTYPE ("Output frequency and precision for .xtc file");
- ITYPE ("nstxtcout", ir->nstxtcout, 0);
- RTYPE ("xtc-precision", ir->xtcprec, 1000.0);
- CTYPE ("This selects the subset of atoms for the .xtc file. You can");
- CTYPE ("select multiple groups. By default all atoms will be written.");
- STYPE ("xtc-grps", xtc_grps, NULL);
+ ITYPE ("nstxout-compressed", ir->nstxout_compressed, 0);
+ RTYPE ("compressed-x-precision", ir->x_compression_precision, 1000.0);
+ CTYPE ("This selects the subset of atoms for the compressed");
+ CTYPE ("trajectory file. You can select multiple groups. By");
+ CTYPE ("default, all atoms will be written.");
+ STYPE ("compressed-x-grps", is->x_compressed_groups, NULL);
CTYPE ("Selection of energy groups");
- STYPE ("energygrps", energy, NULL);
+ STYPE ("energygrps", is->energy, NULL);
/* Neighbor searching */
CCTYPE ("NEIGHBORSEARCHING PARAMETERS");
- CTYPE ("cut-off scheme (group: using charge groups, Verlet: particle based cut-offs)");
+ CTYPE ("cut-off scheme (Verlet: particle based cut-offs, group: using charge groups)");
EETYPE("cutoff-scheme", ir->cutoff_scheme, ecutscheme_names);
CTYPE ("nblist update frequency");
ITYPE ("nstlist", ir->nstlist, 10);
CTYPE ("ns algorithm (simple or grid)");
EETYPE("ns-type", ir->ns_type, ens_names);
- /* set ndelta to the optimal value of 2 */
- ir->ndelta = 2;
CTYPE ("Periodic boundary conditions: xyz, no, xy");
EETYPE("pbc", ir->ePBC, epbc_names);
EETYPE("periodic-molecules", ir->bPeriodicMols, yesno_names);
- CTYPE ("Allowed energy drift due to the Verlet buffer in kJ/mol/ps per atom,");
+ CTYPE ("Allowed energy error due to the Verlet buffer in kJ/mol/ps per atom,");
CTYPE ("a value of -1 means: use rlist");
- RTYPE("verlet-buffer-drift", ir->verletbuf_drift, 0.005);
+ RTYPE("verlet-buffer-tolerance", ir->verletbuf_tol, 0.005);
CTYPE ("nblist cut-off");
RTYPE ("rlist", ir->rlist, 1.0);
CTYPE ("long-range cut-off for switched potentials");
CTYPE ("Extension of the potential lookup tables beyond the cut-off");
RTYPE ("table-extension", ir->tabext, 1.0);
CTYPE ("Separate tables between energy group pairs");
- STYPE ("energygrp-table", egptable, NULL);
+ STYPE ("energygrp-table", is->egptable, NULL);
CTYPE ("Spacing for the PME/PPPM FFT grid");
RTYPE ("fourierspacing", ir->fourier_spacing, 0.12);
CTYPE ("FFT grid size, when a value is 0 fourierspacing will be used");
CTYPE ("EWALD/PME/PPPM parameters");
ITYPE ("pme-order", ir->pme_order, 4);
RTYPE ("ewald-rtol", ir->ewald_rtol, 0.00001);
+ RTYPE ("ewald-rtol-lj", ir->ewald_rtol_lj, 0.001);
+ EETYPE("lj-pme-comb-rule", ir->ljpme_combination_rule, eljpme_names);
EETYPE("ewald-geometry", ir->ewald_geometry, eewg_names);
RTYPE ("epsilon-surface", ir->epsilon_surface, 0.0);
- EETYPE("optimize-fft", ir->bOptFFT, yesno_names);
CCTYPE("IMPLICIT SOLVENT ALGORITHM");
EETYPE("implicit-solvent", ir->implicit_solvent, eis_names);
CTYPE ("Temperature coupling");
EETYPE("tcoupl", ir->etc, etcoupl_names);
ITYPE ("nsttcouple", ir->nsttcouple, -1);
- ITYPE("nh-chain-length", ir->opts.nhchainlength, NHCHAINLENGTH);
+ ITYPE("nh-chain-length", ir->opts.nhchainlength, 10);
EETYPE("print-nose-hoover-chain-variables", ir->bPrintNHChains, yesno_names);
CTYPE ("Groups to couple separately");
- STYPE ("tc-grps", tcgrps, NULL);
+ STYPE ("tc-grps", is->tcgrps, NULL);
CTYPE ("Time constant (ps) and reference temperature (K)");
- STYPE ("tau-t", tau_t, NULL);
- STYPE ("ref-t", ref_t, NULL);
+ STYPE ("tau-t", is->tau_t, NULL);
+ STYPE ("ref-t", is->ref_t, NULL);
CTYPE ("pressure coupling");
EETYPE("pcoupl", ir->epc, epcoupl_names);
EETYPE("pcoupltype", ir->epct, epcoupltype_names);
CCTYPE ("OPTIONS FOR QMMM calculations");
EETYPE("QMMM", ir->bQMMM, yesno_names);
CTYPE ("Groups treated Quantum Mechanically");
- STYPE ("QMMM-grps", QMMM, NULL);
+ STYPE ("QMMM-grps", is->QMMM, NULL);
CTYPE ("QM method");
- STYPE("QMmethod", QMmethod, NULL);
+ STYPE("QMmethod", is->QMmethod, NULL);
CTYPE ("QMMM scheme");
EETYPE("QMMMscheme", ir->QMMMscheme, eQMMMscheme_names);
CTYPE ("QM basisset");
- STYPE("QMbasis", QMbasis, NULL);
+ STYPE("QMbasis", is->QMbasis, NULL);
CTYPE ("QM charge");
- STYPE ("QMcharge", QMcharge, NULL);
+ STYPE ("QMcharge", is->QMcharge, NULL);
CTYPE ("QM multiplicity");
- STYPE ("QMmult", QMmult, NULL);
+ STYPE ("QMmult", is->QMmult, NULL);
CTYPE ("Surface Hopping");
- STYPE ("SH", bSH, NULL);
+ STYPE ("SH", is->bSH, NULL);
CTYPE ("CAS space options");
- STYPE ("CASorbitals", CASorbitals, NULL);
- STYPE ("CASelectrons", CASelectrons, NULL);
- STYPE ("SAon", SAon, NULL);
- STYPE ("SAoff", SAoff, NULL);
- STYPE ("SAsteps", SAsteps, NULL);
+ STYPE ("CASorbitals", is->CASorbitals, NULL);
+ STYPE ("CASelectrons", is->CASelectrons, NULL);
+ STYPE ("SAon", is->SAon, NULL);
+ STYPE ("SAoff", is->SAoff, NULL);
+ STYPE ("SAsteps", is->SAsteps, NULL);
CTYPE ("Scale factor for MM charges");
RTYPE ("MMChargeScaleFactor", ir->scalefactor, 1.0);
CTYPE ("Optimization of QM subsystem");
- STYPE ("bOPT", bOPT, NULL);
- STYPE ("bTS", bTS, NULL);
+ STYPE ("bOPT", is->bOPT, NULL);
+ STYPE ("bTS", is->bTS, NULL);
/* Simulated annealing */
CCTYPE("SIMULATED ANNEALING");
CTYPE ("Type of annealing for each temperature group (no/single/periodic)");
- STYPE ("annealing", anneal, NULL);
+ STYPE ("annealing", is->anneal, NULL);
CTYPE ("Number of time points to use for specifying annealing in each group");
- STYPE ("annealing-npoints", anneal_npoints, NULL);
+ STYPE ("annealing-npoints", is->anneal_npoints, NULL);
CTYPE ("List of times at the annealing points for each group");
- STYPE ("annealing-time", anneal_time, NULL);
+ STYPE ("annealing-time", is->anneal_time, NULL);
CTYPE ("Temp. at each annealing point, for each group.");
- STYPE ("annealing-temp", anneal_temp, NULL);
+ STYPE ("annealing-temp", is->anneal_temp, NULL);
/* Startup run */
CCTYPE ("GENERATE VELOCITIES FOR STARTUP RUN");
EETYPE("gen-vel", opts->bGenVel, yesno_names);
RTYPE ("gen-temp", opts->tempi, 300.0);
- ITYPE ("gen-seed", opts->seed, 173529);
+ ITYPE ("gen-seed", opts->seed, -1);
/* Shake stuff */
CCTYPE ("OPTIONS FOR BONDS");
/* Energy group exclusions */
CCTYPE ("ENERGY GROUP EXCLUSIONS");
CTYPE ("Pairs of energy groups for which all non-bonded interactions are excluded");
- STYPE ("energygrp-excl", egpexcl, NULL);
+ STYPE ("energygrp-excl", is->egpexcl, NULL);
/* Walls */
CCTYPE ("WALLS");
ITYPE ("nwall", ir->nwall, 0);
EETYPE("wall-type", ir->wall_type, ewt_names);
RTYPE ("wall-r-linpot", ir->wall_r_linpot, -1);
- STYPE ("wall-atomtype", wall_atomtype, NULL);
- STYPE ("wall-density", wall_density, NULL);
+ STYPE ("wall-atomtype", is->wall_atomtype, NULL);
+ STYPE ("wall-density", is->wall_density, NULL);
RTYPE ("wall-ewald-zfac", ir->wall_ewald_zfac, 3);
/* COM pulling */
if (ir->ePull != epullNO)
{
snew(ir->pull, 1);
- pull_grp = read_pullparams(&ninp, &inp, ir->pull, &opts->pull_start, wi);
+ is->pull_grp = read_pullparams(&ninp, &inp, ir->pull, &opts->pull_start, wi);
}
/* Enforced rotation */
if (ir->bRot)
{
snew(ir->rot, 1);
- rot_grp = read_rotparams(&ninp, &inp, ir->rot, wi);
+ is->rot_grp = read_rotparams(&ninp, &inp, ir->rot, wi);
+ }
+
+ /* Interactive MD */
+ ir->bIMD = FALSE;
+ CCTYPE("Group to display and/or manipulate in interactive MD session");
+ STYPE ("IMD-group", is->imd_grp, NULL);
+ if (is->imd_grp[0] != '\0')
+ {
+ snew(ir->imd, 1);
+ ir->bIMD = TRUE;
}
/* Refinement */
CTYPE ("Orientation restraints force constant and tau for time averaging");
RTYPE ("orire-fc", ir->orires_fc, 0.0);
RTYPE ("orire-tau", ir->orires_tau, 0.0);
- STYPE ("orire-fitgrp", orirefitgrp, NULL);
+ STYPE ("orire-fitgrp", is->orirefitgrp, NULL);
CTYPE ("Output frequency for trace(SD) and S to energy file");
ITYPE ("nstorireout", ir->nstorireout, 100);
/* free energy variables */
CCTYPE ("Free energy variables");
EETYPE("free-energy", ir->efep, efep_names);
- STYPE ("couple-moltype", couple_moltype, NULL);
+ STYPE ("couple-moltype", is->couple_moltype, NULL);
EETYPE("couple-lambda0", opts->couple_lam0, couple_lam);
EETYPE("couple-lambda1", opts->couple_lam1, couple_lam);
EETYPE("couple-intramol", opts->bCoupleIntra, yesno_names);
ITYPE ("init-lambda-state", fep->init_fep_state, -1);
RTYPE ("delta-lambda", fep->delta_lambda, 0.0);
ITYPE ("nstdhdl", fep->nstdhdl, 50);
- STYPE ("fep-lambdas", fep_lambda[efptFEP], NULL);
- STYPE ("mass-lambdas", fep_lambda[efptMASS], NULL);
- STYPE ("coul-lambdas", fep_lambda[efptCOUL], NULL);
- STYPE ("vdw-lambdas", fep_lambda[efptVDW], NULL);
- STYPE ("bonded-lambdas", fep_lambda[efptBONDED], NULL);
- STYPE ("restraint-lambdas", fep_lambda[efptRESTRAINT], NULL);
- STYPE ("temperature-lambdas", fep_lambda[efptTEMPERATURE], NULL);
+ STYPE ("fep-lambdas", is->fep_lambda[efptFEP], NULL);
+ STYPE ("mass-lambdas", is->fep_lambda[efptMASS], NULL);
+ STYPE ("coul-lambdas", is->fep_lambda[efptCOUL], NULL);
+ STYPE ("vdw-lambdas", is->fep_lambda[efptVDW], NULL);
+ STYPE ("bonded-lambdas", is->fep_lambda[efptBONDED], NULL);
+ STYPE ("restraint-lambdas", is->fep_lambda[efptRESTRAINT], NULL);
+ STYPE ("temperature-lambdas", is->fep_lambda[efptTEMPERATURE], NULL);
ITYPE ("calc-lambda-neighbors", fep->lambda_neighbors, 1);
- STYPE ("init-lambda-weights", lambda_weights, NULL);
- EETYPE("dhdl-print-energy", fep->bPrintEnergy, yesno_names);
+ STYPE ("init-lambda-weights", is->lambda_weights, NULL);
+ EETYPE("dhdl-print-energy", fep->edHdLPrintEnergy, edHdLPrintEnergy_names);
RTYPE ("sc-alpha", fep->sc_alpha, 0.0);
ITYPE ("sc-power", fep->sc_power, 1);
RTYPE ("sc-r-power", fep->sc_r_power, 6.0);
/* Non-equilibrium MD stuff */
CCTYPE("Non-equilibrium MD stuff");
- STYPE ("acc-grps", accgrps, NULL);
- STYPE ("accelerate", acc, NULL);
- STYPE ("freezegrps", freeze, NULL);
- STYPE ("freezedim", frdim, NULL);
+ STYPE ("acc-grps", is->accgrps, NULL);
+ STYPE ("accelerate", is->acc, NULL);
+ STYPE ("freezegrps", is->freeze, NULL);
+ STYPE ("freezedim", is->frdim, NULL);
RTYPE ("cos-acceleration", ir->cos_accel, 0);
- STYPE ("deform", deform, NULL);
+ STYPE ("deform", is->deform, NULL);
/* simulated tempering variables */
CCTYPE("simulated tempering variables");
CCTYPE("Electric fields");
CTYPE ("Format is number of terms (int) and for all terms an amplitude (real)");
CTYPE ("and a phase angle (real)");
- STYPE ("E-x", efield_x, NULL);
- STYPE ("E-xt", efield_xt, NULL);
- STYPE ("E-y", efield_y, NULL);
- STYPE ("E-yt", efield_yt, NULL);
- STYPE ("E-z", efield_z, NULL);
- STYPE ("E-zt", efield_zt, NULL);
+ STYPE ("E-x", is->efield_x, NULL);
+ STYPE ("E-xt", is->efield_xt, NULL);
+ STYPE ("E-y", is->efield_y, NULL);
+ STYPE ("E-yt", is->efield_yt, NULL);
+ STYPE ("E-z", is->efield_z, NULL);
+ STYPE ("E-zt", is->efield_zt, NULL);
+
+ CCTYPE("Ion/water position swapping for computational electrophysiology setups");
+ CTYPE("Swap positions along direction: no, X, Y, Z");
+ EETYPE("swapcoords", ir->eSwapCoords, eSwapTypes_names);
+ if (ir->eSwapCoords != eswapNO)
+ {
+ snew(ir->swap, 1);
+ CTYPE("Swap attempt frequency");
+ ITYPE("swap-frequency", ir->swap->nstswap, 1);
+ CTYPE("Two index groups that contain the compartment-partitioning atoms");
+ STYPE("split-group0", splitgrp0, NULL);
+ STYPE("split-group1", splitgrp1, NULL);
+ CTYPE("Use center of mass of split groups (yes/no), otherwise center of geometry is used");
+ EETYPE("massw-split0", ir->swap->massw_split[0], yesno_names);
+ EETYPE("massw-split1", ir->swap->massw_split[1], yesno_names);
+
+ CTYPE("Group name of ions that can be exchanged with solvent molecules");
+ STYPE("swap-group", swapgrp, NULL);
+ CTYPE("Group name of solvent molecules");
+ STYPE("solvent-group", solgrp, NULL);
+
+ CTYPE("Split cylinder: radius, upper and lower extension (nm) (this will define the channels)");
+ CTYPE("Note that the split cylinder settings do not have an influence on the swapping protocol,");
+ CTYPE("however, if correctly defined, the ion permeation events are counted per channel");
+ RTYPE("cyl0-r", ir->swap->cyl0r, 2.0);
+ RTYPE("cyl0-up", ir->swap->cyl0u, 1.0);
+ RTYPE("cyl0-down", ir->swap->cyl0l, 1.0);
+ RTYPE("cyl1-r", ir->swap->cyl1r, 2.0);
+ RTYPE("cyl1-up", ir->swap->cyl1u, 1.0);
+ RTYPE("cyl1-down", ir->swap->cyl1l, 1.0);
+
+ CTYPE("Average the number of ions per compartment over these many swap attempt steps");
+ ITYPE("coupl-steps", ir->swap->nAverage, 10);
+ CTYPE("Requested number of anions and cations for each of the two compartments");
+ CTYPE("-1 means fix the numbers as found in time step 0");
+ ITYPE("anionsA", ir->swap->nanions[0], -1);
+ ITYPE("cationsA", ir->swap->ncations[0], -1);
+ ITYPE("anionsB", ir->swap->nanions[1], -1);
+ ITYPE("cationsB", ir->swap->ncations[1], -1);
+ CTYPE("Start to swap ions if threshold difference to requested count is reached");
+ RTYPE("threshold", ir->swap->threshold, 1.0);
+ }
/* AdResS defined thingies */
CCTYPE ("AdResS parameters");
/* User defined thingies */
CCTYPE ("User defined thingies");
- STYPE ("user1-grps", user1, NULL);
- STYPE ("user2-grps", user2, NULL);
+ STYPE ("user1-grps", is->user1, NULL);
+ STYPE ("user2-grps", is->user2, NULL);
ITYPE ("userint1", ir->userint1, 0);
ITYPE ("userint2", ir->userint2, 0);
ITYPE ("userint3", ir->userint3, 0);
}
opts->couple_moltype = NULL;
- if (strlen(couple_moltype) > 0)
+ if (strlen(is->couple_moltype) > 0)
{
if (ir->efep != efepNO)
{
- opts->couple_moltype = strdup(couple_moltype);
+ opts->couple_moltype = gmx_strdup(is->couple_moltype);
if (opts->couple_lam0 == opts->couple_lam1)
{
warning(wi, "The lambda=0 and lambda=1 states for coupling are identical");
}
else
{
- warning(wi, "Can not couple a molecule with free_energy = no");
+ warning_note(wi, "Free energy is turned off, so we will not decouple the molecule listed in your input.");
}
}
/* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
}
}
- if (ir->bSimTemp)
+ if (fep->edHdLPrintEnergy == edHdLPrintEnergyYES)
{
- fep->bPrintEnergy = TRUE;
- /* always print out the energy to dhdl if we are doing expanded ensemble, since we need the total energy
- if the temperature is changing. */
+ fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
+ warning_note(wi, "Old option for dhdl-print-energy given: "
+ "changing \"yes\" to \"total\"\n");
+ }
+
+ if (ir->bSimTemp && (fep->edHdLPrintEnergy == edHdLPrintEnergyNO))
+ {
+ /* always print out the energy to dhdl if we are doing
+ expanded ensemble, since we need the total energy for
+ analysis if the temperature is changing. In some
+ conditions one may only want the potential energy, so
+ we will allow that if the appropriate mdp setting has
+ been enabled. Otherwise, total it is:
+ */
+ fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
}
if ((ir->efep != efepNO) || ir->bSimTemp)
{
ir->bExpanded = TRUE;
}
- do_fep_params(ir, fep_lambda, lambda_weights);
+ do_fep_params(ir, is->fep_lambda, is->lambda_weights);
if (ir->bSimTemp) /* done after fep params */
{
do_simtemp_params(ir);
/* WALL PARAMETERS */
- do_wall_params(ir, wall_atomtype, wall_density, opts);
+ do_wall_params(ir, is->wall_atomtype, is->wall_density, opts);
/* ORIENTATION RESTRAINT PARAMETERS */
- if (opts->bOrire && str_nelem(orirefitgrp, MAXPTR, NULL) != 1)
+ if (opts->bOrire && str_nelem(is->orirefitgrp, MAXPTR, NULL) != 1)
{
warning_error(wi, "ERROR: Need one orientation restraint fit group\n");
}
{
dumdub[0][i] = 0;
}
- m = sscanf(deform, "%lf %lf %lf %lf %lf %lf",
+ m = sscanf(is->deform, "%lf %lf %lf %lf %lf %lf",
&(dumdub[0][0]), &(dumdub[0][1]), &(dumdub[0][2]),
&(dumdub[0][3]), &(dumdub[0][4]), &(dumdub[0][5]));
for (i = 0; i < 3; i++)
}
}
+ /* Ion/water position swapping checks */
+ if (ir->eSwapCoords != eswapNO)
+ {
+ if (ir->swap->nstswap < 1)
+ {
+ warning_error(wi, "swap_frequency must be 1 or larger when ion swapping is requested");
+ }
+ if (ir->swap->nAverage < 1)
+ {
+ warning_error(wi, "coupl_steps must be 1 or larger.\n");
+ }
+ if (ir->swap->threshold < 1.0)
+ {
+ warning_error(wi, "Ion count threshold must be at least 1.\n");
+ }
+ }
+
sfree(dumstr[0]);
sfree(dumstr[1]);
}
-static int search_QMstring(char *s, int ng, const char *gn[])
+static int search_QMstring(const char *s, int ng, const char *gn[])
{
/* same as normal search_string, but this one searches QM strings */
int i;
} /* search_QMstring */
-
-int search_string(char *s, int ng, char *gn[])
+/* We would like gn to be const as well, but C doesn't allow this */
+int search_string(const char *s, int ng, char *gn[])
{
int i;
t_grpopts *opts;
gmx_groups_t *groups;
t_pull *pull;
- int natoms, ai, aj, i, j, d, g, imin, jmin, nc;
+ int natoms, ai, aj, i, j, d, g, imin, jmin;
t_iatom *ia;
int *nrdf2, *na_vcm, na_tot;
double *nrdf_tc, *nrdf_vcm, nrdf_uc, n_sub = 0;
* to determine the optimal nrdf assignment.
*/
pull = ir->pull;
- if (pull->eGeom == epullgPOS)
+
+ for (i = 0; i < pull->ncoord; i++)
{
- nc = 0;
- for (i = 0; i < DIM; i++)
+ imin = 1;
+
+ for (j = 0; j < 2; j++)
{
- if (pull->dim[i])
+ const t_pull_group *pgrp;
+
+ pgrp = &pull->group[pull->coord[i].group[j]];
+
+ if (pgrp->nat > 0)
{
- nc++;
+ /* Subtract 1/2 dof from each group */
+ ai = pgrp->ind[0];
+ nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
+ nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
+ if (nrdf_tc[ggrpnr(groups, egcTC, ai)] < 0)
+ {
+ gmx_fatal(FARGS, "Center of mass pulling constraints caused the number of degrees of freedom for temperature coupling group %s to be negative", gnames[groups->grps[egcTC].nm_ind[ggrpnr(groups, egcTC, ai)]]);
+ }
}
- }
- }
- else
- {
- nc = 1;
- }
- for (i = 0; i < pull->ngrp; i++)
- {
- imin = 2*nc;
- if (pull->grp[0].nat > 0)
- {
- /* Subtract 1/2 dof from the reference group */
- ai = pull->grp[0].ind[0];
- if (nrdf_tc[ggrpnr(groups, egcTC, ai)] > 1)
+ else
{
- nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5;
- nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5;
- imin--;
+ /* We need to subtract the whole DOF from group j=1 */
+ imin += 1;
}
}
- /* Subtract 1/2 dof from the pulled group */
- ai = pull->grp[1+i].ind[0];
- nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
- nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
- if (nrdf_tc[ggrpnr(groups, egcTC, ai)] < 0)
- {
- gmx_fatal(FARGS, "Center of mass pulling constraints caused the number of degrees of freedom for temperature coupling group %s to be negative", gnames[groups->grps[egcTC].nm_ind[ggrpnr(groups, egcTC, ai)]]);
- }
}
}
double a, phi;
int i;
- t = strdup(s);
+ t = gmx_strdup(s);
trim(t);
cosine->n = 0;
return bSet;
}
+
+static void make_swap_groups(
+ t_swapcoords *swap,
+ char *swapgname,
+ char *splitg0name,
+ char *splitg1name,
+ char *solgname,
+ t_blocka *grps,
+ char **gnames)
+{
+ int ig = -1, i = 0, j;
+ char *splitg;
+
+
+ /* Just a quick check here, more thorough checks are in mdrun */
+ if (strcmp(splitg0name, splitg1name) == 0)
+ {
+ gmx_fatal(FARGS, "The split groups can not both be '%s'.", splitg0name);
+ }
+
+ /* First get the swap group index atoms */
+ ig = search_string(swapgname, grps->nr, gnames);
+ swap->nat = grps->index[ig+1] - grps->index[ig];
+ if (swap->nat > 0)
+ {
+ fprintf(stderr, "Swap group '%s' contains %d atoms.\n", swapgname, swap->nat);
+ snew(swap->ind, swap->nat);
+ for (i = 0; i < swap->nat; i++)
+ {
+ swap->ind[i] = grps->a[grps->index[ig]+i];
+ }
+ }
+ else
+ {
+ gmx_fatal(FARGS, "You defined an empty group of atoms for swapping.");
+ }
+
+ /* Now do so for the split groups */
+ for (j = 0; j < 2; j++)
+ {
+ if (j == 0)
+ {
+ splitg = splitg0name;
+ }
+ else
+ {
+ splitg = splitg1name;
+ }
+
+ ig = search_string(splitg, grps->nr, gnames);
+ swap->nat_split[j] = grps->index[ig+1] - grps->index[ig];
+ if (swap->nat_split[j] > 0)
+ {
+ fprintf(stderr, "Split group %d '%s' contains %d atom%s.\n",
+ j, splitg, swap->nat_split[j], (swap->nat_split[j] > 1) ? "s" : "");
+ snew(swap->ind_split[j], swap->nat_split[j]);
+ for (i = 0; i < swap->nat_split[j]; i++)
+ {
+ swap->ind_split[j][i] = grps->a[grps->index[ig]+i];
+ }
+ }
+ else
+ {
+ gmx_fatal(FARGS, "Split group %d has to contain at least 1 atom!", j);
+ }
+ }
+
+ /* Now get the solvent group index atoms */
+ ig = search_string(solgname, grps->nr, gnames);
+ swap->nat_sol = grps->index[ig+1] - grps->index[ig];
+ if (swap->nat_sol > 0)
+ {
+ fprintf(stderr, "Solvent group '%s' contains %d atoms.\n", solgname, swap->nat_sol);
+ snew(swap->ind_sol, swap->nat_sol);
+ for (i = 0; i < swap->nat_sol; i++)
+ {
+ swap->ind_sol[i] = grps->a[grps->index[ig]+i];
+ }
+ }
+ else
+ {
+ gmx_fatal(FARGS, "You defined an empty group of solvent. Cannot exchange ions.");
+ }
+}
+
+
+void make_IMD_group(t_IMD *IMDgroup, char *IMDgname, t_blocka *grps, char **gnames)
+{
+ int ig, i;
+
+
+ ig = search_string(IMDgname, grps->nr, gnames);
+ IMDgroup->nat = grps->index[ig+1] - grps->index[ig];
+
+ if (IMDgroup->nat > 0)
+ {
+ fprintf(stderr, "Group '%s' with %d atoms can be activated for interactive molecular dynamics (IMD).\n",
+ IMDgname, IMDgroup->nat);
+ snew(IMDgroup->ind, IMDgroup->nat);
+ for (i = 0; i < IMDgroup->nat; i++)
+ {
+ IMDgroup->ind[i] = grps->a[grps->index[ig]+i];
+ }
+ }
+}
+
+
void do_index(const char* mdparin, const char *ndx,
gmx_mtop_t *mtop,
gmx_bool bVerbose,
set_warning_line(wi, mdparin, -1);
- ntau_t = str_nelem(tau_t, MAXPTR, ptr1);
- nref_t = str_nelem(ref_t, MAXPTR, ptr2);
- ntcg = str_nelem(tcgrps, MAXPTR, ptr3);
+ ntau_t = str_nelem(is->tau_t, MAXPTR, ptr1);
+ nref_t = str_nelem(is->ref_t, MAXPTR, ptr2);
+ ntcg = str_nelem(is->tcgrps, MAXPTR, ptr3);
if ((ntau_t != ntcg) || (nref_t != ntcg))
{
gmx_fatal(FARGS, "Invalid T coupling input: %d groups, %d ref-t values and "
nstcmin = tcouple_min_integration_steps(ir->etc);
if (nstcmin > 1)
{
- if (tau_min/(ir->delta_t*ir->nsttcouple) < nstcmin)
+ if (tau_min/(ir->delta_t*ir->nsttcouple) < nstcmin - 10*GMX_REAL_EPS)
{
sprintf(warn_buf, "For proper integration of the %s thermostat, tau-t (%g) should be at least %d times larger than nsttcouple*dt (%g)",
ETCOUPLTYPE(ir->etc),
}
/* Simulated annealing for each group. There are nr groups */
- nSA = str_nelem(anneal, MAXPTR, ptr1);
+ nSA = str_nelem(is->anneal, MAXPTR, ptr1);
if (nSA == 1 && (ptr1[0][0] == 'n' || ptr1[0][0] == 'N'))
{
nSA = 0;
if (bAnneal)
{
/* Read the other fields too */
- nSA_points = str_nelem(anneal_npoints, MAXPTR, ptr1);
+ nSA_points = str_nelem(is->anneal_npoints, MAXPTR, ptr1);
if (nSA_points != nSA)
{
gmx_fatal(FARGS, "Found %d annealing-npoints values for %d groups\n", nSA_points, nSA);
k += ir->opts.anneal_npoints[i];
}
- nSA_time = str_nelem(anneal_time, MAXPTR, ptr1);
+ nSA_time = str_nelem(is->anneal_time, MAXPTR, ptr1);
if (nSA_time != k)
{
gmx_fatal(FARGS, "Found %d annealing-time values, wanter %d\n", nSA_time, k);
}
- nSA_temp = str_nelem(anneal_temp, MAXPTR, ptr2);
+ nSA_temp = str_nelem(is->anneal_temp, MAXPTR, ptr2);
if (nSA_temp != k)
{
gmx_fatal(FARGS, "Found %d annealing-temp values, wanted %d\n", nSA_temp, k);
if (ir->ePull != epullNO)
{
- make_pull_groups(ir->pull, pull_grp, grps, gnames);
+ make_pull_groups(ir->pull, is->pull_grp, grps, gnames);
+
+ make_pull_coords(ir->pull);
}
if (ir->bRot)
{
- make_rotation_groups(ir->rot, rot_grp, grps, gnames);
+ make_rotation_groups(ir->rot, is->rot_grp, grps, gnames);
+ }
+
+ if (ir->eSwapCoords != eswapNO)
+ {
+ make_swap_groups(ir->swap, swapgrp, splitgrp0, splitgrp1, solgrp, grps, gnames);
+ }
+
+ /* Make indices for IMD session */
+ if (ir->bIMD)
+ {
+ make_IMD_group(ir->imd, is->imd_grp, grps, gnames);
}
- nacc = str_nelem(acc, MAXPTR, ptr1);
- nacg = str_nelem(accgrps, MAXPTR, ptr2);
+ nacc = str_nelem(is->acc, MAXPTR, ptr1);
+ nacg = str_nelem(is->accgrps, MAXPTR, ptr2);
if (nacg*DIM != nacc)
{
gmx_fatal(FARGS, "Invalid Acceleration input: %d groups and %d acc. values",
}
}
- nfrdim = str_nelem(frdim, MAXPTR, ptr1);
- nfreeze = str_nelem(freeze, MAXPTR, ptr2);
+ nfrdim = str_nelem(is->frdim, MAXPTR, ptr1);
+ nfreeze = str_nelem(is->freeze, MAXPTR, ptr2);
if (nfrdim != DIM*nfreeze)
{
gmx_fatal(FARGS, "Invalid Freezing input: %d groups and %d freeze values",
}
}
- nenergy = str_nelem(energy, MAXPTR, ptr1);
+ nenergy = str_nelem(is->energy, MAXPTR, ptr1);
do_numbering(natoms, groups, nenergy, ptr1, grps, gnames, egcENER,
restnm, egrptpALL_GENREST, bVerbose, wi);
add_wall_energrps(groups, ir->nwall, symtab);
ir->opts.ngener = groups->grps[egcENER].nr;
- nvcm = str_nelem(vcm, MAXPTR, ptr1);
+ nvcm = str_nelem(is->vcm, MAXPTR, ptr1);
bRest =
do_numbering(natoms, groups, nvcm, ptr1, grps, gnames, egcVCM,
restnm, nvcm == 0 ? egrptpALL_GENREST : egrptpPART, bVerbose, wi);
}
}
- nuser = str_nelem(user1, MAXPTR, ptr1);
+ nuser = str_nelem(is->user1, MAXPTR, ptr1);
do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser1,
restnm, egrptpALL_GENREST, bVerbose, wi);
- nuser = str_nelem(user2, MAXPTR, ptr1);
+ nuser = str_nelem(is->user2, MAXPTR, ptr1);
do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser2,
restnm, egrptpALL_GENREST, bVerbose, wi);
- nuser = str_nelem(xtc_grps, MAXPTR, ptr1);
- do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcXTC,
+ nuser = str_nelem(is->x_compressed_groups, MAXPTR, ptr1);
+ do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcCompressedX,
restnm, egrptpONE, bVerbose, wi);
- nofg = str_nelem(orirefitgrp, MAXPTR, ptr1);
+ nofg = str_nelem(is->orirefitgrp, MAXPTR, ptr1);
do_numbering(natoms, groups, nofg, ptr1, grps, gnames, egcORFIT,
restnm, egrptpALL_GENREST, bVerbose, wi);
/* QMMM input processing */
- nQMg = str_nelem(QMMM, MAXPTR, ptr1);
- nQMmethod = str_nelem(QMmethod, MAXPTR, ptr2);
- nQMbasis = str_nelem(QMbasis, MAXPTR, ptr3);
+ nQMg = str_nelem(is->QMMM, MAXPTR, ptr1);
+ nQMmethod = str_nelem(is->QMmethod, MAXPTR, ptr2);
+ nQMbasis = str_nelem(is->QMbasis, MAXPTR, ptr3);
if ((nQMmethod != nQMg) || (nQMbasis != nQMg))
{
gmx_fatal(FARGS, "Invalid QMMM input: %d groups %d basissets"
eQMbasis_names);
}
- nQMmult = str_nelem(QMmult, MAXPTR, ptr1);
- nQMcharge = str_nelem(QMcharge, MAXPTR, ptr2);
- nbSH = str_nelem(bSH, MAXPTR, ptr3);
+ nQMmult = str_nelem(is->QMmult, MAXPTR, ptr1);
+ nQMcharge = str_nelem(is->QMcharge, MAXPTR, ptr2);
+ nbSH = str_nelem(is->bSH, MAXPTR, ptr3);
snew(ir->opts.QMmult, nr);
snew(ir->opts.QMcharge, nr);
snew(ir->opts.bSH, nr);
ir->opts.bSH[i] = (gmx_strncasecmp(ptr3[i], "Y", 1) == 0);
}
- nCASelec = str_nelem(CASelectrons, MAXPTR, ptr1);
- nCASorb = str_nelem(CASorbitals, MAXPTR, ptr2);
+ nCASelec = str_nelem(is->CASelectrons, MAXPTR, ptr1);
+ nCASorb = str_nelem(is->CASorbitals, MAXPTR, ptr2);
snew(ir->opts.CASelectrons, nr);
snew(ir->opts.CASorbitals, nr);
for (i = 0; i < nr; i++)
}
/* special optimization options */
- nbOPT = str_nelem(bOPT, MAXPTR, ptr1);
- nbTS = str_nelem(bTS, MAXPTR, ptr2);
+ nbOPT = str_nelem(is->bOPT, MAXPTR, ptr1);
+ nbTS = str_nelem(is->bTS, MAXPTR, ptr2);
snew(ir->opts.bOPT, nr);
snew(ir->opts.bTS, nr);
for (i = 0; i < nr; i++)
ir->opts.bOPT[i] = (gmx_strncasecmp(ptr1[i], "Y", 1) == 0);
ir->opts.bTS[i] = (gmx_strncasecmp(ptr2[i], "Y", 1) == 0);
}
- nSAon = str_nelem(SAon, MAXPTR, ptr1);
- nSAoff = str_nelem(SAoff, MAXPTR, ptr2);
- nSAsteps = str_nelem(SAsteps, MAXPTR, ptr3);
+ nSAon = str_nelem(is->SAon, MAXPTR, ptr1);
+ nSAoff = str_nelem(is->SAoff, MAXPTR, ptr2);
+ nSAsteps = str_nelem(is->SAsteps, MAXPTR, ptr3);
snew(ir->opts.SAon, nr);
snew(ir->opts.SAoff, nr);
snew(ir->opts.SAsteps, nr);
nr = groups->grps[egcENER].nr;
snew(ir->opts.egp_flags, nr*nr);
- bExcl = do_egp_flag(ir, groups, "energygrp-excl", egpexcl, EGP_EXCL);
+ bExcl = do_egp_flag(ir, groups, "energygrp-excl", is->egpexcl, EGP_EXCL);
if (bExcl && ir->cutoff_scheme == ecutsVERLET)
{
warning_error(wi, "Energy group exclusions are not (yet) implemented for the Verlet scheme");
warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
}
- bTable = do_egp_flag(ir, groups, "energygrp-table", egptable, EGP_TABLE);
+ bTable = do_egp_flag(ir, groups, "energygrp-table", is->egptable, EGP_TABLE);
if (bTable && !(ir->vdwtype == evdwUSER) &&
!(ir->coulombtype == eelUSER) && !(ir->coulombtype == eelPMEUSER) &&
!(ir->coulombtype == eelPMEUSERSWITCH))
gmx_fatal(FARGS, "Can only have energy group pair tables in combination with user tables for VdW and/or Coulomb");
}
- decode_cos(efield_x, &(ir->ex[XX]));
- decode_cos(efield_xt, &(ir->et[XX]));
- decode_cos(efield_y, &(ir->ex[YY]));
- decode_cos(efield_yt, &(ir->et[YY]));
- decode_cos(efield_z, &(ir->ex[ZZ]));
- decode_cos(efield_zt, &(ir->et[ZZ]));
+ decode_cos(is->efield_x, &(ir->ex[XX]));
+ decode_cos(is->efield_xt, &(ir->et[XX]));
+ decode_cos(is->efield_y, &(ir->ex[YY]));
+ decode_cos(is->efield_yt, &(ir->et[YY]));
+ decode_cos(is->efield_z, &(ir->ex[ZZ]));
+ decode_cos(is->efield_zt, &(ir->et[ZZ]));
if (ir->bAdress)
{
return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
}
+static void
+check_combination_rule_differences(const gmx_mtop_t *mtop, int state,
+ gmx_bool *bC6ParametersWorkWithGeometricRules,
+ gmx_bool *bC6ParametersWorkWithLBRules,
+ gmx_bool *bLBRulesPossible)
+{
+ int ntypes, tpi, tpj, thisLBdiff, thisgeomdiff;
+ int *typecount;
+ real tol;
+ double geometricdiff, LBdiff;
+ double c6i, c6j, c12i, c12j;
+ double c6, c6_geometric, c6_LB;
+ double sigmai, sigmaj, epsi, epsj;
+ gmx_bool bCanDoLBRules, bCanDoGeometricRules;
+ const char *ptr;
+
+ /* A tolerance of 1e-5 seems reasonable for (possibly hand-typed)
+ * force-field floating point parameters.
+ */
+ tol = 1e-5;
+ ptr = getenv("GMX_LJCOMB_TOL");
+ if (ptr != NULL)
+ {
+ double dbl;
+
+ sscanf(ptr, "%lf", &dbl);
+ tol = dbl;
+ }
+
+ *bC6ParametersWorkWithLBRules = TRUE;
+ *bC6ParametersWorkWithGeometricRules = TRUE;
+ bCanDoLBRules = TRUE;
+ bCanDoGeometricRules = TRUE;
+ ntypes = mtop->ffparams.atnr;
+ snew(typecount, ntypes);
+ gmx_mtop_count_atomtypes(mtop, state, typecount);
+ geometricdiff = LBdiff = 0.0;
+ *bLBRulesPossible = TRUE;
+ for (tpi = 0; tpi < ntypes; ++tpi)
+ {
+ c6i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c6;
+ c12i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c12;
+ for (tpj = tpi; tpj < ntypes; ++tpj)
+ {
+ c6j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c6;
+ c12j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c12;
+ c6 = mtop->ffparams.iparams[ntypes * tpi + tpj].lj.c6;
+ c6_geometric = sqrt(c6i * c6j);
+ if (!gmx_numzero(c6_geometric))
+ {
+ if (!gmx_numzero(c12i) && !gmx_numzero(c12j))
+ {
+ sigmai = pow(c12i / c6i, 1.0/6.0);
+ sigmaj = pow(c12j / c6j, 1.0/6.0);
+ epsi = c6i * c6i /(4.0 * c12i);
+ epsj = c6j * c6j /(4.0 * c12j);
+ c6_LB = 4.0 * pow(epsi * epsj, 1.0/2.0) * pow(0.5 * (sigmai + sigmaj), 6);
+ }
+ else
+ {
+ *bLBRulesPossible = FALSE;
+ c6_LB = c6_geometric;
+ }
+ bCanDoLBRules = gmx_within_tol(c6_LB, c6, tol);
+ }
+
+ if (FALSE == bCanDoLBRules)
+ {
+ *bC6ParametersWorkWithLBRules = FALSE;
+ }
+
+ bCanDoGeometricRules = gmx_within_tol(c6_geometric, c6, tol);
+
+ if (FALSE == bCanDoGeometricRules)
+ {
+ *bC6ParametersWorkWithGeometricRules = FALSE;
+ }
+ }
+ }
+ sfree(typecount);
+}
+
+static void
+check_combination_rules(const t_inputrec *ir, const gmx_mtop_t *mtop,
+ warninp_t wi)
+{
+ char err_buf[256];
+ gmx_bool bLBRulesPossible, bC6ParametersWorkWithGeometricRules, bC6ParametersWorkWithLBRules;
+
+ check_combination_rule_differences(mtop, 0,
+ &bC6ParametersWorkWithGeometricRules,
+ &bC6ParametersWorkWithLBRules,
+ &bLBRulesPossible);
+ if (ir->ljpme_combination_rule == eljpmeLB)
+ {
+ if (FALSE == bC6ParametersWorkWithLBRules || FALSE == bLBRulesPossible)
+ {
+ warning(wi, "You are using arithmetic-geometric combination rules "
+ "in LJ-PME, but your non-bonded C6 parameters do not "
+ "follow these rules.");
+ }
+ }
+ else
+ {
+ if (FALSE == bC6ParametersWorkWithGeometricRules)
+ {
+ if (ir->eDispCorr != edispcNO)
+ {
+ warning_note(wi, "You are using geometric combination rules in "
+ "LJ-PME, but your non-bonded C6 parameters do "
+ "not follow these rules. "
+ "This will introduce very small errors in the forces and energies in "
+ "your simulations. Dispersion correction will correct total energy "
+ "and/or pressure for isotropic systems, but not forces or surface tensions.");
+ }
+ else
+ {
+ warning_note(wi, "You are using geometric combination rules in "
+ "LJ-PME, but your non-bonded C6 parameters do "
+ "not follow these rules. "
+ "This will introduce very small errors in the forces and energies in "
+ "your simulations. If your system is homogeneous, consider using dispersion correction "
+ "for the total energy and pressure.");
+ }
+ }
+ }
+}
+
void triple_check(const char *mdparin, t_inputrec *ir, gmx_mtop_t *sys,
warninp_t wi)
{
- char err_buf[256];
- int i, m, g, nmol, npct;
+ char err_buf[STRLEN];
+ int i, m, c, nmol, npct;
gmx_bool bCharge, bAcc;
real gdt_max, *mgrp, mt;
rvec acc;
set_warning_line(wi, mdparin, -1);
+ if (ir->cutoff_scheme == ecutsVERLET &&
+ ir->verletbuf_tol > 0 &&
+ ir->nstlist > 1 &&
+ ((EI_MD(ir->eI) || EI_SD(ir->eI)) &&
+ (ir->etc == etcVRESCALE || ir->etc == etcBERENDSEN)))
+ {
+ /* Check if a too small Verlet buffer might potentially
+ * cause more drift than the thermostat can couple off.
+ */
+ /* Temperature error fraction for warning and suggestion */
+ const real T_error_warn = 0.002;
+ const real T_error_suggest = 0.001;
+ /* For safety: 2 DOF per atom (typical with constraints) */
+ const real nrdf_at = 2;
+ real T, tau, max_T_error;
+ int i;
+
+ T = 0;
+ tau = 0;
+ for (i = 0; i < ir->opts.ngtc; i++)
+ {
+ T = max(T, ir->opts.ref_t[i]);
+ tau = max(tau, ir->opts.tau_t[i]);
+ }
+ if (T > 0)
+ {
+ /* This is a worst case estimate of the temperature error,
+ * assuming perfect buffer estimation and no cancelation
+ * of errors. The factor 0.5 is because energy distributes
+ * equally over Ekin and Epot.
+ */
+ max_T_error = 0.5*tau*ir->verletbuf_tol/(nrdf_at*BOLTZ*T);
+ if (max_T_error > T_error_warn)
+ {
+ sprintf(warn_buf, "With a verlet-buffer-tolerance of %g kJ/mol/ps, a reference temperature of %g and a tau_t of %g, your temperature might be off by up to %.1f%%. To ensure the error is below %.1f%%, decrease verlet-buffer-tolerance to %.0e or decrease tau_t.",
+ ir->verletbuf_tol, T, tau,
+ 100*max_T_error,
+ 100*T_error_suggest,
+ ir->verletbuf_tol*T_error_suggest/max_T_error);
+ warning(wi, warn_buf);
+ }
+ }
+ }
+
+ if (ETC_ANDERSEN(ir->etc))
+ {
+ int i;
+
+ for (i = 0; i < ir->opts.ngtc; i++)
+ {
+ sprintf(err_buf, "all tau_t must currently be equal using Andersen temperature control, violated for group %d", i);
+ CHECK(ir->opts.tau_t[0] != ir->opts.tau_t[i]);
+ sprintf(err_buf, "all tau_t must be postive using Andersen temperature control, tau_t[%d]=%10.6f",
+ i, ir->opts.tau_t[i]);
+ CHECK(ir->opts.tau_t[i] < 0);
+ }
+
+ for (i = 0; i < ir->opts.ngtc; i++)
+ {
+ int nsteps = (int)(ir->opts.tau_t[i]/ir->delta_t);
+ sprintf(err_buf, "tau_t/delta_t for group %d for temperature control method %s must be a multiple of nstcomm (%d), as velocities of atoms in coupled groups are randomized every time step. The input tau_t (%8.3f) leads to %d steps per randomization", i, etcoupl_names[ir->etc], ir->nstcomm, ir->opts.tau_t[i], nsteps);
+ CHECK((nsteps % ir->nstcomm) && (ir->etc == etcANDERSENMASSIVE));
+ }
+ }
+
if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD &&
ir->comm_mode == ecmNO &&
- !(absolute_reference(ir, sys, FALSE, AbsRef) || ir->nsteps <= 10))
+ !(absolute_reference(ir, sys, FALSE, AbsRef) || ir->nsteps <= 10) &&
+ !ETC_ANDERSEN(ir->etc))
{
warning(wi, "You are not using center of mass motion removal (mdp option comm-mode), numerical rounding errors can lead to build up of kinetic energy of the center of mass");
}
}
}
+ /* Check if combination rules used in LJ-PME are the same as in the force field */
+ if (EVDW_PME(ir->vdwtype))
+ {
+ check_combination_rules(ir, sys, wi);
+ }
+
/* Generalized reaction field */
if (ir->opts.ngtc == 0)
{
CHECK((ir->coulombtype == eelGRF) && (ir->opts.ref_t[0] <= 0));
}
- if (ir->eI == eiSD1 &&
- (gmx_mtop_ftype_count(sys, F_CONSTR) > 0 ||
- gmx_mtop_ftype_count(sys, F_SETTLE) > 0))
+ if (ir->eI == eiSD2)
{
- sprintf(warn_buf, "With constraints integrator %s is less accurate, consider using %s instead", ei_names[ir->eI], ei_names[eiSD2]);
+ sprintf(warn_buf, "The stochastic dynamics integrator %s is deprecated, since\n"
+ "it is slower than integrator %s and is slightly less accurate\n"
+ "with constraints. Use the %s integrator.",
+ ei_names[ir->eI], ei_names[eiSD1], ei_names[eiSD1]);
warning_note(wi, warn_buf);
}
if (ir->ePull != epullNO)
{
- if (ir->pull->grp[0].nat == 0)
+ gmx_bool bPullAbsoluteRef;
+
+ bPullAbsoluteRef = FALSE;
+ for (i = 0; i < ir->pull->ncoord; i++)
+ {
+ bPullAbsoluteRef = bPullAbsoluteRef ||
+ ir->pull->coord[i].group[0] == 0 ||
+ ir->pull->coord[i].group[1] == 0;
+ }
+ if (bPullAbsoluteRef)
{
absolute_reference(ir, sys, FALSE, AbsRef);
for (m = 0; m < DIM; m++)
if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
ir->deform[i][m] != 0)
{
- for (g = 1; g < ir->pull->ngrp; g++)
+ for (c = 0; c < ir->pull->ncoord; c++)
{
- if (ir->pull->grp[g].vec[m] != 0)
+ if (ir->pull->coord[c].vec[m] != 0)
{
gmx_fatal(FARGS, "Can not have dynamic box while using pull geometry '%s' (dim %c)", EPULLGEOM(ir->pull->eGeom), 'x'+m);
}
}
}
+ if (bConstr && ir->epc == epcMTTK)
+ {
+ warning_note(wi, "MTTK with constraints is deprecated, and will be removed in GROMACS 5.1");
+ }
+
if (ir->LincsWarnAngle > 90.0)
{
sprintf(warn_buf, "lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
* since user defined interactions might purposely
* not be zero at the cut-off.
*/
- if ((EVDW_IS_ZERO_AT_CUTOFF(ir->vdwtype) ||
- ir->vdw_modifier != eintmodNONE) &&
+ if (ir_vdw_is_zero_at_cutoff(ir) &&
rvdw1 + rvdw2 > ir->rlistlong - ir->rvdw)
{
sprintf(warn_buf, "The sum of the two largest charge group "
warning_note(wi, warn_buf);
}
}
- if ((EEL_IS_ZERO_AT_CUTOFF(ir->coulombtype) ||
- ir->coulomb_modifier != eintmodNONE) &&
+ if (ir_coulomb_is_zero_at_cutoff(ir) &&
rcoul1 + rcoul2 > ir->rlistlong - ir->rcoulomb)
{
sprintf(warn_buf, "The sum of the two largest charge group radii (%f) is larger than %s (%f) - rcoulomb (%f).\n"
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff