Molecular Orbitals & ESP

[sbembenek]

Is there a way to get both the Molecular Orbitals & ESP in a single run?

After optimizing my system, I usually do separate runs for each of these calculations starting from the optimized coordinates. However, since I still need to go through an iteration to get the needed energy and wave function convergence, I’d like to just do this step once and get both the MOs and ESP, if possible.

thanks!

[jmisiewicz]

What commands are you using to get the MOs and ESP currently?

Also, I’ve removed the “geometry optimization” category because that’s not relevant to your question.

[sbembenek]

Thanks for moving it.

For ESP:
e, wfn = energy(‘b3lyp-d3bj/aug-cc-pVTZ’, return_wfn=True)
set cubeprop_tasks [‘esp’]
cubeprop(wfn)

For MOs;
e, wfn = energy(‘b3lyp-d3bj/aug-cc-pVTZ’, return_wfn=True)
cubeprop(wfn)

[jmisiewicz]

See the cubeprop section of the manual. Adding multiple tasks to cubeprop is supported, and having orbitals as a cubeprop task is also supported. If that doesn’t answer your question, I’m going to need more detail.

For future reference, please include a minimal input file whenever you have a question about particular input. When you said “get the molecular orbitals,” I thought you meant the Ca coefficients and was very confused.

[sbembenek]

Thanks for the info!

With the goal being to get the MOs, ESP, and Properties (“DIPOLE_POLARIZABILITIES”, “MULLIKEN_CHARGES”, “DIPOLE”) without have to set up the density each time, would something like this work then:

molecule mol {
0 1
O          0.44260       -0.88720       -0.29490
H          0.73940       -0.26640       -0.96690
H          0.79240       -0.54630        0.53350

no_reorient
no_com
}

set scf_type df
set basis aug-cc-pVTZ
set reference rks
set s_tolerance 1e-9
set_num_threads(64)

e, wfn = energy('b3lyp-d3bj/aug-cc-pVTZ', return_wfn=True)
set cubeprop_tasks ['esp']
set cubeprop_tasks ['orbitals']
cubeprop(wfn)

psi4.properties("b3lyp-d3bj/aug-cc-pVTZ", properties=["DIPOLE_POLARIZABILITIES", "MULLIKEN_CHARGES", "DIPOLE"])
mol.print_out()
print_variables()

[jmisiewicz]

I don’t deal with those kinds of calculations much, but it looks plausible, with a few things to note:

1. Your second set cubeprop_tasks overrides the first. Use a single list, ['esp', 'orbitals']
2. Why do you need a separate energy and properties call? The properties function with return_wfn=True will return the energy and the wavefunction, just like the energy function does.

[sbembenek]

thanks!

I’m still learning much of the nomenclature for psi4. So this is why I have both a property and energy call – I’m not totally clear on all of it yet

Overall then, are you saying it would look like:

psi4.properties("b3lyp-d3bj/aug-cc-pVTZ","return_wfn=True", properties=["DIPOLE_POLARIZABILITIES", "MULLIKEN_CHARGES", "DIPOLE"])
set cubeprop_tasks ['esp','orbitals']
cubeprop(wfn)

Also, are the double quotes and single quotes in the above interchangeable?

thanks again!

[jmisiewicz]

Single quotes and double quotes are inter-changeable.

You’ll still need to set return_wfn=True and actually get the wfn from the properties call.

[sbembenek]

Do you mean (remove the quotes) like this:

psi4.properties("b3lyp-d3bj/aug-cc-pVTZ", return_wfn=True, properties=["DIPOLE_POLARIZABILITIES", "MULLIKEN_CHARGES", "DIPOLE"])

[jmisiewicz]

You also need to actually store the return values to a variable, but that’s the part about supplying the return_wfn keyword, yes.

If you haven’t seen it already, please read the tutorial.

[sbembenek]

Great. Pulling it altogether then, I used:

props, wfn = psi4.properties('b3lyp-d3bj/aug-cc-pVTZ', return_wfn=True, properties=["DIPOLE_POLARIZABILITIES", "MULLIKEN_CHARGES", "DIPOLE"])
set cubeprop_tasks ['esp','orbitals']
cubeprop(wfn)

Which works well.

Thanks again!