*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2016,2017,2018,2019, by the GROMACS development team, led by
+ * Copyright (c) 2013-2019, 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 "gromacs/utility/textwriter.h"
#define MAXPTR 254
-#define NOGID 255
+#define NOGID 255
/* Resource parameters
* Do not change any of these until you read the instruction
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 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];
+ 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];
} gmx_inputrec_strings;
-static gmx_inputrec_strings *is = nullptr;
+static gmx_inputrec_strings* is = nullptr;
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.");
+ 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);
}
}
-enum {
+enum
+{
egrptpALL, /* All particles have to be a member of a group. */
egrptpALL_GENREST, /* A rest group with name is generated for particles *
* that are not part of any group. */
* make a rest group for the remaining particles. */
};
-static const char *constraints[eshNR+1] = {
- "none", "h-bonds", "all-bonds", "h-angles", "all-angles", nullptr
-};
+static const char* constraints[eshNR + 1] = { "none", "h-bonds", "all-bonds",
+ "h-angles", "all-angles", nullptr };
-static const char *couple_lam[ecouplamNR+1] = {
- "vdw-q", "vdw", "q", "none", nullptr
-};
+static const char* couple_lam[ecouplamNR + 1] = { "vdw-q", "vdw", "q", "none", nullptr };
-static void GetSimTemps(int ntemps, t_simtemp *simtemp, double *temperature_lambdas)
+static void GetSimTemps(int ntemps, t_simtemp* simtemp, double* temperature_lambdas)
{
int i;
/* simple linear scaling -- allows more control */
if (simtemp->eSimTempScale == esimtempLINEAR)
{
- simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*temperature_lambdas[i];
+ simtemp->temperatures[i] =
+ simtemp->simtemp_low
+ + (simtemp->simtemp_high - simtemp->simtemp_low) * temperature_lambdas[i];
}
- else if (simtemp->eSimTempScale == esimtempGEOMETRIC) /* should give roughly equal acceptance for constant heat capacity . . . */
+ else if (simtemp->eSimTempScale
+ == esimtempGEOMETRIC) /* should give roughly equal acceptance for constant heat capacity . . . */
{
- simtemp->temperatures[i] = simtemp->simtemp_low * std::pow(simtemp->simtemp_high/simtemp->simtemp_low, static_cast<real>((1.0*i)/(ntemps-1)));
+ simtemp->temperatures[i] = simtemp->simtemp_low
+ * std::pow(simtemp->simtemp_high / simtemp->simtemp_low,
+ static_cast<real>((1.0 * i) / (ntemps - 1)));
}
else if (simtemp->eSimTempScale == esimtempEXPONENTIAL)
{
- simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*(std::expm1(temperature_lambdas[i])/std::expm1(1.0));
+ simtemp->temperatures[i] = simtemp->simtemp_low
+ + (simtemp->simtemp_high - simtemp->simtemp_low)
+ * (std::expm1(temperature_lambdas[i]) / std::expm1(1.0));
}
else
{
}
-
-static void _low_check(bool b, const char *s, warninp_t wi)
+static void _low_check(bool b, const char* s, warninp_t wi)
{
if (b)
{
}
}
-static void check_nst(const char *desc_nst, int nst,
- const char *desc_p, int *p,
- warninp_t wi)
+static void check_nst(const char* desc_nst, int nst, const char* desc_p, int* p, warninp_t wi)
{
char buf[STRLEN];
if (*p > 0 && *p % nst != 0)
{
/* Round up to the next multiple of nst */
- *p = ((*p)/nst + 1)*nst;
- sprintf(buf, "%s should be a multiple of %s, changing %s to %d\n",
- desc_p, desc_nst, desc_p, *p);
+ *p = ((*p) / nst + 1) * nst;
+ sprintf(buf, "%s should be a multiple of %s, changing %s to %d\n", desc_p, desc_nst, desc_p, *p);
warning(wi, buf);
}
}
}
//! Convert legacy mdp entries to modern ones.
-static void process_interaction_modifier(int *eintmod)
+static void process_interaction_modifier(int* eintmod)
{
if (*eintmod == eintmodPOTSHIFT_VERLET_UNSUPPORTED)
{
}
}
-void check_ir(const char *mdparin, const gmx::MdModulesNotifier &mdModulesNotifier,
- t_inputrec *ir, t_gromppopts *opts, warninp_t wi)
+void check_ir(const char* mdparin,
+ const gmx::MdModulesNotifier& mdModulesNotifier,
+ t_inputrec* ir,
+ t_gromppopts* opts,
+ warninp_t wi)
/* 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
char err_buf[256], warn_buf[STRLEN];
int i, j;
real dt_pcoupl;
- t_lambda *fep = ir->fepvals;
- t_expanded *expand = ir->expandedvals;
+ t_lambda* fep = ir->fepvals;
+ t_expanded* expand = ir->expandedvals;
set_warning_line(wi, mdparin, -1);
if (ir->coulombtype == eelRF_NEC_UNSUPPORTED)
{
- sprintf(warn_buf, "%s electrostatics is no longer supported",
- eel_names[eelRF_NEC_UNSUPPORTED]);
+ sprintf(warn_buf, "%s electrostatics is no longer supported", eel_names[eelRF_NEC_UNSUPPORTED]);
warning_error(wi, warn_buf);
}
{
warning_error(wi, "rvdw should be >= 0");
}
- if (ir->rlist < 0 &&
- !(ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0))
+ if (ir->rlist < 0 && !(ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0))
{
warning_error(wi, "rlist should be >= 0");
}
- sprintf(err_buf, "nstlist can not be smaller than 0. (If you were trying to use the heuristic neighbour-list update scheme for efficient buffering for improved energy conservation, please use the Verlet cut-off scheme instead.)");
+ sprintf(err_buf,
+ "nstlist can not be smaller than 0. (If you were trying to use the heuristic "
+ "neighbour-list update scheme for efficient buffering for improved energy "
+ "conservation, please use the Verlet cut-off scheme instead.)");
CHECK(ir->nstlist < 0);
process_interaction_modifier(&ir->coulomb_modifier);
{
// Since we have PME coulomb + LJ cut-off kernels with rcoulomb>rvdw
// for PME load balancing, we can support this exception.
- bool bUsesPmeTwinRangeKernel = (EEL_PME_EWALD(ir->coulombtype) &&
- ir->vdwtype == evdwCUT &&
- ir->rcoulomb > ir->rvdw);
+ bool bUsesPmeTwinRangeKernel = (EEL_PME_EWALD(ir->coulombtype) && ir->vdwtype == evdwCUT
+ && ir->rcoulomb > ir->rvdw);
if (!bUsesPmeTwinRangeKernel)
{
- warning_error(wi, "With Verlet lists rcoulomb!=rvdw is not supported (except for rcoulomb>rvdw with PME electrostatics)");
+ warning_error(wi,
+ "With Verlet lists rcoulomb!=rvdw is not supported (except for "
+ "rcoulomb>rvdw with PME electrostatics)");
}
}
if (ir->vdwtype == evdwSHIFT || ir->vdwtype == evdwSWITCH)
{
- if (ir->vdw_modifier == eintmodNONE ||
- ir->vdw_modifier == eintmodPOTSHIFT)
+ 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]);
+ 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]);
+ 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 and PME 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) ||
- EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD))
+ if (!(ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype) || EEL_PME(ir->coulombtype)
+ || ir->coulombtype == eelEWALD))
{
- warning_error(wi, "With Verlet lists only cut-off, reaction-field, PME and Ewald electrostatics are supported");
+ 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))
+ if (!(ir->coulomb_modifier == eintmodNONE || ir->coulomb_modifier == eintmodPOTSHIFT))
{
sprintf(warn_buf, "coulomb_modifier=%s is not supported", eintmod_names[ir->coulomb_modifier]);
warning_error(wi, warn_buf);
if (EEL_USER(ir->coulombtype))
{
- sprintf(warn_buf, "Coulomb type %s is not supported with the verlet scheme", eel_names[ir->coulombtype]);
+ sprintf(warn_buf, "Coulomb type %s is not supported with the verlet scheme",
+ eel_names[ir->coulombtype]);
warning_error(wi, warn_buf);
}
if (ir->nstlist < 10)
{
- warning_note(wi, "With Verlet lists the optimal nstlist is >= 10, with GPUs >= 20. Note that with the Verlet scheme, nstlist has no effect on the accuracy of your simulation.");
+ warning_note(wi,
+ "With Verlet lists the optimal nstlist is >= 10, with GPUs >= 20. Note "
+ "that with the Verlet scheme, nstlist has no effect on the accuracy of "
+ "your simulation.");
}
rc_max = std::max(ir->rvdw, ir->rcoulomb);
if (ir->rlist < rc_max)
{
- warning_error(wi, "With verlet lists rlist can not be smaller than rvdw or rcoulomb");
+ warning_error(wi,
+ "With verlet lists rlist can not be smaller than rvdw or rcoulomb");
}
if (ir->rlist == rc_max && ir->nstlist > 1)
{
- warning_note(wi, "rlist is equal to rvdw and/or rcoulomb: there is no explicit Verlet buffer. The cluster pair list does have a buffering effect, but choosing a larger rlist might be necessary for good energy conservation.");
+ warning_note(
+ wi,
+ "rlist is equal to rvdw and/or rcoulomb: there is no explicit Verlet "
+ "buffer. The cluster pair list does have a buffering effect, but choosing "
+ "a larger rlist might be necessary for good energy conservation.");
}
}
else
{
if (ir->rlist > rc_max)
{
- 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.");
+ 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 (inputrec2nboundeddim(ir) < 3)
{
- 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.");
+ 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.0 + verlet_buffer_ratio_nodynamics)*rc_max;
+ ir->rlist = (1.0 + verlet_buffer_ratio_nodynamics) * rc_max;
}
}
}
{
if (EI_RANDOM(ir->eI))
{
- sprintf(warn_buf, "Setting tcoupl from '%s' to 'no'. %s handles temperature coupling implicitly. See the documentation for more information on which parameters affect temperature for %s.", etcoupl_names[ir->etc], ei_names[ir->eI], ei_names[ir->eI]);
+ sprintf(warn_buf,
+ "Setting tcoupl from '%s' to 'no'. %s handles temperature coupling "
+ "implicitly. See the documentation for more information on which "
+ "parameters affect temperature for %s.",
+ etcoupl_names[ir->etc], ei_names[ir->eI], ei_names[ir->eI]);
}
else
{
- sprintf(warn_buf, "Setting tcoupl from '%s' to 'no'. Temperature coupling does not apply to %s.", etcoupl_names[ir->etc], ei_names[ir->eI]);
+ sprintf(warn_buf,
+ "Setting tcoupl from '%s' to 'no'. Temperature coupling does not apply to "
+ "%s.",
+ etcoupl_names[ir->etc], ei_names[ir->eI]);
}
warning_note(wi, warn_buf);
}
}
if (ir->eI == eiVVAK)
{
- sprintf(warn_buf, "Integrator method %s is implemented primarily for validation purposes; for molecular dynamics, you should probably be using %s or %s", ei_names[eiVVAK], ei_names[eiMD], ei_names[eiVV]);
+ sprintf(warn_buf,
+ "Integrator method %s is implemented primarily for validation purposes; for "
+ "molecular dynamics, you should probably be using %s or %s",
+ ei_names[eiVVAK], ei_names[eiMD], ei_names[eiVV]);
warning_note(wi, warn_buf);
}
if (!EI_DYNAMICS(ir->eI))
{
if (ir->epc != epcNO)
{
- sprintf(warn_buf, "Setting pcoupl from '%s' to 'no'. Pressure coupling does not apply to %s.", epcoupl_names[ir->epc], ei_names[ir->eI]);
+ sprintf(warn_buf,
+ "Setting pcoupl from '%s' to 'no'. Pressure coupling does not apply to %s.",
+ epcoupl_names[ir->epc], ei_names[ir->eI]);
warning_note(wi, warn_buf);
}
ir->epc = epcNO;
}
}
}
- else if ( (ir->nstenergy > 0 && ir->nstcalcenergy > ir->nstenergy) ||
- (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
- (ir->nstcalcenergy > ir->fepvals->nstdhdl) ) )
+ else if ((ir->nstenergy > 0 && ir->nstcalcenergy > ir->nstenergy)
+ || (ir->efep != efepNO && ir->fepvals->nstdhdl > 0
+ && (ir->nstcalcenergy > ir->fepvals->nstdhdl)))
{
- const char *nsten = "nstenergy";
- const char *nstdh = "nstdhdl";
- const char *min_name = nsten;
+ const char* nsten = "nstenergy";
+ const char* nstdh = "nstdhdl";
+ const char* min_name = nsten;
int min_nst = ir->nstenergy;
/* find the smallest of ( nstenergy, nstdhdl ) */
- if (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
- (ir->nstenergy == 0 || ir->fepvals->nstdhdl < ir->nstenergy))
+ if (ir->efep != efepNO && ir->fepvals->nstdhdl > 0
+ && (ir->nstenergy == 0 || ir->fepvals->nstdhdl < ir->nstenergy))
{
min_nst = ir->fepvals->nstdhdl;
min_name = nstdh;
}
/* If the user sets nstenergy small, we should respect that */
- sprintf(warn_buf,
- "Setting nstcalcenergy (%d) equal to %s (%d)",
- ir->nstcalcenergy, min_name, min_nst);
+ sprintf(warn_buf, "Setting nstcalcenergy (%d) equal to %s (%d)", ir->nstcalcenergy,
+ min_name, min_nst);
warning_note(wi, warn_buf);
ir->nstcalcenergy = min_nst;
}
if (ir->efep != efepNO)
{
/* nstdhdl should be a multiple of nstcalcenergy */
- check_nst("nstcalcenergy", ir->nstcalcenergy,
- "nstdhdl", &ir->fepvals->nstdhdl, wi);
+ check_nst("nstcalcenergy", ir->nstcalcenergy, "nstdhdl", &ir->fepvals->nstdhdl, wi);
}
if (ir->bExpanded)
{
/* nstexpanded should be a multiple of nstcalcenergy */
- check_nst("nstcalcenergy", ir->nstcalcenergy,
- "nstexpanded", &ir->expandedvals->nstexpanded, wi);
+ check_nst("nstcalcenergy", ir->nstcalcenergy, "nstexpanded",
+ &ir->expandedvals->nstexpanded, wi);
}
/* for storing exact averages nstenergy should be
* a multiple of nstcalcenergy
*/
- check_nst("nstcalcenergy", ir->nstcalcenergy,
- "nstenergy", &ir->nstenergy, wi);
+ check_nst("nstcalcenergy", ir->nstcalcenergy, "nstenergy", &ir->nstenergy, wi);
}
// Inquire all MdModules, if their parameters match with the energy
mdModulesNotifier.notifier_.notify(&energyCalculationFrequencyErrors);
// Emit all errors from the energy calculation frequency checks
- for (const std::string &energyFrequencyErrorMessage : energyCalculationFrequencyErrors.errorMessages())
+ for (const std::string& energyFrequencyErrorMessage :
+ energyCalculationFrequencyErrors.errorMessages())
{
warning_error(wi, energyFrequencyErrorMessage);
}
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.");
+ 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)
+ if ((EI_SD(ir->eI) || ir->eI == eiBD) && ir->bContinuation && ir->ld_seed != -1)
{
- warning_note(wi, "You are doing a continuation with SD or BD, make sure that ld_seed is different from the previous run (using ld_seed=-1 will ensure this)");
+ warning_note(wi,
+ "You are doing a continuation with SD or BD, make sure that ld_seed is "
+ "different from the previous run (using ld_seed=-1 will ensure this)");
}
/* TPI STUFF */
}
/* SHAKE / LINCS */
- if ( (opts->nshake > 0) && (opts->bMorse) )
+ if ((opts->nshake > 0) && (opts->bMorse))
{
- sprintf(warn_buf,
- "Using morse bond-potentials while constraining bonds is useless");
+ sprintf(warn_buf, "Using morse bond-potentials while constraining bonds is useless");
warning(wi, warn_buf);
}
- if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
- ir->bContinuation && ir->ld_seed != -1)
+ if ((EI_SD(ir->eI) || ir->eI == eiBD) && ir->bContinuation && ir->ld_seed != -1)
{
- warning_note(wi, "You are doing a continuation with SD or BD, make sure that ld_seed is different from the previous run (using ld_seed=-1 will ensure this)");
+ warning_note(wi,
+ "You are doing a continuation with SD or BD, make sure that ld_seed is "
+ "different from the previous run (using ld_seed=-1 will ensure this)");
}
/* verify simulated tempering options */
bool bAllTempZero = TRUE;
for (i = 0; i < fep->n_lambda; i++)
{
- sprintf(err_buf, "Entry %d for %s must be between 0 and 1, instead is %g", i, efpt_names[efptTEMPERATURE], fep->all_lambda[efptTEMPERATURE][i]);
+ sprintf(err_buf, "Entry %d for %s must be between 0 and 1, instead is %g", i,
+ efpt_names[efptTEMPERATURE], fep->all_lambda[efptTEMPERATURE][i]);
CHECK((fep->all_lambda[efptTEMPERATURE][i] < 0) || (fep->all_lambda[efptTEMPERATURE][i] > 1));
if (fep->all_lambda[efptTEMPERATURE][i] > 0)
{
if (ir->etc == etcNOSEHOOVER)
{
- sprintf(warn_buf, "Nose-Hoover based temperature control such as [%s] my not be entirelyconsistent with simulated tempering", etcoupl_names[ir->etc]);
+ sprintf(warn_buf,
+ "Nose-Hoover based temperature control such as [%s] my not be "
+ "entirelyconsistent with simulated tempering",
+ etcoupl_names[ir->etc]);
warning_note(wi, warn_buf);
}
/* check that the temperatures make sense */
- sprintf(err_buf, "Higher simulated tempering temperature (%g) must be >= than the simulated tempering lower temperature (%g)", ir->simtempvals->simtemp_high, ir->simtempvals->simtemp_low);
+ sprintf(err_buf,
+ "Higher simulated tempering temperature (%g) must be >= than the simulated "
+ "tempering lower temperature (%g)",
+ ir->simtempvals->simtemp_high, ir->simtempvals->simtemp_low);
CHECK(ir->simtempvals->simtemp_high <= ir->simtempvals->simtemp_low);
- sprintf(err_buf, "Higher simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_high);
+ sprintf(err_buf, "Higher simulated tempering temperature (%g) must be >= zero",
+ ir->simtempvals->simtemp_high);
CHECK(ir->simtempvals->simtemp_high <= 0);
- sprintf(err_buf, "Lower simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_low);
+ sprintf(err_buf, "Lower simulated tempering temperature (%g) must be >= zero",
+ ir->simtempvals->simtemp_low);
CHECK(ir->simtempvals->simtemp_low <= 0);
}
if (ir->efep != efepNO)
{
fep = ir->fepvals;
- sprintf(err_buf, "The soft-core power is %d and can only be 1 or 2",
- fep->sc_power);
+ sprintf(err_buf, "The soft-core power is %d and can only be 1 or 2", fep->sc_power);
CHECK(fep->sc_alpha != 0 && fep->sc_power != 1 && fep->sc_power != 2);
sprintf(err_buf, "The soft-core sc-r-power is %d and can only be 6 or 48",
static_cast<int>(fep->sc_r_power));
CHECK(fep->sc_alpha != 0 && fep->sc_r_power != 6.0 && fep->sc_r_power != 48.0);
- sprintf(err_buf, "Can't use positive delta-lambda (%g) if initial state/lambda does not start at zero", fep->delta_lambda);
- CHECK(fep->delta_lambda > 0 && ((fep->init_fep_state > 0) || (fep->init_lambda > 0)));
+ sprintf(err_buf,
+ "Can't use positive delta-lambda (%g) if initial state/lambda does not start at "
+ "zero",
+ fep->delta_lambda);
+ CHECK(fep->delta_lambda > 0 && ((fep->init_fep_state > 0) || (fep->init_lambda > 0)));
- sprintf(err_buf, "Can't use positive delta-lambda (%g) with expanded ensemble simulations", fep->delta_lambda);
+ sprintf(err_buf, "Can't use positive 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)");
if (fep->n_lambda == 0)
{
/* Clear output in case of no states:*/
- sprintf(err_buf, "init-lambda-state set to %d: no lambda states are defined.", fep->init_fep_state);
+ sprintf(err_buf, "init-lambda-state set to %d: no lambda states are defined.",
+ fep->init_fep_state);
CHECK((fep->init_fep_state >= 0) && (fep->n_lambda == 0));
}
else
{
- sprintf(err_buf, "initial thermodynamic state %d does not exist, only goes to %d", fep->init_fep_state, fep->n_lambda-1);
+ sprintf(err_buf, "initial thermodynamic state %d does not exist, only goes to %d",
+ fep->init_fep_state, fep->n_lambda - 1);
CHECK((fep->init_fep_state >= fep->n_lambda));
}
- sprintf(err_buf, "Lambda state must be set, either with init-lambda-state or with init-lambda");
+ sprintf(err_buf,
+ "Lambda state must be set, either with init-lambda-state or with init-lambda");
CHECK((fep->init_fep_state < 0) && (fep->init_lambda < 0));
- sprintf(err_buf, "init-lambda=%g while init-lambda-state=%d. Lambda state must be set either with init-lambda-state or with init-lambda, but not both",
+ sprintf(err_buf,
+ "init-lambda=%g while init-lambda-state=%d. Lambda state must be set either with "
+ "init-lambda-state or with init-lambda, but not both",
fep->init_lambda, fep->init_fep_state);
CHECK((fep->init_fep_state >= 0) && (fep->init_lambda >= 0));
-
if ((fep->init_lambda >= 0) && (fep->delta_lambda == 0))
{
int n_lambda_terms;
}
if (n_lambda_terms > 1)
{
- sprintf(warn_buf, "If lambda vector states (fep-lambdas, coul-lambdas etc.) are set, don't use init-lambda to set lambda state (except for slow growth). Use init-lambda-state instead.");
+ sprintf(warn_buf,
+ "If lambda vector states (fep-lambdas, coul-lambdas etc.) are set, don't "
+ "use init-lambda to set lambda state (except for slow growth). Use "
+ "init-lambda-state instead.");
warning(wi, warn_buf);
}
if (n_lambda_terms < 2 && fep->n_lambda > 0)
{
warning_note(wi,
- "init-lambda is deprecated for setting lambda state (except for slow growth). Use init-lambda-state instead.");
+ "init-lambda is deprecated for setting lambda state (except for slow "
+ "growth). Use init-lambda-state instead.");
}
}
{
for (i = 0; i < fep->n_lambda; i++)
{
- sprintf(err_buf, "Entry %d for %s must be between 0 and 1, instead is %g", i, efpt_names[j], fep->all_lambda[j][i]);
+ sprintf(err_buf, "Entry %d for %s must be between 0 and 1, instead is %g", i,
+ efpt_names[j], fep->all_lambda[j][i]);
CHECK((fep->all_lambda[j][i] < 0) || (fep->all_lambda[j][i] > 1));
}
}
{
for (i = 0; i < fep->n_lambda; i++)
{
- sprintf(err_buf, "For state %d, vdw-lambdas (%f) is changing with vdw softcore, while coul-lambdas (%f) is nonzero without coulomb softcore: this will lead to crashes, and is not supported.", i, fep->all_lambda[efptVDW][i],
- fep->all_lambda[efptCOUL][i]);
- CHECK((fep->sc_alpha > 0) &&
- (((fep->all_lambda[efptCOUL][i] > 0.0) &&
- (fep->all_lambda[efptCOUL][i] < 1.0)) &&
- ((fep->all_lambda[efptVDW][i] > 0.0) &&
- (fep->all_lambda[efptVDW][i] < 1.0))));
+ sprintf(err_buf,
+ "For state %d, vdw-lambdas (%f) is changing with vdw softcore, while "
+ "coul-lambdas (%f) is nonzero without coulomb softcore: this will lead to "
+ "crashes, and is not supported.",
+ i, fep->all_lambda[efptVDW][i], fep->all_lambda[efptCOUL][i]);
+ CHECK((fep->sc_alpha > 0)
+ && (((fep->all_lambda[efptCOUL][i] > 0.0) && (fep->all_lambda[efptCOUL][i] < 1.0))
+ && ((fep->all_lambda[efptVDW][i] > 0.0) && (fep->all_lambda[efptVDW][i] < 1.0))));
}
}
{
real sigma, lambda, r_sc;
- sigma = 0.34;
+ sigma = 0.34;
/* Maximum estimate for A and B charges equal with lambda power 1 */
lambda = 0.5;
- r_sc = std::pow(lambda*fep->sc_alpha*std::pow(sigma/ir->rcoulomb, fep->sc_r_power) + 1.0, 1.0/fep->sc_r_power);
- sprintf(warn_buf, "With PME there is a minor soft core effect present at the cut-off, proportional to (LJsigma/rcoulomb)^%g. This could have a minor effect on energy conservation, but usually other effects dominate. With a common sigma value of %g nm the fraction of the particle-particle potential at the cut-off at lambda=%g is around %.1e, while ewald-rtol is %.1e.",
- fep->sc_r_power,
- sigma, lambda, r_sc - 1.0, ir->ewald_rtol);
+ r_sc = std::pow(lambda * fep->sc_alpha * std::pow(sigma / ir->rcoulomb, fep->sc_r_power) + 1.0,
+ 1.0 / fep->sc_r_power);
+ sprintf(warn_buf,
+ "With PME there is a minor soft core effect present at the cut-off, "
+ "proportional to (LJsigma/rcoulomb)^%g. This could have a minor effect on "
+ "energy conservation, but usually other effects dominate. With a common sigma "
+ "value of %g nm the fraction of the particle-particle potential at the cut-off "
+ "at lambda=%g is around %.1e, while ewald-rtol is %.1e.",
+ fep->sc_r_power, sigma, lambda, r_sc - 1.0, ir->ewald_rtol);
warning_note(wi, warn_buf);
}
if ((ir->bSimTemp) || (ir->efep == efepEXPANDED))
{
- fep = ir->fepvals;
+ fep = ir->fepvals;
/* checking equilibration of weights inputs for validity */
- sprintf(err_buf, "weight-equil-number-all-lambda (%d) is ignored if lmc-weights-equil is not equal to %s",
+ sprintf(err_buf,
+ "weight-equil-number-all-lambda (%d) is ignored if lmc-weights-equil is not equal "
+ "to %s",
expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
CHECK((expand->equil_n_at_lam > 0) && (expand->elmceq != elmceqNUMATLAM));
- sprintf(err_buf, "weight-equil-number-samples (%d) is ignored if lmc-weights-equil is not equal to %s",
+ sprintf(err_buf,
+ "weight-equil-number-samples (%d) is ignored if lmc-weights-equil is not equal to "
+ "%s",
expand->equil_samples, elmceq_names[elmceqSAMPLES]);
CHECK((expand->equil_samples > 0) && (expand->elmceq != elmceqSAMPLES));
- sprintf(err_buf, "weight-equil-number-steps (%d) is ignored if lmc-weights-equil is not equal to %s",
+ sprintf(err_buf,
+ "weight-equil-number-steps (%d) is ignored if lmc-weights-equil is not equal to %s",
expand->equil_steps, elmceq_names[elmceqSTEPS]);
CHECK((expand->equil_steps > 0) && (expand->elmceq != elmceqSTEPS));
- sprintf(err_buf, "weight-equil-wl-delta (%d) is ignored if lmc-weights-equil is not equal to %s",
+ sprintf(err_buf,
+ "weight-equil-wl-delta (%d) is ignored if lmc-weights-equil is not equal to %s",
expand->equil_samples, elmceq_names[elmceqWLDELTA]);
CHECK((expand->equil_wl_delta > 0) && (expand->elmceq != elmceqWLDELTA));
- sprintf(err_buf, "weight-equil-count-ratio (%f) is ignored if lmc-weights-equil is not equal to %s",
+ sprintf(err_buf,
+ "weight-equil-count-ratio (%f) is ignored if lmc-weights-equil is not equal to %s",
expand->equil_ratio, elmceq_names[elmceqRATIO]);
CHECK((expand->equil_ratio > 0) && (expand->elmceq != elmceqRATIO));
- sprintf(err_buf, "weight-equil-number-all-lambda (%d) must be a positive integer if lmc-weights-equil=%s",
+ sprintf(err_buf,
+ "weight-equil-number-all-lambda (%d) must be a positive integer if "
+ "lmc-weights-equil=%s",
expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
CHECK((expand->equil_n_at_lam <= 0) && (expand->elmceq == elmceqNUMATLAM));
- sprintf(err_buf, "weight-equil-number-samples (%d) must be a positive integer if lmc-weights-equil=%s",
+ sprintf(err_buf,
+ "weight-equil-number-samples (%d) must be a positive integer if "
+ "lmc-weights-equil=%s",
expand->equil_samples, elmceq_names[elmceqSAMPLES]);
CHECK((expand->equil_samples <= 0) && (expand->elmceq == elmceqSAMPLES));
- sprintf(err_buf, "weight-equil-number-steps (%d) must be a positive integer if lmc-weights-equil=%s",
+ sprintf(err_buf,
+ "weight-equil-number-steps (%d) must be a positive integer if lmc-weights-equil=%s",
expand->equil_steps, elmceq_names[elmceqSTEPS]);
CHECK((expand->equil_steps <= 0) && (expand->elmceq == elmceqSTEPS));
CHECK((expand->minvarmin <= 0));
sprintf(err_buf, "weight-c-range (%d) must be greater or equal to 0", expand->c_range);
CHECK((expand->c_range < 0));
- sprintf(err_buf, "init-lambda-state (%d) must be zero if lmc-forced-nstart (%d)> 0 and lmc-move != 'no'",
+ sprintf(err_buf,
+ "init-lambda-state (%d) must be zero if lmc-forced-nstart (%d)> 0 and lmc-move != "
+ "'no'",
fep->init_fep_state, expand->lmc_forced_nstart);
- CHECK((fep->init_fep_state != 0) && (expand->lmc_forced_nstart > 0) && (expand->elmcmove != elmcmoveNO));
+ CHECK((fep->init_fep_state != 0) && (expand->lmc_forced_nstart > 0)
+ && (expand->elmcmove != elmcmoveNO));
sprintf(err_buf, "lmc-forced-nstart (%d) must not be negative", expand->lmc_forced_nstart);
CHECK((expand->lmc_forced_nstart < 0));
- sprintf(err_buf, "init-lambda-state (%d) must be in the interval [0,number of lambdas)", fep->init_fep_state);
+ sprintf(err_buf, "init-lambda-state (%d) must be in the interval [0,number of lambdas)",
+ fep->init_fep_state);
CHECK((fep->init_fep_state < 0) || (fep->init_fep_state >= fep->n_lambda));
sprintf(err_buf, "init-wl-delta (%f) must be greater than or equal to 0", expand->init_wl_delta);
CHECK((expand->wl_scale <= 0) || (expand->wl_scale >= 1));
/* if there is no temperature control, we need to specify an MC temperature */
- if (!integratorHasReferenceTemperature(ir) && (expand->elmcmove != elmcmoveNO) && (expand->mc_temp <= 0.0))
+ if (!integratorHasReferenceTemperature(ir) && (expand->elmcmove != elmcmoveNO)
+ && (expand->mc_temp <= 0.0))
{
- sprintf(err_buf, "If there is no temperature control, and lmc-mcmove!='no', mc_temp must be set to a positive number");
+ sprintf(err_buf,
+ "If there is no temperature control, and lmc-mcmove!='no', mc_temp must be set "
+ "to a positive number");
warning_error(wi, err_buf);
}
if (expand->nstTij > 0)
// Avoid modulus by zero in the case that already triggered an error exit.
if (ir->nstlog != 0)
{
- sprintf(err_buf, "nst-transition-matrix (%d) must be an integer multiple of nstlog (%d)",
+ sprintf(err_buf,
+ "nst-transition-matrix (%d) must be an integer multiple of nstlog (%d)",
expand->nstTij, ir->nstlog);
CHECK((expand->nstTij % ir->nstlog) != 0);
}
}
else
{
- sprintf(err_buf, "Can not have pressure coupling with pbc=%s",
- epbc_names[ir->ePBC]);
+ sprintf(err_buf, "Can not have pressure coupling with pbc=%s", epbc_names[ir->ePBC]);
CHECK(ir->epc != epcNO);
}
sprintf(err_buf, "Can not have Ewald with pbc=%s", epbc_names[ir->ePBC]);
CHECK(EEL_FULL(ir->coulombtype));
- sprintf(err_buf, "Can not have dispersion correction with pbc=%s",
- epbc_names[ir->ePBC]);
+ sprintf(err_buf, "Can not have dispersion correction with pbc=%s", epbc_names[ir->ePBC]);
CHECK(ir->eDispCorr != edispcNO);
}
if (ir->rlist == 0.0)
{
- sprintf(err_buf, "can only have neighborlist cut-off zero (=infinite)\n"
+ sprintf(err_buf,
+ "can only have neighborlist cut-off zero (=infinite)\n"
"with coulombtype = %s or coulombtype = %s\n"
"without periodic boundary conditions (pbc = %s) and\n"
"rcoulomb and rvdw set to zero",
eel_names[eelCUT], eel_names[eelUSER], epbc_names[epbcNONE]);
- CHECK(((ir->coulombtype != eelCUT) && (ir->coulombtype != eelUSER)) ||
- (ir->ePBC != epbcNONE) ||
- (ir->rcoulomb != 0.0) || (ir->rvdw != 0.0));
+ CHECK(((ir->coulombtype != eelCUT) && (ir->coulombtype != eelUSER))
+ || (ir->ePBC != epbcNONE) || (ir->rcoulomb != 0.0) || (ir->rvdw != 0.0));
if (ir->nstlist > 0)
{
- warning_note(wi, "Simulating without cut-offs can be (slightly) faster with nstlist=0, nstype=simple and only one MPI rank");
+ warning_note(wi,
+ "Simulating without cut-offs can be (slightly) faster with nstlist=0, "
+ "nstype=simple and only one MPI rank");
}
}
// helpful for a few years, we should reject such input,
// lest we have to support every historical decision
// forever.
- warning(wi, "If you want to remove the rotation around the center of mass, you should set comm_mode = Angular instead of setting nstcomm < 0. nstcomm is modified to its absolute value");
+ warning(wi,
+ "If you want to remove the rotation around the center of mass, you should set "
+ "comm_mode = Angular instead of setting nstcomm < 0. nstcomm is modified to "
+ "its absolute value");
ir->nstcomm = abs(ir->nstcomm);
}
if (ir->nstcalcenergy > 0 && ir->nstcomm < ir->nstcalcenergy)
{
- warning_note(wi, "nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting nstcomm to nstcalcenergy");
+ warning_note(wi,
+ "nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting "
+ "nstcomm to nstcalcenergy");
ir->nstcomm = ir->nstcalcenergy;
}
if (ir->comm_mode == ecmANGULAR)
{
- sprintf(err_buf, "Can not remove the rotation around the center of mass with periodic molecules");
+ sprintf(err_buf,
+ "Can not remove the rotation around the center of mass with periodic "
+ "molecules");
CHECK(ir->bPeriodicMols);
if (ir->ePBC != epbcNONE)
{
- warning(wi, "Removing the rotation around the center of mass in a periodic system, this can lead to artifacts. Only use this on a single (cluster of) molecules. This cluster should not cross periodic boundaries.");
+ warning(wi,
+ "Removing the rotation around the center of mass in a periodic system, "
+ "this can lead to artifacts. Only use this on a single (cluster of) "
+ "molecules. This cluster should not cross periodic boundaries.");
}
}
}
if (EI_STATE_VELOCITY(ir->eI) && !EI_SD(ir->eI) && ir->ePBC == epbcNONE && ir->comm_mode != ecmANGULAR)
{
- sprintf(warn_buf, "Tumbling and flying ice-cubes: We are not removing rotation around center of mass in a non-periodic system. You should probably set comm_mode = ANGULAR or use integrator = %s.", ei_names[eiSD1]);
+ sprintf(warn_buf,
+ "Tumbling and flying ice-cubes: We are not removing rotation around center of mass "
+ "in a non-periodic system. You should probably set comm_mode = ANGULAR or use "
+ "integrator = %s.",
+ ei_names[eiSD1]);
warning_note(wi, warn_buf);
}
if (ir->etc == etcYES)
{
ir->etc = etcBERENDSEN;
- warning_note(wi, "Old option for temperature coupling given: "
+ warning_note(wi,
+ "Old option for temperature coupling given: "
"changing \"yes\" to \"Berendsen\"\n");
}
{
if (ir->opts.nhchainlength < 1)
{
- sprintf(warn_buf, "number of Nose-Hoover chains (currently %d) cannot be less than 1,reset to 1\n", ir->opts.nhchainlength);
+ sprintf(warn_buf,
+ "number of Nose-Hoover chains (currently %d) cannot be less than 1,reset to "
+ "1\n",
+ ir->opts.nhchainlength);
ir->opts.nhchainlength = 1;
warning(wi, warn_buf);
}
if (ir->etc == etcNOSEHOOVER && !EI_VV(ir->eI) && ir->opts.nhchainlength > 1)
{
- warning_note(wi, "leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1");
+ warning_note(
+ wi,
+ "leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1");
ir->opts.nhchainlength = 1;
}
}
if (ir->eI == eiVVAK)
{
- sprintf(err_buf, "%s implemented primarily for validation, and requires nsttcouple = 1 and nstpcouple = 1.",
+ sprintf(err_buf,
+ "%s implemented primarily for validation, and requires nsttcouple = 1 and "
+ "nstpcouple = 1.",
ei_names[eiVVAK]);
CHECK((ir->nsttcouple != 1) || (ir->nstpcouple != 1));
}
if (ETC_ANDERSEN(ir->etc))
{
- sprintf(err_buf, "%s temperature control not supported for integrator %s.", etcoupl_names[ir->etc], ei_names[ir->eI]);
+ sprintf(err_buf, "%s temperature control not supported for integrator %s.",
+ etcoupl_names[ir->etc], ei_names[ir->eI]);
CHECK(!(EI_VV(ir->eI)));
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(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]);
warning_note(wi, warn_buf);
}
- 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]);
+ 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));
}
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(warn_buf,
+ "The %s thermostat does not generate the correct kinetic energy distribution. You "
+ "might want to consider using the %s thermostat.",
ETCOUPLTYPE(ir->etc), ETCOUPLTYPE(etcVRESCALE));
warning_note(wi, warn_buf);
}
- if ((ir->etc == etcNOSEHOOVER || ETC_ANDERSEN(ir->etc))
- && ir->epc == epcBERENDSEN)
+ if ((ir->etc == etcNOSEHOOVER || ETC_ANDERSEN(ir->etc)) && ir->epc == epcBERENDSEN)
{
- sprintf(warn_buf, "Using Berendsen pressure coupling invalidates the "
+ sprintf(warn_buf,
+ "Using Berendsen pressure coupling invalidates the "
"true ensemble for the thermostat");
warning(wi, warn_buf);
}
if (ir->epc == epcISOTROPIC)
{
ir->epc = epcBERENDSEN;
- warning_note(wi, "Old option for pressure coupling given: "
+ warning_note(wi,
+ "Old option for pressure coupling given: "
"changing \"Isotropic\" to \"Berendsen\"\n");
}
if (ir->epc != epcNO)
{
- dt_pcoupl = ir->nstpcouple*ir->delta_t;
+ dt_pcoupl = ir->nstpcouple * ir->delta_t;
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) - 10*GMX_REAL_EPS)
+ 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)",
+ 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);
warning(wi, warn_buf);
}
- sprintf(err_buf, "compressibility must be > 0 when using pressure"
- " coupling %s\n", EPCOUPLTYPE(ir->epc));
- CHECK(ir->compress[XX][XX] < 0 || ir->compress[YY][YY] < 0 ||
- ir->compress[ZZ][ZZ] < 0 ||
- (trace(ir->compress) == 0 && ir->compress[YY][XX] <= 0 &&
- ir->compress[ZZ][XX] <= 0 && ir->compress[ZZ][YY] <= 0));
+ sprintf(err_buf,
+ "compressibility must be > 0 when using pressure"
+ " coupling %s\n",
+ EPCOUPLTYPE(ir->epc));
+ CHECK(ir->compress[XX][XX] < 0 || ir->compress[YY][YY] < 0 || ir->compress[ZZ][ZZ] < 0
+ || (trace(ir->compress) == 0 && ir->compress[YY][XX] <= 0 && ir->compress[ZZ][XX] <= 0
+ && ir->compress[ZZ][YY] <= 0));
if (epcPARRINELLORAHMAN == ir->epc && opts->bGenVel)
{
if (ir->coulombtype == eelSWITCH)
{
- sprintf(warn_buf, "coulombtype = %s is only for testing purposes and can lead to serious "
+ sprintf(warn_buf,
+ "coulombtype = %s is only for testing purposes and can lead to serious "
"artifacts, advice: use coulombtype = %s",
- eel_names[ir->coulombtype],
- eel_names[eelRF_ZERO]);
+ eel_names[ir->coulombtype], eel_names[eelRF_ZERO]);
warning(wi, warn_buf);
}
if (EEL_RF(ir->coulombtype) && ir->epsilon_rf == 1 && ir->epsilon_r != 1)
{
- sprintf(warn_buf, "epsilon-r = %g and epsilon-rf = 1 with reaction field, proceeding assuming old format and exchanging epsilon-r and epsilon-rf", ir->epsilon_r);
+ sprintf(warn_buf,
+ "epsilon-r = %g and epsilon-rf = 1 with reaction field, proceeding assuming old "
+ "format and exchanging epsilon-r and epsilon-rf",
+ ir->epsilon_r);
warning(wi, warn_buf);
ir->epsilon_rf = ir->epsilon_r;
ir->epsilon_r = 1.0;
if (ir->epsilon_r == 0)
{
sprintf(err_buf,
- "It is pointless to use long-range electrostatics with infinite relative permittivity."
- "Since you are effectively turning of electrostatics, a plain cutoff will be much faster.");
+ "It is pointless to use long-range electrostatics with infinite relative "
+ "permittivity."
+ "Since you are effectively turning of electrostatics, a plain cutoff will be much "
+ "faster.");
CHECK(EEL_FULL(ir->coulombtype));
}
if (ir->coulombtype == eelRF_ZERO && ir->epsilon_rf != 0)
{
- sprintf(warn_buf, "With coulombtype = %s, epsilon-rf must be 0, assuming you meant epsilon_rf=0",
+ sprintf(warn_buf,
+ "With coulombtype = %s, epsilon-rf must be 0, assuming you meant epsilon_rf=0",
eel_names[ir->coulombtype]);
warning(wi, warn_buf);
ir->epsilon_rf = 0.0;
}
sprintf(err_buf, "epsilon-rf must be >= epsilon-r");
- CHECK((ir->epsilon_rf < ir->epsilon_r && ir->epsilon_rf != 0) ||
- (ir->epsilon_r == 0));
+ CHECK((ir->epsilon_rf < ir->epsilon_r && ir->epsilon_rf != 0) || (ir->epsilon_r == 0));
if (ir->epsilon_rf == ir->epsilon_r)
{
sprintf(warn_buf, "Using epsilon-rf = epsilon-r with %s does not make sense",
if (ir_coulomb_switched(ir))
{
sprintf(err_buf,
- "With coulombtype = %s rcoulomb_switch must be < rcoulomb. Or, better: Use the potential modifier options!",
+ "With coulombtype = %s rcoulomb_switch must be < rcoulomb. Or, better: Use the "
+ "potential modifier options!",
eel_names[ir->coulombtype]);
CHECK(ir->rcoulomb_switch >= ir->rcoulomb);
}
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);
+ "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);
+ "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)
+ if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT || ir->vdwtype == evdwSWITCH
+ || ir->vdwtype == evdwSHIFT)
{
sprintf(warn_buf,
- "The switch/shift interaction settings are just for compatibility; you will get better "
+ "The switch/shift interaction settings are just for compatibility; you will get "
+ "better "
"performance from applying potential modifiers to your interactions!\n");
warning_note(wi, warn_buf);
}
if (ir->coulombtype == eelPMESWITCH || ir->coulomb_modifier == eintmodPOTSWITCH)
{
- if (ir->rcoulomb_switch/ir->rcoulomb < 0.9499)
+ if (ir->rcoulomb_switch / ir->rcoulomb < 0.9499)
{
- 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.",
+ 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->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.");
+ 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_FULL(ir->coulombtype))
{
- if (ir->coulombtype == eelPMESWITCH || ir->coulombtype == eelPMEUSER ||
- ir->coulombtype == eelPMEUSERSWITCH)
+ if (ir->coulombtype == eelPMESWITCH || ir->coulombtype == eelPMEUSER
+ || ir->coulombtype == eelPMEUSERSWITCH)
{
sprintf(err_buf, "With coulombtype = %s, rcoulomb must be <= rlist",
eel_names[ir->coulombtype]);
if (ir->pme_order < orderMin || ir->pme_order > orderMax)
{
- sprintf(warn_buf, "With coulombtype = %s, you should have %d <= pme-order <= %d", eel_names[ir->coulombtype], orderMin, orderMax);
+ sprintf(warn_buf, "With coulombtype = %s, you should have %d <= pme-order <= %d",
+ eel_names[ir->coulombtype], orderMin, orderMax);
warning_error(wi, warn_buf);
}
}
{
if (ir->ewald_geometry == eewg3D)
{
- sprintf(warn_buf, "With pbc=%s you should use ewald-geometry=%s",
- epbc_names[ir->ePBC], eewg_names[eewg3DC]);
+ sprintf(warn_buf, "With pbc=%s you should use ewald-geometry=%s", epbc_names[ir->ePBC],
+ eewg_names[eewg3DC]);
warning(wi, warn_buf);
}
/* This check avoids extra pbc coding for exclusion corrections */
sprintf(err_buf, "wall-ewald-zfac should be >= 2");
CHECK(ir->wall_ewald_zfac < 2);
}
- if ((ir->ewald_geometry == eewg3DC) && (ir->ePBC != epbcXY) &&
- EEL_FULL(ir->coulombtype))
+ if ((ir->ewald_geometry == eewg3DC) && (ir->ePBC != epbcXY) && EEL_FULL(ir->coulombtype))
{
sprintf(warn_buf, "With %s and ewald_geometry = %s you should use pbc = %s",
eel_names[ir->coulombtype], eewg_names[eewg3DC], epbc_names[epbcXY]);
warning_note(wi, warn_buf);
sprintf(warn_buf,
"With epsilon_surface > 0 you can only use domain decomposition "
- "when there are only small molecules with all bonds constrained (mdrun will check for this).");
+ "when there are only small molecules with all bonds constrained (mdrun will check "
+ "for this).");
warning_note(wi, warn_buf);
}
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)
+ 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.",
+ 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);
}
if (!(ir->vdw_modifier == eintmodNONE || ir->vdw_modifier == eintmodPOTSHIFT))
{
sprintf(err_buf, "With vdwtype = %s, the only supported modifiers are %s and %s",
- evdw_names[ir->vdwtype],
- eintmod_names[eintmodPOTSHIFT],
- eintmod_names[eintmodNONE]);
+ evdw_names[ir->vdwtype], eintmod_names[eintmodPOTSHIFT], eintmod_names[eintmodNONE]);
warning_error(wi, err_buf);
}
}
if (ir->vdwtype == evdwUSER && ir->eDispCorr != edispcNO)
{
- warning_note(wi, "You have selected user tables with dispersion correction, the dispersion will be corrected to -C6/r^6 beyond rvdw_switch (the tabulated interaction between rvdw_switch and rvdw will not be double counted). Make sure that you really want dispersion correction to -C6/r^6.");
+ warning_note(wi,
+ "You have selected user tables with dispersion correction, the dispersion "
+ "will be corrected to -C6/r^6 beyond rvdw_switch (the tabulated interaction "
+ "between rvdw_switch and rvdw will not be double counted). Make sure that you "
+ "really want dispersion correction to -C6/r^6.");
}
- if (ir->eI == eiLBFGS && (ir->coulombtype == eelCUT || ir->vdwtype == evdwCUT)
- && ir->rvdw != 0)
+ if (ir->eI == eiLBFGS && (ir->coulombtype == eelCUT || ir->vdwtype == evdwCUT) && ir->rvdw != 0)
{
warning(wi, "For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
}
/* IMPLICIT SOLVENT */
if (ir->coulombtype == eelGB_NOTUSED)
{
- sprintf(warn_buf, "Invalid option %s for coulombtype",
- eel_names[ir->coulombtype]);
+ sprintf(warn_buf, "Invalid option %s for coulombtype", eel_names[ir->coulombtype]);
warning_error(wi, warn_buf);
}
str = the input string
n = the (pre-allocated) number of doubles read
r = the output array of doubles. */
-static void parse_n_real(char *str, int *n, real **r, warninp_t wi)
+static void parse_n_real(char* str, int* n, real** r, warninp_t wi)
{
auto values = gmx::splitString(str);
- *n = values.size();
+ *n = values.size();
snew(*r, *n);
for (int i = 0; i < *n; i++)
{
(*r)[i] = gmx::fromString<real>(values[i]);
}
- catch (gmx::GromacsException &)
+ catch (gmx::GromacsException&)
{
- warning_error(wi, "Invalid value " + values[i] + " in string in mdp file. Expected a real number.");
+ warning_error(wi, "Invalid value " + values[i]
+ + " in string in mdp file. Expected a real number.");
}
}
}
-static void do_fep_params(t_inputrec *ir, char fep_lambda[][STRLEN], char weights[STRLEN], warninp_t wi)
+static void do_fep_params(t_inputrec* ir, char fep_lambda[][STRLEN], char weights[STRLEN], warninp_t wi)
{
int i, j, max_n_lambda, nweights, nfep[efptNR];
- t_lambda *fep = ir->fepvals;
- t_expanded *expand = ir->expandedvals;
- real **count_fep_lambdas;
+ t_lambda* fep = ir->fepvals;
+ t_expanded* expand = ir->expandedvals;
+ real** count_fep_lambdas;
bool bOneLambda = TRUE;
snew(count_fep_lambdas, efptNR);
{
if (nfep[i] > max_n_lambda)
{
- max_n_lambda = nfep[i]; /* here's a nonzero one. All of them
- must have the same number if its not zero.*/
+ max_n_lambda = nfep[i]; /* here's a nonzero one. All of them
+ must have the same number if its not zero.*/
break;
}
}
}
else if (nfep[i] == max_n_lambda)
{
- if (i != efptTEMPERATURE) /* we treat this differently -- not really a reason to compute the derivative with
- respect to the temperature currently */
+ if (i != efptTEMPERATURE) /* we treat this differently -- not really a reason to compute
+ the derivative with respect to the temperature currently */
{
ir->fepvals->separate_dvdl[i] = TRUE;
}
}
else
{
- gmx_fatal(FARGS, "Number of lambdas (%d) for FEP type %s not equal to number of other types (%d)",
+ gmx_fatal(FARGS,
+ "Number of lambdas (%d) for FEP type %s not equal to number of other types "
+ "(%d)",
nfep[i], efpt_names[i], max_n_lambda);
}
}
for (i = 0; i < efptNR; i++)
{
snew(fep->all_lambda[i], fep->n_lambda);
- if (nfep[i] > 0) /* if it's zero, then the count_fep_lambda arrays
- are zero */
+ if (nfep[i] > 0) /* if it's zero, then the count_fep_lambda arrays
+ are zero */
{
for (j = 0; j < fep->n_lambda; j++)
{
}
sfree(count_fep_lambdas);
- /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for internal
- bookkeeping -- for now, init_lambda */
+ /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for
+ internal bookkeeping -- for now, init_lambda */
if ((nfep[efptFEP] == 0) && (fep->init_lambda >= 0))
{
{
if (fep->sc_alpha > 0.1)
{
- gmx_fatal(FARGS, "sc_alpha (%f) for sc_r_power = 48 should usually be between 0.001 and 0.004", fep->sc_alpha);
+ gmx_fatal(FARGS,
+ "sc_alpha (%f) for sc_r_power = 48 should usually be between 0.001 and 0.004",
+ fep->sc_alpha);
}
}
}
-static void do_simtemp_params(t_inputrec *ir)
+static void do_simtemp_params(t_inputrec* ir)
{
snew(ir->simtempvals->temperatures, ir->fepvals->n_lambda);
GetSimTemps(ir->fepvals->n_lambda, ir->simtempvals, ir->fepvals->all_lambda[efptTEMPERATURE]);
}
-static void
-convertYesNos(warninp_t /*wi*/, gmx::ArrayRef<const std::string> inputs, const char * /*name*/, gmx_bool *outputs)
+static void convertYesNos(warninp_t /*wi*/,
+ gmx::ArrayRef<const std::string> inputs,
+ const char* /*name*/,
+ gmx_bool* outputs)
{
int i = 0;
- for (const auto &input : inputs)
+ for (const auto& input : inputs)
{
outputs[i] = gmx::equalCaseInsensitive(input, "Y", 1);
++i;
}
}
-template <typename T> void
-convertInts(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char *name, T *outputs)
+template<typename T>
+void convertInts(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char* name, T* outputs)
{
int i = 0;
- for (const auto &input : inputs)
+ for (const auto& input : inputs)
{
try
{
outputs[i] = gmx::fromStdString<T>(input);
}
- catch (gmx::GromacsException &)
+ catch (gmx::GromacsException&)
{
- auto message = gmx::formatString("Invalid value for mdp option %s. %s should only consist of integers separated by spaces.",
- name, name);
+ auto message = gmx::formatString(
+ "Invalid value for mdp option %s. %s should only consist of integers separated "
+ "by spaces.",
+ name, name);
warning_error(wi, message);
}
++i;
}
}
-static void
-convertReals(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char *name, real *outputs)
+static void convertReals(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char* name, real* outputs)
{
int i = 0;
- for (const auto &input : inputs)
+ for (const auto& input : inputs)
{
try
{
outputs[i] = gmx::fromString<real>(input);
}
- catch (gmx::GromacsException &)
+ catch (gmx::GromacsException&)
{
- auto message = gmx::formatString("Invalid value for mdp option %s. %s should only consist of real numbers separated by spaces.",
- name, name);
+ auto message = gmx::formatString(
+ "Invalid value for mdp option %s. %s should only consist of real numbers "
+ "separated by spaces.",
+ name, name);
warning_error(wi, message);
}
++i;
}
}
-static void
-convertRvecs(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char *name, rvec *outputs)
+static void convertRvecs(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char* name, rvec* outputs)
{
int i = 0, d = 0;
- for (const auto &input : inputs)
+ for (const auto& input : inputs)
{
try
{
outputs[i][d] = gmx::fromString<real>(input);
}
- catch (gmx::GromacsException &)
+ catch (gmx::GromacsException&)
{
- auto message = gmx::formatString("Invalid value for mdp option %s. %s should only consist of real numbers separated by spaces.",
- name, name);
+ auto message = gmx::formatString(
+ "Invalid value for mdp option %s. %s should only consist of real numbers "
+ "separated by spaces.",
+ name, name);
warning_error(wi, message);
}
++d;
}
}
-static void do_wall_params(t_inputrec *ir,
- char *wall_atomtype, char *wall_density,
- t_gromppopts *opts,
- warninp_t wi)
+static void do_wall_params(t_inputrec* ir, char* wall_atomtype, char* wall_density, t_gromppopts* opts, warninp_t wi)
{
opts->wall_atomtype[0] = nullptr;
opts->wall_atomtype[1] = nullptr;
auto wallAtomTypes = gmx::splitString(wall_atomtype);
if (wallAtomTypes.size() != size_t(ir->nwall))
{
- gmx_fatal(FARGS, "Expected %d elements for wall_atomtype, found %zu",
- ir->nwall, wallAtomTypes.size());
+ gmx_fatal(FARGS, "Expected %d elements for wall_atomtype, found %zu", ir->nwall,
+ wallAtomTypes.size());
}
for (int i = 0; i < ir->nwall; i++)
{
auto wallDensity = gmx::splitString(wall_density);
if (wallDensity.size() != size_t(ir->nwall))
{
- gmx_fatal(FARGS, "Expected %d elements for wall-density, found %zu", ir->nwall, wallDensity.size());
+ gmx_fatal(FARGS, "Expected %d elements for wall-density, found %zu", ir->nwall,
+ wallDensity.size());
}
convertReals(wi, wallDensity, "wall-density", ir->wall_density);
for (int i = 0; i < ir->nwall; i++)
}
}
-static void add_wall_energrps(SimulationGroups *groups, int nwall, t_symtab *symtab)
+static void add_wall_energrps(SimulationGroups* groups, int nwall, t_symtab* symtab)
{
if (nwall > 0)
{
- AtomGroupIndices *grps = &(groups->groups[SimulationAtomGroupType::EnergyOutput]);
+ AtomGroupIndices* grps = &(groups->groups[SimulationAtomGroupType::EnergyOutput]);
for (int i = 0; i < nwall; i++)
{
- groups->groupNames.emplace_back(
- put_symtab(
- symtab,
- gmx::formatString("wall%d", i).c_str()));
+ groups->groupNames.emplace_back(put_symtab(symtab, gmx::formatString("wall%d", i).c_str()));
grps->emplace_back(groups->groupNames.size() - 1);
}
}
}
-static void read_expandedparams(std::vector<t_inpfile> *inp,
- t_expanded *expand, warninp_t wi)
+static void read_expandedparams(std::vector<t_inpfile>* inp, t_expanded* expand, warninp_t wi)
{
/* read expanded ensemble parameters */
printStringNewline(inp, "expanded ensemble variables");
expand->equil_wl_delta = get_ereal(inp, "weight-equil-wl-delta", -1, wi);
expand->equil_ratio = get_ereal(inp, "weight-equil-count-ratio", -1, wi);
printStringNewline(inp, "Seed for Monte Carlo in lambda space");
- expand->lmc_seed = get_eint(inp, "lmc-seed", -1, wi);
- expand->mc_temp = get_ereal(inp, "mc-temperature", -1, wi);
- expand->lmc_repeats = get_eint(inp, "lmc-repeats", 1, wi);
- expand->gibbsdeltalam = get_eint(inp, "lmc-gibbsdelta", -1, wi);
- expand->lmc_forced_nstart = get_eint(inp, "lmc-forced-nstart", 0, wi);
- expand->bSymmetrizedTMatrix = (get_eeenum(inp, "symmetrized-transition-matrix", yesno_names, wi) != 0);
- expand->nstTij = get_eint(inp, "nst-transition-matrix", -1, wi);
- expand->minvarmin = get_eint(inp, "mininum-var-min", 100, wi); /*default is reasonable */
- expand->c_range = get_eint(inp, "weight-c-range", 0, wi); /* default is just C=0 */
- expand->wl_scale = get_ereal(inp, "wl-scale", 0.8, wi);
- expand->wl_ratio = get_ereal(inp, "wl-ratio", 0.8, wi);
- expand->init_wl_delta = get_ereal(inp, "init-wl-delta", 1.0, wi);
- expand->bWLoneovert = (get_eeenum(inp, "wl-oneovert", yesno_names, wi) != 0);
+ expand->lmc_seed = get_eint(inp, "lmc-seed", -1, wi);
+ expand->mc_temp = get_ereal(inp, "mc-temperature", -1, wi);
+ expand->lmc_repeats = get_eint(inp, "lmc-repeats", 1, wi);
+ expand->gibbsdeltalam = get_eint(inp, "lmc-gibbsdelta", -1, wi);
+ expand->lmc_forced_nstart = get_eint(inp, "lmc-forced-nstart", 0, wi);
+ expand->bSymmetrizedTMatrix =
+ (get_eeenum(inp, "symmetrized-transition-matrix", yesno_names, wi) != 0);
+ expand->nstTij = get_eint(inp, "nst-transition-matrix", -1, wi);
+ expand->minvarmin = get_eint(inp, "mininum-var-min", 100, wi); /*default is reasonable */
+ expand->c_range = get_eint(inp, "weight-c-range", 0, wi); /* default is just C=0 */
+ expand->wl_scale = get_ereal(inp, "wl-scale", 0.8, wi);
+ expand->wl_ratio = get_ereal(inp, "wl-ratio", 0.8, wi);
+ expand->init_wl_delta = get_ereal(inp, "init-wl-delta", 1.0, wi);
+ expand->bWLoneovert = (get_eeenum(inp, "wl-oneovert", yesno_names, wi) != 0);
}
/*! \brief Return whether an end state with the given coupling-lambda
*/
static bool couple_lambda_has_vdw_on(int couple_lambda_value)
{
- return (couple_lambda_value == ecouplamVDW ||
- couple_lambda_value == ecouplamVDWQ);
+ return (couple_lambda_value == ecouplamVDW || couple_lambda_value == ecouplamVDWQ);
}
namespace
class MdpErrorHandler : public gmx::IKeyValueTreeErrorHandler
{
- public:
- explicit MdpErrorHandler(warninp_t wi)
- : wi_(wi), mapping_(nullptr)
- {
- }
+public:
+ explicit MdpErrorHandler(warninp_t wi) : wi_(wi), mapping_(nullptr) {}
- void setBackMapping(const gmx::IKeyValueTreeBackMapping &mapping)
- {
- mapping_ = &mapping;
- }
+ void setBackMapping(const gmx::IKeyValueTreeBackMapping& mapping) { mapping_ = &mapping; }
- bool onError(gmx::UserInputError *ex, const gmx::KeyValueTreePath &context) override
- {
- ex->prependContext(gmx::formatString("Error in mdp option \"%s\":",
- getOptionName(context).c_str()));
- std::string message = gmx::formatExceptionMessageToString(*ex);
- warning_error(wi_, message.c_str());
- return true;
- }
+ bool onError(gmx::UserInputError* ex, const gmx::KeyValueTreePath& context) override
+ {
+ ex->prependContext(
+ gmx::formatString("Error in mdp option \"%s\":", getOptionName(context).c_str()));
+ std::string message = gmx::formatExceptionMessageToString(*ex);
+ warning_error(wi_, message.c_str());
+ return true;
+ }
- private:
- std::string getOptionName(const gmx::KeyValueTreePath &context)
+private:
+ std::string getOptionName(const gmx::KeyValueTreePath& context)
+ {
+ if (mapping_ != nullptr)
{
- if (mapping_ != nullptr)
- {
- gmx::KeyValueTreePath path = mapping_->originalPath(context);
- GMX_ASSERT(path.size() == 1, "Inconsistent mapping back to mdp options");
- return path[0];
- }
- GMX_ASSERT(context.size() == 1, "Inconsistent context for mdp option parsing");
- return context[0];
+ gmx::KeyValueTreePath path = mapping_->originalPath(context);
+ GMX_ASSERT(path.size() == 1, "Inconsistent mapping back to mdp options");
+ return path[0];
}
+ GMX_ASSERT(context.size() == 1, "Inconsistent context for mdp option parsing");
+ return context[0];
+ }
- warninp_t wi_;
- const gmx::IKeyValueTreeBackMapping *mapping_;
+ warninp_t wi_;
+ const gmx::IKeyValueTreeBackMapping* mapping_;
};
} // namespace
-void get_ir(const char *mdparin, const char *mdparout,
- gmx::MDModules *mdModules, t_inputrec *ir, t_gromppopts *opts,
- WriteMdpHeader writeMdpHeader, warninp_t wi)
+void get_ir(const char* mdparin,
+ const char* mdparout,
+ gmx::MDModules* mdModules,
+ t_inputrec* ir,
+ t_gromppopts* opts,
+ WriteMdpHeader writeMdpHeader,
+ warninp_t wi)
{
- char *dumstr[2];
- double dumdub[2][6];
- int i, j, m;
- char warn_buf[STRLEN];
- t_lambda *fep = ir->fepvals;
- t_expanded *expand = ir->expandedvals;
+ char* dumstr[2];
+ double dumdub[2][6];
+ int i, j, m;
+ char warn_buf[STRLEN];
+ t_lambda* fep = ir->fepvals;
+ t_expanded* expand = ir->expandedvals;
- const char *no_names[] = { "no", nullptr };
+ const char* no_names[] = { "no", nullptr };
init_inputrec_strings();
gmx::TextInputFile stream(mdparin);
printStringNewline(&inp, "VARIOUS PREPROCESSING OPTIONS");
printStringNoNewline(&inp, "Preprocessor information: use cpp syntax.");
printStringNoNewline(&inp, "e.g.: -I/home/joe/doe -I/home/mary/roe");
- setStringEntry(&inp, "include", opts->include, nullptr);
- printStringNoNewline(&inp, "e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
- setStringEntry(&inp, "define", opts->define, nullptr);
+ setStringEntry(&inp, "include", opts->include, nullptr);
+ printStringNoNewline(
+ &inp, "e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
+ setStringEntry(&inp, "define", opts->define, nullptr);
printStringNewline(&inp, "RUN CONTROL PARAMETERS");
- ir->eI = get_eeenum(&inp, "integrator", ei_names, wi);
+ ir->eI = get_eeenum(&inp, "integrator", ei_names, wi);
printStringNoNewline(&inp, "Start time and timestep in ps");
ir->init_t = get_ereal(&inp, "tinit", 0.0, wi);
- ir->delta_t = get_ereal(&inp, "dt", 0.001, wi);
- ir->nsteps = get_eint64(&inp, "nsteps", 0, wi);
+ ir->delta_t = get_ereal(&inp, "dt", 0.001, wi);
+ ir->nsteps = get_eint64(&inp, "nsteps", 0, wi);
printStringNoNewline(&inp, "For exact run continuation or redoing part of a run");
- ir->init_step = get_eint64(&inp, "init-step", 0, wi);
- printStringNoNewline(&inp, "Part index is updated automatically on checkpointing (keeps files separate)");
+ ir->init_step = get_eint64(&inp, "init-step", 0, wi);
+ printStringNoNewline(
+ &inp, "Part index is updated automatically on checkpointing (keeps files separate)");
ir->simulation_part = get_eint(&inp, "simulation-part", 1, wi);
printStringNoNewline(&inp, "mode for center of mass motion removal");
- ir->comm_mode = get_eeenum(&inp, "comm-mode", ecm_names, wi);
+ ir->comm_mode = get_eeenum(&inp, "comm-mode", ecm_names, wi);
printStringNoNewline(&inp, "number of steps for center of mass motion removal");
- ir->nstcomm = get_eint(&inp, "nstcomm", 100, wi);
+ ir->nstcomm = get_eint(&inp, "nstcomm", 100, wi);
printStringNoNewline(&inp, "group(s) for center of mass motion removal");
- setStringEntry(&inp, "comm-grps", is->vcm, nullptr);
+ setStringEntry(&inp, "comm-grps", is->vcm, nullptr);
printStringNewline(&inp, "LANGEVIN DYNAMICS OPTIONS");
printStringNoNewline(&inp, "Friction coefficient (amu/ps) and random seed");
- ir->bd_fric = get_ereal(&inp, "bd-fric", 0.0, wi);
- ir->ld_seed = get_eint64(&inp, "ld-seed", -1, wi);
+ ir->bd_fric = get_ereal(&inp, "bd-fric", 0.0, wi);
+ ir->ld_seed = get_eint64(&inp, "ld-seed", -1, wi);
/* Em stuff */
printStringNewline(&inp, "ENERGY MINIMIZATION OPTIONS");
printStringNoNewline(&inp, "Force tolerance and initial step-size");
- ir->em_tol = get_ereal(&inp, "emtol", 10.0, wi);
+ ir->em_tol = get_ereal(&inp, "emtol", 10.0, wi);
ir->em_stepsize = get_ereal(&inp, "emstep", 0.01, wi);
printStringNoNewline(&inp, "Max number of iterations in relax-shells");
- ir->niter = get_eint(&inp, "niter", 20, wi);
+ ir->niter = get_eint(&inp, "niter", 20, wi);
printStringNoNewline(&inp, "Step size (ps^2) for minimization of flexible constraints");
ir->fc_stepsize = get_ereal(&inp, "fcstep", 0, wi);
printStringNoNewline(&inp, "Frequency of steepest descents steps when doing CG");
ir->nstcgsteep = get_eint(&inp, "nstcgsteep", 1000, wi);
- ir->nbfgscorr = get_eint(&inp, "nbfgscorr", 10, wi);
+ ir->nbfgscorr = get_eint(&inp, "nbfgscorr", 10, wi);
printStringNewline(&inp, "TEST PARTICLE INSERTION OPTIONS");
- ir->rtpi = get_ereal(&inp, "rtpi", 0.05, wi);
+ ir->rtpi = get_ereal(&inp, "rtpi", 0.05, wi);
/* Output options */
printStringNewline(&inp, "OUTPUT CONTROL OPTIONS");
printStringNoNewline(&inp, "Output frequency for coords (x), velocities (v) and forces (f)");
- ir->nstxout = get_eint(&inp, "nstxout", 0, wi);
- ir->nstvout = get_eint(&inp, "nstvout", 0, wi);
- ir->nstfout = get_eint(&inp, "nstfout", 0, wi);
+ ir->nstxout = get_eint(&inp, "nstxout", 0, wi);
+ ir->nstvout = get_eint(&inp, "nstvout", 0, wi);
+ ir->nstfout = get_eint(&inp, "nstfout", 0, wi);
printStringNoNewline(&inp, "Output frequency for energies to log file and energy file");
ir->nstlog = get_eint(&inp, "nstlog", 1000, wi);
ir->nstcalcenergy = get_eint(&inp, "nstcalcenergy", 100, wi);
- ir->nstenergy = get_eint(&inp, "nstenergy", 1000, wi);
+ ir->nstenergy = get_eint(&inp, "nstenergy", 1000, wi);
printStringNoNewline(&inp, "Output frequency and precision for .xtc file");
- ir->nstxout_compressed = get_eint(&inp, "nstxout-compressed", 0, wi);
+ ir->nstxout_compressed = get_eint(&inp, "nstxout-compressed", 0, wi);
ir->x_compression_precision = get_ereal(&inp, "compressed-x-precision", 1000.0, wi);
printStringNoNewline(&inp, "This selects the subset of atoms for the compressed");
printStringNoNewline(&inp, "trajectory file. You can select multiple groups. By");
printStringNoNewline(&inp, "default, all atoms will be written.");
setStringEntry(&inp, "compressed-x-grps", is->x_compressed_groups, nullptr);
printStringNoNewline(&inp, "Selection of energy groups");
- setStringEntry(&inp, "energygrps", is->energy, nullptr);
+ setStringEntry(&inp, "energygrps", is->energy, nullptr);
/* Neighbor searching */
printStringNewline(&inp, "NEIGHBORSEARCHING PARAMETERS");
- printStringNoNewline(&inp, "cut-off scheme (Verlet: particle based cut-offs, group: using charge groups)");
- ir->cutoff_scheme = get_eeenum(&inp, "cutoff-scheme", ecutscheme_names, wi);
+ printStringNoNewline(
+ &inp, "cut-off scheme (Verlet: particle based cut-offs, group: using charge groups)");
+ ir->cutoff_scheme = get_eeenum(&inp, "cutoff-scheme", ecutscheme_names, wi);
printStringNoNewline(&inp, "nblist update frequency");
- ir->nstlist = get_eint(&inp, "nstlist", 10, wi);
+ ir->nstlist = get_eint(&inp, "nstlist", 10, wi);
printStringNoNewline(&inp, "Periodic boundary conditions: xyz, no, xy");
- ir->ePBC = get_eeenum(&inp, "pbc", epbc_names, wi);
+ ir->ePBC = get_eeenum(&inp, "pbc", epbc_names, wi);
ir->bPeriodicMols = get_eeenum(&inp, "periodic-molecules", yesno_names, wi) != 0;
- printStringNoNewline(&inp, "Allowed energy error due to the Verlet buffer in kJ/mol/ps per atom,");
+ printStringNoNewline(&inp,
+ "Allowed energy error due to the Verlet buffer in kJ/mol/ps per atom,");
printStringNoNewline(&inp, "a value of -1 means: use rlist");
- ir->verletbuf_tol = get_ereal(&inp, "verlet-buffer-tolerance", 0.005, wi);
+ ir->verletbuf_tol = get_ereal(&inp, "verlet-buffer-tolerance", 0.005, wi);
printStringNoNewline(&inp, "nblist cut-off");
- ir->rlist = get_ereal(&inp, "rlist", 1.0, wi);
+ ir->rlist = get_ereal(&inp, "rlist", 1.0, wi);
printStringNoNewline(&inp, "long-range cut-off for switched potentials");
/* Electrostatics */
printStringNewline(&inp, "OPTIONS FOR ELECTROSTATICS AND VDW");
printStringNoNewline(&inp, "Method for doing electrostatics");
- ir->coulombtype = get_eeenum(&inp, "coulombtype", eel_names, wi);
- ir->coulomb_modifier = get_eeenum(&inp, "coulomb-modifier", eintmod_names, wi);
+ ir->coulombtype = get_eeenum(&inp, "coulombtype", eel_names, wi);
+ ir->coulomb_modifier = get_eeenum(&inp, "coulomb-modifier", eintmod_names, wi);
printStringNoNewline(&inp, "cut-off lengths");
- ir->rcoulomb_switch = get_ereal(&inp, "rcoulomb-switch", 0.0, wi);
- ir->rcoulomb = get_ereal(&inp, "rcoulomb", 1.0, wi);
- printStringNoNewline(&inp, "Relative dielectric constant for the medium and the reaction field");
- ir->epsilon_r = get_ereal(&inp, "epsilon-r", 1.0, wi);
+ ir->rcoulomb_switch = get_ereal(&inp, "rcoulomb-switch", 0.0, wi);
+ ir->rcoulomb = get_ereal(&inp, "rcoulomb", 1.0, wi);
+ printStringNoNewline(&inp,
+ "Relative dielectric constant for the medium and the reaction field");
+ ir->epsilon_r = get_ereal(&inp, "epsilon-r", 1.0, wi);
ir->epsilon_rf = get_ereal(&inp, "epsilon-rf", 0.0, wi);
printStringNoNewline(&inp, "Method for doing Van der Waals");
- ir->vdwtype = get_eeenum(&inp, "vdw-type", evdw_names, wi);
- ir->vdw_modifier = get_eeenum(&inp, "vdw-modifier", eintmod_names, wi);
+ ir->vdwtype = get_eeenum(&inp, "vdw-type", evdw_names, wi);
+ ir->vdw_modifier = get_eeenum(&inp, "vdw-modifier", eintmod_names, wi);
printStringNoNewline(&inp, "cut-off lengths");
- ir->rvdw_switch = get_ereal(&inp, "rvdw-switch", 0.0, wi);
- ir->rvdw = get_ereal(&inp, "rvdw", 1.0, wi);
+ ir->rvdw_switch = get_ereal(&inp, "rvdw-switch", 0.0, wi);
+ ir->rvdw = get_ereal(&inp, "rvdw", 1.0, wi);
printStringNoNewline(&inp, "Apply long range dispersion corrections for Energy and Pressure");
- ir->eDispCorr = get_eeenum(&inp, "DispCorr", edispc_names, wi);
+ ir->eDispCorr = get_eeenum(&inp, "DispCorr", edispc_names, wi);
printStringNoNewline(&inp, "Extension of the potential lookup tables beyond the cut-off");
ir->tabext = get_ereal(&inp, "table-extension", 1.0, wi);
printStringNoNewline(&inp, "Separate tables between energy group pairs");
- setStringEntry(&inp, "energygrp-table", is->egptable, nullptr);
+ setStringEntry(&inp, "energygrp-table", is->egptable, nullptr);
printStringNoNewline(&inp, "Spacing for the PME/PPPM FFT grid");
ir->fourier_spacing = get_ereal(&inp, "fourierspacing", 0.12, wi);
printStringNoNewline(&inp, "FFT grid size, when a value is 0 fourierspacing will be used");
- ir->nkx = get_eint(&inp, "fourier-nx", 0, wi);
- ir->nky = get_eint(&inp, "fourier-ny", 0, wi);
- ir->nkz = get_eint(&inp, "fourier-nz", 0, wi);
+ ir->nkx = get_eint(&inp, "fourier-nx", 0, wi);
+ ir->nky = get_eint(&inp, "fourier-ny", 0, wi);
+ ir->nkz = get_eint(&inp, "fourier-nz", 0, wi);
printStringNoNewline(&inp, "EWALD/PME/PPPM parameters");
- ir->pme_order = get_eint(&inp, "pme-order", 4, wi);
+ ir->pme_order = get_eint(&inp, "pme-order", 4, wi);
ir->ewald_rtol = get_ereal(&inp, "ewald-rtol", 0.00001, wi);
ir->ewald_rtol_lj = get_ereal(&inp, "ewald-rtol-lj", 0.001, wi);
ir->ljpme_combination_rule = get_eeenum(&inp, "lj-pme-comb-rule", eljpme_names, wi);
/* Coupling stuff */
printStringNewline(&inp, "OPTIONS FOR WEAK COUPLING ALGORITHMS");
printStringNoNewline(&inp, "Temperature coupling");
- ir->etc = get_eeenum(&inp, "tcoupl", etcoupl_names, wi);
- ir->nsttcouple = get_eint(&inp, "nsttcouple", -1, wi);
+ ir->etc = get_eeenum(&inp, "tcoupl", etcoupl_names, wi);
+ ir->nsttcouple = get_eint(&inp, "nsttcouple", -1, wi);
ir->opts.nhchainlength = get_eint(&inp, "nh-chain-length", 10, wi);
- ir->bPrintNHChains = (get_eeenum(&inp, "print-nose-hoover-chain-variables", yesno_names, wi) != 0);
+ ir->bPrintNHChains = (get_eeenum(&inp, "print-nose-hoover-chain-variables", yesno_names, wi) != 0);
printStringNoNewline(&inp, "Groups to couple separately");
- setStringEntry(&inp, "tc-grps", is->tcgrps, nullptr);
+ setStringEntry(&inp, "tc-grps", is->tcgrps, nullptr);
printStringNoNewline(&inp, "Time constant (ps) and reference temperature (K)");
- setStringEntry(&inp, "tau-t", is->tau_t, nullptr);
- setStringEntry(&inp, "ref-t", is->ref_t, nullptr);
+ setStringEntry(&inp, "tau-t", is->tau_t, nullptr);
+ setStringEntry(&inp, "ref-t", is->ref_t, nullptr);
printStringNoNewline(&inp, "pressure coupling");
- ir->epc = get_eeenum(&inp, "pcoupl", epcoupl_names, wi);
- ir->epct = get_eeenum(&inp, "pcoupltype", epcoupltype_names, wi);
- ir->nstpcouple = get_eint(&inp, "nstpcouple", -1, wi);
+ ir->epc = get_eeenum(&inp, "pcoupl", epcoupl_names, wi);
+ ir->epct = get_eeenum(&inp, "pcoupltype", epcoupltype_names, wi);
+ ir->nstpcouple = get_eint(&inp, "nstpcouple", -1, wi);
printStringNoNewline(&inp, "Time constant (ps), compressibility (1/bar) and reference P (bar)");
- ir->tau_p = get_ereal(&inp, "tau-p", 1.0, wi);
- setStringEntry(&inp, "compressibility", dumstr[0], nullptr);
- setStringEntry(&inp, "ref-p", dumstr[1], nullptr);
+ ir->tau_p = get_ereal(&inp, "tau-p", 1.0, wi);
+ setStringEntry(&inp, "compressibility", dumstr[0], nullptr);
+ setStringEntry(&inp, "ref-p", dumstr[1], nullptr);
printStringNoNewline(&inp, "Scaling of reference coordinates, No, All or COM");
ir->refcoord_scaling = get_eeenum(&inp, "refcoord-scaling", erefscaling_names, wi);
printStringNewline(&inp, "OPTIONS FOR QMMM calculations");
ir->bQMMM = (get_eeenum(&inp, "QMMM", yesno_names, wi) != 0);
printStringNoNewline(&inp, "Groups treated Quantum Mechanically");
- setStringEntry(&inp, "QMMM-grps", is->QMMM, nullptr);
+ setStringEntry(&inp, "QMMM-grps", is->QMMM, nullptr);
printStringNoNewline(&inp, "QM method");
- setStringEntry(&inp, "QMmethod", is->QMmethod, nullptr);
+ setStringEntry(&inp, "QMmethod", is->QMmethod, nullptr);
printStringNoNewline(&inp, "QMMM scheme");
- ir->QMMMscheme = get_eeenum(&inp, "QMMMscheme", eQMMMscheme_names, wi);
+ ir->QMMMscheme = get_eeenum(&inp, "QMMMscheme", eQMMMscheme_names, wi);
printStringNoNewline(&inp, "QM basisset");
- setStringEntry(&inp, "QMbasis", is->QMbasis, nullptr);
+ setStringEntry(&inp, "QMbasis", is->QMbasis, nullptr);
printStringNoNewline(&inp, "QM charge");
- setStringEntry(&inp, "QMcharge", is->QMcharge, nullptr);
+ setStringEntry(&inp, "QMcharge", is->QMcharge, nullptr);
printStringNoNewline(&inp, "QM multiplicity");
- setStringEntry(&inp, "QMmult", is->QMmult, nullptr);
+ setStringEntry(&inp, "QMmult", is->QMmult, nullptr);
printStringNoNewline(&inp, "Surface Hopping");
- setStringEntry(&inp, "SH", is->bSH, nullptr);
+ setStringEntry(&inp, "SH", is->bSH, nullptr);
printStringNoNewline(&inp, "CAS space options");
- setStringEntry(&inp, "CASorbitals", is->CASorbitals, nullptr);
- setStringEntry(&inp, "CASelectrons", is->CASelectrons, nullptr);
+ setStringEntry(&inp, "CASorbitals", is->CASorbitals, nullptr);
+ setStringEntry(&inp, "CASelectrons", is->CASelectrons, nullptr);
setStringEntry(&inp, "SAon", is->SAon, nullptr);
setStringEntry(&inp, "SAoff", is->SAoff, nullptr);
setStringEntry(&inp, "SAsteps", is->SAsteps, nullptr);
/* Simulated annealing */
printStringNewline(&inp, "SIMULATED ANNEALING");
printStringNoNewline(&inp, "Type of annealing for each temperature group (no/single/periodic)");
- setStringEntry(&inp, "annealing", is->anneal, nullptr);
- printStringNoNewline(&inp, "Number of time points to use for specifying annealing in each group");
+ setStringEntry(&inp, "annealing", is->anneal, nullptr);
+ printStringNoNewline(&inp,
+ "Number of time points to use for specifying annealing in each group");
setStringEntry(&inp, "annealing-npoints", is->anneal_npoints, nullptr);
printStringNoNewline(&inp, "List of times at the annealing points for each group");
- setStringEntry(&inp, "annealing-time", is->anneal_time, nullptr);
+ setStringEntry(&inp, "annealing-time", is->anneal_time, nullptr);
printStringNoNewline(&inp, "Temp. at each annealing point, for each group.");
- setStringEntry(&inp, "annealing-temp", is->anneal_temp, nullptr);
+ setStringEntry(&inp, "annealing-temp", is->anneal_temp, nullptr);
/* Startup run */
printStringNewline(&inp, "GENERATE VELOCITIES FOR STARTUP RUN");
- opts->bGenVel = (get_eeenum(&inp, "gen-vel", yesno_names, wi) != 0);
- opts->tempi = get_ereal(&inp, "gen-temp", 300.0, wi);
- opts->seed = get_eint(&inp, "gen-seed", -1, wi);
+ opts->bGenVel = (get_eeenum(&inp, "gen-vel", yesno_names, wi) != 0);
+ opts->tempi = get_ereal(&inp, "gen-temp", 300.0, wi);
+ opts->seed = get_eint(&inp, "gen-seed", -1, wi);
/* Shake stuff */
printStringNewline(&inp, "OPTIONS FOR BONDS");
- opts->nshake = get_eeenum(&inp, "constraints", constraints, wi);
+ opts->nshake = get_eeenum(&inp, "constraints", constraints, wi);
printStringNoNewline(&inp, "Type of constraint algorithm");
ir->eConstrAlg = get_eeenum(&inp, "constraint-algorithm", econstr_names, wi);
printStringNoNewline(&inp, "Do not constrain the start configuration");
ir->bContinuation = (get_eeenum(&inp, "continuation", yesno_names, wi) != 0);
- printStringNoNewline(&inp, "Use successive overrelaxation to reduce the number of shake iterations");
+ printStringNoNewline(&inp,
+ "Use successive overrelaxation to reduce the number of shake iterations");
ir->bShakeSOR = (get_eeenum(&inp, "Shake-SOR", yesno_names, wi) != 0);
printStringNoNewline(&inp, "Relative tolerance of shake");
ir->shake_tol = get_ereal(&inp, "shake-tol", 0.0001, wi);
/* Energy group exclusions */
printStringNewline(&inp, "ENERGY GROUP EXCLUSIONS");
- printStringNoNewline(&inp, "Pairs of energy groups for which all non-bonded interactions are excluded");
- setStringEntry(&inp, "energygrp-excl", is->egpexcl, nullptr);
+ printStringNoNewline(
+ &inp, "Pairs of energy groups for which all non-bonded interactions are excluded");
+ setStringEntry(&inp, "energygrp-excl", is->egpexcl, nullptr);
/* Walls */
printStringNewline(&inp, "WALLS");
- printStringNoNewline(&inp, "Number of walls, type, atom types, densities and box-z scale factor for Ewald");
+ printStringNoNewline(
+ &inp, "Number of walls, type, atom types, densities and box-z scale factor for Ewald");
ir->nwall = get_eint(&inp, "nwall", 0, wi);
- ir->wall_type = get_eeenum(&inp, "wall-type", ewt_names, wi);
+ ir->wall_type = get_eeenum(&inp, "wall-type", ewt_names, wi);
ir->wall_r_linpot = get_ereal(&inp, "wall-r-linpot", -1, wi);
setStringEntry(&inp, "wall-atomtype", is->wall_atomtype, nullptr);
- setStringEntry(&inp, "wall-density", is->wall_density, nullptr);
+ setStringEntry(&inp, "wall-density", is->wall_density, nullptr);
ir->wall_ewald_zfac = get_ereal(&inp, "wall-ewald-zfac", 3, wi);
/* COM pulling */
/* Refinement */
printStringNewline(&inp, "NMR refinement stuff");
printStringNoNewline(&inp, "Distance restraints type: No, Simple or Ensemble");
- ir->eDisre = get_eeenum(&inp, "disre", edisre_names, wi);
- printStringNoNewline(&inp, "Force weighting of pairs in one distance restraint: Conservative or Equal");
+ ir->eDisre = get_eeenum(&inp, "disre", edisre_names, wi);
+ printStringNoNewline(
+ &inp, "Force weighting of pairs in one distance restraint: Conservative or Equal");
ir->eDisreWeighting = get_eeenum(&inp, "disre-weighting", edisreweighting_names, wi);
printStringNoNewline(&inp, "Use sqrt of the time averaged times the instantaneous violation");
ir->bDisreMixed = (get_eeenum(&inp, "disre-mixed", yesno_names, wi) != 0);
- ir->dr_fc = get_ereal(&inp, "disre-fc", 1000.0, wi);
+ ir->dr_fc = get_ereal(&inp, "disre-fc", 1000.0, wi);
ir->dr_tau = get_ereal(&inp, "disre-tau", 0.0, wi);
printStringNoNewline(&inp, "Output frequency for pair distances to energy file");
ir->nstdisreout = get_eint(&inp, "nstdisreout", 100, wi);
printStringNoNewline(&inp, "Orientation restraints: No or Yes");
- opts->bOrire = (get_eeenum(&inp, "orire", yesno_names, wi) != 0);
+ opts->bOrire = (get_eeenum(&inp, "orire", yesno_names, wi) != 0);
printStringNoNewline(&inp, "Orientation restraints force constant and tau for time averaging");
- ir->orires_fc = get_ereal(&inp, "orire-fc", 0.0, wi);
+ ir->orires_fc = get_ereal(&inp, "orire-fc", 0.0, wi);
ir->orires_tau = get_ereal(&inp, "orire-tau", 0.0, wi);
- setStringEntry(&inp, "orire-fitgrp", is->orirefitgrp, nullptr);
+ setStringEntry(&inp, "orire-fitgrp", is->orirefitgrp, nullptr);
printStringNoNewline(&inp, "Output frequency for trace(SD) and S to energy file");
ir->nstorireout = get_eint(&inp, "nstorireout", 100, wi);
/* free energy variables */
printStringNewline(&inp, "Free energy variables");
ir->efep = get_eeenum(&inp, "free-energy", efep_names, wi);
- setStringEntry(&inp, "couple-moltype", is->couple_moltype, nullptr);
+ setStringEntry(&inp, "couple-moltype", is->couple_moltype, nullptr);
opts->couple_lam0 = get_eeenum(&inp, "couple-lambda0", couple_lam, wi);
opts->couple_lam1 = get_eeenum(&inp, "couple-lambda1", couple_lam, wi);
opts->bCoupleIntra = (get_eeenum(&inp, "couple-intramol", yesno_names, wi) != 0);
- fep->init_lambda = get_ereal(&inp, "init-lambda", -1, wi); /* start with -1 so
- we can recognize if
- it was not entered */
+ fep->init_lambda = get_ereal(&inp, "init-lambda", -1, wi); /* start with -1 so
+ we can recognize if
+ it was not entered */
fep->init_fep_state = get_eint(&inp, "init-lambda-state", -1, wi);
fep->delta_lambda = get_ereal(&inp, "delta-lambda", 0.0, wi);
fep->nstdhdl = get_eint(&inp, "nstdhdl", 50, wi);
/* Non-equilibrium MD stuff */
printStringNewline(&inp, "Non-equilibrium MD stuff");
- setStringEntry(&inp, "acc-grps", is->accgrps, nullptr);
- setStringEntry(&inp, "accelerate", is->acc, nullptr);
- setStringEntry(&inp, "freezegrps", is->freeze, nullptr);
- setStringEntry(&inp, "freezedim", is->frdim, nullptr);
+ setStringEntry(&inp, "acc-grps", is->accgrps, nullptr);
+ setStringEntry(&inp, "accelerate", is->acc, nullptr);
+ setStringEntry(&inp, "freezegrps", is->freeze, nullptr);
+ setStringEntry(&inp, "freezedim", is->frdim, nullptr);
ir->cos_accel = get_ereal(&inp, "cos-acceleration", 0, wi);
- setStringEntry(&inp, "deform", is->deform, nullptr);
+ setStringEntry(&inp, "deform", is->deform, nullptr);
/* simulated tempering variables */
printStringNewline(&inp, "simulated tempering variables");
- ir->bSimTemp = (get_eeenum(&inp, "simulated-tempering", yesno_names, wi) != 0);
+ ir->bSimTemp = (get_eeenum(&inp, "simulated-tempering", yesno_names, wi) != 0);
ir->simtempvals->eSimTempScale = get_eeenum(&inp, "simulated-tempering-scaling", esimtemp_names, wi);
- ir->simtempvals->simtemp_low = get_ereal(&inp, "sim-temp-low", 300.0, wi);
- ir->simtempvals->simtemp_high = get_ereal(&inp, "sim-temp-high", 300.0, wi);
+ ir->simtempvals->simtemp_low = get_ereal(&inp, "sim-temp-low", 300.0, wi);
+ ir->simtempvals->simtemp_high = get_ereal(&inp, "sim-temp-high", 300.0, wi);
/* expanded ensemble variables */
if (ir->efep == efepEXPANDED || ir->bSimTemp)
{
gmx::KeyValueTreeObject convertedValues = flatKeyValueTreeFromInpFile(inp);
gmx::KeyValueTreeTransformer transform;
- transform.rules()->addRule()
- .keyMatchType("/", gmx::StringCompareType::CaseAndDashInsensitive);
+ transform.rules()->addRule().keyMatchType("/", gmx::StringCompareType::CaseAndDashInsensitive);
mdModules->initMdpTransform(transform.rules());
- for (const auto &path : transform.mappedPaths())
+ for (const auto& path : transform.mappedPaths())
{
GMX_ASSERT(path.size() == 1, "Inconsistent mapping back to mdp options");
mark_einp_set(inp, path[0].c_str());
}
- MdpErrorHandler errorHandler(wi);
- auto result
- = transform.transform(convertedValues, &errorHandler);
- ir->params = new gmx::KeyValueTreeObject(result.object());
+ MdpErrorHandler errorHandler(wi);
+ auto result = transform.transform(convertedValues, &errorHandler);
+ ir->params = new gmx::KeyValueTreeObject(result.object());
mdModules->adjustInputrecBasedOnModules(ir);
errorHandler.setBackMapping(result.backMapping());
mdModules->assignOptionsToModules(*ir->params, &errorHandler);
}
/* Ion/water position swapping ("computational electrophysiology") */
- printStringNewline(&inp, "Ion/water position swapping for computational electrophysiology setups");
+ printStringNewline(&inp,
+ "Ion/water position swapping for computational electrophysiology setups");
printStringNoNewline(&inp, "Swap positions along direction: no, X, Y, Z");
ir->eSwapCoords = get_eeenum(&inp, "swapcoords", eSwapTypes_names, wi);
if (ir->eSwapCoords != eswapNO)
{
snew(ir->swap->grp[i].molname, STRLEN);
}
- printStringNoNewline(&inp, "Two index groups that contain the compartment-partitioning atoms");
+ printStringNoNewline(&inp,
+ "Two index groups that contain the compartment-partitioning atoms");
setStringEntry(&inp, "split-group0", ir->swap->grp[eGrpSplit0].molname, nullptr);
setStringEntry(&inp, "split-group1", ir->swap->grp[eGrpSplit1].molname, nullptr);
- printStringNoNewline(&inp, "Use center of mass of split groups (yes/no), otherwise center of geometry is used");
+ printStringNoNewline(&inp,
+ "Use center of mass of split groups (yes/no), otherwise center of "
+ "geometry is used");
ir->swap->massw_split[0] = (get_eeenum(&inp, "massw-split0", yesno_names, wi) != 0);
ir->swap->massw_split[1] = (get_eeenum(&inp, "massw-split1", yesno_names, wi) != 0);
printStringNoNewline(&inp, "Name of solvent molecules");
setStringEntry(&inp, "solvent-group", ir->swap->grp[eGrpSolvent].molname, nullptr);
- printStringNoNewline(&inp, "Split cylinder: radius, upper and lower extension (nm) (this will define the channels)");
- printStringNoNewline(&inp, "Note that the split cylinder settings do not have an influence on the swapping protocol,");
- printStringNoNewline(&inp, "however, if correctly defined, the permeation events are recorded per channel");
+ printStringNoNewline(&inp,
+ "Split cylinder: radius, upper and lower extension (nm) (this will "
+ "define the channels)");
+ printStringNoNewline(&inp,
+ "Note that the split cylinder settings do not have an influence on "
+ "the swapping protocol,");
+ printStringNoNewline(
+ &inp,
+ "however, if correctly defined, the permeation events are recorded per channel");
ir->swap->cyl0r = get_ereal(&inp, "cyl0-r", 2.0, wi);
ir->swap->cyl0u = get_ereal(&inp, "cyl0-up", 1.0, wi);
ir->swap->cyl0l = get_ereal(&inp, "cyl0-down", 1.0, wi);
ir->swap->cyl1u = get_ereal(&inp, "cyl1-up", 1.0, wi);
ir->swap->cyl1l = get_ereal(&inp, "cyl1-down", 1.0, wi);
- printStringNoNewline(&inp, "Average the number of ions per compartment over these many swap attempt steps");
+ printStringNoNewline(
+ &inp,
+ "Average the number of ions per compartment over these many swap attempt steps");
ir->swap->nAverage = get_eint(&inp, "coupl-steps", 10, wi);
- printStringNoNewline(&inp, "Names of the ion types that can be exchanged with solvent molecules,");
- printStringNoNewline(&inp, "and the requested number of ions of this type in compartments A and B");
+ printStringNoNewline(
+ &inp, "Names of the ion types that can be exchanged with solvent molecules,");
+ printStringNoNewline(
+ &inp, "and the requested number of ions of this type in compartments A and B");
printStringNoNewline(&inp, "-1 means fix the numbers as found in step 0");
for (i = 0; i < nIonTypes; i++)
{
ir->swap->grp[ig].nmolReq[1] = get_eint(&inp, buf, -1, wi);
}
- printStringNoNewline(&inp, "By default (i.e. bulk offset = 0.0), ion/water exchanges happen between layers");
- printStringNoNewline(&inp, "at maximum distance (= bulk concentration) to the split group layers. However,");
- printStringNoNewline(&inp, "an offset b (-1.0 < b < +1.0) can be specified to offset the bulk layer from the middle at 0.0");
- printStringNoNewline(&inp, "towards one of the compartment-partitioning layers (at +/- 1.0).");
+ printStringNoNewline(
+ &inp,
+ "By default (i.e. bulk offset = 0.0), ion/water exchanges happen between layers");
+ printStringNoNewline(
+ &inp,
+ "at maximum distance (= bulk concentration) to the split group layers. However,");
+ printStringNoNewline(&inp,
+ "an offset b (-1.0 < b < +1.0) can be specified to offset the bulk "
+ "layer from the middle at 0.0");
+ printStringNoNewline(&inp,
+ "towards one of the compartment-partitioning layers (at +/- 1.0).");
ir->swap->bulkOffset[0] = get_ereal(&inp, "bulk-offsetA", 0.0, wi);
ir->swap->bulkOffset[1] = get_ereal(&inp, "bulk-offsetB", 0.0, wi);
if (!(ir->swap->bulkOffset[0] > -1.0 && ir->swap->bulkOffset[0] < 1.0)
- || !(ir->swap->bulkOffset[1] > -1.0 && ir->swap->bulkOffset[1] < 1.0) )
+ || !(ir->swap->bulkOffset[1] > -1.0 && ir->swap->bulkOffset[1] < 1.0))
{
warning_error(wi, "Bulk layer offsets must be > -1.0 and < 1.0 !");
}
- printStringNoNewline(&inp, "Start to swap ions if threshold difference to requested count is reached");
+ printStringNoNewline(
+ &inp, "Start to swap ions if threshold difference to requested count is reached");
ir->swap->threshold = get_ereal(&inp, "threshold", 1.0, wi);
}
/* User defined thingies */
printStringNewline(&inp, "User defined thingies");
- setStringEntry(&inp, "user1-grps", is->user1, nullptr);
- setStringEntry(&inp, "user2-grps", is->user2, nullptr);
- ir->userint1 = get_eint(&inp, "userint1", 0, wi);
- ir->userint2 = get_eint(&inp, "userint2", 0, wi);
- ir->userint3 = get_eint(&inp, "userint3", 0, wi);
- ir->userint4 = get_eint(&inp, "userint4", 0, wi);
- ir->userreal1 = get_ereal(&inp, "userreal1", 0, wi);
- ir->userreal2 = get_ereal(&inp, "userreal2", 0, wi);
- ir->userreal3 = get_ereal(&inp, "userreal3", 0, wi);
- ir->userreal4 = get_ereal(&inp, "userreal4", 0, wi);
+ setStringEntry(&inp, "user1-grps", is->user1, nullptr);
+ setStringEntry(&inp, "user2-grps", is->user2, nullptr);
+ ir->userint1 = get_eint(&inp, "userint1", 0, wi);
+ ir->userint2 = get_eint(&inp, "userint2", 0, wi);
+ ir->userint3 = get_eint(&inp, "userint3", 0, wi);
+ ir->userint4 = get_eint(&inp, "userint4", 0, wi);
+ ir->userreal1 = get_ereal(&inp, "userreal1", 0, wi);
+ ir->userreal2 = get_ereal(&inp, "userreal2", 0, wi);
+ ir->userreal3 = get_ereal(&inp, "userreal3", 0, wi);
+ ir->userreal4 = get_ereal(&inp, "userreal4", 0, wi);
#undef CTYPE
{
/* Process options if necessary */
for (m = 0; m < 2; m++)
{
- for (i = 0; i < 2*DIM; i++)
+ for (i = 0; i < 2 * DIM; i++)
{
dumdub[m][i] = 0.0;
}
case epctISOTROPIC:
if (sscanf(dumstr[m], "%lf", &(dumdub[m][XX])) != 1)
{
- warning_error(wi, "Pressure coupling incorrect number of values (I need exactly 1)");
+ warning_error(
+ wi,
+ "Pressure coupling incorrect number of values (I need exactly 1)");
}
dumdub[m][YY] = dumdub[m][ZZ] = dumdub[m][XX];
break;
case epctSURFACETENSION:
if (sscanf(dumstr[m], "%lf%lf", &(dumdub[m][XX]), &(dumdub[m][ZZ])) != 2)
{
- warning_error(wi, "Pressure coupling incorrect number of values (I need exactly 2)");
+ warning_error(
+ wi,
+ "Pressure coupling incorrect number of values (I need exactly 2)");
}
dumdub[m][YY] = dumdub[m][XX];
break;
case epctANISOTROPIC:
- if (sscanf(dumstr[m], "%lf%lf%lf%lf%lf%lf",
- &(dumdub[m][XX]), &(dumdub[m][YY]), &(dumdub[m][ZZ]),
- &(dumdub[m][3]), &(dumdub[m][4]), &(dumdub[m][5])) != 6)
+ if (sscanf(dumstr[m], "%lf%lf%lf%lf%lf%lf", &(dumdub[m][XX]), &(dumdub[m][YY]),
+ &(dumdub[m][ZZ]), &(dumdub[m][3]), &(dumdub[m][4]), &(dumdub[m][5]))
+ != 6)
{
- warning_error(wi, "Pressure coupling incorrect number of values (I need exactly 6)");
+ warning_error(
+ wi,
+ "Pressure coupling incorrect number of values (I need exactly 6)");
}
break;
default:
ir->ref_p[YY][ZZ] = dumdub[1][5];
if (ir->ref_p[XX][YY] != 0 && ir->ref_p[XX][ZZ] != 0 && ir->ref_p[YY][ZZ] != 0)
{
- warning(wi, "All off-diagonal reference pressures are non-zero. Are you sure you want to apply a threefold shear stress?\n");
+ warning(wi,
+ "All off-diagonal reference pressures are non-zero. Are you sure you want to "
+ "apply a threefold shear stress?\n");
}
ir->compress[XX][YY] = dumdub[0][3];
ir->compress[XX][ZZ] = dumdub[0][4];
{
warning(wi, "The lambda=0 and lambda=1 states for coupling are identical");
}
- if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE ||
- opts->couple_lam1 == ecouplamNONE))
+ if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE || opts->couple_lam1 == ecouplamNONE))
{
- warning(wi, "For proper sampling of the (nearly) decoupled state, stochastic dynamics should be used");
+ warning(wi,
+ "For proper sampling of the (nearly) decoupled state, stochastic dynamics "
+ "should be used");
}
}
else
{
- warning_note(wi, "Free energy is turned off, so we will not decouple the molecule listed in your input.");
+ 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 (fep->edHdLPrintEnergy == edHdLPrintEnergyYES)
{
fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
- warning_note(wi, "Old option for dhdl-print-energy given: "
+ warning_note(wi,
+ "Old option for dhdl-print-energy given: "
"changing \"yes\" to \"total\"\n");
}
* If the (advanced) user does FEP through manual topology changes,
* this check will not be triggered.
*/
- if (ir->efep != efepNO && ir->fepvals->n_lambda == 0 &&
- ir->fepvals->sc_alpha != 0 &&
- (couple_lambda_has_vdw_on(opts->couple_lam0) &&
- couple_lambda_has_vdw_on(opts->couple_lam1)))
+ if (ir->efep != efepNO && ir->fepvals->n_lambda == 0 && ir->fepvals->sc_alpha != 0
+ && (couple_lambda_has_vdw_on(opts->couple_lam0) && couple_lambda_has_vdw_on(opts->couple_lam1)))
{
- warning(wi, "You are using soft-core interactions while the Van der Waals interactions are not decoupled (note that the sc-coul option is only active when using lambda states). Although this will not lead to errors, you will need much more sampling than without soft-core interactions. Consider using sc-alpha=0.");
+ warning(wi,
+ "You are using soft-core interactions while the Van der Waals interactions are "
+ "not decoupled (note that the sc-coul option is only active when using lambda "
+ "states). Although this will not lead to errors, you will need much more "
+ "sampling than without soft-core interactions. Consider using sc-alpha=0.");
}
}
else
}
double gmx_unused canary;
- int ndeform = sscanf(is->deform, "%lf %lf %lf %lf %lf %lf %lf",
- &(dumdub[0][0]), &(dumdub[0][1]), &(dumdub[0][2]),
- &(dumdub[0][3]), &(dumdub[0][4]), &(dumdub[0][5]), &canary);
+ int ndeform = sscanf(is->deform, "%lf %lf %lf %lf %lf %lf %lf", &(dumdub[0][0]), &(dumdub[0][1]),
+ &(dumdub[0][2]), &(dumdub[0][3]), &(dumdub[0][4]), &(dumdub[0][5]), &canary);
if (strlen(is->deform) > 0 && ndeform != 6)
{
- warning_error(wi, gmx::formatString("Cannot parse exactly 6 box deformation velocities from string '%s'", is->deform).c_str());
+ warning_error(
+ wi, gmx::formatString(
+ "Cannot parse exactly 6 box deformation velocities from string '%s'", is->deform)
+ .c_str());
}
for (i = 0; i < 3; i++)
{
{
if (ir->compress[m][j] != 0)
{
- sprintf(warn_buf, "An off-diagonal box element has deform set while compressibility > 0 for the same component of another box vector, this might lead to spurious periodicity effects.");
+ sprintf(warn_buf,
+ "An off-diagonal box element has deform set while "
+ "compressibility > 0 for the same component of another box "
+ "vector, this might lead to spurious periodicity effects.");
warning(wi, warn_buf);
}
}
sfree(dumstr[1]);
}
-static int search_QMstring(const 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;
/* TODO this is utility functionality (search for the index of a
string in a collection), so should be refactored and located more
centrally. */
-int search_string(const char *s, int ng, char *gn[])
+int search_string(const char* s, int ng, char* gn[])
{
int i;
s);
}
-static bool do_numbering(int natoms, SimulationGroups *groups,
+static bool do_numbering(int natoms,
+ SimulationGroups* groups,
gmx::ArrayRef<std::string> groupsFromMdpFile,
- t_blocka *block, char *gnames[],
- SimulationAtomGroupType gtype, int restnm,
- int grptp, bool bVerbose,
- warninp_t wi)
+ t_blocka* block,
+ char* gnames[],
+ SimulationAtomGroupType gtype,
+ int restnm,
+ int grptp,
+ bool bVerbose,
+ warninp_t wi)
{
- unsigned short *cbuf;
- AtomGroupIndices *grps = &(groups->groups[gtype]);
- int j, gid, aj, ognr, ntot = 0;
- const char *title;
- bool bRest;
- char warn_buf[STRLEN];
+ unsigned short* cbuf;
+ AtomGroupIndices* grps = &(groups->groups[gtype]);
+ int j, gid, aj, ognr, ntot = 0;
+ const char* title;
+ bool bRest;
+ char warn_buf[STRLEN];
title = shortName(gtype);
}
/* Now go over the atoms in the group */
- for (j = block->index[gid]; (j < block->index[gid+1]); j++)
+ for (j = block->index[gid]; (j < block->index[gid + 1]); j++)
{
aj = block->a[j];
ognr = cbuf[aj];
if (ognr != NOGID)
{
- gmx_fatal(FARGS, "Atom %d in multiple %s groups (%d and %d)",
- aj+1, title, ognr+1, i+1);
+ gmx_fatal(FARGS, "Atom %d in multiple %s groups (%d and %d)", aj + 1, title,
+ ognr + 1, i + 1);
}
else
{
{
if (grptp == egrptpALL)
{
- gmx_fatal(FARGS, "%d atoms are not part of any of the %s groups",
- natoms-ntot, title);
+ gmx_fatal(FARGS, "%d atoms are not part of any of the %s groups", natoms - ntot, title);
}
else if (grptp == egrptpPART)
{
- sprintf(warn_buf, "%d atoms are not part of any of the %s groups",
- natoms-ntot, title);
+ sprintf(warn_buf, "%d atoms are not part of any of the %s groups", natoms - ntot, title);
warning_note(wi, warn_buf);
}
/* Assign all atoms currently unassigned to a rest group */
{
if (bVerbose)
{
- fprintf(stderr,
- "Making dummy/rest group for %s containing %d elements\n",
- title, natoms-ntot);
+ fprintf(stderr, "Making dummy/rest group for %s containing %d elements\n", title,
+ natoms - ntot);
}
/* Add group name "rest" */
grps->emplace_back(restnm);
return (bRest && grptp == egrptpPART);
}
-static void calc_nrdf(const gmx_mtop_t *mtop, t_inputrec *ir, char **gnames)
+static void calc_nrdf(const gmx_mtop_t* mtop, t_inputrec* ir, char** gnames)
{
- t_grpopts *opts;
- pull_params_t *pull;
- int natoms, imin, jmin;
- int *nrdf2, *na_vcm, na_tot;
- double *nrdf_tc, *nrdf_vcm, nrdf_uc, *nrdf_vcm_sub;
- ivec *dof_vcm;
- int as;
+ t_grpopts* opts;
+ pull_params_t* pull;
+ int natoms, imin, jmin;
+ int * nrdf2, *na_vcm, na_tot;
+ double * nrdf_tc, *nrdf_vcm, nrdf_uc, *nrdf_vcm_sub;
+ ivec* dof_vcm;
+ int as;
/* Calculate nrdf.
* First calc 3xnr-atoms for each group
opts = &ir->opts;
- const SimulationGroups &groups = mtop->groups;
- natoms = mtop->natoms;
+ const SimulationGroups& groups = mtop->groups;
+ natoms = mtop->natoms;
/* Allocate one more for a possible rest group */
/* We need to sum degrees of freedom into doubles,
* since floats give too low nrdf's above 3 million atoms.
*/
- snew(nrdf_tc, groups.groups[SimulationAtomGroupType::TemperatureCoupling].size()+1);
- snew(nrdf_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size()+1);
- snew(dof_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size()+1);
- snew(na_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size()+1);
- snew(nrdf_vcm_sub, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size()+1);
+ snew(nrdf_tc, groups.groups[SimulationAtomGroupType::TemperatureCoupling].size() + 1);
+ snew(nrdf_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size() + 1);
+ snew(dof_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size() + 1);
+ snew(na_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size() + 1);
+ snew(nrdf_vcm_sub, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size() + 1);
for (gmx::index i = 0; i < gmx::ssize(groups.groups[SimulationAtomGroupType::TemperatureCoupling]); i++)
{
nrdf_tc[i] = 0;
}
- for (gmx::index i = 0; i < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval])+1; i++)
+ for (gmx::index i = 0;
+ i < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval]) + 1; i++)
{
- nrdf_vcm[i] = 0;
+ nrdf_vcm[i] = 0;
clear_ivec(dof_vcm[i]);
na_vcm[i] = 0;
nrdf_vcm_sub[i] = 0;
snew(nrdf2, natoms);
for (const AtomProxy atomP : AtomRange(*mtop))
{
- const t_atom &local = atomP.atom();
+ const t_atom& local = atomP.atom();
int i = atomP.globalAtomNumber();
- nrdf2[i] = 0;
+ nrdf2[i] = 0;
if (local.ptype == eptAtom || local.ptype == eptNucleus)
{
int g = getGroupType(groups, SimulationAtomGroupType::Freeze, i);
if (opts->nFreeze[g][d] == 0)
{
/* Add one DOF for particle i (counted as 2*1) */
- nrdf2[i] += 2;
+ nrdf2[i] += 2;
/* VCM group i has dim d as a DOF */
- dof_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, i)][d] = 1;
+ dof_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, i)][d] =
+ 1;
}
}
- nrdf_tc [getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, i)] += 0.5*nrdf2[i];
- nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, i)] += 0.5*nrdf2[i];
+ nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, i)] +=
+ 0.5 * nrdf2[i];
+ nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, i)] +=
+ 0.5 * nrdf2[i];
}
}
as = 0;
- for (const gmx_molblock_t &molb : mtop->molblock)
+ for (const gmx_molblock_t& molb : mtop->molblock)
{
- const gmx_moltype_t &molt = mtop->moltype[molb.type];
- const t_atom *atom = molt.atoms.atom;
+ const gmx_moltype_t& molt = mtop->moltype[molb.type];
+ const t_atom* atom = molt.atoms.atom;
for (int mol = 0; mol < molb.nmol; mol++)
{
for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
{
gmx::ArrayRef<const int> ia = molt.ilist[ftype].iatoms;
- for (int i = 0; i < molt.ilist[ftype].size(); )
+ for (int i = 0; i < molt.ilist[ftype].size();)
{
/* Subtract degrees of freedom for the constraints,
* if the particles still have degrees of freedom left.
*/
int ai = as + ia[i + 1];
int aj = as + ia[i + 2];
- if (((atom[ia[i + 1]].ptype == eptNucleus) ||
- (atom[ia[i + 1]].ptype == eptAtom)) &&
- ((atom[ia[i + 2]].ptype == eptNucleus) ||
- (atom[ia[i + 2]].ptype == eptAtom)))
+ if (((atom[ia[i + 1]].ptype == eptNucleus) || (atom[ia[i + 1]].ptype == eptAtom))
+ && ((atom[ia[i + 2]].ptype == eptNucleus) || (atom[ia[i + 2]].ptype == eptAtom)))
{
if (nrdf2[ai] > 0)
{
{
imin = 2;
}
- imin = std::min(imin, nrdf2[ai]);
- jmin = std::min(jmin, nrdf2[aj]);
+ imin = std::min(imin, nrdf2[ai]);
+ jmin = std::min(jmin, nrdf2[aj]);
nrdf2[ai] -= imin;
nrdf2[aj] -= jmin;
- nrdf_tc [getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -= 0.5*imin;
- nrdf_tc [getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, aj)] -= 0.5*jmin;
- nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -= 0.5*imin;
- nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, aj)] -= 0.5*jmin;
+ nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -=
+ 0.5 * imin;
+ nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, aj)] -=
+ 0.5 * jmin;
+ nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -=
+ 0.5 * imin;
+ nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, aj)] -=
+ 0.5 * jmin;
}
- i += interaction_function[ftype].nratoms+1;
+ i += interaction_function[ftype].nratoms + 1;
}
}
gmx::ArrayRef<const int> ia = molt.ilist[F_SETTLE].iatoms;
- for (int i = 0; i < molt.ilist[F_SETTLE].size(); )
+ for (int i = 0; i < molt.ilist[F_SETTLE].size();)
{
/* Subtract 1 dof from every atom in the SETTLE */
for (int j = 0; j < 3; j++)
{
- int ai = as + ia[i + 1 + j];
- imin = std::min(2, nrdf2[ai]);
+ int ai = as + ia[i + 1 + j];
+ imin = std::min(2, nrdf2[ai]);
nrdf2[ai] -= imin;
- nrdf_tc [getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -= 0.5*imin;
- nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -= 0.5*imin;
+ nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -=
+ 0.5 * imin;
+ nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -=
+ 0.5 * imin;
}
- i += 4;
+ i += 4;
}
as += molt.atoms.nr;
}
for (int j = 0; j < 2; j++)
{
- const t_pull_group *pgrp;
+ const t_pull_group* pgrp;
pgrp = &pull->group[pull->coord[i].group[j]];
{
/* Subtract 1/2 dof from each group */
int ai = pgrp->ind[0];
- nrdf_tc [getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -= 0.5*imin;
- nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -= 0.5*imin;
+ nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -=
+ 0.5 * imin;
+ nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -=
+ 0.5 * imin;
if (nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, 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.groups[SimulationAtomGroupType::TemperatureCoupling][getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)]]);
+ 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.groups[SimulationAtomGroupType::TemperatureCoupling][getGroupType(
+ groups, SimulationAtomGroupType::TemperatureCoupling, ai)]]);
}
}
else
* the number of degrees of freedom in each vcm group when COM
* translation is removed and 6 when rotation is removed as well.
*/
- for (gmx::index j = 0; j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval])+1; j++)
+ for (gmx::index j = 0;
+ j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval]) + 1; j++)
{
switch (ir->comm_mode)
{
}
}
break;
- case ecmANGULAR:
- nrdf_vcm_sub[j] = 6;
- break;
- default:
- gmx_incons("Checking comm_mode");
+ case ecmANGULAR: nrdf_vcm_sub[j] = 6; break;
+ default: gmx_incons("Checking comm_mode");
}
}
- for (gmx::index i = 0; i < gmx::ssize(groups.groups[SimulationAtomGroupType::TemperatureCoupling]); i++)
+ for (gmx::index i = 0;
+ i < gmx::ssize(groups.groups[SimulationAtomGroupType::TemperatureCoupling]); i++)
{
/* Count the number of atoms of TC group i for every VCM group */
- for (gmx::index j = 0; j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval])+1; j++)
+ for (gmx::index j = 0;
+ j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval]) + 1; j++)
{
na_vcm[j] = 0;
}
*/
nrdf_uc = nrdf_tc[i];
nrdf_tc[i] = 0;
- for (gmx::index j = 0; j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval])+1; j++)
+ for (gmx::index j = 0;
+ j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval]) + 1; j++)
{
if (nrdf_vcm[j] > nrdf_vcm_sub[j])
{
- nrdf_tc[i] += nrdf_uc*(static_cast<double>(na_vcm[j])/static_cast<double>(na_tot))*
- (nrdf_vcm[j] - nrdf_vcm_sub[j])/nrdf_vcm[j];
+ nrdf_tc[i] += nrdf_uc * (static_cast<double>(na_vcm[j]) / static_cast<double>(na_tot))
+ * (nrdf_vcm[j] - nrdf_vcm_sub[j]) / nrdf_vcm[j];
}
}
}
{
opts->nrdf[i] = 0;
}
- fprintf(stderr,
- "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
+ fprintf(stderr, "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
gnames[groups.groups[SimulationAtomGroupType::TemperatureCoupling][i]], opts->nrdf[i]);
}
sfree(nrdf_vcm_sub);
}
-static bool do_egp_flag(t_inputrec *ir, SimulationGroups *groups,
- const char *option, const char *val, int flag)
+static bool do_egp_flag(t_inputrec* ir, SimulationGroups* groups, const char* option, const char* val, int flag)
{
/* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
* But since this is much larger than STRLEN, such a line can not be parsed.
* The real maximum is the number of names that fit in a string: STRLEN/2.
*/
-#define EGP_MAX (STRLEN/2)
- int j, k, nr;
- bool bSet;
+#define EGP_MAX (STRLEN / 2)
+ int j, k, nr;
+ bool bSet;
- auto names = gmx::splitString(val);
+ auto names = gmx::splitString(val);
if (names.size() % 2 != 0)
{
gmx_fatal(FARGS, "The number of groups for %s is odd", option);
{
// TODO this needs to be replaced by a solution using std::find_if
j = 0;
- while ((j < nr) &&
- gmx_strcasecmp(names[2*i].c_str(), *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][j]])))
+ while ((j < nr)
+ && gmx_strcasecmp(
+ names[2 * i].c_str(),
+ *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][j]])))
{
j++;
}
if (j == nr)
{
- gmx_fatal(FARGS, "%s in %s is not an energy group\n",
- names[2*i].c_str(), option);
+ gmx_fatal(FARGS, "%s in %s is not an energy group\n", names[2 * i].c_str(), option);
}
k = 0;
- while ((k < nr) &&
- gmx_strcasecmp(names[2*i+1].c_str(), *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][k]])))
+ while ((k < nr)
+ && gmx_strcasecmp(
+ names[2 * i + 1].c_str(),
+ *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][k]])))
{
k++;
}
if (k == nr)
{
- gmx_fatal(FARGS, "%s in %s is not an energy group\n",
- names[2*i+1].c_str(), option);
+ gmx_fatal(FARGS, "%s in %s is not an energy group\n", names[2 * i + 1].c_str(), option);
}
if ((j < nr) && (k < nr))
{
- ir->opts.egp_flags[nr*j+k] |= flag;
- ir->opts.egp_flags[nr*k+j] |= flag;
+ ir->opts.egp_flags[nr * j + k] |= flag;
+ ir->opts.egp_flags[nr * k + j] |= flag;
bSet = TRUE;
}
}
}
-static void make_swap_groups(
- t_swapcoords *swap,
- t_blocka *grps,
- char **gnames)
+static void make_swap_groups(t_swapcoords* swap, t_blocka* grps, char** gnames)
{
int ig = -1, i = 0, gind;
- t_swapGroup *swapg;
+ t_swapGroup* swapg;
/* Just a quick check here, more thorough checks are in mdrun */
{
swapg = &swap->grp[ig];
gind = search_string(swap->grp[ig].molname, grps->nr, gnames);
- swapg->nat = grps->index[gind+1] - grps->index[gind];
+ swapg->nat = grps->index[gind + 1] - grps->index[gind];
if (swapg->nat > 0)
{
fprintf(stderr, "%s group '%s' contains %d atoms.\n",
- ig < 3 ? eSwapFixedGrp_names[ig] : "Swap",
- swap->grp[ig].molname, swapg->nat);
+ ig < 3 ? eSwapFixedGrp_names[ig] : "Swap", swap->grp[ig].molname, swapg->nat);
snew(swapg->ind, swapg->nat);
for (i = 0; i < swapg->nat; i++)
{
- swapg->ind[i] = grps->a[grps->index[gind]+i];
+ swapg->ind[i] = grps->a[grps->index[gind] + i];
}
}
else
}
-static void make_IMD_group(t_IMD *IMDgroup, char *IMDgname, t_blocka *grps, char **gnames)
+static void make_IMD_group(t_IMD* IMDgroup, char* IMDgname, t_blocka* grps, char** gnames)
{
- int ig, i;
+ int ig, i;
ig = search_string(IMDgname, grps->nr, gnames);
- IMDgroup->nat = grps->index[ig+1] - grps->index[ig];
+ 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",
+ 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];
+ IMDgroup->ind[i] = grps->a[grps->index[ig] + i];
}
}
}
-void do_index(const char* mdparin, const char *ndx,
- gmx_mtop_t *mtop,
- bool bVerbose,
- const gmx::MdModulesNotifier ¬ifier,
- t_inputrec *ir,
- warninp_t wi)
+void do_index(const char* mdparin,
+ const char* ndx,
+ gmx_mtop_t* mtop,
+ bool bVerbose,
+ const gmx::MdModulesNotifier& notifier,
+ t_inputrec* ir,
+ warninp_t wi)
{
- t_blocka *defaultIndexGroups;
+ t_blocka* defaultIndexGroups;
int natoms;
- t_symtab *symtab;
+ t_symtab* symtab;
t_atoms atoms_all;
char warnbuf[STRLEN], **gnames;
int nr;
defaultIndexGroups = init_index(ndx, &gnames);
}
- SimulationGroups *groups = &mtop->groups;
- natoms = mtop->natoms;
- symtab = &mtop->symtab;
+ SimulationGroups* groups = &mtop->groups;
+ natoms = mtop->natoms;
+ symtab = &mtop->symtab;
for (int i = 0; (i < defaultIndexGroups->nr); i++)
{
groups->groupNames.emplace_back(put_symtab(symtab, gnames[i]));
}
groups->groupNames.emplace_back(put_symtab(symtab, "rest"));
- restnm = groups->groupNames.size() - 1;
+ restnm = groups->groupNames.size() - 1;
GMX_RELEASE_ASSERT(restnm == defaultIndexGroups->nr, "Size of allocations must match");
- srenew(gnames, defaultIndexGroups->nr+1);
- gnames[restnm] = *(groups->groupNames.back());
+ srenew(gnames, defaultIndexGroups->nr + 1);
+ gnames[restnm] = *(groups->groupNames.back());
set_warning_line(wi, mdparin, -1);
auto temperatureCouplingTauValues = gmx::splitString(is->tau_t);
auto temperatureCouplingReferenceValues = gmx::splitString(is->ref_t);
auto temperatureCouplingGroupNames = gmx::splitString(is->tcgrps);
- if (temperatureCouplingTauValues.size() != temperatureCouplingGroupNames.size() ||
- temperatureCouplingReferenceValues.size() != temperatureCouplingGroupNames.size())
+ if (temperatureCouplingTauValues.size() != temperatureCouplingGroupNames.size()
+ || temperatureCouplingReferenceValues.size() != temperatureCouplingGroupNames.size())
{
- gmx_fatal(FARGS, "Invalid T coupling input: %zu groups, %zu ref-t values and "
+ gmx_fatal(FARGS,
+ "Invalid T coupling input: %zu groups, %zu ref-t values and "
"%zu tau-t values",
- temperatureCouplingGroupNames.size(),
- temperatureCouplingReferenceValues.size(),
+ temperatureCouplingGroupNames.size(), temperatureCouplingReferenceValues.size(),
temperatureCouplingTauValues.size());
}
const bool useReferenceTemperature = integratorHasReferenceTemperature(ir);
do_numbering(natoms, groups, temperatureCouplingGroupNames, defaultIndexGroups, gnames,
- SimulationAtomGroupType::TemperatureCoupling,
- restnm, useReferenceTemperature ? egrptpALL : egrptpALL_GENREST, bVerbose, wi);
+ SimulationAtomGroupType::TemperatureCoupling, restnm,
+ useReferenceTemperature ? egrptpALL : egrptpALL_GENREST, bVerbose, wi);
nr = groups->groups[SimulationAtomGroupType::TemperatureCoupling].size();
ir->opts.ngtc = nr;
snew(ir->opts.nrdf, nr);
if (ir->etc != etcVRESCALE && ir->opts.tau_t[i] == 0)
{
- warning_note(wi, "tau-t = -1 is the value to signal that a group should not have temperature coupling. Treating your use of tau-t = 0 as if you used -1.");
+ warning_note(
+ wi,
+ "tau-t = -1 is the value to signal that a group should not have "
+ "temperature coupling. Treating your use of tau-t = 0 as if you used -1.");
}
if (ir->opts.tau_t[i] >= 0)
{
if ((ir->etc == etcNOSEHOOVER) && (ir->epc == epcBERENDSEN))
{
- gmx_fatal(FARGS, "Cannot do Nose-Hoover temperature with Berendsen pressure control with md-vv; use either vrescale temperature with berendsen pressure or Nose-Hoover temperature with MTTK pressure");
+ gmx_fatal(FARGS,
+ "Cannot do Nose-Hoover temperature with Berendsen pressure control with "
+ "md-vv; use either vrescale temperature with berendsen pressure or "
+ "Nose-Hoover temperature with MTTK pressure");
}
if (ir->epc == epcMTTK)
{
if (ir->etc != etcNOSEHOOVER)
{
- gmx_fatal(FARGS, "Cannot do MTTK pressure coupling without Nose-Hoover temperature control");
+ gmx_fatal(FARGS,
+ "Cannot do MTTK pressure coupling without Nose-Hoover temperature "
+ "control");
}
else
{
if (ir->nstpcouple != ir->nsttcouple)
{
- int mincouple = std::min(ir->nstpcouple, ir->nsttcouple);
+ int mincouple = std::min(ir->nstpcouple, ir->nsttcouple);
ir->nstpcouple = ir->nsttcouple = mincouple;
- sprintf(warn_buf, "for current Trotter decomposition methods with vv, nsttcouple and nstpcouple must be equal. Both have been reset to min(nsttcouple,nstpcouple) = %d", mincouple);
+ sprintf(warn_buf,
+ "for current Trotter decomposition methods with vv, nsttcouple and "
+ "nstpcouple must be equal. Both have been reset to "
+ "min(nsttcouple,nstpcouple) = %d",
+ mincouple);
warning_note(wi, warn_buf);
}
}
if (ir->nsttcouple != 1)
{
ir->nsttcouple = 1;
- sprintf(warn_buf, "Andersen temperature control methods assume nsttcouple = 1; there is no need for larger nsttcouple > 1, since no global parameters are computed. nsttcouple has been reset to 1");
+ sprintf(warn_buf,
+ "Andersen temperature control methods assume nsttcouple = 1; there is no "
+ "need for larger nsttcouple > 1, since no global parameters are computed. "
+ "nsttcouple has been reset to 1");
warning_note(wi, warn_buf);
}
}
nstcmin = tcouple_min_integration_steps(ir->etc);
if (nstcmin > 1)
{
- if (tau_min/(ir->delta_t*ir->nsttcouple) < nstcmin - 10*GMX_REAL_EPS)
+ 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),
- tau_min, nstcmin,
- ir->nsttcouple*ir->delta_t);
+ 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), tau_min, nstcmin, ir->nsttcouple * ir->delta_t);
warning(wi, warn_buf);
}
}
/* Simulated annealing for each group. There are nr groups */
auto simulatedAnnealingGroupNames = gmx::splitString(is->anneal);
- if (simulatedAnnealingGroupNames.size() == 1 &&
- gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[0], "N", 1))
+ if (simulatedAnnealingGroupNames.size() == 1
+ && gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[0], "N", 1))
{
simulatedAnnealingGroupNames.resize(0);
}
- if (!simulatedAnnealingGroupNames.empty() &&
- gmx::ssize(simulatedAnnealingGroupNames) != nr)
+ if (!simulatedAnnealingGroupNames.empty() && gmx::ssize(simulatedAnnealingGroupNames) != nr)
{
gmx_fatal(FARGS, "Wrong number of annealing values: %zu (for %d groups)\n",
simulatedAnnealingGroupNames.size(), nr);
{
if (ir->opts.anneal_npoints[i] == 1)
{
- gmx_fatal(FARGS, "Please specify at least a start and an end point for annealing\n");
+ gmx_fatal(
+ FARGS,
+ "Please specify at least a start and an end point for annealing\n");
}
snew(ir->opts.anneal_time[i], ir->opts.anneal_npoints[i]);
snew(ir->opts.anneal_temp[i], ir->opts.anneal_npoints[i]);
std::vector<real> allSimulatedAnnealingTimes(numSimulatedAnnealingFields);
std::vector<real> allSimulatedAnnealingTemperatures(numSimulatedAnnealingFields);
- convertReals(wi, simulatedAnnealingTimes, "anneal-time", allSimulatedAnnealingTimes.data());
- convertReals(wi, simulatedAnnealingTemperatures, "anneal-temp", allSimulatedAnnealingTemperatures.data());
+ convertReals(wi, simulatedAnnealingTimes, "anneal-time",
+ allSimulatedAnnealingTimes.data());
+ convertReals(wi, simulatedAnnealingTemperatures, "anneal-temp",
+ allSimulatedAnnealingTemperatures.data());
for (i = 0, k = 0; i < nr; i++)
{
for (j = 0; j < ir->opts.anneal_npoints[i]; j++)
ir->opts.anneal_temp[i][j] = allSimulatedAnnealingTemperatures[k];
if (j == 0)
{
- if (ir->opts.anneal_time[i][0] > (ir->init_t+GMX_REAL_EPS))
+ if (ir->opts.anneal_time[i][0] > (ir->init_t + GMX_REAL_EPS))
{
gmx_fatal(FARGS, "First time point for annealing > init_t.\n");
}
else
{
/* j>0 */
- if (ir->opts.anneal_time[i][j] < ir->opts.anneal_time[i][j-1])
+ if (ir->opts.anneal_time[i][j] < ir->opts.anneal_time[i][j - 1])
{
- gmx_fatal(FARGS, "Annealing timepoints out of order: t=%f comes after t=%f\n",
- ir->opts.anneal_time[i][j], ir->opts.anneal_time[i][j-1]);
+ gmx_fatal(FARGS,
+ "Annealing timepoints out of order: t=%f comes after "
+ "t=%f\n",
+ ir->opts.anneal_time[i][j], ir->opts.anneal_time[i][j - 1]);
}
}
if (ir->opts.anneal_temp[i][j] < 0)
{
- gmx_fatal(FARGS, "Found negative temperature in annealing: %f\n", ir->opts.anneal_temp[i][j]);
+ gmx_fatal(FARGS, "Found negative temperature in annealing: %f\n",
+ ir->opts.anneal_temp[i][j]);
}
k++;
}
ir->opts.anneal_npoints[i]);
fprintf(stderr, "Time (ps) Temperature (K)\n");
/* All terms except the last one */
- for (j = 0; j < (ir->opts.anneal_npoints[i]-1); j++)
+ for (j = 0; j < (ir->opts.anneal_npoints[i] - 1); j++)
{
- fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
+ fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j],
+ ir->opts.anneal_temp[i][j]);
}
/* Finally the last one */
- j = ir->opts.anneal_npoints[i]-1;
+ j = ir->opts.anneal_npoints[i] - 1;
if (ir->opts.annealing[i] == eannSINGLE)
{
- fprintf(stderr, "%9.1f- %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
+ fprintf(stderr, "%9.1f- %5.1f\n", ir->opts.anneal_time[i][j],
+ ir->opts.anneal_temp[i][j]);
}
else
{
- fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
- if (std::fabs(ir->opts.anneal_temp[i][j]-ir->opts.anneal_temp[i][0]) > GMX_REAL_EPS)
+ fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j],
+ ir->opts.anneal_temp[i][j]);
+ if (std::fabs(ir->opts.anneal_temp[i][j] - ir->opts.anneal_temp[i][0]) > GMX_REAL_EPS)
{
- warning_note(wi, "There is a temperature jump when your annealing loops back.\n");
+ warning_note(wi,
+ "There is a temperature jump when your annealing "
+ "loops back.\n");
}
}
}
accelerationGroupNames.size(), accelerations.size());
}
do_numbering(natoms, groups, accelerationGroupNames, defaultIndexGroups, gnames,
- SimulationAtomGroupType::Acceleration,
- restnm, egrptpALL_GENREST, bVerbose, wi);
+ SimulationAtomGroupType::Acceleration, restnm, egrptpALL_GENREST, bVerbose, wi);
nr = groups->groups[SimulationAtomGroupType::Acceleration].size();
snew(ir->opts.acc, nr);
ir->opts.ngacc = nr;
freezeGroupNames.size(), freezeDims.size());
}
do_numbering(natoms, groups, freezeGroupNames, defaultIndexGroups, gnames,
- SimulationAtomGroupType::Freeze,
- restnm, egrptpALL_GENREST, bVerbose, wi);
+ SimulationAtomGroupType::Freeze, restnm, egrptpALL_GENREST, bVerbose, wi);
nr = groups->groups[SimulationAtomGroupType::Freeze].size();
ir->opts.ngfrz = nr;
snew(ir->opts.nFreeze, nr);
{
if (!gmx::equalCaseInsensitive(freezeDims[k], "N", 1))
{
- sprintf(warnbuf, "Please use Y(ES) or N(O) for freezedim only "
- "(not %s)", freezeDims[k].c_str());
+ sprintf(warnbuf,
+ "Please use Y(ES) or N(O) for freezedim only "
+ "(not %s)",
+ freezeDims[k].c_str());
warning(wi, warn_buf);
}
}
auto energyGroupNames = gmx::splitString(is->energy);
do_numbering(natoms, groups, energyGroupNames, defaultIndexGroups, gnames,
- SimulationAtomGroupType::EnergyOutput,
- restnm, egrptpALL_GENREST, bVerbose, wi);
+ SimulationAtomGroupType::EnergyOutput, restnm, egrptpALL_GENREST, bVerbose, wi);
add_wall_energrps(groups, ir->nwall, symtab);
- ir->opts.ngener = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
+ ir->opts.ngener = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
auto vcmGroupNames = gmx::splitString(is->vcm);
- bRest =
- do_numbering(natoms, groups, vcmGroupNames, defaultIndexGroups, gnames,
- SimulationAtomGroupType::MassCenterVelocityRemoval,
- restnm, vcmGroupNames.empty() ? egrptpALL_GENREST : egrptpPART, bVerbose, wi);
+ bRest = do_numbering(natoms, groups, vcmGroupNames, defaultIndexGroups, gnames,
+ SimulationAtomGroupType::MassCenterVelocityRemoval, restnm,
+ vcmGroupNames.empty() ? egrptpALL_GENREST : egrptpPART, bVerbose, wi);
if (bRest)
{
- warning(wi, "Some atoms are not part of any center of mass motion removal group.\n"
+ warning(wi,
+ "Some atoms are not part of any center of mass motion removal group.\n"
"This may lead to artifacts.\n"
"In most cases one should use one group for the whole system.");
}
auto user1GroupNames = gmx::splitString(is->user1);
do_numbering(natoms, groups, user1GroupNames, defaultIndexGroups, gnames,
- SimulationAtomGroupType::User1,
- restnm, egrptpALL_GENREST, bVerbose, wi);
+ SimulationAtomGroupType::User1, restnm, egrptpALL_GENREST, bVerbose, wi);
auto user2GroupNames = gmx::splitString(is->user2);
do_numbering(natoms, groups, user2GroupNames, defaultIndexGroups, gnames,
- SimulationAtomGroupType::User2,
- restnm, egrptpALL_GENREST, bVerbose, wi);
+ SimulationAtomGroupType::User2, restnm, egrptpALL_GENREST, bVerbose, wi);
auto compressedXGroupNames = gmx::splitString(is->x_compressed_groups);
do_numbering(natoms, groups, compressedXGroupNames, defaultIndexGroups, gnames,
- SimulationAtomGroupType::CompressedPositionOutput,
- restnm, egrptpONE, bVerbose, wi);
+ SimulationAtomGroupType::CompressedPositionOutput, restnm, egrptpONE, bVerbose, wi);
auto orirefFitGroupNames = gmx::splitString(is->orirefitgrp);
do_numbering(natoms, groups, orirefFitGroupNames, defaultIndexGroups, gnames,
- SimulationAtomGroupType::OrientationRestraintsFit,
- restnm, egrptpALL_GENREST, bVerbose, wi);
+ SimulationAtomGroupType::OrientationRestraintsFit, restnm, egrptpALL_GENREST,
+ bVerbose, wi);
/* QMMM input processing */
auto qmGroupNames = gmx::splitString(is->QMMM);
auto qmBasisSets = gmx::splitString(is->QMbasis);
if (ir->eI != eiMimic)
{
- if (qmMethods.size() != qmGroupNames.size() ||
- qmBasisSets.size() != qmGroupNames.size())
+ if (qmMethods.size() != qmGroupNames.size() || qmBasisSets.size() != qmGroupNames.size())
{
- gmx_fatal(FARGS, "Invalid QMMM input: %zu groups %zu basissets"
- " and %zu methods\n", qmGroupNames.size(),
- qmBasisSets.size(), qmMethods.size());
+ gmx_fatal(FARGS,
+ "Invalid QMMM input: %zu groups %zu basissets"
+ " and %zu methods\n",
+ qmGroupNames.size(), qmBasisSets.size(), qmMethods.size());
}
/* group rest, if any, is always MM! */
do_numbering(natoms, groups, qmGroupNames, defaultIndexGroups, gnames,
- SimulationAtomGroupType::QuantumMechanics,
- restnm, egrptpALL_GENREST, bVerbose, wi);
+ SimulationAtomGroupType::QuantumMechanics, restnm, egrptpALL_GENREST, bVerbose, wi);
nr = qmGroupNames.size(); /*atoms->grps[egcQMMM].nr;*/
ir->opts.ngQM = qmGroupNames.size();
snew(ir->opts.QMmethod, nr);
/* input consists of strings: RHF CASSCF PM3 .. These need to be
* converted to the corresponding enum in names.c
*/
- ir->opts.QMmethod[i] = search_QMstring(qmMethods[i].c_str(),
- eQMmethodNR,
- eQMmethod_names);
- ir->opts.QMbasis[i] = search_QMstring(qmBasisSets[i].c_str(),
- eQMbasisNR,
- eQMbasis_names);
-
+ ir->opts.QMmethod[i] = search_QMstring(qmMethods[i].c_str(), eQMmethodNR, eQMmethod_names);
+ ir->opts.QMbasis[i] = search_QMstring(qmBasisSets[i].c_str(), eQMbasisNR, eQMbasis_names);
}
auto qmMultiplicities = gmx::splitString(is->QMmult);
auto qmCharges = gmx::splitString(is->QMcharge);
}
/* group rest, if any, is always MM! */
do_numbering(natoms, groups, qmGroupNames, defaultIndexGroups, gnames,
- SimulationAtomGroupType::QuantumMechanics,
- restnm, egrptpALL_GENREST, bVerbose, wi);
+ SimulationAtomGroupType::QuantumMechanics, restnm, egrptpALL_GENREST, bVerbose, wi);
ir->opts.ngQM = qmGroupNames.size();
}
for (auto group : gmx::keysOf(groups->groups))
{
fprintf(stderr, "%-16s has %zu element(s):", shortName(group), groups->groups[group].size());
- for (const auto &entry : groups->groups[group])
+ for (const auto& entry : groups->groups[group])
{
fprintf(stderr, " %s", *(groups->groupNames[entry]));
}
}
nr = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
- snew(ir->opts.egp_flags, nr*nr);
+ snew(ir->opts.egp_flags, nr * nr);
bExcl = do_egp_flag(ir, groups, "energygrp-excl", is->egpexcl, EGP_EXCL);
if (bExcl && ir->cutoff_scheme == ecutsVERLET)
}
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))
+ 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");
+ gmx_fatal(FARGS,
+ "Can only have energy group pair tables in combination with user tables for VdW "
+ "and/or Coulomb");
}
/* final check before going out of scope if simulated tempering variables
*/
if ((ir->expandedvals->nstexpanded < 0) && ir->bSimTemp)
{
- ir->expandedvals->nstexpanded = 2*static_cast<int>(ir->opts.tau_t[0]/ir->delta_t);
- warning(wi, gmx::formatString("the value for nstexpanded was not specified for "
- " expanded ensemble simulated tempering. It is set to 2*tau_t (%d) "
- "by default, but it is recommended to set it to an explicit value!",
- ir->expandedvals->nstexpanded));
+ ir->expandedvals->nstexpanded = 2 * static_cast<int>(ir->opts.tau_t[0] / ir->delta_t);
+ warning(wi, gmx::formatString(
+ "the value for nstexpanded was not specified for "
+ " expanded ensemble simulated tempering. It is set to 2*tau_t (%d) "
+ "by default, but it is recommended to set it to an explicit value!",
+ ir->expandedvals->nstexpanded));
}
for (i = 0; (i < defaultIndexGroups->nr); i++)
{
sfree(gnames);
done_blocka(defaultIndexGroups);
sfree(defaultIndexGroups);
-
}
-
-static void check_disre(const gmx_mtop_t *mtop)
+static void check_disre(const gmx_mtop_t* mtop)
{
if (gmx_mtop_ftype_count(mtop, F_DISRES) > 0)
{
- const gmx_ffparams_t &ffparams = mtop->ffparams;
+ const gmx_ffparams_t& ffparams = mtop->ffparams;
int ndouble = 0;
int old_label = -1;
for (int i = 0; i < ffparams.numTypes(); i++)
}
if (ndouble > 0)
{
- gmx_fatal(FARGS, "Found %d double distance restraint indices,\n"
+ gmx_fatal(FARGS,
+ "Found %d double distance restraint indices,\n"
"probably the parameters for multiple pairs in one restraint "
- "are not identical\n", ndouble);
+ "are not identical\n",
+ ndouble);
}
}
}
-static bool absolute_reference(t_inputrec *ir, gmx_mtop_t *sys,
- bool posres_only,
- ivec AbsRef)
+static bool absolute_reference(t_inputrec* ir, gmx_mtop_t* sys, bool posres_only, ivec AbsRef)
{
- int d, g, i;
- gmx_mtop_ilistloop_t iloop;
- int nmol;
- t_iparams *pr;
+ int d, g, i;
+ gmx_mtop_ilistloop_t iloop;
+ int nmol;
+ t_iparams* pr;
clear_ivec(AbsRef);
/* Check for position restraints */
iloop = gmx_mtop_ilistloop_init(sys);
- while (const InteractionLists *ilist = gmx_mtop_ilistloop_next(iloop, &nmol))
+ while (const InteractionLists* ilist = gmx_mtop_ilistloop_next(iloop, &nmol))
{
- if (nmol > 0 &&
- (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
+ if (nmol > 0 && (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
{
for (i = 0; i < (*ilist)[F_POSRES].size(); i += 2)
{
{
switch (pr->fbposres.geom)
{
- case efbposresSPHERE:
- AbsRef[XX] = AbsRef[YY] = AbsRef[ZZ] = 1;
- break;
- case efbposresCYLINDERX:
- AbsRef[YY] = AbsRef[ZZ] = 1;
- break;
- case efbposresCYLINDERY:
- AbsRef[XX] = AbsRef[ZZ] = 1;
- break;
+ case efbposresSPHERE: AbsRef[XX] = AbsRef[YY] = AbsRef[ZZ] = 1; break;
+ case efbposresCYLINDERX: AbsRef[YY] = AbsRef[ZZ] = 1; break;
+ case efbposresCYLINDERY: AbsRef[XX] = AbsRef[ZZ] = 1; break;
case efbposresCYLINDER:
/* efbposres is a synonym for efbposresCYLINDERZ for backwards compatibility */
- case efbposresCYLINDERZ:
- AbsRef[XX] = AbsRef[YY] = 1;
- break;
+ case efbposresCYLINDERZ: AbsRef[XX] = AbsRef[YY] = 1; break;
case efbposresX: /* d=XX */
case efbposresY: /* d=YY */
case efbposresZ: /* d=ZZ */
AbsRef[d] = 1;
break;
default:
- gmx_fatal(FARGS, " Invalid geometry for flat-bottom position restraint.\n"
- "Expected nr between 1 and %d. Found %d\n", efbposresNR-1,
- pr->fbposres.geom);
+ gmx_fatal(FARGS,
+ " Invalid geometry for flat-bottom position restraint.\n"
+ "Expected nr between 1 and %d. Found %d\n",
+ efbposresNR - 1, pr->fbposres.geom);
}
}
}
return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
}
-static void
-check_combination_rule_differences(const gmx_mtop_t *mtop, int state,
- bool *bC6ParametersWorkWithGeometricRules,
- bool *bC6ParametersWorkWithLBRules,
- bool *bLBRulesPossible)
+static void check_combination_rule_differences(const gmx_mtop_t* mtop,
+ int state,
+ bool* bC6ParametersWorkWithGeometricRules,
+ bool* bC6ParametersWorkWithLBRules,
+ bool* bLBRulesPossible)
{
- int ntypes, tpi, tpj;
- int *typecount;
- real tol;
- double c6i, c6j, c12i, c12j;
- double c6, c6_geometric, c6_LB;
- double sigmai, sigmaj, epsi, epsj;
- bool bCanDoLBRules, bCanDoGeometricRules;
- const char *ptr;
+ int ntypes, tpi, tpj;
+ int* typecount;
+ real tol;
+ double c6i, c6j, c12i, c12j;
+ double c6, c6_geometric, c6_LB;
+ double sigmai, sigmaj, epsi, epsj;
+ bool bCanDoLBRules, bCanDoGeometricRules;
+ const char* ptr;
/* A tolerance of 1e-5 seems reasonable for (possibly hand-typed)
* force-field floating point parameters.
ptr = getenv("GMX_LJCOMB_TOL");
if (ptr != nullptr)
{
- double dbl;
+ double dbl;
double gmx_unused canary;
if (sscanf(ptr, "%lf%lf", &dbl, &canary) != 1)
{
- gmx_fatal(FARGS, "Could not parse a single floating-point number from GMX_LJCOMB_TOL (%s)", ptr);
+ gmx_fatal(FARGS,
+ "Could not parse a single floating-point number from GMX_LJCOMB_TOL (%s)", ptr);
}
tol = dbl;
}
- *bC6ParametersWorkWithLBRules = TRUE;
- *bC6ParametersWorkWithGeometricRules = TRUE;
- bCanDoLBRules = TRUE;
- ntypes = mtop->ffparams.atnr;
+ *bC6ParametersWorkWithLBRules = TRUE;
+ *bC6ParametersWorkWithGeometricRules = TRUE;
+ bCanDoLBRules = TRUE;
+ ntypes = mtop->ffparams.atnr;
snew(typecount, ntypes);
gmx_mtop_count_atomtypes(mtop, state, typecount);
- *bLBRulesPossible = TRUE;
+ *bLBRulesPossible = TRUE;
for (tpi = 0; tpi < ntypes; ++tpi)
{
c6i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c6;
{
if (!gmx_numzero(c12i) && !gmx_numzero(c12j))
{
- sigmai = gmx::sixthroot(c12i / c6i);
- sigmaj = gmx::sixthroot(c12j / c6j);
- epsi = c6i * c6i /(4.0 * c12i);
- epsj = c6j * c6j /(4.0 * c12j);
- c6_LB = 4.0 * std::sqrt(epsi * epsj) * gmx::power6(0.5 * (sigmai + sigmaj));
+ sigmai = gmx::sixthroot(c12i / c6i);
+ sigmaj = gmx::sixthroot(c12j / c6j);
+ epsi = c6i * c6i / (4.0 * c12i);
+ epsj = c6j * c6j / (4.0 * c12j);
+ c6_LB = 4.0 * std::sqrt(epsi * epsj) * gmx::power6(0.5 * (sigmai + sigmaj));
}
else
{
sfree(typecount);
}
-static void
-check_combination_rules(const t_inputrec *ir, const gmx_mtop_t *mtop,
- warninp_t wi)
+static void check_combination_rules(const t_inputrec* ir, const gmx_mtop_t* mtop, warninp_t wi)
{
bool bLBRulesPossible, bC6ParametersWorkWithGeometricRules, bC6ParametersWorkWithLBRules;
- check_combination_rule_differences(mtop, 0,
- &bC6ParametersWorkWithGeometricRules,
- &bC6ParametersWorkWithLBRules,
- &bLBRulesPossible);
+ check_combination_rule_differences(mtop, 0, &bC6ParametersWorkWithGeometricRules,
+ &bC6ParametersWorkWithLBRules, &bLBRulesPossible);
if (ir->ljpme_combination_rule == eljpmeLB)
{
if (!bC6ParametersWorkWithLBRules || !bLBRulesPossible)
{
- warning(wi, "You are using arithmetic-geometric combination rules "
+ warning(wi,
+ "You are using arithmetic-geometric combination rules "
"in LJ-PME, but your non-bonded C6 parameters do not "
"follow these rules.");
}
{
if (ir->eDispCorr != edispcNO)
{
- warning_note(wi, "You are using geometric combination rules in "
+ 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.");
+ "and/or pressure for isotropic systems, but not forces or surface "
+ "tensions.");
}
else
{
- warning_note(wi, "You are using geometric combination rules in "
+ 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 "
+ "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)
+void triple_check(const char* mdparin, t_inputrec* ir, gmx_mtop_t* sys, warninp_t wi)
{
char err_buf[STRLEN];
int i, m, c, nmol;
bool bCharge, bAcc;
- real *mgrp, mt;
+ real * mgrp, mt;
rvec acc;
gmx_mtop_atomloop_block_t aloopb;
ivec AbsRef;
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)))
+ 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.
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;
+ const real nrdf_at = 2;
real T, tau, max_T_error;
int i;
* 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);
+ 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);
+ 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);
}
}
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);
+ 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 positive using Andersen temperature control, tau_t[%d]=%10.6f",
+ sprintf(err_buf,
+ "all tau_t must be positive 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 = gmx::roundToInt(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);
+ int nsteps = gmx::roundToInt(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 != 0);
}
}
}
- 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) &&
- !ETC_ANDERSEN(ir->etc))
+ 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) && !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");
+ 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");
}
- if (ir->epc == epcPARRINELLORAHMAN &&
- ir->etc == etcNOSEHOOVER)
+ if (ir->epc == epcPARRINELLORAHMAN && ir->etc == etcNOSEHOOVER)
{
real tau_t_max = 0;
for (int g = 0; g < ir->opts.ngtc; g++)
{
tau_t_max = std::max(tau_t_max, ir->opts.tau_t[g]);
}
- if (ir->tau_p < 1.9*tau_t_max)
+ if (ir->tau_p < 1.9 * tau_t_max)
{
- std::string message =
- gmx::formatString("With %s T-coupling and %s p-coupling, "
- "%s (%g) should be at least twice as large as %s (%g) to avoid resonances",
- etcoupl_names[ir->etc],
- epcoupl_names[ir->epc],
- "tau-p", ir->tau_p,
- "tau-t", tau_t_max);
+ std::string message = gmx::formatString(
+ "With %s T-coupling and %s p-coupling, "
+ "%s (%g) should be at least twice as large as %s (%g) to avoid resonances",
+ etcoupl_names[ir->etc], epcoupl_names[ir->epc], "tau-p", ir->tau_p, "tau-t",
+ tau_t_max);
warning(wi, message.c_str());
}
}
{
if (AbsRef[m] && norm2(ir->compress[m]) > 0)
{
- warning(wi, "You are using pressure coupling with absolute position restraints, this will give artifacts. Use the refcoord_scaling option.");
+ warning(wi,
+ "You are using pressure coupling with absolute position restraints, "
+ "this will give artifacts. Use the refcoord_scaling option.");
break;
}
}
bCharge = FALSE;
aloopb = gmx_mtop_atomloop_block_init(sys);
- const t_atom *atom;
+ const t_atom* atom;
while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
{
if (atom->q != 0 || atom->qB != 0)
{
sprintf(err_buf,
"You are using full electrostatics treatment %s for a system without charges.\n"
- "This costs a lot of performance for just processing zeros, consider using %s instead.\n",
+ "This costs a lot of performance for just processing zeros, consider using %s "
+ "instead.\n",
EELTYPE(ir->coulombtype), EELTYPE(eelCUT));
warning(wi, err_buf);
}
/* Generalized reaction field */
if (ir->coulombtype == eelGRF_NOTUSED)
{
- warning_error(wi, "Generalized reaction-field electrostatics is no longer supported. "
- "You can use normal reaction-field instead and compute the reaction-field constant by hand.");
+ warning_error(wi,
+ "Generalized reaction-field electrostatics is no longer supported. "
+ "You can use normal reaction-field instead and compute the reaction-field "
+ "constant by hand.");
}
bAcc = FALSE;
snew(mgrp, sys->groups.groups[SimulationAtomGroupType::Acceleration].size());
for (const AtomProxy atomP : AtomRange(*sys))
{
- const t_atom &local = atomP.atom();
+ const t_atom& local = atomP.atom();
int i = atomP.globalAtomNumber();
mgrp[getGroupType(sys->groups, SimulationAtomGroupType::Acceleration, i)] += local.m;
}
{
for (m = 0; (m < DIM); m++)
{
- acc[m] += ir->opts.acc[i][m]*mgrp[i];
+ acc[m] += ir->opts.acc[i][m] * mgrp[i];
}
mt += mgrp[i];
}
{
if (fabs(acc[m]) > 1e-6)
{
- const char *dim[DIM] = { "X", "Y", "Z" };
- fprintf(stderr,
- "Net Acceleration in %s direction, will %s be corrected\n",
- dim[m], ir->nstcomm != 0 ? "" : "not");
+ const char* dim[DIM] = { "X", "Y", "Z" };
+ fprintf(stderr, "Net Acceleration in %s direction, will %s be corrected\n", dim[m],
+ ir->nstcomm != 0 ? "" : "not");
if (ir->nstcomm != 0 && m < ndof_com(ir))
{
acc[m] /= mt;
- for (i = 0; (i < gmx::ssize(sys->groups.groups[SimulationAtomGroupType::Acceleration])); i++)
+ for (i = 0;
+ (i < gmx::ssize(sys->groups.groups[SimulationAtomGroupType::Acceleration])); i++)
{
ir->opts.acc[i][m] -= acc[m];
}
sfree(mgrp);
}
- if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0 &&
- !gmx_within_tol(sys->ffparams.reppow, 12.0, 10*GMX_DOUBLE_EPS))
+ if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0
+ && !gmx_within_tol(sys->ffparams.reppow, 12.0, 10 * GMX_DOUBLE_EPS))
{
gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
}
bWarned = FALSE;
for (i = 0; i < ir->pull->ncoord && !bWarned; i++)
{
- if (ir->pull->coord[i].group[0] == 0 ||
- ir->pull->coord[i].group[1] == 0)
+ if (ir->pull->coord[i].group[0] == 0 || ir->pull->coord[i].group[1] == 0)
{
absolute_reference(ir, sys, FALSE, AbsRef);
for (m = 0; m < DIM; m++)
{
if (ir->pull->coord[i].dim[m] && !AbsRef[m])
{
- warning(wi, "You are using an absolute reference for pulling, but the rest of the system does not have an absolute reference. This will lead to artifacts.");
+ warning(wi,
+ "You are using an absolute reference for pulling, but the rest of "
+ "the system does not have an absolute reference. This will lead to "
+ "artifacts.");
bWarned = TRUE;
break;
}
{
for (m = 0; m <= i; m++)
{
- if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
- ir->deform[i][m] != 0)
+ if ((ir->epc != epcNO && ir->compress[i][m] != 0) || ir->deform[i][m] != 0)
{
for (c = 0; c < ir->pull->ncoord; c++)
{
- if (ir->pull->coord[c].eGeom == epullgDIRPBC &&
- ir->pull->coord[c].vec[m] != 0)
+ if (ir->pull->coord[c].eGeom == epullgDIRPBC && 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->coord[c].eGeom), 'x'+m);
+ gmx_fatal(FARGS,
+ "Can not have dynamic box while using pull geometry '%s' "
+ "(dim %c)",
+ EPULLGEOM(ir->pull->coord[c].eGeom), 'x' + m);
}
}
}
check_disre(sys);
}
-void double_check(t_inputrec *ir, matrix box,
- bool bHasNormalConstraints,
- bool bHasAnyConstraints,
- warninp_t wi)
+void double_check(t_inputrec* ir, matrix box, bool bHasNormalConstraints, bool bHasAnyConstraints, warninp_t wi)
{
char warn_buf[STRLEN];
- const char *ptr;
+ const char* ptr;
ptr = check_box(ir->ePBC, box);
if (ptr)
{
if (ir->shake_tol <= 0.0)
{
- sprintf(warn_buf, "ERROR: shake-tol must be > 0 instead of %g\n",
- ir->shake_tol);
+ sprintf(warn_buf, "ERROR: shake-tol must be > 0 instead of %g\n", ir->shake_tol);
warning_error(wi, warn_buf);
}
}
- if ( (ir->eConstrAlg == econtLINCS) && bHasNormalConstraints)
+ if ((ir->eConstrAlg == econtLINCS) && bHasNormalConstraints)
{
/* If we have Lincs constraints: */
- if (ir->eI == eiMD && ir->etc == etcNO &&
- ir->eConstrAlg == econtLINCS && ir->nLincsIter == 1)
+ if (ir->eI == eiMD && ir->etc == etcNO && ir->eConstrAlg == econtLINCS && ir->nLincsIter == 1)
{
- sprintf(warn_buf, "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
+ sprintf(warn_buf,
+ "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
warning_note(wi, warn_buf);
}
if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder < 8))
{
- sprintf(warn_buf, "For accurate %s with LINCS constraints, lincs-order should be 8 or more.", ei_names[ir->eI]);
+ sprintf(warn_buf,
+ "For accurate %s with LINCS constraints, lincs-order should be 8 or more.",
+ ei_names[ir->eI]);
warning_note(wi, warn_buf);
}
if (ir->epc == epcMTTK)
{
if (ir->nstlist == 0)
{
- warning(wi, "With nstlist=0 atoms are only put into the box at step 0, therefore drifting atoms might cause the simulation to crash.");
+ warning(wi,
+ "With nstlist=0 atoms are only put into the box at step 0, therefore drifting "
+ "atoms might cause the simulation to crash.");
}
if (gmx::square(ir->rlist) >= max_cutoff2(ir->ePBC, box))
{
- sprintf(warn_buf, "ERROR: The cut-off length is longer than half the shortest box vector or longer than the smallest box diagonal element. Increase the box size or decrease rlist.\n");
+ sprintf(warn_buf,
+ "ERROR: The cut-off length is longer than half the shortest box vector or "
+ "longer than the smallest box diagonal element. Increase the box size or "
+ "decrease rlist.\n");
warning_error(wi, warn_buf);
}
}
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff