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)
```python
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:

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

Additional change:

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