====================== Solution to Exercise 1 ====================== --------------------------- Getting started with Dacapo --------------------------- .. contents:: .. section-numbering:: .. |a0| replace:: *a*\ :sub:`0` .. |angst| unicode:: U+0212B .. ANGSTROM SIGN .. |infin| unicode:: U+0221E .. INFINITY .. |simeq| unicode:: U+02243 .. ASYMPTOTICALLY EQUAL TO .. |sigma| unicode:: U+003C3 .. GREEK SMALL LETTER SIGMA .. |Delta| unicode:: U+00394 .. GREEK CAPITAL LETTER DELTA .. |mu| unicode:: U+003BC .. GREEK SMALL LETTER MU .. |beta| unicode:: U+003B2 .. GREEK SMALL LETTER BETA .. |pi| unicode:: U+003C0 .. GREEK SMALL LETTER PI .. |alpha| unicode:: U+003B1 .. GREEK SMALL LETTER ALPHA .. |nu| unicode:: U+003BD .. GREEK SMALL LETTER NU .. |deg| unicode:: U+000B0 .. DEGREE SIGN Bulk Al ======= Running a calculation --------------------- * The default exchange correlation functional used by Dacapo is the Perdew-Wang 91 GGA functional. * Taking the weighted average over all **k**-points one would expect that one band is fully occupied and one band is partly occupied, as Aluminum has three valence electrons. Of course, the occupation numbers vary over the different **k**-points. Convergence in **k**-points and planewave energy cutoff ------------------------------------------------------- For the convergence with respect to **k**-points, one gets the diagram shown here: .. figure:: conv_kpoints.png :width: 400 :align: center :alt: k-point convergence :figwidth: 400 Convergence of the total energy of bulk fcc Aluminum with respect to the **k**-point sampling From this one can see that a **k**-point sampling of 6 x 6 x 6 (6 **k**-points in each direction of the reciprocal unit cell) should be sufficient. As the last part of the exercise (relative stability of the fcc and the bcc lattice structures) requires a good **k**-point sampling, we will use a sampling of 10 x 10 x10 in the following. For the convergence with respect to the kinetic energy cutoff of the plane waves, one gets this diagram: .. figure:: conv_cutoff.png :width: 400 :align: center :alt: k-point convergence :figwidth: 400 Convergence of the total energy of bulk fcc Aluminum with respect to the kinetic energy cutoff of the plane waves From this we see that a kinetic energy cutoff of 200-250 eV is sufficient for Aluminum. In the following we use a cutoff of 250 eV. * The necessary energy cutoff depends on the softness of the pseudopotential and is therefore the same for all Aluminium structures. For other chemical elements, other pseudopotentials with different properties are employed, and therefore, a different energy cutoff may be required. For example, transition metals have often pseudopotentials that require an energy cutoff of 25-30 Ry = 340 - 400 eV. The necessary **k**-point sampling depends on the complexity of the Fermi surface and the bandstructure. This can be assumed to be similar for different Aluminum bulk structures, and therefore, the same **k**-point sampling can be used. For other chemical elements, the required **k**-point sampling is different. In general, metals require a higher **k**-point sampling than semiconductors or insulators. * The energy cutoff determines the number of basis functions, which represent the density and the Kohn-Sham wavefunctions. Therefore, the energy cutoff is subject to the variational energy principles, and the total energy of a system has to decrease if the cutoff energy is increased. The necessary **k**-point sampling depends on the complexity of the Fermi surface. Therefore, one cannot predict, whether the total energy decreases or increases for a finer **k**-point sampling. * Al is an easy element to work with, because its pseudopotential is softer than many other metals. It is not magnetic, and therefore spin-polarized calculations are not required. Furthermore, its electronic structure resembles the free electron gas, and therefore convergence is easy and fast. Equilibrium lattice properties ============================== First, we map out the cohesive curve *E*\ (*a*) for Al(fcc), i.e. the total energy as function of lattice constant *a*, around the experimental equilibrium value of |a0| = 4.05 |angst|. We obtain the curve shown in the following figure, which we fit to a third-order polynomial. .. figure:: fcc_lattice.png :width: 550 :align: center :alt: energy versus lattice constant :figwidth: 550 Total energy versus the lattice constant for bulk fcc Aluminum. The solid line is a third-order polynomial fit. The next step is to take the first and the second derivative of the third order polynomial fit. Now we can easily determine the equilibrium lattice constant |a0| by reading off, at which lattice constant the first derivative is zero. Then the value of the second derivative at |a0| has to be read off. The results are: |a0| = 4.05 |angst| and *d*\ :sup:`2`\ *E*/\ *da*\ :sup:`2` = 3.84 eV/\ |angst|\ :sup:`2` For calculating the bulk modulus using the expression *B* = *M*/(9\ |a0|\ )\ *d*\ :sup:`2`\ *E*/*da*\ :sup:`2` we note that *M*, which is the number of atoms per cubic unit cell, is *M* = 4 for an fcc lattice. Remembering that 1 GPa = 10\ :sup:`9` J/m\ :sup:`3` we have to convert the second derivation from eV/\ |angst|\ :sup:`2` to J/m\ :sup:`2`. This results in: *d*\ :sup:`2`\ *E*/\ *da*\ :sup:`2` = 3.84 eV/\ |angst|\ :sup:`2` = 3.84 x 1.602 x 10\ :sup:`-19` J / (10\ :su p:`-10` m)\ :sup:`2` = 61.52 J/m\ :sup:`2` Keeping in mind that *giga* corresponds to the prefactor 10\ :sup:`9`, we obtain for the bulk modulus: *B* = 4/(9 x 4.05 x 10\ :sup:`-10` m) x 61.52 J/m\ :sup:`2` = 67.51 GPa These results compare well to the experimental values |a0| = 4.05 |angst| and *B* = 76 GPa. The bulk modulus compares not as well as the lattice constant. As it depends on the second derivative of the fit, it is a rather sensitive quantity. Equilibrium lattice properties for bcc ====================================== First, we try to find a qualified starting guess for *a*\ :sub:`bcc` from the value of *a*\ :sub:`fcc`. One possibility is to assume that the primitive unit cell volumes of the fcc and bcc structure are equal. The unit cell volume for the cubic unit cell of fcc is *a*\ :sub:`fcc`\ :sup:`3` and there are four atoms in the cubic unit cell. Therefore the unit cell volume of the fcc primitive unit cell is *a*\ :sub:`fcc`\ :sup:`3`/4. In the same way, one obtains the unit cell volume for the bcc primitive unit cell, which is *a*\ :sub:`bcc`\ :sup:`3`/2. Thus *a*\ :sub:`bcc` = *a*\ :sub:`fcc` 2\ :sup:`-1/3` = 3.21 |angst| Another possibility is to assume that the nearest neighbor distances in the fcc and the bcc structures are equal. The nearest neighbor distance in fcc is half a face diagonal, *a*\ :sub:`fcc`/2\ :sup:`1/2`. The nearest neighbor distance in bcc is half a body diagonal, (3\ :sup:`1/2`/2)\ *a*\ :sub:`bcc`. From this we obtain *a*\ :sub:`bcc` = (2/3)\ :sup:`1/2`\ *a*\ :sub:`fcc` = 3.31 |angst| As the lattice constant for the bcc lattice is smaller, one should use a higher **k**-point sampling to keep the product of lattice constant and number of **k**-points constant. The cohesive energy curve for the bcc Aluminum structure is shown here: .. figure:: bcc_lattice.png :width: 450 :align: center :alt: energy versus lattice constant :figwidth: 450 Total energy versus the lattice constant for bulk bcc Aluminum. The solid line is a third-order polynomial fit. From a third order polynomial fit, we obtain: |a0| = 3.25 |angst| and *d*\ :sup:`2`\ *E*/\ *da*\ :sup:`2` = 5.56 eV/\ |angst|\ :sup:`2` This gives a bulk modulus of *B* = 60.9 GPa The bcc bulk modulus is very similar to the fcc one. One would expect this, as both structures have approximately the same density, and therefore one would expect the elastic properties to be fairly similar. Now, we calculate the fcc/bcc total energies at their respective lattice constants and for two different plane wave cutoffs: 200 eV and 450 eV. We obtain these values: .. |Efcc| replace:: *E*\ :sub:`fcc` .. |Ebcc| replace:: *E*\ :sub:`bcc` =========== =========== =========== =============== cutoff [eV] |Ebcc| [eV] |Efcc| [eV] |Ebcc| - |Efcc| =========== =========== =========== =============== 200 -56.342 -56.448 0.107 450 -56.398 -56.509 0.111 =========== =========== =========== =============== One can see that the energy difference between the bcc and the fcc structure is well converged at 200 eV, whereas the total energy is not yet fully converged. In general, energy differences converge much faster than total energies. One can use that fact to reduce the cutoff energy, because convergence of energy differences is all that is needed. The absolute total energy values depend on the pseudopotentials and do not have any physical significance.