Failing to generate .cube files

Dear Community,
I am new to Psi4 and want to do optimised geometry ground state HOMO-LUMO calculations and visualise them.

When I try to generate a cube file, the script runs without any errors but does not generate a cube file. When generating a .molden file, it works great (although I don’t know how to tell the optimiser to only calculate the frontier orbitals+/-n, such as: “set cubeprop_tasks [‘frontier_orbitals’]”).
Adding to it, I need to to operations on the orbitals, such as to convolve it with a kernel and I think a .cube file is the best for that and I cannot import .molden files into MatLAB or Mathematica.

Searching through the documentation and examples, I think this should be everything needed to generate these files:

memory 6 GB
set basis 6-31G

molecule {pubchem:benzene}

set freeze_core True
set cubeprop_tasks [‘frontier_orbitals’]
E, wfn = optimize(‘b3lyp’, return_wfn=True)
cubeprop(wfn)

Thank you for your help!

P.S. I did not know how to categorise this. If you have suggestions please let me know and I will change it.

Are you using the latest version, 1.5?

When I run your output I get 2 cube files.

Perhaps. Keep in mind that cube files are not necessarily the best representation of the orbitals. Good enough to visualize. Especially if you need to represent core electrons the typically coarse cubic grids are often insufficient. You should integrate the cube file and see if you the electron count is correct.
Depends what you want to do, of course.

You can get the MO orbital coefficient as numpy matrix as well.

Hey hokru,
thanks for the quick reply. When I run the file as posted in the latest psi4 1.5 in the anaconda psi4-conda environment using psi4 input.in -n 4 output.out. I actually get the same output. I shortened the script to feature pubchem:benzene to keep it short instead of my original input (as I figured it would not make a difference):

molecule {
0 1
   C       -6.39157        4.17697       -0.00000
   C       -5.15443        3.53275       -0.00000
   C       -5.09377        2.13924       -0.00000
   C       -6.27026        1.38996        0.00000
   C       -7.50740        2.03418        0.00000
   C       -7.56806        3.42768        0.00000
   H       -8.53193        3.92961        0.00000
   H       -6.43884        5.26267       -0.00000
   H       -4.23782        4.11653       -0.00000
   H       -4.12990        1.63732       -0.00000
   H       -6.22300        0.30425        0.00000
   H       -8.42401        1.45040        0.00000
}

And when I run this, I only the .out and geom.xyz files.

Changing the molecule gives me Psi_a_22_4-A"_LUMO.cube and Psi_a_21_3-A"_HOMO.cube

Alright, that is very strange. I made sure that I have all updates installed in anaconda. Do you have any idea what to do? For me it doesn’t make sense that it runs on the pubchem molecule but not the other one defined in Avogadro and exported using it’s extension.

Yes, this is strange. Could you upload your input and output files?

Yes, of course. I noticed that at the end of the working.out there is a success message, which is absent on the issue.out.
issue_in.txt (959 Bytes)
issue_out.txt (92.7 KB)

working_in.txt (364 Bytes)
workingout.txt (92.1 KB)

Hey Hokru,
is there maybe a specific format that the pubchem: command uses that I can mimick while specifying. Or can I use cdxml files locally to define something similar as pubchem?