How to use the matscipy.fracture_mechanics.crack function in matscipy

reftip_y = tip_y
                else:
                    a.set_cell(refcell)

                # Shift tip position so all systems are exactly centered at the same spot
                a.positions[:, 0] += reftip_x - tip_x
                a.positions[:, 1] += reftip_y - tip_y

                refpositions = a.positions.copy()

                # Move reference crystal by same amount
                cryst.set_cell(a.cell)
                cryst.set_pbc([False, False, True])
                cryst.translate(a[0].position - oldr)

                bond1, bond2 = crack.find_tip_coordination(a, bondlength=bondlength*1.2)

                # Groups mark the fixed region and the region use for fitting the crack tip.
                g = a.get_array('groups')
                gcryst = cryst.get_array('groups')

                ase.io.write('cryst_{}.xyz'.format(i), cryst)

                a.set_calculator(calc)
                a.set_constraint(FixAtoms(mask=g==0))
                FIRE(a, logfile=None).run(fmax=1e-6)

                dpos = np.sqrt(((a.positions[:, 0]-refpositions[:, 0])/ux)**2 + ((a.positions[:, 1]-refpositions[:, 1])/uy)**2)
                a.set_array('dpos', dpos)

                distance_from_tip = np.sqrt((a.positions[:, 0]-reftip_x)**2 + (a.positions[:, 1]-reftip_y)**2)

###

a, cryst, crk, k1g, tip_x, tip_y, bond1, bond2, boundary_mask, \
    boundary_mask_bulk, tip_mask = setup_crack(logger=logger)
ase.io.write('notch.xyz', a, format='extxyz')

# Get general parameters

basename = parameter('basename', 'crack_tip')
calc = parameter('calc')
fmax = parameter('fmax', 0.01)

# Get parameter used for fitting crack tip position

residual_func = parameter('residual_func', crack.displacement_residual)
_residual_func = residual_func
tip_tol = parameter('tip_tol', 1e-4)
tip_mixing_alpha = parameter('tip_mixing_alpha', 1.0)
write_trajectory_during_optimization = parameter('write_trajectory_during_optimization', False)

tip_x = (a.positions[bond1, 0] + a.positions[bond2, 0])/2
tip_y = (a.positions[bond1, 1] + a.positions[bond2, 1])/2
logger.pr('Optimizing tip position -> initially centering tip bond. '
          'Tip positions = {} {}'.format(tip_x, tip_y))

# Check if there is a request to restart from a positions file
restart_from = parameter('restart_from', 'N/A')
if restart_from != 'N/A':
    logger.pr('Restarting from {0}'.format(restart_from))
    a = ase.io.read(restart_from)
    # remove any marker atoms

###

sys.path += [ "." ]
import params

###

# Atom types used for outputting the crack tip position.
ACTUAL_CRACK_TIP = 'Au'
FITTED_CRACK_TIP = 'Ag'

###

cryst = params.cryst.copy()
crk = crack.CubicCrystalCrack(params.C11, params.C12, params.C44,
                              params.crack_surface, params.crack_front)

# Get Griffith's k1.
k1g = crk.k1g(params.surface_energy)
parprint('Griffith k1 = %f' % k1g)

# Crack tip position.
if hasattr(params, 'tip_x'):
    tip_x = params.tip_x
else:
    tip_x = cryst.cell.diagonal()[0]/2
if hasattr(params, 'tip_y'):
    tip_y = params.tip_y
else:
    tip_y = cryst.cell.diagonal()[1]/2

elastic_symmetry = parameter('elastic_symmetry', 'triclinic')
    elastic_optimizer = parameter('elastic_optimizer', ase.optimize.FIRE)

    if compute_elastic_constants:
        cryst.set_calculator(calc)
        log_file = open('elastic_constants.log', 'w')
        C, C_err = fit_elastic_constants(cryst, verbose=False,
                                         symmetry=elastic_symmetry,
                                         optimizer=elastic_optimizer,
                                         logfile=log_file,
                                         fmax=elastic_fmax)
        log_file.close()
        logger.pr('Measured elastic constants (in GPa):')
        logger.pr(np.round(C*10/GPa)/10)
    
        crk = crack.CubicCrystalCrack(parameter('crack_surface'),
                                      parameter('crack_front'),
                                      Crot=C/GPa)
    else:
        if has_parameter('C'):
            crk = crack.CubicCrystalCrack(parameter('crack_surface'),
                                          parameter('crack_front'),
                                          C=parameter('C'))
        else:    
            crk = crack.CubicCrystalCrack(parameter('crack_surface'),
                                          parameter('crack_front'),
                                          parameter('C11'), parameter('C12'),
                                          parameter('C44'))
    
    
    logger.pr('Elastic constants used for boundary condition (in GPa):')
    logger.pr(np.round(crk.C*10)/10)

residual_func = parameter('residual_func', crack.displacement_residual)
_residual_func = residual_func

basename = parameter('basename', 'neb')
calc = parameter('calc')
fmax_neb = parameter('fmax_neb', 0.1)
maxsteps_neb = parameter('maxsteps_neb', 100)
Nimages = parameter('Nimages', 7)
k_neb = parameter('k_neb', 1.0)

a, cryst, crk, k1g, tip_x, tip_y, bond1, bond2, boundary_mask, \
    boundary_mask_bulk, tip_mask = setup_crack(logger=logger)

# Deformation gradient residual needs full Atoms object and therefore
# special treatment here.
if _residual_func == crack.deformation_gradient_residual:
    residual_func = lambda r0, crack, x, y, ref_x, ref_y, k, mask=None:\
        _residual_func(r0, crack, x, y, a, ref_x, ref_y, cryst, k,
                       params.cutoff, mask)

# Groups
g = a.get_array('groups')

dirname = os.path.basename(os.getcwd())
initial = ase.io.read('../../initial_state/{0}/initial_state.xyz'.format(dirname))

final_pos = ase.io.read('../../final_state/{0}/final_state.xyz'.format(dirname))
final = initial.copy()
final.set_positions(final_pos.positions)

# remove marker atoms
mask = np.logical_and(initial.numbers != atomic_numbers[ACTUAL_CRACK_TIP],

# Assign calculator.
a.set_calculator(calc)

log_file = open('{0}.log'.format(basename), 'w')
if write_trajectory_during_optimization:
    traj_file = ase.io.NetCDFTrajectory('{0}.nc'.format(basename), mode='w',
                                        atoms=a)
    traj_file.write()
else:
    traj_file = None

# Deformation gradient residual needs full Atoms object and therefore
# special treatment here.
if _residual_func == crack.deformation_gradient_residual:
    residual_func = lambda r0, crack, x, y, ref_x, ref_y, k, mask=None:\
        _residual_func(r0, crack, x, y, a, ref_x, ref_y, cryst, k,
                       params.cutoff, mask)



old_x = tip_x
old_y = tip_y
converged = False
while not converged:
    #b = cryst.copy()
    u0x, u0y = crk.displacements(cryst.positions[:,0],
                                 cryst.positions[:,1],
                                 old_x, old_y, params.k1*k1g)
    ux, uy = crk.displacements(cryst.positions[:,0],
                               cryst.positions[:,1],

tip_y += a[0].position[1] - oldr[1]
cryst.set_cell(a.cell)
cryst.translate(a[0].position - oldr)
ase.io.write('notch.xyz', a, format='extxyz')

parprint('Initial tip position {0} {1}'.format(tip_x, tip_y))

# Groups mark the fixed region and the region use for fitting the crack tip.
g = a.get_array('groups')

# Assign calculator.
a.set_calculator(params.calc)

for j in range(params.n_bonds):

    bond1, bond2 = crack.find_tip_coordination(a)
    print 'Breaking tip bond {0}--{1}'.format(bond1, bond2)

    # Run crack calculation.
    for i, bond_length in enumerate(params.bond_lengths):
        parprint('=== bond_length = {0} ==='.format(bond_length))
        xyz_file = 'bond_%d_%4d.xyz' % (j+1, int(bond_length*1000))
        if os.path.exists(xyz_file):
            parprint('%s found, skipping' % xyz_file)
            a = ase.io.read(xyz_file)
            del a[np.logical_or(a.numbers == atomic_numbers[ACTUAL_CRACK_TIP],
                                a.numbers == atomic_numbers[FITTED_CRACK_TIP])]
            a.set_calculator(params.calc)
        else:
            a.set_constraint(None)
            a.set_distance(bond1, bond2, bond_length)
            bond_length_constraint = ase.constraints.FixBondLength(bond1, bond2)

ux, uy = crk.displacements(cryst.positions[:,0], cryst.positions[:,1],
                               tip_x, tip_y, k1*k1g)
    a.positions[:len(cryst),0] += ux
    a.positions[:len(cryst),1] += uy
       
    # Center notched configuration in simulation cell and ensure enough vacuum.
    oldr = a[0].position.copy()
    vacuum = parameter('vacuum')
    a.center(vacuum=vacuum, axis=0)
    a.center(vacuum=vacuum, axis=1)
    tip_x += a[0].x - oldr[0]
    tip_y += a[0].y - oldr[1]
    
    # Choose which bond to break.
    bond1, bond2 = \
        parameter('bond', crack.find_tip_coordination(a, bondlength=bondlength))

    if parameter('center_crack_tip_on_bond', False):
        tip_x, tip_z, dummy = (a.positions[bond1]+a.positions[bond2])/2
    
    # Hydrogenate?
    coord = np.bincount(neighbour_list('i', a, bondlength), minlength=len(a))
    a.set_array('coord', coord)
    
    if parameter('optimize_full_crack_face', False):
        g = a.get_array('groups')
        gcryst = cryst.get_array('groups')
        coord = a.get_array('coord')
        g[coord!=4] = -1
        gcryst[coord!=4] = -1
        a.set_array('groups', g)
        cryst.set_array('groups', gcryst)

Optimizer = ase.optimize.FIRE
#Optimizer = ase.optimize.precon.LBFGS

# Atom types used for outputting the crack tip position.
ACTUAL_CRACK_TIP = 'Au'
FITTED_CRACK_TIP = 'Ag'

###

logger = screen

###

# Get general parameters

residual_func = parameter('residual_func', crack.displacement_residual)
_residual_func = residual_func

basename = parameter('basename', 'neb')
calc = parameter('calc')
fmax_neb = parameter('fmax_neb', 0.1)
maxsteps_neb = parameter('maxsteps_neb', 100)
Nimages = parameter('Nimages', 7)
k_neb = parameter('k_neb', 1.0)

a, cryst, crk, k1g, tip_x, tip_y, bond1, bond2, boundary_mask, \
    boundary_mask_bulk, tip_mask = setup_crack(logger=logger)

# Deformation gradient residual needs full Atoms object and therefore
# special treatment here.
if _residual_func == crack.deformation_gradient_residual:
    residual_func = lambda r0, crack, x, y, ref_x, ref_y, k, mask=None:\

print(len(a), len(b))
            print(a.numbers)
            print(b.numbers)
            assert np.all(a.numbers == b.numbers)
            a = b
            a.set_calculator(calc)
            tip_x, tip_y, dummy = a.info['actual_crack_tip']
            a.set_cell(original_cell, scale_atoms=True)

    a.set_constraint(None)
    a.set_distance(bond1, bond2, bond_length)
    bond_length_constraint = FixBondLength(bond1, bond2)

    # Deformation gradient residual needs full Atoms object and therefore
    # special treatment here.
    if _residual_func == crack.deformation_gradient_residual:
        residual_func = lambda r0, crack, x, y, ref_x, ref_y, k, mask=None:\
            _residual_func(r0, crack, x, y, a, ref_x, ref_y, cryst, k,
                           params.cutoff, mask)

    # Optimize x and z position of crack tip.
    if optimize_tip_position:
        old_x = tip_x
        old_y = tip_y
        converged = False
        while not converged:
            #b = cryst.copy()
            u0x, u0y = crk.displacements(cryst.positions[:,0],
                                         cryst.positions[:,1],
                                         old_x, old_y, params.k1*k1g)
            ux, uy = crk.displacements(cryst.positions[:,0],
                                       cryst.positions[:,1],