Problem with generating cube density files

Hello everyone!

I’m trying to generate cube files for ground state and delta-SCF state densities ($rem and $plot sections of the jobs below). However, when I later subtract them using gabedit, the difference density is zero. What could be the reason for such behaviour?

$rem
jobtype sp
method camb3lyp
dft_d d3_bj
basis gen
basis2 sto-3g
purecart 11
print_general_basis true
gen_scfman false
max_scf_cycles 500
scf_convergence 9
scf_final_print true
mem_static 2000
mem_total 80000
solvent_method pcm
make_cube_files true
plots true
$end

$plots
grid_points 200 200 200
total_density 0
$end

@@@

$rem
jobtype sp
method camb3lyp
dft_d d3_bj
basis gen
purecart 11
print_general_basis true
gen_scfman false
max_scf_cycles 500
scf_guess read
scf_convergence 9
scf_final_print true
mem_static 2000
mem_total 80000
solvent_method pcm
molden_format true
state_analysis true
iqmol_fchk true
plots true
make_cube_files true
$end

$plots
grid_points 200 200 200
total_density 0
$end

I will be thankful for any explanations or tips.

I don’t use gabedit, usually just subtract with a command-line awk script. Have you looked at the cube files? They are human-readable ASCI. Are they identical? Maybe run a test with a smaller grid to make it easier to see what’s going on.

Looking at your input, however, I don’t understand why the 2nd job should produce anything different from the first job. What are you trying to do here?

Thank you for replying. As I mentioned, I am performing deltaSCF calculation - I change the occupation manually so that it corresponds to HOMO-LUMO exctitation. My bad that I didn’t include the $occupied section of the second job above, I’ll append it below.

The cube files are indeed identical, but I don’t know why that is - given that the orbitals for two distinct states that I’m working with are actually different (I checked that through fchk file).

$occupied
1:148
1:146
$end

Are the SCF energies the same? My guess is that the 2nd calculation undergoes variational collapse to the ground state. To avoid that, you need a SCF convergence algorithm that is designed to converge non-Aufbau solutions. Choices in Q-Chem are MOM, IMOM, STEP, or SGM. In Q-Chem manual, see here:
https://manual.q-chem.com/latest/sec_DeltaSCF.html
and in the literature, see here:
https://doi.org/10.1021/acs.jctc.0c00502

I recommend starting with IMOM (not regular MOM). If that fails, try STEP or SGM, which are more robust but a bit more expensive.

I don’t think that is the case here; the SCF energies are different (the difference is 3,24 eV).

It’s a bit hard to believe that this just works without some alternative convergence algorithm. What about the orbital energies? For successful delta-SCF, you should be able to see a virtual level that is below the HOMO, corresponding to the hole.

Well, this is a system with many excited states in close proximity to each other. Followng your advice I have run additional calculations with 1) IMOM, 2) STEP and 3) MOM. Using IMOM I get the same excitation energy as without any additional procedures, while with STEP and MOM I obtain a slightly lower lying state, with E(exc) = 3,23 eV.

As to the orbital energies, in all the cases the unoccupied 147 (hole) beta orbital is in fact higher than 148 alpha, so this part is fine.

However, I figured out my problem: the cube files weren’t generated for dSCF state at all, it worked only after I had removed “state_analysis true” from the second job.

That makes some sense in hindsight, sorry I didn’t catch that. The original cube file code (my code) is sort of orthogonal to the new libwfa code (STATE_ANALYSIS, Felix Plasser’s code).