integrate

Metropolis Monte Carlo

The hard particle Monte Carlo (HPMC) integrator HPMCIntegrator samples equilibrium system states using the Metropolis Monte Carlo method. In this method, HPMCIntegrator takes the existing system state in the configuration \(C = (\vec{r}_0, \vec{r}_1, \ldots \vec{r}_{N_\mathrm{particles}-1}, \mathbf{q}_0, \mathbf{q}_2, \ldots \mathbf{q}_{N_\mathrm{particles}-1})\) with potential energy \(U\) and perturbs it to a trial configuration \(C^t\) with potential energy \(U^t\) leading to an energy difference \(\Delta U = U^t - U\). The trial move is accepted with the probability:

\[\begin{split}p_\mathrm{accept} = \begin{cases} \exp(-\beta \Delta U) & \Delta U > 0 \\ 1 & \Delta U \le 0 \\ \end{cases}\end{split}\]

When the trial move is accepted, the system state is set to the the trial configuration. When it is not accepted, the move is rejected and the state is not modified.

Temperature

HPMCIntegrator.kT (and the related \(\beta = \frac{1}{kT}\)) are used throughout all HPMC operations. Set HPMCIntegrator.kT to control the temperature in systems with finite interactions (see Energy evaluation below). Use the default \(kT = 1\) for systems of purely hard particles where \(\Delta U\) is either 0 or \(\infty\).

Local trial moves

HPMCIntegrator generates local trial moves for a single particle \(i\) at a time. The move is either a translation move or a rotation move, selected randomly with the probability of a translation move set by HPMCIntegrator.translation_move_probability (\(p_\mathrm{translation}\)).

The form of the trial move depends on the dimensionality of the system. Let \(u\) be a random value in the interval \([0,1]\), \(\vec{v}\) be a random vector uniformly distributed within the ball of radius 1, and \(\mathbf{w}\) be a random unit quaternion from the set of uniformly distributed rotations. Then the 3D trial move for particle \(i\) is:

\[\begin{split}\begin{cases} \left( \begin{array}{l} \vec{r}^t_i = \vec{r}_i + d_i \vec{v}, \\ \mathbf{q}^t_i = \mathbf{q}_i \end{array} \right) & u \le p_\mathrm{translation} \\ \left( \begin{array}{l} \vec{r}^t_i = \vec{r}_i, \\ \mathbf{q}^t_i = \frac{\mathbf{q}_i + a_i \mathbf{w}} {\vert \mathbf{q}_i + a_i \mathbf{w} \vert} \end{array} \right) & u > p_\mathrm{translation} \\ \end{cases}\end{split}\]

where \(d_i\) is the translation move size for particle \(i\) (set by particle type with HPMCIntegrator.d) and \(a_i\) is the rotation move size (set by particle type with HPMCIntegrator.a).

In 2D boxes, let \(\vec{v}\) be a random vector uniformly distributed within the disk of radius 1 in the x,y plane and \(\alpha\) be a random angle in radians in the interval \([-a_i,a_i]\). Form a quaternion that rotates about the z axis by \(\alpha\): \(\mathbf{w} = (\cos(\alpha/2), 0, 0, \sin(\alpha/2))\). The 2D trial move for particle \(i\) is:

\[\begin{split}\begin{cases} \left( \begin{array}{l} \vec{r}^t_i = \vec{r}_i + d_i \vec{v}, \\ \mathbf{q}^t_i = \mathbf{q}_i \end{array} \right) & u \le p_\mathrm{translation} \\ \left( \begin{array}{l} \vec{r}^t_i = \vec{r}_i, \\ \mathbf{q}^t_i = \frac{\mathbf{q}_i \cdot \mathbf{w}} {\vert \mathbf{q}_i \cdot \mathbf{w} \vert} \end{array} \right) & u > p_\mathrm{translation} \\ \end{cases}\end{split}\]

Note

For non-orientable spheres, \(p_\mathrm{translation} = 1\).

Timesteps

In the serial CPU implementation, HPMCIntegrator performs nselect trial moves per particle in each timestep (which defaults to 4). To achieve detailed balance at the level of a timestep, HPMCIntegrator randomly chooses with equal probability to loop through particles in forward index or reverse index order (random selection severely degrades performance due to cache incoherency). In the GPU and MPI implementations, trial moves are performed in parallel for particles in active domains while leaving particles on the border fixed (see Anderson 2016 for a full description). As a consequence, a single timestep may perform more or less than nselect trial moves per particle when using the parallel code paths. Monitor the number of trial moves performed with HPMCIntegrator.translate_moves and HPMCIntegrator.rotate_moves.

Random numbers

HPMCIntegrator uses a pseudorandom number stream to generate the trial moves. Set the seed using hoomd.Simulation.seed. Given the same seed, the same initial configuration, and the same execution configuration (device and MPI configuration), HPMCIntegrator, will produce exactly the same trajectory.

Note

Due to limited floating point precision, full trajectory reproducibility is only possible with the same binary installation running on the same hardware device. Compiler optimizations, changes to the HOOMD source code, and machine specific code paths may lead to different trajectories.

Energy evaluation

HPMCIntegrator evaluates the energy of a configuration from a number of terms:

\[U = U_{\mathrm{pair}} + U_{\mathrm{shape}} + U_{\mathrm{external}}\]

To enable simulations of small systems, the pair and shape energies evaluate interactions between pairs of particles in multiple box images:

\[\begin{split}U_{\mathrm{pair}} = & \sum_{i=0}^{N_\mathrm{particles}-1} \sum_{j=i+1}^{N_\mathrm{particles}-1} U_{\mathrm{pair},ij}(\vec{r}_j - \vec{r}_i, \mathbf{q}_i, \mathbf{q}_j) \\ + & \sum_{i=0}^{N_\mathrm{particles}-1} \sum_{j=i}^{N_\mathrm{particles}-1} \sum_{\vec{A} \in B_\mathrm{images}, \vec{A} \ne \vec{0}} U_{\mathrm{pair},ij}(\vec{r}_j - (\vec{r}_i + \vec{A}), \mathbf{q}_i, \mathbf{q}_j)\end{split}\]

where \(\vec{A} = h\vec{a}_1 + k\vec{a}_2 + l\vec{a}_3\) is a vector that translates by periodic box images and the set of box images includes all image vectors necessary to find interactions between particles in the primary image with particles in periodic images The first sum evaluates interactions between particle \(i\) with other particles (not itself) in the primary box image. The second sum evaluates interactions between particle \(i\) and all potentially interacting periodic images of all particles (including itself). HPMCIntegrator computes \(U_{\mathrm{shape}}\) similarly (see below).

External potentials apply to each particle individually:

\[U_\mathrm{external} = \sum_{i=0}^\mathrm{N_particles-1} U_{\mathrm{external},i}(\vec{r}_i, \mathbf{q}_i)\]

Potential classes in hoomd.hpmc.pair evaluate \(U_{\mathrm{pair},ij}\). HPMC sums all the Pair potentials in pair_potentials during integration.

Similarly, potential classes in hoomd.hpmc.external evaluate \(U_{\mathrm{external},i}\). Add instances of these classes to external_potentials to apply during integration.

Shape overlap tests

HPMCIntegrator performs shape overlap tests to evaluate \(U_{\mathrm{shape}}\). Let \(S\) be the set of all points inside the shape in the local coordinate system of the shape:

\[S = \{ \vec{a} \in \mathbb{R}^3 : \vec{a} \enspace \mathrm{inside\ the\ shape} \}\]

See the subclasses of HPMCIntegrator for formal definitions of the shapes, whose parameters are set by particle type. Let \(S_i\) refer specifically to the shape for particle \(i\).

The quaternion \(\mathbf{q}\) represents a rotation of the shape from its local coordinate system to the given orientation:

\[S(\mathbf{q}) = \{ \mathbf{q}\vec{a}\mathbf{q}^* : \vec{a} \in S \}\]

The full transformation from the local shape coordinate system to the simulation box coordinate system includes a rotation and translation:

\[S(\mathbf{q}, \vec{r}) = \{ \mathbf{q}\vec{a}\mathbf{q}^* + \vec{r} : \vec{a} \in S \}\]

HPMCIntegrator defines the shape overlap test for two shapes:

\[\mathrm{overlap}(S_1, S_2) = S_1 \bigcap S_2\]

To check for overlaps between two particles in the box, rotating both shapes from their local frame to the box frame, and translate \(S_2\) relative to particle 1:

\[\mathrm{overlap}(S_1(\mathbf{q}_1), S_2(\mathbf{q}_2, \vec{r}_2 - \vec{r}_1))\]

The complete hard shape interaction energy for a given configuration is:

\[\begin{split}U_\mathrm{shape} = \quad & \infty \cdot \sum_{i=0}^{N_\mathrm{particles}-1} \sum_{j=i+1}^{N_\mathrm{particles}-1} \left[ \mathrm{overlap}\left( S_i(\mathbf{q}_i), S_j(\mathbf{q}_j, \vec{r}_j - \vec{r}_i) \right) \ne \emptyset \right] \\ + & \infty \cdot \sum_{i=0}^{N_\mathrm{particles}-1} \sum_{j=i}^{N_\mathrm{particles}-1} \sum_{\vec{A} \in B_\mathrm{images}, \vec{A} \ne \vec{0}} \left[ \mathrm{overlap}\left( S_i(\mathbf{q}_i), S_j(\mathbf{q}_j, \vec{r}_j - (\vec{r}_i + \vec{A})) \right) \ne \emptyset \right]\end{split}\]

where the square brackets denote the Iverson bracket.

Note

While this notation is written in as sums over all particles HPMCIntegrator uses spatial data structures to evaluate these calculations efficiently. Similarly, while the overlap test is notated as a set intersection, HPMCIntegrator employs efficient computational geometry algorithms to determine whether there is or is not an overlap.

Classes

• ConvexPolygon
• ConvexPolyhedron
• ConvexSpheropolygon
• ConvexSpheropolyhedron
• ConvexSpheropolyhedronUnion
• Ellipsoid
• FacetedEllipsoid
• FacetedEllipsoidUnion
• HPMCIntegrator
• Polyhedron
• SimplePolygon
• Sphere
• SphereUnion
• Sphinx