1. The Si lattice constant
In this tutorial the inpgen and fleur input files have to be modified. Please consult the documentation for the inpgen input and the inp.xml file for relevant aspects of these files you are not yet familiar with. You can also have a look at the inpgen example page. Note that the example inputs are far more complex than what is needed in this tutorial.
1.1. Modifying inp.xml for different test lattice constants
The first predictions we want to obtain from the Fleur code are the lattice constant and the bulk modulus of Si. We already created an inp.xml file from a simple inpgen input last time. In this tutorial we will adapt this inp.xml file to perform calculations for different test lattice constants. It is important to adapt the inp.xml file and not the inpgen input since different test lattice constants in the inpgen input will cause a different parametrization in inp.xml. As a consequence the calculations would not be comparable.
The first thing we do is to adjust the maximum number of iterations we are about to perform to 40. This should be enough for obtaining converged results for this system. The XML attribute to be adapted for this is calculationSetup/scfLoop/@itmax.
We want to test lattice constants in 1% steps from 97% to 103% of the lattice constant we initially provided to inpgen. To assure that all of these calculations work we have to guarantee that for none of them the MT spheres overlap. We already saw that inpgen yields valid input. Therefore it is enough to reduce the Si MT radius to 97% of the already available default value. We can calculate the new radius on our own or just write the mathematical expression into to respective XML attribute atomSpecies/species/mtSphere/@radius for the Si-1 species. The radius has to be the same for all calculations.
Next we generate the subdirectories latt-0-97 to latt-1-03 and copy inp.xml and sym.out into each of them. In each directory we can easily adjust a scaling factor in the respective inp.xml file to obtain the desired lattice constant. This scaling factor is found in the XML attribute cell/bulkLattice/@scale. We perform the 7 calculations by calling fleur in each of the subdirectories.
Check whether all calculations are converged, e.g., by looking at the direct Fleur output or by calling "grep dist out". The total energies can be found for each iteration in the out file (grep "total energy" out) or in the out.xml file (grep totalEnergy out.xml). We collect them in a file totalEnergies.txt with two columns: 1. scaling factor, 2. converged total energy in Htr.
Plot totalEnergies.txt with some plotting tool like xmgrace, gnuplot, or Matplotlib.
1.2. The Birch-Murnaghan equation of state
The calculated total energy values should be describable by the Birch-Murnaghan equation of state. As the lattice constant and the bulk modulus are parameters of this equation we fit our data to it to extract these quantities. A program that can be used to fit the data can be found here.
To use the program we copy totalEnergies.txt to totalEnergies-murn.txt and provide some more information in the first lines of the new file. An example for the additional lines is
3 # Hartree 270.107 # unit cell ratio 0.960 1.040 50 # min, max unit cell, no points 7 # no energy val
The first line sets the unit for the provided energy values (3 for Hartree). The second line is the unit cell volume for the scaling factor 1.0. You can obtain it by calling "grep volume latt-1-00/out.xml". The third line controls the output of the program. In this example beyond the desired parameters it will provide 50 energy values obtained from the fitted equation between a scaling factor of 0.96 and 1.04. The fourth line is the number of energy values we provide as input to the fitting program.
We store the Murnaghan fitting program in a file murn.f and compile it with a Fortran compiler like ifort by invoking "ifort -o murn.x murn.f".
We execute the program with "murn.x < totalEnergies-murn.txt".
The program writes out the equilibrium lattice constant in terms of the optimal scaling factor. Convert it to Angstrom or a_0. The bulk modulus is written out in MBar. Convert that to GPa. Finally compare the obtained values to experimental values availaible at wikipedia, WebElements, or other freely available sources (normally we would compare to experimental or theoretical results from scientific papers). Plot the energy values obtained through the fitting procedure together with those obtained from the DFT calculation.
1.3. The lattice constant with LDA
By default inpgen sets up the inp.xml file such that the PBE GGA XC functional is used. We want to compare the results we obtained with this functional with analogous LDA calculations. For this we change the the xcFunctional/@name XML attribute in the inp.xml file to "vwn" and perform the steps described above again for this Vosko-Wilk-Nusair parametrization of the LDA XC functional.
Compare the lattice constant obtained from the PBE calculations with the one from the VWN calculation. Which one is larger and why?
2. Exercise: Some more lattice constants
For many materials lattice constants predicted by DFT calculations with local or semilocal XC functionals match the experimental values with a typical deviation of only a few percent. The bulk modulus typically deviates by 10-20%. Thus, many DFT structure predictions are very good. In the excercises we want to test some limiting cases. For each example make sure that you calculate at least 2 test lattice constants on each side of the total energy minimum.
2.1. fcc vs. hcp Cu
The fcc lattice features an ABCABC stacking of the close-packed atomic monolayers. The hcp lattice is very similar to this but features an ABABAB stacking. This similarity yields tiny energetic differences between the respective structures. We compare these structures for Cu which features an fcc ground-state structure with a lattice constant of 3.61 Angstrom. Calculate the PBE equilibrium nearest-neighbor distance for the atoms in fcc and hcp Cu by performing calculations with different test lattice constants. For hcp Cu we assume the optimal c/a ratio for the packing of spheres. Plot the "energy per atom" vs. "nearest neighbor distance" values for both calculations in one plot. According to your plot: Which is the ground-state structure? How large is the energetic difference per atom in meV?
NOTE1: The choice of the kPoint set for the Brillouin zone integration is an important aspect to obtain comparable results. To obtain comparable kPoint sets we add to the end of the inpgen input for the fcc calculation
&kpt div1=14 div2=14 div3=14 /
and to the one for the hcp calculation
&kpt div1=14 div2=14 div3=7 /
Alternatively we can also modify the specification of the kPoint set in the inp.xml file.
NOTE2: To assure that the calculations are comparable we also have to set some numerical parameters to the same values. These are provided in the XML attributes:
calculationSetup/cutoffs@Kmax calculationSetup/cutoffs@Gmax calculationSetup/cutoffs@GmaxXC atomSpecies/species/mtSphere@radius atomSpecies/species/mtSphere@gridPoints atomSpecies/species/mtSphere@logIncrement atomSpecies/species/atomicCutoffs@lmax atomSpecies/species/atomicCutoffs@lnonsphr
Choose one of the calculations as reference for the values and set the attributes in the other calculations to the same values.
2.2. fcc Ar
An Argon noble gas crystal in fcc structure has an experimental lattice constant of 5.26 Angstrom. Which equilibrium lattice constant do you obtain with the PBE functional?
2.3. Lattice constant of your choice
Choose a chemical element and find out in the literature the ground-state structure or another common structure in which this element can be observed. If the respective unit cell contains more than 2 atoms, the material is magnetic, or describes a molecule crystal choose again. For the final material determine the PBE lattice constant and the bulk modulus. Compare the obtained values to the literature.