From ddc7a13127f298dc091328cd7ace1aa09cdcd3e7 Mon Sep 17 00:00:00 2001 From: lakshenoy Date: Thu, 13 Jul 2023 10:47:23 +0100 Subject: [PATCH] rearrange functions --- .../sinclair_crack/params_Fe.py | 3 +- matscipy/fracture_mechanics/clusters.py | 53 ----- matscipy/fracture_mechanics/crack.py | 198 +----------------- matscipy/fracture_mechanics/solvers.py | 197 +++++++++++++++++ matscipy/surface.py | 54 +++++ 5 files changed, 254 insertions(+), 251 deletions(-) create mode 100644 matscipy/fracture_mechanics/solvers.py diff --git a/examples/fracture_mechanics/sinclair_crack/params_Fe.py b/examples/fracture_mechanics/sinclair_crack/params_Fe.py index 9d7ca68e..184670b1 100644 --- a/examples/fracture_mechanics/sinclair_crack/params_Fe.py +++ b/examples/fracture_mechanics/sinclair_crack/params_Fe.py @@ -25,7 +25,8 @@ crack_front = [0, 0, 1] # crack surface energy -from matscipy.fracture_mechanics.clusters import find_surface_energy +#from matscipy.fracture_mechanics.clusters import find_surface_energy +from matscipy.surface import find_surface_energy surface_energy = find_surface_energy(el,calc,a0,'bcc100',unit='0.1J/m^2')[0] # simulation control diff --git a/matscipy/fracture_mechanics/clusters.py b/matscipy/fracture_mechanics/clusters.py index 5af1ad40..1b8d5c78 100644 --- a/matscipy/fracture_mechanics/clusters.py +++ b/matscipy/fracture_mechanics/clusters.py @@ -171,56 +171,3 @@ def sc(*args, **kwargs): kwargs['lattice'] = SimpleCubic return cluster(*args, **kwargs) - -def find_surface_energy(symbol,calc,a0,surface,size=(8,1,1),vacuum=10,fmax=0.0001,unit='0.1J/m^2'): - - # Import required lattice builder - if surface.startswith('bcc'): - from ase.lattice.cubic import BodyCenteredCubic as lattice_builder - elif surface.startswith('fcc'): - from ase.lattice.cubic import FaceCenteredCubic as lattice_builder #untested - elif surface.startswith('diamond'): - from ase.lattice.cubic import Diamond as lattice_builder #untested - ## Append other lattice builders here - else: - print('Error: Unsupported lattice ordering.') - - # Set orthogonal directions for cell axes - if surface.endswith('100'): - directions=[[1,0,0], [0,1,0], [0,0,1]] #tested for bcc - elif surface.endswith('110'): - directions=[[1,1,0], [-1,1,0], [0,0,1]] #tested for bcc - elif surface.endswith('111'): - directions=[[1,1,1], [-2,1,1],[0,-1,1]] #tested for bcc - ## Append other cell axis options here - else: - print('Error: Unsupported surface orientation.') - - # Make bulk and slab with same number of atoms (size) - bulk = lattice_builder(directions=directions, size=size, symbol=symbol, latticeconstant=a0, pbc=(1,1,1)) - cell = bulk.get_cell() ; cell[0,:] *=2 # vacuum along x axis (surface normal) - slab = bulk.copy() ; slab.set_cell(cell) - - # Optimize the geometries - from ase.optimize import LBFGSLineSearch - bulk.calc = calc ; opt_bulk = LBFGSLineSearch(bulk) ; opt_bulk.run(fmax=fmax) - slab.calc = calc ; opt_slab = LBFGSLineSearch(slab) ; opt_slab.run(fmax=fmax) - - # Find surface energy - import numpy as np - Ebulk = bulk.get_potential_energy() ; Eslab = slab.get_potential_energy() - area = np.linalg.norm(np.cross(slab.get_cell()[1,:],slab.get_cell()[2,:])) - gamma_ase = (Eslab - Ebulk)/(2*area) - - # Convert to required units - if unit == 'ASE': - return [gamma_ase,'ase_units'] - else: - from ase import units - gamma_SI = (gamma_ase / units.J ) * (units.m)**2 - if unit =='J/m^2': - return [gamma_SI,'J/m^2'] - elif unit == '0.1J/m^2': - return [10*gamma_SI,'0.1J/m^2'] # units required for the fracture code - else: - print('Error: Unsupported unit of surface energy.') \ No newline at end of file diff --git a/matscipy/fracture_mechanics/crack.py b/matscipy/fracture_mechanics/crack.py index 56106f35..6ea23ee9 100644 --- a/matscipy/fracture_mechanics/crack.py +++ b/matscipy/fracture_mechanics/crack.py @@ -43,8 +43,8 @@ rotate_cubic_elastic_constants, Voigt_6_to_full_3x3_stress) from matscipy.surface import MillerDirection, MillerPlane - from matscipy.neighbours import neighbour_list +from matscipy.fracture_mechanics.solvers import ode12r import matplotlib.pyplot as plt from ase.optimize.precon import Exp, PreconLBFGS @@ -61,202 +61,6 @@ MPa_sqrt_m = 1e6*units.Pascal*np.sqrt(units.m) -### -# Temporary function. Eventually meant to be pushed to ase.optimize.ode - -import numpy as np -from ase.optimize.sciopt import SciPyOptimizer, OptimizerConvergenceError -def ode12r(f, X0, args=None, h=None, verbose=1, fmax=1e-6, maxtol=1e3, steps=100, - rtol=1e-1, C1=1e-2, C2=2.0, hmin=1e-10, extrapolate=3, - callback=None, apply_precon=None, converged=None, residual=None): - """ - Adaptive ODE solver, which uses 1st and 2nd order approximations to - estimate local error and choose a new step length. - - This optimizer is described in detail in: - - S. Makri, C. Ortner and J. R. Kermode, J. Chem. Phys. - 150, 094109 (2019) - https://dx.doi.org/10.1063/1.5064465 - - Parameters - ---------- - - f : function - function returning driving force on system - X0 : 1-dimensional array - initial value of degrees of freedom - h : float - step size, if None an estimate is used based on ODE12 - verbose: int - verbosity level. 0=no output, 1=log output (default), 2=debug output. - fmax : float - convergence tolerance for residual force - maxtol: float - terminate if reisdual exceeds this value - rtol : float - relative tolerance - C1 : float - sufficient contraction parameter - C2 : float - residual growth control (Inf means there is no control) - hmin : float - minimal allowed step size - extrapolate : int - extrapolation style (3 seems the most robust) - callback : function - optional callback function to call after each update step - apply_precon: function - apply a apply_preconditioner to the optimisation - converged: function - alternative function to check convergence, rather than - using a simple norm of the forces. - residual: function - compute the residual from the current forces - - Returns - ------- - - X: array - final value of degrees of freedom - """ - - X = X0 - if args is None: - Fn = f(X) - else: - #x1, xdot1, ds = args - #Fn = f(X, x1, xdot1, ds) - Fn = f(X, *args) - # For debugging - print('Last 4 forces:', Fn[-4:]) - - if callback is None: - def callback(X): - pass - callback(X) - - if residual is None: - def residual(F, X): - return np.linalg.norm(F, np.inf) - Rn = residual(Fn, X) - - if apply_precon is None: - def apply_precon(F, X): - return F, residual(F, X) - Fp, Rp = apply_precon(Fn, X) - # For debugging - if args is not None: - print('Last 4 precon forces:',Fp[-4:]) - - def log(*args): - if verbose >= 1: - print(*args) - - def debug(*args): - if verbose >= 2: - print(*args) - - if converged is None: - def converged(F, X): - return residual(F, X) <= fmax - - if converged(Fn, X): - log("ODE12r terminates successfully after 0 iterations") - return X, 0 - if Rn >= maxtol: - raise OptimizerConvergenceError(f"ODE12r: Residual {Rn} is too large " - "at iteration 0") - - # computation of the initial step - r = residual(Fp, X) # pick the biggest force - if h is None: - h = 0.5 * rtol ** 0.5 / r # Chose a stepsize based on that force - h = max(h, hmin) # Make sure the step size is not too big - - for nit in range(1, steps): - Xnew = X + h * Fp # Pick a new position - #Fn_new = f(Xnew) if args is None else f(Xnew, x1, xdot1, ds) # Calculate the new forces at this position - Fn_new = f(Xnew) if args is None else f(Xnew, *args) # Calculate the new forces at this position - # For debugging - if args is not None: - print('Last 4 forces:',Fn_new[-4:]) - Rn_new = residual(Fn_new, Xnew) - Fp_new, Rp_new = apply_precon(Fn_new, Xnew) - # For debugging - if args is not None: - print('Last 4 precon forces:',Fp_new[-4:]) - - e = 0.5 * h * (Fp_new - Fp) # Estimate the area under the forces curve - err = np.linalg.norm(e, np.inf) # Error estimate - - # Accept step if residual decreases sufficiently and/or error acceptable - accept = ((Rp_new <= Rp * (1 - C1 * h)) or - ((Rp_new <= Rp * C2) and err <= rtol)) - - # Pick an extrapolation scheme for the system & find new increment - y = Fp - Fp_new - if extrapolate == 1: # F(xn + h Fp) - h_ls = h * (Fp @ y) / (y @ y) - elif extrapolate == 2: # F(Xn + h Fp) - h_ls = h * (Fp @ Fp_new) / (Fp @ y + 1e-10) - elif extrapolate == 3: # min | F(Xn + h Fp) | - h_ls = h * (Fp @ y) / (y @ y + 1e-10) - else: - raise ValueError(f'invalid extrapolate value: {extrapolate}. ' - 'Must be 1, 2 or 3') - if np.isnan(h_ls) or h_ls < hmin: # Rejects if increment is too small - h_ls = np.inf - - h_err = h * 0.5 * np.sqrt(rtol / err) - - # Accept the step and do the update - if accept: - X = Xnew - Rn = Rn_new - Fn = Fn_new - Fp = Fp_new - Rp = Rp_new - callback(X) - - # We check the residuals again - if Rn >= maxtol: - raise OptimizerConvergenceError( - f"ODE12r: Residual {Rn} is too " - f"large at iteration number {nit}") - - if converged(Fn, X): - log(f"ODE12r: terminates successfully " - f"after {nit} iterations.") - return X, nit - - # Compute a new step size. - # Based on the extrapolation and some other heuristics - h = max(0.25 * h, - min(4 * h, h_err, h_ls)) # Log steep-size analytic results - - debug(f"ODE12r: accept: new h = {h}, |F| = {Rp}") - debug(f"ODE12r: hls = {h_ls}") - debug(f"ODE12r: herr = {h_err}") - else: - # Compute a new step size. - h = max(0.1 * h, min(0.25 * h, h_err, - h_ls)) - debug(f"ODE12r: reject: new h = {h}") - debug(f"ODE12r: |Fnew| = {Rp_new}") - debug(f"ODE12r: |Fold| = {Rp}") - debug(f"ODE12r: |Fnew|/|Fold| = {Rp_new/Rp}") - - # abort if step size is too small - if abs(h) <= hmin: - raise OptimizerConvergenceError('ODE12r terminates unsuccessfully' - f' Step size {h} too small') - - raise OptimizerConvergenceError(f'ODE12r terminates unsuccessfully after ' - f'{steps} iterations.') - - -### #wrapper to count function calls for writing def counted(f): def wrapped(*args): diff --git a/matscipy/fracture_mechanics/solvers.py b/matscipy/fracture_mechanics/solvers.py new file mode 100644 index 00000000..60a4b382 --- /dev/null +++ b/matscipy/fracture_mechanics/solvers.py @@ -0,0 +1,197 @@ + +# This function is temporarily here, for use in NCFlex. +# Eventually meant to be pushed to ase.optimize.ode + +import numpy as np +from ase.optimize.sciopt import SciPyOptimizer, OptimizerConvergenceError + +def ode12r(f, X0, args=None, h=None, verbose=1, fmax=1e-6, maxtol=1e3, steps=100, + rtol=1e-1, C1=1e-2, C2=2.0, hmin=1e-10, extrapolate=3, + callback=None, apply_precon=None, converged=None, residual=None): + """ + Adaptive ODE solver, which uses 1st and 2nd order approximations to + estimate local error and choose a new step length. The original + implementation in ase.optimize.ode, and has been modified to this + matscipy version to ensure compatibility with the NCFlex routines. + + This optimizer is described in detail in: + + S. Makri, C. Ortner and J. R. Kermode, J. Chem. Phys. + 150, 094109 (2019) + https://dx.doi.org/10.1063/1.5064465 + + Parameters + ---------- + + f : function + function returning driving force on system + X0 : 1-dimensional array + initial value of degrees of freedom + h : float + step size, if None an estimate is used based on ODE12 + verbose: int + verbosity level. 0=no output, 1=log output (default), 2=debug output. + fmax : float + convergence tolerance for residual force + maxtol: float + terminate if reisdual exceeds this value + rtol : float + relative tolerance + C1 : float + sufficient contraction parameter + C2 : float + residual growth control (Inf means there is no control) + hmin : float + minimal allowed step size + extrapolate : int + extrapolation style (3 seems the most robust) + callback : function + optional callback function to call after each update step + apply_precon: function + apply a apply_preconditioner to the optimisation + converged: function + alternative function to check convergence, rather than + using a simple norm of the forces. + residual: function + compute the residual from the current forces + + Returns + ------- + + X: array + final value of degrees of freedom + """ + + X = X0 + if args is None: + Fn = f(X) + else: + #x1, xdot1, ds = args + #Fn = f(X, x1, xdot1, ds) + Fn = f(X, *args) + # For debugging + print('Last 4 forces:', Fn[-4:]) + + if callback is None: + def callback(X): + pass + callback(X) + + if residual is None: + def residual(F, X): + return np.linalg.norm(F, np.inf) + Rn = residual(Fn, X) + + if apply_precon is None: + def apply_precon(F, X): + return F, residual(F, X) + Fp, Rp = apply_precon(Fn, X) + # For debugging + if args is not None: + print('Last 4 precon forces:',Fp[-4:]) + + def log(*args): + if verbose >= 1: + print(*args) + + def debug(*args): + if verbose >= 2: + print(*args) + + if converged is None: + def converged(F, X): + return residual(F, X) <= fmax + + if converged(Fn, X): + log("ODE12r terminates successfully after 0 iterations") + return X, 0 + if Rn >= maxtol: + raise OptimizerConvergenceError(f"ODE12r: Residual {Rn} is too large " + "at iteration 0") + + # computation of the initial step + r = residual(Fp, X) # pick the biggest force + if h is None: + h = 0.5 * rtol ** 0.5 / r # Chose a stepsize based on that force + h = max(h, hmin) # Make sure the step size is not too big + + for nit in range(1, steps): + Xnew = X + h * Fp # Pick a new position + #Fn_new = f(Xnew) if args is None else f(Xnew, x1, xdot1, ds) # Calculate the new forces at this position + Fn_new = f(Xnew) if args is None else f(Xnew, *args) # Calculate the new forces at this position + # For debugging + if args is not None: + print('Last 4 forces:',Fn_new[-4:]) + Rn_new = residual(Fn_new, Xnew) + Fp_new, Rp_new = apply_precon(Fn_new, Xnew) + # For debugging + if args is not None: + print('Last 4 precon forces:',Fp_new[-4:]) + + e = 0.5 * h * (Fp_new - Fp) # Estimate the area under the forces curve + err = np.linalg.norm(e, np.inf) # Error estimate + + # Accept step if residual decreases sufficiently and/or error acceptable + accept = ((Rp_new <= Rp * (1 - C1 * h)) or + ((Rp_new <= Rp * C2) and err <= rtol)) + + # Pick an extrapolation scheme for the system & find new increment + y = Fp - Fp_new + if extrapolate == 1: # F(xn + h Fp) + h_ls = h * (Fp @ y) / (y @ y) + elif extrapolate == 2: # F(Xn + h Fp) + h_ls = h * (Fp @ Fp_new) / (Fp @ y + 1e-10) + elif extrapolate == 3: # min | F(Xn + h Fp) | + h_ls = h * (Fp @ y) / (y @ y + 1e-10) + else: + raise ValueError(f'invalid extrapolate value: {extrapolate}. ' + 'Must be 1, 2 or 3') + if np.isnan(h_ls) or h_ls < hmin: # Rejects if increment is too small + h_ls = np.inf + + h_err = h * 0.5 * np.sqrt(rtol / err) + + # Accept the step and do the update + if accept: + X = Xnew + Rn = Rn_new + Fn = Fn_new + Fp = Fp_new + Rp = Rp_new + callback(X) + + # We check the residuals again + if Rn >= maxtol: + raise OptimizerConvergenceError( + f"ODE12r: Residual {Rn} is too " + f"large at iteration number {nit}") + + if converged(Fn, X): + log(f"ODE12r: terminates successfully " + f"after {nit} iterations.") + return X, nit + + # Compute a new step size. + # Based on the extrapolation and some other heuristics + h = max(0.25 * h, + min(4 * h, h_err, h_ls)) # Log steep-size analytic results + + debug(f"ODE12r: accept: new h = {h}, |F| = {Rp}") + debug(f"ODE12r: hls = {h_ls}") + debug(f"ODE12r: herr = {h_err}") + else: + # Compute a new step size. + h = max(0.1 * h, min(0.25 * h, h_err, + h_ls)) + debug(f"ODE12r: reject: new h = {h}") + debug(f"ODE12r: |Fnew| = {Rp_new}") + debug(f"ODE12r: |Fold| = {Rp}") + debug(f"ODE12r: |Fnew|/|Fold| = {Rp_new/Rp}") + + # abort if step size is too small + if abs(h) <= hmin: + raise OptimizerConvergenceError('ODE12r terminates unsuccessfully' + f' Step size {h} too small') + + raise OptimizerConvergenceError(f'ODE12r terminates unsuccessfully after ' + f'{steps} iterations.') diff --git a/matscipy/surface.py b/matscipy/surface.py index 2dd85be6..4fbbd271 100644 --- a/matscipy/surface.py +++ b/matscipy/surface.py @@ -281,3 +281,57 @@ def make_unit_slab(unit_cell, axes): sup.set_scaled_positions(sup.get_scaled_positions()) return sup + + +def find_surface_energy(symbol,calc,a0,surface,size=(8,1,1),vacuum=10,fmax=0.0001,unit='0.1J/m^2'): + + # Import required lattice builder + if surface.startswith('bcc'): + from ase.lattice.cubic import BodyCenteredCubic as lattice_builder + elif surface.startswith('fcc'): + from ase.lattice.cubic import FaceCenteredCubic as lattice_builder #untested + elif surface.startswith('diamond'): + from ase.lattice.cubic import Diamond as lattice_builder #untested + ## Append other lattice builders here + else: + print('Error: Unsupported lattice ordering.') + + # Set orthogonal directions for cell axes + if surface.endswith('100'): + directions=[[1,0,0], [0,1,0], [0,0,1]] #tested for bcc + elif surface.endswith('110'): + directions=[[1,1,0], [-1,1,0], [0,0,1]] #tested for bcc + elif surface.endswith('111'): + directions=[[1,1,1], [-2,1,1],[0,-1,1]] #tested for bcc + ## Append other cell axis options here + else: + print('Error: Unsupported surface orientation.') + + # Make bulk and slab with same number of atoms (size) + bulk = lattice_builder(directions=directions, size=size, symbol=symbol, latticeconstant=a0, pbc=(1,1,1)) + cell = bulk.get_cell() ; cell[0,:] *=2 # vacuum along x axis (surface normal) + slab = bulk.copy() ; slab.set_cell(cell) + + # Optimize the geometries + from ase.optimize import LBFGSLineSearch + bulk.calc = calc ; opt_bulk = LBFGSLineSearch(bulk) ; opt_bulk.run(fmax=fmax) + slab.calc = calc ; opt_slab = LBFGSLineSearch(slab) ; opt_slab.run(fmax=fmax) + + # Find surface energy + import numpy as np + Ebulk = bulk.get_potential_energy() ; Eslab = slab.get_potential_energy() + area = np.linalg.norm(np.cross(slab.get_cell()[1,:],slab.get_cell()[2,:])) + gamma_ase = (Eslab - Ebulk)/(2*area) + + # Convert to required units + if unit == 'ASE': + return [gamma_ase,'ase_units'] + else: + from ase import units + gamma_SI = (gamma_ase / units.J ) * (units.m)**2 + if unit =='J/m^2': + return [gamma_SI,'J/m^2'] + elif unit == '0.1J/m^2': + return [10*gamma_SI,'0.1J/m^2'] # units required for the fracture code + else: + print('Error: Unsupported unit of surface energy.')