Computing Hubbard parameters#
In this tutorial you will learn how to calculate the Hubbard parameters step by step using aiida-hubbard.
We can divide this goal in three phases:
Define the manifolds: define the target Hubbard manifolds via the
HubbardStructureDataSCF ground-state: calculate the ground-state using the
PwBaseWorkChainDFPT calculation: use the
HpBaseWorkChainto do a self-consistent perturbation calculation to predict the Hubbard parameters.
In this tutorial we will make use of the silicon structure to give you an overall understanding of the usage of the package. If you are interested in more advanced features, please have a look at the next tutorial or to the how tos.
Let’s get started!
Defining the target manifold through the HubbardStructureData#
The Hubbard correction is a corrective term that is added to the Hamiltonian of a system which suffers from great self-interaction errors. This is usually the case for transition metals on their d manifolds. An extra correction to account for the hybridization can be accounted for with the ligands, typically belonging to the p element group. Such interaction needs to be localized in space. This is the reason why we need to define the projectors. Quantum ESPRESSO allows you to define different type of projections \(| \phi^I_m \rangle\) (\(m\) orbital quantum number, \(I\) atom in cell). Currently, the ortho-atomic projectors are the most accurate ones implemented.
Still, we need to ask the program on which atoms \(I\) and which manifolds \(m\) to project and correct for this self-interaction.
Since manifolds and atoms belong to the structure, then you need to definet them together as an HubbardStructureData.
In the following, we take LiCoO2 as example, and we suppose we want to target the 3d orbitals of cobalt and the intersite interaction between 2p of oxygen and 3d of cobalt.
Note
By default we set ortho-atomic projectors and we use the Dudarev formulation.
Let’s define the HubbardStructureData:
Note
If you already have a aiida.orm.StructureData, you can load the structure information in HubbardStructureData as follows:
my_structure = load_node(IDENTIFIER)
hubbard_structure = HubbardStructureData.from_structure(my_structure)
from aiida_quantumespresso.data.hubbard_structure import HubbardStructureData
a, b, c, d = 1.40803, 0.81293, 4.68453, 1.62585
cell = [[a, -b, c], [0.0, d, c], [-a, -b, c]]
sites = [
['Co', 'Co', (0, 0, 0)],
['O', 'O', (0, 0, 3.6608)],
['O', 'O', (0, 0, 10.392)],
['Li', 'Li', (0, 0, 7.0268)],
]
hubbard_structure = HubbardStructureData(cell=cell, sites=sites)
hubbard_structure.initialize_onsites_hubbard("Co", "3d")
hubbard_structure.initialize_intersites_hubbard("Co", "3d", "O", "2p")
hubbard_structure.store()
<HubbardStructureData: uuid: e7f67116-a0d2-476c-a35b-265691ed6b3c (pk: 106)>
Let’s visualize what will be print in the Hubbard card of Quantum ESPRESSO.
from aiida_quantumespresso.utils.hubbard import HubbardUtils
print(HubbardUtils(hubbard_structure).get_hubbard_card())
HUBBARD ortho-atomic
V Co-3d Co-3d 1 1 1e-08
V Co-3d O-2p 1 46 1e-08
As you can see, the desired interactions have been initialized correctly.
This is important because hp.x needs to know which atoms need to be perturbed.
As you will see later, hp.x will take care of adding the remaining interactions with neighbouring atoms.
Important
When you use your own structures, make sure to have your ‘Hubbard atoms’ first in the list of atoms.
This is due to the way the hp.x routine works internally, requiring those to be first.
You can simply do this with the following snippet (IF THE NODE IS YET NOT STORED!):
from aiida_quantumespresso.utils.hubbard import HubbardUtils
HubbardUtils(hubbard_structure).reorder_atoms
Calculating the SCF ground-state#
Now that we have defined the structure, we can calculate its ground-state via an SCF using the PwBaseWorkChain.
We can fill the inputs of the builder of the PwBaseWorkChain through the get_builder_from_protocol() method.
from aiida.engine import run_get_node
from aiida.orm import KpointsData
from aiida_quantumespresso.workflows.pw.base import PwBaseWorkChain
from aiida_quantumespresso.common.types import ElectronicType
kpoints = KpointsData()
kpoints.set_kpoints_mesh([1,1,1])
builder = PwBaseWorkChain.get_builder_from_protocol(
code=data.pw_code, # modify here if you downloaded the notebook
structure=hubbard_structure,
protocol="fast",
electronic_type=ElectronicType.INSULATOR,
overrides={"kpoints":kpoints, "clean_workdir":False}
)
results, pw_node = run_get_node(builder)
results
As you can notice from the results, the workchain (actually, the PwCalculation!) has a remote_folder output.
This is what we need in order to run the HpBaseWorkChain.
DFPT calculation of Hubbard parameters#
We can perturb the ground-state previously found to compute the Hubbard parameters.
Here we will need to use the HpBaseWorkChain, and link the remote_folder previously produced via the parent_scf input.
from aiida.orm import Dict
from aiida_hubbard.workflows.hp.base import HpBaseWorkChain
qpoints = KpointsData()
qpoints.set_kpoints_mesh([1,1,1])
builder = HpBaseWorkChain.get_builder()
builder.hp.code = data.hp_code
builder.hp.hubbard_structure = data.structure
builder.hp.parameters = Dict({"INPUTHP":{"conv_thr_chi": 1e-4}})
builder.hp.qpoints = qpoints
builder.hp.parent_scf = pw_node.outputs.remote_folder
Or via the get_builder_from_protocol:
from aiida_hubbard.workflows.hp.base import HpBaseWorkChain
builder = HpBaseWorkChain.get_builder_from_protocol(
code=data.hp_code, # modify here if you downloaded the notebook
protocol="fast",
parent_scf_folder=pw_node.outputs.remote_folder,
overrides={'hp':{'hubbard_structure':hubbard_structure}},
)
results, hp_node = run_get_node(builder)
results
🚀 Let’s inspect the results!
print(HubbardUtils(results['hubbard_structure']).get_hubbard_card())
HUBBARD ortho-atomic
V Co-3d Co-3d 1 1 9.8969
V Co-3d O-2p 1 11 3.3429
V Co-3d O-2p 1 22 3.3407
Final considerations#
We managed to compute the Hubbard parameters of LiCoO2 fully ab initio! 🎉 However, we had to execute quite a few steps manually, which can be tedious and error prone. Moreover, there are the following considerations:
For larger and more complex structures you will need to perturb many more atoms. Moreover, to get converged results you will need more than one q point. Click here to learn how to parallelize over atoms and q points.
To do a full self-consistent calculation of these parameters, you should relax your structure with the Hubbard parameters from the
hp.xrun, repeat the steps of this tutorial, relax again, and do this procedure over and over till convergence. Learn the automated way here!
Learn more and in details
To learn the full sets of inputs, to use proficiently the get_builder_from_protocol and more, have a look at the following sections:
Note
We suggest to proceed first with the tutorial for point (1) and then the one for point (2). Nevertheless, tutorial (1) is not strictly necessary for (2).