MAINT: Refactored fracture mechanics scripts to CLI

docs/cli/interatomic.md → docs/cli/calculators.md

@@ -4,7 +4,7 @@
 
 This command generates an averaged EAM potential for an alloy. The basic theory behind averaged alloy EAM potentials is
 described in
-[Varvenne, Luque,  Nöhring, Curtin, Phys. Rev. B 93, 104201 (2016)](https://doi.org/10.1103/PhysRevB.93.104201). The
+[Varvenne, Luque, Nöhring, Curtin, Phys. Rev. B 93, 104201 (2016)](https://doi.org/10.1103/PhysRevB.93.104201). The
 basic usage of the command is
 
 ```bash

docs/cli/fracture_mechanics.md

@@ -0,0 +1,67 @@
+# Fracture mechanics
+
+## `matscipy-quasistatic-crack`
+
+This command carries out a quasistatic K-field-controlled crack calculation on a small cluster. The outer boundaries of
+the cluster are fixed according to the near field solution of linear elastic fracture mechanics, considering elastic
+anisotropy, as described in [Sih, Paris, Irwin, Int. J. Fract. Mech. 1, 189 (1965)](https://doi.org/10.1007/BF00186854).
+The crack advances by stepwise updates of $K_\textrm{I}$ and and crack tip position.
+
+The command does not take any arguments but can be configured using a `params.py` file in the current directory.
+An example `params.py` file follows:
+```Python
+import numpy as np
+
+import ase.io
+
+from matscipy.fracture_mechanics.clusters import diamond, set_groups
+
+from atomistica import TersoffScr, Tersoff_PRB_39_5566_Si_C__Scr
+
+###
+# Interaction potential
+calc = TersoffScr(**Tersoff_PRB_39_5566_Si_C__Scr)
+
+# Fundamental material properties
+el = 'C'
+a0 = 3.57
+
+# Elastic constants: either specify explicitly or compute automatically
+compute_elastic_constants = True
+# C11 = 1220.   # GPa
+# C12 = -3.    # GPa
+# C44 = 535.    # GPa
+
+# Surface energy
+surface_energy = 2.7326 * 10  # GPa*A = 0.1 J/m^2
+
+# Crack system
+crack_surface = [1, 1, 1]  # Normal of crack face
+crack_front = [1, -1, 0]  # Direction of crack front
+# bond = (10, 11)
+bondlength = 1.7  # Bond length for crack tip detection
+bulk_nn = 4  # Number of nearest neighbors in the bulk
+
+vacuum = 6.0  # Vacuum surrounding the cracked cluster
+
+# Simulation control
+nsteps = 31
+# Increase stress intensity factor
+k1 = np.linspace(0.8, 1.2, nsteps)
+# Don't move crack tip
+tip_dx = np.zeros_like(k1)
+tip_dz = np.zeros_like(k1)
+
+fmax = 0.05  # Tolerance for quasistatic optimizer
+
+# The actual crack system
+n = [1, 1, 1]
+skin_x, skin_y = 1, 1
+cryst = diamond(el, a0, n, crack_surface, crack_front)
+set_groups(cryst, n, skin_x, skin_y)  # Outer fixed atoms
+ase.io.write('cryst.xyz', cryst, format='extxyz')  # Dump initial crack system (without notch)
+```
+
+## `matscipy-sinclair-continuation`
+
+## `matscipy-sinclair-crack`

docs/cli/index.rst

@@ -14,5 +14,6 @@ whether a default parameter was used, to the log.
 
    diffusion
    electrochemistry
-   interatomic
+   fracture_mechanics
+   calculators
    structure

matscipy/tool/meson.build → matscipy/cli/fracture_mechanics/meson.build

@@ -3,11 +3,14 @@
 # Pure Python sources
 python_sources = [
     '__init__.py',
-    'plot.py'
+    'quasistatic_crack.py',
+    'setup_crack.py',
+    'sinclair_continuation.py',
+    'sinclair_crack.py'
 ]
 
 # Install pure Python
 python.install_sources(
     python_sources,
-    subdir: 'matscipy/tool'
+    subdir: 'matscipy/cli/fracture_mechanics'
 )

matscipy/cli/fracture_mechanics/quasistatic_crack.py

@@ -0,0 +1,177 @@
+#
+# Copyright 2014-2015, 2021 Lars Pastewka (U. Freiburg)
+#           2015 James Kermode (Warwick U.)
+#
+# matscipy - Materials science with Python at the atomic-scale
+# https://github.com/libAtoms/matscipy
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 2 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+
+# USAGE:
+#
+# Code imports the file 'params.py' from current working directory. params.py
+# contains simulation parameters. Some parameters can be omitted, see below.
+#
+# Parameters
+# ----------
+# calc : ase.Calculator
+#     Calculator object for energy and force computation.
+# tip_x0 : float
+#     Initial x-position of crack tip.
+# tip_y0 : float
+#     Initial y-position of crack tip.
+# tip_dx : array-like
+#     Displacement of tip in x-direction during run. x- and y-positions will be
+#     optimized self-consistently if omitted.
+# tip_dy : array-like
+#     Displacement of tip in y-direction during run. Position will be optimized
+#     self-consistently if omitted.
+
+import sys
+
+import numpy as np
+
+import ase.io
+import ase.constraints
+import ase.optimize
+
+from matscipy import parameter
+from matscipy.logger import screen
+
+from setup_crack import setup_crack
+from ase.optimize.precon import PreconLBFGS
+
+###
+
+def main():
+    # Atom types used for outputting the crack tip position.
+    ACTUAL_CRACK_TIP = 'Au'
+    FITTED_CRACK_TIP = 'Ag'
+
+    ###
+
+    logger = screen
+
+    ###
+
+    a, cryst, crk, k1g, tip_x0, tip_y0, bond1, bond2, boundary_mask, \
+        boundary_mask_bulk, tip_mask = setup_crack(logger=logger)
+    ase.io.write('notch.xyz', a, format='extxyz')
+
+    # Global parameters
+    basename = parameter('basename', 'quasistatic_crack')
+    calc = parameter('calc')
+    fmax = parameter('fmax', 0.01)
+
+    # Determine simulation control
+    k1_list = parameter('k1')
+    old_k1 = k1_list[0]
+    nsteps = len(k1_list)
+    tip_dx_list = parameter('tip_dx', np.zeros(nsteps))
+    tip_dy_list = parameter('tip_dy', np.zeros(nsteps))
+
+    # Run crack calculation.
+    tip_x = tip_x0
+    tip_y = tip_y0
+    a.set_calculator(calc)
+    for i, (k1, tip_dx, tip_dy) in enumerate(zip(k1_list, tip_dx_list,
+                                                 tip_dy_list)):
+        logger.pr('=== k1 = {0}*k1g, tip_dx = {1}, tip_dy = {2} ===' \
+                  .format(k1, tip_dx, tip_dy))
+        if tip_dx is None or tip_dy is None:
+            #
+            # Optimize crack tip position
+            #
+            old_y = tip_y + 1.0
+            old_x = tip_x + 1.0
+            while abs(tip_x - old_x) > 1e-6 and abs(tip_y - old_y) > 1e-6:
+                b = cryst.copy()
+                ux, uy = crk.displacements(cryst.positions[:, 0], cryst.positions[:, 1],
+                                           tip_x, tip_y, k * k1g)
+                b.positions[:, 0] += ux
+                b.positions[:, 1] += uy
+                a.set_constraint(None)
+                a.positions[boundary_mask] = b.positions[boundary_mask]
+                a.set_constraint(ase.constraints.FixAtoms(mask=boundary_mask))
+                logger.pr('Optimizing positions...')
+                opt = ase.optimize.FIRE(a, logfile=None)
+                opt.run(fmax=fmax)
+                logger.pr('...done. Converged within {0} steps.' \
+                          .format(opt.get_number_of_steps()))
+
+                old_x = tip_x
+                old_y = tip_y
+                tip_x, tip_y = crk.crack_tip_position(a.positions[:, 0],
+                                                      a.positions[:, 1],
+                                                      cryst.positions[:, 0],
+                                                      cryst.positions[:, 1],
+                                                      tip_x, tip_y, k * k1g,
+                                                      mask=mask)
+        else:
+            #
+            # Do not optimize tip position.
+            #
+            tip_x = tip_x0 + tip_dx
+            tip_y = tip_y0 + tip_dy
+            logger.pr('Setting crack tip position to {0} {1}' \
+                      .format(tip_x, tip_y))
+            # Scale strain field and optimize crack
+            a.set_constraint(None)
+            x, y = crk.scale_displacements(a.positions[:len(cryst), 0],
+                                           a.positions[:len(cryst), 1],
+                                           cryst.positions[:, 0],
+                                           cryst.positions[:, 1],
+                                           old_k1, k1)
+            a.positions[:len(cryst), 0] = x
+            a.positions[:len(cryst), 1] = y
+            # Optimize atoms in center
+            a.set_constraint(ase.constraints.FixAtoms(mask=boundary_mask))
+            logger.pr('Optimizing positions...')
+            opt = PreconLBFGS(a)
+            opt.run(fmax=fmax)
+            logger.pr('...done. Converged within {0} steps.' \
+                      .format(opt.get_number_of_steps()))
+            old_k1 = k1
+
+        # Output optimized configuration ot file. Include the mask array in
+        # output, so we know which atoms were used in fitting.
+        a.set_array('atoms_used_for_fitting_crack_tip', tip_mask)
+
+        # Output a file that contains the target crack tip (used for
+        # the displacmenets of the boundary atoms) and the fitted crack tip
+        # positions. The target crack tip is marked by a Hydrogen atom.
+        b = a.copy()
+        b += ase.Atom(ACTUAL_CRACK_TIP, (tip_x, tip_y, b.cell[2, 2] / 2))
+
+        # Measure the true (fitted) crack tip position.
+        try:
+            measured_tip_x, measured_tip_y = \
+                crk.crack_tip_position(a.positions[:, 0], a.positions[:, 1],
+                                       cryst.positions[:, 0], cryst.positions[:, 1],
+                                       tip_x, tip_y, k1 * k1g, mask=tip_mask)
+            measured_tip_x %= a.cell[0][0]
+            measured_tip_y %= a.cell[0][0]
+        except:
+            measured_tip_x = 0.0
+            measured_tip_y = 0.0
+
+        # The fitted crack tip is marked by a Helium atom.
+        b += ase.Atom(FITTED_CRACK_TIP, (measured_tip_x, measured_tip_y,
+                                         b.cell[2, 2] / 2))
+
+        b.info['bond_length'] = a.get_distance(bond1, bond2)
+        b.info['energy'] = a.get_potential_energy()
+        b.info['cell_origin'] = [0, 0, 0]
+        ase.io.write('%s_%4.4i.xyz' % (basename, i), b, format='extxyz')