Rewrite b2 lmp

[YuuuXie]

The pair style with summation order reorganized. Now has passed unit test.
Some timing results with SiC:

• Our previous implementation (unit: s)

function single_bond [time=1.76481]

function B2_descriptor [time=13.68479]
    for n1, n2, l, m 
        B(n1, n2, l) += A(n1, l, m) * A(n2, l, m)  
        for neigh_k
             dB(n1, n2, l)/dr_ik

beta_matrix * B (matrix-vector multiplication) [time=6.23888]

[total time=21.68848]

• This PR (with more details below)

function single_bond [time=1.636844]

function B2_descriptor [time=7.245533]
    for n1, n2, l, m
	B(n1,n2,l) += A(n1,l,m) * A(n2,l,m)

    w = beta_matrix * B (matrix-vector multiplication) # originally outside of the function B2_descriptor

    for n1, n2, l, m
	u(n1,l,m) += w(n1,n2,l) * A(n2,l,m)

compute partial forces [time=1.807217]
    for n1, l, m, neigh_k # In practice the loop is only a matrix-vector multiplication u * dA
	f_ik += 2 * u(n1,l,m) * dA(n1,l,m) / dr_ik

[total time=10.689594]

TODO:

• Fix variance compute command with the new B2_descriptor function
• remove time stamps in pair style

Additional change:

• Divide the total energy kernel by the number of atoms for normalization

[jonpvandermause]

We should have Github actions run the Python tests too (might involve compiling lammps)