IntroductionWe 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.
Unit ConversionsImportant note! PWSCF reports distances in bohr and energy in rydbergs:
Getting StartedFor 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 1Remember 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,  sums over all (at least in principle) reciprocal lattice vectors, and ’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.
Theoretical Background: Problem 2Because 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   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 BackgroundFor 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 FilesWe 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 (line numbers are added for reference): (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) 10 10 10 0 0 0
Interpreting OutputLook 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. PWSCF : 0.98s CPU time
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 ...
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. |
 |
|
|
|
