Fourier Surface

Added in version 2.0.0.

With surface_method='fourier', MembraneCurvature fits a smooth periodic surface to atomistic or bead height data using a truncated 2D Fourier series. The Fourier method comprises six steps:

Choose a non-redundant set of Fourier modes:
- fourier_mode_list()
- n_fourier_parameters()
Compute the corresponding wave numbers:
- _compute_wavevector()
Build a design matrix from cosine and sine basis functions
- _build_fourier_matrix()
Solve the least-squares problem for the Fourier coefficients
Reconstruct the continuous surface on a grid:
- _bin_centre_mesh()
- _eval_fourier_surface() with derivatives=False
Evaluate partial derivatives analytically:
- _eval_fourier_surface() with derivatives=True, or
- fourier_height_derivatives_from_atoms()

Public entry points fourier_height_from_atoms() and fourier_height_derivatives_from_atoms() chain these steps on bin centres.

Note that the use of a non-redundant set of Fourier modes is justified by the conjugate symmetry of Fourier coefficients for real-valued functions:

\[F(-m, -n) = F(m, n)^*\]

We therefore keep only the non-redundant half of the mode grid. The \(m = 0\) row has the additional symmetry

\[F(0, -n) = F(0, n)^*\]

so we keep only \(n \ge 0\) there and store the mean term \(A_{00}\) separately.

A diagram of a non-redundant Fourier mode grid looks like:

n
^
N  |  ✓  |  ✓  ✓  ✓ ... ✓  |
:  |  ✓  |  ✓  ✓  ✓ ... ✓  |
2  |  ✓  |  ✓  ✓  ✓ ... ✓  |
1  |  ✓  |  ✓  ✓  ✓ ... ✓  |
0  | A00 |  ✓  ✓  ✓ ... ✓  |   <- right block includes n=0; left cell is A00
-1 |  x  |  ✓  ✓  ✓ ... ✓  |
-2 |  x  |  ✓  ✓  ✓ ... ✓  |
:  |  x  |  ✓  ✓  ✓ ... ✓  |
-N |  x  |  ✓  ✓  ✓ ... ✓  |
   +-----+-------------------+--> m
      0     1   2   3  ...  M
      ^           ^
      |           |
 modes_with_     modes_with_
 m_equal_zero    m_positive
 (n=1..N only)   (all n, -N..N)

where the ✓ symbol denotes that the wave includes 2 parameters each: \(cos\), and \(sin\) coefficients, and where the x symbols denotes that the wave is excluded to avoid redundancy by \(F(0,-n) = F(0,n)^*\), only occurs at \(m=0\).

The associated wavevector for each mode is

\[\mathbf{k}_{mn} = \left(\frac{2\pi m}{L_x}, \frac{2\pi n}{L_y}\right)\]

Important

The simulation box defines the fundamental periodic domain and therefore sets the largest allowed wavelength in the system. The surface is represented as a sum of periodic waves that are compatible with these boundary conditions, meaning all basis functions must fit exactly within the box without discontinuities at the boundaries.

As a result, only discrete Fourier modes are allowed, corresponding to integer wave numbers:

\[k_x = \frac{2\pi m}{L_x}, \qquad k_y = \frac{2\pi n}{L_y}\]

with \(m, n \in \mathbb{Z}\).

The box size \(L_x\) and \(L_y\) therefore determines the spacing of allowed wavevectors, while higher modes represent progressively shorter wavelengths.

For more details on the Fourier method, refer to the 2.1 Fourier method section in the Algorithm page.

References

The Fourier method in MembraneCurvature uses the Fourier surface modeling [CAG2009] and molecular Fourier shape descriptors [JMG1988] as conceptual references:

[CAG2009]

Shen et al., Fourier method for large-scale surface modeling and registration, Computers & Graphics (2009), doi: 10.1016/j.cag.2009.03.002.

[JMG1988]

Leicester et al., Description of molecular surface shape using Fourier descriptors, Journal of Molecular Graphics (1988), doi: 10.1016/0263-7855(88)85008-2.

Functions

membrane_curvature.fourier_surface._bin_centre_mesh(Lx: float, Ly: float, n_x_bins: int, n_y_bins: int) → tuple[ndarray, ndarray][source]

Relative-coordinate mesh at bin centres for surface reconstruction (step 5).

Bin centres lie at (i + 1/2) L / n_bins for each axis, compatible with the periodic domain [0, L).

Parameters:

Lx (float) – Period length in \(x\).
Ly (float) – Period length in \(y\).
n_x_bins (int) – Number of bins along \(x\).
n_y_bins (int) – Number of bins along \(y\).

Returns:

X_rel (ndarray) – Mesh of relative \(x\) coordinates, shape (n_x_bins, n_y_bins).
Y_rel (ndarray) – Mesh of relative \(y\) coordinates, same shape as X_rel.

Raises:

ValueError – If n_x_bins or n_y_bins is not a positive integer. The public entry points (fourier_height_from_atoms(), fourier_height_derivatives_from_atoms(), MembraneCurvature) already validate via validate_positive_bin_counts().

membrane_curvature.fourier_surface._build_fourier_matrix(x_positions: ndarray, y_positions: ndarray, box_length_x: float, box_length_y: float, fourier_modes: list[tuple[int, int]]) → ndarray[source]

Build the least-squares design matrix: basis functions evaluated at sample points.

Implements step 3 (and uses _compute_wavevector() for step 2 inside the loop). Each row is an atom position; columns are 1, then cos(k\cdot r), sin(k\cdot r) pairs for each mode, so that fitting recovers \(z(x,y)\) at the atoms.

In the matrix, each row corresponds to a spatial point \((x_i, y_i)\), and each column corresponds to a basis function:

\[[1, cos(k·r), sin(k·r), cos(k·r), sin(k·r), ...]\]

where the wave vectors are given by \(k=(2\pi m/L_x, 2\pi n/L_y)\).

Parameters:

x_positions (ndarray) – \(x\) coordinates of points within one periodic box (same shape as y_positions).
y_positions (ndarray) – \(y\) coordinates of points within one periodic box.
box_length_x (float) – Period length \(L_x\).
box_length_y (float) – Period length \(L_y\).
fourier_modes (list of tuple of int) – Non-redundant (m, n) indices, e.g. from fourier_mode_list().

Returns:

design_matrix – Basis evaluated at each point; dtype is numpy.float64.

Return type:

ndarray, shape (n_points, 1 + 2 * n_modes)

membrane_curvature.fourier_surface._compute_wavevector(mode_x, mode_y, kx_base, ky_base)[source]

Wavevector components for a single mode (step 2 of the module workflow).

Builds \((k_x, k_y) = (m\, k_{x,\mathrm{base}}, n\, k_{y,\mathrm{base}})\) used in the phase \(k_x x + k_y y\) for that mode.

Parameters:

mode_x (int) – Mode index \(m\).
mode_y (int) – Mode index \(n\).
kx_base (float) – Reciprocal step \(2\pi / L_x\).
ky_base (float) – Reciprocal step \(2\pi / L_y\).

Returns:

kx (float) – mode_x * kx_base.
ky (float) – mode_y * ky_base.

membrane_curvature.fourier_surface._eval_fourier_surface(X_rel: ndarray, Y_rel: ndarray, Lx: float, Ly: float, A00: float, coeffs: dict[tuple[int, int], tuple[float, float]], *, derivatives: Literal[False]) → tuple[ndarray][source]

membrane_curvature.fourier_surface._eval_fourier_surface(X_rel: ndarray, Y_rel: ndarray, Lx: float, Ly: float, A00: float, coeffs: dict[tuple[int, int], tuple[float, float]], *, derivatives: Literal[True]) → tuple[ndarray, ndarray, ndarray, ndarray, ndarray, ndarray]

Evaluate the fitted Fourier sum on a relative-coordinate grid (steps 5-6).

Step 5: reconstruct \(Z(x,y) = A_{00} + \sum_{(m,n)} A_{mn}\cos\phi_{mn} + B_{mn}\sin\phi_{mn}\) with \(\phi_{mn} = k_x x + k_y y\). Step 6 (optional): accumulate first and second partial derivatives for Monge-gauge curvature using the analytic derivatives of the same sum.

Parameters:

X_rel (ndarray) – Relative \(x\) coordinates on the periodic domain (same shape as Y_rel).
Y_rel (ndarray) – Relative \(y\) coordinates on the periodic domain.
Lx (float) – Period length in \(x\) (same units as X_rel).
Ly (float) – Period length in \(y\).
A00 (float) – Mean term \(A_{00}\).
coeffs (dict[tuple[int, int], tuple[float, float]]) – Mode amplitudes (A, B) from _fourier_fit_from_atoms() / _unpack_coefficients().
derivatives (bool) – If False, return (Z,). If True, return (Z, fx, fy, fxx, fyy, fxy).

Returns:

If derivatives is False, the 1-tuple (Z,) where Z matches the shape of X_rel. If derivatives is True, the 6-tuple with each grid matching X_rel.

Return type:

tuple of ndarray

Notes

The derivatives=False branch returns (Z,) — a 1-tuple, not a bare ndarray — to keep a single return type (tuple[ndarray, ...]) regardless of the flag. Callers therefore unwrap with [0].

membrane_curvature.fourier_surface._fourier_fit_from_atoms(positions: ndarray, x_range: tuple[float, float], y_range: tuple[float, float], M: int, N: int, rcond: float | None = None) → tuple[float, float, float, dict[tuple[int, int], tuple[float, float]]][source]

Fit Fourier coefficients to atom heights (steps 1-4; no bin grid).

Builds the design matrix, solves for \(A_{00}\) and the mode amplitudes, and returns domain sizes and coefficients for subsequent grid evaluation (steps 5-6).

Parameters:

positions (ndarray, shape (n_atoms, 3)) – Atom coordinates; the third column (positions[:, 2]) is the height \(z\).
x_range, y_range (tuple of float) – (min, max) extents of the periodic domain.
M, N (int) – Maximum Fourier mode indices (see fourier_mode_list()).
rcond (float or None, optional) – Relative cutoff for small singular values in _solve_design_least_squares_svd().

Returns:

Lx (float) – Period length \(L_x = x_{\max} - x_{\min}\).
Ly (float) – Period length \(L_y = y_{\max} - y_{\min}\).
A00 (float) – Mean height coefficient \(A_{00}\).
coeffs (dict[tuple[int, int], tuple[float, float]]) – Mapping (m, n) -> (A, B) cosine and sine coefficients per mode.

Raises:

ValueError – If positions cannot be converted to a float64 NumPy array, or has the wrong shape (see validate_positions()); if x_range / y_range have non-positive width (see validate_positive_domain_widths()); or if M / N are negative (see fourier_mode_list()).

Warns:

UserWarning – If the design matrix is rank-deficient or underdetermined (same condition as fourier_height_from_atoms() and fourier_height_derivatives_from_atoms()).

membrane_curvature.fourier_surface._harmonic_height_and_phi_derivatives(A: float, B: float, cos_phase: ndarray, sin_phase: ndarray) → tuple[ndarray, ndarray, ndarray][source]

Height and \(\phi\)-derivatives for one mode \(h(\phi)=A\cos\phi+B\sin\phi\).

Used internally by _eval_fourier_surface() to accumulate fx, fy, fxx, fyy, fxy using the chain rule. For example, \(\partial_x h = (\partial_\phi h)\, \partial_x \phi\).

Parameters:

A, B (float) – Cosine and sine coefficients for this mode.
cos_phase (ndarray) – cos(\phi), same shape as sin_phase.
sin_phase (ndarray) – sin(\phi).

Notes

This helper returns derivatives with respect to phase \(\phi\) (h_phi, h_phiphi), not spatial derivatives. Spatial derivatives (fx, fy, fxx, fyy, fxy) are accumulated in _eval_fourier_surface() via the chain rule using \(k_x\) and \(k_y\).