I recently got interested in printing orbitals and densities as cube files by using the new version of the $plots section, which is invoked with MAKE_CUBE_FILES = TRUE and PLOTS = TRUE.
Given the following input:
$molecule
0 1
O -2.1904701 0.5202006 0.0237374
H -1.2011437 0.4799213 0.0090967
H -2.4798356 -0.3698583 -0.2997841
$end
$rem
BASIS = 6-31G*
MAKE_CUBE_FILES = TRUE
PLOTS = TRUE
!STATE_ANALYSIS = TRUE
METHOD = PBE
CIS_N_ROOTS = 2
CIS_TRIPLETS = 0
CIS_SINGLETS = 1
NTO_PAIRS = 4
$end
$plots
grid_range (-8,8) (-8,8) (-8,8)
grid_points 150 150 150
total_density 0-2
transition_density 1-2
natural_transition_orbital 1-2
$end
There are 3 issues I would like to discuss when running it with Q-Chem:
- A minor thing, I think there is a typo in the Q-Chem manual since the “natural_transition_orbitals” keyword in the “$plot” section is not recognized while “natural_transition_orbital” is.
- If I run the input with STATE_ANALYSIS = TRUE (commented in the example above) the job crashes with an error involving the FM_READ function. Is this keyword incompatible with printing cube files as described before?
- If NTOs are successfully printed. How do I understand which pair corresponds to the highest singular value? Files are named: NTO.TDA.state1.2.cube, NTO.TDA.state1.3.cube, … , up to NTO.TDA.state1.9.cube but this doesn’t help with the ordering criteria.
Thank you for the questions, and for your patience!
- You’re correct that this is a typo, and I’ll be sure to correct this in the manual. Thanks for pointing it out!
- After a bit of experimentation, the issue seems to be that STATE_ANALYSIS can’t read the $plots section in the new format, as you suggested. I’ve submitted a ticket about this in our system. In the meantime, I would suggest splitting the job into two separate jobs, when possible, using the old $plots formatting for the STATE_ANALYSIS job.
- I believe the ordering is the same as described in the “Direct Generation of Cube Files” section (see 11.5.4 Direct Generation of “Cube” Files‣ 11.5 Visualizing and Plotting Orbitals, Densities, and Other Volumetric Data ‣ Chapter 11 Molecular Properties and Analysis ‣ Q-Chem 5.2 User’s Manual). For a system with N_a alpha-spin MOs, the occupied NTOs (1, 2, 3… N_a) are stored in order of increasing amplitude, so the last one (N_a) is that with the highest amplitude. Virtual NTOs are stored next, in order of decreasing importance. Based on this, the NTO.TDA.state1.2.cube should be the occupied NTO with the the lowest amplitude, while NTO.TDA.state1.9.cube is the virtual with the lowest amplitude. NTO.TDA.state1.N_a and NTO.TDA.state1.(N_a+1) should represent the most important pair (I believe this should be 5 and 6 in this case).
Please let me know if there’s anything here that seems incorrect or needs clarification.
1 Like
Regarding point (3): what Shannon says sounds right to me. Note that you can always check this with IQmol, which is set up to visualize NTOs directly from a TDDFT calculation. If the look the same as your cube files, then you’ve figured out the order.
1 Like