Tables for
Volume B
Reciprocal space
Edited by U. Shmueli

International Tables for Crystallography (2010). Vol. B, ch. 3.5, pp. 471-474

Section 3.5.3. Generalization to Gaussian- and Hermite-based continuous charge distributions

T. A. Dardena*

aLaboratory of Structural Biology, National Institute of Environmental Health Sciences, 111 T. W. Alexander Drive, Research Triangle Park, NC 27709, USA
Correspondence e-mail:

3.5.3. Generalization to Gaussian- and Hermite-based continuous charge distributions

| top | pdf |

In this section we develop the generalization of Ewald sums to interacting molecular densities modelled by linear combinations of spherical or Hermite Gaussians. In Coppens (1997[link]) the electrostatic potential and energy for crystal densities modelled by linear combinations of Slater-like density elements is discussed. As with quantum-chemical programs, there are trade-offs in the choice of density elements. The use of Slater-like densities provides accuracy with the use of fewer elements, but the Gaussian and Hermite Gaussian elements are far more convenient analytically. In particular the Ewald summation, as well as the fast-Fourier-transform-based methods discussed in Section 3.5.4[link], generalize straightforwardly to these elements. We note ahead of time that the case of lattice sums of point multipoles are obtained by specializing the results for Hermite Gaussians, just as the results for spherical Gaussians specialize to the above lattice summation results for point charges.

A somewhat different approach to Ewald summation over Hermite Gaussian elements, focused on application to periodic density-functional calculations, is given in Trickey et al. (2004[link]). The continuous fast multipole method (Challacombe et al., 1997[link]) provides yet a different but highly efficient approach to lattice summation of continuous charge densities. Lattice sums of interacting spherical Gaussians

| top | pdf |

We begin by discussing the electrostatic interaction energy between a normalized spherical Gaussian charge distribution [\rho_1], centred at a point [{\bf R}_1 \in U] and having exponent [\theta_1], and a second normalized spherical Gaussian charge distribution [\rho_2], centred at [{\bf R}_2 \in U ] and having exponent [\theta_2], together with the latter's periodic images, centred at [{\bf R} + {\bf n}], [{\bf n} \in \Omega(S,K) ]. Their interaction energy can be written[\eqalignno{E_{12} &= \textstyle\sum\limits_{{\bf n} \in \Omega(S,K)} \textstyle\int\limits_{\Re^3}\textstyle\int\limits_{\Re^3} \rho_1({\bf r}_1 - {\bf R}_1) \rho_2({\bf r}_2 -{\bf R}_2 - {\bf n})&\cr &\quad\times {{1}/{|{\bf r}_2 - {\bf r}_1|}}\,{\rm d}^3 {\bf r}_1 \,{\rm d}^3 {\bf r}_2,&(} ]where [\rho_j({\bf r}) = (\theta_j/\pi)^{3/2} \exp(-\theta_j {\bf r}^2) ], j = 1, 2. We now recall equation ([link], applying its two right-hand components separately, to get [E_{12} = E_{12}^{\,\rm d} + E_{12}^{\,\rm r} ], where[\eqalign{E_{12}^{\,\rm d} &= \sum\limits_{{\bf n} \in \Omega(S,K)} \int\limits_{\Re^3} \int\limits_{\Re^3} \rho_1({\bf r}_1 - {\bf R}_1) \rho_2({\bf r}_2 -{\bf R}_2 - {\bf n}) \cr &\quad\times {{{\rm erfc}(\alpha^{1/2}|{\bf r}_2 - {\bf r}_1|)}\over{|{\bf r}_2 - {\bf r}_1|}} \,{\rm d}^3{\bf r}_1 \,{\rm d}^3{\bf r}_2} ]and[\eqalign{E_{12}^{\,\rm r} &= {{1}\over{\pi}} \sum\limits_{{\bf n} \in \Omega(S,K)} \int\limits_{\Re^3} \int\limits_{\Re^3} \rho_1({\bf r}_1 - {\bf R}_1) \rho_2({\bf r}_2 -{\bf R}_2 - {\bf n}) \cr &\quad \times \int\limits_{\Re^3} {{\exp(-\pi^2 {\bf u}^2 / \alpha)}\over{{\bf u}^2}} \exp[-2 \pi i {\bf u}\cdot ( {\bf r}_2 - {\bf r}_1)]\,{\rm d}^3 {\bf u} \, {\rm d}^3{\bf r}_1 \,{\rm d}^3{\bf r}_2.} ]

The direct pair energy [E_{12}^{\,\rm d}] involves the damped Coulomb interaction, which has been discussed by Gill & Adamson (1996[link]). We can calculate [E_{12}^{\,\rm d}] using the following standard results from quantum chemistry, which are derived in Helgaker et al. (2000[link]): Let [\rho_p({\bf r})], [\rho_q({\bf r}) ] be two normalized spherical Gaussian charge distributions having exponents p and q and centred at points P and Q, i.e. [\rho_p({\bf r}) = (p / \pi)^{3/2} \exp[-p({\bf r} - {\bf P})^2] ] and [\rho_q({\bf r})] = [(q / \pi)^{3/2} \exp[-q({\bf r} - {\bf Q})^2] ]. Then the electrostatic potential [V_p({\bf R})] at a point R due to the charge distribution [\rho_p({\bf r})] is given by[V_p({\bf R}) = \left({{4 p}\over{\pi}}\right)^{1/2} F_0[p ({\bf P} - {\bf R})^2], \eqno( ]where [F_0] is the zeroth-order Boy's function, with [F_0(x) =] [ \textstyle\int_0^1 \exp(-x t^2) \,{\rm d}t ] for x > 0. The Coulomb interaction energy [E_{pq}] between [\rho_p] and [\rho_q] is given by[\eqalignno{E_{pq} &= \int\limits_{\Re^3} \int\limits_{\Re^3} {{\rho_p({\bf r}_1 - {\bf P}) \rho_q({\bf r}_2 - {\bf Q})}\over{|{\bf r}_1 - {\bf r}_2|}} \,{\rm d}^3{\bf r}_1 \,{\rm d}^3{\bf r}_2 &\cr &= \int\limits_{\Re^3} V_p({\bf r}_2) \rho_q({\bf r}_2 - {\bf Q}) \,{\rm d}^3{\bf r}_2&\cr&= \left({{4 \mu_{pq}}\over{\pi}}\right)^{1/2} F_0[\mu_{pq} ({\bf P} - {\bf Q})^2],&(} ]where [ 1 / \mu_{pq} = 1 / p + 1 / q]. The Boy's function can be related to the error function:[{{2}\over{\pi^{1/2}}} F_0(x^2) = {{{\rm erf}(|x|)}\over{|x|}}.\eqno( ]Substituting equation ([link] into ([link], we see that[V_p({\bf R}) = {{{\rm erf}(p^{1/2} |{\bf R} - {\bf P}|) }\over{ |{\bf R} - {\bf P}|}},\eqno( ]a result which we had derived above by another method [see equation ([link]]. From equations ([link] and ([link], substituting [\alpha] for p and [\rho_p ] for [\rho_q], we have for any point [{\bf r} \in \Re^3][\int\limits_{\Re^3} {{ {\rm erf}(\alpha^{1/2} |{\bf r} - {\bf r}_1|)}\over{|{\bf r} - {\bf r}_1|}} \rho_p({\bf r}_1 - {\bf P})\,{\rm d}^3{\bf r}_1 = {{{\rm erf}(\mu_{p\alpha}^{1/2}\ |{\bf r} - {\bf P}|)}\over{|{\bf r} - {\bf P}|}}, ]where [1/\mu_{p \alpha} = 1 / p + 1 / \alpha]. Iterating this result with [\mu_{p \alpha}] replacing [\alpha] and [\rho_q ] replacing [\rho_p] we get the damped interaction energy between [\rho_p] and [\rho_q],[\eqalignno{ E_{pq}^{\,\rm d} &= \int\limits_{\Re^3} \int\limits_{\Re^3} \rho_p({\bf r}_1 - {\bf P}) \rho_q({\bf r}_2 - {\bf Q}) &\cr &\quad\times {{{\rm erfc}(\alpha^{1/2}|{\bf r}_2 - {\bf r}_1|)}\over{|{\bf r}_2 - {\bf r}_1|}} \,{\rm d}^3{\bf r}_1 \,{\rm d}^3{\bf r}_2&\cr =& \int\limits_{\Re^3} \int\limits_{\Re^3} {{\rho_p({\bf r}_1 - {\bf P}) \rho_q({\bf r}_2 - {\bf Q})}\over{|{\bf r}_1 - {\bf r}_2|}} \,{\rm d}^3{\bf r}_1 \,{\rm d}^3{\bf r}_2 &\cr&\quad- \int\limits_{\Re^3} \int\limits_{\Re^3} \rho_p({\bf r}_1 - {\bf P}) \rho_q({\bf r}_2 - {\bf Q}) &\cr &\quad\times {{{\rm erf}(\alpha^{1/2}|{\bf r}_2 - {\bf r}_1|)}\over{|{\bf r}_2 - {\bf r}_1|}} \,{\rm d}^3{\bf r}_1 \,{\rm d}^3{\bf r}_2&\cr &= {{{\rm erf}(\mu_{pq}^{1/2} |{\bf P} - {\bf Q}|)}\over{|{\bf P} - {\bf Q}|}} - {{{\rm erf}(\mu_{pq\alpha}^{1/2} |{\bf P} - {\bf Q}|)}\over{|{\bf P} - {\bf Q}|}} &\cr &= {{{\rm erfc}(\mu_{pq\alpha}^{1/2} |{\bf P} - {\bf Q}|)}\over{|{\bf P} - {\bf Q}|}} - {{{\rm erfc}(\mu_{pq}^{1/2} |{\bf P} - {\bf Q}|)}\over{|{\bf P} - {\bf Q}|}},&(} ]where [1 / \mu_{pq} = 1 / p + 1 / q] and [1 / \mu_{pq\alpha} = 1 / p + 1 / q + 1 / \alpha ]. Notice that in the limit [\alpha \to \infty], [E_{pq,\alpha} \to 0 ] as expected, since [{\rm erfc}(x) \to 0] as [x \to \infty ], and that in the limit [p \to \infty], [q \to \infty], [E_{pq,\alpha} \to {\rm erfc}(\alpha^{1/2}|{\bf P} - {\bf Q}|) / |{\bf P} - {\bf Q}| ], the damped interaction between two unit point charges at P and Q, as in the standard Ewald direct sum.

To calculate the reciprocal pair energy [E_{12}^{\,\rm r}] we first apply the first right-hand component of equation ([link] together with the Fourier transforms of the normalized Gaussians [\rho_p ] and [\rho_q] from equation ([link] to obtain[\eqalignno{E_{pq}^{\,\rm r} &= {{1}\over{\pi}} \int\limits_{\Re^3} \int\limits_{\Re^3} \int\limits_{\Re^3} \rho_p({\bf r}_1 - {\bf P}) \rho_q({\bf r}_2 - {\bf Q}) {{\exp(-\pi^2 {\bf u}^2 / \alpha)}\over{{\bf u}^2}}&\cr &\quad \times \exp[-2 \pi i {\bf u}\cdot ({\bf r}_1 - {\bf r}_2)] \,{\rm d}^3 {\bf u} \,{\rm d}^3{\bf r}_1 \,{\rm d}^3{\bf r}_2 &\cr&= {{1}\over{\pi}} \int\limits_{\Re^3} \int\limits_{\Re^3} \int\limits_{\Re^3} {{1}\over{{\bf u}^2}} \exp[-\pi^2 ({\bf u}^2/\alpha + {\bf v}^2/p + {\bf w}^2/q)] &\cr&\quad\times \exp(2 \pi i {\bf v} \cdot {\bf P}) \exp(2 \pi i {\bf w} \cdot {\bf Q}) \,{\rm d}^3{\bf u} \,{\rm d}^3{\bf v}\,{\rm d}^3{\bf w}&\cr&\quad\times \int\limits_{\Re^3} \exp[-2 \pi i ({\bf u}+{\bf v})\cdot{\bf r}_1]\,{\rm d}^3{\bf r}_1&\cr &\quad\times \int\limits_{\Re^3}\exp[2 \pi i ({\bf u}-{\bf w})\cdot{\bf r}_2] \,{\rm d}^3{\bf r}_2 &\cr &= {{1}\over{\pi}} \int\limits_{\Re^3} {{ \exp(-\pi^2 {\bf u}^2 / \mu_{pq\alpha}) }\over{{\bf u}^2}} \exp[-2 \pi i {\bf u} \cdot ({\bf P} - {\bf Q})] \,{\rm d}^3{\bf u},&\cr &&(} ]where we have used the fact that [\textstyle\int \exp[-2 \pi i ({\bf u}+{\bf v})\cdot{\bf r}_1] \,{\rm d}^3{\bf r}_1] = [\delta^3({\bf u}+{\bf v}) ], with [\delta^3] the three-dimensional delta function, and similarly for the integral over [{\bf r}_2]. Note that in the limit [p \to \infty ], [q \to \infty], [E_{pq}^{\,\rm r}] reduces to the equivalent expression for unit point charges at P and Q, and in the limit [\alpha \to \infty], [E_{pq}^{\,\rm r}] reduces to the Fourier-space expression for the direct Coulomb inter­action between the Gaussians [\rho_p ] and [\rho_q].

Recalling the developments leading to equation ([link], and using equations ([link] and ([link], we can write the inter­action energy [E_{12}] from equation ([link] as[\eqalignno{E_{12} &= \sum_{\bf n} {{{\rm erfc}(\mu_{12\alpha}^{1/2} |{\bf R}_{12} - {\bf n}|) - {\rm erfc}(\mu_{12}^{1/2} |{\bf R}_{12} - {\bf n}|)}\over{|{\bf R}_{12} - {\bf n}|}}&\cr &\quad+ {{1}\over{\pi V}} \sum_{{\bf m} \ne {\bf 0}} {{ \exp(-\pi^2 {\bf m}^2 / \mu_{12\alpha}) }\over{{\bf m}^2}}\exp[-2 \pi i {\bf m} \cdot ({\bf R}_{12})] &\cr&\quad- {{\pi}\over{\mu_{12\alpha} V}} + {{1}\over{\pi}} H_{S,K}({\bf R}_{12}) + {\varepsilon}(K),&(} ]where [{\bf R}_{12} = {\bf R}_1 - {\bf R}_2], [1 / \mu_{12} = 1 / \theta_1 + 1 / \theta_2 ] and [1 / \mu_{12\alpha}= 1 / \theta_1 ] [+ \ 1 / \theta_2 + 1 / \alpha ].

Note that although the energy of a finite-exponent Gaussian charge distribution interacting with itself is finite, when considering a generalization of lattice sums involving point charges it is necessary to remove the direct Coulomb interaction of the Gaussian with itself. From equations ([link] and ([link], if p = q, then as [{\bf Q} \to {\bf P}], [E_{pq}^{\,\rm d} - E_{pq} \to -2 (\mu_{pq\alpha} / \pi)^{1/2} ]. We can now generalize equation ([link]. Let [\rho_1, \ldots,\rho_N ] be normalized spherical Gaussian charge distributions with Gaussian exponents [\theta_1,\ldots,\theta_N], centred at [\lbrace{\bf R}_1,\ldots,{\bf R}_N\rbrace = \underline{\bf R}^{\lbrace\bf N\rbrace}\in U ] (point charges are included by taking [\theta \to \infty]), and let [q_1,\ldots,q_N] satisfy [q_1+\ldots+q_N = 0] (neutral unit cell). Then the energy of the central unit cell [U] within a large spherical crystal, due to interactions of the Gaussian charge distributions [q_i \rho_i] with each other and with their periodic images within the crystal, is given by[\eqalign{E_{S,K}(\underline{\bf R}^{\lbrace\bf N\rbrace}) &= {\textstyle{{1}\over{2}}} \sum_{\bf n}{}^{'}\sum_{i,j=1}^N q_i q_j \left \lbrace {{{\rm erfc}(\mu_{ij\alpha}^{1/2} |{\bf R}_{ij} - {\bf n}|)}\over{|{\bf R}_{ij} - {\bf n}|}} \right. \cr &\quad- \left. \ {{{\rm erfc}(\mu_{ij}^{1/2} |{\bf R}_{ij} - {\bf n}|)}\over{|{\bf R}_{ij} - {\bf n}|}} \right \rbrace \cr &\quad+ {{1}\over{2 \pi V}} \sum_{{\bf m} \ne {\bf 0}} \sum_{i,j=1}^N q_iq_j {{ \exp(-\pi^2 {\bf m}^2 / \mu_{ij\alpha})}\over{{\bf m}^2}} \cr &\quad\times \exp(-2 \pi i {\bf m} \cdot {\bf R}_{ij}) \cr &\quad- {{\pi}\over{2 V}} \sum_{i,j = 1}^N {{q_i q_j}\over{\mu_{ij\alpha}}} - \sum_{i=1}^N q_i^2 \left({{\mu_{ii\alpha}}\over{\pi}}\right)^{1/2} \cr &+\quad {{2 \pi {\bf D}^2}\over{3 V}} + {\varepsilon}(K),} ]where [{\bf R}_{ij} = {\bf R}_i - {\bf R}_j], [1 / \mu_{ij} = 1 / \theta_i + 1 / \theta_j ], [1 / \mu_{ij\alpha}= 1 / \theta_i +  1 / \theta_j] [ +\ 1 / \alpha] and as before [{\bf D} = q_1{\bf R}_1 +\ldots + q_N {\bf R}_N] is the unit-cell dipole.

Note that from the above derivation, the [\alpha] in [\mu_{ij\alpha} ] is allowed to be different for each pair ij. One choice that leads to an efficient algorithm (particularly when combined with the fast-Fourier-transform-based methods discussed below) is to separate the Gaussian charge distributions [q_j \rho_j] into compact and diffuse sets, based on their exponents [\theta_j], and choose [\alpha] to be infinite for ij pairs where at least one of the two Gaussians is diffuse. These pair interactions are then evaluated entirely in reciprocal space. Next choose [\alpha] so that [\mu_{ij\alpha}] is constant for all compact pairs ij. More specifically, given [\beta > 0], a Gaussian [q_j \rho_j] is classified as compact if [\theta_j \ge 2 \beta ] and diffuse otherwise. We denote these cases as [j \in C] and [j \in D], respectively. Then for [i,j \in C], choose [\alpha ] so that [1/\mu_{ij\alpha}] = [1/\theta_i + 1/\theta_j + 1/\alpha = 1/\beta ]. Otherwise, choose [\alpha] to be infinite. Then the Coulomb energy of the central unit cell [E_{S,K}(\underline{\bf R}^{\lbrace\bf N\rbrace}) ] can be re-expressed as[\eqalignno{E_{S,K}(\underline{\bf R}^{\lbrace\bf N\rbrace}) &= {\textstyle{{1}\over{2}}} \sum_{\bf n}{}^{'}\sum_{(i,j)\in C\times C} q_i q_j \left \lbrace {{{\rm erfc}(\beta^{1/2} |{\bf R}_{ij} - {\bf n}|)}\over{|{\bf R}_{ij} - {\bf n}|}} \right. &\cr&\quad- \left. \ {{{\rm erfc}(\mu_{ij}^{1/2} |{\bf R}_{ij} - {\bf n}|)}\over{|{\bf R}_{ij} - {\bf n}|}} \right \rbrace &\cr &\quad+ {{1}\over{2 \pi V}} \sum_{{\bf m} \ne {\bf 0}}\ \sum_{(i,j)\in C\times C} q_iq_j \ {{ \exp(-\pi^2 {\bf m}^2 / \beta ) }\over{{\bf m}^2}}&\cr&\quad\times \exp(-2 \pi i {\bf m} \cdot {\bf R}_{ij})&\cr &\quad+ {{1}\over{2\pi V}} \sum_{{\bf m} \ne {\bf 0}}\ \sum_{(i,j) \not\in C\times C}q_iq_j \ {{ \exp(-\pi^2 {\bf m}^2 / \mu_{ij})}\over{{\bf m}^2}} &\cr &\quad\times \exp(-2 \pi i {\bf m} \cdot {\bf R}_{ij}) &\cr&\quad- {{\pi}\over{2 V}} \sum_{(i,j)\in C\times C} q_i q_j\left ({{1}\over{\beta}} - {{1}\over{\theta_i}} - {{1}\over{\theta_j}}\right) &\cr&\quad- \sum_{i\in C} q_i^2 \left({{\beta}\over{\pi}}\right)^{1/2} - \sum_{i \not\in C}q_i^2 \left({{\theta_i}\over{2\pi}}\right)^{1/2} + {{2 \pi {\bf D}^2}\over{3 V}} &\cr &\quad+ {\varepsilon}(K).&(} ]

Note that in the limit that [\theta_i \to \infty] for all i, all of the Gaussians become compact and this formula reduces to equation ([link] for a large spherical crystal of point charges. The volume-dependent term [({{\pi}/{2 V}} )\textstyle\sum_{(i,j)\in C\times C} q_i q_j ({{1}/{\beta}} - {{1}/{\theta_i}} - {{1}/{\theta_j}}) ] is somewhat surprising, and would not appear in a more straightforward Ewald derivation such as our first Ewald derivation. It vanishes by unit-cell neutrality for the case of point charges or multipoles, and would also vanish if we chose a single [\alpha] for all Gaussian pairs. Higher angular momentum: Hermite Gaussians

| top | pdf |

More complex crystal charge distributions can be realized by considering normalized Hermite Gaussian distributions in place of the normalized spherical Gaussian distributions discussed above. Let [\rho] be a normalized spherical Gaussian charge distribution with exponent [\theta], centred at [{\bf R} = ({\bf R}_x,{\bf R}_y,{\bf R}_z)], i.e. [\rho({\bf r}) = (\theta/\pi)^{3/2} \exp[-\theta({\bf r} - {\bf R})^2] ]. The associated normalized Hermite Gaussians [\Lambda_{tuv}], where t, u, v are non-negative integers, are defined by[\eqalign{\Lambda_{tuv}({\bf r},\theta,{\bf R}) &= \left({{\theta}\over{\pi}}\right)^{3/2} \left(\partial / \partial{\bf R}_x\right)^t \left(\partial / \partial{\bf R}_y\right)^u \left(\partial / \partial{\bf R}_z\right)^v \cr &\quad\times \exp[-\theta({\bf r} - {\bf R})^2].} ]Note that this definition of Hermite Gaussians differs from that found in some textbooks, e.g. Helgaker et al. (2000[link]), in that we have normalized the Gaussian [\rho] it is derived from in order to allow exponents to tend to infinity to arrive at ideal point charges or point multipoles. We will refer to these Hermite Gaussians as normalized Hermite Gaussians to emphasize this distinction. The Hermite definition reduces to the spherical Gaussian when t, u and v are all zero, i.e. [\Lambda_{000}({\bf r},\theta,{\bf R}) = \rho({\bf r})].

Above we discussed the interaction between spherical Gaussian charge distributions [\rho]. Here we generalize these to higher angular momentum charge distributions [\varrho] given by[\varrho({\bf r},\theta,{\bf R},{\bf c}) = \textstyle\sum\limits_{tuv} {\bf c}_{tuv}\Lambda_{tuv}({\bf r},\theta,{\bf R}), ]where c is a vector of coefficients, a finite number of which are nonzero. For example, if [{\bf c}_{000}] is the only nonzero coefficient, we are reduced to the spherical charge distributions discussed above.

The fact that the partial derivatives within [\Lambda_{tuv}] are with respect to the position of the centre of the distribution allows us to easily calculate the Coulomb interaction energy between two normalized Hermite Gaussian distributions. The partial derivatives can be pulled out of the double integral, so that the interaction energy is given by the corresponding partial derivatives of the interaction energy of the two spherical Gaussian charge distributions:[\eqalign{E_{12} &= \int\limits_{\Re^3} \int\limits_{\Re^3} {{\Lambda_{t_1u_1v_1}({\bf r}_1,\theta_1,{\bf R}_1) \Lambda_{t_2 u_2 v_2}({\bf r}_2,\theta_2,{\bf R}_2)}\over{|{\bf r}_1 - {\bf r}_2|}} \, {\rm d}^3{\bf r}_1 \,{\rm d}^3{\bf r}_2 \cr &= \left({{\partial }\over{\partial{\bf R}_{1x}}}\right)^{t_1} \left({{\partial}\over{ \partial{\bf R}_{1y}}}\right)^{u_1} \left({{\partial}\over{ \partial{\bf R}_{1z}}}\right)^{v_1} \left({{\partial }\over{ \partial{\bf R}_{2x}}}\right)^{t_2} \cr &\quad\times \left({{\partial}\over{ \partial{\bf R}_{2y}}}\right)^{u_2} \left({{\partial }\over{\partial{\bf R}_{2z}}}\right)^{v_2} \left({{\theta_1\theta_2}\over{\pi^2}}\right)^{3/2} \int\limits_{\Re^3} \int\limits_{\Re^3} {{1}\over{|{\bf r}_1 - {\bf r}_2|}} \cr &\quad\times \exp[-\theta_1({\bf r}_1 - {\bf R}_1)^2] \exp[-\theta_2({\bf r}_2 - {\bf R}_2)^2]\,{\rm d}^3{\bf r}_1\,{\rm d}^3{\bf r}_2.} ]

Next we consider the interaction energy between a generalized charge distribution [\varrho_1] centred at [{\bf R}_1 \in U] and another generalized charge distribution [\varrho_2] centred at [{\bf R}_2 \in U], together with all the latter's periodic images, centred at [{\bf R}_2 + {\bf n} ], [{\bf n} \in \Omega(S,K)]. Examining equation ([link], we note that the terms to be differentiated all depend on [{\bf R}_{12} = {\bf R}_1 - {\bf R}_2 ]. Thus[\eqalignno{E_{12} &= \sum_{\bf n} \sum_{t_1 u_1 v_1} \sum_{t_2 u_2 v_2} (-1)^{(t_2+u_2+v_2)} {\bf c}_{1,t_1u_1v_1} {\bf c}_{2,t_2u_2v_2} &\cr &\times \left({{\partial }\over{\partial{\bf R}_{12x}}}\right)^{t_1+t_2} \left({{\partial}\over{ \partial{\bf R}_{12y}}}\right)^{u_1+u_2} \left({{\partial}\over{ \partial{\bf R}_{12z}}} \right)^{v_1+v_2} &\cr&\quad\times \left \lbrace {{{\rm erfc}(\mu_{12\alpha}^{1/2} |{\bf R}_{12} - {\bf n}|)}\over {|{\bf R}_{12} - {\bf n}|}} - {{{\rm erfc}(\mu_{12}^{1/2} |{\bf R}_{12} - {\bf n}|)}\over {|{\bf R}_{12} - {\bf n}|}} \right \rbrace &\cr &\quad+ {{1}\over{\pi V}} \sum_{{\bf m} \ne {\bf 0}} {{ \exp(-\pi^2 {\bf m}^2 / \mu_{12\alpha}) }\over{{\bf m}^2}} C_1({\bf m}) C_2(-{\bf m})&\cr &\quad- {\bf c}_{1,000} {\bf c}_{2,000}{{\pi}\over{\mu_{p12\alpha} V}} + {{1}\over{\pi}} H_{S,K}^{12}({\bf R}_{12}) &\cr &\quad+ {\varepsilon}(K),&(} ]where [C_i({\bf m})], i =1, 2 are given by[\eqalign{C_i({\bf m}) &= \sum_{t_i u_i v_i}{\bf c}_{i,t_iu_iv_i} \cr &\quad\times \left({{\partial }\over{\partial{\bf R}_{ix}}}\right)^{t_i} \left({{\partial}\over{ \partial{\bf R}_{iy}}}\right)^{u_i} \left({{\partial}\over{ \partial{\bf R}_{iz}}} \right)^{v_i} \exp(-2 \pi i {\bf m}\cdot{\bf R}_i) } ]and[\eqalign{H_{S,K}^{12}({\bf R}_{12}) &= \sum_{t_1 u_1 v_1} \sum_{t_2 u_2 v_2} (-1)^{(t_2+u_2+v_2)} {\bf c}_{1,t_1u_1v_1} {\bf c}_{2,t_2u_2v_2} \cr &\quad\times \left({{\partial }\over{\partial{\bf R}_{12x}}}\right)^{t_1+t_2} \left({{\partial}\over{ \partial{\bf R}_{12y}}}\right)^{u_1+u_2} \left({{\partial}\over{ \partial{\bf R}_{12z}}} \right)^{v_1+v_2} \cr &\quad\times H_{S,K}({\bf R}_{12}).} ]Note that by the arguments leading to equation ([link], [H_{S,K}({\bf r}) ] is quadratic in r, i.e. derivatives of higher than second order vanish.

We now can discuss the Coulomb energy of the central unit cell U within a large spherical crystal when the charge density in the unit cell is described by a basis [\Theta] of normalized Hermite Gaussians and a set of expansion points [{\bf R}_1,\ldots,{\bf R}_N] within U that are repeated periodically in the crystal (for example, some of the expansion points could be atomic nuclei). Let [\theta_1,\ldots,\theta_L ] be positive Gaussian exponents (or infinite, in the case of nuclear charge or point charges or ideal multipoles) and let the charge density about the expansion point [{\bf R}_i] be given by[\varrho_i({\bf r},{\bf R}_i,\Theta) = \textstyle\sum\limits_{l = 1}^L \textstyle\sum\limits_{tuv} {\bf c}_{i,l,tuv} \Lambda_{tuv}({\bf r},\theta_l,{\bf R}_i), ]where, as above, for each i, l the coefficients [{\bf c}_{i,l,tuv}] are nonzero for a finite number of tuv. As above we can separate the exponents [\theta_l] into compact and diffuse, i.e. [l \in C] or [l \in D], according to whether [\theta_l \ge 2 \beta ] or [\theta \,\lt\, 2 \beta], respectively. Let [q_i = \textstyle\sum_{l = 1}^L {\bf c}_{i,l,000} ] denote the net charge at the expansion point [{\bf R}_i].

We are interested in the Coulomb energy of charge densities [\varrho_i({\bf r},{\bf R}_i,\Theta) ], [i = 1,\ldots,N], interacting with each other and their periodic images at expansion points [{\bf R}_i + {\bf n}], [{\bf n} \in \Omega(S,K) ] for large K. As above, the direct Coulomb interaction of [\varrho_i({\bf r},{\bf R}_i,\Theta)] with itself is not permitted. This can be modified for more sophisticated treatments e.g. in the context of the Coulomb integrals in a periodic density-functional-theory approach, but clearly a nuclear charge cannot interact with itself. The unit cell is assumed to be neutral, i.e. [\textstyle\sum_{i=1}^N q_i = 0]. Then, using equations ([link] and ([link] we can write the Coulomb energy of [\underline\varrho^{\lbrace\bf N\rbrace}] = [\lbrace \varrho_1,\ldots,\varrho_N\rbrace] within the spherical crystal as[\eqalignno{ E_{S,K}(\underline\varrho^{\lbrace\bf N\rbrace}) &= {\textstyle{{1}\over{2}}} \sum_{\bf n}{}^{'} \sum_{i=1}^N \sum_{l_i \in C} \sum_{t_i u_i v_i} {\bf c}_{i,l_i,t_i u_i v_i} &\cr &\quad\times \sum_{j=1}^N \sum_{l_j \in C} \sum_{t_j u_j v_j} (-1)^{(t_j+u_j+v_j)} {\bf c}_{j,l_j,t_j u_j v_j} &\cr &\quad\times \left({{\partial }\over{\partial{\bf R}_{ijx}}}\right)^{t_i+t_j} \left({{\partial}\over{ \partial{\bf R}_{ijy}}}\right)^{u_i+u_j} \left({{\partial}\over{ \partial{\bf R}_{ijz}}} \right)^{v_i+v_j} &\cr &\quad\times \left \lbrace {{{\rm erfc}(\beta^{1/2} |{\bf R}_{ij} - {\bf n}|)}\over {|{\bf R}_{ij} - {\bf n}|}} - {{{\rm erfc}(\mu_{l_i l_j}^{1/2} |{\bf R}_{ij} - {\bf n}|)}\over {|{\bf R}_{ij} - {\bf n}|}} \right\rbrace&\cr &\quad+ {{1}\over{2 \pi V}} \sum_{{\bf m} \ne {\bf 0}}{{1}\over{{\bf m}^2}} \exp(-\pi^2 {\bf m}^2 /2 \beta) \sum_{l_1\in C}S_{l_1}({\bf m}) &\cr &\quad\times \exp(-\pi^2 {\bf m}^2 /2 \beta) \sum_{l_2\in C} S_{l_2}(-{\bf m})&\cr &\quad+ {{1}\over{2 \pi V}} \sum_{{\bf m} \ne {\bf 0}} {{1}\over{{\bf m}^2}} \sum_{(l_1,l_2)\not\in C\times C} \exp(-\pi^2 {\bf m}^2 / \theta_{l_1}) &\cr &\quad\times \exp(-\pi^2 {\bf m}^2 / \theta_{l_2}) S_{l_1}({\bf m}) S_{l_2}(-{\bf m}) &\cr &\quad- {{\pi}\over{2 V}} \sum_{l_1 \in C} \sum_{l_2 \in C} \sum_{i=1}^N \sum_{j=1}^N {\bf c}_{i,l_1,000} {\bf c}_{j,l_2,000} &\cr &\quad\times \left({{1}\over{\beta}} - {{1}\over{\theta_{l_1}}} - {{1}\over{\theta_{l_2}}}\right) &\cr&\quad- \sum_{i=1}^N E_{\rm self}(\varrho_i) + {{2 \pi {\bf D}^2}\over{3 V}} + {\varepsilon}(K), &(} ]where the structure factors [S_{l}({\bf m})] are given by[\eqalignno{S_{l}({\bf m}) &= \sum_{i=1}^N \sum_{tuv} {\bf c}_{i,l,tuv} \left({{\partial }\over{\partial{\bf R}_{ix}}}\right)^t \left({{\partial }\over{\partial{\bf R}_{iy}}}\right)^u \left({{\partial }\over{\partial{\bf R}_{iz}}}\right)^v &\cr &\quad\times\exp(2 \pi i {\bf m}\cdot{\bf R}_i).&(} ][E_{\rm self}(\varrho_i)] is given by[\eqalign{E_{\rm self}(\varrho_i) &= \lim_{{\bf R} \to {\bf 0}} \sum_{t_1u_1v_1} \sum_{t_2 u_2 v_2} (-1)^{(t_2+u_2+v_2)} \cr &\quad\times \left({{\partial }\over{\partial{\bf R_x}}}\right)^{t_1+t_2} \left({{\partial}\over{ \partial{\bf R_y}}}\right)^{u_1+u_2} \left({{\partial}\over{ \partial{\bf R_z}}} \right)^{v_1+v_2} \cr &\quad\times \left\lbrace \sum_{(l_1,l_2) \in C\times C}{\bf c}_{i,l_1,t_1u_1v_1}{\bf c}_{i,l_2,t_2u_2v_2} {{{\rm erf}(\beta^{1/2} |{\bf R}|)}\over{|{\bf R}|}} \right. \cr &\quad+ \left. \sum_{(l_1,l_2) \not \in C\times C}{\bf c}_{i,l_1,t_1u_1v_1}{\bf c}_{i,l_2,t_2u_2v_2} {{{\rm erf}(\mu_{12}^{1/2} |{\bf R}|)}\over{|{\bf R}|}} \right\rbrace,} ]where [1 / \mu_{12} = 1 / \theta_{l_1} + 1 / \theta_{l_2}], and the unit-cell dipole components are given by[\eqalign{{\bf D}_x &= \textstyle\sum\limits_{i=1}^N q_i {\bf R}_x + \textstyle\sum\limits_{i=1}^N \textstyle\sum\limits_{l = 1}^L {\bf c}_{i,l,100}, \cr {\bf D}_y &= \textstyle\sum\limits_{i=1}^N q_i {\bf R}_y + \textstyle\sum\limits_{i=1}^N \textstyle\sum\limits_{l = 1}^L {\bf c}_{i,l,010}, \cr {\bf D}_z &= \textstyle\sum\limits_{i=1}^N q_i {\bf R}_z + \textstyle\sum\limits_{i=1}^N \textstyle\sum\limits_{l = 1}^L {\bf c}_{i,l,001}.} ]

In the limit that all Gaussians are compact with [\theta_l \to \infty ], [E_{S,K}(\underline\varrho^{\lbrace\bf N\rbrace})] reduces to the Coulomb energy of ideal multipoles positioned at the expansion points.


Challacombe, M., White, C. & Head-Gordon, M. (1997). Periodic boundary conditions and the fast multipole method. J. Chem. Phys. 107, 10131–10140.
Coppens, P. (1997). X-ray Charge Densities and Chemical Bonding. New York: Oxford University Press.
Gill, P. M. W. & Adamson, R. D. (1996). A family of attenuated Coulomb operators. Chem. Phys. Lett. 261, 105–110.
Helgaker, T., Jorgensen, P. & Olsen, J. (2000). Molecular Electronic-Structure Theory. Chichester: John Wiley and Sons.
Trickey, S. B., Alford, J. A. & Boettger, J. C. (2004). Methods and implementation of robust, high-precision Gaussian basis DFT calculations for periodic systems: the GTOFF code. In Theoretical and Computational Chemistry, edited by J. Leszcynski, Vol. 15, ch. 6. Amsterdam: Elsevier.

to end of page
to top of page