Fix: various matrix elements calculations in orbital generation

[jinzx10]

Reminder

• [x] Have you linked an issue with this pull request?
• [x] Have you added adequate unit tests and/or case tests for your pull request?
• [x] Have you noticed possible changes of behavior below or in the linked issue?
• [x] Have you explained the changes of codes in core modules of ESolver, HSolver, ElecState, Hamilt, Operator or Psi? (ignore if not applicable)

Linked Issue

Fix #3618 Resolve #3616

Unit Tests and/or Case Tests for my changes

• A unit test for two-center integral <jY|jY> is added

What's changed?

• The algorithm for computing <jY|jY> in orbital generation is switched from a summation in plane-wave basis to two-center integration. This change is expected to bring 2-3 more significant digits. A side effect of this change is that the previous code missed a factor of 2*pi/lat0 in several places, which resulted in an effective 10x weight for the gradient term in spillage. After this change, the same spillage optimization would need to manually set the weight to reproduce previous calculations.
• a complex phase in the plane-wave explansion in the <mo|jY> calculation is corrected;
• bugs in the polynomial interpolation in bessel_basis.cpp are fixed
• removes the usage of ucell via GlobalC and simplifies several function interfaces in numerical_basis

Any changes of core modules? (ignore if not applicable)

No.

[jinzx10]

This PR has been converted to a draft since its consequence to subsequent spillage optimization is not fully studied yet. It might be better to merge this PR after sufficient tests on the final orbitals. @kirk0830