Quantum Espresso
Let us use the package "Quantum espresso" to learn what we can do with DFT.
Why QE?
Quantum espresso was called PWscf before. With PWscf, one can do a DFT calculation using the plane wave basis set and the pseudopotential. This software package has been extended by including many subroutines and was renamed as Quantum espresso. There have been continuous efforts to distribute this package to general users and this has become one of the most popular packages. Detailed documents and tutorials are available; from this viewpoint, QE is the most suitable for learners.
Let us learn from the tutorial
Let me acknowledge Prof. Piccinin
Let me acknowledge Prof. Piccinin
QE can be obtained from this site.
Suppose we want to calculate the crystalline silicon of diamond structure.
One need to prepare the input file indicating
-
How to approach the problem
-
Information on the cell
-
Number of atoms and elements
-
Number of basis functions
-
How to achieve self-consistency
and then
-
Element information (mass, pseudopotential)
-
Position of each atom
-
Mesh used for the Brillouine zone integration
QE solves the Schroedinger type equation, called Kohn-Sham equation. Because of the nonlinearlity, it needs to be solved iteratively by assuming a form for the electron density n(r) and correcting the assumption afterwards.
Coming back to the input file, we notice that there are several ways for approach:
-
scf: self-consistently solve the KS equation
-
nscf: non-self-consistently solve it using a guess for n(r)
-
relax: optimize the nuclear position according to the atomic force obtained from the solution of the KS equation.
-
md: run a molecular dynamics simulation using the atomic force
There are some comments on constructing the Colomb potential originated from the nuclear charge.
In QE, we apply the periodic boundary condition unless otherwise specified. So, the atoms in a unit cell is repeated three-dimensionally.
Therefore we need to specify the unit cell. For this, we have to provide the type Bravais lattice (ibrav), length of the cell, and so on.
For silicon, which has a fcc unit cell, we need to specify the Bravais lattice type and the cell dimension a.
We have two silicon atoms in the cell.
We have only one element, "silicon".
We have to specify the species, mass, and the pseudopotential (PP). The PP is specified by the name of the file that contains information on PP.
PP will be explained later.
The PP file is stored in a directory specified by "pseudo_dir"
Position of atoms is given in unit of a. It can be alternatively given in unit of Bohr or Angstrom.
First of all, we need to assume a form for the electron density n(r). It is automatically assumed on the basis of the atomic charge stored in PP file, but you can alternatively give a randomly generated density or specify a file that store the information of n(r).
As written above, we normally use the plane wave as the basis function. We then need to specify the number of the functions. This is done by indicating the threshold value. Normally we use the cut-off energy for the kinetic energy.
Here I comment on the pseudopotential (PP). PP is necessary for the plane wave because the localized nature of the core electrons cannot be described using unrealistically large number of plane waves. So, effect of core is described by introducing an effective potential for valence electrons; core electrons are not treated explicitly. The valence electrons are then made node-less in the core region. PPs are constructed by smoothly extrapolating the potential inside the cut-off radius.
There are various methods for constructing PP and the constructed PPs are available on the internet.
You can construct your PP using the code provided by QE. Here is the input file fot the PP generation calculation.
Here we show a download site for the PP provided by QE.
The next topic is about the k-point mesh used for the Brillouine zone (BZ) integration.
We discretize inside the BZ and apply the trapezoidal rule or the tetrahedron method for the BZ integration.
This figure assumes 4 4 4 1 1 1, differently from the sample input. If you specify 4 4 4 0 0 0, then the k-points include the origin, called gamma point.
We normally use different BZ integration method for metals and insulators. The BZ for metals is demanding. The simplest one is to apply the trapezoidal rule with the Fermi surface smeared by introducing the Gaussian broadening.
We need to specify the functional form for the exchange-correlation. This is put implicitly in the PP file. In the sample input, a semi-local functional of generalized gradient approximation (GGA) type is used. PBE refers to those who developed the functional.
Popular GGA functionals include RPBE, PZ, B3LYP, and so on.
There are a few ways to solve the Shroedinger type differential equation. Davidson method and the conjugate gradient (CG) are the most popular.
If the input density and the output density are the same within the given threshold, self-consistency is achieved and QE proceeds to calculate the total energy and the force.
It is known for a long time that the self-consistency can be generally achieved by mixing the input density with the output density. This parameter depends significantly on the electronic structure. In some cases we need take 0.05!
To relax the atomic structure, we need to specify other parameters as well.
Report (1) deadline is 13th December
The Quantum Espresso (QE) community provides excellent tutorial not only for beginners but also more advanced users. Here we follow the tutorial given by Prof. Madhura Marathe to solve the Kohn-Sham equation of crystalline silicon and plot the band dispersion along lines connecting symmetric points in the Brillouine zone.
Using the obtained plot, let us estimate the effective mass of the electrons and holes. By comparing the result with experiment, let us consider the reason for the discrepancy, which is due to approximations made in the calculation.
Report the results. This is a first step toward more complicated calculations. Try to understand what you are doing.