[BUG] Inconsistent contractions between gbasis inner workings and iodata wrapper?

[marco-2023]

Describe the bug
It is not clear what are the objects returned by the from_iodata interface. If the returned lists contains the contracted shells, then the overlap and other integrals should return square matrix with the number of rows/columns
equal to the returned value from from_iodata.

To Reproduce

from iodata import load_one
from gbasis.wrappers import from_iodata
import numpy as np
from gbasis.integrals.overlap import overlap_integral

mol = load_one("./data/c2h4_q0.fchk")
ao_basis = from_iodata(mol)

print("Number of contracted basis functions: ", len(ao_basis))

# compute overlap of atomic orbitals (AO)
int1e_overlap_ao = overlap_integral(ao_basis)
print("Overlap matrix (S) of atomic orbitals (AO) shape: ", int1e_overlap_ao.shape)

returns:

Number of contracted basis functions:  64
Overlap matrix (S) of atomic orbitals (AO) shape:  (184, 184)

Expected behaviour
The overlap matrix should be of dimension NxN where N is the number of contracted shells that from_iodata returns

Screenshots

System Information:

• OS:
• Python version:
• NumPy version:
• SciPy version:

Additional context

[leila-pujal]

Hey @marco-2023, if I am not mistaken and understood your message properly, the length of generalized contractions shells is not equal to the number of AO. Each shell represents all the contractions that share the same set of exponents, centers, and angular momentum (AOs). For the overlap, this is taken care of in base_two_symm.py the specific construct_array function calls construct_array_contraction in overlap.py and it passes the angular momentum of the two pair contractions integrated.

gbasis/gbasis/integrals/overlap.py

Line 101 in bd7bd7d

contractions_one.angmom_components_cart,

This for a p orbital for example looks like:

[[1 0 0]
 [0 1 0]
 [0 0 1]]

And this done for all the angular momentums generates the (184,184). Let me know if this makes sense to you and explains your bug.

[leila-pujal]

If you try the code below you should get the 184 as total number of atomic orbitals:

total_ao = 0
for s in basis:
    if s.coord_type == 'cartesian':
        total_ao += s.angmom_components_cart.shape[0]
    elif s.coord_type == 'spherical':
        total_ao += len(s.angmom_components_sph)

print('TOTAL NUMBER AOs: ', total_ao)

[marco-2023]

Thank you very much @leila-pujal, your answer helps a lot. I was not clear on what the from_iodata function was returning.

[leila-pujal]

No worries, I know the shell contractions are confusing because it looks like they should relate to the number of AOs. Let me know if you have any more questions related to this. Feel free to close the issue if you think everything is clear.

[PaulWAyers]

@leila-pujal and @marco-2023 is this a case where there is a failure of documentation? We should probably be really clear about documentation, especially the dimension of objects. We probably also need to define, for those who are not among the QC-Devs cognoscenti, our verbiage (AOs = atomic basis functions; shells = basis functions with the same set of exponents; etc.) Having an example tutorial (or including it in an example tutorial) or something like that is helpful too.

In some sense, it may go into the scientific documentation, but we need to make it feasible/usable for people who haven't seen all the documentation too.

[leila-pujal]

hi, I agree we may have not been really conscious about our notation/dimensions. In my opinion, I think the part about covering it in the documentation is there. I think we are including this in the paper and the docstring for the contractions base class. I was trying to look at IOData earlier to see if they give more information about this I only found this. On the other hand, I agree that reading that a contractions groups together all AO for that angular momentum does not give you a clear idea of what the shells mean. I had to spend some time myself at the begining of using Gbasis. I agree this is one of the most important things to show in the tutorial. I am mentioning things from the email and issue #162 (sorry about mixing) about the units. if I am not mistaken (I will double check) Gbasis does not follow a particular convention meaning it just takes the convention of the parser or wrapper the user is using. Indeed this is missing in the documentation and we should state it clearly early in the documentation that the user should know what are the units is using. For example if using IOData it internally uses Bohr (I think). PySCF wrapper I checked and Gbasis just reads whatever is stored there so it could be Bohr or Angstrom. Anyway, I'll check the documentation and comment on issue #162 about where we could add this. @PaulWAyers I know you mentioned there to check point-sets passed in, basis functions, and the properties/integrals, I ll check those.