fast atom-pair decomposition of the SELF (Pauli) descriptor, 2/2 partition More...

#include <vector>

Include dependency graph for SELFATOMIC.h:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Detailed Description

fast atom-pair decomposition of the SELF (Pauli) descriptor, 2/2 partition

The quantity partitioned here is

 c(r) = 8.rho.G - grad(rho)^2
      = 4 x Sum_(a,b,c,d) [Dab.Dcd - Dac.Dbd] x (a'.b') x c x d

The 2/2 partition keeps only the quadruplets having two AOs on FRAG1 and two AOs on FRAG2. The historical implementation (NCISolver.cpp, ~line 3604) walks those quadruplets explicitly over the PRIMITIVES: six quadruple loops, i.e. O(N^4) work per grid point.

Here the same quantity is obtained in O(N^2) by first contracting the density matrix into three ATOM-PAIR resolved quantities (I and J are ATOM indices, x is a cartesian component, a runs over the AOs borne by atom I).

Note that the AO products a.b, a'.b' and a'.b are still needed – they are the bricks of R, S and T below. What goes away is their TABULATION into the OV, DOV and COV matrices: each product is now formed inside the contraction loop and consumed on the spot. Those tables were sized on the full fragNbPrim^2, rebuilt at every grid point, and the primRmin pruning did not apply to them.

 R_IJ    = Sum_(a in I, b in J) Dab .  a  .  b        (symmetric)
 S_x,IJ  = Sum_(a in I, b in J) Dab . da  .  b        (NOT symmetric)
 T_x,IJ  = Sum_(a in I, b in J) Dab . da  . db        (symmetric)

so that

 c(r) = 4 Sum_x Sum_IJKL [ T_x,IJ . R_KL - S_x,IK . S_x,JL ]

Beware of the index pairing: the direct term pairs (I,J) and (K,L) whereas the exchange term pairs (I,K) and (J,L). This follows from Dab.Dcd versus Dac.Dbd and must be respected when binning into atom pairs.

Each 2/2 quadruplet holds two atoms per fragment, hence two INTRA-FRAGMENT atom pairs; its value is split half and half between them. This is exactly what the historical code does, so the pauli[][] array is reproduced as is.

See the design note for the full derivation.