Get water orbitals

Hi all,

I have problems understanding how to get the orbitals of the water molecule.
I have the water molecule:
h2o = psi4.geometry("""
H 1 0.96
H 1 0.96 2 104.5
and compute the energy:
eng, wfn =“scf/sto-3g”, molecule=h2o, return_wfn=True)

I want to get the orbitals of all the 10 electrons. Since for water, the alpha ad beta orbitals are equal
and I use the sto-3g basis set (2 x s + px + py + pz for O and 1 x s for each H).
I expected to obtain a 5x7 matrix. Now, when I run
I obtain a 4x4, 0x0, 1x1 and 2x2 matrix. How do I get the orbitals of the corresponding electrons?

Thanks in advance,

PSI4 is making use of point group symmetry of water, i.e C2V which has 4 irreps, A1, A2, B1, B2. The following is the occupancy of different irreps: A1 - 4, A2 - 0, B1 - 1, B2 - 2
Since the wavefunction transforms as fully symmetric irrep, i.e A1, you will get a 4x4, 0x0, 1x1 and 2x2 matrix, corresponding to A1 * A1, A2 * A2, B1 * B1, B2 * B2 combinations.
You can get the 5x7 matrix if you set the symmetry as c1 in your input file:

h2o = psi4.geometry(""“
H 1 0.96
H 1 0.96 2 104.5
symmetry c1
and then use
Cocc = wfn.Ca_subset("AO", "OCC")

Thanks for the quick answer! Just two followup questions:
1.Does the line
Cocc = wfn.Ca_subset(“AO”, “OCC”)
work for all molecules?
2. How to I assign the electrons to the atoms in the psi.geometry structure?

Yes, since Ca_subset(“AO”, “OCC”) is associated with the wavefunction.

Maybe, I didn’t understand this question properly, but you can specify the total charge and spin multiplicity in psi4.geometry.

h2o = psi4.geometry("""
0 1 # total charge and spin multiplicity
H 1 0.96
H 1 0.96 2 104.5
symmetry c1

Also, you could fragment your system and then can put charges on individual fragments as well. Let me know if this answers your question.

I meant:
How do I get the basis functions (along with the centers) corresponding to the rows of Cocc?

sorry for the late reply. I think you can read this:
Some examples:


rhf_e, rhf_wfn = energy('SCF', return_wfn=True)
basis = rhf_wfn.basisset()
# Returns the atomic center for the 0th shell
a = basis.shell_to_center(0)
# Return the number of shells on atom center 0
b = basis.nshell_on_center(0)
# Returns the first basis function of shell 0
c = basis.shell_to_basis_function(0)

Also, the atomic centers should be ordered the same as in input. hope this helps.