PWSCF Guide
Introduction
We will be using the Quantum-Espresso package as our first-principles code. Quantum- Espresso is a full ab initio package implementing electronic structure and energy calculations, linear response methods (to calculate phonon dispersion curves, dielectric constants, and Born effective charges) and third-order anharmonic perturbation theory. Quantym-Espresso also contains two molecular-dynamics codes, CPMD (Car-Parrinello Molecular Dynamics) and FPMD (First-Principles Molecular Dynamics). Inside this package, PWSCF is the code we will use to perform total energy calculations. PWSCF uses both norm-conserving pseudopotentials (PP) and ultrasoft pseudopotentials (US-PP), within density functional theory (DFT). This is a tutorial on how to get energies using the PWSCF code in Quantum-Espresso.
Quantum-Espresso is free under the conditions of the GNU GPL. Further information (including online manual) can be found at the Quantum-Espresso website http://www.quantum-espresso.org/. Quantum-Espresso is under very active development. There are many other first-principles codes that one can use. These have been listed in the lecture slides, and include:
ABINIT (http://www.abinit.org): DFT Plane wave pseudo-potential (PP) code. ABINIT is also distributed under the GPL license.
CPMD (http://www.cpmd.org): DFT Plane wave pseudopotential code. Car-Parrinello molecular dynamics. GPL.
SIESTA (http://www.uam.es/departamentos/ciencias/fismateriac/siesta/): DFT in a localized atomic base.
VASP (http://cms.mpi.univie.ac.at/vasp/) DFT Ultrasoft pseudo-potential (US-PP) and PAW code. US PP’s are faster and more complex than corresponding PP codes, with similar accuracy. Moderate cost for academics (~$2000)
WIEN2k (http://www.wien2k.at/): DFT Full-Potential Linearized Augmented Plane-Wave (FLAPW). FLAPW is the most accurate implementation of DFT, but the slowest. Small cost for academics (~$500)
Gaussian (http://www.gaussian.com): Quantum Chemistry Code - includes Hartree-Fock and higher-order correlated electrons approaches. Moderate cost for academics (~$3000). Has recently been extended to solids.
Crystal (http://www.cse.clrc.ac.uk/cmg/CRYSTAL): HF and DFT code. Small cost for academics (~$500).
Unit Conversions
Important note! PWSCF reports distances in bohr and energy in rydbergs:
1 bohr = 1 a.u. (atomic unit) = 0.529177249 angstroms.
1 Rydberg (R ) = 13.6056981 eV
1 eV =1.60217733 x 10-19 Joules
Getting Started
For the first half of Lab 2, you only need two files. Navigate to /projects/e20438/lab2/ and copy C.scf.in and pwscf.q to your directory of choice.
Theoretical Background: Problem 1
Remember that we are dealing with infinite systems using periodic boundary conditions. This means that we can use the Bloch theorem to help us solve the Schrödinger equation. The Bloch theorem says:
with
In these equations, is the wavefunction (Psi_n_k), is a function that is periodic in the same way as the lattice (G), G sums over all (at least in principle) reciprocal lattice vectors, and c_G’s are coefficients in the expansion. In this case, our basis functions (ie what we expand in) are plane waves. They are called “plane waves” because surfaces of constant phase are parallel planes perpendicular to the direction of propagation. Remember that limiting the plane wave expansion to the infinite, but numerable and discrete set of G vectors that are integer multiples of the three primitive lattice vectors, we are selecting plane waves that have a periodicity compatible with the periodic boundary conditions of our direct lattice.
In actual calculations, we have to limit the plane wave expansion at some point (i.e. stop taking more G 's). This is called the planewave cutoff. Cutoffs are always given in energy units (such as Rydberg or eV) corresponding to the kinetic energy of the highest G. Note: The units of reciprocal lattice are the inverse of the direct lattice, or 1/length. However, we can convert 1/length to energy units (Remember λν = c, and hν = E. λ is wavelength, and is in units of 1/length). Problem 1 is designed to test cutoff convergence issues. You can always take a higher cutoff than you need, but the calculations will take longer.
Theoretical Background: Problem 2
Because of the Bloch theorem, we need to solve a Schrödinger-like Kohn-Sham equation (i.e. iterate selfconsistently the diagonalization of a MxM matrix, where M is the number of basis functions) everywhere in the Brillouin Zone. In practice, we do it for a finite number of k wavevector values (e.g. the coarse grained Monkhorst-Pack mesh), and get a value for E at each k . This is seen in the schematic below.
The picture goes over the first Brillouin zone. (If you don’t know what a Brillouin zone is, it’s time to go over the concepts of direct and reciprocal lattice). To obtain a value for E, the energy of the crystal, we need to integrate over the first Brillouin zone, where the bands are occupied (and divide by the volume). Thus, summing over a finite number of k points is an approximation to performing an integral. You will need to make sure you have enough k -points to have a converged value for the energy.
Summary of Background
For all first-principles calculations, you must pay attention to two convergence issues. The first is energy cutoffs, which is the cutoff for the wave-function expansion. The second is number of k -points, which measures how well your discrete grid has approximated the continuous integral.
Input Files
We will look at the C.scf.in input file, which is contained in the script that we will use later for the problems. This is an input file for carbon in its diamond form. Diamond is a face-centered cubic structure with two C atoms at 0 0 0 and 0.25 0.25 0.25 (of course, there would be no difference in having those two atoms in 0.13, 0.20, 0.22 and 0.38, 0.45, 047. Why ?). Diamond looks like this (a is the lattice parameter).
Look at the input file, which you have just copied over to your directory. The file will look something like this:
(1) &control
(2) calculation = 'scf'
(3) restart_mode='from_scratch'
(4) prefix='diamond'
(5) tstress = .true.
(6) tprnfor = .true.
(7) pseudo_dir = '/usr/local/espresso-4.0/pseudo'
(8) outdir = './'
(9) /
(10) &system
(11) ibrav= 2, celldm(1) =6.60, nat= 2, ntyp= 1
(12) ecutwfc =110
(13) /
(14) &electrons
(15) diagonalization='david'
(16) mixing_mode = 'plain'
(17) mixing_beta = 0.7
(18) conv_thr = 1.0d-8
(19) /
(20) ATOMIC_SPECIES
(21) C 12.011 C.pz-vbc.UPF
(22) ATOMIC_POSITIONS
(23) C 0.00 0.00 0.00
(24) C 0.25 0.25 0.25
(25) K_POINTS {automatic}
(26) 4 4 4 0 0 0
Note: k-point grid has been changed to 4x4x4 and line numbers are added for reference:
Line 1: &control declares the control block.
Line 2: calculation = 'scf' tells PWSCF that this will be a self-consistent field calculation.
Line 3: restart_mode = 'from_scratch', declares that we will be generating a new structure.
Line 4: prefix='diamond', declares the filename prefix to be used for temporary files.
Line 5: tstress = .true. is a flag to calculate the stresses.
Line 6: tprnfor = .true. is a flag to calculate the forces.
Line 7: pseudo_dir defines the location of the directory where you store the pseudo-potentials.
Line 8: outdir='./' defines the location of the temporary files. This should always be a local scratch disk so that large I/O operations do not occur across the network.
Line 9: / denotes the end of a block.
Lines 10-13: the system block
ibrav – gives the crystal system. ibrav=2 is a face-centered cubic structure. This is used because the symmetry of the structure can reduce the number of calculations you need to do. If you need other crystal systems, consult the PWSCF manual.
celldm – defines the dimensions of the cell. You will be changing this parameter. celldm is in atomic units, or bohrs. Remember that 1 bohr = 0.529177249 angstroms. The value will depend on the Bravais lattice of the structure (see below for a key). For FCC, celldm(1)=a. In cubic systems, a=b=c.
nat – number of atoms (each individual unique atom). Note that diamond has two unique atoms in the smallest asymmetric unit.
ntyp – number of types of atoms
ecutwfc – Energy cutoff for pseudo-potentials. This one is important; you will be changing this parameter.
Lines 14-19: The electrons block
diagonalization – diagonalization method. Use the default for now.
mixing_beta - Mixing factor. Don’t worry about this for now.
mixing_mode – Mixing method. Don't worry about this for now.
conv_thr – Convergence threshold. Use the default for now.
Lines 20-21: Atomic species declaration. After the keyword ATOMIC_SPECIES, for each ntyp enter: "atomic symbol atomic weight pseudo-potential"
Lines 22-24: Atomic positions. After the keyword ATOMIC_POSITIONS, for each nat enter "atomic symbol x y z" where x,y,z are given as fractional coordinates of the conventional cell.
Lines 25-26: k-point selection. After the keywork K_POINTS, “automatic” tells PWSCF to automatically generate a k-point grid. The format of the next line is nkx nky nkz offx offy offz where nk* is the number of intervals in a direction and off* is the offset of the origin of the grid.
Interpreting Output
Look at one of your .out files generated by issuing a qsub command with the provided pwscf.q file. Scroll through the output file. The beginning will just recap the configuration that is being calculated. Then there is some information about the pseudo-potentials that PWSCF just read in. The next part tells you about intermediate energies that PWSCF calculates, before the calculation is converged. Near the end, there will be something like:
! total energy = -22.51880584 ryd
(your number may be slightly different). Note, the final occurrence of “total energy” will have an exclamation point by it, something you can use to hunt for it. You can skip to this right away by using search functions (in vi, type esc, /total energy, enter), or just scroll down a lot. This is your total energy, as calculated by PWSCF. Or you can type
$ grep “\!” *.out
Which will grep all total energies from all files terminated by “.out” in the current directory. The “!” is a special character, which is negated by using “\”. At the end of the file, it will tell you how long your program took to run, and how much memory it used.
init_run : 0.14s CPU
electrons : 0.81s CPU
forces : 0.00s CPU
stress : 0.02s CPU
electrons : 0.81s CPU (6 calls, 0.100 s avg)
c_bands : 0.60s CPU ( 6 calls, 0.023 s avg)
sum_band : 0.14s CPU (13 calls, 0.005 s avg)
v_of_rho : 0.06s CPU
... etc ...
PWSCF : 0.98s CPU time
Your numbers may be different from these. You should start to develop a feel for how long your runs take, and how much memory they will use.
There will also be lines which say:
Cartesian axes
site n. atom positions (a_0 units)
1 C tau( 1) = ( 0.0000000 0.0000000 0.0000000 )
2 C tau( 2) = ( 0.2500000 0.2500000 0.2500000 )
number of k points= 8
cart. coord. in units 2pi/a_0
k( 1) = ( 0.0000000 0.0000000 0.0000000), wk = 0.0312500
k( 2) = ( -0.2500000 0.2500000 -0.2500000), wk = 0.2500000
k( 3) = ( 0.5000000 -0.5000000 0.5000000), wk = 0.1250000
k( 4) = ( 0.0000000 0.5000000 0.0000000), wk = 0.1875000
k( 5) = ( 0.7500000 -0.2500000 0.7500000), wk = 0.7500000
k( 6) = ( 0.5000000 0.0000000 0.5000000), wk = 0.3750000
k( 7) = ( 0.0000000 -1.0000000 0.0000000), wk = 0.0937500
k( 8) = ( -0.5000000 -1.0000000 0.0000000), wk = 0.1875000
This records the number of unique k-points that were calculated. The input here has a 4x4x4=64 k-point mesh, however some of these k-points have the same energy because of the symmetry of the crystal. They are then weighted differently. The weights add up to 2, an idiosyncrasy of this program. Your numbers will be different if you use a different k-point mesh.