HOMO-LUMO energies

kpinto · April 18, 2017, 10:22am

Hello,

Just a silly question. Anyone can help me to retrieve HOMO and LUMO energies? I cannot see which keyword I should use to do that. Any template. I would like to use also semi empirical methods.

Thanks in advance,

Best Warmest Regards,

Kevin

dgasmith · April 18, 2017, 11:26pm

We do not currently have semiempirical methods in Psi4. The HOMO-LUMO is not explicitly printed, but part of the eigenvalue printing at the end of a SCF program. Did you want to obtain these programmatically?

kpinto · April 20, 2017, 1:30pm

Hello again,

Thanks for your answer. What I would like to extract is the average energy of HOMO and LUMO for example using semiempirical methods, and then calculate HOMO-LUMO difference. Due to the fact that Psi4 does not have semiempirical, DFT would be perfect. Any suggestion of how to get that value?

Regads,

Kevin

dgasmith · April 20, 2017, 2:20pm

Sure ,you can extract this data from a Wavefunction object like so:

molecule {
He
}

scf_e, scf_wfn = energy("B3LYP/cc-pVDZ", return_wfn=True)

HOMO = scf_wfn.epsilon_a_subset("AO", "ALL")[scf_wfn.nalpha()]
LUMO = scf_wfn.epsilon_a_subset("AO", "ALL")[scf_wfn.nalpha() + 1]

print("\nThe HOMO - LUMO gap is %16.8f hartree" % (LUMO - HOMO))

kpinto · April 21, 2017, 10:43am

Thank you so much for your answer.

Would you mind telling me which energy units are used?

Best Regards,

Kevin

dgasmith · April 21, 2017, 12:11pm

It will be in hartrees.

kpinto · April 24, 2017, 6:45am

Thank you so much for answering me.

I got the following error when running the program for this molecule (COC(=O)C@@HCCSC):

To start psi4 job, run: psi4 psi4-mol00000004.opt -n 2

The HOMO - LUMO gap is 0.02278216 hartree
1 molecule converted
24 audit log messages
To start psi4 job, run: psi4 psi4-mol00000000.opt -n 2
An error has occurred python-side. Traceback (most recent call last):

File “”, line 54, in

File “/opt/psi4conda/share/psi4/python/driver.py”, line 405, in energy
return driver_cbs._cbs_gufunc(energy, name, ptype=‘energy’, **kwargs)

File “/opt/psi4conda/share/psi4/python/driver_cbs.py”, line 1557, in _cbs_gufunc
ptype_value, wfn = func(method_name, return_wfn=True, molecule=molecule, **kwargs)

File “/opt/psi4conda/share/psi4/python/driver.py”, line 446, in energy
wfn = procedures[‘energy’][lowername](lowername, molecule=molecule, **kwargs)

File “/opt/psi4conda/share/psi4/python/procedures/proc.py”, line 2446, in run_dft
scf_wfn = run_scf(name, **kwargs)

File “/opt/psi4conda/share/psi4/python/procedures/proc.py”, line 1689, in run_scf
scf_wfn = scf_helper(name, **kwargs)

File “/opt/psi4conda/share/psi4/python/procedures/proc.py”, line 1170, in scf_helper
scf_wfn = psi4.scf(ref_wfn, precallback, postcallback)

RuntimeError:
Fatal Error: RHF: RHF reference is only for singlets.
Error occurred in file: /scratch/cdsgroup/conda-builds/work/src/lib/libscf_solver/rhf.cc on line: 71
The most recent 5 function calls were:

psi::PsiException::PsiException(std::string, char const*, int)
psi::scf::RHF::common_init()
psi::scf::RHF::RHF(boost::shared_ptrpsi::Wavefunction, psi::Options&, boost::shared_ptrpsi::PSIO)
psi::scf::RKS::RKS(boost::shared_ptrpsi::Wavefunction, psi::Options&, boost::shared_ptrpsi::PSIO)
psi::scf::scf(boost::shared_ptrpsi::Wavefunction, psi::Options&, _object*, _object*)

Which keyword should I use to run UHF for that molecules ( radicals, …) so that I can obtain HOMO and LUMO energies?

Best Warmest Regards,

Kevin

dgasmith · April 24, 2017, 6:56pm

You need to place set reference uhf in the options to run the computation. Not too familiar with HOMO/LUMO gaps for open shell systems. If you want to the HOMO/LUMO in terms of a singlet excitation, the previous script will work fine. If you could describe exactly what you want it would be helpful.

kpinto · April 25, 2017, 7:14am

I would like to calculate several molecular descriptors for building a photolysis model. Then I need to calculate HOMO/LUMO gap to obtain the absorbance for different kind of molecules (radicals, neutral, …). I tried to set up reference as uhf and it worked perfectly.

Best Regards,

Kevin

kpinto · April 25, 2017, 8:18am

Hi again,

I tried with another molecule and it showed this error:

Failed to converge.
Iterations did not converge.An error has occurred python-side. Traceback (most recent call last):

File “”, line 55, in

File “/opt/psi4conda/share/psi4/python/driver.py”, line 405, in energy
return driver_cbs._cbs_gufunc(energy, name, ptype=‘energy’, **kwargs)

File “/opt/psi4conda/share/psi4/python/driver_cbs.py”, line 1557, in _cbs_gufunc
ptype_value, wfn = func(method_name, return_wfn=True, molecule=molecule, **kwargs)

File “/opt/psi4conda/share/psi4/python/driver.py”, line 446, in energy
wfn = procedures[‘energy’][lowername](lowername, molecule=molecule, **kwargs)

File “/opt/psi4conda/share/psi4/python/procedures/proc.py”, line 2446, in run_dft
scf_wfn = run_scf(name, **kwargs)

File “/opt/psi4conda/share/psi4/python/procedures/proc.py”, line 1689, in run_scf
scf_wfn = scf_helper(name, **kwargs)

File “/opt/psi4conda/share/psi4/python/procedures/proc.py”, line 1170, in scf_helper
scf_wfn = psi4.scf(ref_wfn, precallback, postcallback)

RuntimeError:
Fatal Error: Iterations did not converge.
Error occurred in file: /scratch/cdsgroup/conda-builds/work/src/lib/libparallel/process.cc on line: 207
The most recent 5 function calls were:

psi::PsiException::PsiException(std::string, char const*, int)
psi::die_if_not_converged()
psi::scf::HF::finalize_E()
psi::scf::HF::compute_energy()
psi::scf::scf(boost::shared_ptrpsi::Wavefunction, psi::Options&, _object*, _object*)

dgasmith · April 25, 2017, 7:02pm

Sure, not all molecules converge well (especially radicals with DFT methods) without some work. You can try posting your molecule/method to receive feedback on convergence manipulation.

kpinto · April 26, 2017, 6:42am

calculate homo - lumo gap energy with SCF in hartree

memory 8 GB # use as much as you’ve got

For Cartesian coordinates, use XYZ coordinates

COORDINATES will be replaced with XYZ or ZMatric coordinates

molecule mol00000005 {
C -2.957900 -3.165400 0.000000
C -3.404400 -1.733400 0.000000
C -4.781900 -2.327100 0.000000
S -1.933600 -1.439000 0.000000
C -1.759000 0.050800 0.000000
C -1.132600 1.413700 0.000000
N 0.274100 1.934500 0.000000
C 1.428500 0.976700 0.000000
O 1.176200 -0.502000 0.000000
C 2.835200 1.497500 0.000000
N 3.087500 2.976100 0.000000
C 3.989600 0.539600 0.000000
C 3.737200 -0.939000 0.000000
C 4.891600 -1.896800 0.000000
C 6.298300 -1.376000 0.000000
O 7.452700 -2.333800 0.000000
C 6.550600 0.102600 0.000000
C 5.396200 1.060400 0.000000
C -2.495500 2.040100 0.000000
O -3.016300 3.446800 0.000000
N -3.121900 0.677200 0.000000
C -4.138800 -0.425500 0.000000
C -5.628700 -0.250900 0.000000
O -6.222400 1.126600 0.000000
O -6.524800 -1.453800 0.000000
}

set reference uhf
scf_e, scf_wfn = energy(“B3LYP/cc-pVDZ”, return_wfn=True)

HOMO = scf_wfn.epsilon_a_subset(“AO”, “ALL”)[scf_wfn.nalpha()]
LUMO = scf_wfn.epsilon_a_subset(“AO”, “ALL”)[scf_wfn.nalpha() + 1]

HL_nrg = HOMO - LUMO

print("\nThe HOMO - LUMO gap is %16.8f hartree" % (LUMO - HOMO))

XYZ_FILE_TO_SAVE and SD_FILE_TO_SAVE will be replaced with filenames

mol00000005.save_xyz_file(“psi4-mol00000005.xyz”, False)

command="babel -ixyz psi4-mol00000005.xyz -osd psi4-mol00000005.sdf --property ‘HL_nrg(hartree)’ %.4f " % HL_nrg
os.system(command)
command=‘echo "Psi4 results for psi4-mol00000005.sdf " | mail -s “Psi4 results” -a psi4-mol00000005.sdf kpinto.gil@gmail.com’
os.system(command)

loriab · April 26, 2017, 5:40pm

You’ll want to look at the output file and see how the SCF iterations are (not) converging – try grepping @ in the output and posting it. Haven’t run your molecule, but some things to try altering are set guess gwh or switching from b3lyp to hf and trying set soscf true.

cisert · September 24, 2020, 10:34am

Dear all, I’m wondering about the indexing into the numpy array in the line

HOMO = scf_wfn.epsilon_a_subset("AO", "ALL")[scf_wfn.nalpha()]

To my understanding, scf_wfn.nalpha() provides the number of occupied orbitals. Indexing into a numpy array should be done with scf_wfn.nalpha()-1 in order to retrieve the nalpha-th element. E.g., for water:

print(scf_wfn.epsilon_a_subset('AO', 'ALL').np)
[-19.273, -1.1025, -0.628, -0.479, -0.408, 0.075, 0.134 ...] # rounded for simplicity

print(scf_wfn.nalpha())
5

print(scf_wfn.epsilon_a_subset('AO', 'ALL').np[scf_wfn.nalpha()])
0.075 # instead of the expected 5-th element of the array (-0.408)

Am I missing something or should HOMO and LUMO rather be accessed via

HOMO = scf_wfn.epsilon_a_subset("AO", "ALL")[scf_wfn.nalpha()-1]
LUMO = scf_wfn.epsilon_a_subset("AO", "ALL")[scf_wfn.nalpha()]

Thanks for your help!
Best, cisert

jmisiewicz · September 24, 2020, 12:06pm

You’re correct. The HOMO is scf_wfn.nalpha()-1 and the LUMO is scf_wfn.nalpha() for a closed-shell system. For an open-shell system, determining what you mean by HOMO and LUMO is more complicated.

As this is a sufficiently long topic that the context of previous posts hurts rather than helps, I’m going to lock it. Let’s put further discussion in a new topic.