Plotting orbital densities¶
In this tutorial we will plot the HOMO density of the water molecule. We will do this in two steps: first we will run the SCF and save the DFCOEF file, then we will restart from this and calculate the orbital density on a grid of points and write this into a cube file which can be opened by your favorite molecular visualization program (the author of this tutorial likes Avogadro).
Getting the wave function¶
Let us use the following molecule input (h2o.mol):
DIRAC C 2 1. 2 H -1.4523499293 .0000000000 .8996235720 H 1.4523499293 .0000000000 .8996235720 LARGE BASIS cc-pVDZ 8. 1 O .0000000000 0.0000000000 -.2249058930 LARGE BASIS cc-pVDZ FINISH
Together with a simple job input (scf.inp):
**DIRAC .WAVE FUNCTION **HAMILTONIAN .DFT B3LYP **WAVE FUNCTION .SCF *END OF INPUT
We can get the wave function with the following run script:
pam --inp=scf.inp --mol=h2o.mol --get=DFCOEF
After running the script you should see the file DFCOEF in your submit directory.
Plotting the orbital density¶
Water has 5 orbitals and we will plot the HOMO, orbital nr. 5 (density.inp):
**DIRAC #.WAVE FUNCTION **HAMILTONIAN .DFT B3LYP **WAVE FUNCTION .SCF **VISUAL .3D 40 40 40 # will generate 40x40x40 cubes, the more points the better looking .DENSITY DFCOEF .OCCUPATION 1 1 5 1.0 .3D_INT # to double-check that density integrates to 2 electrons *END OF INPUT
The orbital occupation is controlled by the following keyword:
.OCCUPATION 1 1 5 1.0
This means that we have one occupation set (first “1”), this occupation set is for Fermion irrep 1 (second “1”), the range is orbital 5 (“5”) and it is fully occupied. To get the total density we could remove this keyword or alternatively:
.OCCUPATION 1 1 1-5 1.0
We will use the following run script to generate the cube file:
pam --inp=density.inp --mol=h2o.mol --put=DFCOEF --get=plot.3d.cube
Finally you can open plot.3d.cube with your favorite program.