hoomd.hpmc.compute

Overview

FreeVolume

Compute the free volume available to a test particle.

SDF

Compute the scale distribution function.

Details

Compute properties of hard particle configurations.

The HPMC compute classes analyze the system configuration and provide results as loggable quantities for use with hoomd.logging.Logger or by direct access via the Python API. FreeVolume computes the free volume available to small particles, such as depletants, and SDF computes the pressure in system of convex particles with a fixed box size.

class hoomd.hpmc.compute.FreeVolume(test_particle_type, num_samples)

Compute the free volume available to a test particle.

Parameters
  • test_particle_type (str) – Test particle type.

  • num_samples (int) – Number of samples to evaluate.

FreeVolume computes the free volume in the simulation state available to a given test particle shape using Monte Carlo integration. Use it in combination with hoomd.hpmc.integrate.HPMCIntegrator, which defines the particle shape parameters. Particles of test_particle_type may or may not be present in the simulation state.

FreeVolume generates num_samples (\(n_\mathrm{samples}\)) trial particle configurations with positions \(\vec{r}^t_j\) uniformly distributed in the simulation box, and orientations \(\mathbf{q}^t_j\) uniformly distributed among rotations matching the box dimensionality. FreeVolume counts the number of successful samples that do not overlap particles in the simulation state:

\[n_\mathrm{success} = \sum_{j=1}^{n_\mathrm{samples}} \prod_{i=0}^{N_\mathrm{particles}-1} \prod_{\vec{A} \in B_\mathrm{images}} \left[ \mathrm{overlap}\left( S_i(\mathbf{q}_i), S_t(\mathbf{q}^t_j, \vec{r}^t_j - (\vec{r}_i + \vec{A})) \right) = \emptyset \right]\]

where \(\mathrm{overlap}\) is the shape overlap function defined in hoomd.hpmc.integrate, \(S_i\) is the shape of particle \(i\), \(S_t\) is the shape of the test particle, \(\vec{A} = h\vec{a}_1 + k\vec{a}_2 + l\vec{a}_3\) is a vector that translates by periodic box images, the set of box images includes all image vectors necessary to find overlaps between particles in the primary image with particles in periodic images, and the square brackets denote the Iverson bracket.

The free volume \(V_\mathrm{free}\) is given by:

\[V_\mathrm{free} = \frac{n_\mathrm{success}} {n_\mathrm{samples}} V_\mathrm{box}\]

where \(V_\mathrm{box}\) is the volume of the simulation box (or area in 2D).

Note

FreeVolume respects the HPMC integrator’s interaction_matrix.

Mixed precision

FreeVolume uses reduced precision floating point arithmetic when checking for particle overlaps in the local particle reference frame.

Box images

On CPU devices, FreeVolume does not apply the minimum image convention. It supports small boxes where particles may overlap with non-primary images of other particles, including self overlap. On GPU devices, FreeVolume applies the minimum image convention.

Examples:

fv = hoomd.hpmc.compute.FreeVolume(test_particle_type='B',
                                   num_samples=1000)
test_particle_type

Test particle type.

Type

str

num_samples

Number of samples to evaluate.

Type

int

property free_volume

Free volume available to the test particle \([\mathrm{length}^{2}]\) in 2D and \([\mathrm{length}^{3}]\) in 3D.

(Loggable: category=”scalar”)

class hoomd.hpmc.compute.SDF(xmax, dx)

Compute the scale distribution function.

Parameters
  • xmax (float) – Maximum x value at the right hand side of the rightmost bin \([\mathrm{length}]\).

  • dx (float) – Bin width \([\mathrm{length}]\).

SDF computes the proability distribution \(s(x)\) of particles overlapping as a function of separation. It estimates \(s(x)\) numerically by computing a histogram with \(\lfloor x_\mathrm{max}/ \delta x \rfloor\) bins of width dx (\(\delta x\)).

See also

Anderson 2016 describes the theory relating SDF to the system pressure.

Implementation

For each pair of particles \(i\) and \(j\) SDF scales the particle separation vector by the factor \((1-x)\) and finds the smallest positive value of \(x\) leading to either an overlap of the particle shapes (a “hard overlap”) or a discontinuous change in the pair energy \(U_{\mathrm{pair},ij}\) (a “soft overlap”):

\[\begin{split}x_{ij}(\vec{A}) = \min \{ & x \in \mathbb{R}_{> 0} : \\ & \mathrm{overlap}\left( S_i(\mathbf{q}_i), S_j(\mathbf{q}_j, (1-x)(\vec{r}_j - (\vec{r}_i + \vec{A}))) \right) \ne \emptyset \\ &\lor \\ & U_{\mathrm{pair},ij}((1-x)(\vec{r}_j - (\vec{r}_i + \vec{A})), \mathbf{q}_i, \mathbf{q}_j) \ne U_{\mathrm{pair},ij}(\vec{r}_j - (\vec{r}_i + \vec{A}), \mathbf{q}_i, \mathbf{q}_j) \\ \} &\end{split}\]

where \(\mathrm{overlap}\) is the shape overlap function defined in hoomd.hpmc.integrate, \(S_i\) is the shape of particle \(i\), and \(\vec{A} = h\vec{a}_1 + k\vec{a}_2 + l\vec{a}_3\) is a vector that translates by periodic box images.

\(x_i\) is the minimum value of \(x_{ij}\) for a single particle:

\[x_i = \min \{ x_{ij} : \vec{A} \in B_\mathrm{images}, j \in [0,N_\mathrm{particles}) \}\]

where the set of box images includes all image vectors necessary to find overlaps between particles in the primary image with particles in periodic images.

SDF adds a single count to the histogram for each particle \(i\), weighted by a factor that is a function of the change in energy upon overlap:

\[s(x + \delta x/2) = \frac{1}{N_\mathrm{particles} \cdot \delta x} \sum_{i=0}^{N_\mathrm{particles}-1} [x \le x_i < x + \delta x] \cdot (1 - \exp(-\beta \Delta U_{i}))\]

where \(\Delta U_{i}\) is the change in energy associated with the first overlap involving particle \(i\) (\(\infty\) for hard overlaps), the square brackets denote the Iverson bracket, and \(s(x + \delta x/2)\) is evaluated for \(\{ x \in \mathbb{R}, 0 \le x < x_\mathrm{max}, x = k \cdot \delta x, k \in \mathbb{Z}^* \}\).

Pressure

The extrapolation of \(s(x)\) to \(x = 0\), \(s(0+)\) is related to the pressure \(P\):

\[\beta P = \rho \left(1 + \frac{s(0+)}{2d} \right)\]

where \(d\) is the dimensionality of the system, \(\rho\) is the number density, and \(\beta = \frac{1}{kT}\). This measurement of the pressure is inherently noisy due to the nature of the sampling. Average betaP over many timesteps to obtain accurate results.

Assuming particle diameters are ~1, these paramater values typically achieve good results:

  • xmax = 0.02

  • dx = 1e-4

In systems near densest packings, dx=1e-5 may be needed along with smaller xmax. Check that \(\sum_k s(x_k) \cdot dx \approx 0.5\).

Warning

SDF only considers negative volume perturbations, and therefore does not compute the correct pressure in simulations where positive volume perturbations may change the system’s potential energy, e.g., systems of concave particles or with non-monotonic enthalpic interactions.

Warning

Because SDF samples pair configurations at discrete separations, the computed pressure is correct only for potentials with constant values and step discontinuities, e.g., square well potentials.

Note

SDF always runs on the CPU.

Mixed precision

SDF uses reduced precision floating point arithmetic when checking for particle overlaps in the local particle reference frame.

Box images

SDF does not apply the minimum image convention. It supports small boxes where particles may overlap with non-primary images of other particles, including self overlap.

xmax

Maximum x value at the right hand side of the rightmost bin \([\mathrm{length}]\).

Type

float

dx

Bin width \([\mathrm{length}]\).

Type

float

property betaP

Beta times pressure in NVT simulations \(\left[ \mathrm{length}^{-d} \right]\).

Uses a polynomial curve fit of degree 5 to estimate \(s(0+)\) and computes the pressure via:

\[\beta P = \rho \left(1 + \frac{s(0+)}{2d} \right)\]

where \(d\) is the dimensionality of the system, \(\rho\) is the number density, and \(\beta = \frac{1}{kT}\).

Attention

In MPI parallel execution, betaP is available on rank 0 only. betaP is None on ranks >= 1.

(Loggable: category=”scalar”)

Type

float

property sdf

\(s[k]\) - The scale distribution function \([\mathrm{probability\ density}]\).

The \(x\) at the center of bin \(k\) is: \(x = k \cdot \delta x + \delta x/2\).

Attention

In MPI parallel execution, the array is available on rank 0 only. sdf is None on ranks >= 1.

(Loggable: category=”sequence”)

Type

(N_bins,) numpy.ndarray of float)