gbasis icon indicating copy to clipboard operation
gbasis copied to clipboard
verified publisher icon theochem

Evauation of angular momentum in GBasis vs. Libcint bindings

Open msricher opened this issue 2 years ago • 8 comments

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 https://github.com/theochem/gbasis/issues/137#issuecomment-1885115839

angularmomentum

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 avatar Jan 10 '24 19:01 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

msricher avatar Jan 18 '24 15:01 msricher

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.

Ali-Tehrani avatar Jan 18 '24 16:01 Ali-Tehrani

@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?

msricher avatar Jan 19 '24 22:01 msricher

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?

PaulWAyers avatar Jan 20 '24 18:01 PaulWAyers

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).

msricher avatar Jan 20 '24 19:01 msricher

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

PaulWAyers avatar Jan 20 '24 21:01 PaulWAyers

I compiled it the with flag to set PyPzPx.

msricher avatar Jan 20 '24 23:01 msricher

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.

PaulWAyers avatar Jan 21 '24 02:01 PaulWAyers