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

<p>Reference data for the Perram and Wertheim (1985) contact function of ellipsoids</p> <p>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&ndash;416. <a href="https://doi.org/10.1016/0021-9991(85)90171-8">DOI:10.1016/0021-9991(85)90171-8</a>). This paper will be referred to as PW85 in what follows.</p> <p>Reference values of the <code>F</code> function</p> <p>The data is shared as a HDF5 file <code>pw85_ref_data-YYYYMMDD.h5</code>, which contains the following datasets (to be described below)</p> <ul> <li><code>directions</code>: a 12&times;3 array,</li> <li><code>F</code>: a 108&times;108&times;12&times;9 array,</li> <li><code>lambdas</code>: a length-9 array,</li> <li><code>radii</code>: a length-3 array,</li> <li><code>spheroids</code>: a 108&times;6 array.</li> </ul> <p>The attached Python script <code>pw85_gen_ref_data.py</code> was used to generate the data; it uses the <a href="http://mpmath.org/">mpmath</a> library.</p> <p>Mathematical definition of the contact function</p> <p>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]</p> <pre><code>F(λ) = λ(1-λ)r₁₂ᵀ⋅Q⁻¹⋅r₁₂,</code></pre> <p>where <code>0 &le; &lambda; &le; 1</code> is a scalar, <code>r₁₂</code> is the center-to-center vector. <code>Q</code> is the matrix defined as follows</p> <pre><code>Q = (1-λ)Q₁ + λQ₂,</code></pre> <p>where <code>Qᵢ</code> is the symmetric, positive definite matrix that defines ellipsoid <code>&Omega;ᵢ</code> through</p> <pre><code>m ∈ Ωᵢ iff (m-cᵢ)ᵀ⋅Qᵢ⁻¹⋅(m-cᵢ) ≤ 1,</code></pre> <p>where <code>cᵢ</code> is the center of <code>&Omega;ᵢ</code>. Then, the contact function <code>F₁₂</code> is defined as the maximum of <code>F</code> [see Eq. (3.8) in PW85]</p> <pre><code>F₁₂(r₁₂, Q₁, Q₂) = max{ F(λ), 0 ≤ λ ≤ 1 }.</code></pre> <p>Parametrization</p> <p>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>)</p> <pre><code>Qᵢ = aᵢ²I + (cᵢ²-aᵢ²)nᵢᵀ⋅nᵢ,</code></pre> <p>(<code>I</code>: identity matrix). The radii take the following values</p> <pre><code>aᵢ, cᵢ ∈ {0.01999, 1.999, 9.999}.</code></pre> <p>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</p> <pre><code>nᵢ = [0, ±u, ±v]ᵀ or nᵢ = [±v, 0, ±u]ᵀ or nᵢ = [±u, ±v, 0]ᵀ,</code></pre> <p>where</p> <pre><code> 1 φ 1+√5 u = ───────, v = ─────── and φ = ────. √(1+φ²) √(1+φ²) 2</code></pre> <p>The orientations are stored in the <code>directions</code> dataset as a 12&times;3 array. The matrices <code>Qᵢ</code> are precomputed and stored in the <code>spheroids</code> dataset as a 108&times;6 array (note: 108 = 12 orientations &times; 3 equatorial radii &times; 3 polar radii). <code>spheroids[i, :]</code> stores the upper triangular part of the corresponding matrix in row-major order</p> <pre><code>⎡ spheroids[i, 0] spheroids[i, 1] spheroids[i, 2] ⎤ ⎢ spheroids[i, 3] spheroids[i, 4] ⎥. ⎣ sym. spheroids[i, 5] ⎦</code></pre> <p>The scalar <code>&lambda;</code> takes tabulated values (see the <code>lambdas</code> dataset)</p> <pre><code>λ ∈ {0.1, 0.2, …, 0.9}.</code></pre> <p>Note that <code>&lambda; = 0.0</code> and <code>&lambda; = 1.0</code> are excluded, since <code>F</code> is uniformly 0 in that case.</p> <p>Reference values of the <code>F</code> function</p> <p>The reference values of the function <code>F</code> are stored in the <code>F</code> dataset, which is a 108&times;108&times;12&times;9, such that <code>F[i, j, h, k]</code> is the value of <code>F</code> for</p> <pre><code>Q₁ = spheroids[i], Q₂ = spheroids[j], r₁₂ = directions[h] and λ = lambdas[k].</code></pre> <p>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.</p> <p>Reference values of the contact function</p> <p>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>.</p> <blockquote> <p>The reference values of the contact function <code>F₁₂</code> are stored in the <code>contact_function</code> dataset, which is a 108&times;108&times;12&times;3 array, such that <code>contact_function[i, j, h, k]</code> is the value of <code>F₁₂</code> for</p> <pre><code>Q₁ = spheroids[i], Q₂ = spheroids[j] and r₁₂ = radii[h] * directions[k].</code></pre> <p>Note that the <code>r₁₂</code> vector is not normed, here.</p> </blockquote>