\subsection{Electrostatic interaction energy}
\label{sec:distributed_multipoles}
\index{distributed multipoles}
\index{site energy!distributed multipoles}

We represent the molecular charge density by choosing multiple expansion sites (``polar sites'') per molecule in such a way as to accurately reproduce the molecular electrostatic potential (ESP), with a set of suitably chosen multipole moments $\{Q_{lk}^a\}$ (in spherical-tensor notation) allocated to each site. The expression for the electrostatic interaction energy between two molecules $A$ and $B$ in the multi-point expansion includes an implicit sum over expansion sites $a\epsilon A$ and $b\epsilon B$,
\begin{align}
 U_{AB} & = \sum_{a\epsilon A} \sum_{b\epsilon B} \hat{Q}_{l_1k_1}^a T_{l_1k_1l_2k_2}^{a,b} \hat{Q}_{l_2k_2}^b \equiv  \hat{Q}_{l_1k_1}^a T_{l_1k_1l_2k_2}^{a,b} \hat{Q}_{l_2k_2}^b,
 \label{equ:mol_distributed_U}
\end{align}
where we have used the Einstein sum convention for the site indices $a$ and $b$ on the right-hand side of the equation, in addition to the sum convention that is in place for the multipole-moment components $t\equiv l_1 k_1$ and $u\equiv l_2 k_2$. The $T_{l_1k_1l_2k_2}^{a,b}$ are tensors that mediate the interaction between a multipole component $l_1 k_1$ on site $a$ with the moment $l_2 k_2$ on site $b$. If we include the molecular environment into a perturbative term $W$ to enter in the single-molecule Hamiltonian, the above expression is exactly the first-order correction to the energy where the quantum-mechanical detail has been absorbed in classical multipole moments.

The are a number of strategies how to arrive at such a collection of {\em distributed multipoles}. They can be classified according to whether the multipoles are derived (a) from the electrostatic potential generated by the SCF charge density or (b) from a decomposition of the wavefunction itself. Here, we will only draft two of those approaches, CHELPG~\cite{breneman_determining_1990} from category (a) and DMA~\cite{stone_distributed_1985} from category (b).

The CHELPG (CHarges from ELectrostatic Potentials, Grid-based) method relies on performing a least-squares fit of atom-placed charges to reproduce the electrostatic potential as evaluated from the SCF density on a regularly spaced grid~\cite{breneman_determining_1990}. The fitted charges result from minimizing the Lagrangian function~\cite{chirlian_atomic_1987}
\begin{align}
 z(\{q_i\}) = \sum_{k=1}^M \left( \phi(\vec{r}_k) - \sum_{i=1}^N \frac{1}{4\pi\varepsilon_0} \frac{q_i}{|\vec{r}_i-\vec{r}_k|} \right) + \lambda \left( q_\textrm{mol} - \sum_{i=1}^N q_i \right),
\end{align}
with $M$ grid points, $N$ atomic sites, the set of atomic partial charges $\{q_i\}$ and the SCF potential $\phi$. The Lagrange multiplier $\lambda$ constrains the sum of the fitted charges to the molecular charge $q_\textrm{mol}$. The main difference from other fitting schemes~\cite{singh_approach_1984} is the algorithm that selects the positions at which the potential is evaluated (we note that the choice of grid points can have substantial effects especially for bulky molecules).
Clearly, the CHELPG method can be (and has been) extended to include higher atomic multipoles. It should be noted, however, how already the inclusion of atomic dipoles hardly improves the parametrization, and can in fact be harmful to its conformational stability.

The Distributed-Multipole-Analysis (DMA) approach~\cite{stone_distributed_1985, stone_distributed_2005}, developed by A. Stone, operates directly on the quantum-mechanical density matrix, expanded in terms of atom- and bond-centered Gaussian functions $\chi_\alpha = R_{LK}(\vec{x}-\vec{s}_\alpha) \exp[-\zeta(\vec{x}-\vec{s}_\alpha)^2]$,
\begin{align}
 \rho(\vec{x}) = \sum_{\alpha,\beta} \rho_{\alpha\beta} \chi_\alpha(\vec{x}-\vec{s}_\alpha) \chi_\beta(\vec{x}-\vec{s}_\beta). 
\end{align}
The aim is to compute multipole moments according in a distributed fashion: If we use that the overlap product $\chi_\alpha \chi_\beta$ of two Gaussian basis functions yields itself a Gaussian centered at $\vec{P} = (\zeta_\alpha \vec{s}_\alpha + \zeta_\beta \vec{s}_\beta) / (\zeta_\alpha + \zeta_\beta)$, it is possible to proceed in two steps: First, we compute the multipole moments associated with a specific summand in the density matrix, referred to the overlap center $\vec{P}$:
\begin{align}
 Q_{LK}[\vec{P}] = - \int R_{LK}(\vec{x}-\vec{P}) \rho_{\alpha\beta} \chi_\alpha \chi_\beta d^3\!x.
\end{align}
Second, we transfer the resulting $Q_{lk}[\vec{P}]$ to the position $\vec{S}$ of a polar site according to the rule~\cite{stone_distributed_1985}
\begin{align}
 Q_{nm}[\vec{S}] = \sum_{l=0}^L \sum_{k=-l}^l \left[ \left(\begin{array}{c} n+m \\ l+k \end{array}\right)\left(\begin{array}{c} n-m \\ l-k \end{array}\right) \right]^{1/2} R_{n-l,m-k}(\vec{S}-\vec{P})\cdot Q_{lk}[\vec{P}].
\end{align}
Note how this requires a rule for the choice of the expansion site to which the multipole moment should be transferred. In the near past~\cite{stone_distributed_2005}, the nearest-site algorithm, which allocates the multipole moments to the site closest to the overlap center, was replaced for diffuse functions by an algorithm based on a sxtpth weighting function in conjunction with grid-based integration methods in order to decrease the basis-set dependence of the resulting set of distributed multipoles.

One important advantage of the DMA approach over fitting algorithms such as CHELPG or Merz-Kollman (MK) is that higher-order moments can also be derived without too large an ambiguity.

The `mps' file format used by VOTCA for the definition of distributed multipoles (as well as point polarizabilities, see subsequent section) is based on the GDMA punch format of A. Stone's GDMA program~\cite{stone_distributed_2005} (the punch output file can be immediately plugged into VOTCA without any conversions to be applied). Furthermore the log-file of different QM packages (currently \gaussian, \turbomole and \nwchem) may be fed into the \toolref{log2mps} \tool, which will subsequently generate the appropriate mps-file.

\votcacommand{Read in ESP charges from a QM log file}{\cmdlogmps}