# How to obtain the solid formation of atomic orbitals

Dear Developers,
Is there a way to obtain the solid formation of atomic orbitals used in a calculation? For example, following the H2O example in Psi4,

import psi4

h2o = psi4.geometry("""
O
H 1 0.96
H 1 0.96 2 104.5
symmetry c1
units angstrom
""")
psi4.set_options({'reference': 'rhf', 'basis': 'def2-SVP'})
scf_e, scf_wfn = psi4.energy('scf', molecule=h2o, return_wfn=True)


scf_wfn.get_basisset('ORBITAL').print_detail_out() will return the following information,

  -AO BASIS SET INFORMATION:
Name                   = DEF2-SVP
Blend                  = DEF2-SVP
Total number of shells = 12
Number of primitives   = 22
Number of AO           = 25
Number of SO           = 24
Maximum AM             = 2
Spherical Harmonics    = TRUE

-Contraction Scheme:
Atom   Type   All Primitives // Shells:
------ ------ --------------------------
1     O     7s 4p 1d // 3s 2p 1d
2     H     4s 1p // 2s 1p
3     H     4s 1p // 2s 1p

==> AO Basis Functions <==

[ DEF2-SVP ]
spherical
****
O   1
S   5 1.00
2266.17677850          -0.00534318
340.87010191          -0.03989004
77.36313517          -0.17853912
21.47964494          -0.46427685
6.65894331          -0.44309745
S   1 1.00
0.80975976           1.00000000
S   1 1.00
0.25530772           1.00000000
P   3 1.00
17.72150432           0.04339457
3.86355054           0.23094121
1.04809209           0.51375311
P   1 1.00
0.27641544           1.00000000
D   1 1.00
1.20000000           1.00000000
****
H   2
S   3 1.00
13.01070100           0.01968216
1.96225720           0.13796524
0.44453796           0.47831935
S   1 1.00
0.12194962           1.00000000
P   1 1.00
0.80000000           1.00000000
****
H   3
S   3 1.00
13.01070100           0.01968216
1.96225720           0.13796524
0.44453796           0.47831935
S   1 1.00
0.12194962           1.00000000
P   1 1.00
0.80000000           1.00000000
****


It is obvious that 25 AOs were used, while scf_wfn.basisset().compute_phi(0, 0, 0).shape will return (24,), which seems to be the number of SO. Is there a way to obtain the details of atomic orbitals?

I don’t understand what “the solid formation” is. Are you asking for the change-of-basis matrix between AO and SO?

My naive guess is the “discrepancy” is due to Cartesian vs spherical.

Thanks for your reply. I want to obtain an array of AOs when a position (x,y,z) is given, just like what the function scf_wfn.basisset().compute_phi does. Since the number of AOs is 25, the function should return a (25, ) array rather than a (24, ) array.

This depends on how you define an AO. Does the d shell have 5 AOs (including x^2-y^2, which is the more physically sensible one), or does it have 6 AOs (including x^2 and y^2, which is less physically sensible)?

When the printout says “Number of AO = 25”, it is using the version where d shells have 6 AOs. Is this what you want?

What assumptions are you making about the normalization of each AO?

Yes, you’re right. However, I think which AOs of d shell are selected by the psi4.set_options('basis': 'def2-SVP'). If I chose another basis, the number of AOs will also change. The form of AOs is determined by the selection of basis and I want to know if psi4 could outpout this function.

I do not understand what the “form” is. Can you give me an example?

Psi4 always includes either all AOs of a given shell on a given atom or none. wfn.get_basisset('ORBITAL').print_detail_out() will always determine “how many AOs are in this shell?” assuming the Cartesian convention. Therefore, the only reason why that command gives different numbers of AOs for different basis sets on the same molecule is because different basis sets include different shells. So for example, the cc-pVDZ basis set may include the 1s, 2s, and 2p shells for hydrogen atoms, where cc-pVTZ includes the 3s, 3p, and 2d shells, leading to 1+3+6=10 additional AOs per hydrogen atom.

I also don’t understand how this relates to your original question about 25 AOs vs 24 SOs.

For example, maybe the Cartesian Gaussian-type functions is used as AOs, then the form is $\chi_i=N(x-X)^k(y-Y)^l(z-Z)^m\exp{-\zeta_i(r-R)^2}$. I’m not sure whether the form is this and how the coefficients is determined. In principle, the coefficents should be found in the basis settings in psi4 but I don’t know the format of the basis set files.
And the question w.r.t 25 AOs and 24 SOs is the shape of matrix obtained from psi4, such as density matrix, which is represented under the AO basis, is (24, 24). However, scf_wfn.get_basisset('ORBITAL').print_detail_out() tells us that the number of AOs is 25, which is not 24. I understand your statements that d shells have 6 AOs, because for the form of AOs above, k+l+m=2 will return 6 results. Then how the shape of density matrix represented in AO basis is (24, 24)?

I can see how it can be confusing since:

-AO BASIS SET INFORMATION:
Name                   = DEF2-SVP
Blend                  = DEF2-SVP
Total number of shells = 12
Number of primitives   = 22
Number of AO           = 25   # <-- always cartesian, even if spherical harmonics =true
Number of SO           = 24
Maximum AM             = 2
Spherical Harmonics    = TRUE


You have 24 AOs using spherical harmonics which can be seen in the basis set output and which gives the shape of the compute_phi matrix.

  ==> Primary Basis <==

Basis Set: DEF2-SVP
Blend: DEF2-SVP
Number of shells: 12
Number of basis functions: 24
Number of Cartesian functions: 25
Spherical Harmonics?: true
Max angular momentum: 2

1 Like

Yes. So the question is how to determine the coefficients of the spherical harmonics for atom A and orbital $\mu$ in psi4.

Yeah the number of AOs is definitely reported incorrectly…

If I’m understanding what you’re looking for, you can go through each shell in the basisset. For example:

basis this_basis {
assign    def2-SVPD
}

molecule {
0 1
C        0.60720        0.00000       -0.00040
O       -0.60040        0.00000        0.00010
H        1.14720        0.93530        0.00160
H        1.14720       -0.93530        0.00160
}

nrg, wfn = energy('M06-2X', return_wfn=True)

basis_set = wfn.basisset()
for i in range(0, basis_set.nshell()):
shell = basis_set.shell(i)
print("shell number:", i + 1)
print("functions", shell.nfunction)
print("angular momentum:", shell.amchar)
for j in range(0, shell.nprimitive):
print("primitive:", j + 1)
print("exponent", shell.exp(j))
print("coefficient", shell.coef(j))
print("\n")


This goes through each shell based on the order, but the documentation says you can also get a shell from a center: BasisSet

The GaussianShell documentation is here: GaussianShell. There’s more data available than what my example will print, so check that out.

If you are interested in the format of basis set files, Psi4 uses the Gaussian basis file format (Basis Sets). This format is described in Gaussian’s documentation (see the Input tab): Gen and GenECP | Gaussian.com

Thank you vey much. Another related question, I found a link. It seems that psi4 uses gau2grid to return the value of the spherical AOs. So I can invoke the gau2grid.collocation function to obtain what I want. Is there a misunderstanding?

I still cannot identify what you want.

So basically, you want a function where the input is a point in space, and the output is the function value of each individual Gaussian-type atomic orbital at that point? And you want this to work for pure AM basis sets where applicable?

I’m not exactly knowledgeable in the programming side of Psi4, so I’m afraid I can’t be too helpful here. I can’t tell what that collocation function is doing by looking at it.

If you wanted to make your own function for this, you would need to know the functional form that Psi4 uses (e.g., does Psi4 use the standard f orbital functions, or the cubic f functions, and what exactly the standard/cubic functions are).

As @ajs99778 said, I need the function whose input is a point in space, and output is the function value of each individual Gaussian-type atomic orbital at that point, but I don’t know what the functional form that Psi4 uses.

Okay, that helps. Two more technical clarifying questions:

1. Is an “individual” Gaussian-type orbital a linear combination of Gaussians used to form an AO, or a single Gaussian even though many of them are needed to form an AO?
2. When you say “AO,” are you talking about a Cartesian or a spherical? Again, Cartesians use x/y/z, and there are 6 d Cartesians, which is unphysical. Spherical are what you’d see for the exact eigenstates of atomic hydrogen, and there are 5 d sphericals. You can express either as a linear combination of the other. “Whichever one Psi actually uses when doing its SCF for def2-SVP” is a perfectly fine answer to this question.

For the first question, I think the Gaussian-type orbital should be the linear combination of Gaussian used to form an AO. For the second question, I think it should be the spherical AO, since the shape of Hamiltonian or the density matrix obtained from Psi4 is the (number of SO, number of SO).

Thanks for your patience - I’ve been getting a lot of things back in order the last few days, so I’ve let the Psi forum slip.

It sounds like you want the compute_phi function but for a single AO rather than summed over all AOs. Is that about right?

Yep, that’s what I want.

Got it, thanks.

I’ve put up an issue to get this request added to our list. This isn’t something that Psi currently supports. If you wanted to add it yourself (the joy of Psi being an open-source code), I could give you some assistance.