This is a collection of raw input and output files that were used to produce the following publication:
Eric Macke, Iurii Timrov, Nicola Marzari, and Lucio Colombi Ciacchi
"Orbital-resolved DFT+U for molecules and solids"
Calculations were performed using the open-source Quantum ESPRESSO distribution, version 7.0,
which can be downloaded from www.quantum-espresso.org. The capability to perform orbital-resolved DFT+U calculations
was implemented by us, and will be made available in a future release. Note that the input format of Hubbard
parameters has changed for versions of QE >=7.1.
The following codes of Quantum ESPRESSO were used:
pw.x - the code that performs self-consistent-field (SCF) ground-state calculations
(to obtain total energy, forces, stress, and other properties) using
Hubbard-corrected density-functional theory (DFT+U and DFT+U+V); (Usage : pw.x < scf.in > scf.out)
hp.x - the code that computes the Hubbard U and V parameters using density-functional
perturbation theory (DFPT); (Usage: hp.x < hp.in > hp.out)
projwfc.x - the code to compute the projected density of states (PDOS) (Usage : projwfc.x < pdos.in > pdos.out).
pp.x - the code for postprocessing the charge density file, used here for plotting KS orbitals
(Usage : pp.x < pp.in > pp.out)
sumpdos.x - a tool to sum up over PDOS output files of projwfc.x
(Usage : sumpdos.x *\(Fe\)*\(d\) > Fe_d.dat).
_____________________________
Description of the contents:
_____________________________
The directory structure of the data folder can be visualized as follows:
-- DATA_FOLDER
-- PSEUDOS
-- CORRECTION_TYPE
-- COMPOUND
[IF CORRECTION_TYPE == DFT-orbitalResolvedU-LR-cDFT]
-- HUBBARD_MANIFOLD
-- CALCFOLDERS
[ELSE]
-- CALCFOLDERS
The systems in this work are studied with up to five types of Hubbard corrections. Aside from the single folder PSEUDOS that contains all
pseudopotentials needed to carry out the calculations (SSSP precision v. 1.1.2), the CORRECTION_TYPE forms the first level of subdirectories:
1) DFT-UNCORRECTED - Uncorrected (bare) PBE/PBEsol calculations
2) DFT-U-EMPICIRAL - DFT+U calculations with empirical U values
3) DFT-U-DFPT - DFT+U calculations with U values computed using DFPT (hp.x)
4) DFT-UV-DFPT - DFT+U+V calculations with U and V values computed using DFPT (hp.x)
5) DFT-orbitalResolvedU-LR-cDFT - Orbital-resolved DFT+U calculations with orbital-resolved U parameters computed using LR-cDFT
Inside these folders, there is another level of subdirectories given by the COMPOUND.
Note that not all compounds are studied with all types of corrections.
In particular, we study the bulk solids pyrite (FeS2) and pyrolusite (β-MnO2), as well as the six Fe(II) molecular hexacomplexes Fe[H2O]6,
Fe[NH3]6, Fe[NCH]6, Fe[PH3]6, Fe[CO]6 and Fe[CNH]6 in their +2 charge state. To calculate the adiabatic spin energies,
the Fe(II) molecular hexacomplexes are computed in both their low-spin states (LS, 0.0 µB total magnetization) and in their
high-spin states (HS, 4.0 µB total magnetization), respectively.
If the correction type is orbital-resolved DFT+U, the next level of directories is the HUBBARD_MANIFOLD, for example ``t2g''. Otherwise, the next
level is directly comprised of the CALCFOLDERS, which contain the input and output files explained below. For the DFPT and LR-cDFT corrections,
some of the calcfolders are named "iter_N", where N are integers starting from 1. ``N'' refers to the self-consistency cycle required to converge
Hubbard parameters and structures (see manuscript).
Finally, in each of these CALCFOLDERS the following files can be found:
- *.relax.in - input files for (vc-)relax calculations using pw.x
- *.relax.out - output files for (vc-)relax calculation
- *.scf.in - input files for SCF ground-state calculations using pw.x
- *.scf.out - output files for SCF ground-state calculations
- *.hp.in - input file for hp.x
- *.hp.out - output file from hp.x
- *.Hubbard_parameters.dat - output file from hp.x containing Hubbard U (and V)
- Hubbard_parameters.dat - output file containing orbital-resolved Hubbard U parameters
- parameters.in - input file containing Hubbard U and V (DFT+U+V only)
- parameters.out - output file containing Hubbard U and V (DFT+U+V only)
- *.pdos.in - input file for the PDOS calculation using projwfc.x
- *.pdos.out - output file from the PDOS calculation using projwfc.x
- *.pdos_tot - output file containing the DOS and summed PDOS
- *.pdos_atm#*(*)_wfc#*(*) - output files from projwfc.x with the PDOS for each orbital of each atom
- *.pdos_*_* - output files of sumpdos.x, used to sum up over the PDOS of symmetry equivalent atoms
- *.band.in - input files for bands calculations using pw.x
- *.band.out - output files for bands calculations using pw.x
- *.pp.in - input file for isosurface calculations using pp.x
- *.pp.out - output file for isosurface calculations using pp.x
- *.cube - output file containing the contribution of selected wavefunction(s) to the (pseudo-)charge density.
_________________________________________________
On the practical evaluation of orbital-resolved
Hubbard U parameters from LR-cDFT
_________________________________________________
The folder ``DFT-orbitalResolvedU-LR-cDFT'' contains a subdirectory ``EXAMPLES'' where we demonstrate how
orbital-resolved Hubbard U parameters can be evaluated by applying linear perturbations to the Hubbard manifold.
CAVEAT: The formatting of pw.x's input and output files is specific to the in-house version of QE that was used to produce
this work. Future (publicly avilable) releases of QE featuring orbital-resolved Hubbard U are likely to use different
(and more user friendly) formats. Thus, the input files and scripts provided in the ``EXAMPLE'' folder might require
slight adaptions in order to function properly with future versions of QE.
MANUAL EVALULATION:
EXAMPLE_1 shows the manual calculation of the orbital-resolved U_{t2g} parameter of low-spin Fe[CO]6. To compute U, we start from
an unperturbed SCF calculation (feco6.scf.in /.out) with an intial U parameter of 3.5eV that serves as a starting guess.
Then, we perform two perturbative calculations in which we apply linear perturbations alpha of strengths +0.05eV and -0.05eV to the eigenstates
representing the t2g manifold (3,4, and 5 of Fe-d). The input and output files are named feco6.05.1.in /out and feco6.min05.1.in /out, respectively.
Since this system contains only one Hubbard atom and only one Hubbard manifold (Fe-t2g), we can easily extract the
unperturbed and the perturbed responses of the t2g manifold by looking into the output files.
For convenience, pw.x computes the atomwise sum over the eigenvalues of all eigenstates with nonzero Hubbard alpha values and
outputs it starting with the ``@'' symbol after every SCF iteration.
Thus, to obtain the unscreened occupation numbers of t2g (n0) that are obtained after the first SCF iteration, it suffices to type:
grep @ feco6.min05.1.out | head -n 1
grep @ feco6.05.1.out | head -n 1
This yields the values 5.20273873 for the -0.05eV perturbation and 5.18051626 for the +0.05eV perturbation.
The unscreened response chi0 equals the slope Δn0/Δalpha, i.e.:
chi_0= (5.18051626 - 5.20273873) / (0.05eV - (-0.05eV)) = −0.2222247 ev^(-1)
In the same fashion, we can obtain the screened occupation numbers following the perturbation (n). Since these are
found when the SCF calulation converges, we extract them using the ``tail'' utility:
grep @ feco6.min05.1.out | tail -n 1
grep @ feco6.05.1.out | tail -n 1
We obtain 5.19725870 (+0.05eV) and 5.18615574 (-0.05eV) and proceed to compute the screened response chi:
chi= (5.18615574 - 5.19725870) / (0.05eV - (-0.05eV)) = −0.111029 ev^(-1)
In the final step, we invert both numbers to obtain U_{t2g}:
U_{t2g} = [chi0_{t2g}^(-1) - chi_{t2g}^(-1)] = (1/−0.2222247) - (1/−0.111029) = 4.51eV
SCRIPT-BASED EVALUATION:
The manual calculation shown in EXAMPLE_1 is good for illustrative purposes and suffices for simple molecules with a single Hubbard atom.
However, working with the responses of multiple Hubbard manifolds in supercells of bulk solids can become quite tedious.
In EXAMPLE_2, a script is used to extract the responses of Mn-t2g and O-p_z from 48 atom (2x2x2) supercell calculations
of bulk β-MnO2. The output files of the perturbed runs are:
mno2.05.1.out ---> +0.05eV perturbation of t2g in atom 1 (Mn)
mno2.min05.1.out ---> -0.05eV perturbation of t2g in atom 1 (Mn)
mno2.05.17.out ---> +0.05eV perturbation of pz in atom 17 (O)
mno2.min05.17.out ---> -0.05eV perturbation of pz in atom 17 (O)
The script can be used by typing (in the console, with bash available):
bash extractor_script.bash
When executed,``extractor_script.bash'' runs over the output files fechting the unscreened and screened occupation numbers of the perturbed states
and writes them into files named dn0.RESPONSE_ATOM.da.PERTURBED_ATOM.dat for the unscreened case and dn.RESPONSE_ATOM.da.PERTURBED_ATOM.dat
for the screened case, respectively. For example, dn.34.da.1.dat contains the screened responses of the p_z manifold of atom 34 (oxygen) to
the positive and negative perturbations to the t2g manifold of atom 1 (manganese).
Moreover, ``extractor_script.bash'' outputs a file containing a list of all produced data files called ``file.mno2''. This file is needed to run
the small program ``resp_mat.f90'', which performs the inversion of the chi0 and chi response matrices. ``resp_mat.f90'' was obtained from
https://doi.org/10.24435/materialscloud:vp-wm and was used to generate the data of Phys. Rev. B 103, 045141 (2021).
For the final step, we first need to compile resp_mat.f90, which requires a working Fortran compiler as well as BLAS and LAPACK libraries.
You can download LAPACK (and BLAS, included) from the official GitHub Repository, currently at https://github.com/Reference-LAPACK/lapack/releases .
Asuming the all requirements are fulfiled, and using the GNU Fortran compiler, we can compile the utility by typing:
gfortran -o resp_mat.x resp_mat.f90 -lblas -llapack
The resulting executable ``resp_mat.x'' needs to be made executable:
chmod +x resp_mat.x
Finally, it can be called and provided with an input file:
resp_mat.x < resp_mat.in
``resp_mat.x'' will read this input file, the response data files (dn.*) and a file containing information on the geometry of the system (``pos_mno2''),
assemble the response matrices chi0 and chi and finally invert them to evaluate the Hubbard U parameters, which can be found at the end of the output
file ``Hubbard_parameters.dat''. We find U_{t2g} = 1.77eV and U_{p-z}=4.47eV.
GENERAL RECOMENDATIONSW WHEN COMPUTING (ORBITAL-RESOLVED) HUBBARD PARAMETERS:
1. Perturbative calculations (any calculation with nonzero Hubbard alpha values) MUST be restarted from a converged,
unperturbed SCF charge density (startingpot = 'file').
2. Be careful NOT to overwrite the converged unperturbed charge density while performing the perturbative runs. This is best achieved
by setting disk_io='none' in the &CONTROL section.
3. To obtain meaningful unscreened occupation numbers (n0), set the convergence threshold for iterative diagonalization to tight values, e.g., 1.D-9 or tighter.
This is done using the ``diago_thr_init'' keyword in the &ELECTRONS section.
4. During perturbative calculations, the Hubbard potential should be kept fixed to its self-consistent (unperturbed) value to ensure that U is evaluated
using only the DFT part of the curvature. This is achived using the currently undocumented pw.x input keyword ``hub_pot_fix=.true.'' in the &SYSTEM section.
Note that these points will be simplified in the publicly avilable release of QE that supports orbital-resolved Hubbard U calculations.