# Permanent Multipole¶

## Definitions and Units¶

The electrostatic potential at r due to the charge distribution nearby is

(1)$\phi(r) = \frac{1}{4\pi\epsilon_0} \int ds\frac{\rho(s)}{|r-s|},\ ds=\mathrm{d}x \mathrm{d}y \mathrm{d}z.$

Tinker uses a variable electric (in chgpot module) to represent the the factor $$1/(4\pi\epsilon_0)$$. Its default magnitude is 332.063713, which is a constant defined by variable coulomb (in units module), and its units are kcal/mol Å/e2. The default value can be modified by the ELECTIRIC keyword.

Note

Should the value of coulomb documented here be out-dated and become inconsistent with our code, please send us a pull request.

Expanding $$1/|r-s|$$ in Taylor series, $$4\pi\epsilon_0\phi(r)$$ can be rewritten as

(2)$\left[\int ds\rho(s)\right]\frac{1}{r} +\sum_i\left[\int ds\rho(s)s_i\right]\nabla_i\frac{1}{r} +\sum_{ij}\left[\frac{1}{2}\int ds\rho(s)s_i s_j\right]\nabla_i\nabla_j\frac{1}{r} +\cdots,$

where three pairs of square brackets give a set of definitions of monopole (charge, C), dipole (D), and quadrupole moments (Q*), respectively. The units of the multipole moments used in Tinker parameter files and internal calculation are different.

Multipole

Parameter Units

Internal Units

Charge

e

e

Dipole

e Bohr

e Å

e Bohr2

e Å2

In addition to different units, the quadrupole moments in Tinker parameter files use what is traditionally called traceless quadrupole $$\Theta$$ that has a different definition than Q*. The third term in (2) can be rewritten as

$\sum_{ij}\left[\frac{1}{2}\int ds\rho(s)(3s_i s_j - s^2\delta_{ij})\right] \frac{r_i r_j}{r^5},$

hence the traceless quadrupole can be defined as

$\Theta_{ij} = \frac{1}{2}\int ds\rho(s)(3s_i s_j - s^2\delta_{ij}).$

It is easy to confirm that $$\sum_k^{x,y,z}(3 s_k s_k - s^2)=0$$, thus,

$\Theta_{ij} = 3Q_{ij}^* - \delta_{ij}\sum_k^{x,y,z}Q_{kk}^*.$

Internally, Tinker scales $$\Theta$$ by 1/3

$Q = \Theta/3,$

so that the energy expression is the same as if we were using Q*.

Potential energy

(3)$U = \frac{1}{4\pi\epsilon_0}\int ds\rho(s)\phi(s).$

Potential energy with discretized charge distribution in (3)

(4)$U(r) = \phi(r) C(r) + \phi'(r) D(r) + \phi''(r) Q(r) + \cdots.$

Distance

$(r_x,r_y,r_z)=\boldsymbol{r}=r_2-r_1.$

Pairwise (atoms 1 and 2) quadrupole energy

$U_{12} = M_1^T T_{12} M_2.$

Multipoles

$\begin{split}M_1 = \begin{pmatrix} C_1 \\ D_1 \\ Q_1 \end{pmatrix},\ M_2 = \begin{pmatrix} C_2 \\ D_2 \\ Q_2 \end{pmatrix}.\end{split}$

T matrix

$\begin{split}T_{12} = \begin{pmatrix} 1 & \nabla_2 & \nabla_2^2 \\ \nabla_1 & \nabla_1\nabla_2 & \nabla_1\nabla_2^2 \\ \nabla_1^2 & \nabla_1^2\nabla_2 & \nabla_1^2\nabla_2^2 \end{pmatrix}\frac{1}{r}.\end{split}$

The upper left 4×4 elements of $$T_{12}$$

$\begin{split}T_{12}^{4 \times 4} = \begin{pmatrix} 1/r & -r_x/r^3 & -r_y/r^3 & -r_z/r^3 \\ r_x/r^3 & -3r_x^2/r^5 +1/r^3 & -3r_xr_y/r^5 & -3r_xr_z/r^5 \\ r_y/r^3 & -3r_xr_y/r^5 & -3r_y^2/r^5 +1/r^3 & -3r_yr_z/r^5 \\ r_z/r^3 & -3r_xr_z/r^5 & -3r_yr_z/r^5 & -3r_z^2/r^5 +1/r^3 \end{pmatrix}.\end{split}$

In the EWALD summation, $$1/r^k$$ terms will have different forms (Bn). Neverthelss, they are still connected through derivatives.

Non-EWALD

EWALD

$$1/r$$

$$B_0=\mathrm{erfc}(\alpha r)/r$$

$$1/r^3$$

$$B_1$$

$$3/r^5$$

$$B_2$$

$$15/r^7$$

$$B_3$$

$$105/r^9$$

$$B_4$$

$$945/r^{11}$$

$$B_5$$

The Bn terms are related to the (complementary) Boys functions and (complementary) error functions. For $$x>0$$ and $$n\ge 0$$,

$\begin{split}\frac{\mathrm{erf}(x)}{x} &= \frac{2}{\sqrt{\pi}}F_0(x^2), \\ \frac{\mathrm{erfc}(x)}{x} &= \frac{2}{\sqrt{\pi}}F_0^C(x^2), \\ F_n(x) &= \int_0^1 \exp(-xt^2)t^{2n} dt, \\ F_n^C(x) &= \int_1^\infty \exp(-xt^2)t^{2n} dt.\end{split}$

The Boys functions can be generated through upward and downward recursions

$\begin{split}F_n(x) = \frac{2xF_{n+1}(x) + \exp(-x)}{2n+1}, \\ F_n^C(x) = \frac{2xF_{n+1}^C(x) - \exp(-x)}{2n+1}.\end{split}$

Energy, torque, and force

Terms

Energy

Torque

Force

C

$$\phi C$$

N/A

$$\phi' C$$

D

$$\phi' D$$

$$\phi'\times D$$

$$\phi'' D$$

Q

$$\phi'' Q$$

$$\phi''\times Q$$

$$\phi''' Q$$

$\begin{split}\tau(D) &= \phi'\times D = D \times E, \\ \tau_{i}(Q) &= -2 \sum_{jk}\sum_{l}^{xyz} \epsilon_{ijk}Q_{jl}\phi''_{kl},\end{split}$

where $$\epsilon_{ijk}$$ is the Levi-Civita symbol.

Reference: [1].