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., &amp; 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 <code>F</code> function The data is shared as a HDF5 file <code>pw85_ref_data-YYYYMMDD.h5</code>, which contains the following datasets (to be described below) <code>directions</code>: a 12×3 array, <code>F</code>: a 108×108×12×9 array, <code>lambdas</code>: a length-9 array, <code>radii</code>: a length-3 array, <code>spheroids</code>: a 108×6 array. The attached Python script <code>pw85_gen_ref_data.py</code> 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 <code>(0, 1)</code> of the <code>F</code> function which is defined as follows [see Eq. (3.7) in PW85, with slightly different notations] <pre><code>F(λ) = λ(1-λ)r₁₂ᵀ⋅Q⁻¹⋅r₁₂,</code></pre> where <code>0 ≤ λ ≤ 1</code> is a scalar, <code>r₁₂</code> is the center-to-center vector. <code>Q</code> is the matrix defined as follows <pre><code>Q = (1-λ)Q₁ + λQ₂,</code></pre> where <code>Qᵢ</code> is the symmetric, positive definite matrix that defines ellipsoid <code>Ωᵢ</code> through <pre><code>m ∈ Ωᵢ iff (m-cᵢ)ᵀ⋅Qᵢ⁻¹⋅(m-cᵢ) ≤ 1,</code></pre> where <code>cᵢ</code> is the center of <code>Ωᵢ</code>. Then, the contact function <code>F₁₂</code> is defined as the maximum of <code>F</code> [see Eq. (3.8) in PW85] <pre><code>F₁₂(r₁₂, Q₁, Q₂) = max{ F(λ), 0 ≤ λ ≤ 1 }.</code></pre> Parametrization The reference data is restricted to spheroids (equatorial radius: <code>aᵢ</code>; polar radius: <code>cᵢ</code>; direction of axis of revolution: <code>nᵢ</code>) <pre><code>Qᵢ = aᵢ²I + (cᵢ²-aᵢ²)nᵢᵀ⋅nᵢ,</code></pre> (<code>I</code>: identity matrix). The radii take the following values <pre><code>aᵢ, cᵢ ∈ {0.01999, 1.999, 9.999}.</code></pre> These values of the radii are stored in the <code>radii</code> dataset of the HDF5 file. The orientations <code>nᵢ</code> coincide with the vertices of an icosahedron <pre><code>nᵢ = [0, ±u, ±v]ᵀ or nᵢ = [±v, 0, ±u]ᵀ or nᵢ = [±u, ±v, 0]ᵀ,</code></pre> where <pre><code> 1 φ 1+√5 u = ───────, v = ─────── and φ = ────. √(1+φ²) √(1+φ²) 2</code></pre> The orientations are stored in the <code>directions</code> dataset as a 12×3 array. The matrices <code>Qᵢ</code> are precomputed and stored in the <code>spheroids</code> dataset as a 108×6 array (note: 108 = 12 orientations × 3 equatorial radii × 3 polar radii). <code>spheroids[i, :]</code> stores the upper triangular part of the corresponding matrix in row-major order <pre><code>⎡ spheroids[i, 0] spheroids[i, 1] spheroids[i, 2] ⎤ ⎢ spheroids[i, 3] spheroids[i, 4] ⎥. ⎣ sym. spheroids[i, 5] ⎦</code></pre> The scalar <code>λ</code> takes tabulated values (see the <code>lambdas</code> dataset) <pre><code>λ ∈ {0.1, 0.2, …, 0.9}.</code></pre> Note that <code>λ = 0.0</code> and <code>λ = 1.0</code> are excluded, since <code>F</code> is uniformly 0 in that case. Reference values of the <code>F</code> function The reference values of the function <code>F</code> are stored in the <code>F</code> dataset, which is a 108×108×12×9, such that <code>F[i, j, h, k]</code> is the value of <code>F</code> for <pre><code>Q₁ = spheroids[i], Q₂ = spheroids[j], r₁₂ = directions[h] and λ = lambdas[k].</code></pre> Note that the <code>r₁₂</code> vector takes values in the <code>directions</code> dataset. In other words, only unit-length center-to-center vectors are considered here. Indeed, <code>F</code> trivially depends on the norm of <code>r₁₂</code>, 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 <em>not</em> implemented yet, as reference values of the contact function were not deemed useful. Indeed, once <code>F</code> is validated, it is straightforward to check that the implementation of <code>F₁₂</code> to be tested indeed maximizes <code>F</code>. The reference values of the contact function <code>F₁₂</code> are stored in the <code>contact_function</code> dataset, which is a 108×108×12×3 array, such that <code>contact_function[i, j, h, k]</code> is the value of <code>F₁₂</code> for <pre><code>Q₁ = spheroids[i], Q₂ = spheroids[j] and r₁₂ = radii[h] * directions[k].</code></pre> Note that the <code>r₁₂</code> vector is not normed, here.