matscipy.fracture_mechanics package

Submodules

matscipy.fracture_mechanics.clusters module

matscipy.fracture_mechanics.clusters.set_groups(a, n, skin_x, skin_y, central_x=-0.5, central_y=-0.5, invert_central=False)
matscipy.fracture_mechanics.clusters.set_regions(cryst, r_I, cutoff, r_III)
matscipy.fracture_mechanics.clusters.cluster(el, a0, n, crack_surface=[1, 1, 0], crack_front=[0, 0, 1], lattice=None, shift=None)
matscipy.fracture_mechanics.clusters.diamond(*args, **kwargs)
matscipy.fracture_mechanics.clusters.fcc(*args, **kwargs)
matscipy.fracture_mechanics.clusters.bcc(*args, **kwargs)
matscipy.fracture_mechanics.clusters.sc(*args, **kwargs)

matscipy.fracture_mechanics.crack module

class matscipy.fracture_mechanics.crack.RectilinearAnisotropicCrack

Bases: object

Near field solution for a crack in a rectilinear anisotropic elastic medium. See: G. C. Sih, P. C. Paris and G. R. Irwin, Int. J. Frac. Mech. 1, 189 (1965)

__init__()
set_plane_stress(a11, a22, a12, a16, a26, a66)
set_plane_strain(b11, b22, b33, b12, b13, b23, b16, b26, b36, b66)
displacements(r, theta, k)

Displacement field in mode I fracture. Positions are passed in cylinder coordinates.

Parameters:

  • r (array_like) – Distances from the crack tip.

  • theta (array_like) – Angles with respect to the plane of the crack.

  • k (float) – Stress intensity factor.

Returns:

  • u (array) – Displacements parallel to the plane of the crack.

  • v (array) – Displacements normal to the plane of the crack.

deformation_gradient(r, theta, k)

Deformation gradient tensor in mode I fracture. Positions are passed in cyclinder coordinates.

Parameters:

  • r (array_like) – Distances from the crack tip.

  • theta (array_like) – Angles with respect to the plane of the crack.

  • k (float) – Stress intensity factor.

Returns:

  • du_dx (array) – Derivatives of displacements parallel to the plane within the plane.

  • du_dy (array) – Derivatives of displacements parallel to the plane perpendicular to the plane.

  • dv_dx (array) – Derivatives of displacements normal to the plane of the crack within the plane.

  • dv_dy (array) – Derivatives of displacements normal to the plane of the crack perpendicular to the plane.

stresses(r, theta, k)

Stress field in mode I fracture. Positions are passed in cylinder coordinates.

Parameters:

  • r (array_like) – Distances from the crack tip.

  • theta (array_like) – Angles with respect to the plane of the crack.

  • k (float) – Stress intensity factor.

Returns:

  • sig_x (array) – Diagonal component of stress tensor parallel to the plane of the crack.

  • sig_y (array) – Diagonal component of stress tensor normal to the plane of the crack.

  • sig_xy (array) – Off-diagonal component of the stress tensor.

rtheta(u, v, k)

Invert displacement field in mode I fracture, i.e. compute r and theta from displacements.

k1g(surface_energy)

K1G, Griffith critical stress intensity in mode I fracture

k1gsqG()
matscipy.fracture_mechanics.crack.displacement_residuals(r0, crack, x, y, ref_x, ref_y, k, power=1)

Return actual displacement field minus ideal displacement field, divided by r**alpha.

matscipy.fracture_mechanics.crack.displacement_residual(r0, crack, x, y, ref_x, ref_y, k, mask=None, power=1)
matscipy.fracture_mechanics.crack.deformation_gradient_residuals(r0, crack, x, y, cur, ref_x, ref_y, ref, k, cutoff)

Return actual displacement field minus ideal displacement field.

matscipy.fracture_mechanics.crack.deformation_gradient_residual(r0, crack, x, y, cur, ref_x, ref_y, ref, k, cutoff, mask=None)
class matscipy.fracture_mechanics.crack.CubicCrystalCrack(crack_surface, crack_front, C11=None, C12=None, C44=None, stress_state='plane strain', C=None, Crot=None)

Bases: object

Crack in a cubic crystal.

__init__(crack_surface, crack_front, C11=None, C12=None, C44=None, stress_state='plane strain', C=None, Crot=None)

Initialize a crack in a cubic crystal with elastic constants C11, C12 and C44 (or optionally a full 6x6 elastic constant matrix C). The crack surface is given by crack_surface, the cracks runs in the plane given by crack_front.

k1g(surface_energy)

Compute Griffith critical stress intensity in mode I fracture.

Parameters:

surface_energy (float) – Surface energy of the respective crystal surface.

Returns:

k1g – Stress intensity factor.

Return type:

float

k1gsqG()
displacements_from_cylinder_coordinates(r, theta, k)

Displacement field in mode I fracture from cylindrical coordinates.

displacements_from_cartesian_coordinates(dx, dy, k)

Displacement field in mode I fracture from cartesian coordinates.

displacements(ref_x, ref_y, x0, y0, k)

Displacement field for a list of cartesian positions.

Parameters:

  • ref_x (array_like) – x-positions of the reference crystal.

  • ref_y (array_like) – y-positions of the reference crystal.

  • x0 (float) – x-coordinate of the crack tip.

  • y0 (float) – y-coordinate of the crack tip.

  • k (float) – Stress intensity factor.

Returns:

  • ux (array_like) – x-displacements.

  • uy (array_like) – y-displacements.

deformation_gradient_from_cylinder_coordinates(r, theta, k)

Displacement field in mode I fracture from cylindrical coordinates.

deformation_gradient_from_cartesian_coordinates(dx, dy, k)

Displacement field in mode I fracture from cartesian coordinates.

deformation_gradient(ref_x, ref_y, x0, y0, k)

Deformation gradient for a list of cartesian positions.

Parameters:

  • ref_x (array_like) – x-positions of the reference crystal.

  • ref_y (array_like) – y-positions of the reference crystal.

  • x0 (float) – x-coordinate of the crack tip.

  • y0 (float) – y-coordinate of the crack tip.

  • k (float) – Stress intensity factor.

Returns:

  • ux (array_like) – x-displacements.

  • uy (array_like) – y-displacements.

crack_tip_position(x, y, ref_x, ref_y, x0, y0, k, mask=None, residual_func=<function displacement_residual>, method='Powell', return_residuals=False)

Return an estimate of the real crack tip position by minimizing the mean square error of the current configuration relative to the displacement field obtained for a stress intensity factor k from linear elastic fracture mechanics.

Parameters:

  • x (array_like) – x-positions of the atomic system containing the crack.

  • y (array_like) – y-positions of the atomic system containing the crack.

  • ref_x (array_like) – x-positions of the reference crystal.

  • ref_y (array_like) – y-positions of the reference crystal.

  • x0 (float) – Initial guess for the x-coordinate of the crack tip.

  • y0 (float) – Initial guess for the y-coordinate of the crack tip.

  • k (float) – Stress intensity factor.

  • mask (array_like, optional) – Marks the atoms to use for this calculation.

  • residual_func (function) – Function returning per-atom residuals to be minimized.

  • method (str) – Optimization method. See method argument of scipy.optimize.minimize. Additionally, ‘leastsq’ invokes scipy.optimize.leastsq.

  • return_residuals (bool) – Function returns residuals if set to True.

Returns:

  • x0 (float) – x-coordinate of the crack tip.

  • y0 (float) – y-coordinate of the crack tip.

  • residuals (array, optional) – Per-atom residuals at end of optimization.

crack_tip_position_y(x, y, ref_x, ref_y, x0, y0, k, mask=None)

Return an estimate of the y-coordinate of the real crack tip position assuming the stress intensity factor is k.

Parameters:

  • x (array_like) – x-positions of the atomic system containing the crack.

  • y (array_like) – y-positions of the atomic system containing the crack.

  • ref_x (array_like) – x-positions of the reference crystal.

  • ref_y (array_like) – y-positions of the reference crystal.

  • x0 (float) – Initial guess for the x-coordinate of the crack tip.

  • y0 (float) – Initial guess for the y-coordinate of the crack tip.

  • k (float) – Stress intensity factor.

  • mask (array_like, optional) – Marks the atoms to use for this calculation.

Returns:

y0 – y-coordinate of the crack tip

Return type:

float

scale_displacements(x, y, ref_x, ref_y, old_k, new_k)

Rescale atomic positions from stress intensity factor old_k to the new stress intensity factor new_k. This is useful for extrapolation of relaxed positions.

Parameters:

  • x (array_like) – x-positions of the atomic system containing the crack.

  • y (array_like) – y-positions of the atomic system containing the crack.

  • ref_x (array_like) – x-positions of the reference crystal.

  • ref_y (array_like) – y-positions of the reference crystal.

  • old_k (float) – Current stress intensity factor.

  • new_k (float) – Stress intensity factor for the output positions.

Returns:

  • x (array_like) – New x-positions of the atomic system containing the crack.

  • y (array_like) – New y-positions of the atomic system containing the crack.

stresses_from_cylinder_coordinates(r, theta, k)

Stress field in mode I fracture from cylindrical coordinates

stresses_from_cartesian_coordinates(dx, dy, k)

Stress field in mode I fracture from cartesian coordinates

stresses(ref_x, ref_y, x0, y0, k)

Stress field for a list of cartesian positions

Parameters:

  • ref_x (array_like) – x-positions of the reference crystal.

  • ref_y (array_like) – y-positions of the reference crystal.

  • x0 (float) – x-coordinate of the crack tip.

  • y0 (float) – y-coordinate of the crack tip.

  • k (float) – Stress intensity factor.

Returns:

  • sig_x (array_like) – xx-component of the stress tensor

  • sig_y (array_like) – yy-component of the stress tensor

  • sig_xy (array_like) – xy-component of the stress tensor

class matscipy.fracture_mechanics.crack.SinclairCrack(crk, cryst, calc, k, alpha=0.0, vacuum=6.0, variable_alpha=True, variable_k=False, alpha_scale=None, k_scale=None, extended_far_field=False)

Bases: object

Flexible boundary conditions for a Mode I crack as described in

Sinclair, J. E. The Influence of the Interatomic Force Law and of Kinks on the

Propagation of Brittle Cracks. Philos. Mag. 31, 647–671 (1975)

__init__(crk, cryst, calc, k, alpha=0.0, vacuum=6.0, variable_alpha=True, variable_k=False, alpha_scale=None, k_scale=None, extended_far_field=False)
Parameters:

  • CubicCrystalCrack (crk - instance of) –

  • crystal (cryst - Atoms object representing undeformed) –

  • calculator (calc - ASE-compatible) –

  • factor (k - stress intensity) –

  • position (alpha - crack tip) –

  • cell (vacuum - amount of vacuum to add to unit) –

  • True (extended_far_field - if) –

  • BCs) (crack tip position can vary (flexible) –

  • True – (needed for arc-length continuation)

  • vary (stress intensity factor can) – (needed for arc-length continuation)

  • True –

  • contrib (crack tip force includes region III) –

pack(u, alpha, k)
unpack(x, reshape=False, defaults=None)
get_dofs()
set_dofs(x)
u_cle(alpha=None)

Returns CLE displacement solution at current (alpha, K)

Note that this is NOT pre-multipled by the stress intensity factor k

fit_cle(r_fit=20.0, variable_alpha=True, variable_k=True, x0=None, grid=None)
update_atoms()

Update self.atoms from degrees of freedom (self.u, self.alpha, self.k)

set_atoms(atoms)
get_crack_tip_force(forces=None, mask=None)
get_xdot(x1, x2, ds=None)
get_k_force(x1, xdot1, ds)
get_forces(x1=None, xdot1=None, ds=None, forces=None, mask=None)
update_precon(x, F=None)
get_precon(x, F)
optimize(ftol=0.001, steps=20, dump=False, args=None, precon=False, method='krylov', check_grad=True, dump_interval=10)
get_potential_energy()
rescale_k(new_k)
arc_length_continuation(x0, x1, N=10, ds=0.01, ftol=0.01, direction=1, steps=100, continuation=False, traj_file='x_traj.h5', traj_interval=1, precon=False)
plot(ax=None, regions='1234', styles=None, bonds=None, cutoff=2.8, tip=False, atoms_args=None, bonds_args=None, tip_args=None)
animate(x, k1g, regions='12', cutoff=2.8, frames=None, callback=None)
matscipy.fracture_mechanics.crack.isotropic_modeI_crack_tip_stress_field(K, r, t, xy_only=True, nu=0.5, stress_state='plane strain')

Compute Irwin singular crack tip stress field for mode I fracture.

Parameters:

  • K (float) – Mode I stress intensity factor. Units should match units of r.

  • r (array_like) – Radial distances from crack tip. Can be a multidimensional array to evaluate stress field on a grid.

  • t (array_like) – Angles from horzontal line y=0 ahead of crack tip, measured anticlockwise. Should have same shape as r.

  • xy_only (bool) – If True (default) only xx, yy, xy and yx components will be set.

  • nu (float) – Poisson ratio. Used only when xy_only=False, to determine zz stresses

  • stress_state (str) – One of”plane stress” or “plane strain”. Used if xyz_only=False to determine zz stresses.

Returns:

sigma

Return type:

array with shape r.shape + (3,3)

matscipy.fracture_mechanics.crack.isotropic_modeI_crack_tip_displacement_field(K, G, nu, r, t, stress_state='plane strain')

Compute Irwin singular crack tip displacement field for mode I fracture.

Parameters:

  • K (float) – Mode I stress intensity factor. Units should match units of G and r.

  • G (float) – Shear modulus. Units should match units of K and r.

  • nu (float) – Poisson ratio.

  • r (array_like) – Radial distances from crack tip. Can be a multidimensional array to evaluate stress field on a grid.

  • t (array_like) – Angles from horizontal line y=0 ahead of crack tip, measured anticlockwise. Should have same shape as r.

  • stress_state (str) – One of”plane stress” or “plane strain”. Used if xyz_only=False to determine zz stresses.

Returns:

  • u (array)

  • v (array) – Displacements. Same shape as r and t.

class matscipy.fracture_mechanics.crack.IsotropicStressField(K=None, x0=None, y0=None, sxx0=0.0, syy0=0.0, sxy0=0.0, nu=0.5, stress_state='plane strain')

Bases: object

Calculator to return Irwin near-tip stress field at atomic sites

__init__(K=None, x0=None, y0=None, sxx0=0.0, syy0=0.0, sxy0=0.0, nu=0.5, stress_state='plane strain')
get_stresses(atoms)
matscipy.fracture_mechanics.crack.strain_to_G(strain, E, nu, orig_height)

Convert from strain to energy release rate G for thin strip geometry

Parameters:

  • strain (float) – Dimensionless ratio (current_height - orig_height)/orig_height

  • E (float) – Young’s modulus relevant for a pull in y direction sigma_yy/eps_yy

  • nu (float) – Poission ratio -eps_yy/eps_xx

  • orig_height (float) – Unstrained height of slab

Returns:

G – Energy release rate in units consistent with input (i.e. in eV/A**2 if eV/A/fs units used)

Return type:

float

matscipy.fracture_mechanics.crack.G_to_strain(G, E, nu, orig_height)

Convert from energy release rate G to strain for thin strip geometry

Parameters:

  • G (float) – Energy release rate in units consistent with E and orig_height

  • E (float) – Young’s modulus relevant for a pull in y direction sigma_yy/eps_yy

  • nu (float) – Poission ratio -eps_yy/eps_xx

  • orig_height (float) – Unstrained height of slab

Returns:

strain – Dimensionless ratio (current_height - orig_height)/orig_height

Return type:

float

matscipy.fracture_mechanics.crack.get_strain(atoms)

Return the current strain on thin strip configuration atoms

Requires unstrained height of slab to be stored as OrigHeight key in atoms.info dictionary.

Also updates value stored in atoms.info.

matscipy.fracture_mechanics.crack.get_energy_release_rate(atoms)

Return the current energy release rate G for atoms

Result is computed assuming thin strip geometry, and using stored Young’s modulus and Poission ratio and original slab height from atoms.info dictionary.

Also updates G value stored in atoms.info dictionary.

matscipy.fracture_mechanics.crack.get_stress_intensity_factor(atoms, stress_state='plane strain')

Compute stress intensity factor K_I

Calls get_energy_release_rate() to compute G, then uses stored YoungsModulus and PoissionRatio_yz values from atoms.info dictionary to compute K_I.

Also updates value stored in atoms.info dictionary.

matscipy.fracture_mechanics.crack.fit_crack_stress_field(atoms, r_range=(0.0, 50.0), initial_params=None, fix_params=None, sigma=None, avg_sigma=None, avg_decay=0.005, calc=None, verbose=False)

Perform a least squares fit of near-tip stress field to isotropic solution

Stresses on the atoms are fit to the Irwin K-field singular crack tip solution, allowing the crack position, stress intensity factor and far-field stress components to vary during the fit.

Parameters:

  • atoms (Atoms object) –

    Crack system. For the initial fit, the following keys are used from the info dictionary:

    • YoungsModulus

    • PossionRatio_yx

    • G — current energy release rate

    • strain — current applied strain

    • CrackPos — initial guess for crack tip position

    The initial guesses for the stress intensity factor K and far-field stress sigma0 are computed from YoungsModulus, PoissonRatio_yx, G and strain, assuming plane strain in thin strip boundary conditions.

    On exit, new K, sigma0 and CrackPos entries are set in the info dictionary. These values are then used as starting guesses for subsequent fits.

  • r_range (sequence of two floats, optional) – If present, restrict the stress fit to an annular region r_range[0] <= r < r_range[1], centred on the previous crack position (from the CrackPos entry in atoms.info). If r_range is None, fit is carried out for all atoms.

  • initial_params (dict) – Names and initial values of parameters. Missing initial values are guessed from Atoms object.

  • fix_params (dict) – Names and values of parameters to fix during the fit, e.g. {y0: 0.0} to constrain the fit to the line y=0

  • sigma (None or array with shape (len(atoms), 3, 3)) – Explicitly provide the per-atom stresses. Avoids calling Atoms’ calculators get_stresses() method.

  • avg_sigma (None or array with shape (len(atoms), 3, 3)) – If present, use this array to accumulate the time-averaged stress field. Useful when processing a trajectory.

  • avg_decay (real) – Factor by which average stress is attenuated at each step. Should be set to dt/tau where dt is MD time-step and tau is a characteristic averaging time.

  • calc (Calculator object, optional) –

    If present, override the calculator used to compute stresses on the atoms. Default is atoms.get_calculator.

    To use the atom resolved stress tensor pass an instance of the AtomResolvedStressField class.

  • verbose (bool, optional) – If set to True, print additional information about the fit.

Returns:

params – Fitted parameters, in a form suitable for passing to IsotropicStressField constructor. These are the stress intensity factor K, the centre of the stress field (x0, y0), and the far field contribution to the stress (sxx0, syy0, sxy0).

Return type:

dict with keys [K, x0, y0, sxx0, syy0, sxy0]

matscipy.fracture_mechanics.crack.find_tip_coordination(a, bondlength=2.6, bulk_nn=4)

Find position of tip in crack cluster from coordination

matscipy.fracture_mechanics.crack.find_tip_broken_bonds(atoms, cutoff, bulk_nn=4, boundary_thickness=None)

Find position of the tip from the atom coordination, i.e. broken bonds. Using the C implementation of ‘neighbour_list’. Returns the tip’s position in cartesian coordinates.

Parameters:

  • atoms (ase.Atoms) – Atomic configuration.

  • cutoff (float) – Cutoff distance for neighbour search.

  • bulk_nn (integer) – Number of nearest neighbours for the standard bulk configuration.

  • boundary_buffer (float) – Thickness of the boundaries. Defaults to cutoff distance.

Returns:

tip_position – The x and y values are found. The z value is calculated as the midpoint of the depth.

Return type:

numpy array

matscipy.fracture_mechanics.crack.find_tip_stress_field(atoms, r_range=None, initial_params=None, fix_params=None, sigma=None, avg_sigma=None, avg_decay=0.005, calc=None)

Find the position of crack tip by fitting to the isotropic K-field stress

Fit is carried out using fit_crack_stress_field(), and parameters have the same meaning as there.

See also

fit_crack_stress_field

matscipy.fracture_mechanics.crack.plot_stress_fields(atoms, r_range=None, initial_params=None, fix_params=None, sigma=None, avg_sigma=None, avg_decay=0.005, calc=None)

Fit and plot atomistic and continuum stress fields

Firstly a fit to the Irwin K-field solution is carried out using fit_crack_stress_field(), and parameters have the same meaning as for that function. Then plots of the \(\sigma_{xx}\), \(\sigma_{yy}\), \(\sigma_{xy}\) fields are produced for atomistic and continuum cases, and for the residual error after fitting.

matscipy.fracture_mechanics.crack.thin_strip_displacement_y(x, y, strain, a, b)

Return vertical displacement ramp used to apply initial strain to slab

Strain is increased from 0 to strain over distance \(a <= x <= b\). Region \(x < a\) is rigidly shifted up/down by strain*height/2.

Here is an example of how to use this function on an artificial 2D square atomic lattice. The positions are plotted before (left) and after (right) applying the displacement, and the horizontal and vertical lines show the strain (red), a (green) and b (blue) parameters.

import matplotlib.pyplot as plt
import numpy as np

w = 1; h = 1; strain = 0.1; a = -0.5; b =  0.0
x = np.linspace(-w, w, 20)
y = np.linspace(-h, h, 20)
X, Y = np.meshgrid(x, y)
u_y = thin_strip_displacement_y(X, Y, strain, a, b)

for i, disp in enumerate([0, u_y]):
    plt.subplot(1,2,i+1)
    plt.scatter(X, Y + disp, c='k', s=5)
    for y in [-h, h]:
        plt.axhline(y, color='r', linewidth=2, linestyle='dashed')
        plt.axhline(y*(1+strain), color='r', linewidth=2)
    for x, c in zip([a, b], ['g', 'b']):
        plt.axvline(x, color=c, linewidth=2)
thin-strip-displacement-y.png
Parameters:

  • x (array) –

  • y (array) – Atomic positions in unstrained slab, centered on origin x=0,y=0

  • strain (float) – Far field strain to apply

  • a (float) – x coordinate for beginning of strain ramp

  • b (float) – x coordinate for end of strain ramp

matscipy.fracture_mechanics.crack.print_crack_system(directions)

Pretty printing of crack crystallographic coordinate system

Specified by list of Miller indices for crack_direction (x), cleavage_plane (y) and crack_front (z), each of which should be a sequence of three floats

class matscipy.fracture_mechanics.crack.ConstantStrainRate(orig_height, delta_strain, mask=None)

Bases: object

Constraint which increments epsilon_yy at a constant strain rate

Rescaling is applied only to atoms where mask is True (default is all atoms)

__init__(orig_height, delta_strain, mask=None)
adjust_forces(atoms, forces)
adjust_positions(atoms, newpos)
copy()
apply_strain(atoms, rigid_constraints=False)

Applies a constant strain to the system.

Parameters:

  • atoms (ASE.atoms) – Atomic configuration.

  • rigid_constraints (boolean) – Apply (or not apply) strain to every atom. i.e. allow constrainted atoms to move during strain application

matscipy.fracture_mechanics.crackpathsel module

function coordination(r, cutoff, transition_width)
if r > cutoff

f = 0.0 df = 0.0

elseif r > cutoff-transition_width

f = 0.5 * ( cos(pi*(r-cutoff+transition_width)/transition_width) + 1.0 ) df = - 0.5 * pi * sin(pi*(r-cutoff+transition_width)/transition_width) / transition_width

else

f = 1.0 df = 0.0

end f

end

function dcoordination(r, cutoff, transition_width)
if r > cutoff

df = 0.0

elseif r > cutoff-transition_width

df = - 0.5 * pi * sin(pi*(r-cutoff+transition_width)/transition_width) / transition_width

else

df = 0.0

end df

end

function energy(pos, neighb_j, neighb_rij, cutoff, transition_width, epsilon)

N = size(pos, 2)

n = zeros(Float64, N) energies = zeros(Float64, N)

for i = 1:N
for (m, j) in enumerate(neighb_j[i])

r_ij = neighb_rij[i][m] #@printf(“i %d j %d r_ij %f

“, i, j, r_ij)

r_ij > cutoff && continue

f_ij = coordination(r_ij, cutoff, transition_width) n[i] += f_ij

end energies[i] += (n[i] - 3)^2

end

for i = 1:N

sum_B_ij = 0.0

for (m, j) in enumerate(neighb_j[i])

r_ij = neighb_rij[i][m] r_ij > cutoff && continue

f_ij = coordination(r_ij, cutoff, transition_width) B_ij = (n[j] - 3.0)^2*f_ij sum_B_ij += B_ij

end

Eb_i = epsilon*(n[i] - 3.0)^2*sum_B_ij energies[i] += Eb_i

end

E = sum(energies) return (E, energies, n)

end

function force(pos, neighb_j, neighb_rij, cutoff, transition_width, epsilon, dx)

N = size(pos, 2) f = zeros(Float64, (3, N)) p = zeros(Float64, (3, N))

p[:, :] = pos for i = 1:N

for j = 1:3

p[j, i] += dx ep, local_e_p, n_p = energy(p, neighb_j, neighb_rij, cutoff, transition_width, epsilon) p[j, i] -= 2dx em, local_e_m, n_m = energy(p, neighb_j, neighb_rij, cutoff, transition_width, epsilon) f[j, i] = -(ep - em)/(2dx) p[j, i] += dx

end

end f

end

matscipy.fracture_mechanics.energy_release module

matscipy.fracture_mechanics.energy_release.J_integral(a, deformation_gradient, virial, epot, e0, tip_x, tip_y, r1, r2, mask=None)

Compute the energy release rate from the J-integral. Converts contour integral into a domain integral. See: Li, Shih, Needleman, Eng. Fract. Mech. 21, 405 (1985); Jin, Yuan, J. Nanosci. Nanotech. 5, 2099 (2005) Domain function is currently fixed: q(r) = (r-r1)/(r2-r1)

Parameters:

  • a (ase.Atoms) – Relaxed atomic configuration of the crack.

  • deformation_gradient (array_like) – len(a) x 3x3 array of atomic deformation gradients.

  • virial (array_like) – len(a) x 3x3 array of atomic virials.

  • epot (float) – Potential energy per atom of the cracked configuration.

  • e0 (float) – Reference energy (cohesive energy per atom of the bulk crystal at equilibrium).

  • tip_x (float) – Position of the crack tip.

  • tip_y (float) – Position of the crack tip.

  • r1 (float) – Volume integration is carried out in region at a distance between r1 and r2 from the crack tip.

  • r2 (float) – Volume integration is carried out in region at a distance between r1 and r2 from the crack tip.

  • mask (array_like) – Include only a subset of all atoms into J-integral computation.

Returns:

J – Value of the J-integral.

Return type:

float

matscipy.fracture_mechanics.idealbrittlesolid module

matscipy.fracture_mechanics.idealbrittlesolid.triangular_lattice_slab(a, n, m)
matscipy.fracture_mechanics.idealbrittlesolid.find_triangles_2d(atoms, cutoff, minangle=0.5235987755982988, maxangle=2.0943951023931953, xdim=0, ydim=1)

Return a list of all triangles of a triangular lattice sitting in the x-y plane.

class matscipy.fracture_mechanics.idealbrittlesolid.IdealBrittleSolid(*args, **kwargs)

Bases: Calculator

Implementation of force field for an ideal brittle solid

Described in Marder, Int. J. Fract. 130, 517-555 (2004)

implemented_properties: List[str] = ['energy', 'free_energy', 'energies', 'stress', 'forces']

Properties calculator can handle (energy, forces, …)

default_parameters: Dict[str, Any] = {'a': 1.0, 'b': 0.01, 'beta': 0.01, 'k': 1.0, 'linear': False, 'rc': 1.01}

Default parameters

__init__(*args, **kwargs)

Basic calculator implementation.

restart: str

Prefix for restart file. May contain a directory. Default is None: don’t restart.

ignore_bad_restart_file: bool

Deprecated, please do not use. Passing more than one positional argument to Calculator() is deprecated and will stop working in the future. Ignore broken or missing restart file. By default, it is an error if the restart file is missing or broken.

directory: str or PurePath

Working directory in which to read and write files and perform calculations.

label: str

Name used for all files. Not supported by all calculators. May contain a directory, but please use the directory parameter for that instead.

atoms: Atoms object

Optional Atoms object to which the calculator will be attached. When restarting, atoms will get its positions and unit-cell updated from file.

set_reference_crystal(crystal)
calculate(atoms, properties, system_changes)

Do the calculation.

properties: list of str

List of what needs to be calculated. Can be any combination of ‘energy’, ‘forces’, ‘stress’, ‘dipole’, ‘charges’, ‘magmom’ and ‘magmoms’.

system_changes: list of str

List of what has changed since last calculation. Can be any combination of these six: ‘positions’, ‘numbers’, ‘cell’, ‘pbc’, ‘initial_charges’ and ‘initial_magmoms’.

Subclasses need to implement this, but can ignore properties and system_changes if they want. Calculated properties should be inserted into results dictionary like shown in this dummy example:

self.results = {'energy': 0.0,
                'forces': np.zeros((len(atoms), 3)),
                'stress': np.zeros(6),
                'dipole': np.zeros(3),
                'charges': np.zeros(len(atoms)),
                'magmom': 0.0,
                'magmoms': np.zeros(len(atoms))}

The subclass implementation should first call this implementation to set the atoms attribute and create any missing directories.

get_wave_speeds(atoms)

Return longitudinal, shear and Rayleigh wave speeds

get_elastic_moduli()

Return Lam’e constants lambda and mu

get_youngs_modulus()
get_poisson_ratio()
band_structure()

Create band-structure object for plotting.

calculate_numerical_forces(atoms, d=0.001)

Calculate numerical forces using finite difference.

All atoms will be displaced by +d and -d in all directions.

calculate_numerical_stress(atoms, d=1e-06, voigt=True)

Calculate numerical stress using finite difference.

calculate_properties(atoms, properties)

This method is experimental; currently for internal use.

calculation_required(atoms, properties)
check_state(atoms, tol=1e-15)

Check for any system changes since last calculation.

property directory: str
discard_results_on_any_change = False

Whether we purge the results following any change in the set() method.

export_properties()
get_atoms()
get_charges(atoms=None)
get_default_parameters()
get_dipole_moment(atoms=None)
get_forces(atoms=None)
get_magnetic_moment(atoms=None)
get_magnetic_moments(atoms=None)

Calculate magnetic moments projected onto atoms.

get_potential_energies(atoms=None)
get_potential_energy(atoms=None, force_consistent=False)
get_property(name, atoms=None, allow_calculation=True)

Get the named property.

get_stress(atoms=None)
get_stresses(atoms=None)

the calculator should return intensive stresses, i.e., such that stresses.sum(axis=0) == stress

ignored_changes: Set[str] = {}

Properties of Atoms which we ignore for the purposes of cache

property label
read(label)

Read atoms, parameters and calculated properties from output file.

Read result from self.label file. Raise ReadError if the file is not there. If the file is corrupted or contains an error message from the calculation, a ReadError should also be raised. In case of succes, these attributes must set:

atoms: Atoms object

The state of the atoms from last calculation.

parameters: Parameters object

The parameter dictionary.

results: dict

Calculated properties like energy and forces.

The FileIOCalculator.read() method will typically read atoms and parameters and get the results dict by calling the read_results() method.

classmethod read_atoms(restart, **kwargs)
reset()

Clear all information from old calculation.

set(**kwargs)

Set parameters like set(key1=value1, key2=value2, …).

A dictionary containing the parameters that have been changed is returned.

Subclasses must implement a set() method that will look at the chaneged parameters and decide if a call to reset() is needed. If the changed parameters are harmless, like a change in verbosity, then there is no need to call reset().

The special keyword ‘parameters’ can be used to read parameters from a file.

set_label(label)

Set label and convert label to directory and prefix.

Examples:

  • label=’abc’: (directory=’.’, prefix=’abc’)

  • label=’dir1/abc’: (directory=’dir1’, prefix=’abc’)

  • label=None: (directory=’.’, prefix=None)

todict(skip_default=True)
matscipy.fracture_mechanics.idealbrittlesolid.find_crack_tip(atoms, dt=None, store=True, results=None)

Return atom at the crack tip and its x-coordinate

Crack tip is defined to be location of rightmost atom whose nearest neighbour is at distance > 2.5*a

matscipy.fracture_mechanics.idealbrittlesolid.set_initial_velocities(c)

Initialise a dynamical state by kicking some atoms behind tip

matscipy.fracture_mechanics.idealbrittlesolid.set_constraints(c, a)
matscipy.fracture_mechanics.idealbrittlesolid.extend_strip(atoms, a, N, M, vacuum)

Module contents