@@ -0,0 +1,2602 @@
+import numpy as np%0Afrom numba import njit%0Afrom math import exp%0A%0A%0A@njit%0Adef two_body_energy(bond_array_1, c1, etypes1,%0A                    bond_array_2, c2, etypes2,%0A                    sig, ls, r_cut, cutoff_func):%0A    %22%22%22Compute the 2-body covariance between the local energy of one%0A        environment and the local energy, forces, and partial stresses (times%0A        the structure volume) of a second environment.%22%22%22%0A%0A    energy_kernel = 0.%0A    force_kernels = np.array(%5B0., 0., 0.%5D)%0A    stress_kernels = np.array(%5B0., 0., 0., 0., 0., 0.%5D)%0A%0A    ls1 = 1 / (2 * ls * ls)%0A    ls2 = 1 / (ls * ls)%0A    sig2 = sig * sig%0A%0A    for m in range(bond_array_1.shape%5B0%5D):%0A        ri = bond_array_1%5Bm, 0%5D%0A        fi, _ = cutoff_func(r_cut, ri, 1)%0A        e1 = etypes1%5Bm%5D%0A%0A        for n in range(bond_array_2.shape%5B0%5D):%0A            e2 = etypes2%5Bn%5D%0A%0A            # check if bonds agree%0A            if (c1 == c2 and e1 == e2) or (c1 == e2 and c2 == e1):%0A                rj = bond_array_2%5Bn, 0%5D%0A                xrel = bond_array_2%5Bn, 1%5D%0A                yrel = bond_array_2%5Bn, 2%5D%0A                zrel = bond_array_2%5Bn, 3%5D%0A                xval = xrel * rj%0A                yval = yrel * rj%0A                zval = zrel * rj%0A                fj, fdj = cutoff_func(r_cut, rj, 1)%0A%0A                r11 = ri - rj%0A%0A                # local energy kernel%0A                cons1 = r11 * r11%0A                cons2 = exp(-cons1 * ls1)%0A                energy_kernel += fi * fj * sig2 * cons2 / 4%0A%0A                # helper constants%0A                cons3 = cons2 * ls2 * fi * fj * r11%0A                cons4 = cons2 * fi * fdj%0A%0A                # fx + exx, exy, exz stress components%0A                fx = xrel * cons3 - cons4 * xrel%0A                force_kernels%5B0%5D += sig2 * fx / 2%0A                stress_kernels%5B0%5D -= sig2 * fx * xval / 4%0A                stress_kernels%5B1%5D -= sig2 * fx * yval / 4%0A                stress_kernels%5B2%5D -= sig2 * fx * zval / 4%0A%0A                # fy + eyy, eyz stress components%0A                fy = yrel * cons3 - cons4 * yrel%0A                force_kernels%5B1%5D += sig2 * fy / 2%0A                stress_kernels%5B3%5D -= sig2 * fy * yval / 4%0A                stress_kernels%5B4%5D -= sig2 * fy * zval / 4%0A%0A                # fz + ezz stress component%0A                fz = zrel * cons3 - cons4 * zrel%0A                force_kernels%5B2%5D += sig2 * fz / 2%0A                stress_kernels%5B5%5D -= sig2 * fz * zval / 4%0A%0A    return energy_kernel, force_kernels, stress_kernels%0A%0A%0A@njit%0Adef two_body_force(bond_array_1, c1, etypes1,%0A                   bond_array_2, c2, etypes2,%0A                   sig, ls, r_cut, cutoff_func):%0A    pass%0A

add two body energy kernel