How to use matscipy - 10 common examples

def test_out_of_bounds(self):
        nat = 10
        atoms = ase.Atoms(numbers=range(nat),
                          cell=[(0.2, 1.2, 1.4),
                                (1.4, 0.1, 1.6),
                                (1.3, 2.0, -0.1)])
        atoms.set_scaled_positions(3 * np.random.random((nat, 3)) - 1)
        
        for p1 in range(2):
            for p2 in range(2):
                for p3 in range(2):
                    atoms.set_pbc((p1, p2, p3))
                    i, j, d, D, S = neighbour_list("ijdDS", atoms, atoms.numbers * 0.2 + 0.5)
                    c = np.bincount(i)
                    atoms2 = atoms.repeat((p1 + 1, p2 + 1, p3 + 1))
                    i2, j2, d2, D2, S2 = neighbour_list("ijdDS", atoms2, atoms2.numbers * 0.2 + 0.5)
                    c2 = np.bincount(i2)
                    c2.shape = (-1, nat)
                    dd = d.sum() * (p1 + 1) * (p2 + 1) * (p3 + 1) - d2.sum()
                    dr = np.linalg.solve(atoms.cell.T, (atoms.positions[1]-atoms.positions[0]).T).T+np.array([0,0,3])
                    self.assertTrue(abs(dd) < 1e-10)
                    self.assertTrue(not (c2 - c).any())

def test_small_cell(self):
        a = ase.Atoms('C', positions=[[0.5, 0.5, 0.5]], cell=[1, 1, 1],
                      pbc=True)
        i, j, dr, shift = neighbour_list("ijDS", a, 1.1)
        assert np.bincount(i)[0] == 6
        assert (dr == shift).all()

        i, j = neighbour_list("ij", a, 1.5)
        assert np.bincount(i)[0] == 18

        a.set_pbc(False)
        i = neighbour_list("i", a, 1.1)
        assert len(i) == 0

        a.set_pbc([True, False, False])
        i = neighbour_list("i", a, 1.1)
        assert np.bincount(i)[0] == 2

        a.set_pbc([True, False, True])
        i = neighbour_list("i", a, 1.1)
        assert np.bincount(i)[0] == 4

def test_small_cell(self):
        a = ase.Atoms('C', positions=[[0.5, 0.5, 0.5]], cell=[1, 1, 1],
                      pbc=True)
        i, j, dr, shift = neighbour_list("ijDS", a, 1.1)
        assert np.bincount(i)[0] == 6
        assert (dr == shift).all()

        i, j = neighbour_list("ij", a, 1.5)
        assert np.bincount(i)[0] == 18

        a.set_pbc(False)
        i = neighbour_list("i", a, 1.1)
        assert len(i) == 0

        a.set_pbc([True, False, False])
        i = neighbour_list("i", a, 1.1)
        assert np.bincount(i)[0] == 2

        a.set_pbc([True, False, True])
        i = neighbour_list("i", a, 1.1)
        assert np.bincount(i)[0] == 4

def test_hydrogenate(self):
        a = Diamond('Si', size=[2,2,1])
        b = hydrogenate(a, 2.85, 1.0, mask=[True,True,False], vacuum=5.0)
        # Check if every atom is fourfold coordinated
        syms = np.array(b.get_chemical_symbols())
        c = np.bincount(neighbour_list('i', b, 2.4))
        self.assertTrue((c[syms!='H']==4).all())

def test_consistency_of_deformation_gradient_and_stress(self):
        for C11, C12, C44, surface_energy, k1 in self.materials:
            crack = CubicCrystalCrack([1,0,0], [0,1,0], C11, C12, C44)
            k = crack.k1g(surface_energy)*k1
            for i in range(10):
                x = np.random.uniform(-10, 10)
                y = np.random.uniform(-10, 10)
                F = crack.deformation_gradient(x, y, 0, 0, k)
                eps = (F+F.T)/2-np.eye(2)
                sig_x, sig_y, sig_xy = crack.stresses(x, y, 0, 0, k)
                eps_xx = crack.crack.a11*sig_x + \
                         crack.crack.a12*sig_y + \
                         crack.crack.a16*sig_xy
                eps_yy = crack.crack.a12*sig_x + \
                         crack.crack.a22*sig_y + \
                         crack.crack.a26*sig_xy
                # Factor 1/2 because elastic constants and matrix product are
                # expressed in Voigt notation.
                eps_xy = (crack.crack.a16*sig_x + \

refcell = None
            reftip_x = None
            reftip_y = None
            #[41, 39, 1], 
            for i, n in enumerate([[21, 19, 1], [11, 9, 1], [6, 5, 1]]):
                #print(n)
                cryst = diamond(el, a0, n, crack_surface, crack_front)
                set_groups(cryst, n, skin_x, skin_y)
                cryst.set_pbc(True)
                cryst.set_calculator(calc)
                FIRE(UnitCellFilter(cryst), logfile=None).run(fmax=1e-6)
                cryst.set_cell(cryst.cell.diagonal(), scale_atoms=True)

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

                crk = crack.CubicCrystalCrack(crack_surface,
                                              crack_front,
                                              Crot=C/units.GPa)
                k1g = crk.k1g(surface_energy)

                tip_x = cryst.cell.diagonal()[0]/2
                tip_y = cryst.cell.diagonal()[1]/2

                a = cryst.copy()
                a.set_pbc([False, False, True])

                k1 = 1.0
                ux, uy = crk.displacements(cryst.positions[:,0], cryst.positions[:,1],
                                           tip_x, tip_y, k1*k1g)

                a.positions[:, 0] += ux
                a.positions[:, 1] += uy

def test_elastostatics(self):
        eps = 1e-3
        for C11, C12, C44, surface_energy, k1 in self.materials:
            crack = CubicCrystalCrack([1,0,0], [0,1,0], C11, C12, C44)
            k = crack.k1g(surface_energy)*k1
            for i in range(10):
                x = i+1
                y = i+1
                # Finite difference approximation of the stress divergence
                sxx0x, syy0x, sxy0x = crack.stresses(x-eps, y, 0, 0, k)
                sxx0y, syy0y, sxy0y = crack.stresses(x, y-eps, 0, 0, k)
                sxx1x, syy1x, sxy1x = crack.stresses(x+eps, y, 0, 0, k)
                sxx1y, syy1y, sxy1y = crack.stresses(x, y+eps, 0, 0, k)
                divsx = (sxx1x-sxx0x)/(2*eps) + (sxy1y-sxy0y)/(2*eps)
                divsy = (sxy1x-sxy0x)/(2*eps) + (syy1y-syy0y)/(2*eps)
                # Check that divergence of stress is zero (elastostatic
                # equilibrium)
                self.assertAlmostEqual(divsx, 0.0, places=3)
                self.assertAlmostEqual(divsy, 0.0, places=3)

return

        nx = 128
        for calc, a0, C11, C12, C44, surface_energy, bulk_coordination in [
            ( atomistica.DoubleHarmonic(k1=1.0, r1=1.0, k2=1.0,
                                        r2=math.sqrt(2), cutoff=1.6),
              clusters.sc('He', 1.0, [nx,nx,1], [1,0,0], [0,1,0]),
              3, 1, 1, 0.05, 6 ),
#            ( atomistica.Harmonic(k=1.0, r0=1.0, cutoff=1.3, shift=False),
#              clusters.fcc('He', math.sqrt(2.0), [nx/2,nx/2,1], [1,0,0],
#                           [0,1,0]),
#              math.sqrt(2), 1.0/math.sqrt(2), 1.0/math.sqrt(2), 0.05, 12)
            ]:
            print '{} atoms.'.format(len(a0))

            crack = CubicCrystalCrack(C11, C12, C44, [1,0,0], [0,1,0])

            x, y, z = a0.positions.T
            r2 = min(np.max(x)-np.min(x), np.max(y)-np.min(y))/4
            r1 = r2/2

            a = a0.copy()
            a.center(vacuum=20.0, axis=0)
            a.center(vacuum=20.0, axis=1)
            ref = a.copy()
            r0 = ref.positions

            a.set_calculator(calc)
            print 'epot = {}'.format(a.get_potential_energy())

            sx, sy, sz = a.cell.diagonal()
            tip_x = sx/2

def test_J_integral(self):
        if quadrature is None:
            print('No scipy.integrate.quadrature. Skipping J-integral test.')
            return

        for C11, C12, C44, surface_energy, k1 in self.materials:
            crack = CubicCrystalCrack([1,0,0], [0,1,0], C11, C12, C44)
            k = crack.k1g(surface_energy)*k1

            def polar_path(theta, r=1, x0=0, y0=0):
                nx = np.cos(theta)
                ny = np.sin(theta)
                n = np.transpose([nx, ny])
                return r*n-np.array([x0, y0]), n, r

            def elliptic_path(theta, r=1, x0=0, y0=0):
                rx, ry = r
                x = rx*np.cos(theta)
                y = ry*np.sin(theta)
                nx = ry*np.cos(theta)
                ny = rx*np.sin(theta)
                ln = np.sqrt(nx**2+ny**2)
                nx /= ln

def test_consistency_of_deformation_gradient_and_displacement(self):
        eps = 1e-3
        for C11, C12, C44, surface_energy, k1 in self.materials:
            crack = CubicCrystalCrack([1,0,0], [0,1,0], C11, C12, C44)
            k = crack.k1g(surface_energy)*k1
            for i in range(10):
                x = i+1
                y = i+1
                F = crack.deformation_gradient(x, y, 0, 0, k)
                # Finite difference approximation of deformation gradient
                u, v = crack.displacements(x, y, 0, 0, k)
                ux0, vx0 = crack.displacements(x-eps, y, 0, 0, k)
                uy0, vy0 = crack.displacements(x, y-eps, 0, 0, k)
                ux1, vx1 = crack.displacements(x+eps, y, 0, 0, k)
                uy1, vy1 = crack.displacements(x, y+eps, 0, 0, k)
                du_dx = (ux1-ux0)/(2*eps)
                du_dy = (uy1-uy0)/(2*eps)
                dv_dx = (vx1-vx0)/(2*eps)
                dv_dy = (vy1-vy0)/(2*eps)
                du_dx += np.ones_like(du_dx)