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