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:
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,
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
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/Psps_for_tests/ # 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 .
# Dynamic quadrupoles calculation for AlAs # Tutorial based on tests/v9/Input/t145.abi made by M. Royo ndtset 5 #Set 1 : ground state self-consistency #************************************* getwfk1 0 kptopt1 1 nqpt1 0 tolvrs1 1.0d-10 #Set 2: Response function calculation of d/dk wave function #********************************************************** iscf2 -3 rfelfd2 2 tolwfr2 1.0d-22 rfdir2 1 1 1 #Set 3: Response function calculation of d2/dkdk wave function #************************************************************* getddk3 2 iscf3 -3 rf2_dkdk3 3 tolwfr3 1.0d-22 #Set 4 : Response function calculation of Q=0 phonons and electric field #*********************************************************************** getddk4 2 rfelfd4 3 rfphon4 1 rfatpol4 1 2 rfdir4 1 1 1 tolvrs4 1.0d-10 prepalw4 2 #Set 5 : Dynamic Quadrupoles calculation #*************************************** optdriver5 10 get1wf5 4 get1den5 4 getddk5 2 getdkdk5 3 lw_qdrpl5 1 ####################################################################### #Common input variables getwfk 1 useylm 1 kptopt 2 #Definition of the unit cell acell 3*10.61 rprim 0.0 0.5 0.5 0.5 0.0 0.5 0.5 0.5 0.0 #Definition of the atom types ntypat 2 znucl 13 33 #Definition of the atoms natom 2 typat 1 2 xred 0.0 0.0 0.0 0.25 0.25 0.25 #Gives the number of bands nband 4 #Definition of the planewave basis set ecut 5.0 # Maximal kinetic energy cut-off, in Hartree #Definition of the k-point grid ngkpt 6 6 6 nshiftk 4 shiftk 0.5 0.5 0.5 0.5 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.5 #Definition of the SCF procedure nstep 100 diemac 9.0 pp_dirpath "$ABI_PSPDIR" pseudos "13al.981214.fhi, As.fhi" ############################################################## # This section is used only for regression testing of ABINIT # ############################################################## ## After modifying the following section, one might need to regenerate the pickle database with runtests.py -r #%%<BEGIN TEST_INFO> #%% [setup] #%% executable = abinit #%% [files] #%% files_to_test = #%% tquad_1.abo, tolnlines= 12, tolabs= 3.e-04, tolrel= 5.00e-04, fld_options= -easy #%% [paral_info] #%% max_nprocs = 4 #%% [extra_info] #%% authors = S. Ponce #%% keywords = DFPT, LONGWAVE #%% description = Dynamic Quadrupoles Calculation for AlAs #%% topics = longwave #%%<END TEST_INFO>
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
# Author: Samuel Ponc\'e # Date: 02/07/2024 # Converter of quadrupole tensor from Abinit to EPW format (Quadrupole.fmt). import numpy as np quad = np.zeros((4,3,3,3)) # atoms, direction, efield, qgrad # Read the input file user_input = input('Enter name of Abinit .abo file\n') AB_file = user_input.strip() with open(AB_file,'r') as N: for lines in N: tmp = lines.split() if (len(tmp) < 1): continue if (tmp[0] == 'natom' and len(tmp) == 2): natom = int(tmp[1]) exit quad = np.zeros((natom,3,3,3)) # atoms, direction, efield, qgrad Found = 0 with open(AB_file,'r') as N: for lines in N: tmp = lines.split() if (len(tmp) < 1): continue if (tmp[0] == 'Quadrupole'): Found = 1 continue if (tmp[0] == 'efidir' and Found == 1): Found = 2 continue if (Found == 2 and tmp[0] == 'Electronic'): Found = 0 continue if (Found == 2): at = int(tmp[1]) - 1 dirc = int(tmp[2]) - 1 efield = int(tmp[0]) - 1 qgrad = int(tmp[3]) - 1 val = np.float(tmp[4]) if np.abs(val) < 1E-6: val = 0.0 quad[at,dirc,efield,qgrad] = val print(' atom dir Qxx Qyy Qzz Qyz Qxz Qxy') for ii in np.arange(natom): at = ii dirc = 0 print(' %3i %3i %12.8f %12.8f %12.8f %12.8f %12.8f %12.8f'%(at+1, dirc+1, quad[at,dirc,0,0], quad[at,dirc,1,1], quad[at,dirc,2,2], quad[at,dirc,1,2], quad[at,dirc,0,2], quad[at,dirc,0,1],)) at = ii dirc = 1 print(' %3i %3i %12.8f %12.8f %12.8f %12.8f %12.8f %12.8f'%(at+1, dirc+1, quad[at,dirc,0,0], quad[at,dirc,1,1], quad[at,dirc,2,2], quad[at,dirc,1,2], quad[at,dirc,0,2], quad[at,dirc,0,1],)) at = ii dirc = 2 print(' %3i %3i %12.8f %12.8f %12.8f %12.8f %12.8f %12.8f'%(at+1, dirc+1, quad[at,dirc,0,0], quad[at,dirc,1,1], quad[at,dirc,2,2], quad[at,dirc,1,2], quad[at,dirc,0,2], quad[at,dirc,0,1],))
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.