Visualize canonical orbitals from CPMD

For different reasons, one might be interested in looking at the Kohn-Sham canonical orbitals. In my case it was to identify the main orbitals involved TDDFT electronic transitions.

Step 1: Get the CPMD files for the orbitals to plot

The CPMD keyword used to generate the required file is RHOOUT BANDS in the &CPMD section.

I have written a script that reads a CPMD TDDFT output file and extract information regarding the orbital involved in the TDDFT electronic transitions. The script is called read_orbs_info and is part of the comp_chem_py package. The script prints some of the required lines to plot orbitals with CPMD. For example:

RESTART WAVEFUNCTION COORDINATES LINRES LATEST
KOHN-SHAM ENERGIES
    20
RHOOUT BANDS
    4
    -138 -139 -140 -141

Those lines have to be pasted in the &CPMD section of a CPMD input file. The restart information and number of Kohn-Sham energies is set to match the TDDFT output file that is read by the script.

Finally, the lines after RHOOUT BANDS are set such that a grid file will be generated for all the orbitals involved in the TDDFT electronic transitions.

Of course all the information can be modified to suit your purpose.

Step 2: Generate cube files

After the first step we should have a bunch of WAVEFUNCTION.index files available. 8 in our case. Those files have to be converted to the Gaussian cube format. This can be done using the cpmd2cube.x utility. This script is provided with the CPMD package and is located at:

/share/lcbcpc35/software/cpmd2cube/cpmd2cube.x

On the LCBC cluster.

To use it just run:

cpmd2cube.x WAVEFUNCTION.*

Or use the -h option for help. This operation should produce a cube file and a pdb file for each input, e.g. WAVEFUNCTION.8 leads to WAVEFUNCTION.8.cube and WAVEFUNCTION.8.pdb.

The cube files can then be visualized using e.g. the VMD package:

module load vmd
vmd WAVEFUNCTION.8.pdb

Then in the VMD terminal window:

vmd > mol addfile WAVEFUNCTION.8.cube
vmd > mol addrep 0
vmd > mol modmaterial 1 0 Opaque
vmd > mol modstyle 1 0 Isosurface  0.01 0 0 0 1 1
vmd > mol modcolor 1 0 ColorID 0
vmd > mol addrep 0
vmd > mol modstyle 2 0 Isosurface -0.01 0 0 0 1 1

Or using the vmd_plot_cube script from the comp_chem_py package:

vmd_plot_cube GEOMETRY.xyz -c WAVEFUNCTION.*.cube