Q-Chem 4.4 User’s Manual

10.9 Intracules

The many dimensions of electronic wave functions makes them difficult to analyze and interpret. It is often convenient to reduce this large number of dimensions, yielding simpler functions that can more readily provide chemical insight. The most familiar of these is the one-electron density $\rho (\ensuremath{\mathbf{r}})$, which gives the probability of an electron being found at the point $\ensuremath{\mathbf{r}}$. Analogously, the one-electron momentum density $\pi (\ensuremath{\mathbf{p}})$ gives the probability that an electron will have a momentum of $\ensuremath{\mathbf{p}}$. However, the wave function is reduced to the one-electron density much information is lost. In particular, it is often desirable to retain explicit two-electron information. Intracules are two-electron distribution functions and provide information about the relative position and momentum of electrons. A detailed account of the different type of intracules can be found in Ref. Gill:2003. Q-Chem’s intracule package was developed by Aaron Lee and Nick Besley, and can compute the following intracules for or HF wave functions:

10.9.1 Position Intracules

The intracule density, $I(\ensuremath{\mathbf{u}})$, represents the probability for the inter-electronic vector $\ensuremath{\mathbf{u}} = \ensuremath{\mathbf{u}}_1-\ensuremath{\mathbf{u}}_2$:

  \begin{equation} \label{eq1005} I(\ensuremath{\mathbf{u}}) = \int \rho (\ensuremath{\mathbf{r}}_1 \ensuremath{\mathbf{r}}_2 )\; \delta (\ensuremath{\mathbf{r}}_{12} - \ensuremath{\mathbf{u}}) \; d\ensuremath{\mathbf{r}}_1 \, \mbox{d} \ensuremath{\mathbf{r}}_2 \end{equation}   (10.13)

where $\rho (\ensuremath{\mathbf{r}}_1,\ensuremath{\mathbf{r}}_2)$ is the two-electron density. A simpler quantity is the spherically averaged intracule density,

  \begin{equation} \label{eq1006} P(u) = \int I(\ensuremath{\mathbf{u}}) \mbox{d} \Omega _{\ensuremath{\mathbf{u}}} \; , \end{equation}   (10.14)

where $\Omega _{\ensuremath{\mathbf{u}}}$ is the angular part of $\ensuremath{\mathbf{v}}$, measures the probability that two electrons are separated by a scalar distance $u = \vert \ensuremath{\mathbf{u}}\vert $. This intracule is called a position intracule [563]. If the molecular orbitals are expanded within a basis set

  \begin{equation} \label{eq1007} \psi _ a (\ensuremath{\mathbf{r}}) = \sum _\mu c_{\mu a} \phi _\mu (\ensuremath{\mathbf{r}}) \end{equation}   (10.15)

The quantity $P(u)$ can be expressed as

  \begin{equation} \label{eq1008} P(u) = \sum \limits _{\mu \nu \lambda \sigma } \Gamma _{\mu \nu \lambda \sigma } (\mu \nu \lambda \sigma )_\ensuremath{\mathrm{}}{P} \end{equation}   (10.16)

where $ \Gamma _{\mu \nu \lambda \sigma }$ is the two-particle density matrix and $(\mu \nu \lambda \sigma )_{\ensuremath{\mathrm{P}}}$ is the position integral

  \begin{equation} (\mu \nu \lambda \sigma )_\ensuremath{\mathrm{}}{P} = \int \phi _\mu ^\ast (\ensuremath{\mathbf{r}})\; \phi _\nu (\ensuremath{\mathbf{r}})\; \phi _\lambda ^\ast (\ensuremath{\mathbf{r}}+\ensuremath{\mathbf{u}}) \phi _\sigma (\ensuremath{\mathbf{r}}+\ensuremath{\mathbf{u}}) \; d\ensuremath{\mathbf{r}}\; d\Omega \end{equation}   (10.17)

and $\phi _{\mu }(\ensuremath{\mathbf{r}})$, $\phi _{\nu }(\ensuremath{\mathbf{r}})$, $\phi _{\lambda }(\ensuremath{\mathbf{r}})$ and $\phi _{\sigma }(\ensuremath{\mathbf{r}})$ are basis functions. For HF wave functions, the position intracule can be decomposed into a Coulomb component,

  \begin{equation} P_{\ensuremath{\mathrm{J}}}(u) = \frac{1}{2} \sum _{\mu \nu \lambda \sigma } D_{\mu \nu } D_{\lambda \sigma } ({\mu \nu \lambda \sigma })_\ensuremath{\mathrm{}}{P} \end{equation}   (10.18)

and an exchange component,

  \begin{equation} P_\ensuremath{\mathrm{}}{K}(u) = -\frac{1}{2} \sum _{\mu \nu \lambda \sigma } \left[ D_{\mu \lambda }^\alpha D_{\nu \sigma }^\alpha + D_{\mu \lambda }^\beta D_{\nu \sigma }^\beta \right] ({\mu \nu \lambda \sigma })_\ensuremath{\mathrm{}}{P} \end{equation}   (10.19)

where $D_{\mu \nu }$ etc. are density matrix elements. The evaluation of $P(u)$, $P_{\ensuremath{\mathrm{J}}}(u)$ and $P_{\ensuremath{\mathrm{K}}}(u)$ within Q-Chem has been described in detail in Ref. Lee:1999.

Some of the moments of $P(u)$ are physically significant [565], for example

  $\displaystyle \int \limits _0^\infty u^0 P(u) du $ $\displaystyle = $ $\displaystyle \frac{n(n-1)}{2} $   (10.20)
  $\displaystyle \int \limits _0^\infty u^0 P_{\ensuremath{\mathrm{J}}}(u) du $ $\displaystyle = $ $\displaystyle \frac{n^2}{2} $   (10.21)
  $\displaystyle \int \limits _0^\infty u^2 P_{\ensuremath{\mathrm{J}}}(u) du $ $\displaystyle = $ $\displaystyle n Q- \mu ^2 $   (10.22)
  $\displaystyle \int \limits _0^\infty u^0 P_{\ensuremath{\mathrm{K}}}(u) du $ $\displaystyle = $ $\displaystyle -\frac{n}{2} $   (10.23)

where $n$ is the number of electrons and, $\mu $ is the electronic dipole moment and $Q$ is the trace of the electronic quadrupole moment tensor. Q-Chem can compute both moments and derivatives of position intracules.

10.9.2 Momentum Intracules

Analogous quantities can be defined in momentum space; $\bar{I}(\ensuremath{\mathbf{v}})$, for example, represents the probability density for the relative momentum $\ensuremath{\mathbf{v}} = \ensuremath{\mathbf{p}}_1 - \ensuremath{\mathbf{p}}_2$:

  \begin{equation} \bar{I}(\ensuremath{\mathbf{v}}) = \int \pi (\ensuremath{\mathbf{p}}_1,\ensuremath{\mathbf{p}}_2 )\; \delta (\ensuremath{\mathbf{p}}_{12}-\ensuremath{\mathbf{v}}) d\ensuremath{\mathbf{p}}_1 d\ensuremath{\mathbf{p}}_2 \end{equation}   (10.24)

where $\pi (\ensuremath{\mathbf{p}}_1,\ensuremath{\mathbf{p}}_2)$ momentum two-electron density. Similarly, the spherically averaged intracule

  \begin{equation} M(v) = \int \bar{I}(\ensuremath{\mathbf{v}}) \mbox{d}\Omega _{\ensuremath{\mathbf{v}}} \end{equation}   (10.25)

where $\Omega _{\ensuremath{\mathbf{v}}}$ is the angular part of $\ensuremath{\mathbf{v}}$, is a measure of relative momentum $v = \vert \ensuremath{\mathbf{v}}\vert $ and is called the momentum intracule. The quantity $M(v)$ can be written as

  \begin{equation} M(v)=\sum \limits _{\mu \nu \lambda \sigma } \Gamma _{\mu \nu \lambda \sigma } \left( \mu \nu \lambda \sigma \right)_{\ensuremath{\mathrm{M}}} \end{equation}   (10.26)

where $ \Gamma _{\mu \nu \lambda \sigma } $ is the two-particle density matrix and $(\mu \nu \lambda \sigma )_{\ensuremath{\mathrm{M}}}$ is the momentum integral [566]

  \begin{equation} (\mu \nu \lambda \sigma )_{\ensuremath{\mathrm{M}}} = \frac{v^2}{2\pi ^2} \int \phi _\mu ^\ast (\ensuremath{\mathbf{r}}) \phi _\nu (\ensuremath{\mathbf{r}} +\ensuremath{\mathbf{q}}) \phi _\lambda ^\ast (\ensuremath{\mathbf{u}}+\ensuremath{\mathbf{q}}) \phi _\sigma (\ensuremath{\mathbf{u}}) j_0(q v)\; d\ensuremath{\mathbf{r}}\; d\ensuremath{\mathbf{q}}\; d\ensuremath{\mathbf{u}} \end{equation}   (10.27)

The momentum integrals only possess four-fold permutational symmetry, i.e.,

  $\displaystyle (\mu \nu \lambda \sigma )_{\ensuremath{\mathrm{M}}} = (\nu \mu \lambda \sigma )_{\ensuremath{\mathrm{M}}} = (\sigma \lambda \nu \mu )_{\ensuremath{\mathrm{M}}} = (\lambda \sigma \mu \nu )_{\ensuremath{\mathrm{M}}} $   (10.28)
  $\displaystyle (\nu \mu \lambda \sigma )_{\ensuremath{\mathrm{M}}} = (\mu \nu \sigma \lambda )_{\ensuremath{\mathrm{M}}} = (\lambda \sigma \nu \mu )_{\ensuremath{\mathrm{M}}} = (\sigma \lambda \mu \nu )_{\ensuremath{\mathrm{M}}} $   (10.29)

and therefore generation of $M(v)$ is roughly twice as expensive as $P(u)$. Momentum intracules can also be decomposed into Coulomb $M_{\ensuremath{\mathrm{J}}}(v)$ and exchange $M_{\ensuremath{\mathrm{K}}}(v)$ components:

  \begin{equation} M_ J(v) = \frac{1}{2}\sum \limits _{\mu \nu \lambda \sigma } D_{\mu \nu } D_{\lambda \sigma }(\mu \nu \lambda \sigma )_{\ensuremath{\mathrm{M}}} \end{equation}   (10.30)
  \begin{equation} M_{\ensuremath{\mathrm{K}}}(v) = -\frac{1}{2} \sum \limits _{\mu \nu \lambda \sigma } \left[ D_{\mu \lambda }^\alpha D_{\nu \sigma }^\alpha + D_{\mu \lambda }^\beta D_{\nu \sigma }^\beta \right] (\mu \nu \lambda \sigma )_{\ensuremath{\mathrm{M}}} \end{equation}   (10.31)

Again, the even-order moments are physically significant [566]:

  \begin{equation} \int \limits _0^\infty v^0 M(v) dv = \frac{n(n-1)}{2} \end{equation}   (10.32)
  \begin{equation} \int \limits _0^\infty u^0 M_{\ensuremath{\mathrm{J}}}(v) dv = \frac{n^2}{2} \end{equation}   (10.33)
  \begin{equation} \int \limits _0^\infty v^2 P_{\ensuremath{\mathrm{J}}}(v) dv = 2 n E_{\ensuremath{\mathrm{T}}} \end{equation}   (10.34)
  \begin{equation} \int \limits _0^\infty v^0 M_{\ensuremath{\mathrm{K}}}(v) dv = -\frac{n}{2} \end{equation}   (10.35)

where $n$ is the number of electrons and $E_{\ensuremath{\mathrm{T}}}$ is the total electronic kinetic energy. Currently, Q-Chem can compute $M(v)$, $M_{\ensuremath{\mathrm{J}}}(v)$ and $M_{\ensuremath{\mathrm{K}}}(v)$ using $s$ and $p$ basis functions only. Moments are generated using quadrature and consequently for accurate results $M(v)$ must be computed over a large and closely spaced $v$ range.

10.9.3 Wigner Intracules

The intracules $P(u)$ and $M(v)$ provide a representation of an electron distribution in either position or momentum space but neither alone can provide a complete description. For a combined position and momentum description an intracule in phase space is required. Defining such an intracule is more difficult since there is no phase space second-order reduced density. However, the second-order Wigner distribution [567],

  \begin{equation} W_2(\ensuremath{\mathbf{r}}_1,\ensuremath{\mathbf{p}}_1,\ensuremath{\mathbf{r}}_2,\ensuremath{\mathbf{p}}_2) = \frac{1}{\pi ^6} \int \rho _2(\ensuremath{\mathbf{r}}_1+\ensuremath{\mathbf{q}}_1, \ensuremath{\mathbf{r}}_1-\ensuremath{\mathbf{q}}_1, \ensuremath{\mathbf{r}}_2+\ensuremath{\mathbf{q}}_2, \ensuremath{\mathbf{r}}_2-\ensuremath{\mathbf{q}}_2) e^{-2i(\ensuremath{\mathbf{p}}_1 \cdot \ensuremath{\mathbf{q}}_1 + \ensuremath{\mathbf{p}}_2 \cdot \ensuremath{\mathbf{q}}_2)} d\ensuremath{\mathbf{q}}_1 d\ensuremath{\mathbf{q}}_2 \end{equation}   (10.36)

can be interpreted as the probability of finding an electron at $\ensuremath{\mathbf{r}}_1$ with momentum $\ensuremath{\mathbf{p}}_1$ and another electron at $\ensuremath{\mathbf{r}}_2$ with momentum $\ensuremath{\mathbf{p}}_2$. [The quantity $W_2(\ensuremath{\mathbf{r}}_1,\ensuremath{\mathbf{r}}_2,\ensuremath{\mathbf{p}}_1,\ensuremath{\mathbf{p}}_2$ is often referred to as “quasi-probability distribution” since it is not positive everywhere.]

The Wigner distribution can be used in an analogous way to the second order reduced densities to define a combined position and momentum intracule. This intracule is called a Wigner intracule, and is formally defined as

  \begin{equation} W(u,v) = \int W_2(\ensuremath{\mathbf{r}}_1, \ensuremath{\mathbf{p}}_1, \ensuremath{\mathbf{r}}_2, \ensuremath{\mathbf{p}}_2) \delta (\ensuremath{\mathbf{r}}_{12}-\ensuremath{\mathbf{u}}) \delta (\ensuremath{\mathbf{p}}_{12}-\ensuremath{\mathbf{v}}) d\ensuremath{\mathbf{r}}_1\, d\ensuremath{\mathbf{r}}_2\, d\ensuremath{\mathbf{p}}_1\, d\ensuremath{\mathbf{p}}_2\, d\Omega _{\ensuremath{\mathbf{u}}}\, d\Omega _{\ensuremath{\mathbf{v}}} \end{equation}   (10.37)

If the orbitals are expanded in a basis set, then $W(u,v)$ can be written as

  \begin{equation} W(u,v) = \sum \limits _{\mu \nu \lambda \sigma } \Gamma _{\mu \nu \lambda \sigma } \left( {\mu \nu \lambda \sigma } \right)_{\ensuremath{\mathrm{W}}} \end{equation}   (10.38)

where ($\mu \nu \lambda \sigma )_{\ensuremath{\mathrm{W}}}$ is the Wigner integral

  \begin{equation} \label{eq1009} (\mu \nu \lambda \sigma )_{\ensuremath{\mathrm{W}}} = \frac{v^2}{2\pi ^2} \int \int \phi _{\mu }^\ast (\ensuremath{\mathbf{r}}) \phi _{\nu }(\ensuremath{\mathbf{r}}+\ensuremath{\mathbf{q}}) \phi _\lambda ^\ast (\ensuremath{\mathbf{r}}+\ensuremath{\mathbf{q}}+\ensuremath{\mathbf{u}}) \phi _{\sigma }(\ensuremath{\mathbf{r}}+\ensuremath{\mathbf{u}}) j_0(q\, v)\; d\ensuremath{\mathbf{r}}\; d{q}\; d\Omega _{\ensuremath{\mathbf{u}}} \end{equation}   (10.39)

Wigner integrals are similar to momentum integrals and only have four-fold permutational symmetry. Evaluating Wigner integrals is considerably more difficult that their position or momentum counterparts. The fundamental $\left[ssss\right]_{\ensuremath{\mathrm{w}}}$ integral,

  $\displaystyle \left[ssss\right]_{\ensuremath{\mathrm{W}}} $ $\displaystyle = $ $\displaystyle \frac{u^2v^2}{2\pi ^2}\; \int \int \exp \left[-\alpha |\ensuremath{\mathbf{r}}\! -\! \ensuremath{\mathbf{A}}|^2 -\! \beta |\ensuremath{\mathbf{r}}\! +\! \ensuremath{\mathbf{q}}\! -\! \ensuremath{\mathbf{B}}|^2 -\! \gamma |\ensuremath{\mathbf{r}}\! +\! \ensuremath{\mathbf{q}}\! +\! \ensuremath{\mathbf{u}}\! -\! \ensuremath{\mathbf{C}}|^2 -\! \delta |\ensuremath{\mathbf{r}}\! +\! \ensuremath{\mathbf{u}}\! -\! \ensuremath{\mathbf{D}}|^2 \right] \times \nonumber $    
  $\displaystyle $ $\displaystyle $ $\displaystyle j_0(q v)\; d\ensuremath{\mathbf{r}}\; d\ensuremath{\mathbf{q}}\; d\Omega _{\ensuremath{\mathbf{u}}} $   (10.40)

can be expressed as

  \begin{equation} \left[ssss\right]_{\ensuremath{\mathrm{W}}} = \frac{\pi u^2 v^2\; e^{-(R+\lambda ^2 u^2 +\mu ^2 v^2)}}{ 2(\alpha +\delta )^{3/2}(\beta +\gamma )^{3/2}}\int {e^{-\ensuremath{\mathbf{P}}\cdot \ensuremath{\mathbf{u}}}} j_0 \left(|\ensuremath{\mathbf{Q}}+\eta \ensuremath{\mathbf{u}}|v \right)\; d\Omega _ u \end{equation}   (10.41)

or alternatively

  \begin{equation} \left[ssss\right]_{\ensuremath{\mathrm{W}}} = \frac{2\pi ^2 u^2 v^2 e^{-(R+\lambda ^2 u^2+\mu ^2 v^2)}}{(\alpha +\delta )^{3/2}(\beta +\gamma )^{3/2}} \sum \limits _{n=0}^\infty (2n+1) i_ n(P\, u) j_ n(\eta u v) j_ n(Q v) P_ n \left( {\frac{\ensuremath{\mathbf{P}}\cdot \ensuremath{\mathbf{Q}}}{P\; Q}} \right) \end{equation}   (10.42)

Two approaches for evaluating $(\mu \nu \lambda \sigma )_{\ensuremath{\mathrm{W}}}$ have been implemented in Q-Chem, full details can be found in Ref. Wigner:1932. The first approach uses the first form of $\left[ssss\right]_{\ensuremath{\mathrm{W}}}$ and used Lebedev quadrature to perform the remaining integrations over $\Omega _{\ensuremath{\mathbf{u}}}$. For high accuracy large Lebedev grids [148, 146, 149] should be used, grids of up to 5294 points are available in Q-Chem. Alternatively, the second form can be adopted and the integrals evaluated by summation of a series. Currently, both methods have been implemented within Q-Chem for $s$ and $p$ basis functions only.

When computing intracules it is most efficient to locate the loop over $u$ and/or $v$ points within the loop over shell-quartets [569]. However, for $W(u,v)$ this requires a large amount of memory to store all the integrals arising from each $(u,v)$ point. Consequently, an additional scheme, in which the $u$ and $v$ points loop is outside the shell-quartet loop, is available. This scheme is less efficient, but substantially reduces the memory requirements.

10.9.4 Intracule Job Control

The following $rem variables can be used to control the calculation of intracules.

INTRACULE

Controls whether intracule properties are calculated (see also the $intracule section).


TYPE:

LOGICAL


DEFAULT:

FALSE


OPTIONS:

FALSE

No intracule properties.

TRUE

Evaluate intracule properties.


RECOMMENDATION:

None


WIG_MEM

Reduce memory required in the evaluation of $W(u,v)$.


TYPE:

LOGICAL


DEFAULT:

FALSE


OPTIONS:

FALSE

Do not use low memory option.

TRUE

Use low memory option.


RECOMMENDATION:

The low memory option is slower, so use the default unless memory is limited.


WIG_LEB

Use Lebedev quadrature to evaluate Wigner integrals.


TYPE:

LOGICAL


DEFAULT:

FALSE


OPTIONS:

FALSE

Evaluate Wigner integrals through series summation.

TRUE

Use quadrature for Wigner integrals.


RECOMMENDATION:

None


WIG_GRID

Specify angular Lebedev grid for Wigner intracule calculations.


TYPE:

INTEGER


DEFAULT:

194


OPTIONS:

Lebedev grids up to 5810 points.


RECOMMENDATION:

Larger grids if high accuracy required.


N_WIG_SERIES

Sets summation limit for Wigner integrals.


TYPE:

INTEGER


DEFAULT:

10


OPTIONS:

$n<100$


RECOMMENDATION:

Increase $n$ for greater accuracy.


N_I_SERIES

Sets summation limit for series expansion evaluation of $i_ n(x)$.


TYPE:

INTEGER


DEFAULT:

40


OPTIONS:

$n>0$


RECOMMENDATION:

Lower values speed up the calculation, but may affect accuracy.


N_J_SERIES

Sets summation limit for series expansion evaluation of $j_ n(x)$.


TYPE:

INTEGER


DEFAULT:

40


OPTIONS:

$n>0$


RECOMMENDATION:

Lower values speed up the calculation, but may affect accuracy.


10.9.5 Format for the $intracule Section

int_type

0

Compute $P(u)$ only

 

1

Compute $M(v)$ only

 

2

Compute $W(u,v)$ only

 

3

Compute $P(u)$, $M(v)$ and $W(u,v)$

 

4

Compute $P(u)$ and $M(v)$

 

5

Compute $P(u)$ and $W(u,v)$

 

6

Compute $M(v)$ and $W(u,v)$

u_points

 

Number of points, start, end.

v_points

 

Number of points, start, end.

moments

0–4

Order of moments to be computed ($P(u)$ only).

derivs

0–4

order of derivatives to be computed ($P(u)$ only).

accuracy

$n$

($10^{-n})$ specify accuracy of intracule interpolation table ($P(u)$ only).

Example 10.226  Compute HF/STO-3G $P(u)$, $M(v)$ and $W(u,v)$ for Ne, using Lebedev quadrature with 974 point grid.

$molecule
   0  1
   Ne
$end

$rem
   METHOD      hf
   BASIS       sto-3g
   INTRACULE   true
   WIG_LEB     true
   WIG_GRID    974
$end

$intracule
   int_type   3
   u_points  10   0.0  10.0
   v_points   8   0.0   8.0
   moments    4
   derivs     4
   accuracy   8
$end

Example 10.227  Compute HF/6-31G $W(u,v)$ intracules for H$_{2}$O using series summation up to $n$=25 and 30 terms in the series evaluations of $j_ n(x)$ and $i_ n(x)$.

$molecule
   0  1
   H1
   O   H1  r
   H2  O   r  H1  theta

   r = 1.1
   theta = 106
$end

$rem
   METHOD         hf
   BASIS          6-31G
   INTRACULE      true
   WIG_MEM        true
   N_WIG_SERIES   25
   N_I_SERIES     40
   N_J_SERIES     50
$end

$intracule
   int_type    2
   u_points   30   0.0   15.0
   v_points   20   0.0   10.0
$end