Plot densities of states with castep

Castepdos produces densities of states by interpolation of a Monkhorst-Pack k-point set used in castep. The necessary input file is the seedname.bands file. If your castep run used a symmetry-reduced set of k-points, you will also need to create a .symm file. Written by Graeme Ackland, please acknowledge authorship if you use it, and report bugs to gjackland@ed.ac.uk . I'll write up the method in my copious spare time.

Released 4.11.2010

Castepdos is distributed in a fortran files, in this directory castepdos.f

To compile castepdos you simply need a vanilla fortran compiler, e.g.

gfortran -o castepdos castepdos.f

I haven't tested it across many compilers, but there's nothing very sophisticated in the code.

Not using symmetry

To run the code, you first need to do a standard castep run with a Monkhorst-Pack k-point grid and no symmetry which produces a .bands file. This contains almost all the information castepdos needs. If you then run castepdos it will first prompt you for the seedname e.g.

/Home/gja/castep/Ti>castepdos
Enter castep seedname

Type in enter the seedname: this enables castepdos to look for the .bands file and (optionally) the .symm file as described later.

castepdos will then prompt for the mp k-point grid (annoyingly, this isn't contained in the .bands file). Enter the three numbers as requested, the same as in your .cell file (or .castep file under MP grid size for SCF calculation

The program will produce a file called seedname.dos. This contains three columns:

the energy in eV, with zero set to the Fermi Energy
the density of states at the energy
the cumulative density of states.

Which you can plot using your favourite graphics package. For magnetic systems spin_polarized : true there is another file called seedname.dos2.

If you are working in the Northern Hemisphere, this is the spin down band, the other being spin up. Researchers in the Southern Hemisphere will find these designations reversed.

Using symmetry

Castepdos can also use files created using symmetry, which is a good thing to do because it saves you computing time. Unfortunately, CASTEP doesn't write out the symmetry in any sensible place, so you have to do a bit of work to extract it.

The slow, tedious way

First, you need to run with high verbosity, include in the .param file:

iprint : 2

and in the .cell file...

symmetry_generate

This will produce lots of verbosity in the .castep file (sorry). You now need to create a .symm file for yourself.

First open the seedname.castep file in an editor, then copy the chunk containing the symmetry elements starting with the string ``1 rotation'' into a file called seedname.symm. An example for fcc is attached.

Now run castepdos as above, it should read both the .bands and the .symm files and produce a seedname.dos as before.

Using symmetry_snap

An better alternative to creating a .symm file is to use the symmetry_snap utility supplied with CASTEP5.5. This creates a file called seedname-symm.cell. If such a file exists, castepdos will use it to extract the symmetry information.

This method does not require you to use iprint : 2 , it needs only the .cell and .bands files. You still need to input the k-point grid at runtime.

Known Bugs.

castepdos is standard Fortran 77, except that it uses the F90 fortran intrinsic "trim" which is not implemented on some old compilers. It has no other known bugs. If you report unknown ones I'll fix them.

Known non-features:

There is no manual in a nice word processing format.

castepdos does not produce partial densities of states.

castepdos does not work with other DFT codes, although interfaces should be fairly easy to write, there is the usual pronoun problem.

Resolution

The accuracy of the DOS depends on the number of k-points you use. There are also two magic numbers in the code, NHIST and EFINE.

NHIST determines the size of the bins in the histogram used for the DOS plot. The default is 2000 which is fine for a single band. If your calculation includes semi-core electrons you might end up with a very wide energy range, in which case changing NHIST might be a good idea.

EFINE sets the trapezium width for the xy integration (z is done explicitly). If you increase it, you'll get a better DOS but the code will run slower.