# Induced Dipole¶

## Energy¶

$$\mu$$ is the induced dipole in the external field E. The induced field due to the induced dipole is $$E^u=-T\mu$$, and the induced dipole is proportional to the total field $$E^t$$:

$\mu = \alpha E^t = \alpha(E+E^u),$

where $$\alpha$$ is the polarizability. Defining $$\tilde{T}=\alpha^{-1}+T$$, the induced dipole is the solution to the linear equation

(5)$\tilde{T}\mu = E.$

The polarization energy is given by

(6)$\begin{split}U &= -\mu E +\int_0^\mu d\mu\ \tilde{T}\mu \\ &= -\mu E +\frac{1}{2}\mu\tilde{T}\mu.\end{split}$

On the right-hand side of (6):

• the 1st term is the contribution from the external field;

• the 2nd term is the mutual and self polarization energy.

Finally, the polarization energy is

(7)$U = -\frac{1}{2}\mu E.$

With limited numerical precision, the solution to linear equation (5) cannot be fully precise:

(8)$\tilde{T}\mu = \epsilon + E.$

The gradient of the induced dipole can be written in

$\mu' = \tilde{T}^{-1}(\epsilon' + E' - \tilde{T}'\mu),$

and the polarization gradient is

$\begin{split}U' &= -\frac{1}{2} (E\mu' + \mu E') \\ &= -\frac{1}{2} [(-\epsilon+\mu\tilde{T})\tilde{T}^{-1}(\epsilon' +E' -\tilde{T}'\mu) +\mu E'] \\ &\approx -\frac{1}{2} (\mu E' -\mu\tilde{T}'\mu +\mu E'),\end{split}$

only if the convergence of (8) is tight that $$\epsilon$$ and $$\epsilon'$$ terms will drop.

Tinker uses the following Conjugate Gradient algorithm (C.G.) with a sparse matrix preconditioner (denoted as M) [2] to obtain the induced dipoles. Related Tinker variables and routines are tabulated.

C.G. Terms

Tinker variables and routines

$$\gamma$$

a

$$\beta$$

b

$$r$$

rsd

$$M r$$

zrsd

$$p$$

conj

$$\tilde{T} p$$

vec

$$-T$$

ufield()

$$M$$

uscale()

## Polarization Model: AMOEBA (Thole Damping 2)¶

AMOEBA force field adopts two polarization schemes, d and p, for the external field due to the permanent multipoles, and a third scheme u for mutual induced dipole interactions. Both d and u schemes are group-based. The p scheme is atomic connectivity-based. Tinker uses C.G. iterations to solve the following linear equations

$\begin{split}(1/\alpha+T^u)\mu_d &= E_d \\ (1/\alpha+T^u)\mu_p &= E_p,\end{split}$

and defines the polarization energy as

(9)$U = -\frac{1}{2}\mu_d E_p.$

From an optimizational perspective, (9) is the minimum of the target function

$f_1(\mu_d,\mu_p)=\frac{1}{2}\left(\frac{1}{2}\mu_d\tilde{T}\mu_p +\frac{1}{2}\mu_p\tilde{T}\mu_d -E_d\mu_p-E_p\mu_d\right),$

whereas the way C.G. coded in Tinker is to solve the minimum of another target function

$f_2(\mu_d,\mu_p)=\frac{1}{2}\left(\frac{1}{2}\mu_d\tilde{T}\mu_d +\frac{1}{2}\mu_p\tilde{T}\mu_p -E_d\mu_d-E_p\mu_p\right).$

The difference in two target functions is usually negligible unless other loose convergence methods are used to compute the induced dipoles.

In the Thole damping model, a charge distribution $$\rho$$ is used as a replacement for the point dipole model. AMOEBA adopts the second functional form

$\rho = \frac{3a}{4\pi}\exp(-au^3)$

from paper [3], where u is the polarizability-scaled distance. The electrostatic field and potential at distance r can be obtained from Gauss’s law,

$E(r) = -\phi'(r) = \frac{1}{r^2} \int_0^u du\ 4\pi u^2 \rho = \frac{1-\exp(-au^3)}{r^2},$
$\phi(r) = \int_r^\infty dr\ E(r) = \frac{\lambda_1}{r} = \frac{1}{r}\left[ 1-\frac{(au^3)^\frac{1}{3}}{3}\Gamma(-\frac{1}{3},au^3)\right],$

where $$\lambda_1$$ serves as the $$B_0$$ term in EWALD quadrupole interactions. $$\lambda_n$$ terms are also related via derivatives

$\begin{split}\phi'' &= \frac{1}{r^3}\left[2-(2+3au^3)\exp(-au^3)\right], \\ \phi''' &= \frac{3}{r^4}\left[-2+(2+2au^3+3a^2u^6)\exp(-au^3)\right], \\ \phi'''' &= \frac{3}{r^5}\left[8-(8+8au^3+9a^3u^9)\exp(-au^3)\right],\end{split}$
$\begin{split}\phi'_i &= \phi'\frac{r_i}{r}, \\ \phi''_{ij} &= \left(\phi''-\frac{\phi'}{r}\right)\frac{r_i r_j}{r^2} +\frac{\phi'}{r}\delta_{ij}, \\ \phi'''_{ijk} &= \left(\phi'''-\frac{3\phi''}{r}+\frac{3\phi'}{r^2}\right)\frac{r_i r_j r_k}{r^3} +\left(\frac{\phi''}{r}-\frac{\phi'}{r^2}\right)\frac{\sum r_k \delta_{ij}}{r}, \\ \phi''''_{ijkl}&= \left(\phi''''-\frac{6\phi'''}{r}+\frac{15\phi''}{r^2}-\frac{15\phi'}{r^3}\right) \frac{r_i r_j r_k r_l}{r^4} \\ &+\left(\frac{\phi'''}{r}-\frac{3\phi''}{r^2}+\frac{3\phi'}{r^3}\right)\frac{\sum r_k r_l\delta_{ij}}{r^2} +\left(\frac{\phi''}{r^2}-\frac{\phi'}{r^3}\right)\sum\delta_{kl}\delta_{ij}.\end{split}$

Thus,

$\begin{split}-\lambda_3/r^3 &= \phi'/r \Rightarrow \\ &\lambda_3 = 1 - \exp(-au^3), \\ 3\lambda_5/r^5 &= (\phi''-\phi'/r)/r^2 \Rightarrow \\ &\lambda_5 = 1 - (1+au^3)\exp(-au^3), \\ -15\lambda_7/r^7 &= (\phi'''-3\phi''/r+3\phi'/r^2)/r^3 \Rightarrow \\ &\lambda_7 = 1 - \left(1+au^3+\frac{3}{5}a^2 u^6\right)\exp(-au^3), \\ 105\lambda_9/r^9 &= (\phi''''-6\phi'''/r+15\phi''/r^2-15\phi'/r^3)/r^4 \Rightarrow \\ &\lambda_9 = 1 - \left(1+au^3+\frac{18}{35}a^2 u^6+\frac{9}{35}a^3 u^9\right)\exp(-au^3).\end{split}$