# Plan: Pairwise Contact Area Potential
## Context
The faunus crate has a many-body `ContactTessellationEnergy` using radical tessellation (voronota-ltr) to compute inter-body contact areas weighted by surface tension γ. This is a **pairwise-additive approximation** for the `interatomic` crate: an analytical contact disc area on the radical (Voronoi) plane between two probe-expanded spheres, weighted by γ.
## Physics
For two spheres with probe-expanded radii R = σ_a/2 + p and S = σ_b/2 + p, separated by distance d, the radical plane divides space where the power w.r.t. both spheres is equal. The contact disc area on this plane is:
```
A(q) = π × [4qR² − (q + R² − S²)²] / (4q) for q ≤ (R + S)²
A(q) = 0 otherwise
```
where q = d² = distance_squared (no sqrt needed).
**Energy:** `U(q) = γ × A(q)`
**Force** (−∂U/∂(r²)):
```
F(q) = γπ(q² − (R² − S²)²) / (4q²) for q < (R + S)²
F(q) = 0 otherwise
```
**Equal spheres** (R = S) simplify to:
- `U(q) = γπ(R² − q/4)`
- `F = γπ/4` (constant)
**Properties:**
- Natural finite cutoff at d = R + S = (σ_a + σ_b)/2 + 2p
- Energy continuous (zero at cutoff)
- Force discontinuous at cutoff (acceptable for MC; can be splined for MD)
- Returns 0 if one sphere contains the other (d ≤ |R − S|)
- Sign of γ determines attraction (γ < 0) or repulsion (γ > 0)
## Implementation
### New file: `interatomic/src/twobody/contact.rs`
**Internal struct:**
```rust
pub struct ContactArea {
gamma_pi: f64, // γ × π (precomputed)
radius_a_sq: f64, // R² = (σ_a/2 + p)²
radius_b_sq: f64, // S² = (σ_b/2 + p)²
cutoff_squared: f64, // (R + S)²
delta_sq: f64, // R² − S² (precomputed for force)
}
```
**Serde helper struct** (like AshbaughHatch pattern):
```rust
#[cfg(feature = "serde")]
#[derive(Deserialize, Serialize)]
#[serde(deny_unknown_fields)]
struct ContactAreaSerde {
gamma: f64, // γ
sigma_a: f64, // σ_a (diameter of particle a)
sigma_b: f64, // σ_b (diameter of particle b)
#[serde(default)]
probe_radius: f64, // p (default: 0.0)
}
```
`From<ContactAreaSerde> for ContactArea` computes R = σ_a/2 + p, S = σ_b/2 + p, then precomputed fields.
`From<ContactArea> for ContactAreaSerde` reverses (sqrt, ÷π, etc.).
**Constructors:**
- `new(gamma, radius_a, radius_b)` — takes probe-expanded radii directly
- `from_sigma(gamma, sigma_a, sigma_b, probe_radius)` — convenience
**Trait impls:**
- `IsotropicTwobodyEnergy` — energy and analytical force (sqrt-free)
- `Cutoff` — finite cutoff at R + S
- `Debug, Clone, PartialEq`
### Modify: `interatomic/src/twobody/mod.rs`
Add `mod contact;` and `pub use contact::ContactArea;`
### Tests (in `contact.rs`)
1. Equal spheres: energy at r=0 equals `γπR²`, at cutoff equals 0, beyond cutoff equals 0
2. Unequal spheres: energy at cutoff equals 0, energy positive/negative matches γ sign
3. Contained sphere (d < |R−S|): energy = 0
4. Analytical force matches numerical central difference
5. Force constant `γπ/4` for equal spheres
6. Serde round-trip (serialize → deserialize → same energy)
## Verification
```bash
cd /Users/mikael/github/faunus-rs/interatomic
cargo test contact
cargo clippy
```