Phase separation & droplet physics

The pimms.lemonade.phase_separation and pimms.lemonade.surface_tension modules quantify liquid-liquid phase separation of a PIMMS system: the coexistence (binodal) densities, the amount of condensed material, cluster-size order parameters, droplet shape and the interfacial tension.

import pimms.lemonade as lemonade
from pimms.lemonade import phase_separation as ps
from pimms.lemonade import surface_tension as st

traj = lemonade.load(xtc="traj.xtc", pdb="START.pdb", keyfile="KEYFILE.kf")

Densities throughout are occupied lattice-site fractions in [0, 1], so they are directly comparable across box sizes.

Clusters and condensates

Clusters are connected groups of chains, found per frame by Frame.clusters (largest first). The largest cluster - the condensate - is frame.droplet. Each cluster carries its own geometry (The object hierarchy): radius_of_gyration, asphericity, sphericity, volume, surface_area, density, radial_density_profile() and its composition.

Order parameters

Simple, per-frame measures of how phase separated the system is:

ps.condensed_fraction(traj)          # (n_frames,) fraction of beads in the largest cluster
ps.largest_cluster_size(traj)        # (n_frames,) beads in the largest cluster
ps.number_of_clusters(traj, min_beads=2)
ps.cluster_size_distribution(traj)   # all cluster sizes pooled across frames

condensed_fraction is the basic order parameter: near zero when the system is well mixed, approaching one when most material collects into a single condensate.

The binodal (coexistence densities)

The dense- and dilute-phase densities are read off a density profile fit to a hyperbolic tangent. Two geometries are supported.

Droplet (spherical). A radial profile about the condensate centre of mass - occupied fraction as a function of distance - falls from the dense core to the dilute background:

r, rho = ps.radial_density_profile(traj)
fit = ps.fit_radial_profile(r, rho)
fit.rho_dense, fit.rho_dilute        # coexistence densities
fit.radius, fit.interface_width      # droplet radius and interface width

The fit is \(\rho(r) = \tfrac12(\rho_d + \rho_v) - \tfrac12(\rho_d - \rho_v)\tanh((r-R)/w)\).

Slab. For a slab condensate that spans the periodic plane and is bounded along one axis (the geometry of the slab_phase_separation demo), a 1D profile along the long axis - with the slab re-centred each frame so it does not smear - is fit to a two-interface tanh:

z, rho = ps.slab_density_profile(traj, axis=2)   # axis defaults to the longest
fit = ps.fit_slab_profile(z, rho)
fit.rho_dense, fit.rho_dilute, fit.interface_width

One call for everything

analyze() runs the whole pipeline and auto-detects the geometry (slab if one box axis is much longer than the others, else spherical):

result = ps.analyze(traj)

result.geometry                      # 'sphere' or 'slab'
result.rho_dense, result.rho_dilute  # binodal
result.condensed_fraction            # time-averaged
result.binodal.interface_width
result.shape                         # {'radius_of_gyration', 'sphericity', ...}
result.is_phase_separated            # heuristic: density gap AND most material condensed
result.profile                       # (coordinate, density) for plotting

Surface tension from undulations

Both surface-tension estimators use capillary-wave theory - the interface’s fluctuation spectrum. Because PIMMS uses \(\exp(-\Delta E/T)\) with \(k_B = 1\), the trajectory temperature is \(k_B T\), and the returned \(\gamma\) is in reduced units (interaction energy per lattice area).

st.surface_tension(traj)             # auto-dispatch by geometry
st.slab_surface_tension(traj)        # planar capillary waves  (robust)
st.droplet_surface_tension(traj)     # spherical-harmonic shape fluctuations

  • Slab - the two flat interfaces of a box-spanning condensate have a height field \(h(x,y)\) obeying \(\langle|h(q)|^2\rangle = k_BT/(\gamma A q^2)\); fitting the low-\(q\) spectrum gives \(\gamma\). This is the reliable method.

  • Droplet - the radius \(R(\theta,\phi)\) fluctuates in spherical-harmonic modes with \(\langle|u_{lm}|^2\rangle = k_BT/(\gamma R_0^2 (l-1)(l+2))\) for \(l \ge 2\). Best-effort: it needs a single, compact, reasonably large droplet sampled over many frames.

Each returns a SurfaceTension with the estimate gamma, a per-mode spread gamma_std (an uncertainty proxy), the number of modes used, and the raw spectrum for inspection:

result = st.slab_surface_tension(traj)
result.gamma, result.gamma_std
q, power = result.spectrum           # inspect the capillary spectrum

Warning

Surface tension is intrinsically noisy for small lattice condensates (few long-wavelength modes, rough interfaces). Always check gamma_std and plot the spectrum; use a large interface (a big slab cross-section, or a single large droplet) and many frames for a precise value. If the temperature is unknown (loaded without a keyfile), pass temperature= explicitly.

Worked example

import pimms.lemonade as lemonade
from pimms.lemonade import phase_separation as ps
from pimms.lemonade import surface_tension as st

traj = lemonade.load(xtc="traj.xtc", pdb="START.pdb", keyfile="KEYFILE.kf")

result = ps.analyze(traj)
if result.is_phase_separated:
    print(f"{result.geometry}: rho_dense = {result.rho_dense:.3f}, "
          f"rho_dilute = {result.rho_dilute:.3f}, "
          f"condensed fraction = {result.condensed_fraction:.2f}")
    gamma = st.surface_tension(traj)
    print(f"surface tension = {gamma.gamma:.2f} +/- {gamma.gamma_std:.2f} (kT units)")