1. Convergence with respect to parameters
There are several parameters that control the precision of DFT calculation. Some of them depend on the concrete approach used to implement DFT, others are relevant for all approaches. In this tutorial we want to systematically converge our results with respect to these parameters.
1.1. fcc Cu and Kmax
The LAPW basis set size is mainly controlled by . This is a reciprocal plane-wave cutoff parameter used to specify which plane waves are used to represent the wave functions in the interstitial region of the unit cell. Similar parameters are also present in plane-wave pseudopotential or projector-augmented-wave codes. In some codes this parameter is specified in terms of an energy. Assume the kinetic energy operator (here in Hartree atomic units). What kinetic energy does a plane wave with a of have (in Hartree and in Rydberg).
In Fleur we also have the parameters and which are used to limit the expansion of the density and the XC potential in terms of plane waves. A condition the three plane-wave cutoff parameters have to fulfill is .
We want to test for one of our test systems (choose fcc Cu with a lattice constant near the equilibrium lattice constant) the dependence of the total energy on . For this we perform self-convergent total energy calculations with increasing in steps of starting with the default until the total energy seems to be stable. Make sure to set fitting (and same) , for all calculations (assume that over the calculations you might increase by up to ). How large is the total energy difference between the calculations with the default parameters and the one with the largest . Plot ((total energy)-(min(total energies))) vs. .
Note: The cutoffs can be set in the XML element calculationSetup/cutoffs.
Note2: The basis set size depends cubically on . On the other hand for large unit cells the algorithm scales cubically with the number of basis functions. For large unit cells this implies a scaling of the runtime with . However, does not depend on the unit cell size. It does, however, vary from material to material.
1.2. lmax and lnonsph
For the calculation with the converged Kmax we now converge the lmax and lnonsph cutoffs for each atom (if you have multiple atom species in principle these are separate optimizations but in such a case we perform the optimizations simultaneously). The lmax cutoff parameters limit the spherical harmonics expansion of the basis functions in the MT spheres. Beyond this cutoff the basis functions are discontinuous at the respective MT boundary. As such discontinuities are related to large kinetic energy contributions care has to be taken to obtain reasonably smooth basis functions, even if physical effects are related to smaller angular momentum quantum numbers. A rule of thumb for the choice of lmax is , where is the respective MT radius.
Simultaneously with lmax we can also increase lnonsph, though the only connection of this parameter to lmax is . lnonsph has a stronger connection to the actual physics of the material under investigation as it determines up to which of the basis functions the contributions originating from the nonspherical part of the potential in the MT sphere are considered. Typically we choose . For many materials lnonsph can be choosen smaller than lmax but for certain calculations it may be more reasonable to have . Note that increasing lnonsph is computationally much more expensive than increasing lmax.
We perform total energy calculations for increasing lmax, lnonsph values in steps of 2, starting with the default value. Plot ((total energy)-(min(total energies))) vs. .
1.3. k point set and Fermi smearing
All electronic structure calculations for crystals require a k point sampling of the reciprocal unit cell. This is primarily needed for the Brillouin zone integration for the construction of the electron density. Similarly to the previous two tasks we test the convergence of the total energy with respect to the k point set. For this we increase the k point sampling in each direction simultaneously in steps of 1 until the result is stable. Plot ((total energy)-(min(total energies))) vs. numkPointsInOneDirection.
Note that the k point sampling is more demanding for metals than it is for insulators. This is due to the complexity of the electronic structure in the near of the Fermi energy in combination with the occupation of states. For materials with a very complex electronic structure around the Fermi energy a too small k point sampling may even impede the convergence of the DFT calculation, i.e., the iterative procedure does not find a self-consistent solution. One reason for this can be the shifting of states around the Fermi energy from iteration to iteration. This can lead to large changes in the occupation of states. To soften such effects one typically introduces a smooth occupation function. Often this is a Fermi-Dirac distribution, parametrized by an energy or temperature. In Fleur we provide the fermiSmearingEnergy in calculationSetup/bzIntegration/@fermiSmearingEnergy (in Hartree). Alternatively we can also provide a temperature to specify the smearing.
After we reached a self-consistent solution with optimized Kmax, lmax, k point set we reduce the fermiSmearing starting with the self-consistent result for the default smearing. Plot the convergence behavior.
2. Excercises
NOTE: REFERENCE RESULTS FOR LAST WEEKS EXERCISES CAN BE FOUND BELOW. THEY MAY BE NEEDED TO PERFORM THIS WEEKS EXERCISES.
2.1. fcc vs. hcp Cu revisited (CANCELED, you don't have to do this due to time demands!)
As practiced above we also converge the parameters for hcp Cu. When finished decide on a common set of parameters for the fcc and hcp calculations. For fcc optimize the lattice parameter and perform a Birch-Murnaghan fit. How large are the deviations from the previous results. For hcp Cu we also recalculate the results, perform a Birch-Murnaghan fit and compare the results. For the optimized 'a' parameter we then relax the predefined c/a ratio. What is the optimal c? To what extent is the energy difference between fcc and hcp Cu affected by these more sophisticated calculations.
2.2. Si revisited
From your PBE Si calculations use a test lattice constant near the equilibrium lattice constant and converge Kmax, lmax, lnonsph, and the k point set (ignore the Fermi smearing, it is irrelevant for insulators). With the optimized parameters determine the lattice constant and the bulk modulus again. Compare to the previous results.
Note: If the k point set is defined by a kPointCount XML element the provided number does not exactly define the number of k points. It is only a rough suggestion to Fleur how many k points should be used. Different provided numbers can yield the same k point set.