Reference data for the Perram and Wertheim (1985) contact function of ellipsoids

Reference data for the Perram and Wertheim (1985) contact function of ellipsoids

This dataset provides reference values of the contact function of two ellipsoids, as defined by Perram and Wertheim (Perram, J. W., & Wertheim, M. S. (1985). Statistical mechanics of hard ellipsoids. I. Overlap algorithm and the contact function. Journal of Computational Physics, 58(3), 409–416. DOI:10.1016/0021-9991(85)90171-8). This paper will be referred to as PW85 in what follows.

Reference values of the F function

The data is shared as a HDF5 file pw85_ref_data-YYYYMMDD.h5, which contains the following datasets (to be described below)

directions: a 12×3 array,

F: a 108×108×12×9 array,

lambdas: a length-9 array,

radii: a length-3 array,

spheroids: a 108×6 array.

The attached Python script pw85_gen_ref_data.py was used to generate the data; it uses the mpmath library.

Mathematical definition of the contact function

The contact function is defined in PW85 as the maximum over (0, 1) of the F function which is defined as follows [see Eq. (3.7) in PW85, with slightly different notations]

F(λ) = λ(1-λ)r₁₂ᵀ⋅Q⁻¹⋅r₁₂,

where 0 ≤ λ ≤ 1 is a scalar, r₁₂ is the center-to-center vector. Q is the matrix defined as follows

Q = (1-λ)Q₁ + λQ₂,

where Qᵢ is the symmetric, positive definite matrix that defines ellipsoid Ωᵢ through

m ∈ Ωᵢ iff (m-cᵢ)ᵀ⋅Qᵢ⁻¹⋅(m-cᵢ) ≤ 1,

where cᵢ is the center of Ωᵢ. Then, the contact function F₁₂ is defined as the maximum of F [see Eq. (3.8) in PW85]

F₁₂(r₁₂, Q₁, Q₂) = max{ F(λ), 0 ≤ λ ≤ 1 }.

Parametrization

The reference data is restricted to spheroids (equatorial radius: aᵢ; polar radius: cᵢ; direction of axis of revolution: nᵢ)

Qᵢ = aᵢ²I + (cᵢ²-aᵢ²)nᵢᵀ⋅nᵢ,

(I: identity matrix). The radii take the following values

aᵢ, cᵢ ∈ {0.01999, 1.999, 9.999}.

These values of the radii are stored in the radii dataset of the HDF5 file. The orientations nᵢ coincide with the vertices of an icosahedron

nᵢ = [0, ±u, ±v]ᵀ or nᵢ = [±v, 0, ±u]ᵀ or nᵢ = [±u, ±v, 0]ᵀ,

where

1 φ 1+√5
u = ───────, v = ─────── and φ = ────.
√(1+φ²) √(1+φ²) 2

The orientations are stored in the directions dataset as a 12×3 array. The matrices Qᵢ are precomputed and stored in the spheroids dataset as a 108×6 array (note: 108 = 12 orientations × 3 equatorial radii × 3 polar radii). spheroids[i, :] stores the upper triangular part of the corresponding matrix in row-major order

⎡ spheroids[i, 0] spheroids[i, 1] spheroids[i, 2] ⎤
⎢ spheroids[i, 3] spheroids[i, 4] ⎥.
⎣ sym. spheroids[i, 5] ⎦

The scalar λ takes tabulated values (see the lambdas dataset)

λ ∈ {0.1, 0.2, …, 0.9}.

Note that λ = 0.0 and λ = 1.0 are excluded, since F is uniformly 0 in that case.

Reference values of the F function

The reference values of the function F are stored in the F dataset, which is a 108×108×12×9, such that F[i, j, h, k] is the value of F for

Q₁ = spheroids[i], Q₂ = spheroids[j], r₁₂ = directions[h] and λ = lambdas[k].

Note that the r₁₂ vector takes values in the directions dataset. In other words, only unit-length center-to-center vectors are considered here. Indeed, F trivially depends on the norm of r₁₂, which is therefore not considered here in order to reduce the size of the dataset.

Reference values of the contact function

Note: the following is not implemented yet, as reference values of the contact function were not deemed useful. Indeed, once F is validated, it is straightforward to check that the implementation of F₁₂ to be tested indeed maximizes F.

The reference values of the contact function F₁₂ are stored in the contact_function dataset, which is a 108×108×12×3 array, such that contact_function[i, j, h, k] is the value of F₁₂ for

Q₁ = spheroids[i], Q₂ = spheroids[j] and r₁₂ = radii[h] * directions[k].

Note that the r₁₂ vector is not normed, here.