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 –xvf 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’s 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—the charge density file Silicon.save—a save file with general info on the run Silicon.wfc—the 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—lattice type, celldm—lattice constant in bohrs, nat—# of atoms, ntyp—number 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) |