Sequential calculations in single python function

I’m trying to do a geometry optimization with b3lyp/6-31g(d,p), and then get the HF orbitals at that geometry with, say, the cc-pvdz basis. My attempt is given below, but I think it’s simply assigning the KS orbitals to the wfn object. Is there something besides clean() that I need to do here?

    psi4.core.set_output_file('output.dat', False)
    psi4.set_options({'scf_type': 'pk', 'd_convergence': 1e-12})
    psi4.core.set_output_file('output.dat', False)
    psi4.set_options({'basis': str(basis), 'scf_type': 'pk', 'reference': 'rhf', 'd_convergence': 1e-12})
    if molecule.multiplicity != 1:
        psi4.set_options({'reference': 'uhf'})
    self.hf_energy, wfn ='scf', return_wfn=True)

When I try modifying your script to get something that runs, wfn.Ca() returns orbitals in the new basis, not the optimization basis. Please post a minimal working example, and a demonstration that the orbitals you get are incorrect.

Thanks for taking a look; I was unable to reproduce the error in a minimal example, so it must be a bug in my code. Isolated working example for completeness:

import psi4
psi4.set_options({‘scf_type’: ‘pk’, ‘reference’: ‘rhf’, ‘d_convergence’: 1e-12})
geometry = “”"
0 1
H 0 0 -1.249957094407
Cl 0 0 0.036024574790
molecule = psi4.geometry(geometry)
psi4.set_options({‘scf_type’: ‘pk’, ‘reference’: ‘rhf’, ‘basis’: ‘sto-3g’, ‘d_convergence’: 1e-12})
hf_energy, wfn =‘scf’, return_wfn = True)

with Ca() printout:

(array([[ 5.42518555e-04, 1.06618800e-02, 1.44155335e-02,
1.84047944e-01, 5.01509845e-01, 1.15140659e+00],
[ 9.94504338e-01, 3.77058344e-01, -5.01841637e-03,
8.23316225e-02, -3.11008416e-02, -2.52963093e-02],
[ 1.57387517e-02, -1.05359763e+00, 1.42078692e-02,
-2.73128059e-01, 9.51501779e-02, 6.64916480e-02],
[-1.24231787e-04, 1.13384921e-02, 9.88437588e-01,
4.26989079e-02, 2.00777618e-01, -2.20267433e-01],
[-1.84369795e-03, -4.51056306e-02, -6.80266191e-03,
9.14080471e-01, -4.70595132e-01, -5.37773644e-01],
[ 2.64843001e-04, 6.19057862e-03, 4.73626720e-02,
-1.29484589e-01, -6.76314203e-01, 9.62101210e-01]]), array([], shape=(0, 0), dtype=float64), array([[ 0.99024376, 0.27677981],
[ 0.03887301, -1.02746224]]), array([[ 0.99024376, 0.27677981],
[ 0.03887301, -1.02746224]]))

Good luck with the debugging!

It sounds like that’s all you needed from Psi developers, so I’m going to mark this topic solved. Post back if there’s still a question.

Just confirmed this is not a psi4 issue. What I wasn’t capturing in my working example was that I was attempting to transform both the alpha and beta fock matrices, but psi4 has these referencing the same matrix for rhf. So I was trying to run code for MO’s on untransformed AO’s.