hoomd.hpmc.compute#

Overview

FreeVolume

Compute the free volume available to a test particle.

SDF

Compute the scale distribution function via volume perturbations.

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 via volume perturbations.

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 probability distributions \(s_{\mathrm{comp}}(x)\) and \(s_{\mathrm{exp}}(x)\) of particles overlapping as a function of separation for compressive and expansive perturbations, respectively. It estimates \(s_{\mathrm{comp}}(x)\) and \(s_{\mathrm{exp}}(x)\) numerically by computing histograms 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 and its implementation in HOOMD-blue. Eppenga and Frenkel 1984 present a derivation relating the scale distribution function to the system pressure for hard, convex particles. Allen 2006 describes the theory for calculating the pressure in systems with discontinuous potential energy functions. The expansive perturbations are based on theory described in de Miguel and Jackson.

Implementation

SDF constructs two histograms, one for compressive volume perturbations and one for expansive volume perturbations. These histograms represent \(s_{\mathrm{comp}}(x)\) and \(s_{\mathrm{exp}}(x)\). The following discussion applies to compressive volume perturbations and the computation of \(s_{\mathrm{comp}}(x)\); the computation of \(s_{\mathrm{exp}}(x)\) proceeds similarly as noted throughout the description.

For each pair of particles \(i\) and \(j\) SDF scales the particle separation vector by the factor \((1 \pm x)\) (\(+\) for expansive perturbations, \(-\) for compressive perturbations) 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”). For compressive perturbations:

\[\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. For expansive perturbations,

\[\begin{split}x_{ij}(\vec{A}) = \max \{ & 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}\]

For particle \(i\), SDF finds the the minimum (maximum for expansive perturbations) value \(x_i\). For compressive perturbations:

\[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. For expansive perturbations:

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

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

\[s_{\mathrm{comp}}(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_{\mathrm{comp}}(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}^* \}\) for compressive perturbations. For expansive perturbations,

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

where \(s_{\mathrm{exp}}(x - \delta x/2)\) is evaluated for \(\{ x \in \mathbb{R}, -|x_\mathrm{max}| < x \le 0, x = (k - \lfloor x_\mathrm{max} / \delta x \rfloor + 1) \cdot \delta x, k \in \mathbb{Z}^* \}\).

Pressure

The pressure \(P\) is related to the one-sided limits \(s_{\mathrm{comp}}(0+)\) and \(s_{\mathrm{exp}}(0-)\), computed by fitting and extrapolating \(s_{\mathrm{comp}}\) and \(s_{\mathrm{exp}}\) to \(x = 0\).

\[\beta P = \rho \left( 1 + \frac{s_{\mathrm{comp}}(0+) - s_{\mathrm{exp}}(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 parameter 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_\mathrm{comp}(x_k) \cdot dx \approx 0.5\).

Important

SDF samples pair configurations at discrete separations. Therefore, the computed pressure is correct only for potentials with constant values and step discontinuities.

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_\mathrm{comp}(0+)\) (and \(s_\mathrm{exp}(0-)\) if required) and computes the pressure via:

\[\beta P = \rho \left(1 + \frac{s_\mathrm{comp}(0+)}{2d} + \frac{s_\mathrm{exp}(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_compression#

\(s_\mathrm{comp}[k]\) - The scale distribution function for compression moves \([\mathrm{probability\ density}]\).

See also

x_compression defines the bin center locations.

Attention

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

(Loggable: category=”sequence”)

Type:

(N_bins,) numpy.ndarray of float)

property sdf_expansion#

\(s_\mathrm{exp}[k]\) - The scale distribution function for the expansion moves \([\mathrm{probability\ density}]\).

See also

x_expansion defines the bin center locations..

Attention

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

(Loggable: category=”sequence”)

Type:

(N_bins,) numpy.ndarray of float)

property x_compression#

The x values at the center of each bin corresponding to the scale distribution function for the compressive perturbations \([\mathrm{length}]\).

(Loggable: category=”sequence”)

Type:

(N_bins,) numpy.ndarray of float)

property x_expansion#

The x values at the center of each bin corresponding to the scale distribution function for the expansion moves \([\mathrm{length}]\).

(Loggable: category=”sequence”)

Type:

(N_bins,) numpy.ndarray of float)