How to use the matscipy.elasticity.fit_elastic_constants function in matscipy

def testtriclinic_relaxed(self):
            C, C_err = fit_elastic_constants(self.at0, 'triclinic',
                                             optimizer=FIRE, fmax=self.fmax,
                                             verbose=False, graphics=False)
            self.assertArrayAlmostEqual(C/units.GPa, self.C_ref_relaxed, tol=0.2)

def testtriclinic_unrelaxed(self):
            C, C_err = fit_elastic_constants(self.at0, 'triclinic',
                                             verbose=False, graphics=False)
            self.assertArrayAlmostEqual(C/units.GPa, self.C_ref, tol=0.2)

def test_Grochola(self):
        a = FaceCenteredCubic('Au', size=[2,2,2])
        calc = EAM('Au-Grochola-JCP05.eam.alloy')
        a.set_calculator(calc)
        FIRE(StrainFilter(a, mask=[1,1,1,0,0,0]), logfile=None).run(fmax=0.001)
        a0 = a.cell.diagonal().mean()/2
        self.assertTrue(abs(a0-4.0701)<2e-5)
        self.assertTrue(abs(a.get_potential_energy()/len(a)+3.924)<0.0003)
        C, C_err = fit_elastic_constants(a, symmetry='cubic', verbose=False)
        C11, C12, C44 = Voigt_6x6_to_cubic(C)
        self.assertTrue(abs((C11-C12)/GPa-32.07)<0.7)
        self.assertTrue(abs(C44/GPa-45.94)<0.5)

def testmonoclinic_unrelaxed(self):
            C, C_err = fit_elastic_constants(self.at0, 'monoclinic',
                                             verbose=False, graphics=False)
            self.assertArrayAlmostEqual(C/units.GPa, self.C_ref, tol=0.2)

def testorthorhombic_relaxed(self):
            C, C_err = fit_elastic_constants(self.at0, 'orthorhombic',
                                             optimizer=FIRE, fmax=self.fmax,
                                             verbose=False, graphics=False)
            self.assertArrayAlmostEqual(C/units.GPa, self.C_ref_relaxed, tol=0.1)

def testcubic_relaxed(self):
            C, C_err = fit_elastic_constants(self.at0, 'cubic',
                                             optimizer=FIRE, fmax=self.fmax,
                                             verbose=False, graphics=False)
            self.assertArrayAlmostEqual(C/units.GPa, self.C_ref_relaxed, tol=1e-2)
            if abs(C_err).max() > 0.:
                self.assertArrayAlmostEqual(C_err/units.GPa, self.C_err_ref_relaxed, tol=1e-2)

def testmonoclinic_relaxed(self):
            C, C_err = fit_elastic_constants(self.at0, 'monoclinic',
                                             optimizer=FIRE, fmax=self.fmax,
                                             verbose=False, graphics=False)
            self.assertArrayAlmostEqual(C/units.GPa, self.C_ref_relaxed, tol=0.2)

def testorthorhombic_unrelaxed(self):
            C, C_err = fit_elastic_constants(self.at0, 'orthorhombic',
                                             verbose=False, graphics=False)
            self.assertArrayAlmostEqual(C/units.GPa, self.C_ref, tol=0.1)

cryst = parameter('cryst').copy()
    cryst.set_pbc(True)
    
    # Double check elastic constants. We're just assuming this is really a periodic
    # system. (True if it comes out of the cluster routines.)
    
    compute_elastic_constants = parameter('compute_elastic_constants', False)
    elastic_fmax = parameter('elastic_fmax', 0.01)
    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'))

atom_types={'W': 1}, keep_alive=True)
        calculator = lammps

    unit_cell.set_calculator(calculator)

#   simple calculation to get the lattice constant and cohesive energy
#    alat0 = W.cell[0][1] - W.cell[0][0]
    sf = StrainFilter(unit_cell)  # or UnitCellFilter(W) -> to minimise wrt pos, cell
    opt = FIRE(sf)
    opt.run(fmax=1e-4)  # max force in eV/A
    alat = unit_cell.cell[0][1] - unit_cell.cell[0][0]
#    print("a0 relaxation %.4f --> %.4f" % (a0, a))
#    e_coh = W.get_potential_energy()
#    print("Cohesive energy %.4f" % e_coh)

    Cij, Cij_err = fit_elastic_constants(unit_cell,
                                         symmetry="cubic",
                                         delta=delta)

    Cij = Cij/GPa  # unit conversion to GPa

    elasticMatrix3x3 = Cij[:3, :3]
    # average of diagonal elements: C11, C22, C33
    C11 = elasticMatrix3x3.diagonal().mean()
    # make mask to extract non diagonal elements
    mask = np.ones((3, 3), dtype=bool)
    np.fill_diagonal(mask, False)

    # average of all non diagonal elements from 1 to 3
    C12 = elasticMatrix3x3[mask].mean()

    # average of diagonal elements from 4 till 6: C44, C55, C66,