Procedures

$$\alpha(\mathrm iu)=\frac{\alpha_0}{1+u^2/\omega^2},\qquad \partial\alpha(\mathrm iu)=\alpha(\mathrm iu)\left( \frac{\partial\alpha_0}{\alpha_0}+ \frac2\omega\frac{\partial\omega}{1+\omega^2/u^2} \right)$$

$$\begin{aligned} B(R,\gamma) &=\operatorname{erfc}(\gamma R) +\frac{2\gamma R}{\sqrt\pi}\mathrm e^{-(\gamma R)^2} \\ \partial B(R,\gamma) &=-\frac4{\sqrt\pi}(\gamma R)^2\mathrm e^{-(\gamma R)^2} (R\partial\gamma+\gamma\partial R) \end{aligned}$$

$$\bar C_6=\frac3\pi\int_0^\infty\mathrm du\,\bar\alpha(u)^2,\qquad \partial\bar C_6=\frac6\pi\int_0^\infty\mathrm du \bar\alpha(u)\partial\bar\alpha(u)$$

$$\begin{aligned} C(r,\gamma) &=3\operatorname{erfc}(\gamma R) +\frac{2\gamma R}{\sqrt\pi}(3+2(\gamma R)^2)\mathrm e^{-(\gamma R)^2} \\ \partial C(R,\gamma) &=-\frac8{\sqrt\pi}(\gamma R)^4\mathrm e^{-(\gamma R)^2} (R\partial\gamma+\gamma\partial R) \end{aligned}$$

$$\begin{gathered} f_{(ij)}=\frac1{1+\exp\big({-}a(\eta-1)\big)},\qquad \eta=\frac{R_{(ij)}}{S_{\text{vdW}(ij)}}\equiv \frac{R_{(ij)}}{\beta R_{\text{vdW}(ij)}} \\ \frac{\mathrm df}{\mathrm dR_c}= \frac a{2+2\cosh\big(a(\eta-1)\big)}\frac{\mathrm d\eta}{\mathrm dR_c},\qquad \frac{\mathrm d\eta}{\mathrm dR_c}= \frac{R_c}{RS_\text{vdW}}- \frac{R}{S_\text{vdW}^2}\frac{\mathrm dS_\text{vdW}}{\mathrm dR_c} \end{gathered}$$

Form either a real or a complex dipole matrix.

The real-typed version is equivalent to $\mathbf q=0$.

$$\boldsymbol{\mathcal A}:=[\mathbf a_1\mathbf a_2\mathbf a_3],\qquad\boldsymbol{\mathcal B}:=[\mathbf b_1\mathbf b_2\mathbf b_3] \\ \boldsymbol{\mathcal B}=2\pi(\boldsymbol{\mathcal A}^{-1})^\mathrm T,\qquad \partial\boldsymbol{\mathcal B}=-\big((\partial\boldsymbol{\mathcal A})\boldsymbol{\mathcal A}^{-1}\big)^\mathrm T\boldsymbol{\mathcal B} \\ \mathbf R_\mathbf n=\boldsymbol{\mathcal A}\mathbf n,\qquad\partial\mathbf R_\mathbf n=(\partial\boldsymbol{\mathcal A})\mathbf n, \\ \mathbf G_\mathbf m=\boldsymbol{\mathcal B}\mathbf m,\qquad \partial\mathbf G_\mathbf m=-\big((\partial\boldsymbol{\mathcal A})\boldsymbol{\mathcal A}^{-1}\big)^\mathrm T\mathbf G_\mathbf m, \\ \frac{\partial G_{\mathbf ma}}{\partial A_{bc}}=-\mathcal A^{-1}_{ca}G_{\mathbf mb}$$

$$\begin{gathered} \mathbf T_{ij}(\mathbf q)=\sum_{\mathbf n}\mathbf T(\mathbf R_{\mathbf nij})\mathrm e^{-\mathrm i\mathbf q\cdot\mathbf R_{\mathbf nij}},\quad\mathbf R_{\mathbf nij}=\mathbf R_j+\mathbf R_\mathbf n-\mathbf R_i \\ \frac{\mathrm d\mathbf R_{\mathbf nij}}{\mathrm d\mathbf R_k}=(\delta_{jk}-\delta_{ik})\mathbf I \\ \mathbf{T}_{ij}(\mathbf{q})\approx\mathbf{T}^\text{Ew}_{ij}(\mathbf{q}) =\sum_\mathbf n^{|\mathbf R_{\mathbf nij}|<R_\text c}\mathbf T^\text{erfc}(\mathbf R_{\mathbf nij};\gamma)\mathrm e^{-\mathrm i\mathbf q\cdot\mathbf R_{\mathbf nij}} +\frac{4\pi}{V_\text{uc}}\sum_{\mathbf m}^{0<|\mathbf k_\mathbf m|<k_\text c}\mathbf{\hat k}_\mathbf m\otimes\mathbf{\hat k}_\mathbf m\,\mathrm e^{-\frac{k_\mathbf m^2}{4\gamma^2}-\mathrm i\mathbf G_\mathbf m\cdot\mathbf R_{ij}} \\ -\frac{4\gamma^3}{3\sqrt\pi}\delta_{ij}\mathbf I +\delta(\mathbf q)\frac{4 \pi}{3 V_\text{uc}}\mathbf I,\qquad \mathbf k_\mathbf m=\mathbf G_\mathbf m+\mathbf q \end{gathered}$$

$$\partial\left(\frac{4\pi}{V_\text{uc}}\right)=-\frac{4\pi}{V_\text{uc}}\frac{\partial V_\text{uc}}{V_\text{uc}},\qquad\frac{\partial V_\text{uc}}{\partial\boldsymbol{\mathcal A}}=V_\text{uc}(\boldsymbol{\mathcal A}^{-1})^\mathrm T \\ \partial(k^2)=2\mathbf k\cdot\partial\mathbf k \\ \mathbf{\hat k}\otimes\partial\mathbf{\hat k}=\frac{\mathbf k\otimes\partial\mathbf k}{k^2}-\frac{\mathbf k\otimes\mathbf k}{2k^4}\partial(k^2)$$

$$\gamma:=\frac{2.5}{\sqrt[3]{V_\text{uc}}},\quad R_\text c:=\frac6\gamma,\quad k_\text c:=10\gamma$$

Get MBD energy.

For a nonperiodic system, the method just transforms $C_6$ coefficients to frequencies, and performs a single call to get_mbd_hamiltonian_energy. For a periodic system, the method integrates the energy over the frist Brillouin zone.

$$E=\int_\text{FBZ}\mathrm d\mathbf q\,E(\mathbf q)\approx\frac1{N_k}\sum_i^{N_k}E(\mathbf q_i) \\ \mathbf q_i=\boldsymbol{\mathcal B}\mathbf n_i,\qquad\partial\mathbf q_i=-\big((\partial\boldsymbol{\mathcal A})\boldsymbol{\mathcal A}^{-1}\big)^\mathrm T\mathbf q_i$$

Form and solve either a real or a complex MBD Hamiltonian.

The real-typed version is equivalent to $\mathbf q=0$.

$$\begin{gathered} E_\text{MBD}(\mathbf q)=\frac12\operatorname{Tr}\big(\sqrt{\mathbf Q(\mathbf q)}\big)- 3\sum_i\frac{\omega_i}2,\qquad \mathbf Q_{ij}(\mathbf q)=\omega_i^2\delta_{ij}\mathbf I+ \omega_i\omega_j\sqrt{\alpha_{0,i}\alpha_{0,j}}\mathbf T_{ij}(\mathbf q) \\ \mathbf Q(\mathbf q)\equiv \mathbf C(\mathbf q)\boldsymbol\Lambda(\mathbf q)\mathbf C(\mathbf q)^\dagger,\qquad \boldsymbol\Lambda(\mathbf q) \equiv\operatorname{diag}(\{\tilde\omega_i(\mathbf q)^2\}),\qquad \operatorname{Tr}\big(\sqrt{\mathbf Q(\mathbf q)}\big) =\sum_i\tilde\omega_i(\mathbf q) \end{gathered}$$

$$\begin{aligned} \partial E_\text{MBD}&=\frac14\operatorname{Tr}\big( \mathbf C\boldsymbol\Lambda^{-\frac12}\mathbf C^\dagger \partial\mathbf Q \big)- 3\sum_i\frac{\partial\omega_i}2 \\ \frac{\partial E_\text{MBD}}{\partial X_i}&= \operatorname{Re}\frac12\sum_{p\zeta}( \mathbf C\boldsymbol\Lambda^{-\frac12}\mathbf C^\dagger )_{i\zeta,p} \frac{\partial Q_{p,i\zeta}}{\partial X_i}- \frac32\frac{\partial\omega_i}{\partial X_i} \end{aligned}$$

$$\begin{aligned} \partial\mathbf Q_{ij}=& 2\delta_{ij}\omega_i\partial\omega_i\mathbf I+ \omega_i\omega_j\sqrt{\alpha_{0,i}\alpha_{0,j}}\mathbf T_{ij}\left( \frac{\partial\omega_i}{\omega_i}+ \frac{\partial\omega_j}{\omega_j}+ \frac12\frac{\partial\alpha_{0,i}}{\alpha_{0,i}}+ \frac12\frac{\partial\alpha_{0,j}}{\alpha_{0,j}} \right) \\ &+\omega_i\omega_j\sqrt{\alpha_{0,i}\alpha_{0,j}} \partial\mathbf T_{ij} \end{aligned}$$

Get screened MBD energy.

Get TS energy.

$$\omega=\frac{4C_6}{3\alpha_{0}^2},\qquad \partial\omega=\omega\left( \frac{\partial C_6}{C_6}-\frac{2\partial\alpha_0}{\alpha_0} \right)$$

$$\begin{gathered} \bar\alpha_i=\tfrac13\operatorname{Tr} \big(\textstyle\sum_j\boldsymbol{\bar\alpha}_{ij}\big),\qquad \boldsymbol{\bar\alpha}=(\boldsymbol\alpha^{-1}+\mathbf T_\text{GG})^{-1} \\ \partial\boldsymbol{\bar\alpha}= -\boldsymbol{\bar\alpha}( \partial\boldsymbol\alpha^{-1}+\partial\mathbf T_\text{GG} )\boldsymbol{\bar\alpha},\qquad \frac{\partial\bar\alpha_i}{\partial X_j}= -\frac13\sum_{\zeta\eta}\big( B_{i\zeta,j\eta}\bar\alpha'_{\zeta,j\eta}+ B'_{j\eta,\zeta}\bar\alpha_{j\eta,i\zeta} \big) \\ \mathbf B=\boldsymbol{\bar\alpha}\mathbf A, \quad A_{i\zeta,j\eta}= \frac{\partial(\alpha_i^{-1})}{\partial X_i} \delta_{ij}\delta_{\zeta\eta}+ \frac{\partial T^\text{GG}_{i\zeta,j\eta}}{\partial X_i},\quad \bar\alpha'_{\zeta,p}=\sum_i\bar\alpha_{i\zeta,p},\quad B'_{p,\zeta}=\sum_iB_{p,i\zeta} \end{gathered}$$

$$x'=x\left(\frac{y'}y\right)^q,\qquad \partial x'=x\left( \frac{\partial x}x+ q\frac{\partial y'}{y'}- q\frac{\partial y}{y} \right)$$

$$\begin{gathered} \sigma_i(u)=\left(\frac13\sqrt{\frac2\pi}\alpha_i(u)\right)^{\frac13},\qquad \partial\sigma_i=\sigma_i\frac{\partial\alpha_i}{3\alpha_i} \\ \sigma_{ij}(u)=\sqrt{\sigma_i(u)^2+\sigma_j(u)^2},\qquad \partial\sigma_{ij}= \frac{\sigma_i\partial\sigma_i+\sigma_j\partial\sigma_j}{\sigma_{ij}} \end{gathered}$$

$$T_{ab}(\mathbf r)=\frac{\partial^2}{\partial r_a\partial r_b}\frac1r= \frac{-3r_ar_b+r^2\delta_{ab}}{r^5},\qquad \frac{\partial T_{ab}(\mathbf r)}{\partial r_c}=-3\left( \frac{r_a\delta_{bc}+r_b\delta_{ca}+r_c\delta_{ab}}{r^5}- \frac{5r_ar_br_c}{r^7} \right)$$

$$\begin{aligned} T^\text{GG}_{ab}(\mathbf r,\sigma)&= \frac{\partial^2}{\partial r_a\partial r_b}\frac{\operatorname{erf}(\zeta)}r= \big(\operatorname{erf}(\zeta)-\Theta(\zeta)\big)T_{ab}(\mathbf r)+ 2\zeta^2\Theta(\zeta)\frac{r_ar_b}{r^5} \\ \Theta(\zeta)&=\frac{2\zeta}{\sqrt\pi}\exp(-\zeta^2),\qquad \zeta=\frac r\sigma \\ \frac{\mathrm d T_{ab}^\text{GG}(\mathbf r,\sigma)}{\mathrm dr_c}&= 2\zeta\Theta(\zeta)\left(T_{ab}(\mathbf r)+(3-2\zeta^2)\frac{r_ar_b}{r^5}\right) \frac{\mathrm d\zeta}{\mathrm dr_c} \\ &+\big(\operatorname{erf}(\zeta)-\Theta(\zeta)\big) \frac{\partial T_{ab}(\mathbf r)}{\partial r_c}- 2\zeta^2\Theta(\zeta)\left( \frac13\frac{\partial T_{ab}(\mathbf r)}{\partial r_c}+ \frac{r_c\delta_{ab}}{r^5} \right) \\ \qquad\frac{\mathrm d\zeta}{\mathrm dr_c}&= \frac{r_c}{r\sigma}-\frac r{\sigma^2}\frac{\mathrm d\sigma}{\mathrm dr_c} \end{aligned}$$

$$T_{ab}^\text{erfc}(\mathbf r,\gamma) =-3\frac{r_ar_b}{r^5}C(r,\gamma)+\frac{\delta_{ab}}{r^3}B(r,\gamma)$$

$$\begin{aligned} \frac{\partial T_{ab}^\text{erfc}(\mathbf r,\gamma)}{\partial r_c} &=-\left( \frac{r_a\delta_{bc}+r_b\delta_{ca}}{r^5}- 5\frac{r_ar_br_c}{r^7} \right)C(r,\gamma)-3\frac{r_c\delta_{ab}}{r^5}B(r,\gamma) \\ &-\frac{r_ar_br_c}{r^6}\frac{\partial C(r,\gamma)}{\partial r}+\frac{r_c\delta_{ab}}{r^4}\frac{\partial B(r,\gamma)}{\partial r} \end{aligned}$$