How does pyscf's internal basis format distinguish between 1s and 2s orbitals?

Question

For example: gto.basis.load("6-31g", "Li") will show:

[[0,
  [642.41892, 0.0021426],
  [96.798515, 0.0162089],
  [22.091121, 0.0773156],
  [6.2010703, 0.245786],
  [1.9351177, 0.470189],
  [0.6367358, 0.3454708]],
 [0, [2.3249184, -0.0350917], [0.6324306, -0.1912328], [0.0790534, 1.0839878]],
 [0, [0.035962, 1.0]],
 [1, [2.3249184, 0.0089415], [0.6324306, 0.1410095], [0.0790534, 0.9453637]],
 [1, [0.035962, 1.0]]]

This is a list with 5 elements (which are also lists). The first element of each of these 5 lists represents angular momentum. In this example, there are 3 lists starting with a 0; how does pyscf determine if the list of [exponent, coefficient] that follow are for the 1s or 2s orbitals?

Answer 1

There is no distinction of any basis functions in PySCF. Basis functions are just that: basis functions, a mathematical abstraction. The individual basis functions do not have to represent orbitals; all that matters is the span of the basis set. The 1s and 2s orbitals in atoms, for instance, will be variationally found by optimizing the self-consistent field energy in the defined basis set.

Answer 2

The difference between "S" and "P" is quite significant in that S-type basis functions have a different mathematical form compared to P-type functions (and therefore the shapes of these functions are fundamentally different), but the difference between the "1" and the "2" is not as significant: The 1 just denotes the first orbital of a particular type, whereas the 2 denotes the next one.

So which orbital is labeled with a 1 and which one is labeled with a 2? If you need to label the orbitals numerically like this, what would make the most sense is to label them in a way consistent with the way everyone learns about orbitals starting from high school and early-undergraduate-level chemistry classes: the 1s orbital is closest to the nucleus, the 2s is further out, etc.

Usually basis sets are presented with the S-type orbital values appearing in order, so the first number you see is for the 1s, the second is for the 2s, etc. This happens to be true for the basis set in your question. But in case you doubted that the basis set was presented "in order" like this, you can look at the size of the exponents. 642.41892 is much larger than 96.798515, so if you were to plot the corresponding two basis functions, you'd find the first one to resemble the "1s" vs the one with the smaller number representing the "2s" orbital.

For the record I'll put the basis set from your question in ACESII/CFOUR/MRCC format below, since it's less confusing (at least to me), and also because the numbers in your version seem to be truncated, which would cause you to get energies inconsistent with other 6-31G calculations:

    LI:6-31G
    6-31G valence double-zeta
    
      2
        0    1
        3    2
       10    4
    
    0.6424189150D+03 0.9679851530D+02 0.2209112120D+02 0.6201070250D+01 0.2324918408D+01 
    0.1935117680D+01 0.6367357890D+00 0.6324303556D+00 0.7905343475D-01 0.3596197175D-01 
    
    0.2142607810D-02 0.00000000       0.00000000 
    0.1620887150D-01 0.00000000       0.00000000 
    0.7731557250D-01 0.00000000       0.00000000 
    0.2457860520D+00 0.00000000       0.00000000 
    0.00000000      -0.3509174574D-01 0.00000000 
    0.4701890040D+00 0.00000000       0.00000000 
    0.3454708450D+00 0.00000000       0.00000000 
    0.00000000      -0.1912328431D+00 0.00000000 
    0.00000000       0.1083987795D+01 0.00000000 
    0.00000000       0.00000000       0.1000000000D+01 
    
    0.2324918408D+01 0.6324303556D+00 0.7905343475D-01 0.3596197175D-01 
    
    0.8941508043D-02 0.00000000 
    0.1410094640D+00 0.00000000 
    0.9453636953D+00 0.00000000 
    0.00000000       0.1000000000D+01