}
}
+/* Returns the mass of atom atomIndex or 1 when setMassesToOne=true */
+static real getMass(const t_atoms &atoms,
+ int atomIndex,
+ bool setMassesToOne)
+{
+ if (!setMassesToOne)
+ {
+ return atoms.atom[atomIndex].m;
+ }
+ else
+ {
+ return 1;
+ }
+}
+
static void get_vsite_masses(const gmx_moltype_t *moltype,
const gmx_ffparams_t *ffparams,
+ bool setMassesToOne,
real *vsite_m,
int *n_nonlin_vsite)
{
{
const t_iparams *ip;
real inv_mass, coeff, m_aj;
- int a1, aj;
+ int a1;
ip = &ffparams->iparams[il->iatoms[i]];
assert(maxj <= 5);
for (j = 1; j < maxj; j++)
{
- cam[j] = moltype->atoms.atom[il->iatoms[i+1+j]].m;
+ int aj = il->iatoms[i + 1 + j];
+ cam[j] = getMass(moltype->atoms, aj, setMassesToOne);
if (cam[j] == 0)
{
- cam[j] = vsite_m[il->iatoms[i+1+j]];
+ cam[j] = vsite_m[aj];
}
if (cam[j] == 0)
{
gmx_fatal(FARGS, "In molecule type '%s' %s construction involves atom %d, which is a virtual site of equal or high complexity. This is not supported.",
*moltype->name,
interaction_function[ft].longname,
- il->iatoms[i+1+j]+1);
+ aj + 1);
}
}
inv_mass = 0;
for (j = 0; j < 3*ffparams->iparams[il->iatoms[i]].vsiten.n; j += 3)
{
- aj = il->iatoms[i+j+2];
- coeff = ffparams->iparams[il->iatoms[i+j]].vsiten.a;
+ int aj = il->iatoms[i + j + 2];
+ coeff = ffparams->iparams[il->iatoms[i+j]].vsiten.a;
if (moltype->atoms.atom[aj].ptype == eptVSite)
{
m_aj = vsite_m[aj];
}
static void get_verlet_buffer_atomtypes(const gmx_mtop_t *mtop,
+ bool setMassesToOne,
verletbuf_atomtype_t **att_p,
int *natt_p,
int *n_nonlin_vsite)
ip = &mtop->ffparams.iparams[il->iatoms[i]];
a1 = il->iatoms[i+1];
a2 = il->iatoms[i+2];
- if (atoms->atom[a2].m > prop[a1].con_mass)
+ real mass1 = getMass(*atoms, a1, setMassesToOne);
+ real mass2 = getMass(*atoms, a2, setMassesToOne);
+ if (mass2 > prop[a1].con_mass)
{
- prop[a1].con_mass = atoms->atom[a2].m;
+ prop[a1].con_mass = mass2;
prop[a1].con_len = ip->constr.dA;
}
- if (atoms->atom[a1].m > prop[a2].con_mass)
+ if (mass1 > prop[a2].con_mass)
{
- prop[a2].con_mass = atoms->atom[a1].m;
+ prop[a2].con_mass = mass1;
prop[a2].con_len = ip->constr.dA;
}
}
* If this is not the case, we overestimate the displacement,
* which leads to a larger buffer (ok since this is an exotic case).
*/
- prop[a1].con_mass = atoms->atom[a2].m;
+ prop[a1].con_mass = getMass(*atoms, a2, setMassesToOne);
prop[a1].con_len = ip->settle.doh;
- prop[a2].con_mass = atoms->atom[a1].m;
+ prop[a2].con_mass = getMass(*atoms, a1, setMassesToOne);
prop[a2].con_len = ip->settle.doh;
- prop[a3].con_mass = atoms->atom[a1].m;
+ prop[a3].con_mass = getMass(*atoms, a1, setMassesToOne);
prop[a3].con_len = ip->settle.doh;
}
get_vsite_masses(&moltype,
&mtop->ffparams,
+ setMassesToOne,
vsite_m,
&n_nonlin_vsite_mol);
if (n_nonlin_vsite != nullptr)
}
else
{
- prop[a].mass = atoms->atom[a].m;
+ prop[a].mass = getMass(*atoms, a, setMassesToOne);
}
prop[a].type = atoms->atom[a].type;
prop[a].q = atoms->atom[a].q;
// Returns an (over)estimate of the energy drift for a single atom pair,
// given the kinetic properties, displacement variances and list buffer.
-static real energyDriftAtomPair(const atom_nonbonded_kinetic_prop_t *prop_i,
- const atom_nonbonded_kinetic_prop_t *prop_j,
+static real energyDriftAtomPair(bool isConstrained_i,
+ bool isConstrained_j,
real s2, real s2i_2d, real s2j_2d,
real r_buffer,
const pot_derivatives_t *der)
else
{
/* For constraints: adapt r and scaling for the Gaussian */
- if (prop_i->bConstr)
+ if (isConstrained_i)
{
real sh, sc;
rsh += sh;
sc_fac *= sc;
}
- if (prop_j->bConstr)
+ if (isConstrained_j)
{
real sh, sc;
lj.d2 = c6*ljDisp->d2 + c12*ljRep->d2;
lj.md3 = c6*ljDisp->md3 + c12*ljRep->md3;
- real pot_lj = energyDriftAtomPair(prop_i, prop_j,
+ real pot_lj = energyDriftAtomPair(prop_i->bConstr, prop_j->bConstr,
s2, s2i_2d, s2j_2d,
rlist - rlj,
&lj);
elec_qq.d2 = elec->d2 *prop_i->q*prop_j->q;
elec_qq.md3 = 0;
- real pot_q = energyDriftAtomPair(prop_i, prop_j,
+ real pot_q = energyDriftAtomPair(prop_i->bConstr, prop_j->bConstr,
s2, s2i_2d, s2j_2d,
rlist - rcoulomb,
&elec_qq);
return md3_pot + md3_sw;
}
+/* Returns the maximum reference temperature over all coupled groups */
+static real maxReferenceTemperature(const t_inputrec &ir)
+{
+ if (EI_MD(ir.eI) && ir.etc == etcNO)
+ {
+ /* This case should be handled outside calc_verlet_buffer_size */
+ gmx_incons("calc_verlet_buffer_size called with an NVE ensemble and reference_temperature < 0");
+ }
+
+ real maxTemperature = 0;
+ for (int i = 0; i < ir.opts.ngtc; i++)
+ {
+ if (ir.opts.tau_t[i] >= 0)
+ {
+ maxTemperature = std::max(maxTemperature, ir.opts.ref_t[i]);
+ }
+ }
+
+ return maxTemperature;
+}
+
+/* Returns the variance of the atomic displacement over timePeriod.
+ *
+ * Note: When not using BD with a non-mass dependendent friction coefficient,
+ * the return value still needs to be divided by the particle mass.
+ */
+static real displacementVariance(const t_inputrec &ir,
+ real temperature,
+ real timePeriod)
+{
+ real kT_fac;
+
+ if (ir.eI == eiBD)
+ {
+ /* Get the displacement distribution from the random component only.
+ * With accurate integration the systematic (force) displacement
+ * should be negligible (unless nstlist is extremely large, which
+ * you wouldn't do anyhow).
+ */
+ kT_fac = 2*BOLTZ*temperature*timePeriod;
+ if (ir.bd_fric > 0)
+ {
+ /* This is directly sigma^2 of the displacement */
+ kT_fac /= ir.bd_fric;
+ }
+ else
+ {
+ /* Per group tau_t is not implemented yet, use the maximum */
+ real tau_t = ir.opts.tau_t[0];
+ for (int i = 1; i < ir.opts.ngtc; i++)
+ {
+ tau_t = std::max(tau_t, ir.opts.tau_t[i]);
+ }
+
+ kT_fac *= tau_t;
+ /* This kT_fac needs to be divided by the mass to get sigma^2 */
+ }
+ }
+ else
+ {
+ kT_fac = BOLTZ*temperature*gmx::square(timePeriod);
+ }
+
+ return kT_fac;
+}
+
+/* Returns the largest sigma of the Gaussian displacement over all particle
+ * types. This ignores constraints, so is an overestimate.
+ */
+static real maxSigma(real kT_fac,
+ int natt,
+ const verletbuf_atomtype_t *att)
+{
+ assert(att);
+ real smallestMass = att[0].prop.mass;
+ for (int i = 1; i < natt; i++)
+ {
+ smallestMass = std::min(smallestMass, att[i].prop.mass);
+ }
+
+ 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,
real nb_clust_frac_pairs_not_in_list_at_cutoff;
verletbuf_atomtype_t *att = nullptr;
- int natt = -1, i;
+ int natt = -1;
real elfac;
- real kT_fac, mass_min;
int ib0, ib1, ib;
real rb, rl;
real drift;
if (reference_temperature < 0)
{
- if (EI_MD(ir->eI) && ir->etc == etcNO)
- {
- /* This case should be handled outside calc_verlet_buffer_size */
- gmx_incons("calc_verlet_buffer_size called with an NVE ensemble and reference_temperature < 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 = 0;
- for (i = 0; i < ir->opts.ngtc; i++)
- {
- if (ir->opts.tau_t[i] >= 0)
- {
- reference_temperature = std::max(reference_temperature,
- ir->opts.ref_t[i]);
- }
- }
+ reference_temperature = maxReferenceTemperature(*ir);
}
/* 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);
- get_verlet_buffer_atomtypes(mtop, &att, &natt, n_nonlin_vsite);
+ /* 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);
+ get_verlet_buffer_atomtypes(mtop, setMassesToOne, &att, &natt, n_nonlin_vsite);
assert(att != nullptr && natt >= 0);
if (debug)
* For inertial dynamics (not Brownian dynamics) the mass factor
* is not included in kT_fac, it is added later.
*/
- if (ir->eI == eiBD)
- {
- /* Get the displacement distribution from the random component only.
- * With accurate integration the systematic (force) displacement
- * should be negligible (unless nstlist is extremely large, which
- * you wouldn't do anyhow).
- */
- kT_fac = 2*BOLTZ*reference_temperature*list_lifetime*ir->delta_t;
- if (ir->bd_fric > 0)
- {
- /* This is directly sigma^2 of the displacement */
- kT_fac /= ir->bd_fric;
-
- /* Set the masses to 1 as kT_fac is the full sigma^2,
- * but we divide by m in ener_drift().
- */
- for (i = 0; i < natt; i++)
- {
- att[i].prop.mass = 1;
- }
- }
- else
- {
- real tau_t;
-
- /* Per group tau_t is not implemented yet, use the maximum */
- tau_t = ir->opts.tau_t[0];
- for (i = 1; i < ir->opts.ngtc; i++)
- {
- tau_t = std::max(tau_t, ir->opts.tau_t[i]);
- }
-
- kT_fac *= tau_t;
- /* This kT_fac needs to be divided by the mass to get sigma^2 */
- }
- }
- else
- {
- kT_fac = BOLTZ*reference_temperature*gmx::square(list_lifetime*ir->delta_t);
- }
-
- mass_min = att[0].prop.mass;
- for (i = 1; i < natt; i++)
- {
- mass_min = std::min(mass_min, att[i].prop.mass);
- }
+ const real kT_fac = displacementVariance(*ir, reference_temperature,
+ list_lifetime*ir->delta_t);
if (debug)
{
fprintf(debug, "LJ rep. -V' %9.2e V'' %9.2e -V''' %9.2e\n", ljRep.md1, ljRep.d2, ljRep.md3);
fprintf(debug, "Electro. -V' %9.2e V'' %9.2e\n", elec.md1, elec.d2);
fprintf(debug, "sqrt(kT_fac) %f\n", std::sqrt(kT_fac));
- fprintf(debug, "mass_min %f\n", mass_min);
}
/* Search using bisection */
ib0 = -1;
/* The drift will be neglible at 5 times the max sigma */
- ib1 = static_cast<int>(5*2*std::sqrt(kT_fac/mass_min)/resolution) + 1;
+ ib1 = static_cast<int>(5*maxSigma(kT_fac, natt, att)/resolution) + 1;
while (ib1 - ib0 > 1)
{
ib = (ib0 + ib1)/2;
*rlist = std::max(ir->rvdw, ir->rcoulomb) + ib1*resolution;
}
+
+/* Returns the pairlist buffer size for use as a minimum buffer size
+ *
+ * Note that this is a rather crude estimate. It is ok for a buffer
+ * set for good energy conservation or RF electrostatics. But it is
+ * too small with PME and the buffer set with the default tolerance.
+ */
+static real minCellSizeFromPairlistBuffer(const t_inputrec &ir)
+{
+ return ir.rlist - std::max(ir.rvdw, ir.rcoulomb);
+}
+
+real minCellSizeForAtomDisplacement(const gmx_mtop_t &mtop,
+ const t_inputrec &ir,
+ real chanceRequested)
+{
+ if (!EI_DYNAMICS(ir.eI) || (EI_MD(ir.eI) && ir.etc == etcNO))
+ {
+ return minCellSizeFromPairlistBuffer(ir);
+ }
+
+ /* 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 displacements.
+ */
+ const real temperature = maxReferenceTemperature(ir);
+
+ const bool setMassesToOne = (ir.eI == eiBD && ir.bd_fric > 0);
+
+ verletbuf_atomtype_t *att = nullptr;
+ int natt = -1;
+ get_verlet_buffer_atomtypes(&mtop, setMassesToOne, &att, &natt, nullptr);
+
+ const real kT_fac = displacementVariance(ir, temperature,
+ ir.nstlist*ir.delta_t);
+
+ /* Resolution of the cell size */
+ real resolution = 0.001;
+
+ /* Search using bisection, avoid 0 and start at 1 */
+ int ib0 = 0;
+ /* The chance will be neglible at 10 times the max sigma */
+ int ib1 = int(10*maxSigma(kT_fac, natt, att)/resolution) + 1;
+ real cellSize = 0;
+ while (ib1 - ib0 > 1)
+ {
+ int ib = (ib0 + ib1)/2;
+ cellSize = ib*resolution;
+
+ /* We assumes atom are distributed uniformly over the cell width.
+ * Once an atom has moved by more than the cellSize (as passed
+ * as the buffer argument to energyDriftAtomPair() below),
+ * the chance of crossing the boundary of the neighbor cell
+ * thus increases as 1/cellSize with the additional displacement
+ * on to of cellSize. We thus create a linear interaction with
+ * derivative = -1/cellSize. Using this in the energyDriftAtomPair
+ * function will return the chance of crossing the next boundary.
+ */
+ const pot_derivatives_t boundaryInteraction = { 1/cellSize, 0, 0 };
+
+ real chance = 0;
+ for (int i = 0; i < natt; i++)
+ {
+ const atom_nonbonded_kinetic_prop_t &propAtom = att[i].prop;
+ real s2_2d;
+ real s2_3d;
+ get_atom_sigma2(kT_fac, &propAtom, &s2_2d, &s2_3d);
+
+ real chancePerAtom = energyDriftAtomPair(propAtom.bConstr, false,
+ s2_2d + s2_3d, s2_2d, 0,
+ cellSize,
+ &boundaryInteraction);
+
+ if (propAtom.bConstr)
+ {
+ /* energyDriftAtomPair() uses an unlimited Gaussian displacement
+ * distribution for constrained atoms, whereas they can
+ * actually not move more than the COM of the two constrained
+ * atoms plus twice the distance from the COM.
+ * Use this maximum, limited displacement when this results in
+ * a smaller chance (note that this is still an overestimate).
+ */
+ real massFraction = propAtom.con_mass/(propAtom.mass + propAtom.con_mass);
+ real comDistance = propAtom.con_len*massFraction;
+
+ real chanceWithMaxDistance =
+ energyDriftAtomPair(false, false,
+ s2_3d, 0, 0,
+ cellSize - 2*comDistance,
+ &boundaryInteraction);
+ chancePerAtom = std::min(chancePerAtom, chanceWithMaxDistance);
+ }
+
+ /* Take into account the line density of the boundary */
+ chancePerAtom /= cellSize;
+
+ chance += att[i].n*chancePerAtom;
+ }
+
+ /* Note: chance is for every nstlist steps */
+ if (chance > chanceRequested*ir.nstlist)
+ {
+ ib0 = ib;
+ }
+ else
+ {
+ ib1 = ib;
+ }
+ }
+
+ sfree(att);
+
+ return cellSize;
+}
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff