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

a.positions[:,1],
                                                        r0[:,0], r0[:,1],
                                                        tip_x, tip_y, k1g,
                                                        mask=g!=0)

                print tip_x, tip_y

            # Get atomic strain
            i, j = neighbour_list("ij", a, cutoff=1.3)
            deformation_gradient, residual = \
                atomic_strain(a, ref, neighbours=(i, j))

            # Get atomic stresses
            # Note: get_stresses returns the virial in Atomistica!
            virial = a.get_stresses() 
            virial = Voigt_6_to_full_3x3_stress(virial)
            
            # Compute J-integral
            epot = a.get_potential_energies()
            eref = np.zeros_like(epot)

            for r1, r2 in [(r1, r2), (r1/2, r2/2), (r1/2, r2)]:
                print 'r1 = {}, r2 = {}'.format(r1, r2)

                J = J_integral(a, deformation_gradient, virial, epot, eref,
                               tip_x, tip_y, r1, r2)

                print '2*gamma = {0}, J = {1}'.format(2*surface_energy, J)

# Get atomic strain
i, j = neighbour_list("ij", a, cutoff=2.85)
deformation_gradient, residual = atomic_strain(a, ref, neighbours=(i, j))

# Get atomic stresses
virial = a.get_stresses() # Note: get_stresses returns the virial in Atomistica!
strain = full_3x3_to_Voigt_6_strain(deformation_gradient)
vol_strain, dev_strain, J3_strain = invariants(strain)

a.set_array('strain', strain)
a.set_array('vol_strain', vol_strain)
a.set_array('dev_strain', dev_strain)
a.set_array('virial', virial)

virial = Voigt_6_to_full_3x3_stress(virial)

# Coordination count
coord = np.bincount(i)
coord_coord = np.bincount(i, weights=coord[j])
mask = coord!=4
mask = np.ones_like(mask)

eref = e0*np.ones(len(a))
eref[coord_coord==15] = esurf1
eref[coord!=4] = esurf0

#virial[coord_coord==15] -= vsurf1
#virial[coord!=4] -= vsurf0

a.set_array('coordination', coord)
a.set_array('coord_coord', coord_coord)

if fix_params is None:
       fix_params = {}
    params.update(fix_params)

    x = atoms.positions[:, 0]
    y = atoms.positions[:, 1]
    r = np.sqrt((x - params['x0'])**2 + (y - params['y0'])**2)

    # Get local stresses
    if sigma is None:
       if calc is None:
           calc = atoms.get_calculator()
       sigma = calc.get_stresses(atoms)

    if sigma.shape != (len(atoms), 3, 3):
        sigma = Voigt_6_to_full_3x3_stress(sigma)

    # Update avg_sigma in place
    if avg_sigma is not None:
       avg_sigma[...] = np.exp(-avg_decay)*avg_sigma + (1.0 - np.exp(-avg_decay))*sigma
       sigma = avg_sigma.copy()

    # Zero components out of the xy plane
    sigma[:,2,2] = 0.0
    sigma[:,0,2] = 0.0
    sigma[:,2,0] = 0.0
    sigma[:,1,2] = 0.0
    sigma[:,2,1] = 0.0

    mask = Ellipsis # all atoms
    if r_range is not None:
        rmin, rmax = r_range

assert np.all(np.abs(forces[np.logical_not(mask)]-b.get_forces()) < 1e-6)

        #virial = params.calc.wpot_per_at
        deformation_gradient, residual = atomic_strain(b, ref, cutoff=2.85)
        virial = b.get_stresses()
        strain = full_3x3_to_Voigt_6_strain(deformation_gradient)
        #virial2 = strain.dot(C6)*vol0

        #vols = voropp(b)
        #virial3 = strain.dot(C6)*vols.reshape(-1,1)

        #print virial[175]
        #print virial2[175]
        #print virial3[175]

        virial = Voigt_6_to_full_3x3_stress(virial)
        #vol, dev, J3 = invariants(strain)
        #b.set_array('vol', vol)
        #b.set_array('dev', dev)
        x, y, z = b.positions.T
        r = np.sqrt((x-_tip_x)**2 + (y-_tip_y)**2)
        #b.set_array('J_eval', np.logical_and(r > params.eval_r1,
        #                                     r < params.eval_r2))

        #epot = b.get_potential_energies()
        G = J_integral(b, deformation_gradient, virial, epot_per_at, epotref_per_at,
                       _tip_x, _tip_y, (ref_sx+ref_sy)/8, 3*(ref_sx+ref_sy)/8)
    J_int += [ G ]

    #b.set_array('J_integral', np.logical_and(r > params.eval_r1,
    #                                         r < params.eval_r2))
    #ase.io.write('eval_'+fn, b, format='extxyz')

# Get surface energy and stress
a *= (1,3,1)
a.set_pbc([True, False, True])
ase.optimize.FIRE(a, logfile=None).run(fmax=params.fmax)
coord = np.bincount(neighbour_list("i", a, 2.85))
esurf0 = (a.get_potential_energy() - e0*len(a))
sx, sy, sz = a.cell.diagonal()
print 'surface energy = {} J/m^2'.format(esurf0/(2*sx*sz)/J_m2)
esurf0 = a.get_potential_energies()
esurf1 = esurf0[1]
esurf0 = esurf0[0]
print 'cohesive energy (surface atoms) = {}, {} eV/atom'.format(esurf0, esurf1)
vsurf0 = a.get_stresses()
vsurf1 = Voigt_6_to_full_3x3_stress(vsurf0[1])
vsurf0 = Voigt_6_to_full_3x3_stress(vsurf0[0])

# Reference configuration for strain calculation
a = ase.io.read(sys.argv[1])
ref = ase.io.read(sys.argv[2])

###

# Get crack tip position
tip_x, tip_y, tip_z = a.info['fitted_crack_tip']

# Set calculator
a.set_calculator(params.calc)

# Relax positions
del a[np.logical_or(a.numbers == atomic_numbers[ACTUAL_CRACK_TIP],
                    a.numbers == atomic_numbers[FITTED_CRACK_TIP])]

print 'cohesive energy = {} eV/atom'.format(e0)

# Get surface energy and stress
a *= (1,3,1)
a.set_pbc([True, False, True])
ase.optimize.FIRE(a, logfile=None).run(fmax=params.fmax)
coord = np.bincount(neighbour_list("i", a, 2.85))
esurf0 = (a.get_potential_energy() - e0*len(a))
sx, sy, sz = a.cell.diagonal()
print 'surface energy = {} J/m^2'.format(esurf0/(2*sx*sz)/J_m2)
esurf0 = a.get_potential_energies()
esurf1 = esurf0[1]
esurf0 = esurf0[0]
print 'cohesive energy (surface atoms) = {}, {} eV/atom'.format(esurf0, esurf1)
vsurf0 = a.get_stresses()
vsurf1 = Voigt_6_to_full_3x3_stress(vsurf0[1])
vsurf0 = Voigt_6_to_full_3x3_stress(vsurf0[0])

# Reference configuration for strain calculation
a = ase.io.read(sys.argv[1])
ref = ase.io.read(sys.argv[2])

###

# Get crack tip position
tip_x, tip_y, tip_z = a.info['fitted_crack_tip']

# Set calculator
a.set_calculator(params.calc)

# Relax positions
del a[np.logical_or(a.numbers == atomic_numbers[ACTUAL_CRACK_TIP],