Charge density due to doping atom in Si

This example shows how to use c2x to display the change in charge density caused by replacing an atom of Si with P in bulk silicon. The calculations are performed in Quantum Espresso, but Abinit, Castep, Siesta or VASP could have been used instead. The minimum version of c2x for this example is 2.40.

Running the calculations

The starting point was the eight atom silicon cell, si8.cell, given on the crystals page. This was then expanded to 2x2x2 64-atom supercell.

$ c2x -x=2x2x2 -Q --cell si8.cell si64.cell

By sorting the atomic co-ordinates (-Q) it is easy to edit the file to replace the central Si atom with P.

$ cp si64.cell si63p.cell
$ emacs si63p.cell

Separate subdirectories were prepared for the two QE runs. The bulk Si run is the easiest.

run_si64$ c2x --qe ../si64.cell
run_si64$ curl -o Si.UPF
run_si64$ pw.x -in > si64.log

The doped system is slightly harder, as QE defaults to assuming the system is insulating, and does not switch automatically to partial occupancies if there is an odd number of electrons. So one must edit the .in file produced by c2x and add to the SYSTEM namelist the lines:

  occupations = 'smearing',
  degauss = 0.2,

So the full sequence is

run_si63p$ c2x --qe ../si63p.cell
run_si63p$ emacs
run_si63p$ curl -o Si.UPF
run_si63p$ curl -o P.UPF
run_si63p$ pw.x -in > si63p.log

Calculating the Density Difference

As the same cutoff was used in both calculations, the grid sizes are the same, and it suffices to use

$ c2x -cv --diff run_si63p/pwscf.xml run_si64/pwscf.xml diff.xsf
Found 3D data for Charge
min=-0.0619541 max=0.330954 sum=69.9793 int=1

Had the grid sizes been different, it would have been necessary to add -i=0x0x0 to cause the coarser grid to be interpolated to the finer.

Note that the integral of charge difference is a single electron, as expected, and the there are places where there is less charge than before, as well as places where there is more.

In pictures

(The images below use XCrysDen. In order to show the bonds, it was necessary to increase the "connectivity factor" for Si using Modify|Atomic Radii|Chemical connectivity factor. It was set to 1.15.)

Isosurface of doping charge in Si

The charge density isosurface at 0.045e/A3 for P substitutionally doping Si in bulk Si. The extra charge (red) is localised around the single P atom at the centre of the image, and there is a slight reduction in charge (blue) in the bonds between the P and its Si neighbours.

Another isosurface of doping charge in Si

The charge density isosurface at 0.010e/A3. Now it can be seen that the effect of the addition of the P is not entirely localised, but that all the Si-Si bonds have lost a small amount of charge. It would seem credible that this system is now conducting (which, of course, P-doped Si is).

Not Science

The above is a demonstration of how to use c2x, not of how to do good science. No relaxation of the atoms around the phosphorous was performed, and no thought was given to convergence with respect to k-point density, basis set cut-off, or smearing width.