Using PWscf for density functional calculations

 

Today, I will be going over a few examples of systems that you can calculate with the PWscf code from the Quantum Espresso suite. To obtain your own copy of this code, you can visit their website, www.pwscf.org.

 

Density functional approaches are typically used to find the band structure of a material and also its ground state energy. For simple crystals, the ground state can be found by adjusting the lattice constant and fitting the minimum energy. For more complex systems, it is often important to use forces to relax the system.

 

In these first calculations, we will focus on finding the minimum energy for a particular configuration. PWscf reads in an input file that describes several features of the system under study. This can include things like lattice constant, atom type, pseudopotential used, and parameters related to the approach.

 

For our work today, I will focus on two key parameters, the energy cut-off and the k-point grid used. One of the fundamental tests of DFT calculations is to insure that your ground state energy is converged with respect to these two parameters.

 

Energy cut-off

PWscf uses a plane-wave basis set to describe the ground state wavefunction for the system. The number of plane waves that we use to describe the system is determined by an energy cut-off parameter. The more plane waves we use, the more accurate the calculation. However, this also increases the calculation time.

 

You can think of this plane wave expansion in similar terms to descriptions of Fourier expansions. Relatively smooth curves are easy to describe. Curves with sharp corners, such as a square step, require a large basis set to describe accurately.

 

The same thing applies with the pseudopotentials that are used in density functional calculations. If we look at the real potential around the atom, electron wavefunctions near the nucleus have wide oscillations to insure orthogonality. These oscillations are very difficult to describe in a plane wave basis. However, we are only really interested in the phase shift in the wavefunction that interactions with the atom produces. So if we can find a simpler potential or pseudopotential that generates the same phase shift, we should be able to use fewer plane-waves and speed up our calculations.

 

Pseudopotentials for atoms usually come in two varieties, hard and soft. Hard pseudopotentials are deeper and require more plane waves and a higher energy cut-off. Soft or ultra-soft pseudopotentials uses much lower energy cut-offs but relax certain conditions.

 

K-point grid

Since we are dealing with periodic structures, it is often easier to describe our system in momentum or k-space. The region of available k points for a system is known as the Brillouin zone (BZ). Many of the physical parameters associated with a system, (i.e. total energy, charge, magnetic moment, etc) can be expressed in terms of integrations over k-points in the BZ. For DFT calculations, it is essential that a sufficiently large k-point grid is used in the calculations. Note that the k-point grid does not have any direct relationship with the basis we will use. So k-point integrations will be important for other approaches (LMTO, KKR, etc) as well.

 

The k-grid is usually expressed in terms of in terms of a Monkhorst-Pack grid. These grids describe divisions in the kx, ky, and kz directions. Since semiconductors do not have a Fermi surface, we can get away with a fairly small k-point grid. Metals usually require much better k-point grids to insure convergence (i.e. 16 16 16 or greater).

 

Running calculations

You will need to copy the file pwscf_workshop.tar.gz into your home directory. Then you should create a directory for our work today.

 

mkdir yourname_pwscf

 

Move the tar.gz file into that directory:

 

mv pwscf_workshop.tar.gz yourname_pwscf/

 

Now we need to get the contents out of the file

 

gunzip pwscf_workshop.tar.gz

 

tar 堀vf pwscf_workshop.tar

 

This will create a directory called workshop_runs/

 

This directory contains a number of different input files and the associated pseudopotentials for the atoms. In particular, we have input files for silicon, aluminum, nickel, CO molecule, and an aluminum surface. You will notice that silicon, aluminum and nickel have two input files associated with them. The scf version is used to find the ground state energy and the band version is used to plot the band based on the self-consistent density.

 

Let痴 try running the silicon case.

 

pw.x < si.scf.cg.in > scf.out &

 

If we look at our directory, we will see that a number of files have been created:

 

Silicon.pot - the potential file

Silicon.rho葉he charge density file

Silicon.save預 save file with general info on the run

Silicon.wfc葉he wavefunction file for the system

 

We can also look at the convergence properties of the calculation. DFT calculations perform a self-consistent loop between calculating the potential (Quantum Mechanics side) and the charge (Electrostatics side). We can look at the convergence in energy by grepping for that term.

 

nnin-login-01% grep energy si.scf.out

kinetic-energy cutoff = 18.0000 Ry

total energy = -14.54442731 ryd

total energy = -14.56252129 ryd

total energy = -14.56383710 ryd

total energy = -14.56385997 ryd

total energy = -14.56386162 ryd

! total energy = -14.56386165 ryd

 

You can also look at the convergence by grepping scf

nnin-login-01% grep scf si.scf.out

Current dimensions of program pwscf are:

estimated scf accuracy < 0.34803700 ryd

estimated scf accuracy < 0.01086909 ryd

estimated scf accuracy < 0.00026578 ryd

estimated scf accuracy < 0.00000338 ryd

estimated scf accuracy < 0.00000016 ryd

estimated scf accuracy < 4.0E-09 ryd

 

Looking at the input file si.scf.cg.in

覧覧覧覧覧覧覧覧覧覧覧覧覧覧覧覧覧覧覧覧覧覧

&control Section that describes what type of job we will do (scf, phonon, nscf)

calculation = 'scf' We are doing self consistent calculation.

restart_mode='from_scratch',

prefix='silicon', Prefix for all output files.

tstress = .true. Calculate stress in system.

tprnfor = .true. Print out the forces

/

&system

ibrav= 2, celldm(1) =10.20, nat= 2, ntyp= 1,

ibrav様attice type,

celldm様attice constant in bohrs, nat# of atoms,

ntyp溶umber of different atom types

ecutwfc =18.0, Energy cut-off for the wavefunction in Ryd

/

&electrons

diagonalization='cg'

mixing_beta = 0.7 mixing factor for self consistent cycles

conv_thr = 1.0d-8 convergence threshold in Ryd

/

ATOMIC_SPECIES

Si 28.086 Si.vbc.UPF atom name, atomic weight, and potential

ATOMIC_POSITIONS

Si 0.00 0.00 0.00

Si 0.25 0.25 0.25

K_POINTS

4 4 4 1 1 1 Monkhorst pack grid (nx, ny, nz, shift_x, shift_y, shift_z)