RPA for solids
The Random-Phase Approximation (RPA) is an approach for computing the electronic correlation energy. The RPA correlation energy + exact exchange energy is often referred to as the reference DFT method for total energy calculations and is on the fifth rung of the density functional approximation ladder (Jacob's ladder; 1: LDA, 2: GGA, 3: metGGA, 4: Hybrid functionals, 5: RPA + exact exchange, ...).
Let us motivate why this may be of importance. We performed a stability analysis of rutile-TiO2. In total, we performed four structure optimizations: a symmetry-constrained and symmetry-unconstrained relaxation for PBE and PBEsol. For PBE, we found the following minimum on the potential energy surface:
lattice_vector 4.65014962 0.00000000 0.00000000
lattice_vector -0.00000000 4.65014059 0.00000000
lattice_vector -0.00000000 -0.00000000 2.98377469
atom_frac 0.30471240 0.30472314 0.04256115 O
atom_frac 0.69528758 0.69527687 0.04256115 O
atom_frac 0.19528759 0.80472314 0.54256115 O
atom_frac 0.80471239 0.19527688 0.54256115 O
atom_frac -0.00000001 0.00000001 0.00245810 Ti
atom_frac 0.49999999 0.50000000 0.50245852 Ti
and for PBEsol:
lattice_vector 4.59011462 -0.00000000 0.00000000
lattice_vector -0.00000000 4.59011462 0.00000000
lattice_vector 0.00000000 -0.00000000 2.93865021
atom_frac 0.30447192 0.30447193 0.00000000 O
atom_frac 0.69552808 0.69552807 0.00000000 O
atom_frac 0.19552808 0.80447193 0.50000000 O
atom_frac 0.80447192 0.19552807 0.50000000 O
atom_frac -0.00000000 0.00000000 -0.00000000 Ti
atom_frac 0.50000000 0.50000000 0.50000000 Ti
Despite the slightly different volume of 64.520 Å\(^3\) and 61.915 Å\(^3\) for PBE and PBEsol, respectively, we also find that the symmetry is broken for PBE. You can see it best from the following figure. The central O-Ti-O axis has an angle of 173.17\(^{\circ}\) instead of 180\(^{\circ}\) in case of the correct rutile structure. From the experiment, the answer is clear that the correct structure should be the rutile-TiO2. However, from theory the correct answer cannot be predicted by using only GGA functionals.
By computing the RPA correlation energy, we will determine which GGA functional predicts the correct structure.
Requesting the RPA correlation energy
The RPA correlation energy is available as a total energy method and is computed once the SCF cycle has finished.
RPA calculations are very expensive!!!
The calculation of the RPA correlation energy is computationally very expensive. We do not recommend to repeat the calculations as you'll need for some data point more than 1000 CPUs. The idea of this chapter is to introduce you to the main concepts for periodic RPA calculations. You will learn the most important things by just reading this chapter. You can use the output files that we have computed for you from here: https://gitlab.com/FHI-aims-club/tutorials/rpa-and-gw-for-molecules-and-solids/-/tree/main/Tutorials/Part-2/solutions/RPA-calculations
There are two aspects that have to be covered for a periodic RPA calculation (this applies to probably all beyond-DFT methods):
- Convergence of the Observable of interest with the number of k-points.
- Convergence of the Observable of interest with the number of basis functions.
Both aspects will be discussed in the subsequent subsections. Our Observable of interest will be the energy difference of the perfect and distorted rutile structure.
To address these two points, we created the following directory tree:
RPA-calculations
├── distorted
│ ├── light
│ │ ├── k_grid_2x2x3
│ │ ├── k_grid_4x4x6
│ │ └── k_grid_6x6x9
│ ├── intermediate
│ │ ├── k_grid_2x2x3
│ │ └── k_grid_4x4x6
│ ├── tight
│ │ ├── k_grid_2x2x3
│ │ └── k_grid_4x4x6
│ └── tier2
│ └── k_grid_2x2x3
└── perfect
├── light
│ ├── k_grid_2x2x3
│ ├── k_grid_4x4x6
│ └── k_grid_6x6x9
├── intermediate
│ ├── k_grid_2x2x3
│ └── k_grid_4x4x6
├── tight
│ ├── k_grid_2x2x3
│ └── k_grid_4x4x6
└── tier2
└── k_grid_2x2x3
and copied the corresponding geometry.in file into the corresponding subfolder and adapted the following control.in snippet:
xc pbe
relativistic atomic_zora scalar
k_grid <k_x> <k_y> <k_y>
total_energy_method rpa
frozen_core_postscf 2
[attach the corresponding species defaults]
The calculation of the RPA total energy is activated by the keyword total_energy_method rpa. We also use the frozen-core approximation by setting the keyword frozen_core_postscf. The line frozen_core_postscf 2 means that only 2 highest occupied shells are considered. The remaining shells will get frozen. This especially means that for the oxygen atom all electrons are considered and for the Ti atom only the shells 3 and 4.
For the tier2 folders we used the tight species defaults and uncommented the second tier for the Ti species.
We carried out all the calculations and you can find the output files here: https://gitlab.com/FHI-aims-club/tutorials/rpa-and-gw-for-molecules-and-solids/-/tree/main/Tutorials/Part-2/solutions/RPA-calculations
You can download all of the solutions by simply cloning the whole repository by:
git clone git@gitlab.com:FHI-aims-club/tutorials/rpa-and-gw-for-molecules-and-solids.git
Convergence with number of k-points
Let us first analyze the convergence of the energy difference with the number of k-points. In the main output file, the RPA energy can be found after the SCF cycle has finished. The important energies are summarized in the following output block:
--------------------------------------------------------------------
Exact exchange energy : -115.42030152 Ha, -3140.74620389 eV
DFT/HF total energy : -2013.89413206 Ha, -54800.84756990 eV
Exchange-only total energy : -2010.69376830 Ha, -54713.76124086 eV
RPA total energy : -2012.93989759 Ha, -54774.88152873 eV
RPA+SE total energy : -2013.64773347 Ha, -54794.14272298 eV
RPA+rSE total energy : -2013.19202676 Ha, -54781.74231241 eV
------------------------------------------------------------------------------
The energies we are looking for are DFT/HF total energy and RPA total energy.
You can use the script RPA_energies_k_convergence.py to just extract the relevant RPA and DFT energies from the main output file (aims.out). Feel free to alter the script to check the k-point convergence for the other species defaults. The energy difference should be simply calculated by:
Let us look at the k-point convergence for the light species defaults, which are shown in the picture below. Even though the RPA total energies of the individual structure are not fully converged (right side), the total energy difference appears to be converged using a 4x4x6 k-point grid (left side).
Luckily, we can already find the desired effect of the RPA approach: Indeed, the perfect rutile structure is the more stable structure and the symmetry breaking is only an artifact of the PBE XC functional approximation. In the figure below, positive numbers indicate that the distorted geometry (as predicted by PBE) is more stable and negative numbers indicate that the perfect geometry is more stable (as predicted by PBEsol).
Convergence with the number of basis functions
The k-point convergence turned out to be a rather affordable task. In contrast, the convergence with the number of basis functions is by far more difficult. There are two reasons for that:
- The computational load with the number of basis functions increases rapidly and you will need soon thousands of CPU to run the calculation.
- Systematical convergence with numeric atom-centered Orbitals (NAOs) of the desired property is usually hard to achieve.
To analyze the basis set convergence you can use the script RPA_energies_basis_convergence.py, which will produce the below plot.
The figure shows several curves among which the blue corresponds to the convergence of energy difference for the 2x2x3 k-point grid and the green to the 4x4x6 k-point grid, each using the species default light, intermediate, tight, tier2. Tier2 refers to tight+additional basis function of tier 2 from the Ti species.
The data points clearly show that the perfect structure is stable for all used species defaults. Even though the convergence is not linear, the trend seems clearly to indicate that by adding more and more basis functions, eventually, the perfect rutile structure becomes more and more stable.
However, we also find that the energy difference hardly converges and we cannot really give a final number up to the accuracy of tier2 species defaults.
As conclusion we find that the structure predicted by PBEsol is in agreement with the RPA calculations.