Dielectric properties of Si#

In this tutorial you will learn how to calculate the high-frequency dielectric tensor of silicon using aiida-vibroscopy.

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 other tutorials.

Let’s get started!

Finite electric fields#

In the harmonic and Born-Oppenheimer approximation, the high-frequency dielectric tensor is defined as second order derivative of the polarization of the system in respect to an external electric field \(\mathcal{E}\). The polarization \(\mathbf{P}\) is computed via the modern theory of polarization, using a Berry phase formalism. An external electric field can still be applied, even though the Hamiltonian will result not bounded from below, and translational symmetry is broken (hence, in principle no Bloch’s theorem).

Note

All these difficulties were solved in the early 2000s. See:

Ivo Souza, Jorge Iniguez, and David Vanderbilt, Phys. Rev. Lett, 89, 117602 (2002)
Paolo Umari and Alfredo Pasquarello, Phys. Rev. Lett, 89, 157602 (2002)

Extending the Kohn-Sham Hamiltonian, we can include the external electric field as follows:

(1)#\[\begin{equation} \mathcal{F}[\{ \psi \}] = E_{KS}[\rho] - \Omega \mathbf{P}[\{ \psi \}] \cdot \mathcal{E} \end{equation}\]

This extended functional is usually referred to as electric enthalpy functional.

One subtlelty of this formalism is the slow convergence is respect to k points. In particular, to converge well the polarization operator, one has to sample with many k points the direction of the induced polarization (i.e. the direction of the applied electric field).

In the following workflows, we will define the amount of k point in this direction through what we call parallel distance between k points. This distance is used to define the number of k points. The smaller it is, the more points will be defined, hence more accurate results.

The high-frequency dielectric tensor is computed like frozen phonons: forces become the polarization, atomic displacements become the electric field step. In the DielectricWorkChain, this step can either be automatically chosen or set manually.

Important

It is fundamental that the electric field step is properly chosen. In fact, only values below the critical electric field \(\mathcal{E}_c\) can be chosen. An estimate of such value can be done knowing the amount of k points and the band gap of the material.

(2)#\[\begin{equation} \mathcal{E}_c \approx \frac{E_{gap} }{e \cdot a \cdot N_k } \end{equation}\]

Thus, increasing the number of k points signifies a lower critical field, hence we are allowed to use smaller field steps. As a consequence, trade off with convergence must be made. Nevertheless, for the majority of cases this does not represent a limitation.

The `DielectricWorkChain`#

As for phonons, we need to compute different ground-states with different displacements, in this case different directions of the electric field. We can use the DielectricWorkChain to do this in automated fashion.

Note

As again for phonons, also the dielectric properties can be computed in the linear regime as perturbative respones of the system. Quantum ESPRESSO can carry out such simulation using the PHonon code (ph.x routine). Nevertheless, not all functionals are implemented. For example, DFT+U(+V) using ortho-atomic projectors is not implemented.

Let’s import the WorkChain and run it! We use the get_builder_from_protocol to get a prefilled builder with all inputs.

Note

These inputs should be considered not as converged parameters, but as a good starting point. You may also need to tweak some parameters, e.g. add magnetization etc., depending on your case.

from aiida.plugins import WorkflowFactory
from aiida.engine import run_get_node

DielectricWorkChain = WorkflowFactory("vibroscopy.dielectric")

builder = DielectricWorkChain.get_builder_from_protocol(
    code=data.pw_code,
    structure=structure,
    protocol="fast",
    overrides={"property":"dielectric"}
)

results, calc = run_get_node(builder)

Show code cell output Hide code cell output

05/24/2023 11:02:04 AM <4148633> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [98|DielectricWorkChain|run_base_scf]: launching base scf PwBaseWorkChain<106>
05/24/2023 11:02:04 AM <4148633> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [106|PwBaseWorkChain|run_process]: launching PwCalculation<109> iteration #1
05/24/2023 11:02:11 AM <4148633> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [106|PwBaseWorkChain|sanity_check_insufficient_bands]: PwCalculation<109> run with smearing and highest band is occupied
05/24/2023 11:02:11 AM <4148633> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [106|PwBaseWorkChain|sanity_check_insufficient_bands]: BandsData<112> has invalid occupations: Occupation of 0.00930717805266601 at last band lkn<0,0,20>
05/24/2023 11:02:11 AM <4148633> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [106|PwBaseWorkChain|sanity_check_insufficient_bands]: PwCalculation<109> had insufficient bands
05/24/2023 11:02:11 AM <4148633> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [106|PwBaseWorkChain|sanity_check_insufficient_bands]: Action taken: increased number of bands to 24 and restarting from the previous charge density.
05/24/2023 11:02:11 AM <4148633> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [106|PwBaseWorkChain|inspect_process]: PwCalculation<109> finished successfully but a handler was triggered, restarting
05/24/2023 11:02:11 AM <4148633> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [106|PwBaseWorkChain|run_process]: launching PwCalculation<117> iteration #2
05/24/2023 11:02:18 AM <4148633> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [106|PwBaseWorkChain|results]: work chain completed after 2 iterations
05/24/2023 11:02:18 AM <4148633> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [106|PwBaseWorkChain|on_terminated]: remote folders will not be cleaned
05/24/2023 11:02:19 AM <4148633> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [98|DielectricWorkChain|run_nscf]: launching base scf PwBaseWorkChain<125>
05/24/2023 11:02:19 AM <4148633> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [125|PwBaseWorkChain|run_process]: launching PwCalculation<128> iteration #1
05/24/2023 11:02:34 AM <4148633> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [125|PwBaseWorkChain|results]: work chain completed after 1 iterations
05/24/2023 11:02:34 AM <4148633> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [125|PwBaseWorkChain|on_terminated]: remote folders will not be cleaned
05/24/2023 11:02:35 AM <4148633> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [98|DielectricWorkChain|run_null_field_scfs]: launching PwBaseWorkChain<142> with null electric field 0
05/24/2023 11:02:36 AM <4148633> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [142|PwBaseWorkChain|run_process]: launching PwCalculation<145> iteration #1
05/24/2023 11:03:03 AM <4148633> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [142|PwBaseWorkChain|results]: work chain completed after 1 iterations
05/24/2023 11:03:03 AM <4148633> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [142|PwBaseWorkChain|on_terminated]: remote folders will not be cleaned
05/24/2023 11:03:04 AM <4148633> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [98|DielectricWorkChain|run_electric_field_scfs]: launching PwBaseWorkChain<153> with electric field index 2 and sign 1.0 iteration #0
05/24/2023 11:03:06 AM <4148633> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [153|PwBaseWorkChain|run_process]: launching PwCalculation<156> iteration #1
05/24/2023 11:07:27 AM <4148633> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [153|PwBaseWorkChain|results]: work chain completed after 1 iterations
05/24/2023 11:07:27 AM <4148633> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [153|PwBaseWorkChain|on_terminated]: remote folders will not be cleaned
05/24/2023 11:07:29 AM <4148633> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [98|DielectricWorkChain|run_electric_field_scfs]: launching PwBaseWorkChain<164> with electric field index 2 and sign 1.0 iteration #1
05/24/2023 11:07:31 AM <4148633> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [164|PwBaseWorkChain|run_process]: launching PwCalculation<167> iteration #1
05/24/2023 11:14:32 AM <4148633> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [164|PwBaseWorkChain|results]: work chain completed after 1 iterations
05/24/2023 11:14:32 AM <4148633> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [164|PwBaseWorkChain|on_terminated]: remote folders will not be cleaned
05/24/2023 11:14:35 AM <4148633> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [98|DielectricWorkChain|run_numerical_derivatives]: launching NumericalDerivativesWorkChain<178> for computing numerical derivatives.
05/24/2023 11:14:35 AM <4148633> aiida.engine.processes.functions: [WARNING] function `generate_preprocess_data` has invalid type hints: unsupported operand type(s) for |: 'AbstractNodeMeta' and 'NoneType'
05/24/2023 11:14:35 AM <4148633> aiida.engine.processes.functions: [WARNING] function `generate_phonopy_data` has invalid type hints: unsupported operand type(s) for |: 'AbstractNodeMeta' and 'NoneType'
05/24/2023 11:14:42 AM <4148633> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [98|DielectricWorkChain|on_terminated]: cleaned remote folders of calculations: 109 117 128 145 156 167

These are the results:

results

{'critical_electric_field': <Float: uuid: 38a272fa-7624-470b-95d9-9a8fd72d3c3f (pk: 135) value: 0.00091230365735908>,
 'electric_field_step': <Float: uuid: 4cdbe871-803c-4e00-9aa0-13c2c337ca3f (pk: 139) value: 0.00045615182867954>,
 'fields_data': {'field_index_2': {'0': <TrajectoryData: uuid: fc58409c-8948-4ef6-b025-705ec01d80ac (pk: 160)>,
   '1': <TrajectoryData: uuid: 24d8fb1c-0b1e-4e45-9b18-2f137e06ea51 (pk: 171)>}},
 'tensors': AttributeDict({'numerical_accuracy_2_step_2': <ArrayData: uuid: 2edd1da2-69bd-44a9-b97a-d96cb6bf39de (pk: 182)>, 'numerical_accuracy_2_step_1': <ArrayData: uuid: 5866824b-274c-46d3-8bdf-8dbab1f2ef1d (pk: 183)>, 'numerical_accuracy_4': <ArrayData: uuid: a476ac2a-3992-47a6-80f6-f1f79c9340c3 (pk: 184)>})}

As you can see, we have a tensors output. This is a namespace, meaning it contains more outputs. If we inspect it, we can see we have some other names, all starting with numerical_accuracy_*. This is because the workflow decided to use more points (electric field steps) to have a better accuracy for the numerical differentiation. In this case, up to the 4th order accuracy is given, using the following formula:

(3)#\[\begin{equation} \left . \frac{\partial f(x)}{\partial x} \right|_{x=0} = \frac{1}{12 h} \left [ -f(2h) +8f(h) -8f(-h) +f(-2h) \right ] + \mathcal{O}(h^4) \end{equation}\]

\(f(h)\) will then be \(\mathbf{P}(\delta \mathbf(E))\), a three variable function computing three dimensional output; hence, we will get a matrix out of the differentiation. This expression goes as \(\mathcal{O}(h^4)\), meaning the numerical error due to the finite step will be much smaller then with a second order formula.

Important

Having computed the 4th order accuracy, one can compute also the 2nd order accuracy with the two external and two internal points!

Have a look at the finite difference coefficients to convince yourself.

This is extremely useful, as one can check a posteriori the convergence in respect to finite size step.

Let’s print the results!

print("Numerical accuracy order 4:")
print(calc.outputs.tensors.numerical_accuracy_4.get_array("dielectric").round(5), "\n")
print("Numerical accuracy order 2 (inner points):")
print(calc.outputs.tensors.numerical_accuracy_2_step_1.get_array("dielectric").round(5), "\n")
print("Numerical accuracy order 2 (external points):")
print(calc.outputs.tensors.numerical_accuracy_2_step_2.get_array("dielectric").round(5))

Numerical accuracy order 4:
[[11.08489  0.       0.     ]
 [-0.      11.08489  0.     ]
 [ 0.      -0.      11.08489]] 

Numerical accuracy order 2 (inner points):
[[11.08526  0.       0.     ]
 [-0.      11.08526  0.     ]
 [ 0.      -0.      11.08526]] 

Numerical accuracy order 2 (external points):
[[11.0864 -0.      0.    ]
 [ 0.     11.0864  0.    ]
 [-0.     -0.     11.0864]]

We can conclude the 2nd order accuracy is accurate till the third digit, much more accurate then the convergence error!

print("The estimated critical electric field step: ", calc.outputs.critical_electric_field.value)
print("The electric field step: ", calc.outputs.electric_field_step.value)

The estimated critical electric field step:  0.00091230365735908
The electric field step:  0.00045615182867954

Analysing the workflow#

Many things happened for computing this dielectric tensor. Let’s inspect the process tree.

%verdi process status {calc.pk}

DielectricWorkChain<98> Finished [0] [11:results]
    ├── create_kpoints_from_distance<99> Finished [0]
    ├── create_directional_kpoints<102> Finished [0]
    ├── PwBaseWorkChain<106> Finished [0] [3:results]
    │   ├── PwCalculation<109> Finished [0]
    │   └── PwCalculation<117> Finished [0]
    ├── PwBaseWorkChain<125> Finished [0] [3:results]
    │   └── PwCalculation<128> Finished [0]
    ├── compute_critical_electric_field<134> Finished [0]
    ├── get_accuracy_from_critical_field<136> Finished [0]
    ├── get_electric_field_step<138> Finished [0]
    ├── PwBaseWorkChain<142> Finished [0] [3:results]
    │   └── PwCalculation<145> Finished [0]
    ├── PwBaseWorkChain<153> Finished [0] [3:results]
    │   └── PwCalculation<156> Finished [0]
    ├── PwBaseWorkChain<164> Finished [0] [3:results]
    │   └── PwCalculation<167> Finished [0]
    ├── subtract_residual_forces<175> Finished [0]
    └── NumericalDerivativesWorkChain<178> Finished [0] [None]
        ├── generate_preprocess_data<179> Finished [0]
        └── compute_nac_parameters<181> Finished [0]

Here the explanation:

The k points grids for SCF and electric field are generated.
A base ground-state calculation is performed via a PwBaseWorkChain. This will be used in the following as starting point for all the electric field steps.
An NSCF is performed to evaluate the band gap of the material.
The critical electric field is evaluated.
From this estimate, the numerical accuracy (i.e. number of points) and the electric field step are extrapolated. The accuracy is 4, the value of the field is \(\approx 5\cdot 10^{-4}\) Ry atomic units (1 Ry a.u. \(\approx\) \(36.3609 \cdot 10^{10}\) V/m).
A ground-state elecrtric field with \(\mathcal{E}=0\) is performed for each inequivalent k point grid.
Once finished, SCF with electric field \(\mathcal{E}=\delta \mathbf{E}\) are run, restarting from the previous calculations.
Once finished, SCF with electric field \(\mathcal{E}=2 \delta \mathbf{E}\) are run, restarting from the previous calculations.
The numerical derivatives are performed, and tensors are stored.

Can you explain why only one direction?

We used symmetries to reduce the number of directions to compute. In materials, we can exploit the fact that different operations are connected by symmetry. As an example, if we have inversion symmetry, it holds (for a one dimensional function):

(4)#\[\begin{equation} f(x) = -f(-x) \end{equation}\]

Thus, we would save one calculation. In silicon, we have only one symmetrically inequivalent direction, thanks to its very symmetric point group. The power of symmetries!

You can also inspect this data to convince yourself.

calc.outputs.fields_data

AttributeDict({'field_index_2': AttributeDict({'0': <TrajectoryData: uuid: fc58409c-8948-4ef6-b025-705ec01d80ac (pk: 160)>, '1': <TrajectoryData: uuid: 24d8fb1c-0b1e-4e45-9b18-2f137e06ea51 (pk: 171)>})})

Born effective charges#

With this method, we can also compute the Born effective charges \(Z^*\). They are the derivatives of the forces in respect to electric fields, in contrast with the polarization for the dielectric tensor. In the case of silicon, as this is a homopolar material, they will be zero, as you can verify. In the tutorial on AlAs, you will calculate non zero BEC!

calc.outputs.tensors.numerical_accuracy_4.get_array("born_charges").round(5)

array([[[ 0., -0., -0.],
        [-0., -0., -0.],
        [ 0., -0., -0.]],

       [[ 0., -0.,  0.],
        [-0.,  0., -0.],
        [ 0., -0.,  0.]],

       [[-0., -0., -0.],
        [ 0., -0.,  0.],
        [-0., -0., -0.]],

       [[-0., -0., -0.],
        [ 0.,  0., -0.],
        [-0.,  0.,  0.]],

       [[-0.,  0.,  0.],
        [ 0., -0.,  0.],
        [ 0.,  0.,  0.]],

       [[-0.,  0.,  0.],
        [ 0.,  0., -0.],
        [ 0., -0., -0.]],

       [[ 0.,  0., -0.],
        [-0., -0., -0.],
        [-0., -0., -0.]],

       [[ 0., -0.,  0.],
        [ 0.,  0.,  0.],
        [-0.,  0., -0.]]])

Convergence with k points#

As a final remark, we show here how the convergence with k points is speeded up by sampling in this directional manner, compared to the most ‘traditional’ uniform mesh.

<matplotlib.legend.Legend at 0x7feaa4445d30>

_images/5251f80ba5b7ea199011b96afbf173429cf019ea2ac1930bbb855f6d4669d9ab.png

Final considerations#

We managed to compute the dielectric tensor of silicon fully ab initio! :tada:

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:

Dielectric properties of Si

Contents

Dielectric properties of Si#

Finite electric fields#

The DielectricWorkChain#

Analysing the workflow#

Born effective charges#

Convergence with k points#

Final considerations#

The `DielectricWorkChain`#