1. Band structures, semicore states, and local orbitals

1.1. Silicon band structure

Perform a self consistency calculation of Si with equilibrium lattice constant (or another test lattic constant). Then modify the inp.xml file in the following way:

1. Set output/@band to "T". This enables the band structure calculation mode in fleur. With this only a single iteration is performed and no new charge density is generated. Instead some files related to the band structure of the material are written.

2. For the calculation of a band structure a special $\vec{k}$-point set with points along the desired $\vec{k}$-point path has to be used. Such a $\vec{k}$-point set can only be specified within a 'kPointCount' XML element. The high symmetry path is either defined explicitly or a default $\vec{k}$-point path is generated. By default in the inp.xml file such a $\vec{k}$-point set is specified in the optional 'altKPointSet' XML element if the respective XML attribute 'purpose' is set to 'bands'. Check whether such an entry is present in your inp.xml file. If not replace the specification of the $\vec{k}$-point set by

   <kPointCount count="240" gamma="F"/>

   With this you specify that 240 $\vec{k}$ points along the high symmetry path of the band structure are calculated.

3. Set calculationSetup/cutoffs/@numbands to 25. If numbands is 0 a default number of states is calculated for each $\vec{k}$ point. This number is adapted to the charge density construction. Setting the parameter to 25 explicitly specifies that we want to calculate 25 states at each $\vec{k}$ point. With this we make sure that we also obtain a reasonable number of bands above the Fermi energy.

Call fleur on top of this calculation with the new inp.xml file. This will produce a band structure ouput, i.e., a band.gnu file with a gnuplot script to visualize the band structure and a bands.1 file with a list of eigenenergies for each $\vec{k}$ point. Use the gnuplot script and convert the result into a pdf with

gnuplot < band.gnu > bands.ps
ps2pdf bands.ps

Visualize the pdf with some viewer (e.g. evince). Where is the highest occupied state, where is the lowest unoccupied state? How large is the bandgap (grep for bandgap in the out file)? Compare the result with the experimental value? How large is the difference and why is there a difference?

1.2. Strontium band structure

Use this inp.xml file for Strontium together with the sym.out file for fcc Cu to calculate the band structure for Strontium! Does the result look correct to you? Can it be improved? If so: How? Improve it!

Note: To change the description of a semicore state from a core electron description to a valence electron description you have to perform the following steps:

1. Identify the semicore state, i.e., the main quantum number and the angular momentum quantum number. This is best done by taking a look into the list of "coreStates" in the out.xml file. Semicore states feature energies that are near the valence band. Depending on the angular momentum of the state this may actually be 2 states as the fully-relativistic description for the core electrons yields a splitting of the states due to spin-orbit coupling. The list of core states also provides information on how many electrons there are in this state. This is the "weight".

2. For each atom (not atom type, group, or species) for which the description of the semicore state has to be changed add the number of related electrons to calculationSetup/bzIntegration@valenceElectrons.

3. Add the respective local orbital for the state (type is SCLO; n, l as identified; eDeriv 0) to the related atom species section just below the energy parameters.

4. Reduce the number of core states for the respective atom species by the number of states identified in (1.). This is 1 or 2. The XML attribute to be changed is atomSpecies/species/@coreStates.

2. Exercises

2.1. Proof of orthogonality to core states

Proof $\langle u_l^c(\epsilon_c,r) \vert \dot{u}_l(E_l,r) \rangle = 0$ under the assumptions $u_l^c(\epsilon_c,r) \vert_{R_\text{MT}} = 0$ and $\frac{\partial}{\partial r} u_l^c(\epsilon_c,r) \vert_{R_\text{MT}} = 0$. The function $u_l^c(\epsilon_c,r)$ is the wave function of a core state with eigenenergy $\epsilon_c$. $\dot{u}_l(E_l,r)$ is the energy derivative of $u_l(E_l,r)$ and thus part of the LAPW basis.

Hint 1: $E_l$ is an energy parameter and not an eigenenergy for the Hamilton operator. $u_l(E_l,r)$ and $\dot{u}_l(E_l,r)$ are not eigenfunctions.

Hint 2: The differential equation to calculate $\dot{u}_l(E_l,r)$ is obtained by deriving $\hat{H} u_l(E_l,r) = E_l u_l(E_l,r)$ with respect to $E_l$. The non-relativistic Hamilton operator has no energy dependence.

Hint 3: You may assume that $\langle u_l^c(\epsilon_c,r) \vert u_l(E_l,r) \rangle = 0$ has already been shown.

2.2. Ne band structure

Note: In the following you have to modify the gnuplot script band.gnu such that bands up to 40 eV above the Fermi energy are displayed. The simplest way to do this is to replace in the line that starts with "plot" the last bracket by "[-10:40]" and the line before that by "set ytics -10,5,40". However the plot will not look beautiful with this single adaption.

1. Set up an input for a Ne fcc crystal with a lattice constant of 442.9 pm. Calculate the self-consistent solution and calculate the band structure up to energies 40 eV above the Fermi level.

2. Copy the input from (1.) into a new directory and modify it in the following way: The difference to the reference calculation is supposed to be the addition of a few local orbitals constructed from the second energy derivative. Add the following line after the energyParameters xml element for both atom species:

   <lo type="SCLO" l="0-3" n="2,2,3,4" eDeriv="2"/>

   Calculate the self-consistent solution and calculate the band structure up to energies 40 eV above the Fermi level.

3. Copy the input from (1.) into a new directory and modify it in the following way: The difference to the reference calculation is supposed to be the addition of a different set of local orbitals constructed from radial functions evaluated at higher energy parameters. Add the following line after the energyParameters xml element for both atom species:

   <lo type="HELO" l="0-3" n="3,3,4,5" eDeriv="0"/>

   Calculate the self-consistent solution and calculate the band structure up to energies 40 eV above the Fermi level.

How do the results differ? Grep for "EP" in the out.xml file. You obtain a list of the energy parameters with energies relative to the 0 of the potential. Grep for "Fermi" in the out.xml file to obtain the Fermi energy in each iteration, also relative to the 0 of the potential. How far away from the conduction band are the $s$ and $p$ energy parameters.

3. Last weeks results

3.1. Cu

After performing the convergence tests we assume the following optimized parameters:

• $K_\text{max} = 5.5~a_0^{-1}$
• $G_\text{max} = 16.2~a_0^{-1}$
• $K_\text{max}^\text{XC} = 14.5~a_0^{-1}$
• $l_\text{max} = 10$
• $l_\text{max} = 8$
• k point mesh (fcc): 26 x 26 x 26
• k point mesh (hcp): 26 x 26 x 13
• lower Fermi smearing energy (single shot, no optimization): 0.0001 Htr

	default parameters	refined parameters	refined parameters, low smearing	refined parameters, low smearing, c/a optimized
$a_\text{fcc} (a_0)$	$1.000764 \cdot 6.82$	$1.005883 \cdot 6.82$	$1.005896 \cdot 6.82$
bulk modulus (fcc, GPa)	$142.43$	$140.52$	$140.50$
total energy minimum per atom (fcc, Htr)	$−1655.030758$	$−1655.039326$	$−1655.039341$	$−1655.039341$
$a_\text{hcp} (a_0)$	$1.001666 \cdot 4.82247$	$1.006336 \cdot 4.82247$	$1.006345 \cdot 4.82247$	$1.006345 \cdot 4.82247$
$c_\text{hcp} (a_0)$	$1.001666 \cdot 7.87506$	$1.006336 \cdot 7.87506$	$1.006345 \cdot 7.87506$	${} \approx 1.003 * 1.006345 \cdot 7.87506$
bulk modulus (hcp, GPa)	$145.76$	$139.15$	$139.25$
total energy minimum per atom (hcp, Htr)	$−1655.0304925$	$−1655.039014$	$−1655.0390295$	${} \approx −1655.039034$
hcp-fcc energy difference per atom (meV)	$7.22$	$8.49$	$8.48$	$8.35$

Note 1: Even though the total energy difference per atom between the default parameter calculation and the one with optimized parameters is about $8.5~\text{mHtr} \approx 231~\text{meV}$ the energy difference between the hcp and fcc structure for the default parameters differs by less than $1.5~\text{meV}$ from the optimized calculation. Also the differences for the lattice constant and the bulk modulus is minimal. Lesson learned: total energy differences typically converge much faster than total energies. This implies that comparing total energies from different calculations is only meaningful if cutoff parameters and other numerical parameters are identical. Otherwise one might see the effect of these different parameters and not the total energy changes due to the intended differences in the calculation.

Note 2: On the reference computer an iteration of the default parameter calculation takes about $1~\text{s}$. It takes about $26~\text{s}$ for the optimized parameter calculation. Considering more complex unit cells converging the total energy can be very expensive with respect to the computational demands.