3.37 Response to atomic coordinate perturbations using DFPT for Hessian and phonon calculations

In order to get vibration/phonon frequencies, first we need to get dynamical matrix, Let’s define a lattice vector ${𝐑}_{Im}$ as

$${𝐑}_{Im}={𝐑}_I+{𝐑}_m,$$

(3.133)

whereby ${𝐑}_m$ denotes an arbitrary linear combination of ${𝐚}_1$, ${𝐚}_2$, and ${𝐚}_3$. And the dynamical matrix $D_{IJ}(𝐪)$ is a Fourier transform of harmonic Hessian matrix (Force constant) ${\mathrm{\Phi }}_{Im,J}^{harm}$.

	$D_{IJ}(𝐪)$	$={\displaystyle \frac{1}{\sqrt{M_I{M}_J}}}{\displaystyle \sum _m}{\mathrm{\Phi }}_{Im,J}^{harm}\exp \left(i𝐪\cdot {𝐑}_{𝐦}\right)$
		$={\displaystyle \frac{1}{\sqrt{M_I{M}_J}}}{\displaystyle \sum _m}{\displaystyle \frac{d^2{E}_{KS}}{d{𝐑}_{Im}d{𝐑}_J}}\exp \left(i𝐪\cdot {𝐑}_{𝐦}\right)$		(3.134)

Since the finite ($3N\times 3N$) dynamical matrix $𝐃(𝐪)$ would in principle have to be determined for an infinite number of $𝐪$-points in the Brillouin zone. Its diagonalization would produce a set of $3N$ $𝐪$-dependent eigenfrequencies ${\omega }_{\lambda }(𝐪)$ and -vectors ${𝐞}_{\lambda }(𝐪)$.

In order to derive the second order derivatives for total energy analytically (which is the key idea for the DFPT approach), the ground state total energy and force is derived first and we could get this second order derivatives directly. It is should be noted that, here real-space DFPT method is used, so that we could get the real-space force constant directly.

In FHI-aims, the total energy is calculated by using band-energy, a derivation for cluster systems is a following, and the extending to extended systems is straightforward.

$$\frac{1}{2}\int \int \frac{n(𝐫)n({𝐫}^{\prime })}{|𝐫-{𝐫}^{\prime }|}𝑑𝐫𝑑{𝐫}^{\prime }+\frac{1}{2}\sum _I\sum _{J\ne I}\frac{Z_I{Z}_J}{|{𝐑}_I-{𝐑}_J|}+E_{xc}(n)$$

(3.135)

using

$$\Rightarrow \widehat{h_{ks}}=-\frac{1}{2}{\nabla }^2-\sum _I\frac{Z_I}{|𝐫-{𝐑}_I|}+\int \frac{n({𝐫}^{\prime })}{|𝐫-{𝐫}^{\prime }|}𝑑{𝐫}^{\prime }+v_{xc}(𝐫)$$

we have

$$-\frac{1}{2}\int \int \frac{n(𝐫)n({𝐫}^{\prime })}{|𝐫-{𝐫}^{\prime }|}𝑑𝐫𝑑{𝐫}^{\prime }+\frac{1}{2}\sum _I\sum _{J\ne I}\frac{Z_I{Z}_J}{|{𝐮}_I-{𝐮}_I|}$$

(3.136)

$$=\sum _i{f}_i{ϵ}_i-\int [n(𝐫)v_{xc}(𝐫)]𝑑𝐫+E_{xc}(n)$$

$$-\frac{1}{2}\int \int \frac{n𝐫)n({𝐫}^{\prime })}{|𝐫-{𝐫}^{\prime }|}𝑑𝐫𝑑{𝐫}^{\prime }+\frac{1}{2}\sum _I\sum _{J\ne I}\frac{Z_I{Z}_J}{|{𝐑}_I-{𝐑}_I|}$$

(3.137)

$$=\sum _i{f}_i{ϵ}_i-\int [n(𝐫)v_{xc}(𝐫)]𝑑𝐫+E_{xc}(n)$$

$$-\frac{1}{2}\int \int \frac{n(𝐫)n({𝐫}^{\prime })}{|𝐫-{𝐫}^{\prime }|}𝑑𝐫𝑑{𝐫}^{\prime }+\frac{1}{2}\int n(𝐫)\sum _I\frac{Z_I}{|𝐫-{𝐑}_I|}d𝐫$$

$$+\frac{1}{2}\sum _I\sum _{J\ne I}\frac{Z_I{Z}_J}{|{𝐑}_I-{𝐑}_J|}-\frac{1}{2}\int n(𝐫)\sum _I\frac{Z_I}{|𝐫-{𝐑}_I|}d𝐫$$

(3.138)

$$=\sum _i{f}_i{ϵ}_i-\int [n(𝐫)v_{xc}(𝐫)]𝑑𝐫+E_{xc}(n)$$

$$-\frac{1}{2}\int n(𝐫)[\sum _I{V}_I^{free}(|𝐫-{𝐑}_I|)+\delta V_I(|𝐫-{𝐑}_I|)]𝑑𝐫$$

$$-\frac{1}{2}\sum _I{Z}_I[\sum _J{V}_J^{free}(|{𝐑}_J-{𝐑}_I|)+\sum _{J\ne I}\delta V_J(|{𝐑}_J-{𝐑}_I|)]$$

(3.139)

$$=\sum _i{f}_i{ϵ}_i-\int [n(𝐫)v_{xc}(𝐫)]𝑑𝐫+E_{xc}(n)$$

$$-\frac{1}{2}\int n(𝐫)[\sum _I{V}_I^{es,tot}(|𝐫-{𝐑}_I|)]𝑑𝐫$$

$$-\frac{1}{2}\sum _I{Z}_I[V_I^{es,tot}(0)+\sum _{J\ne I}{V}_J^{es,tot}(|{𝐑}_J-{𝐑}_I|)]$$

(3.140)

Here Eq.(3.140) is exactly the one to calculate Kohn-Sham total energy. An extension expression to extended system is in Eq.3.141

	$E_{tot}=-{\displaystyle \frac{1}{N_k}}{\displaystyle \sum _i^{uc}}{f}_{𝐤i}{ϵ}_{𝐤i}-{\displaystyle {\int }_{uc}}[n(𝐫)v_{xc}(𝐫)]𝑑𝐫+E_{xc}(n)$		(3.141)
	$-{\displaystyle \frac{1}{2}}{\displaystyle {\int }_{uc}}n(𝐫)[{\displaystyle \sum _{I,m}}{V}_I^{es,tot}(|𝐫-{𝐑}_{I,m}|)]𝑑𝐫$
	$-{\displaystyle \frac{1}{2}}{\displaystyle \sum _I}{Z}_I[V_I^{es,tot}(0)+{\displaystyle \sum _{(J,m)\ne (I,0)}}{V}_J^{es,tot}(|{𝐑}_{Jm}-{𝐑}_I|)]$

Then we get the analytical force expression for cluster systems:

$$F_I=-\frac{d{E}_{KS}}{d{𝐑}_I}={𝐅}_I^{HF}+{𝐅}_I^P+{𝐅}_I^M$$

$$-[-\frac{\partial (\int n(𝐫){\sum }_I{\displaystyle \frac{Z_I}{|𝐫-{𝐑}_I|}}d𝐫)}{\partial {𝐑}_I}+\frac{\partial ({\displaystyle \frac{1}{2}}{\sum }_I{\sum }_{J\ne I}{\displaystyle \frac{Z_I{Z}_J}{|{𝐑}_I-{𝐑}_J|}})}{\partial {𝐑}_I}]$$

$$-\sum _{\mu \nu }[P_{\mu \nu }\int (\frac{\partial {\chi }_{\mu }}{\partial {𝐑}_I}{\widehat{h}}_{ks}{\chi }_{\nu }+{\chi }_{\mu }{\widehat{h}}_{ks}\frac{\partial {\chi }_{\nu }}{\partial {𝐑}_I})𝑑𝐫-W_{\mu \nu }\frac{\partial S_{\mu \nu }}{\partial {𝐑}_I}]$$

$$-\int [n(𝐫)-n^{MP}(𝐫)]\frac{\partial V^{MP}(𝐫)}{\partial {𝐑}_I}𝑑𝐫$$

(3.142)

Here $P_{\mu \nu }$ refers to density matrix and $W_{\mu \nu }$ refers to energy density matrix. Here the first term is Hellmann-Feynman force using Hellmann-Feynman theorem ; The second term is Pulay term, it comes from basis set dependence of atomic coordinate; The third term is multipole force, it count the contribution from multipole expansion error for Hartree energy calculation (multipole Poisson solver).

Under periodic boundary conditions, we get the following expression for force for extended systems:

$${𝐅}_J^{HF}=Z_J[\frac{\partial V_J^{es,tot}(0)}{\partial {𝐑}_J}+\sum _{(I,m)\ne (J,0)}(\frac{\partial V_I^{es,tot}(|{𝐑}_J-{𝐑}_{Im})|}{\partial {𝐑}_J})]$$

(3.143)

${𝐅}_J^P=-{\displaystyle \frac{2}{N_k}}{\displaystyle \sum _{𝐤,i,\mu ,\nu }}{\displaystyle \int f_{i𝐤}{C}_{\mu i𝐤}^{\ast }\frac{\partial {\phi }_{\mu 𝐤}^{\ast }(𝐫)}{\partial {𝐑}_J}({\widehat{h}}_{ks}-{ϵ}_{i𝐤})C_{\nu i𝐤}{\phi }_{\nu 𝐤}(𝐫)𝑑𝐫},$

(3.144)

In our real-space DFPT method, this harmonic Hessian matrix is calculated directly, as explained in our paper. The Hessian can be split into two part

$$E_{KS}^{(2)}=\frac{d^2{E}_{KS}}{d{𝐑}_{I}d{𝐑}_J}={\mathrm{\Phi }}_{I,J}^{HF}+{\mathrm{\Phi }}_{I,J}^P+{\mathrm{\Phi }}_{I,J}^{MP}$$

(3.145)

$$=[-\frac{{\partial }^2(\int n(𝐫){\sum }_I{\displaystyle \frac{Z_I}{|𝐫-{𝐑}_I|}}d𝐫)}{\partial {𝐑}_I\partial {𝐮}_J}+\frac{{\partial }^2({\displaystyle \frac{1}{2}}{\sum }_I{\sum }_{J\ne I}{\displaystyle \frac{Z_I{Z}_J}{|{𝐑}_I-{𝐑}_J|}})}{\partial {𝐑}_I\partial {𝐑}_J}]$$

$$+\sum _{\mu \nu }\frac{\partial P_{\mu \nu }}{\partial {𝐑}_J}\int (\frac{\partial {\chi }_{\mu }}{\partial {𝐑}_I}{\widehat{h}}_{ks}{\chi }_{\nu }+{\chi }_{\mu }{\widehat{h}}_{ks}\frac{\partial {\chi }_{\nu }}{\partial {𝐑}_I})𝑑𝐫$$

$$+\sum _{\mu \nu }{P}_{\mu \nu }\int [\frac{{\partial }^2{\chi }_{\mu }}{\partial {𝐑}_I\partial {𝐑}_J}{\widehat{h}}_{ks}{\chi }_{\nu }+\frac{\partial {\chi }_{\mu }}{\partial {𝐑}_I}\frac{\partial {\widehat{h}}_{ks}}{\partial {𝐑}_J}{\chi }_{\nu }+\frac{\partial {\chi }_{\mu }}{\partial {𝐑}_I}{\widehat{h}}_{ks}\frac{\partial {\chi }_{\nu }}{\partial {𝐑}_J}]𝑑𝐫$$

$$+\sum _{\mu \nu }{P}_{\mu \nu }\int [\frac{\partial {\chi }_{\mu }}{\partial {𝐑}_J}{\widehat{h}}_{ks}\frac{\partial {\chi }_{\nu }}{\partial {𝐑}_I}+{\chi }_{\mu }\frac{\partial {\widehat{h}}_{ks}}{\partial {𝐑}_J}\frac{\partial {\chi }_{\nu }}{\partial {𝐑}_I}+{\chi }_{\nu }{\widehat{h}}_{ks}\frac{{\partial }^2{\chi }_{\nu }}{\partial {𝐑}_I\partial {𝐑}_J}]𝑑𝐫$$

$$-(\sum _{\mu \nu }\frac{\partial E_{\mu \nu }}{\partial {𝐑}_J}\frac{\partial S_{\mu \nu }}{\partial {𝐑}_I}+\sum _{\mu \nu }{E}_{\mu \nu }\frac{{\partial }^2{S}_{\mu \nu }}{\partial {𝐑}_I\partial {𝐑}_J})$$

$$+\int \frac{\partial [n(𝐫)-n^{MP}(𝐫)]}{\partial {𝐑}_I}\frac{\partial V^{MP}(𝐫)}{\partial {𝐑}_J}𝑑𝐫$$

$$+\int [n(r)-n^{MP}(r)]\frac{{\partial }^2{V}^{MP}(r)}{\partial {𝐑}_I\partial {𝐑}_J}𝑑𝐫$$

It can also be divided into three part: Hellmann-Feynman Hessian, Pulay Hessian, Multipole Hessian. In real calculation, we drop multipole Hessian due to its value is only ${10}^{-3}$ times small compared with the other two.

Under periodic boundary conditions, we get the following expression for force constants:

$${\mathrm{\Phi }}_{Is,J}^{harm}=\frac{d^2{𝐄}^{KS}}{d{𝐑}_{Is}d{𝐑}_J}=-\frac{d{𝐅}_J}{d{𝐑}_{Is}}=-\frac{d{𝐅}_{Is}}{d{𝐑}_J}={\mathrm{\Phi }}_{Is,J}^{HF}+{\mathrm{\Phi }}_{Is,J}^P.$$

(3.146)

	${\mathrm{\Phi }}_{Is,J}^{HF}$	$=$	$-Z_J\left({\displaystyle \frac{d}{d{𝐑}_J}}{\displaystyle \frac{\partial V_J^{es}(0)}{\partial {𝐑}_J}}\right){\delta }_{Is,J0}$
			$-Z_J\left({\displaystyle \frac{d}{d{𝐑}_J}}{\displaystyle \frac{\partial V_I^{es,tot}(|{𝐑}_J-{𝐑}_{Is})|}{\partial {𝐑}_{Is}}}\right)\left(1-{\delta }_{Is,J0}\right),$

in which ${\delta }_{Is,J0}={\delta }_{IJ}{\delta }_{s0}$ denotes a multi-index Kronecker delta.

For the sake of readability, its total derivative is split into four terms:

$${\mathrm{\Phi }}_{Is,J}^P={\mathrm{\Phi }}_{Is,J}^{P-P}+{\mathrm{\Phi }}_{Is,J}^{P-H}+{\mathrm{\Phi }}_{Is,J}^{P-W}+{\mathrm{\Phi }}_{Is,J}^{P-S}.$$

(3.148)

The first term

$${\mathrm{\Phi }}_{Is,J}^{P-P}=2\sum _{\mu m,\nu n}\left(\frac{d{P}_{\mu m,\nu n}}{d{𝐑}_J}\right)\int \frac{\partial {\chi }_{\mu m}(𝐫)}{\partial {𝐑}_{Is}}{\widehat{h}}_{ks}{\chi }_{\nu n}(𝐫)𝑑𝐫$$

(3.149)

accounts for the response of the density matrix $P_{\mu m,\nu n}^{m,n}$. The second term

	${\mathrm{\Phi }}_{Is,J}^{P-H}$	$=$	$2{\displaystyle \sum _{\mu m,\nu n}}{P}_{\mu m,\nu n}\cdot$		(3.153)
			$({\displaystyle \int }{\displaystyle \frac{{\partial }^2{\chi }_{\mu m}(𝐫)}{\partial {𝐑}_{Is}\partial {𝐑}_J}}{\widehat{h}}_{ks}{\chi }_{\nu n}(𝐫)d𝐫$
			$+{\displaystyle \int \frac{\partial {\chi }_{\mu m}(𝐫)}{\partial {𝐑}_{Is}}\frac{d{\widehat{h}}_{ks}}{d{𝐑}_J}{\chi }_{\nu n}(𝐫)𝑑𝐫}$
			$+{\displaystyle \int }{\displaystyle \frac{\partial {\chi }_{\mu m}(𝐫)}{\partial {𝐑}_{Is}}}{\widehat{h}}_{ks}{\displaystyle \frac{\partial {\chi }_{\nu n}(𝐫)}{\partial {𝐑}_J}}d𝐫)$

accounts for the response of the Hamiltonian ${\widehat{h}}_{ks}(𝐤)$, while the third and fourth term

	${\mathrm{\Phi }}_{Is,J}^{P-W}$	$=$	$-2{\displaystyle \sum _{\mu m,\nu n}}{\displaystyle \frac{d{W}_{\mu m,\nu n}}{d{𝐑}_J}}{\displaystyle \int \frac{\partial {\chi }_{\mu m}(𝐫)}{\partial {𝐑}_{Is}}{\chi }_{\nu n}(𝐫)𝑑𝐫}$		(3.154)
	${\mathrm{\Phi }}_{Is,J}^{P-S}$	$=$	$-2{\displaystyle \sum _{\mu m,\nu n}}{W}_{\mu m,\nu n}{\displaystyle \frac{\partial }{\partial {𝐑}_J}}{\displaystyle \int \frac{\partial {\chi }_{\mu m}(𝐫)}{\partial {𝐑}_{Is}}{\chi }_{\nu n}(𝐫)𝑑𝐫}$		(3.155)

for the response of the energy weighted density matrix $W_{\mu m,\nu n}$ and the overlap matrix $S_{\mu m,\nu n}$, respectively, Please note that in all four contributions many terms vanish due to the fact that the localized atomic orbitals ${\chi }_{\mu m}(𝐫)$ are associated with one specific atom/periodic image ${𝐑}_{J(\mu )m}$, which implies, e.g.,

$$\frac{\partial {\chi }_{\mu m}(𝐫)}{\partial {𝐑}_{Is}}=\frac{\partial {\chi }_{\mu m}(𝐫)}{\partial {𝐑}_{Is}}{\delta }_{J(\mu )m,Is}.$$

(3.156)

In the above force constants, it is clear that the first order density matrix ${\displaystyle \frac{\partial P_{\mu \nu }}{\partial {𝐑}_J}}$ and the first order energy density matrix ${\displaystyle \frac{\partial W_{\mu \nu }}{\partial {𝐑}_J}}$ are needed. These first order qualtities are obtained in the DFPT cycle. The flowchart of our DFPT cycle for lattice dynamics is shown in Fig.3.8.

The atomic perturbation response calculation is triggered by setting the atomic_pert_response keyword to DFPT in the control.in file. The atomic_pert_serial keyword can be used to determine whether the CPSCF cycle is carried out for each atomic coordinate in serial, or for all atomic coordinates at the same time (i.e $3N$ CPSCF cycles or 1 longer CPSCF cycle).