API Reference

Solver

High-level solver API for moire relaxation.

Typical usage:

from moire_metrology import RelaxationSolver
from moire_metrology.interfaces import GRAPHENE_GRAPHENE

solver = RelaxationSolver()
result = solver.solve(
    moire_interface=GRAPHENE_GRAPHENE,
    theta_twist=2.0,
)
result.plot_stacking()
class moire_metrology.solver.SolverConfig(method='trust-ncg', max_iter=200, gtol=1e-06, rtol=0.0001, etol=1e-06, etol_window=5, pixel_size=0.2, n_scale=1, display=True, min_mesh_points=100, beta=1.0, dt0=None, linear_solver='direct', linear_solver_tol=1e-06, linear_solver_maxiter=200, elastic_strain='cauchy', max_iter_discover=300, gtol_discover_factor=10.0)[source]

Bases: object

Configuration for the relaxation solver.

Parameters:

  • method (str) –

    Optimization method.

    • 'newton' – Levenberg-Marquardt damped Newton with a modified per-vertex GSFE Hessian (negative eigenvalues are flipped so the assembled Hessian is positive-definite by construction). Fast and well-behaved for nearly-convex problems (typical at moderate-to-high twist angles), but by design cannot follow negative-curvature directions and therefore cannot cross saddles between basins. In multi-basin landscapes (e.g. low-twist TDBG with significant AB/BA asymmetry, where single-DW and 2DW topologies are both stationary states) the eigenvalue flipping causes the iterate to commit early to whichever basin is locally convex-nearest, with no recourse to escape that basin later. Not recommended as a default; kept available for high-twist / near-convex problems where its fast quadratic convergence is desirable and saddle-crossing is moot.

    • 'trust-ncg' (default) – True trust-region Newton-CG via scipy.optimize.minimize (Steihaug-Toint inner solver) using the true (unmodified) sparse Hessian via matrix-free Hessian- vector products. Detects negative-curvature directions and steps to the trust-region boundary along them — the property that lets it cross saddles, which the 'newton' LM-modified path cannot. Good general-purpose default; supports the same mean-constraint API as the other methods via null-space projection. Caveat: although saddle-crossing is available, from a generic IC like U=0 the inner CG’s quadratic-step preference still biases toward whichever basin is locally quadratic-nearest. For deliberate basin discovery from U=0 in low-twist problems, prefer 'L-BFGS-B' or 'two_phase'.

    • 'pseudo_dynamics' – implicit theta-method on the gradient flow dU/dt = -nabla E. At each step solves (M + beta*dt*H)*dU = -dt*grad with adaptive dt and energy-monitored step rejection. More robust than Newton on stiff/multi-layer systems. Ports the algorithm used in the paper MATLAB code.

    • 'L-BFGS-B' – Gradient-only quasi-Newton via scipy.optimize.minimize. Slower convergence but no Hessian needed. Supports the same mean-constraint API as 'newton' / 'trust-ncg' via null-space projection. Empirically lands in qualitatively different basins than 'trust-ncg' from the same U=0 initial guess in multi-basin landscapes (e.g. low-twist TDBG with AB/BA asymmetry); use this for discovery when you suspect the U=0 IC sits near a saddle.

    • 'two_phase' – L-BFGS-B discovery (using max_iter_discover and gradient tolerance gtol_discover_factor * gtol) followed by trust-ncg polish (using the standard max_iter / gtol / rtol / etol criteria) from the discovery output. Combines L-BFGS-B’s basin-selection behaviour with trust-ncg’s tight 2nd-order convergence. Recommended for low-twist relaxation when both basin choice and high precision matter.

  • max_iter (int) – Maximum number of optimizer iterations.

  • gtol (float) – Absolute gradient norm tolerance for convergence. Convergence is declared when |grad| < gtol. For L-BFGS-B this is passed to scipy as the gtol option (max-norm threshold).

  • rtol (float) – Relative gradient tolerance. Convergence is declared when |grad| / |grad_initial| < rtol, i.e. the gradient has dropped by this factor from its initial value. This criterion matters at low twist angles where the absolute gradient norm at the energy minimum can be O(1-10) due to large total energies, making any reasonable absolute gtol unreachable. For L-BFGS-B the check is applied post-hoc after scipy returns.

  • etol (float) – Energy stagnation tolerance. If the fractional energy improvement over the last etol_window accepted steps falls below etol, convergence is declared (newton and pseudo_dynamics only).

  • etol_window (int) – Number of recent iterations over which to measure energy stagnation.

  • pixel_size (float) – Target mesh element size in nm.

  • n_scale (int) – Number of moire unit cells in each direction.

  • display (bool) – If True, print progress during optimization.

  • min_mesh_points (int) – Floor on the mesh grid resolution per direction. Default 100.

  • beta (float) – Theta-method parameter for the ‘pseudo_dynamics’ solver. 0 = explicit Euler (unstable), 1 = fully implicit backward Euler (default, unconditionally stable). Values in [0.5, 1.0] are stable.

  • dt0 (float or None) – Initial pseudo-time step for the ‘pseudo_dynamics’ solver. If None, chosen automatically as ~1/(K1+G1) so that dt0·|H| ~ 1.

  • linear_solver (str) –

    Linear-system solver used inside each pseudo_dynamics step.

    • 'direct' (default) – build the sparse Hessian explicitly and call spsolve. Per-iteration cost is the LU factorization, which scales like O(n_sol^1.5) on 2D meshes. Best for small/moderate problems.

    • 'iterative' – matrix-free preconditioned MINRES, using only Hessian-vector products and a Jacobi preconditioner from the diagonal of the constant elastic Hessian. Much faster per iteration on large problems (n_sol >> 10^4). Uses MINRES (not CG) so it stays robust when the operator is symmetric but indefinite.

  • linear_solver_tol (float) – Relative-residual tolerance for the iterative linear solver. Only used when linear_solver=’iterative’. Default 1e-6.

  • linear_solver_maxiter (int) – Max iterations for the iterative linear solver per pseudo_dynamics step. Default 200. The outer pseudo_dynamics loop will rebuild the operator and try again with a smaller dt if MINRES fails to converge.

  • elastic_strain (str)

  • max_iter_discover (int)

  • gtol_discover_factor (float)

class moire_metrology.solver.RelaxationSolver(config=None)[source]

Bases: object

Solver for atomic relaxation in twisted 2D heterostructures.

Parameters:

config (SolverConfig | None)

solve(moire_interface=None, *, top_interface=None, bottom_interface=None, n_top=1, n_bottom=1, theta_twist=0.0, delta=None, theta0=0.0, initial_solution=None, constraints=None, fix_top=False, fix_bottom=False, pin_mean=False, mesh=None, mean_constraints=None, **legacy_kwargs)[source]

Solve the relaxation problem.

Parameters:

  • moire_interface (Interface) – The twisted A-B interface between the bottom layer of the top flake and the top layer of the bottom flake. Carries both materials (moire_interface.bottom is the bottom-flake material, moire_interface.top is the top-flake material) as well as the GSFE coefficients for their stacking.

  • top_interface (Interface, optional) – Required iff n_top > 1. The homobilayer interface used between successive layers within the top flake. Must have bottom == top == moire_interface.top.

  • bottom_interface (Interface, optional) – Required iff n_bottom > 1. The homobilayer interface used between successive layers within the bottom flake. Must have bottom == top == moire_interface.bottom.

  • n_top (int) – Number of layers in each flake.

  • n_bottom (int) – Number of layers in each flake.

  • theta_twist (float) – Twist angle in degrees.

  • delta (float or None) – Lattice mismatch. If None, computed from the two materials.

  • theta0 (float) – Lattice orientation angle in degrees.

  • initial_solution (ndarray or None) – Initial guess. If None, starts from zero.

  • constraints (PinnedConstraints or None) – If set, pins certain DOFs to fixed displacements while optimizing the rest. Build via PinningMap.build_constraints(). Mutually exclusive with fix_top / fix_bottom.

  • fix_top (bool) – Pin all DOFs of the topmost layer (top of the top flake) to zero. Use this to clamp the upper free surface of the heterostructure and approximate a semi-infinite top.

  • fix_bottom (bool) – Pin all DOFs of the bottommost layer (bottom of the bottom flake) to zero. Use this to clamp the substrate and approximate a semi-infinite bottom — typical for simulating a twisted flake on a thick substrate (e.g. graphene on graphite).

  • mesh (MoireMesh or None) – Pre-built mesh to use. If None (default), the solver builds a periodic moire-cell mesh from the SolverConfig parameters (pixel_size, n_scale, min_mesh_points). Pass a mesh built via moire_metrology.mesh.generate_finite_mesh (or any other MoireMesh constructor) to relax on a finite, non-periodic domain — typically combined with constraints from PinningMap.build_constraints to pin selected stacking sites in the experimental image.

  • pin_mean (bool)

  • mean_constraints (list | None)

Return type:

RelaxationResult

Notes

Internally the solver still uses the legacy “stack 1” / “stack 2” numbering, where stack 1 is the top flake and stack 2 is the bottom flake. This is an implementation detail of the construction layer; the public API is in terms of moire_interface.top / moire_interface.bottom and the n_top / n_bottom flake sizes. The translation happens once at the top of this method.

Result

Result container for moire relaxation calculations.

class moire_metrology.result.RelaxationResult(mesh, geometry, moire_interface, top_interface, bottom_interface, displacement_x1, displacement_y1, displacement_x2, displacement_y2, total_energy, unrelaxed_energy, gsfe_map, elastic_map1, elastic_map2, solution_vector, optimizer_result)[source]

Bases: object

Container for relaxation results with post-processing and plotting.

Parameters:
mesh

The computational mesh.

Type:

MoireMesh

geometry

Moire geometry used in the calculation.

Type:

MoireGeometry

moire_interface

The twisted A-B interface that drove the calculation. Carries both materials (moire_interface.top, moire_interface.bottom) and the GSFE coefficients of the moire stacking.

Type:

Interface

top_interface

Homobilayer interface used inside the top flake when nlayer1 > 1. None for monolayer top flakes.

Type:

Interface or None

bottom_interface

Homobilayer interface used inside the bottom flake when nlayer2 > 1. None for monolayer bottom flakes.

Type:

Interface or None

displacement_x1, displacement_y1

Displacement fields for the top flake (stack 1).

Type:

ndarray, shape (nlayer1, Nv)

displacement_x2, displacement_y2

Displacement fields for the bottom flake (stack 2).

Type:

ndarray, shape (nlayer2, Nv)

total_energy

Total relaxed energy in meV/uc.

Type:

float

unrelaxed_energy

Total unrelaxed energy in meV/uc.

Type:

float

gsfe_map

GSFE energy density at interface vertices (meV/nm^2).

Type:

ndarray, shape (Nv,)

elastic_map1

Elastic energy density per layer in the top flake (meV/nm^2).

Type:

ndarray, shape (nlayer1, Nv)

elastic_map2

Elastic energy density per layer in the bottom flake (meV/nm^2).

Type:

ndarray, shape (nlayer2, Nv)

solution_vector

Raw optimizer solution vector.

Type:

ndarray

optimizer_result

scipy.optimize result object.

Type:

object

property material1: Material

Top flake material (stack 1).

Convenience accessor — equivalent to self.moire_interface.top. Stack-1 / stack-2 numbering follows the internal solver convention (stack 1 = top flake, stack 2 = bottom flake).

property material2: Material

Bottom flake material (stack 2).

Convenience accessor — equivalent to self.moire_interface.bottom.

property energy_reduction: float

Fractional energy reduction from relaxation.

property converged: bool

Whether the optimizer reported successful convergence.

property convergence_message: str

Human-readable description of the convergence outcome.

local_twist(stack=1, layer=0)[source]

Compute local twist angle (degrees) at each vertex.

local_twist = theta_twist + 0.5*(duy/dx - dux/dy) * 180/pi

Return type:

ndarray

Parameters:
plot_stacking(ax=None, n_tile=2, **kwargs)[source]

Plot the GSFE (stacking energy) map at the interface.

Parameters:

n_tile (int)

plot_elastic_energy(stack=1, layer=0, ax=None, n_tile=2, **kwargs)[source]

Plot elastic energy density for a given layer.

Parameters:
plot_local_twist(stack=1, layer=0, ax=None, n_tile=2, **kwargs)[source]

Plot local twist angle map.

Parameters:
save(path)[source]

Save result to a .npz file.

Parameters:

path (str)

Materials

Per-layer material properties for 2D van der Waals heterostructures.

A Material carries the intra-layer properties of a 2D crystal: lattice constant and 2D elastic moduli. The inter-layer stacking interaction (GSFE) is a property of a pair of materials in contact, not of either material individually, and lives on moire_metrology.interfaces.Interface instead.

Units convention

Lattice constants are in nm. Bulk and shear moduli are stored in meV per unit cell, matching the Carr/Halbertal convention used throughout the package: the elastic energy density per unit cell is

E_elastic_per_uc = (1/2) K (∂_x u_x + ∂_y u_y)²

  • (1/2) G [(∂_x u_x − ∂_y u_y)² + (∂_x u_y + ∂_y u_x)²]

with K = λ + μ (the 2D bulk modulus) and G = μ (the 2D shear modulus). The conversion to literature-standard 2D moduli in N/m is

K [N/m] = K [meV/uc] / (S_uc [m²] × 6.241509e21 meV/J)

where S_uc = (√3/2) × a² is the unit cell area. For graphene this factor is ~329.5 meV/uc per N/m, so K = 69518 meV/uc corresponds to K = 211 N/m — matching the experimental indentation value of Lee et al. Science 2008.

The Material.from_2d_moduli_n_per_m() constructor and the Material.moduli_n_per_m property provide the round-trip in both directions, so user-defined materials can be specified directly in N/m and existing materials can be sanity-checked against the literature.

Sources for the bundled values

  • GRAPHENE: K, G from Carr et al. PRB 98, 224102 (2018), Table I.

  • HBN_AA, HBN_AAP: K, G derived from Falin et al. Nat. Commun. 8, 15815 (2017), monolayer hBN indentation, E_2D ≈ 286 N/m, combined with the literature 2D Poisson ratio ν ≈ 0.21 to give K_2D ≈ 181 N/m, G_2D ≈ 118 N/m. The AA / AA’ distinction is a stacking convention for the hBN/hBN bilayer interface, not a monolayer material property — the Material entries are physically the same; only the paired Interface entries differ.

  • MOSE2, WSE2: K, G from Halbertal et al. Nat. Commun. 12, 242 (2021), SI Table 1, and independently confirmed in the Shabani, Halbertal et al. Nat. Phys. 17, 720 (2021) Methods section.

  • GRAPHENE_BILAYER: a Bernal-stacked graphene bilayer treated as a single effective 2D elastic continuum, with K and G taken as exactly twice the single-layer graphene values per the convention of Halbertal et al. Nat. Commun. 12, 242 (2021), SI Table 1 (the “TDBG” rows). Lattice constant is inherited from graphene.

class moire_metrology.materials.Material(name, lattice_constant, bulk_modulus, shear_modulus)[source]

Bases: object

Per-layer properties of a single 2D material.

Parameters:

  • name (str) – Human-readable name.

  • lattice_constant (float) – In-plane lattice constant in nm.

  • bulk_modulus (float) – 2D bulk modulus K = lambda + mu in meV/unit cell.

  • shear_modulus (float) – 2D shear modulus G = mu in meV/unit cell.

property unit_cell_area: float

Unit cell area in nm^2.

property moduli_n_per_m: tuple[float, float]

Return (K, G) converted to literature-standard units of N/m.

The package stores moduli in meV/unit-cell (the Carr/Halbertal convention used throughout the elastic energy expression). This property converts back to the more familiar 2D N/m units used in the experimental literature, e.g. for sanity-checking a bundled material against published indentation values.

For graphene the conversion factor is ~329.5 meV/uc per N/m, so K = 69518 meV/uc returns ~211 N/m, matching the Lee et al. Science 2008 indentation value.

classmethod from_2d_moduli_n_per_m(name, lattice_constant, bulk_modulus_n_per_m, shear_modulus_n_per_m)[source]

Construct a Material from K, G in literature-standard N/m units.

The package internally stores moduli in meV per unit cell to match the Carr/Halbertal elastic energy convention. Most published 2D elastic constants (Lee et al. Science 2008, Carr et al. PRB 98, etc.) are reported in N/m instead, so this constructor exists to make user-defined materials specifiable in their natural units without forcing the user to redo the conversion arithmetic by hand.

Parameters:

  • name (str) – Human-readable name.

  • lattice_constant (float) – In-plane lattice constant in nm.

  • bulk_modulus_n_per_m (float) – 2D bulk modulus K = λ + μ in N/m.

  • shear_modulus_n_per_m (float) – 2D shear modulus G = μ in N/m.

Return type:

Material

Examples

Reproducing the bundled GRAPHENE from literature values:

>>> g = Material.from_2d_moduli_n_per_m(
...     name="Graphene",
...     lattice_constant=0.247,
...     bulk_modulus_n_per_m=211.0,   # Lee et al. Science 2008
...     shear_modulus_n_per_m=144.0,
... )
>>> int(round(g.bulk_modulus))   # within rounding of paper Table 1
69533
classmethod n_layer_stack(base, n=1, name=None)[source]

Return a Material for n coherently-strained copies of base.

The returned Material has the same lattice constant as base and elastic moduli scaled by n (each additional layer adds linearly to the effective stiffness under the assumption of identical in-plane strain). It is otherwise a first-class Material that can be used anywhere a single-layer one can.

This is the modelling convention used in the literature for twisted double-bilayer graphene (n=2) and twisted double- trilayer graphene (n=3): see e.g. Halbertal et al. Nat. Commun. 12, 242 (2021), SI Table 1, which states K, G were simply taken as twice that of TBG for the TDBG row, i.e. n=2.

Parameters:

  • base (Material) – The per-layer material to stack.

  • n (int, default 1) – Number of coherent layers in the stack. n=1 reproduces base (with a fresh name if one is provided). Must be at least 1.

  • name (str, optional) – Human-readable name for the returned Material. Defaults to f"{base.name} (n={n})".

Return type:

Material

Examples

>>> from moire_metrology import GRAPHENE
>>> bilayer = Material.n_layer_stack(GRAPHENE, n=2)
>>> round(bilayer.bulk_modulus / GRAPHENE.bulk_modulus)
2
>>> bilayer.lattice_constant == GRAPHENE.lattice_constant
True
classmethod from_dict(data)[source]

Build a Material from a plain dict (e.g. parsed from TOML).

Required keys: name, lattice_constant, bulk_modulus, shear_modulus. Extra keys raise a ValueError so that typos in user TOML files surface immediately instead of being silently ignored.

Return type:

Material

Parameters:

data (Mapping[str, Any])

classmethod from_toml(path)[source]

Load a Material from a TOML file.

The file must contain a top-level [material] table with the four required fields:

[material]
name = "MyTMD"
lattice_constant = 0.330      # nm
bulk_modulus = 42000.0        # meV/uc
shear_modulus = 28000.0       # meV/uc

See Material.from_dict() for the field semantics.

Return type:

Material

Parameters:

path (str | Path)

Interfaces

Interlayer interfaces for moire heterostructures.

An Interface describes the registry-dependent stacking energy between two adjacent layers — the GSFE — together with the convention for how the upper layer is offset relative to the lower layer when multi-layer flakes are stacked. GSFE is physically a property of the pair of materials in contact, not of either material individually, so this module separates that concern from the per-material elastic moduli that live in moire_metrology.materials.

The GSFE Fourier expansion follows the Carr et al. convention:

V(v, w) = c0 + c1*(cos(v) + cos(w) + cos(v + w))

  • c2*(cos(v + 2w) + cos(v - w) + cos(2v + w))

  • c3*(cos(2v) + cos(2w) + cos(2v + 2w))

  • c4*(sin(v) + sin(w) - sin(v + w))

  • c5*(sin(2v + 2w) - sin(2v) - sin(2w))

For centrosymmetric homobilayers (e.g. graphene/graphene, hBN-AA) the sin coefficients vanish (c4 = c5 = 0). Heterointerfaces with broken inversion symmetry (e.g. MoSe2/WSe2 H-stacking, hBN AA’) carry non-zero c4, c5.

References

  • Zhou et al., PRB 92, 155438 (2015) — original GSFE parameterizations for graphene and hBN homo/heterointerfaces.

  • Carr et al., PRB 98, 224102 (2018) — Fourier basis convention used here, including the Zhou->Carr coefficient transformation.

  • Shabani, Halbertal et al., Nature Physics 17, 720 (2021) — GSFE parameters for the H-stacked MoSe2/WSe2 heterointerface (already in Carr basis, no transformation needed).

class moire_metrology.interfaces.Interface(name, bottom, top, gsfe_coeffs, stacking_func=None, reference='')[source]

Bases: object

The registry-dependent stacking interaction between two adjacent layers.

GSFE characterizes a pair of materials in contact, not a single material. An Interface couples two Material instances with the GSFE Fourier coefficients that describe their stacking energy landscape, plus (for homobilayer interfaces) the natural Bernal stacking convention used when multiple layers are stacked within a flake.

Parameters:

  • name (str) – Human-readable name (e.g. “Graphene/Graphene”, “MoSe2/WSe2 (H)”).

  • bottom (Material) – The lower layer of the pair. Boundary conditions like fix_bottom clamp the bottommost layer of the bottom flake.

  • top (Material) – The upper layer of the pair. May equal bottom (homobilayer interface) or differ (heterointerface).

  • gsfe_coeffs (tuple[float, ...]) – GSFE Fourier coefficients (c0, c1, c2, c3, c4, c5) in meV per unit cell, in the Carr basis (see module docstring). For centrosymmetric homobilayers, c4 = c5 = 0.

  • stacking_func (callable or None) – (k: int) -> (I, J) returning the natural Bernal stacking offsets for layer k within a flake of more than one layer of this same material. None for heterointerfaces (where no Bernal convention applies) and for monolayer flake usage.

  • reference (str) – Free-form citation string for the GSFE numbers, so the provenance travels with the data.

property is_homobilayer: bool

True iff the bottom and top materials are the same object.

classmethod from_dict(data)[source]

Build an Interface from a plain dict (e.g. parsed from TOML).

The dict must contain the following keys:

  • name (str)

  • bottom (dict — inline Material spec, see Material.from_dict())

  • top (dict — inline Material spec)

  • gsfe_coeffs (sequence of 6 floats, Carr basis, meV/uc)

  • reference (optional str)

stacking_func is intentionally NOT loadable from a plain dict — it is a Python callable, and serializing arbitrary callables is out of scope for the TOML loader. If you need a non-trivial stacking convention, construct the Interface directly in Python or post-process the loaded instance.

Extra keys raise a ValueError so typos in user TOML files surface immediately.

Return type:

Interface

Parameters:

data (Mapping[str, Any])

classmethod from_toml(path)[source]

Load an Interface from a TOML file.

The file must contain a top-level [interface] table with inline Material definitions for both layers:

[interface]
name = "MoSe2/WSe2 (H-stacked)"
gsfe_coeffs = [42.6, 16.0, -2.7, -1.1, 3.7, 0.6]
reference = "Shabani et al., Nat. Phys. 17, 720 (2021)"

[interface.bottom]
name = "WSe2"
lattice_constant = 0.3282
bulk_modulus = 43113.0
shear_modulus = 30770.0

[interface.top]
name = "MoSe2"
lattice_constant = 0.3288
bulk_modulus = 40521.0
shear_modulus = 26464.0

See Interface.from_dict() for the field semantics. The stacking_func field cannot be loaded from TOML; if you need a non-trivial stacking convention for multi-layer flakes, construct the Interface directly in Python instead.

Return type:

Interface

Parameters:

path (str | Path)

Lattice and geometry

Hexagonal lattice geometry and moire superlattice calculations.

Conventions:

  • Angles in degrees externally, radians internally

  • Lattice constants in nm

  • b1, b2 are unit vectors at angles theta0 and theta0 + 60 degrees

moire_metrology.lattice.rotation_matrix(angle_deg)[source]

2x2 rotation matrix for a given angle in degrees.

Return type:

ndarray

Parameters:

angle_deg (float)

class moire_metrology.lattice.HexagonalLattice(alpha, theta0=0.0)[source]

Bases: object

A 2D hexagonal lattice.

Parameters:

  • alpha (float) – Lattice constant in nm.

  • theta0 (float) – Orientation angle of the lattice in degrees (angle of b1 from x-axis).

property b1: ndarray

First lattice unit vector (length 1).

property b2: ndarray

Second lattice unit vector (length 1), at 60 degrees from b1.

property basis_matrix: ndarray

B = alpha * [b1 | b2], shape (2, 2). Columns are lattice vectors.

property reciprocal_matrix: ndarray

M = (2*pi/alpha) * inv([b1 | b2]), shape (2, 2).

Maps real-space coordinates to stacking phase coordinates:

[v; w] = M @ [x; y]

property unit_cell_area: float

Unit cell area in nm^2.

class moire_metrology.lattice.MoireGeometry(lattice, theta_twist, delta=0.0)[source]

Bases: object

Geometry of a moire superlattice from two hexagonal layers.

Parameters:

  • lattice (HexagonalLattice) – The reference lattice (layer 2 / substrate).

  • theta_twist (float) – Twist angle in degrees between the two layers.

  • delta (float) – Lattice mismatch: (alpha_layer1 / alpha_layer2) - 1. For homostructures (e.g., TBLG), delta = 0.

property R_twist: ndarray

Rotation matrix R(-theta_twist).

property moire_matrix: ndarray

Mr = R(+theta) - (1+delta)*I, shape (2,2).

The moire vectors satisfy: [V1|V2] = inv(Mr) * B.

Using R(+theta) (not R(-theta)) ensures that the stacking phases jump by exact integer multiples of 2*pi across V1 and V2, so the GSFE tiles perfectly over the moire unit cell.

property moire_vectors: tuple[ndarray, ndarray]

Moire superlattice vectors V1, V2 in nm.

property wavelength: float

Moire wavelength in nm (geometric mean of |V1| and |V2|).

property Mu1: ndarray

Conversion matrix for layer 1 (twisted layer) stacking phases.

Maps displacement (ux, uy) to stacking phase change via:

delta_v = -Mu1[0,:] @ [ux, uy] delta_w = -Mu1[1,:] @ [ux, uy]

property Mu2: ndarray

Conversion matrix for layer 2 (substrate) stacking phases.

For the substrate layer, displacement shifts stacking phase as:

delta_v = Mu2[0,:] @ [ux, uy] delta_w = Mu2[1,:] @ [ux, uy]

stacking_phases(x, y, ux=None, uy=None)[source]

Compute stacking phase coordinates (v, w) at positions (x, y).

Parameters:

  • x (ndarray) – Spatial coordinates in nm.

  • y (ndarray) – Spatial coordinates in nm.

  • ux (ndarray or None) – Displacement field components. If None, assumed zero (unrelaxed).

  • uy (ndarray or None) – Displacement field components. If None, assumed zero (unrelaxed).

Returns:

v, w – Stacking phase coordinates.

Return type:

ndarray

Mesh

Triangular mesh generation for moire unit cells.

Generates a regular grid on a parallelogram domain spanned by the moire vectors V1, V2 (or multiples thereof), then triangulates with Delaunay.

class moire_metrology.mesh.MoireMesh(points, triangles, V1, V2, ns, nt, n_scale, is_periodic=True, _boundary_info=None)[source]

Bases: object

Triangular mesh on a moire unit cell parallelogram (or finite domain).

Parameters:
points

Vertex coordinates (x, y) in nm.

Type:

ndarray, shape (2, Nv)

triangles

Triangle connectivity (vertex indices).

Type:

ndarray, shape (Nt, 3)

V1, V2

Parallelogram edge vectors in nm. For a periodic mesh these are the moire lattice vectors. For a finite mesh they describe the bounding parallelogram of the domain (used by some plotting helpers).

Type:

ndarray, shape (2,)

ns, nt

Number of grid divisions along V1 and V2.

Type:

int

is_periodic

True if the triangulation uses periodic wrap-around connectivity (the default, produced by MoireMesh.generate). False if the mesh covers a finite domain with open boundaries (produced by generate_finite_mesh). The discretization uses this flag to decide whether to apply lattice-vector wrap corrections to relative vertex positions when assembling differentiation matrices: periodic meshes need them, finite meshes don’t.

Type:

bool

classmethod generate(geometry, pixel_size=0.2, n_scale=1, min_points=100, max_points=500000)[source]

Generate a mesh on the moire unit cell parallelogram.

Parameters:

  • geometry (MoireGeometry) – Defines V1, V2 moire vectors.

  • pixel_size (float) – Target element size in nm.

  • n_scale (int) – Number of moire unit cells along each direction.

  • min_points (int) – Bounds on the number of grid points per direction.

  • max_points (int) – Bounds on the number of grid points per direction.

Return type:

MoireMesh

parametric_coords()[source]

Return the (s, t) parametric coordinates for each vertex.

s, t in [0, 1), where r = s*V1 + t*V2.

Return type:

tuple[ndarray, ndarray]

get_boundary_vertices(tol=1e-10)[source]

Identify vertices on each boundary of the parallelogram.

Returns dict with keys: ‘s0’, ‘s1’, ‘t0’, ‘t1’ (edges at s=0, s~1, t=0, t~1), and ‘corners’ (vertices at corners).

Return type:

dict

Parameters:

tol (float)

moire_metrology.mesh.generate_finite_mesh(geometry, n_cells=3, pixel_size=0.5, *, n_cells_x=None, n_cells_y=None)[source]

Generate a non-periodic mesh covering multiple moire unit cells.

Creates a parallelogram domain spanning n_cells moire vectors in each direction (or n_cells_x along V1 and n_cells_y along V2 when a rectangular extent is needed), with no periodic wrapping. Suitable for constrained relaxation.

Parameters:

  • geometry (MoireGeometry) – Moire geometry.

  • n_cells (int) – Number of moire unit cells along each direction. Used as the default for both axes when n_cells_x / n_cells_y aren’t given.

  • pixel_size (float) – Target element size in nm.

  • n_cells_x (int or None) – Optional overrides for the number of moire cells along V1 and V2 respectively. If either is None, falls back to n_cells. Use these to build a rectangular finite domain that’s wider in one direction than the other (e.g. for a twist-gradient demo where you want many wavelengths along the gradient axis but only a few perpendicular to it).

  • n_cells_y (int or None) – Optional overrides for the number of moire cells along V1 and V2 respectively. If either is None, falls back to n_cells. Use these to build a rectangular finite domain that’s wider in one direction than the other (e.g. for a twist-gradient demo where you want many wavelengths along the gradient axis but only a few perpendicular to it).

Returns:

Non-periodic mesh (boundary triangles do NOT wrap).

Return type:

MoireMesh

moire_metrology.mesh.generate_hex_periodic_mesh(geometry, pixel_size=2.0, eps_boundary_frac=0.05)[source]

Generate a hexagonal Wigner-Seitz periodic mesh centered on the AA stacking site (origin in our convention).

The mesh has exact 6-fold rotational symmetry around the center: vertex positions are placed on a triangular sub-lattice through origin with spacing vectors V1/N, V2/N (so the sub-lattice has the same orientation as the moire lattice), then filtered to those inside the WS hexagon. The 6 WS corners are added explicitly.

Use this together with the boundary identification helper identify_hex_periodic_boundary() and PeriodicPairConstraint (3 pair sets, one for each pair of opposite edges) plus a PinnedConstraints pinning the 6 corners at U=0 — see the online docs/example for a complete recipe.

Compared to generate_finite_mesh() (parallelogram supercell, C2 symmetry only), the hexagonal cell respects C3 symmetry around the AA site and the C3 orbit of AB / BA corner pins. This eliminates the direction-selection artifact that parallelogram supercells exhibit at low twist angles in problems with rich GSFE structure.

Parameters:

  • geometry (MoireGeometry) – Moire geometry defining V1, V2.

  • pixel_size (float) – Target vertex spacing in nm. The actual spacing is set by rounding |V1| / pixel_size to the nearest integer N, then using spacing |V1| / N (so the triangular sub-lattice tiles the cell exactly).

  • eps_boundary_frac (float) – Tolerance for the “inside hex” test, in units of pixel_size. Vertices within this tolerance of the boundary are kept.

Returns:

Hexagonal-cell mesh with is_periodic=False (open boundary; periodicity is enforced via PeriodicPairConstraint after boundary identification). V1, V2 stored as the moire vectors (NOT the rhombus-cell vectors), so downstream code computing stacking phases from positions uses the correct geometry.

Return type:

MoireMesh

moire_metrology.mesh.identify_hex_periodic_boundary(mesh, tol=1e-06)[source]

Identify boundary topology of a hexagonal-cell mesh built by generate_hex_periodic_mesh().

Returns a dict with:

  • "corners" — shape (6,) int array, mesh vertex indices for the 6 WS corners, ordered c1, c2, c3, c4, c5, c6 (counterclockwise).

  • "edges" — dict mapping k ∈ {1, …, 6} to the array of vertex indices on edge k, sorted along the c_{k} → c_{k+1} tangent direction. Each edge includes its two endpoint corners.

  • "pairs" — list of three dicts, one per opposite-edge pair: {"source": k_src, "dest": k_dst, "translation": v, "src_indices": array, "dst_indices": array}, where dst_indices[i] ≈ src_indices[i] + translation in the parent moire lattice.

Edge labelling (counterclockwise, edge k from c_k to c_{k+1}):

k=1: c1→c2, outward normal +V2, pair with k=4 (translation -V2) k=2: c2→c3, outward normal V2-V1, pair with k=5 (-(V2-V1)) k=3: c3→c4, outward normal -V1, pair with k=6 (+V1) k=4: c4→c5, outward normal -V2, k=5: c5→c6, outward normal V1-V2, k=6: c6→c1, outward normal +V1.

The three pairs (1↔4, 2↔5, 3↔6) cover ALL boundary vertices including corners (each corner is on two edges, so appears in two pair lists).

Return type:

dict

Parameters:
moire_metrology.mesh.generate_custom_mesh(geometry, *, outer_boundary='disk', outer_radius=70.0, center=(0.0, 0.0), holes=None, max_area=2.0, min_angle=25.0, boundary_density=None)[source]

Generate a mesh on an arbitrary 2D domain using meshpy.

Produces an unstructured triangulation with smooth boundaries, suitable for constrained relaxation on non-rectangular domains and domains with interior holes.

Parameters:

  • geometry (MoireGeometry) – Moire geometry (provides V1, V2 for metadata).

  • outer_boundary ("disk" or list of (x, y) tuples) – Outer domain boundary. "disk" creates a circle of radius outer_radius centred at center.

  • outer_radius (float) – Radius in nm (used only when outer_boundary="disk").

  • center (tuple of float) – Centre of the disk in nm (default origin).

  • holes (list of dict, optional) – Each dict describes a hole: {"center": (cx, cy), "radius": r} for a circular hole, or {"points": [(x1,y1), ...]} for an arbitrary polygon.

  • max_area (float) – Maximum triangle area in nm² (controls mesh density).

  • min_angle (float) – Minimum triangle angle in degrees (controls quality).

  • boundary_density (int or None) – Number of boundary points per circle. If None, chosen automatically from outer_radius / sqrt(max_area).

Returns:

Non-periodic mesh with is_periodic=False.

Return type:

MoireMesh

Examples

Disk with a central hole:

mesh = generate_custom_mesh(
    geom,
    outer_boundary="disk",
    outer_radius=70.0,
    holes=[{"center": (0, 0), "radius": 10.0}],
    max_area=1.5,
)

Hourglass (two pads connected by a neck):

# Define as a polygon
mesh = generate_custom_mesh(
    geom,
    outer_boundary=hourglass_points,
    max_area=1.0,
)

GSFE surface

Generalized Stacking Fault Energy (GSFE) surface and its derivatives.

The GSFE is parameterized as a 6-coefficient Fourier expansion in stacking phase coordinates (v, w):

V(v,w) = c0 + c1*(cos(v) + cos(w) + cos(v+w))

  • c2*(cos(v+2w) + cos(v-w) + cos(2v+w))

  • c3*(cos(2v) + cos(2w) + cos(2v+2w))

  • c4*(sin(v) + sin(w) - sin(v+w))

  • c5*(sin(2v+2w) - sin(2v) - sin(2w))

All methods operate element-wise on numpy arrays.

class moire_metrology.gsfe.GSFESurface(coeffs)[source]

Bases: object

GSFE energy surface with analytical derivatives.

Parameters:

coeffs (sequence of float) – (c0, c1, c2, c3, c4, c5) in meV/unit cell.

dv(v, w)[source]

dV/dv.

Return type:

ndarray

Parameters:
dw(v, w)[source]

dV/dw.

Return type:

ndarray

Parameters:
d2v2(v, w)[source]

d^2V/dv^2.

Return type:

ndarray

Parameters:
d2w2(v, w)[source]

d^2V/dw^2.

Return type:

ndarray

Parameters:
d2vw(v, w)[source]

d^2V/dvdw.

Return type:

ndarray

Parameters:
property minimum_value: float

Minimum of V over all (v, w). Cached after first call.

saddle_point_energy()[source]

Energy at the saddle point (SP stacking).

For the hexagonal GSFE, the SP lies along v = w. We maximize V(v, v) in the range [2*pi/3, 4*pi/3] to find the saddle point between adjacent AB/BA domains.

Return type:

float

Energy functional

Energy functional assembly for moire relaxation.

Computes the total energy E = E_elastic + E_GSFE, its gradient, and Hessian.

The elastic Hessian is constant (precomputed once as sparse matrix). The GSFE Hessian is vertex-diagonal (cheap per iteration).

For multi-layer systems, GSFE coupling exists at: - The interface between the two stacks (always present) - Between adjacent layers within each stack (intra-stack, when nlayer > 1)

class moire_metrology.energy.RelaxationEnergy(disc, conv, geometry, gsfe_interface, K1, G1, K2, G2, nlayer1=1, nlayer2=1, gsfe_flake1=None, gsfe_flake2=None, I1_vect=None, J1_vect=None, I2_vect=None, J2_vect=None, constraints=None, elastic_strain='cauchy')[source]

Bases: object

Total energy functional for moire relaxation.

Provides: - __call__(U) -> (E, grad) - hessian(U) -> sparse matrix - hessp(U, p) -> Hessian-vector product

Parameters:
hessp(U, p)[source]

Compute Hessian-vector product H(U) @ p.

If constrained, U and p are in free-DOF space; result is also free-DOF space.

Return type:

ndarray

Parameters:
hessian(U, modified=False)[source]

Compute the Hessian. If constrained, returns the free x free subblock.

Parameters:

  • U (ndarray) – Current DOF vector (free-space if constrained).

  • modified (bool) – If True, flip negative eigenvalues of the per-vertex GSFE Hessian blocks so the returned matrix is positive definite (modulo the elastic part, which is always PSD). This enables Newton steps at non-minimum stacking sites (AA, saddle points) where the unmodified GSFE Hessian is indefinite.

Return type:

csr_matrix

energy_maps(U)[source]

Compute spatially resolved energy density maps.

Return type:

dict

Parameters:

U (ndarray)

Discretization

FEM discretization on triangular meshes.

Provides differentiation matrices, area weights, and conversion matrices for mapping the flat solution vector to per-layer displacement fields.

For a periodic mesh, all vertices are independent (periodicity is encoded in the triangulation via wrapping indices). For constrained problems, PinnedConstraints tracks which DOFs are free vs pinned.

class moire_metrology.discretization.PinnedConstraints(free_indices, pinned_indices, pinned_values, n_free, n_full)[source]

Bases: object

Tracks which DOFs are free vs pinned in a constrained relaxation.

The full solution vector has length n_full. The optimizer only sees the n_free free DOFs. This class handles the mapping between the two.

Parameters:
free_indices

Indices into the full U vector that are free (optimized).

Type:

ndarray of int

pinned_indices

Indices into the full U vector that are pinned (fixed).

Type:

ndarray of int

pinned_values

Displacement values at pinned DOFs.

Type:

ndarray

n_free
Type:

int

n_full
Type:

int

expand(U_free)[source]

Map free DOFs to full vector, inserting pinned values.

Return type:

ndarray

Parameters:

U_free (ndarray)

expand_zeros(v_free)[source]

Map free DOFs to full vector with zeros at pinned positions.

Used for perturbation directions (Hessian-vector products) where pinned DOFs should not contribute.

Return type:

ndarray

Parameters:

v_free (ndarray)

project(vec_full)[source]

Extract free-DOF entries from a full-length vector (e.g., gradient).

Return type:

ndarray

Parameters:

vec_full (ndarray)

project_hessian(H_full)[source]

Extract the free x free subblock of a sparse Hessian.

Return type:

csr_matrix

Parameters:

H_full (spmatrix)

moire_metrology.discretization.build_outer_layer_constraints(conv, fix_top, fix_bottom, pin_mean=False)[source]

Build PinnedConstraints clamping the outer (free-surface) layers to zero.

For a multi-layer stack, the “top” outer layer is the topmost layer of stack 1 (the layer farthest from the twisted interface, on the upper side of the heterostructure). The “bottom” outer layer is the deepest layer of stack 2 (the layer farthest from the interface on the lower side, i.e. the substrate’s free surface in a typical setup).

Pinning these to zero displacement simulates the layers being bonded to a rigid bulk — the conventional way to approximate a semi-infinite substrate with a finite stack. This matches the fix_external_layers parameter in the MATLAB paper code.

Parameters:

  • conv (ConversionMatrices) – Conversion matrices describing the multi-layer DOF layout.

  • fix_top (bool) – Pin all DOFs (ux and uy at every vertex) of stack 1, layer 0.

  • fix_bottom (bool) – Pin all DOFs of stack 2, layer (nlayer2 - 1).

  • pin_mean (bool)

Returns:

Suitable to pass to RelaxationSolver.solve(constraints=…).

Return type:

PinnedConstraints

class moire_metrology.discretization.ConversionMatrices(conv_x1, conv_y1, conv_x2, conv_y2, n_sol, n_vertices, nlayer1, nlayer2)[source]

Bases: object

Matrices mapping the flat solution vector U to per-layer displacements.

For a system with nlayer1 layers in stack 1 and nlayer2 layers in stack 2, the solution vector has length 2 * (nlayer1 + nlayer2) * Nv (x and y displacements for each layer at each vertex).

The conversion matrices extract the x/y displacement for each layer and apply the appropriate strain partition factor (eps for stack 1, 1-eps for stack 2).

Parameters:
class moire_metrology.discretization.Discretization(mesh, geometry)[source]

Bases: object

FEM infrastructure for a triangular moire mesh.

Handles BOTH periodic and finite (open) meshes. The mesh’s is_periodic flag controls whether the diff matrix builder applies lattice-vector wrap corrections to relative vertex positions: periodic meshes need them (the structured triangulation contains wrap-around triangles whose vertices are placed on opposite sides of the parallelogram), finite meshes don’t.

Parameters:

  • mesh (MoireMesh) – The triangular mesh.

  • geometry (MoireGeometry) – Moire geometry for computing stacking phases.

property diff_mat_x: csr_matrix

maps vertex values to per-triangle x-derivatives.

Type:

Sparse matrix (Nt, Nv)

property diff_mat_y: csr_matrix

maps vertex values to per-triangle y-derivatives.

Type:

Sparse matrix (Nt, Nv)

property triangle_areas: ndarray

Area of each triangle in nm^2, shape (Nt,).

property vertex_areas: ndarray

Area associated with each vertex (Voronoi dual), shape (Nv,).

Computed by distributing each triangle’s area equally among its 3 vertices.

property total_area: float

Total mesh area in nm^2.

property triangle_to_vertex: csr_matrix

maps triangle values to vertex values by area-weighted average.

Type:

Sparse matrix (Nv, Nt)

build_conversion_matrices(nlayer1=1, nlayer2=1, eps=0.5)[source]

Build matrices mapping solution vector U to per-layer displacements.

The solution vector U has the layout:
[ux_layer1_0, …, ux_layer1_{n1-1},

ux_layer2_0, …, ux_layer2_{n2-1}, uy_layer1_0, …, uy_layer1_{n1-1}, uy_layer2_0, …, uy_layer2_{n2-1}]

Each layer block has Nv entries (one per vertex).

For a bilayer (nlayer1=nlayer2=1), the strain partition factor eps controls how the total relative displacement is shared: - Layer 1 displacement = eps * u - Layer 2 displacement = (1-eps) * u With eps=0.5, both layers move symmetrically.

Parameters:

  • nlayer1 (int) – Number of layers in each stack.

  • nlayer2 (int) – Number of layers in each stack.

  • eps (float) – Strain partition factor (0 to 1). Default 0.5 (symmetric).

Return type:

ConversionMatrices

moire_metrology.discretization.PeriodicDiscretization

alias of Discretization

Pinning

Interactive and programmatic tools for defining pinning sites.

Provides three ways to define where known stacking configurations are pinned:

  1. Programmatic — PinningMap with pin_stacking() calls

  2. Interactive — InteractivePinner with matplotlib click interface

  3. From file — PinningMap.from_csv()

Example:

from moire_metrology.pinning import PinningMap from moire_metrology import GRAPHENE, RelaxationSolver, SolverConfig from moire_metrology.lattice import HexagonalLattice, MoireGeometry from moire_metrology.mesh import MoireMesh from moire_metrology.discretization import Discretization

geometry = MoireGeometry(HexagonalLattice(0.247), theta_twist=0.5) mesh = MoireMesh.generate(geometry, pixel_size=0.5)

pins = PinningMap(mesh, geometry) pins.pin_stacking(x=10.0, y=20.0, stacking=”AB”, radius=3.0) pins.pin_stacking(x=30.0, y=15.0, stacking=”BA”, radius=3.0)

disc = Discretization(mesh, geometry) conv = disc.build_conversion_matrices() constraints = pins.build_constraints(conv)

solver = RelaxationSolver(SolverConfig(display=True)) result = solver.solve(GRAPHENE, GRAPHENE, 0.5, constraints=constraints)

class moire_metrology.pinning.PinSite(x, y, stacking, radius)[source]

Bases: object

A single pinning site with known stacking.

Parameters:
class moire_metrology.pinning.PinningMap(mesh, geometry)[source]

Bases: object

Collection of pinning sites for constrained relaxation.

Parameters:

  • mesh (MoireMesh) – The computational mesh.

  • geometry (MoireGeometry) – Moire geometry for computing stacking phases.

pin_stacking(x, y, stacking, radius=1.0)[source]

Pin vertices near (x, y) to a known stacking configuration.

Parameters:

  • x (float) – Center of the pinning site (nm).

  • y (float) – Center of the pinning site (nm).

  • stacking (str) – Stacking type: ‘AA’, ‘AB’, or ‘BA’.

  • radius (float) – Radius of the pinning region (nm). All mesh vertices within this radius are pinned.

pin_vertices(vertex_indices, stacking)[source]

Pin specific vertices to a known stacking configuration.

Parameters:

  • vertex_indices (ndarray of int) – Mesh vertex indices to pin.

  • stacking (str) – Stacking type: ‘AA’, ‘AB’, or ‘BA’.

get_pinned_vertex_indices()[source]

Get all pinned vertex indices (union of all pin sites).

Return type:

ndarray

build_constraints(conv, nlayer1=1, nlayer2=1)[source]

Build PinnedConstraints from the current pin sites.

Pins are applied to the interface layers (innermost of each stack). Both interface layers are pinned to produce the target stacking.

Parameters:

  • conv (ConversionMatrices) – From disc.build_conversion_matrices().

  • nlayer1 (int) – Number of layers in each stack (needed for DOF offset calculation).

  • nlayer2 (int) – Number of layers in each stack (needed for DOF offset calculation).

Return type:

PinnedConstraints

classmethod from_csv(path, mesh, geometry)[source]

Load pin sites from a CSV file.

Expected format (with header): x, y, stacking, radius

Return type:

PinningMap

Parameters:
save_csv(path)[source]

Save pin sites to a CSV file.

Parameters:

path (str | Path)

class moire_metrology.pinning.InteractivePinner(mesh, geometry, background_image=None, image_extent=None)[source]

Bases: object

Interactive matplotlib tool for placing pinning sites on a moire mesh.

Usage:

pinner = InteractivePinner(mesh, geometry) pinner.show() # opens interactive window # Left-click to place pins, right-click to remove nearest # Press 1/2/3 to select AA/AB/BA stacking type # Press +/- to adjust radius # Close window when done constraints = pinner.get_constraints(conv)

Parameters:

  • mesh (MoireMesh)

  • geometry (MoireGeometry)

  • background_image (str or ndarray or None) – Path to an image file (e.g., PFM/STM data) to display behind the mesh, or a 2D numpy array. If None, shows the unrelaxed GSFE map.

  • image_extent (tuple or None) – (xmin, xmax, ymin, ymax) for the background image in nm.

show()[source]

Open the interactive pinning window.

get_constraints(conv, nlayer1=1, nlayer2=1)[source]

Build constraints from the interactively placed pins.

Return type:

PinnedConstraints

Parameters:

Plotting

Visualization utilities for moire relaxation results.

Uses matplotlib’s Triangulation and tripcolor for efficient rendering of scalar fields on triangular meshes, with optional tiling to show multiple moire unit cells.

moire_metrology.plotting.plot_scalar_field(mesh, values, ax=None, n_tile=1, title='', cmap='viridis', colorbar=True, **kwargs)[source]

Plot a scalar field on the moire mesh.

Parameters:

  • mesh (MoireMesh) – The triangular mesh.

  • values (ndarray, shape (Nv,)) – Scalar values at each vertex.

  • ax (matplotlib Axes or None) – If None, creates a new figure.

  • n_tile (int) – Number of periodic copies in each direction (for visualization).

  • title (str) – Plot title.

  • cmap (str) – Colormap name.

  • colorbar (bool) – Whether to add a colorbar.

  • **kwargs – Passed to tripcolor.

Return type:

Axes

moire_metrology.plotting.plot_displacement_field(mesh, ux, uy, ax=None, n_tile=1, scale=1.0, **kwargs)[source]

Plot displacement vectors on the mesh.

Parameters:

  • mesh (MoireMesh) – The triangular mesh.

  • ux (ndarray, shape (Nv,)) – Displacement components.

  • uy (ndarray, shape (Nv,)) – Displacement components.

  • ax (matplotlib Axes or None)

  • n_tile (int) – Periodic tiling.

  • scale (float) – Arrow scale factor.

Return type:

Axes

Strain extraction

Strain extraction from moire patterns.

Two workflows:

  1. Pointwise extraction — given moire vector observables (lambda1, lambda2, phi1, phi2) for a single point:

    from moire_metrology.strain import get_strain_minimize_compression
result = get_strain_minimize_compression(
    alpha1, alpha2, lambda1, lambda2, phi1_deg, phi2_deg,
)

  2. Spatial strain mapping — from traced moire fringe polylines:

    from moire_metrology.strain import (
    FringeSet, compute_strain_field, compute_displacement_field,
    convex_hull_mask,
)
fringes = FringeSet.from_matlab("points.mat")
I_field, J_field = fringes.fit_registry_fields(order=11)
strain = compute_strain_field(
    x, y, I_field, J_field, alpha1, alpha2, phi0_deg=-65.6,
)
ux, uy = compute_displacement_field(
    x, y, I_field, J_field, geometry, target_stacking="BA",
)
class moire_metrology.strain.StrainResult(theta_twist, phi0, S11, S12, S22, eps_c, eps_s, eps1, eps2, strain_angle)[source]

Bases: object

Result of strain extraction from moire observables.

Parameters:
theta_twist

Twist angle in degrees.

Type:

float

phi0

Substrate lattice orientation in degrees.

Type:

float

S11, S12, S22

Strain tensor components (symmetric: S21 = S12).

Type:

float

eps_c

Volumetric (compression) strain = (eps1 + eps2) / 2.

Type:

float

eps_s

Shear (deviatoric) strain = (eps1 - eps2) / 2.

Type:

float

eps1, eps2

Principal strains (eigenvalues of strain tensor).

Type:

float

strain_angle

Principal strain axis angle in degrees.

Type:

float

moire_metrology.strain.get_strain(alpha1, alpha2, lambda1, lambda2, phi1_deg, phi2_deg, phi0)[source]

Extract twist angle and strain from moire observables at a given phi0.

Parameters:

  • alpha1 (float) – Lattice constants of the two layers (nm).

  • alpha2 (float) – Lattice constants of the two layers (nm).

  • lambda1 (float) – Moire periods along the two moire vectors (nm).

  • lambda2 (float) – Moire periods along the two moire vectors (nm).

  • phi1_deg (float) – Orientations of the two moire vectors (degrees).

  • phi2_deg (float) – Orientations of the two moire vectors (degrees).

  • phi0 (float) – Substrate lattice orientation angle (degrees).

Return type:

StrainResult

moire_metrology.strain.get_strain_minimize_compression(alpha1, alpha2, lambda1, lambda2, phi1_deg, phi2_deg, phi0_guess=0.0)[source]

Extract strain with phi0 chosen to minimize compression strain |eps_c|.

This solves analytically for phi0 such that eps_c = 0 when possible, or finds the closest achievable phi0 when eps_c = 0 is infeasible.

Parameters:

  • alpha1 (float) – Lattice constants (nm).

  • alpha2 (float) – Lattice constants (nm).

  • lambda1 (float) – Moire periods (nm).

  • lambda2 (float) – Moire periods (nm).

  • phi1_deg (float) – Moire vector orientations (degrees).

  • phi2_deg (float) – Moire vector orientations (degrees).

  • phi0_guess (float) – Initial guess for phi0 (degrees), used for disambiguation.

Return type:

StrainResult

moire_metrology.strain.get_strain_axis(S11, S12, S22)[source]

Diagonalize the strain tensor to get principal strains.

Parameters:

  • S11 (float) – Strain tensor components (S is symmetric: S21 = S12).

  • S12 (float) – Strain tensor components (S is symmetric: S21 = S12).

  • S22 (float) – Strain tensor components (S is symmetric: S21 = S12).

Return type:

tuple[float, float, float]

Returns:

  • eps1 (float) – Principal strain 1 (larger).

  • eps2 (float) – Principal strain 2 (smaller).

  • strain_angle (float) – Principal strain axis angle in degrees.

moire_metrology.strain.shear_strain_invariant(alpha1, alpha2, lambda1, lambda2, phi1_deg, phi2_deg)[source]

Compute the phi0-independent shear strain from moire observables.

This is the closed-form formula from the paper:

eps_s = alpha * r_minus / (2 * lambda1 * lambda2 * sin(dphi))

where r_minus = sqrt(l1^2 + l2^2 - 2*l1*l2*cos(dphi - pi/3))

This quantity depends only on the moire observables, not on phi0.

Parameters:

  • alpha1 (float) – Lattice constants (nm).

  • alpha2 (float) – Lattice constants (nm).

  • lambda1 (float) – Moire periods (nm).

  • lambda2 (float) – Moire periods (nm).

  • phi1_deg (float) – Moire vector orientations (degrees).

  • phi2_deg (float) – Moire vector orientations (degrees).

Returns:

Shear strain magnitude (always non-negative).

Return type:

float

moire_metrology.strain.compute_displacement_field(x, y, I_field, J_field, geometry, target_stacking='BA', remove_mean=True)[source]

Build a displacement field from polynomial registry fits.

The pointwise pinning machinery in moire_metrology.pinning solves Mu @ u = -(v_target - v0) at individual sites to put a given mesh vertex at a chosen stacking. This function does the same solve at every query point (x, y), but with v_target coming from a continuous polynomial fit of the registry fields rather than discrete pin assignments. The result is a smooth displacement field suitable as the initial condition for relaxation against an experimentally measured moiré pattern.

The target phase at each point is

v_target(r) = 2π · I(r) + v_offset w_target(r) = 2π · J(r) + w_offset

where (v_offset, w_offset) is the constant phase shift from moire_metrology.pinning.STACKING_PHASES that puts integer (I, J) data sites at the requested target_stacking (e.g. "BA" is the global GSFE minimum for H-stacked TMDs).

Parameters:

  • x (ndarray) – Query coordinates (nm), of any common shape.

  • y (ndarray) – Query coordinates (nm), of any common shape.

  • I_field (RegistryField) – Polynomial fits of the moire registry index fields.

  • J_field (RegistryField) – Polynomial fits of the moire registry index fields.

  • geometry (MoireGeometry) – The average moire geometry. Provides stacking_phases and Mu1 for the per-point linear solve.

  • target_stacking (str) – One of "AA", "AB", "BA". Default "BA".

  • remove_mean (bool) – If True, subtract the mean of the resulting displacement so the global translation gauge is fixed at zero. Default True.

Returns:

ux, uy – Displacement field components (nm), of the same shape as x.

Return type:

ndarray

Notes

The polynomial fits are unconstrained outside the convex hull of their training data and grow rapidly there. Mask the query points with convex_hull_mask() before calling this function on a mesh that extends past the data extent.

moire_metrology.strain.compute_strain_field(x, y, I_field, J_field, alpha1, alpha2, phi0_deg)[source]

Compute the spatially-varying strain field from registry polynomials.

Implements the spatial extension of the ACS Nano paper strain inversion. At each query point, the local moiré reciprocal lattice vectors are obtained from the polynomial fit gradients via eq. 9:

[v1(r), v2(r)] = inv([∇I(r); ∇J(r)])

yielding (λ1, λ2, φ1, φ2) at every point. The pointwise get_strain() is then called with a global phi0_deg to recover (θ, ε_c, ε_s) per point.

Parameters:

  • x (ndarray) – Query coordinates (nm), of any common shape.

  • y (ndarray) – Query coordinates (nm), of any common shape.

  • I_field (RegistryField) – Polynomial fits of the moire registry index fields, with analytic dx / dy evaluators (e.g. moire_metrology.strain.RegistryField).

  • J_field (RegistryField) – Polynomial fits of the moire registry index fields, with analytic dx / dy evaluators (e.g. moire_metrology.strain.RegistryField).

  • alpha1 (float) – Lattice constants of the two layers (nm). Convention follows get_strain() — alpha2 is the larger of the two for the H-MoSe2/WSe2 sample reproduced in the paper.

  • alpha2 (float) – Lattice constants of the two layers (nm). Convention follows get_strain() — alpha2 is the larger of the two for the H-MoSe2/WSe2 sample reproduced in the paper.

  • phi0_deg (float) – Substrate lattice orientation angle (degrees), held constant across the field. The paper uses a global φ0 rather than a per-point one.

Returns:

Each entry has the same shape as x. Keys:

theta, eps_c, eps_s

Local twist angle (deg) and the volumetric / shear strain scalars used in the paper figures.

S11, S12, S22

Components of the symmetric strain tensor in the global (x, y) frame. S21 = S12. These are the natural inputs to a gradient-integration IC builder for the relaxation framework.

eps1, eps2, strain_angle

Principal strains and the principal-axis angle (deg).

lambda1, lambda2, phi1_deg, phi2_deg

The local moire vector lengths and orientations recovered from the eq. 9 inversion (intermediate quantities, exposed for diagnostics).

Return type:

dict

moire_metrology.strain.displacement_from_strain_field(discretization, *, theta_deg, theta_avg_deg, S11, S12, S22, pin_vertex=0, dirichlet_vertices=None, dirichlet_ux=None, dirichlet_uy=None)[source]

Build a displacement field from a target local twist + strain.

Translates the spatially-varying (θ(r), S(r)) output of compute_strain_field() into the relaxation framework’s native u(r) language. The framework holds a single global twist θ_avg fixed, so any local twist deviation δθ(r) = θ(r) − θ_avg and any local strain S(r) must show up as a gradient of u. Specifically (with the framework’s eps=0.5 layer partition, confirmed against energy.py and discretization.py):

∂u_x/∂x = S11(r) ∂u_y/∂y = S22(r) ∂u_x/∂y = S12(r) − δθ(r) ∂u_y/∂x = S12(r) + δθ(r)

These four equations are over-determined (four constraints on the four components of ∂u/∂r per point — exactly the count, but with a global compatibility condition that need not hold for noisy measured data). The implementation reads them as a sparse linear least-squares problem on the FEM mesh: [Dx; Dy] @ u_x = (S11; S12 − δθ) and similarly for u_y, where Dx and Dy are Discretization’s vertex-to-triangle gradient operators.

When dirichlet_vertices is given, those vertices are held at fixed (dirichlet_ux, dirichlet_uy) values and the Poisson solve optimizes the remaining vertices. This is the natural way to combine the gradient IC with pinning constraints: the pinned vertices get the displacement that the PinningMap computes for their target stacking, and the free vertices are interpolated to match the target gradient field. When no Dirichlet vertices are given, a single pin_vertex is pinned to zero to fix the global translation gauge.

Parameters:

  • discretization (Discretization) – FEM discretization carrying the mesh and the gradient operators.

  • theta_deg (ndarray, shape (Nv,)) – Local twist angle at each mesh vertex (degrees). Use np.abs(compute_strain_field(...)["theta"]) since the package convention reports -θ.

  • theta_avg_deg (float) – Global twist angle the relaxation framework holds fixed (deg).

  • S11 (ndarray, shape (Nv,)) – Strain tensor components at each mesh vertex (the compute_strain_field output keys of the same names).

  • S12 (ndarray, shape (Nv,)) – Strain tensor components at each mesh vertex (the compute_strain_field output keys of the same names).

  • S22 (ndarray, shape (Nv,)) – Strain tensor components at each mesh vertex (the compute_strain_field output keys of the same names).

  • pin_vertex (int) – Index of the vertex pinned to u = 0 to fix the global translation gauge. Ignored when dirichlet_vertices is set.

  • dirichlet_vertices (ndarray of int or None) – Vertex indices where u is prescribed. When given, dirichlet_ux and dirichlet_uy must also be provided.

  • dirichlet_ux (ndarray or None) – Prescribed displacement values at the Dirichlet vertices (nm).

  • dirichlet_uy (ndarray or None) – Prescribed displacement values at the Dirichlet vertices (nm).

Returns:

ux, uy – Per-vertex displacement field components (nm), suitable as the initial condition for a relaxation against the average-twist moiré geometry.

Return type:

ndarray, shape (Nv,)

class moire_metrology.strain.FringeLine(x, y, index, family)[source]

Bases: object

A single traced moire fringe line.

Parameters:
x

x-coordinates along the fringe (nm).

Type:

ndarray

y

y-coordinates along the fringe (nm).

Type:

ndarray

index

Integer registry index assigned to this fringe (I or J value).

Type:

int

family

Fringe family: 1 for I-fringes, 2 for J-fringes.

Type:

int

class moire_metrology.strain.FringeSet(fringes=<factory>)[source]

Bases: object

Collection of traced moire fringes for strain extraction.

Parameters:

fringes (list[FringeLine])

fringes

All traced fringe lines.

Type:

list of FringeLine

property i_fringes: list[FringeLine]

Fringes of family 1 (I-type).

property j_fringes: list[FringeLine]

Fringes of family 2 (J-type).

fit_registry_fields(order=11)[source]

Fit polynomial registry fields I(x,y) and J(x,y) to the fringe data.

The two families of fringes are fit independently, each using its raw polyline points (one constraint per traced point) — this matches the approach in the maintainer’s MATLAB pipeline that was validated against ACS Nano 16, 1471 (2022) Fig. 1.

Parameters:

order (int) – Polynomial order for the 2D fit (default 11, matching the paper Methods section).

Returns:

I_field, J_field – Polynomial fits to the I and J registry fields.

Return type:

RegistryField

estimate_moire_wavelength()[source]

Estimate the moire wavelength from fringe spacing.

Return type:

float

classmethod from_csv(path)[source]

Load fringes from a CSV file.

Expected format: x, y, index, family Each row is a point; fringes are separated by blank lines or by changes in the (index, family) pair.

Return type:

FringeSet

Parameters:

path (str | Path)

classmethod from_matlab(path)[source]

Load fringes from a MATLAB .mat file (GUI output format).

Expects keys: xpts_list, ypts_list, line_integer_val, line_type_list.

Return type:

FringeSet

Parameters:

path (str | Path)

class moire_metrology.strain.RegistryField(coeffs, order, x_center, y_center, x_scale, y_scale)[source]

Bases: object

A 2D polynomial fit of a moire registry field (I or J).

The polynomial has the form:

f(x,y) = sum_{i+j <= order} c_{ij} * x^i * y^j

Parameters:
coeffs

Polynomial coefficients in the order used by _poly_terms.

Type:

ndarray

order

Maximum polynomial order.

Type:

int

x_center, y_center

Coordinate centering for numerical stability.

Type:

float

x_scale, y_scale

Coordinate scaling for numerical stability.

Type:

float

dx(x, y)[source]

Evaluate df/dx at (x, y).

Return type:

ndarray

Parameters:
dy(x, y)[source]

Evaluate df/dy at (x, y).

Return type:

ndarray

Parameters:
classmethod fit(x, y, values, order=8)[source]

Fit a 2D polynomial to scattered data.

Parameters:

  • x (ndarray, shape (N,)) – Data point coordinates (nm).

  • y (ndarray, shape (N,)) – Data point coordinates (nm).

  • values (ndarray, shape (N,)) – Field values at each point (integer registry indices).

  • order (int) – Maximum polynomial order (default 8).

Return type:

RegistryField

moire_metrology.strain.convex_hull_mask(data_x, data_y, query_x, query_y)[source]

Boolean mask of query points inside the convex hull of input data.

Parameters:

  • data_x (ndarray) – Coordinates of the points used to fit a registry polynomial.

  • data_y (ndarray) – Coordinates of the points used to fit a registry polynomial.

  • query_x (ndarray) – Query coordinates (any common shape).

  • query_y (ndarray) – Query coordinates (any common shape).

Returns:

inside – Same shape as query_x. True where the point lies inside (or on the boundary of) the convex hull of the input data.

Return type:

ndarray of bool

Multi-layer

Multi-layer relaxation convenience APIs.

Example:

from moire_metrology.multilayer import LayerStack from moire_metrology import GRAPHENE

stack = LayerStack(top=GRAPHENE, n_top=5, bottom=GRAPHENE, n_bottom=3, theta_twist=0.5) result = stack.solve()

class moire_metrology.multilayer.LayerStack(moire_interface, top_interface=None, bottom_interface=None, n_top=1, n_bottom=1, theta_twist=0.0, delta=None, theta0=0.0)[source]

Bases: object

A heterostructure stack of two flakes.

Parameters:

  • moire_interface (Interface) – The twisted interface between the bottom layer of the top flake and the top layer of the bottom flake. Carries both materials (moire_interface.top, moire_interface.bottom) and the GSFE coefficients of the moire stacking.

  • top_interface (Interface or None) – Required iff n_top > 1. The homobilayer interface inside the top flake. Must satisfy top_interface.bottom is top_interface.top is moire_interface.top.

  • bottom_interface (Interface or None) – Required iff n_bottom > 1. The homobilayer interface inside the bottom flake. Must satisfy bottom_interface.bottom is bottom_interface.top is moire_interface.bottom.

  • n_top (int) – Number of layers in each flake.

  • n_bottom (int) – Number of layers in each flake.

  • theta_twist (float) – Twist angle between the two flakes (degrees).

  • delta (float or None) – Lattice mismatch. If None, computed from material lattice constants.

  • theta0 (float) – Lattice orientation angle (degrees).

property top

Top flake material (convenience accessor).

property bottom

Bottom flake material (convenience accessor).

solve(config=None, *, fix_top=False, fix_bottom=False, pin_mean=False)[source]

Run relaxation on this stack.

Parameters:

  • config (SolverConfig or None) – Solver configuration. If None, uses defaults.

  • fix_top (bool) – Pin all DOFs of the topmost layer (the topmost layer of the top flake) to zero. Use this to simulate a rigid bulk above the heterostructure.

  • fix_bottom (bool) – Pin all DOFs of the bottommost layer (the deepest layer of the bottom flake) to zero. Use this to simulate a rigid bulk below the heterostructure — typical for a twisted flake on a thick substrate (e.g. graphene on graphite).

  • pin_mean (bool)

Return type:

RelaxationResult

describe()[source]

Human-readable description of the stack.

Return type:

str