1 MIT 3.320 Atomistic Modeling of Materials Spring 2005 Lab 3: Handout Quantum-ESPRESSO: a first principles code, part 2. In this lab, we will be using Quantum-ESPRESSO as our first-principles code again. In problem 1, we will compare energy between allotropes of a transition metal and also evaluate stacking fault energy. In problem 2, we will examine the energetics of a perovskite structure, and how first- principles calculations can be applied to study a ferroelectric. Some helpful conversions: 1 bohr = 0.529177249 angstrom 1 Rydberg (R ) = 13.6056981 eV 1 eV =1.60217733 x 10-19 Joules We will once again be using hpcbeo2.mit.edu for this lab. Instructions on how to use the hpcbeo2 cluster were given in Lab 2. Please keep in mind that these calculations are more complex than the previous problem set and your calculations will take longer. So do not expect that you can get all the results before this session ends. *If you want to submit several jobs at a time, use a different $PREFIX for each job. Two jobs can go to one node. Create a directory for LAB3 $ cd ~/3.320 $ mkdir LAB3 $ cd LAB32 MIT 3.320 Atomistic Modeling of Materials Spring 2005 Problem 1 Problem 1 will look at calculations involving Cobalt. We will use the LDA exchange-correlation functional. Since Co is a ferromagnetic material, we will do spin-polarized calculations. To set up these calculations, you should have a basic knowledge of crystal structures. Chapter 1 of Introduction to Solid State Physics by Kittel is a good introduction. First, Create a directory for PROBLEM1 $ mkdir PROBLEM1 $ cd PROBLEM1 Next, copy files from /home/lee0su/3.320/LAB3/ $ cp ~lee0su/3.320/LAB3/Co.* . HCP structure c a b a Below is a copy of the file Co.HCP.scf.in, along with explanations. For the sake of brevity, I will not repeat entries which have the same meaning as the previous handout, except for items which you will be required to change.3 MIT 3.320 Atomistic Modeling of Materials Spring 2005 (1) &control (2) calculation = 'scf' (3) restart_mode='from_scratch' (4) prefix='Co.HCP' (5) tstress = .true. (6) tprnfor = .true. (7) outdir = '/state/partition1/lee0su' (8) pseudo_dir = '/state/partition1/lee0su' (9) / (10) &system (11) ibrav= 4 (12) celldm(1) = 1 (13) celldm(3) = 2 (14) nat= 2 (15) ntyp= 1 (16) ecutwfc = 30 (17) ecutrho = 250 (18) starting_magnetization(1) = 0.7 (19) occupations = 'smearing' (20) degauss = 0.03 (21) smearing = 'cold' (22) nspin = 2 (23) / (24) &electrons (25) mixing_beta = 0.7 (26) conv_thr = 1.0d-8 (27) / (28) ATOMIC_SPECIES (29) Co 58.933 Co.pz-nd-rrkjus.UPF (30) ATOMIC_POSITIONS (crystal) (31) Co 0.3333333333 0.6666666667 0.25 (32) Co 0.6666666667 0.3333333333 0.75 (33) K_POINTS {automatic} (34) 2 2 1 0 0 04 MIT 3.320 Atomistic Modeling of Materials Spring 2005 Line 11 This line contains information about the Bravais lattice. ibrav = 4 refers to the hexagonal lattice. Refer to INPUT_PW to see how each bravais lattice is defined. Line 12-14 There are two independent lattice parameters in a hexagonal lattice (as you probably know). Refer to the figure if you have forgotten how to construct an HCP structure. celldm(1) is a in bohr ( not angstrom ! ) and celldm(3) is c/a ( not the absolute value of c ! ). Find an appropriate value of celldm(3) from your knowledge of the ideal HCP structure. Two atoms comprise one unit cell. Line 16-17 We are going to use ultrasoft pseudopotentials here. In the case of norm-conserving pseudopotentials, ecutrho (the charge density cutoff) is automatically determined by 4*ecutwfc. However, in the case of ultrasoft pseudopotentials, we need an augmented charge around the ion core, so ecutrho should be higher than 4*ecutwfc. The usual value is 25-35 Ry for ecutwfc and 200-300 Ry for ecutrho. You might want to do a convergence check. Keep in mind that the value you should look at is the energy difference or force, not the absolute value of the energy (the energy will not converge unless you use very, very high ecutwfc and ecutrho). Lines 18 starting_magnetization is the starting magnetization for the atom. Set this to a value between -1 and +1. Lines 19-21 These keywords are particular details for the Brillioun zone integration for metals. Since there is a discontinuity of the occupation number for the bands around the Fermi energy, total energy with respect to the number of k-points converges very slowly. Adding electronic temperature (degauss) smooth out the abrupt change of the occupation number and as a result total energy converges with fewer number of k-points. Lines 22 nspin=2 turns ON spin polarization while nspin=1 turns it OFF. We will use nspin=2 throughout Problem 1. Lines 30-32 You should be cautious when you set up atomic positions.5 MIT 3.320 Atomistic Modeling of Materials Spring 2005 The default is alat , which means that you are using Cartesian coordinates and everything is scaled by celldm(1). In the case of HCP , it is easier to express coordinates with respect to crystal axes, which are a, b, and c in the figure. This way you don't need to change the atomic coordinates accordingly whenever you change celldm (3). Lines 33 This is k-point information. One thing to keep in mind is since the unit cell is longer in the c direction, a sparser k-mesh can be used in that direction. Problem 1-1 You should find the correct k-point grid yourself. There are two scripts, one for FCC and the other for HCP. It should be straightforward to check k-point convergence using these scripts. Co.FCC.scf.j Co.HCP.scf.j ( fixed c/a ratio ) To get total energy vs. k-point mesh, ( e.g. a = 4.74 ) $ grep ! Co.FCC.4.74.30.250.*.out | cut -c 24-25,63-90 > ksum.30.250 To get total energy vs. lattice parameter, $ grep ! Co.FCC.*.30.250.12.out | cut -c 12-15,63-90 > esum.30.250.12 If you want to do something more complex, refer to /home/lee0su/3.320/LAB3/getE.sh Make any necessary changes for your own purpose. Problem 1-2 To change a and c independently, use Co.HCP.covera.scf.j The script is almost the same as Co.HCP.scf.j and self-explanatory.6 MIT 3.320 Atomistic Modeling of Materials Spring 2005 Problem 1-4 To generate FCC structure in hexagonal cell, Make a small change ( celldm(3), nat ... ) to Co.HCP.scf.j Problem 1-5 Set up and run a calculation using Co.HCP.scf.j for a 5-atom cell ABCAC Check your result with the automatic script Co.ABCAC.j The script will do everything for you, after you specify the lattice