static void get_vsite_masses(const gmx_moltype_t &moltype,
const gmx_ffparams_t &ffparams,
bool setMassesToOne,
- gmx::ArrayRef<real> vsite_m,
- int *n_nonlin_vsite)
+ gmx::ArrayRef<real> vsite_m)
{
- GMX_RELEASE_ASSERT(n_nonlin_vsite, "Expect a valid pointer");
-
- *n_nonlin_vsite = 0;
+ int numNonlinearVsites = 0;
/* Check for virtual sites, determine mass from constructing atoms */
for (const auto &ilist : extractILists(moltype.ilist, IF_VSITE))
{
vsite_m[a1] = std::min(vsite_m[a1], cam[j]);
}
- (*n_nonlin_vsite)++;
+ numNonlinearVsites++;
break;
}
}
}
}
}
+
+ if (debug && numNonlinearVsites > 0)
+ {
+ fprintf(debug, "The molecule type has %d non-linear virtual constructions\n",
+ numNonlinearVsites);
+ }
}
static std::vector<VerletbufAtomtype>
-get_verlet_buffer_atomtypes(const gmx_mtop_t *mtop,
- bool setMassesToOne,
- int *n_nonlin_vsite)
+getVerletBufferAtomtypes(const gmx_mtop_t &mtop,
+ const bool setMassesToOne)
{
std::vector<VerletbufAtomtype> att;
int ft, i, a1, a2, a3, a;
const t_iparams *ip;
- int n_nonlin_vsite_mol;
-
- if (n_nonlin_vsite != nullptr)
- {
- *n_nonlin_vsite = 0;
- }
- for (const gmx_molblock_t &molblock : mtop->molblock)
+ for (const gmx_molblock_t &molblock : mtop.molblock)
{
int nmol = molblock.nmol;
- const gmx_moltype_t &moltype = mtop->moltype[molblock.type];
+ const gmx_moltype_t &moltype = mtop.moltype[molblock.type];
const t_atoms *atoms = &moltype.atoms;
/* Check for constraints, as they affect the kinetic energy.
for (i = 0; i < il.size(); i += 1+NRAL(ft))
{
- ip = &mtop->ffparams.iparams[il.iatoms[i]];
+ ip = &mtop.ffparams.iparams[il.iatoms[i]];
a1 = il.iatoms[i+1];
a2 = il.iatoms[i+2];
real mass1 = getMass(*atoms, a1, setMassesToOne);
for (i = 0; i < il.size(); i += 1+NRAL(F_SETTLE))
{
- ip = &mtop->ffparams.iparams[il.iatoms[i]];
+ ip = &mtop.ffparams.iparams[il.iatoms[i]];
a1 = il.iatoms[i+1];
a2 = il.iatoms[i+2];
a3 = il.iatoms[i+3];
std::vector<real> vsite_m(atoms->nr);
get_vsite_masses(moltype,
- mtop->ffparams,
+ mtop.ffparams,
setMassesToOne,
- vsite_m,
- &n_nonlin_vsite_mol);
- if (n_nonlin_vsite != nullptr)
- {
- *n_nonlin_vsite += nmol*n_nonlin_vsite_mol;
- }
+ vsite_m);
for (a = 0; a < atoms->nr; a++)
{
return 2*std::sqrt(kT_fac/smallestMass);
}
-void calc_verlet_buffer_size(const gmx_mtop_t *mtop, real boxvol,
- const t_inputrec *ir,
- int nstlist,
- int list_lifetime,
- real reference_temperature,
- const VerletbufListSetup *list_setup,
- int *n_nonlin_vsite,
- real *rlist)
+real
+calcVerletBufferSize(const gmx_mtop_t &mtop,
+ const real boxVolume,
+ const t_inputrec &ir,
+ const int nstlist,
+ const int listLifetime,
+ real referenceTemperature,
+ const VerletbufListSetup &listSetup)
{
double resolution;
char *env;
real rb, rl;
real drift;
- if (!EI_DYNAMICS(ir->eI))
+ if (!EI_DYNAMICS(ir.eI))
{
gmx_incons("Can only determine the Verlet buffer size for integrators that perform dynamics");
}
- if (ir->verletbuf_tol <= 0)
+ if (ir.verletbuf_tol <= 0)
{
gmx_incons("The Verlet buffer tolerance needs to be larger than zero");
}
- if (reference_temperature < 0)
+ if (referenceTemperature < 0)
{
/* We use the maximum temperature with multiple T-coupl groups.
* We could use a per particle temperature, but since particles
* interact, this might underestimate the buffer size.
*/
- reference_temperature = maxReferenceTemperature(*ir);
+ referenceTemperature = maxReferenceTemperature(ir);
- GMX_RELEASE_ASSERT(reference_temperature >= 0, "Without T-coupling we should not end up here");
+ GMX_RELEASE_ASSERT(referenceTemperature >= 0, "Without T-coupling we should not end up here");
}
/* Resolution of the buffer size */
*/
/* Worst case assumption: HCP packing of particles gives largest distance */
- particle_distance = std::cbrt(boxvol*std::sqrt(2)/mtop->natoms);
+ particle_distance = std::cbrt(boxVolume*std::sqrt(2)/mtop.natoms);
/* TODO: Obtain masses through (future) integrator functionality
* to avoid scattering the code with (or forgetting) checks.
*/
- const bool setMassesToOne = (ir->eI == eiBD && ir->bd_fric > 0);
+ const bool setMassesToOne = (ir.eI == eiBD && ir.bd_fric > 0);
const auto att =
- get_verlet_buffer_atomtypes(mtop, setMassesToOne, n_nonlin_vsite);
+ getVerletBufferAtomtypes(mtop, setMassesToOne);
GMX_ASSERT(!att.empty(), "We expect at least one type");
if (debug)
pot_derivatives_t ljDisp = { 0, 0, 0 };
pot_derivatives_t ljRep = { 0, 0, 0 };
- real repPow = mtop->ffparams.reppow;
+ real repPow = mtop.ffparams.reppow;
- if (ir->vdwtype == evdwCUT)
+ if (ir.vdwtype == evdwCUT)
{
real sw_range, md3_pswf;
- switch (ir->vdw_modifier)
+ switch (ir.vdw_modifier)
{
case eintmodNONE:
case eintmodPOTSHIFT:
/* -dV/dr of -r^-6 and r^-reppow */
- ljDisp.md1 = -6*std::pow(ir->rvdw, -7.0);
- ljRep.md1 = repPow*std::pow(ir->rvdw, -(repPow + 1));
+ ljDisp.md1 = -6*std::pow(ir.rvdw, -7.0);
+ ljRep.md1 = repPow*std::pow(ir.rvdw, -(repPow + 1));
/* The contribution of the higher derivatives is negligible */
break;
case eintmodFORCESWITCH:
/* At the cut-off: V=V'=V''=0, so we use only V''' */
- ljDisp.md3 = -md3_force_switch(6.0, ir->rvdw_switch, ir->rvdw);
- ljRep.md3 = md3_force_switch(repPow, ir->rvdw_switch, ir->rvdw);
+ ljDisp.md3 = -md3_force_switch(6.0, ir.rvdw_switch, ir.rvdw);
+ ljRep.md3 = md3_force_switch(repPow, ir.rvdw_switch, ir.rvdw);
break;
case eintmodPOTSWITCH:
/* At the cut-off: V=V'=V''=0.
* V''' is given by the original potential times
* the third derivative of the switch function.
*/
- sw_range = ir->rvdw - ir->rvdw_switch;
+ sw_range = ir.rvdw - ir.rvdw_switch;
md3_pswf = 60.0/gmx::power3(sw_range);
- ljDisp.md3 = -std::pow(ir->rvdw, -6.0 )*md3_pswf;
- ljRep.md3 = std::pow(ir->rvdw, -repPow)*md3_pswf;
+ ljDisp.md3 = -std::pow(ir.rvdw, -6.0 )*md3_pswf;
+ ljRep.md3 = std::pow(ir.rvdw, -repPow)*md3_pswf;
break;
default:
gmx_incons("Unimplemented VdW modifier");
}
}
- else if (EVDW_PME(ir->vdwtype))
+ else if (EVDW_PME(ir.vdwtype))
{
- real b = calc_ewaldcoeff_lj(ir->rvdw, ir->ewald_rtol_lj);
- real r = ir->rvdw;
+ real b = calc_ewaldcoeff_lj(ir.rvdw, ir.ewald_rtol_lj);
+ real r = ir.rvdw;
real br = b*r;
real br2 = br*br;
real br4 = br2*br2;
gmx_fatal(FARGS, "Energy drift calculation is only implemented for plain cut-off Lennard-Jones interactions");
}
- elfac = ONE_4PI_EPS0/ir->epsilon_r;
+ elfac = ONE_4PI_EPS0/ir.epsilon_r;
// Determine the 1st and 2nd derivative for the electostatics
pot_derivatives_t elec = { 0, 0, 0 };
- if (ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype))
+ if (ir.coulombtype == eelCUT || EEL_RF(ir.coulombtype))
{
real eps_rf, k_rf;
- if (ir->coulombtype == eelCUT)
+ if (ir.coulombtype == eelCUT)
{
eps_rf = 1;
k_rf = 0;
}
else
{
- eps_rf = ir->epsilon_rf/ir->epsilon_r;
+ eps_rf = ir.epsilon_rf/ir.epsilon_r;
if (eps_rf != 0)
{
- k_rf = (eps_rf - ir->epsilon_r)/( gmx::power3(ir->rcoulomb) * (2*eps_rf + ir->epsilon_r) );
+ k_rf = (eps_rf - ir.epsilon_r)/( gmx::power3(ir.rcoulomb) * (2*eps_rf + ir.epsilon_r) );
}
else
{
/* epsilon_rf = infinity */
- k_rf = 0.5/gmx::power3(ir->rcoulomb);
+ k_rf = 0.5/gmx::power3(ir.rcoulomb);
}
}
if (eps_rf > 0)
{
- elec.md1 = elfac*(1.0/gmx::square(ir->rcoulomb) - 2*k_rf*ir->rcoulomb);
+ elec.md1 = elfac*(1.0/gmx::square(ir.rcoulomb) - 2*k_rf*ir.rcoulomb);
}
- elec.d2 = elfac*(2.0/gmx::power3(ir->rcoulomb) + 2*k_rf);
+ elec.d2 = elfac*(2.0/gmx::power3(ir.rcoulomb) + 2*k_rf);
}
- else if (EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD)
+ else if (EEL_PME(ir.coulombtype) || ir.coulombtype == eelEWALD)
{
real b, rc, br;
- b = calc_ewaldcoeff_q(ir->rcoulomb, ir->ewald_rtol);
- rc = ir->rcoulomb;
+ b = calc_ewaldcoeff_q(ir.rcoulomb, ir.ewald_rtol);
+ rc = ir.rcoulomb;
br = b*rc;
elec.md1 = elfac*(b*std::exp(-br*br)*M_2_SQRTPI/rc + std::erfc(br)/(rc*rc));
elec.d2 = elfac/(rc*rc)*(2*b*(1 + br*br)*std::exp(-br*br)*M_2_SQRTPI + 2*std::erfc(br)/rc);
* For inertial dynamics (not Brownian dynamics) the mass factor
* is not included in kT_fac, it is added later.
*/
- const real kT_fac = displacementVariance(*ir, reference_temperature,
- list_lifetime*ir->delta_t);
+ const real kT_fac = displacementVariance(ir, referenceTemperature,
+ listLifetime*ir.delta_t);
if (debug)
{
{
ib = (ib0 + ib1)/2;
rb = ib*resolution;
- rl = std::max(ir->rvdw, ir->rcoulomb) + rb;
+ rl = std::max(ir.rvdw, ir.rcoulomb) + rb;
/* Calculate the average energy drift at the last step
* of the nstlist steps at which the pair-list is used.
*/
- drift = energyDrift(att, &mtop->ffparams,
+ drift = energyDrift(att, &mtop.ffparams,
kT_fac,
&ljDisp, &ljRep, &elec,
- ir->rvdw, ir->rcoulomb,
- rl, boxvol);
+ ir.rvdw, ir.rcoulomb,
+ rl, boxVolume);
/* Correct for the fact that we are using a Ni x Nj particle pair list
* and not a 1 x 1 particle pair list. This reduces the drift.
*/
/* We don't have a formula for 8 (yet), use 4 which is conservative */
nb_clust_frac_pairs_not_in_list_at_cutoff =
- surface_frac(std::min(list_setup->cluster_size_i, 4),
+ surface_frac(std::min(listSetup.cluster_size_i, 4),
particle_distance, rl)*
- surface_frac(std::min(list_setup->cluster_size_j, 4),
+ surface_frac(std::min(listSetup.cluster_size_j, 4),
particle_distance, rl);
drift *= nb_clust_frac_pairs_not_in_list_at_cutoff;
/* Convert the drift to drift per unit time per atom */
- drift /= nstlist*ir->delta_t*mtop->natoms;
+ drift /= nstlist*ir.delta_t*mtop.natoms;
if (debug)
{
fprintf(debug, "ib %3d %3d %3d rb %.3f %dx%d fac %.3f drift %.1e\n",
ib0, ib, ib1, rb,
- list_setup->cluster_size_i, list_setup->cluster_size_j,
+ listSetup.cluster_size_i, listSetup.cluster_size_j,
nb_clust_frac_pairs_not_in_list_at_cutoff,
drift);
}
- if (std::abs(drift) > ir->verletbuf_tol)
+ if (std::abs(drift) > ir.verletbuf_tol)
{
ib0 = ib;
}
}
}
- *rlist = std::max(ir->rvdw, ir->rcoulomb) + ib1*resolution;
+ return std::max(ir.rvdw, ir.rcoulomb) + ib1*resolution;
}
/* Returns the pairlist buffer size for use as a minimum buffer size
const bool setMassesToOne = (ir.eI == eiBD && ir.bd_fric > 0);
- const auto atomtypes = get_verlet_buffer_atomtypes(&mtop, setMassesToOne, nullptr);
+ const auto atomtypes = getVerletBufferAtomtypes(mtop, setMassesToOne);
const real kT_fac = displacementVariance(ir, temperature,
ir.nstlist*ir.delta_t);
*/
VerletbufListSetup verletbufGetSafeListSetup(ListSetupType listType);
-/* Calculate the non-bonded pair-list buffer size for the Verlet list
+/* Returns the non-bonded pair-list radius including computed buffer
+ *
+ * Calculate the non-bonded pair-list buffer size for the Verlet list
* based on the particle masses, temperature, LJ types, charges
* and constraints as well as the non-bonded force behavior at the cut-off.
* The pair list update frequency and the list lifetime, which is nstlist-1
* for normal pair-list buffering, are passed separately, as in some cases
* we want an estimate for different values than the ones set in the inputrec.
- * If reference_temperature < 0, the maximum coupling temperature will be used.
+ * If referenceTemperature < 0, the maximum coupling temperature will be used.
* The target is a maximum average energy jump per atom of
- * ir->verletbuf_tol*nstlist*ir->delta_t over the lifetime of the list.
- * Returns the number of non-linear virtual sites. For these it's difficult
- * to determine their contribution to the drift exaclty, so we approximate.
- * Returns the pair-list cut-off.
+ * inputrec.verletbuf_tol*nstlist*inputrec.delta_t over the lifetime of the list.
+ *
+ * \note For non-linear virtual sites it can be problematic to determine their
+ * contribution to the drift exaclty, so we approximate.
+ *
+ * \param[in] mtop The system topology
+ * \param[in] boxVolume The volume of the unit cell
+ * \param[in] inputrec The input record
+ * \param[in] nstlist The pair list update frequency in steps (is not taken from \p inputrec)
+ * \param[in] listLifetime The lifetime of the pair-list, usually nstlist-1, but could be different for dynamic pruning
+ * \param[in] referenceTemperature The reference temperature for the ensemble
+ * \param[in] listSetup The pair-list setup
+ * \returns The computed pair-list radius including buffer
*/
-void calc_verlet_buffer_size(const gmx_mtop_t *mtop, real boxvol,
- const t_inputrec *ir,
- int nstlist,
- int list_lifetime,
- real reference_temperature,
- const VerletbufListSetup *list_setup,
- int *n_nonlin_vsite,
- real *rlist);
+real
+calcVerletBufferSize(const gmx_mtop_t &mtop,
+ real boxVolume,
+ const t_inputrec &inputrec,
+ int nstlist,
+ int listLifetime,
+ real referenceTemperature,
+ const VerletbufListSetup &listSetup);
/* Convenience type */
using PartitioningPerMoltype = gmx::ArrayRef<const gmx::RangePartitioning>;
float listfac_ok, listfac_max;
int nstlist_orig, nstlist_prev;
- real rlistWithReferenceNstlist, rlist_inc, rlist_ok, rlist_max;
+ real rlist_inc, rlist_ok, rlist_max;
real rlist_new, rlist_prev;
size_t nstlist_ind = 0;
gmx_bool bBox, bDD, bCont;
*/
nstlist_prev = ir->nstlist;
ir->nstlist = nbnxnReferenceNstlist;
- calc_verlet_buffer_size(mtop, det(box), ir, ir->nstlist, ir->nstlist - 1,
- -1, &listSetup,
- nullptr, &rlistWithReferenceNstlist);
+ const real rlistWithReferenceNstlist =
+ calcVerletBufferSize(*mtop, det(box), *ir, ir->nstlist, ir->nstlist - 1,
+ -1, listSetup);
ir->nstlist = nstlist_prev;
/* Determine the pair list size increase due to zero interactions */
}
/* Set the pair-list buffer size in ir */
- calc_verlet_buffer_size(mtop, det(box), ir, ir->nstlist, ir->nstlist - 1, -1, &listSetup, nullptr, &rlist_new);
+ rlist_new =
+ calcVerletBufferSize(*mtop, det(box), *ir, ir->nstlist, ir->nstlist - 1, -1, listSetup);
/* Does rlist fit in the box? */
bBox = (gmx::square(rlist_new) < max_cutoff2(ir->ePBC, box));
*/
int listLifetime = tunedNstlistPrune - (useGpuList ? 0 : 1);
listParams->nstlistPrune = tunedNstlistPrune;
- calc_verlet_buffer_size(mtop, det(box), ir,
- tunedNstlistPrune, listLifetime,
- -1, &listSetup, nullptr,
- &listParams->rlistInner);
+ listParams->rlistInner =
+ calcVerletBufferSize(*mtop, det(box), *ir,
+ tunedNstlistPrune, listLifetime,
+ -1, listSetup);
/* On the GPU we apply the dynamic pruning in a rolling fashion
* every c_nbnxnGpuRollingListPruningInterval steps,
}
if (supportsDynamicPairlistGenerationInterval(*ir))
{
- VerletbufListSetup listSetup1x1 = { 1, 1 };
- real rlistOuter;
- real rlistInner;
- calc_verlet_buffer_size(mtop, det(box), ir, ir->nstlist, ir->nstlist - 1,
- -1, &listSetup1x1, nullptr, &rlistOuter);
+ const VerletbufListSetup listSetup1x1 = { 1, 1 };
+ const real rlistOuter =
+ calcVerletBufferSize(*mtop, det(box), *ir, ir->nstlist, ir->nstlist - 1,
+ -1, listSetup1x1);
+ real rlistInner = rlistOuter;
if (listParams->useDynamicPruning)
{
int listLifeTime = listParams->nstlistPrune - (useGpuList ? 0 : 1);
- calc_verlet_buffer_size(mtop, det(box), ir, listParams->nstlistPrune, listLifeTime,
- -1, &listSetup1x1, nullptr, &rlistInner);
+ rlistInner =
+ calcVerletBufferSize(*mtop, det(box), *ir, listParams->nstlistPrune, listLifeTime,
+ -1, listSetup1x1);
}
mesg += gmx::formatString("At tolerance %g kJ/mol/ps per atom, equivalent classical 1x1 list would be:\n",
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff