matscipy.fracture_mechanics.continuation

Force-only pseudo-arclength tracing of clamped crack clusters (lattice-trapping S-curves).

A clamped cluster (clusters.clamped_cluster) in a SinclairCrack with variable_alpha=False, energy=”total” is the problem: free atoms move, clamped atoms sit on the K-field. The tracer works in x = (u, k) with k = K/K_G and u the free-atom displacements relative to the field; the arclength metric is ds^2 = |du|^2 + alpha^2 dk^2. Correctors use a surrogate Jacobian (the Hessian of any matscipy many-body potential, e.g. Stillinger-Weber for Si, plus X-H springs, scaled by one calibrated factor) with good-Broyden updates from every force call, then Newton-Krylov.

Functions

`bordered_solve`(J, t, rhs_F, rhs_g)	Solve [J; t] d = [rhs_F; rhs_g] (the pseudo-arclength bordered system).
`trace_s_curve`(problem, surrogate, k0, *[, ...])	Trace the S-curve K(u) of a clamped crack cluster from load k0 = K/K_G.

Classes

`BroydenCorrector`([ftol, max_iter, step_limit])	Quasi-Newton on the bordered system F(y) = 0, t.(y - yp) = 0; B starts from the surrogate Jacobian and gets a good-Broyden rank-1 update from every force call.
`ClampedProblem`(crack, K_G[, atoms])	x = (u flat over free atoms, k = K/K_G) on a clamped SinclairCrack.
`CorrectorResult`(y, F, iters, calls, ok)
`HessianSurrogate`(problem, calc[, ks, kb, scale])	Approximate Jacobian dF/dx: the Hessian of calc (a matscipy many-body calculator with a hessian property) on the heavy atoms, times a calibrated scale, plus anisotropic X-H springs (stretch ks, bend kb, eV/A^2) for H caps (h_host).
`KrylovCorrector`([ftol, max_iter, eps_fd, ...])	Jacobian-free Newton-Krylov on the bordered system (GMRES; finite-difference J v), right preconditioned by an approximate Jacobian; at most max_calls force calls.
`TraceSettings`([ftol, alpha, ds0, ds_min, ...])

class matscipy.fracture_mechanics.continuation.ClampedProblem(crack, K_G, atoms=None)

Bases: object

x = (u flat over free atoms, k = K/K_G) on a clamped SinclairCrack. Free atoms: region 1 (groups != 0); clamped atoms ride the K-field. U is the field per unit k (all atoms).

atoms: the built cluster (clamped_cluster/cluster_like), whose geometry is the start; the crack’s own atoms at u = 0 put every H cap on its host (a cap’s reference is its host’s), so pass it whenever the cluster is passivated. k_build is the load the crack was set up at.

Attributes:

advance_break: Bond length (A) at which a crack-line bond counts as broken (settable).
line_bonds: Heavy-atom bonds crossing the cleavage plane at or ahead of the tip, sorted along the propagation direction: (i, j, reference length).

Methods

`advance`(x)	Crack advance in bonds: the sum over crack-line bonds of the fraction each is stretched from its reference length towards advance_break (clipped to [0, 1]).
`evaluate`(x)	(energy, forces on all atoms) at x; one force call unless x repeats the last one.

atoms_at
bond_between
dr_max
positions
residual
x_from

ADVANCE_BREAK_RATIO = 1.55

__init__(crack, K_G, atoms=None)

positions(x)

x_from(positions, k)

atoms_at(x)

evaluate(x): (energy, forces on all atoms) at x; one force call unless x repeats the last one.

residual(x)

bond_between(x, i, j)

dr_max(x1, x0)

property advance_break: Bond length (A) at which a crack-line bond counts as broken (settable).

property line_bonds: Heavy-atom bonds crossing the cleavage plane at or ahead of the tip, sorted along the propagation direction: (i, j, reference length).

advance(x): Crack advance in bonds: the sum over crack-line bonds of the fraction each is stretched from its reference length towards advance_break (clipped to [0, 1]).

class matscipy.fracture_mechanics.continuation.HessianSurrogate(problem, calc, ks=16.0, kb=1.5, scale=1.0)

Bases: object

Approximate Jacobian dF/dx: the Hessian of calc (a matscipy many-body calculator with a hessian property) on the heavy atoms, times a calibrated scale, plus anisotropic X-H springs (stretch ks, bend kb, eV/A^2) for H caps (h_host).

Methods

calibrate(x[, F0, n, eps, seed, clip])

Fit the scale from n finite-difference Hessian-vector products of the true forces along random unit vectors over the heavy-atom free dofs; clipped to clip.

hessian
jac

__init__(problem, calc, ks=16.0, kb=1.5, scale=1.0)

hessian(x)

jac(x)

calibrate(x, F0=None, n=2, eps=0.005, seed=0, clip=(0.3, 3.0)): Fit the scale from n finite-difference Hessian-vector products of the true forces along random unit vectors over the heavy-atom free dofs; clipped to clip.

class matscipy.fracture_mechanics.continuation.CorrectorResult(y: numpy.ndarray, F: numpy.ndarray, iters: int, calls: int, ok: bool)

Bases: object

y: ndarray

F: ndarray

iters: int

calls: int

ok: bool

__init__(y: ndarray, F: ndarray, iters: int, calls: int, ok: bool) → None

matscipy.fracture_mechanics.continuation.bordered_solve(J, t, rhs_F, rhs_g): Solve [J; t] d = [rhs_F; rhs_g] (the pseudo-arclength bordered system).

class matscipy.fracture_mechanics.continuation.BroydenCorrector(ftol=0.002, max_iter=8, step_limit=0.1)

Bases: object

Quasi-Newton on the bordered system F(y) = 0, t.(y - yp) = 0; B starts from the surrogate Jacobian and gets a good-Broyden rank-1 update from every force call.

Methods

correct
update

__init__(ftol=0.002, max_iter=8, step_limit=0.1)

update(dy, dF)

correct(F_fn, y, Fy, yp, t, dr_fn)

class matscipy.fracture_mechanics.continuation.KrylovCorrector(ftol=0.002, max_iter=8, eps_fd=0.001, eta=0.2, gmres_max=60, step_limit=0.1, max_calls=40)

Bases: object

Jacobian-free Newton-Krylov on the bordered system (GMRES; finite-difference J v), right preconditioned by an approximate Jacobian; at most max_calls force calls.

Methods

correct

__init__(ftol=0.002, max_iter=8, eps_fd=0.001, eta=0.2, gmres_max=60, step_limit=0.1, max_calls=40)

correct(F_fn, y, Fy, yp, t, dr_fn, prec_J=None)

class matscipy.fracture_mechanics.continuation.TraceSettings(ftol: float = 0.002, alpha: float = 2.0, ds0: float = 0.05, ds_min: float = 0.001, ds_max: float = 0.6, dr_cap: float = 0.05, n_target: int = 3, max_iter: int = 8, broyden_retries: int = 2, surrogate_every: int = 5, max_points: int = 400, k_window: tuple = (0.0, 3.0), stop_after_folds: int = 2, extra_after_last_fold: float = 0.25, reject_factor: float = 3.0, max_calls_per_point: int = 30, krylov_max_calls: int = 40, corrector_step_limit: float = 0.1, start_max_iter: int = 200, max_advance: float = None)

Bases: object

Attributes:

max_advance

ftol: float = 0.002

alpha: float = 2.0

ds0: float = 0.05

ds_min: float = 0.001

ds_max: float = 0.6

dr_cap: float = 0.05

n_target: int = 3

max_iter: int = 8

broyden_retries: int = 2

surrogate_every: int = 5

max_points: int = 400

k_window: tuple = (0.0, 3.0)

stop_after_folds: int = 2

extra_after_last_fold: float = 0.25

reject_factor: float = 3.0

max_calls_per_point: int = 30

krylov_max_calls: int = 40

corrector_step_limit: float = 0.1

start_max_iter: int = 200

max_advance: float = None

__init__(ftol: float = 0.002, alpha: float = 2.0, ds0: float = 0.05, ds_min: float = 0.001, ds_max: float = 0.6, dr_cap: float = 0.05, n_target: int = 3, max_iter: int = 8, broyden_retries: int = 2, surrogate_every: int = 5, max_points: int = 400, k_window: tuple = (0.0, 3.0), stop_after_folds: int = 2, extra_after_last_fold: float = 0.25, reject_factor: float = 3.0, max_calls_per_point: int = 30, krylov_max_calls: int = 40, corrector_step_limit: float = 0.1, start_max_iter: int = 200, max_advance: float = None) → None

matscipy.fracture_mechanics.continuation.trace_s_curve(problem, surrogate, k0, *, direction=1, settings=None, callback=None)

Trace the S-curve K(u) of a clamped crack cluster from load k0 = K/K_G.

The start (the problem’s geometry, moved rigidly with the field from its build load to k0) is relaxed at fixed k; then secant predictor (capped to dr_cap and, with max_advance, to that crack advance) and bordered Broyden corrector (Newton-Krylov after broyden_retries halvings), with the surrogate Jacobian refreshed every surrogate_every points. Folds are sign changes of dk/ds; the branch label starts “stable” and alternates at each fold. Stops after stop_after_folds folds plus extra_after_last_fold in k, on leaving k_window moving outward, at max_points, or when a step needs more than max_calls_per_point force calls or ds falls below ds_min (“stalled”). callback(frame) is called for every accepted point.

Returns (frames, status): extxyz-ready Atoms with arrays forces and info K_over_KG, K, energy, branch, point_index, method, ds, calls; status one of “folds”, “window”, “complete”, “stalled”.