Lab 3: Quantum-Espresso - a First-Principles Code (Part 2)

1 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 firstprinciples 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 LAB3 2 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 b a 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 (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= (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 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 0 4 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 normconserving 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 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 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 parameters and the number of layers ( $nlayer ). It is very important for you to make an input file at least once on your own. Spend several minutes trying to understand what is going on in the script. Extra credit problem To get total energy, magnetic moment, and pressure, use /home/lee0su/3.320/LAB3/getM.sh 7 Problem 2 In this problem, we will be looking at BaTiO3 in the cubic phase and tetragonal phase. The atomic positions of the cubic perovskite structure are shown below. Figure 1. Representation of a cubic perovskite structure. In our case, O is present on each of the faces, Ba on each corner and Ti in the body-centered position. Take a look at the file BaTiO3.ion_dynamics_example.in 8 (1) &control (2) calculation = 'relax', (3) restart_mode = 'from_scratch', (4) prefix = 'BaTiO3', (5) tstress = .true., (6) tprnfor = .true., (7) pseudo_dir = '/state/partition1/yourusername/', (8) outdir = '/state/partition1/yourusername/' (9) / (10) &system (11) ibrav = 1, celldm(1) = 7.5, nat = (12) ecutwfc = 30.0, ecutrho = 240.0 (13) / (14) &electrons (15) diagonalization = 'cg', (16) mixing_mode = 'plain', (17) mixing_beta = 0.7, (18) conv_thr = 1.0d­8 (19) / (20) &ions (21) ion_dynamics = 'bfgs' (22) / (23) (24)ATOMIC_SPECIES (25) Ba 137.327 Ba.UPF (26) Ti 47.88 (27) O Ti.UPF 15.9994 O.UPF (28) (29)ATOMIC_POSITIONS (30) Ba O.0 0.0 0.0 0 0 0 (31) Ti 0.5 0.5 0.5 0 0 0 (32) O 0.5 0.5 0.0 (33) O 0.5 0.0 0.5 (34) O 0.0 0.5 0.5 (35) (36)K_POINTS automatic (37) 4 4 4 1 1 1 5, ntyp = 3, 9 Line 2 For parts (b) and (c), instead of a single self-consistent field calculation, we will be doing a 'relax' calculation which includes a series of SCF calculations. Here, the ions are allowed to move in order to reduce the total system energy. This setting is very much like setting the 'opti' flag in GULP. Note that for part (a), you should use 'scf', as in problem1. Lines 7-8 Remember to set your scratch directory correctly. Lines 20-22 Since we will be using ion “dynamics”, we now need the new IONS section. We are not, however, using real dynamics(i.e. there is no time coordinate used in the relaxations), but just searching for the minimum energy relaxations. This section should be omitted for the scf calculations of part (a). Lines 33-34 You will notice that there are now three additional flags at the end of our atomic coordinates. These flags define the degrees of freedom available for those atoms during relaxation (zero = disallow motion in that direction for that atom). In the example shown above, the Ba and Ti ions are fixed, while the O ions are allowed to relax. So the format is: atomic label pos_x pos_y pos_z allow_x allow_y allow_z You will find that using scripts will save you tons of time on this problem set. Not using scripts would mean that you could spend literally days sitting at a computer waiting for runs to finish. You should be able to set up the appropriate scripts based on the ones we've already used for examining Fe and also the example script from last problem set. If you are still uncomfortable with writing your own scripts, make an appointment with your TA to work out some basic scripts for completing this problem set. A sample shell script for submitting a job using this example input file can be found in the ~brandonw/LAB3 directory. You should copy this file and modify it appropriately.