Dynamical quadrupoles with DFPT

This tutorial describes the computation of dynamical quadrupoles using density functional perturbation theory (DFPT), using AlAs as an example.

It is assumed the user has already completed the two tutorials RF1 and RF2, and is comfortable with the concepts of ground-state and response properties, including phonons, Born effective charges and dielectric tensor.

The first-order terms in the long-wavelength expansion of a material’s charge response to atomic displacements are the Born effective charges (dynamical dipoles). The second-order terms in this expansion are known as dynamical quadrupoles. Accordingly, the quadrupoles can be considered as the spatial-dispersion of the Born charges.

Dynamical quadrupoles can be used to go beyond the typical dipole approach to characterize the macroscopic, long-range electrostatic fields generated by atomic displacements. This extension has been shown to be important in the Fourier interpolation of several quantities, including: interatomic force constants [Royo2020], perturbed potentials [Brunin2020], and electron-phonon matrix elements [Ponce2021].

In addition, adaptations of these concepts exist for 2D materials in the context of interatomic force-constants [Royo2021] and electron-phonon matrices [Ponce2023].

Completing this tutorial should take approximately 1 hour.

Formalism

The general formalism to calculate spatial-dispersion properties from DFPT has been reported in [Royo2019]. Therein, the dynamical quadrupoles can be calculated from the first \({\bf q}\)-gradient of the second-order total energy response due to an electric field and an atomic displacement:

\[\begin{equation} Q_{\kappa\beta}^{(2,\gamma\delta)}= -2 i ( E_{\gamma}^{\mathcal{E}_\delta^* \tau_{\kappa\beta}} + E_{\delta}^{\mathcal{E}_\gamma^* \tau_{\kappa\beta}}), \end{equation}\]

where \(\kappa\) is a sublattice index and the rest of Greek indexes indicate Cartesian directions. The first \(\bf q\)-gradient of the mixed energy due to an electric field and an atomic displacement is obtained as,

\[\begin{equation} E_{\gamma}^{\mathcal{E}_\delta^* \tau_{\kappa\beta}} = s \int_{\rm BZ} [d^3k] \sum_{m} E_{m{\bf k},\gamma}^{\mathcal{E}_\delta^* \tau_{\kappa\beta}} + \frac{1}{2} \int_{\Omega} \int K_{\gamma}({\bf r},{\bf r}') n^{\mathcal{E}_\delta} ({\bf r}) n^{\tau_{\kappa\beta}}({\bf r}') d^3r d^3r', \end{equation}\]

where \(s=2\) is the spin multiplicity, \(K_{\gamma}({\bf r},{\bf r}')\) is the \(\bf q\)-derivative of the Hartree and exchange-correlation kernel in the direction \(\gamma\), and \(n^\lambda\) is the first-order perturbed density due to a generic pertubation \(\lambda\). The first term in the above equation is a band-resolved contribution given by

\[\begin{equation} E_{m{\bf k},\gamma}^{\mathcal{E}_\delta^* \tau_{\kappa\beta}} = \langle u_{m{\bf k}}^{\mathcal{E}_{\delta}} | \partial_{\gamma} \hat{H}^{(0)}_{{\bf k}} | u_{m{\bf k}}^{\tau_{\kappa\beta}} \rangle + \langle u_{m{\bf k}}^{\mathcal{E}_{\delta}}| \partial_{\gamma} \hat{Q}_{{\bf k}} \hat{\mathcal{H}}_{{\bf k}}^{\tau_{\kappa\beta}} | u_{m{\bf k}}^{(0)} \rangle + \langle u_{m{\bf k}}^{(0)} | V^{\mathcal{E}_{\delta}} \partial_{\gamma} \hat{Q}_{{\bf k}}| u_{m{\bf k}}^{\tau_{\kappa\beta}} \rangle + \langle u_{m{\bf k}}^{\mathcal{E}_{\delta}}| \hat{H}_{{\bf k},\gamma}^{\tau_{\kappa\beta}} | u_{m{\bf k}}^{(0)} \rangle +i \langle u_{m{\bf k},\gamma}^{A_\delta}| u_{m{\bf k}}^{\tau_{\kappa\beta}} \rangle . \end{equation}\]

where \(\partial_{\gamma} \hat{H}^{(0)}_{{\bf k}}\) is the velocity operator; \(\hat{\mathcal{H}}_{{\bf k}}^{\tau_{\kappa\beta}}\) and \(\hat{H}_{{\bf k},\gamma}^{\tau_{\kappa\beta}}\) are, respectively, the atomic-displacement first-order Hamiltonian (the caligraphic symbol indicates that self-consistent field potentials are included) and its gradient along \(\gamma\); \(V^{\mathcal{E}_{\delta}}\) is the first-order potential due to an electric field; and \(\partial_{\gamma} \hat{Q}_{{\bf k}}\) is the gradient of the conduction-band projector. In addition, the above equation comprises a set of ground-state (\(u_{m{\bf k}}^{(0)}\)) and first-order wave functions. Among the latter one finds the standard response functions due to an electric field (\(u_{m{\bf k}}^{\mathcal{E}_{\delta}}\)) and an atomic displacement (\(u_{m{\bf k}}^{\tau_{\kappa\beta}}\)), along with the response function due to a gradient of a vector potential (\(u_{m{\bf k},\gamma}^{A_\delta}\)). A discussion on how the latter can be calculated from the auxiliary d2_dkdk second-order response functions (see rf2_dkdk) is provided in [Zabalo2022].

Note

Supposing you made your own installation of ABINIT, the input files to run the examples are in the ~abinit/tests/ directory where ~abinit is the absolute path of the abinit top-level directory. If you have NOT made your own install, ask your system administrator where to find the package, especially the executable and test files.

In case you work on your own PC or workstation, to make things easier, we suggest you define some handy environment variables by executing the following lines in the terminal:

export ABI_HOME=Replace_with_absolute_path_to_abinit_top_level_dir # Change this line
export PATH=$ABI_HOME/src/98_main/:$PATH      # Do not change this line: path to executable
export ABI_TESTS=$ABI_HOME/tests/             # Do not change this line: path to tests dir
export ABI_PSPDIR=$ABI_TESTS/Pspdir/  # Do not change this line: path to pseudos dir

Examples in this tutorial use these shell variables: copy and paste the code snippets into the terminal (remember to set ABI_HOME first!) or, alternatively, source the set_abienv.sh script located in the ~abinit directory:

source ~abinit/set_abienv.sh

The ‘export PATH’ line adds the directory containing the executables to your PATH so that you can invoke the code by simply typing abinit in the terminal instead of providing the absolute path.

To execute the tutorials, create a working directory (Work*) and copy there the input files of the lesson.

Most of the tutorials do not rely on parallelism (except specific tutorials on parallelism). However you can run most of the tutorial examples in parallel with MPI, see the topic on parallelism.

Quadrupole calculation in fcc AlAs

Before beginning, you might consider creating a different subdirectory to work in. Why not create Work_quad ?

The file tquad_1.abi is the input file for the first step (GS + DFPT perturbations for all the \(\qq\)-points in the IBZ). Copy it to the working directory with:

cd $ABI_TESTS/tutorespfn/Input
mkdir Work_quad
cd Work_quad
cp ../tquad_1.abi .

View tests/tutorespfn/Input/tquad_1.abi

This step might be quite time-consuming so you may want to immediately start the job in background with:

abinit tquad_1.abi > tquad_1.log 2> err &

The logic behind this multi-dataset calculation is to sequentialy calculate the required ingredients to be plugged in the quadrupole equations shown above. This involves:

• DS1: Perform ground-state calculation.
• DS2 and DS3: Perform ddk and d2_dkdk (rf2_dkdk=3 is mandatory) response function calculations.
• DS4: Perform response function calculations at q =Γ to atomic displacements and electric field, including the option prepalw=2.
• DS5: Perform a longwave DFTP calculation of third-order energy derivatives (optdriver=10 and lw_qdrpl=1).

Upon finalization, open the output and look for the Quadrupole data block:

 Quadrupole tensor, in cartesian coordinates,
 efidir   atom   atdir    qgrdir          real part        imaginary part
    1       1       1       1           -0.0000000566        -0.0000000000
    1       1       2       1           -0.0000000157        -0.0000000000
    1       1       3       1            0.0000000157        -0.0000000000
    1       2       1       1           -0.0000001675        -0.0000000000
    1       2       2       1           -0.0000000042        -0.0000000000
    1       2       3       1            0.0000000042        -0.0000000000
    2       1       1       1           -0.0000000362        -0.0000000000
    2       1       2       1           -0.0000000205        -0.0000000000
    2       1       3       1           13.4866068621        -0.0000000000
    2       2       1       1           -0.0000000859        -0.0000000000
    2       2       2       1           -0.0000000816        -0.0000000000
    2       2       3       1           -6.0008872938        -0.0000000000
    3       1       1       1           -0.0000000362        -0.0000000000
    3       1       2       1           13.4866068621        -0.0000000000
    3       1       3       1           -0.0000000205        -0.0000000000
    3       2       1       1           -0.0000000859        -0.0000000000
    3       2       2       1           -6.0008872938        -0.0000000000
    3       2       3       1           -0.0000000816        -0.0000000000

Since we are in a binary FCC solid, there are only two independent quadrupoles values given by \(Q_{\kappa\beta}^{\gamma\delta} = Q_\kappa |\varepsilon_{\beta\gamma\delta}|\) where \(\varepsilon_{\beta\gamma\delta}\) is the Levi-Cevita symbol.

In case you want to use the quadrupole tensor for the EPW code, you can use the python convertor located in abinit/scripts/post_processing/ab2epw_quad.py

View ../scripts/post_processing/ab2epw_quad.py

If you run the python code, the script will ask for the name of the output. Once provided you should obtain the following result:

  atom   dir       Qxx             Qyy         Qzz          Qyz            Qxz         Qxy
    1     1      0.00000000   0.00000000   0.00000000  13.48660685   0.00000000   0.00000000
    1     2      0.00000000   0.00000000   0.00000000   0.00000000  13.48660686   0.00000000
    1     3      0.00000000   0.00000000   0.00000000   0.00000000   0.00000000  13.48660686
    2     1      0.00000000   0.00000000   0.00000000  -6.00088730   0.00000000   0.00000000
    2     2      0.00000000   0.00000000   0.00000000   0.00000000  -6.00088729   0.00000000
    2     3      0.00000000   0.00000000   0.00000000   0.00000000   0.00000000  -6.00088729

which can be copied in a file name Quadrupole.fmt and used in EPW.