XC-Functionals

Introduction

These exercises are intended to introduce you to some of the XC functionals available in CASTEP. They will show you how to choose and set up calculations using these functionals. They will also demonstrate the effect of chosen functional upon different physical properties.

It is important to remember throughout that no XC functional is perfect (as discussed in the XC functional talk). In choosing a functional, we must be guided by an understanding of the XC functional itself and the physics of the system under consideration. It is also important to realise that sometimes ''inappropriate functionals get the 'right' answer for the wrong reasons!'' Some of these exercises will further illustrate this point.

Many of the post-LDA functionals (HF, sX-LDA) are computationally intensive. For that reason, this practical session focusses upon small systems. However, by the end of this session you should be able to set up more complex and demanding systems.

Where To Find Help

If you want more information about a particular CASTEP keyword, or you want to find if CASTEP has particular functionality, there are a few places you can look. There is information on this website under documentation or you can use the search bar at the top of the page. CASTEP has an in built help option to assist with using particular keywords. Information on using CASTEP can be seen by using:

$ castep -h

To get more information on a particular input file keyword (e.g. kpoints_mp_grid) use:

$ castep -h kpoints_mp_grid

If you don't know the keyword you need to use, then you can search on a particular keyword. This returns a list of keywords that you might be interested in, e.g. to look at all keywords which contain a reference to symmetry.

$ castep -s symmetry

Finally, to list all keywords, use:

$ castep -h all

Example 1 - Si: LDA and GGA

In this first exercise, we shall explore how to use different XC functionals to determine the electronic structure of silicon.

Get the files required for this exercise : Si2 files

We shall carry out these calculations on the VM. Copy the files from the /course_materials directory to somewhere in your home directory and then unpack them. You can use:
```
cp /course_materials/Si2.tar.gz .
tar -xvf Si2.tar.gz
```
We will begin by performing a simple band structure calculation. To do this, we must edit the Si2.param file. Open this file. Find the line that says task:. Amend this to say:
```
task            spectral
spectral_task   bandstructure
```
We can also choose the XC functional that we wish to employ in the param file. If we set
```
xc_functional : lda
```
then we will run a calculation using the LDA. This is the functional that will be employed in the SCF calculation that determines the ground state electronic density. If we do not specify otherwise, then it will also be used for the band structure calculation.
Now we shall examine the .cell file. Open Si2.cell. Later in this tutorial, we shall use non-local functionals. In CASTEP, non-local functionals can only be used with norm-conserving pseudo potentials. (Question: Do you understand as to why this should be so? Ask if you do not.) For this reason, ensure that you are using the norm-conserving pseudo potential library by including:
```
%block species_pot
NCP
%endblock species_pot
```
in your .cell file.
We will also need to set up an appropriate k-point path for the band structure calculation. We will use the following path:
```
%block spectral_kpoint_path
0.5 0.5 0.5
0.0 0.0 0.0
0.5 0.0 0.5
0.625 0.25 0.625
0.0 0.0 0.0
%endblock spectral_kpoint_path
```
Add this block to the .cell file.
Ensure that your calculation is converged with respect to both kinetic energy cut-off and k-point sampling.
Now run CASTEP on the VM using the 2-atom input files.
```
$ mpirun -n Number castep.mpi Si2 
```
where Number should be replaced by the number of processors that you wish to use (up to a maximum of 16).

This should only take a few seconds and produce a readable output file Si2.castep. Examine this file and try to understand the meaning of the various parts. In particular check the section following the header which lists all of the input parameters, both explicit and default. Note what default values of the major parameters CASTEP chose where you did not specify them explicitly. (There will be some whose meaning has not been explained. Don't worry about these.) Note down how long the calculation takes for each k-point.

We will now analyse the band structure. We can produce a plot of this from the Si2.bands file using:
```
dispersion.pl -xg -np -bs Si2.bands | xmgrace -
```
Using the band structure plot, determine the band gap. How does this compare to the experimental value? Can you explain this?

I recommend that you save the output from this calculation as Si2_lda or similar as later on we wish to compare with GGA and non-local calculations.

Repeat steps 1-8 using a GGA functional. To change the functional employed, amend the .param file as follows:
```
xc_functional : pw91 
```
There is no single unique flavour of GGA, and this specifies the Perdew-Wang 91 implementation of the GGA. Other flavours such as PBE would be specified in similar fashion.

How does the use of the GGA affect the band gap? Does it improve upon the LDA value? Why?

Example 2 - Si: Non-local functionals

We shall now move onto examining what happens with a non-local functional. In principle, we should, for a new functional, determine a ground state electronic density that is consistent with the functional employed. However, CASTEP allows us to be more flexible and to specify a different functional for the band structure calculation. We will use this to save some computational effort by using a LDA-derived ground state density.

To begin, carry out a single-point energy calculation on our silicon system using the LDA.

This calculation will be very quick. Alongside the .castep file, CASTEP there will also be a file Si2.check. This is the restart file.

We will use this restart file to perform a band structure calculation using a different functional from that employed for the SCF.
In the Si2.param file, add the following line:
```
continuation : default 
```
This tells CASTEP that this is a re-start, and that the restart file is the default, which in this case is Si2.check. We could specify a different file here, if we wished.
Add the following line to the .param file:
```
spectral_xc_functional : HF
```
This allows us to specify a different functional for the band structure calculation, which in this case will be a Hartree-Fock calculation. However, we are using a LDA-derived ground state density to feed into this calculation.
Now run the band structure calculation as before.

When examining the .castep output look at the band structure timings. How do these compare to the LDA and GGA cases performed earlier? This should give some idea as to the more computationally intensive nature of a non-local calculation. Similarly, how do the memory estimates change? It is important to understand that this increased memory overhead can limit non-local calculations, even on HPC platforms (indeed, it can be more of a problem than the increased computational effort per step of the calculation).
Plot the resulting band structure. How does this compare to the LDA and GGA values? How does the band gap compare? Why is this?
Repeat steps 1-5, but for screened exchange with LDA correlation (sX-LDA). To do this, in Step 3 set
```
spectral_xc_functional : sX-lda
```
Example 3 - FeO and DFT+U

In the XC Functional talk, data were shown illustrating the effect of functional upon the electronic structure of antiferromagnetic FeO. In this example you will learn how to perform DFT+U calculations.

Get the files required for this exercise : FeO Files

We shall carry out these calculations on the VM. Copy the files from the /course_materials/FeO folder to somewhere in your home directory. For example:
```
cp -r /course_materials/FeO .
```
Open the .cell file. We are interested in examining an anti-ferromagnetic structure. Note the following in the block where the atomic positions are specified:
```
Fe 0.0  0.0  0.0 spin=-4.0
Fe 1.5 -0.5 -0.5 spin=4.0
```
This specifies the initial antiferromagnetic spin arrangement, with spin = $N_\text{up}-N_\text{down}$ defining the initial spin polarisation on each Fe atom.

NOTE: the spins specified here should be consistent with the overall spin (i.e. the initial number of unpaired electrons) specified in the .param file using the spin keyword. In this case, as the overall spin is zero (which is the CASTEP default value), we do not need to specify this, but if we examined a ferromagnetic structure then we would.

CONVERGENCE NOTE: I suggest using a 6x6x6 MP grid. While not as converged as one would desire for serious science, this is sufficient for demonstrative purposes, and will allow us to perform several calculations in the time available.
As we saw in the XC presentation, PBE obtains a metallic ground state for antiferromagnetic FeO. We shall therefore run this calculation as a metal. Ensure that the .param file has the following line:
```
fix_occupancy : false
```
which allows CASTEP to vary electronic occupancies. We shall use the default metals method, which is density mixing. Density mixing is the only method in CASTEP that allows us to specify initial magnetic moments on atoms.

The convergence of calculations on metallic systems is often improved by the inclusion of a number of extra empty bands. We can specify this using the nextra_bands keyword. Set this to 14 with:
```
nextra_bands : 14
```
4. As we wish to examine an anti-ferromagnetic structure in which individual Fe ions have non-zero spin, we must allow the system to be spin-polarised. To do this, open up and edit the .param file. Set:
```
spin_polarized : true
```
5. We will set the XC functional to be the PBE flavour of GGA (xc_functional : pbe).
Now run a band structure calculation and produce a band structure plot as before.

There are certain features of the band structure plot that you should understand: the dispersive bands are ''sp'' bands derived from ''s'' and ''p'' electrons. The flat non-dispersive bands are derived from Fe ''d'' electrons. Can you explain why the ''d'' bands are so flat?
When examining magnetic systems, it is also often useful to examine the spin on each atom. This information can be found in the .castep file in the section
```
Atomic Populations (Mulliken)
```
Note the spin on each Fe atom. How does this vary from the initial value specified?

A DFT+U calculation

To do this, open up and edit the .cell file. Add the following block to it:
```
%block hubbard_u 
Fe 1 d: 2.5 
Fe 2 d: 2.5 
%endblock hubbard_u
```
This specifies a Hubbard U parameter of 2.5 eV for each Fe ion in the system. Whilst schemes exist to determine the value of U from first principles, in CASTEP it is simply specified as a parameter. Of course, one could use ab initio calculations to determine a value of U then specify this in the .cell file. Beware, however, of using a value of U value determined for one XC functional or one DFT code with another as they are not transferable between methods.
Run a band structure calculation and plot the results. How does this band structure compare to the PBE+U result? Can you explain which features change and which stay the same? What happens to the band gap?
Now examine the spin on each Fe atom. Compare to the PBE values. Can you explain your observations?
Repeat steps 7-9 for differing U values. Explain your findings.

Example 4 - Graphite and DFT+D

Graphite represents a prototypical example of a layered system in which the layers are weakly bound by van der Waals forces. As is well known, local and semi-local functionals such as the LDA and GGA neglect such interactions, and we could perhaps anticipate that they therefore perform poorly for such a system. We will explore the performance of these functionals. This example will also show you how to run dispersion corrected DFT+D functionals.

Get the files required for this exercise : Graphite files

As before, we shall carry out these calculations on the VMs. Copy the files from the /course_materials directory to somewhere in your home directory and then unpack them. You can use:
```
cp /course_materials/Graphite.tar.gz .
tar -xvf Graphite.tar.gz
```
Open the .param file using a text editor. We will do our first calculation with the PBE GGA functional. Ensure that this is the xc functional specified (see above if you cannot recall how to do this).
Graphite has a semimetallic nature, with a non-zero DOS at the Fermi energy for the K and H points in reciprocal space. We will therefore run this calculation as a metal, as we did in the FeO example. Modify the .param file to ensure that this is the case. In contrast to FeO, however, we do not need to worry about spin polarisation.

The convergence of calculations on metallic systems is often improved by the inclusion of a number of extra empty bands. We can specify this using the nextra_bands keyword. Set this to 14 with:
```
nextra_bands : 14
```
We will run a geometry optimisation. Change the .param file so that this is the task.
Now open up the .cell file using a text editor. We will perform our calculations for the experimental graphite structure. However, our atomic positions are slightly incommensurate with the symmetry of the unit cell; to remedy this, add the following line to the .cell file:
```
snap_to_symmetry
```
which puts the atoms on the symmetry positions.
Now run a variable cell geometry optimisation. If you cannot recall how to do this, consult the earlier tutorial concerning geometry optimisations.

Ensure that the calculation is converged with respect to both plane wave and k-point sampling. As we are interested in a geometry optimisation, examining the convergence with respect to force is a suitable criterion. However, as we have placed the ions on high symmetry positions, the forces on them will be zero by symmetry - can you think of how we can still examine the convergence of forces, given this fact?
Once the calculation has completed, examine the unit cell parameters - how do they compare to the original (experimental) parameters? Can you explain this behaviour?
We shall now investigate what happens with LDA. Repeat the above calculation with the LDA as the XC functional.
What happens to the unit cell parameters in this case? Pay particular attention to the c axis. How does this compare to the GGA value? Is this behaviour expected? Can you explain it? Does this mean that LDA is particularly suitable to use to describe van der Waals bound systems?
We shall now carry out a GGA calculation with a DFT+D correction. For this you can leave the .cell file unaltered. Open up the .param file using an editor. Add the following line:
```
sedc_apply : true
```
which tells CASTEP that it should apply the semi-empirical dispersion correction (i.e. DFT+D).
Now specify the correction scheme. Begin with the G06 scheme (also known as D2):
```
sedc_scheme : g06
```
Ensure that the XC functional to use is the PBE.
Run the geometry optimisation. How do the results compare to the LDA and GGA answers?
Repeat 10-13 but using the D4 correction scheme. How do the results change?

NOTE: the OBS scheme is not compatible with the PBE functional.

You can also compare your findings to the other tutorial here on DFT+D where the layers of graphite are separated. The Jupyter notebook provided there can also be used on the VMs as an alternative to running CASTEP manually.

Further examples

If you have worked through to this point in the tutorial, then feel free to apply these functionals to some more sophisticated systems that interest you. If you do not have particular systems in mind, then feel free to discuss possible choices with a demonstrator.