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 Å/e². 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 Å
Quadrupole	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*.

Energy Torque Gradient

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].

Electrostatics Induced Dipole