Evauation of angular momentum in GBasis vs. Libcint bindings

[msricher]

There's a difference in some positive and negative signs in the outputs.

tests/test_libcint.py::test_integral[STO-6G-Be-Cartesian-AngularMomentum] 
GBASIS
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00-1.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+1.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00-1.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00-1.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+1.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00-1.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
CBASIS
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+1.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00-1.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00-1.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+1.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+1.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00-1.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
FAILED

Originally posted by @msricher in #137 (comment)

Attached image is what the Gbasis notes say angular momentum is, which matches up with the code as far as I can tell.

Libcint doesn't give as much documentation, but it says it's r cross Nabla; does this exclusively imply the formula in the image attached, or could there be a difference of convention? I feel like a difference of convention could explain why there's a factor of -1 difference in some elements, but not others.

Is there anyway I can evaluate this array numerically using Grid, to verify whether GBasis or Libcint is correct?

[msricher]

Single P orbital, cartesian

GBASIS
( 0,  0): 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
( 0,  1): 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+1.000000000e+00j
( 0,  2): 0.000000000e+00+0.000000000e+00j, 0.000000000e+00-1.000000000e+00j, 0.000000000e+00+0.000000000e+00j
( 1,  0): 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00-1.000000000e+00j
( 1,  1): 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
( 1,  2): 0.000000000e+00+1.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
( 2,  0): 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+1.000000000e+00j, 0.000000000e+00+0.000000000e+00j
( 2,  1): 0.000000000e+00-1.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
( 2,  2): 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
CBASIS
( 0,  0): 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
( 0,  1): 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00-1.000000000e+00j
( 0,  2): 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+1.000000000e+00j, 0.000000000e+00+0.000000000e+00j
( 1,  0): 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00-1.000000000e+00j
( 1,  1): 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
( 1,  2): 0.000000000e+00-1.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
( 2,  0): 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+1.000000000e+00j, 0.000000000e+00+0.000000000e+00j
( 2,  1): 0.000000000e+00-1.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
( 2,  2): 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j

[Ali-Tehrani]

Judging from wikipedia, that formula from the image (with the minus sign) match.

You can use grid it to evaluate this integral. The first grid I would recommend is a atomic grid centered at the p-orbital center.

from grid.molgrid import AtomGrid
atcoords = # p-orbital center 
atgrid = AtomGrid.from_preset(atnum=3, preset="fine", center=atcoords)

Then you can interpolate the p-orbital then approximate its derivative

p_orbital_vals = # Evaluated on atgrid.points
func_interpolate_dens = atgrid.interpolate(p_orbital_vals)
deriv_density = func_interpolate_dens(atgrid.points, deriv=1)  # Get the first derivative wrt to Cartesian

Then you can evaluate those integral by running

x, y, z = atgrid.points.T
atgrid.integrate(p_orbital_vals * y * deriv_density[:, 1])  # etc 

The other option is to analytically calculate the derivatives. If the accuracy isn't good, you can play around increasing the preset parameter or adding your own radial grid.

[msricher]

@FarnazH @PaulWAyers

@leila-pujal and I have tested the angular momentum integrals using Grid. It seems Gbasis is correct. Now, since I can't find a proper LibCint function to call, in order to get the right answer, should we just ask someone or raise an Issue in the Libcint repo?

[PaulWAyers]

Let's make sure it isn't a different ordering. I can give some tests for gbasis too.

I know that libcint does a strange ordering, so that can be a factor. One other idea, however: did you confirm that gbasis and libcint are also different for Cartesians, or is this specific to Spherical?

[msricher]

I will drop both spherical/cartesian here.

CARTESIAN

GBASIS
( 0,  0): +0.000000000e+00, +0.000000000e+00, +0.000000000e+00
( 0,  1): +0.000000000e+00, +0.000000000e+00, +1.000000000e+00
( 0,  2): +0.000000000e+00, -1.000000000e+00, +0.000000000e+00
( 1,  0): +0.000000000e+00, +0.000000000e+00, -1.000000000e+00
( 1,  1): +0.000000000e+00, +0.000000000e+00, +0.000000000e+00
( 1,  2): +1.000000000e+00, +0.000000000e+00, +0.000000000e+00
( 2,  0): +0.000000000e+00, +1.000000000e+00, +0.000000000e+00
( 2,  1): -1.000000000e+00, +0.000000000e+00, +0.000000000e+00
( 2,  2): +0.000000000e+00, +0.000000000e+00, +0.000000000e+00
CBASIS
( 0,  0): +0.000000000e+00, +0.000000000e+00, +0.000000000e+00
( 0,  1): +0.000000000e+00, +0.000000000e+00, +1.000000000e+00
( 0,  2): +0.000000000e+00, -1.000000000e+00, +0.000000000e+00
( 1,  0): +0.000000000e+00, +0.000000000e+00, +1.000000000e+00
( 1,  1): +0.000000000e+00, +0.000000000e+00, +0.000000000e+00
( 1,  2): +1.000000000e+00, +0.000000000e+00, +0.000000000e+00
( 2,  0): +0.000000000e+00, -1.000000000e+00, +0.000000000e+00
( 2,  1): +1.000000000e+00, +0.000000000e+00, +0.000000000e+00
( 2,  2): +0.000000000e+00, +0.000000000e+00, +0.000000000e+00
GRID
( 0,  0): +0.000000000e+00, +0.000000000e+00, +0.000000000e+00
( 0,  1): +0.000000000e+00, +0.000000000e+00, +1.000000000e+00
( 0,  2): +0.000000000e+00, -1.000000000e+00, +0.000000000e+00
( 1,  0): +0.000000000e+00, +0.000000000e+00, -1.000000000e+00
( 1,  1): +0.000000000e+00, +0.000000000e+00, +0.000000000e+00
( 1,  2): +1.000000000e+00, +0.000000000e+00, +0.000000000e+00
( 2,  0): +0.000000000e+00, +9.435335430e-01, +0.000000000e+00
( 2,  1): -9.717667715e-01, +0.000000000e+00, +0.000000000e+00
( 2,  2): +0.000000000e+00, +0.000000000e+00, +0.000000000e+00

SPHERICAL

GBASIS
( 0,  0): +0.000000000e+00, +0.000000000e+00, +0.000000000e+00
( 0,  1): +1.000000000e+00, +0.000000000e+00, +0.000000000e+00
( 0,  2): +0.000000000e+00, +0.000000000e+00, -1.000000000e+00
( 1,  0): -1.000000000e+00, +0.000000000e+00, +0.000000000e+00
( 1,  1): +0.000000000e+00, +0.000000000e+00, +0.000000000e+00
( 1,  2): +0.000000000e+00, +1.000000000e+00, +0.000000000e+00
( 2,  0): +0.000000000e+00, +0.000000000e+00, +1.000000000e+00
( 2,  1): +0.000000000e+00, -1.000000000e+00, +0.000000000e+00
( 2,  2): +0.000000000e+00, +0.000000000e+00, +0.000000000e+00
CBASIS
( 0,  0): +0.000000000e+00, -0.000000000e+00, -0.000000000e+00
( 0,  1): +1.000000000e+00, -0.000000000e+00, -0.000000000e+00
( 0,  2): +0.000000000e+00, -0.000000000e+00, -1.000000000e+00
( 1,  0): +1.000000000e+00, -0.000000000e+00, -0.000000000e+00
( 1,  1): -0.000000000e+00, +0.000000000e+00, -0.000000000e+00
( 1,  2): -0.000000000e+00, +1.000000000e+00, -0.000000000e+00
( 2,  0): +0.000000000e+00, -0.000000000e+00, -1.000000000e+00
( 2,  1): -0.000000000e+00, +1.000000000e+00, -0.000000000e+00
( 2,  2): -0.000000000e+00, -0.000000000e+00, +0.000000000e+00
GRID
( 0,  0): +0.000000000e+00, +0.000000000e+00, +0.000000000e+00
( 0,  1): +1.000000000e+00, +0.000000000e+00, +0.000000000e+00
( 0,  2): +0.000000000e+00, +0.000000000e+00, -1.000000000e+00
( 1,  0): -9.717667715e-01, +0.000000000e+00, +0.000000000e+00
( 1,  1): +0.000000000e+00, +0.000000000e+00, +0.000000000e+00
( 1,  2): +0.000000000e+00, +9.435335430e-01, +0.000000000e+00
( 2,  0): +0.000000000e+00, +0.000000000e+00, +1.000000000e+00
( 2,  1): +0.000000000e+00, -1.000000000e+00, +0.000000000e+00
( 2,  2): +0.000000000e+00, +0.000000000e+00, +0.000000000e+00

The answers are consistent for spherical/cartesian. Libcint is wrong here. At least this integral isn't angular momentum, and I can't find one that is correct.

I made sure the P orbital ordering is Py Pz Px for Libcint, which is the same as Gbasis (S1 C0 C1).

[PaulWAyers]

Are you sure that this is the ordering in libcint? Superficially it seems to disagree with the documentation
https://github.com/sunqm/libcint/blob/3c242cc8ea9d8e87d3f6a685aeb8d2ad897c74fc/doc/program_ref.txt#L174

[msricher]

I compiled it the with flag to set PyPzPx.

[PaulWAyers]

In this case, I think we should just pass on using libcint for this. It's not worth the time, and it's not the most common integral. We have an issue for it, and we can craft an issue for libcint to see if we can resolve it there.