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 Å

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: [6].