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.