Skip to content

Commit

Permalink
Merge branch 'adaptivecont' of github.com:libAtoms/matscipy into adap…
Browse files Browse the repository at this point in the history 
…tivecont
  • Loading branch information
@Fraser-Birks
Fraser-Birks committed Jul 13, 2023
2 parents e7c57f6 + 46609ad commit 441a027
Show file tree
Hide file tree
Showing 5 changed files with 254 additions and 251 deletions.
3 changes: 2 additions & 1 deletion examples/fracture_mechanics/sinclair_crack/params_Fe.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
53 changes: 0 additions & 53 deletions matscipy/fracture_mechanics/clusters.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.')
198 changes: 1 addition & 197 deletions matscipy/fracture_mechanics/crack.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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):
Expand Down
Loading

0 comments on commit 441a027

Please sign in to comment.