# 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.