Samedi aprem

algorithms.org

@@ -1,26 +1,26 @@
 #+TITLE: Important algorithms for CIPSI
 #+DATE: 14/12/2020
-#+AUTHOR: Anthony Scemama
+#+AUTHOR: /Anthony Scemama
 
 #+startup: beamer
-#+LaTeX_HEADER: \institute{Laboratoire de Chimie et Physique Quantiques, IRSAMC, UPS/CNRS, Toulouse}
 #+LATEX_CLASS: beamer
 #+LaTeX_CLASS_OPTIONS:[aspectratio=169]
 #+BEAMER_THEME: pterosor
-#+LaTeX_HEADER: \usepackage{minted}
-#+LaTeX_HEADER: \usepackage[utf8]{inputenc}
-#+LaTeX_HEADER: \usepackage[T1]{fontenc}
-#+LaTeX_HEADER: \usepackage{hyperref}
-#+LaTeX_HEADER: \usepackage{mathtools}
-#+LaTeX_HEADER: \usepackage{physics}
-#+LaTeX_HEADER: \definecolor{darkgreen}{rgb}{0.,0.6,0.}
-#+LaTeX_HEADER: \definecolor{darkblue}{rgb}{0.,0.2,0.7}
-#+LaTeX_HEADER: \definecolor{darkred}{rgb}{0.6,0.1,0.1}
-#+LaTeX_HEADER: \definecolor{darkpink}{rgb}{0.7,0.0,0.7}
+#+BEAMER_HEADER_EXTRA: \institute{Laboratoire de Chimie et Physique Quantiques, IRSAMC, UPS/CNRS, Toulouse}
+#+BEAMER_HEADER_EXTRA: \usepackage{minted}
+#+BEAMER_HEADER_EXTRA: \usepackage[utf8]{inputenc}
+#+BEAMER_HEADER_EXTRA: \usepackage[T1]{fontenc}
+#+BEAMER_HEADER_EXTRA: \usepackage{hyperref}
+#+BEAMER_HEADER_EXTRA: \usepackage{mathtools}
+#+BEAMER_HEADER_EXTRA: \usepackage{physics}
 #+EXPORT_EXCLUDE_TAGS: noexport
 #+OPTIONS: H:2 num:t toc:nil \n:nil @:t ::t |:t ^:t -:t f:t *:t <:t
 #+OPTIONS: TeX:t LaTeX:t skip:nil d:nil todo:t pri:nil tags:not-in-toc
 
+#+LATEX: \definecolor{darkgreen}{rgb}{0.,0.6,0.}
+#+LATEX: \definecolor{darkblue}{rgb}{0.,0.2,0.7}
+#+LATEX: \definecolor{darkred}{rgb}{0.6,0.1,0.1}
+#+LATEX: \definecolor{darkpink}{rgb}{0.7,0.0,0.7}
 #+LATEX: \newcommand{\mcenter}[1]{\multicolumn{1}{c}{#1}}
 #+LATEX: \newcommand{\Ndet}{{\textcolor{darkgreen}{N_\text{det}}}}
 #+LATEX: \newcommand{\Evar}{\textcolor{darkgreen}{E_\text{var}}}
@@ -53,10 +53,116 @@
 #+LATEX: \newcommand{\mM}{{\textcolor{blue}{M}} }
 #+LATEX: \newcommand{\epsik}{{\textcolor{darkgreen}{\epsilon_{I_{\textcolor{blue}{k}}}} }}  
 
-* Efficient direct CI
   
-** Sorting
+** Determinant-driven methods
+   - Integral-driven : sequential access to $\mathcal{O}(N^4)$ integrals, indirect access
+     to vectors
+     #+BEGIN_SRC python
+     for (i,j,k,l,integral) in all_integrals:
+         pairs = find_determinant_pairs(i,j,k,l,ijkl)
+         for (d1,d2) in pairs:
+             do_work(d1,d2)
+     #+END_SRC
+   - Determinant-driven : sequential access to vectors, indirect
+     access to integrals
+     #+BEGIN_SRC python
+     for d1 in determinants:
+       for d2 in determinants:
+         i,j,k,l = get_excitation(d1,d2)
+         do_work(d1,d2)
+     #+END_SRC
+   - Integral-driven: $\mathcal{O}(\nmo^4)$
+   - Determinant-driven: $\mathcal{O}(\Ndet^2)$
+   - Efficient CIPSI: How to be efficient within a determinant-driven approach
 
+
+* Data structures
+
+** Slater-Condon's rules
+   \begin{eqnarray*}
+  \langle I | {\cal O}_1 | I \rangle & = & \sum_{i \in D} \langle \varphi_i | {\cal O}_1 | \varphi_i \rangle \\
+  \langle I | {\cal O}_2 | I \rangle & = & \frac{1}{2} \sum_{(i,j) \in D}  
+      \langle \varphi_i \varphi_j | {\cal O}_2 | \varphi_i \varphi_j \rangle 
+      \langle \varphi_i \varphi_j | {\cal O}_2 | \varphi_j \varphi_i \rangle \\
+  \langle I | {\cal O}_1 | \hat{T}_i^j I \rangle & = & \langle \varphi_i | {\cal O}_1 | \varphi_j \rangle \\
+  \langle I | {\cal O}_2 | \hat{T}_i^j I \rangle & = & \sum_{k \in D} 
+      \langle \varphi_i \varphi_k | {\cal O}_2 | \varphi_j \varphi_k \rangle - 
+      \langle \varphi_i \varphi_k | {\cal O}_2 | \varphi_k \varphi_j \rangle \\
+  \langle I | {\cal O}_2 | \hat{T}_{ik}^{jl} I \rangle & = & 
+      \langle \varphi_i \varphi_k | {\cal O}_2 | \varphi_j \varphi_l \rangle -
+      \langle \varphi_i \varphi_k | {\cal O}_2 | \varphi_l \varphi_j \rangle 
+\end{eqnarray*}
+ 
+   Need for functions : $f(I, J) \rightarrow (i,j,k,l,\phi)$
+** Double-determinant representation
+ A Slater determinant can be written as a Waller-Hartree double determinant
+
+ \[ |I \rangle = \hat{I}\,  |\rangle =
+     -1^p \times \hat{I}_\uparrow \, \hat{I}_\downarrow\, |\rangle =
+     -1^p \times \hat{I}_\uparrow \, |\rangle \otimes\ \hat{I}_\downarrow \, |\rangle
+     \]
+ 
+ Storage:
+ - 1 determinant: one integer for $\hat{I}_\uparrow$ and one integer for $\hat{I}_\downarrow$
+ - Set the bit to ~1~ if the orbital is occupied
+ - $>64$ orbitals: $N_{\text{int}}$ integers for $\hat{I}_\uparrow$ and for $\hat{I}_\downarrow$
+   
+** Fast determinant comparisons
+   Bitwise operations (1 CPU cycle):
+   - ~and~, ~or~, ~xor~, ~shl~, ~shr~: logical 
+   - ~shl~, ~shr~: shift left/right
+   - ~lzcnt~, ~tzcnt~ : Number of leading/trailing zero bits
+   - ~popcnt~ : Number of bits set to ~1~
+
+   Example: degree of excitation between $|I\rangle$ and $|J\rangle$:
+   #+BEGIN_SRC fortran
+    integer function degree(det_i, det_j, N_int)
+      integer,   intent(in) :: N_int
+      integer*8, intent(in) :: det_i(N_int,2), det_j(N_int,2)
+      integer               :: two_d, i
+      two_d = 0
+      do i=1,N_int
+         two_d = two_d + popcnt( ieor( det_i(i,1), det_j(i,1) ) ) &
+                       + popcnt( ieor( det_i(i,2), det_j(i,2) ) ) 
+      end do
+      degree = rshift(two_d,1)
+    end function degree
+   #+END_SRC
+
+** Fast determinant comparisons
+   #+ATTR_LATEX: :width 0.7\textwidth
+   [[./slaterRules.pdf]]
+   To get the orbital indices:
+   number of leading/trailing zeros gives the positions of the ~1~'s. 
+** Integrals
+*** Constraints
+   - Integrals require a fast *random access*
+   - 8-fold permtutation symmetry $\langle ij|kl \rangle = \langle kj|il \rangle = \cdots$
+   - Many integrals are zero: need for a sparse data structure
+
+*** Implementation     
+   - Hash table
+   - $f(i,j,k,l) -> K$ gives the same $K$ for all similar permutations
+   - $f(i+1,j,k,l) - f(i,j,k,l)$ is likely to be 1 : locality
+   - Array (cache) for $128^4$ frequently used integrals
+
+** Integrals
+
+   #+Caption: Time to access integrals (in nanoseconds/integral) with different access patterns. The time to generate random numbers (measured as 67~ns/integral) was not counted in the random access results.
+   |-----------+--------+------------|
+   | Access    |  Array | Hash table |
+   |-----------+--------+------------|
+   | $i,j,k,l$ |   9.72 |     125.79 |
+   | $i,j,l,k$ |   9.72 |     120.64 |
+   | $i,k,j,l$ |  10.29 |     144.65 |
+   | $l,k,j,i$ |  88.62 |     125.79 |
+   | $l,k,i,j$ |  88.62 |     120.64 |
+   | Random    | 170.00 |     370.00 |
+   |-----------+--------+------------|
+
+* Efficient direct CI
+
+** Sorting
 *** Popular misconception
 
  *Sorting is /not/ $\mathcal{O}(N \log(N))$* :
@@ -78,8 +184,8 @@
         for (int i=0 ; i<N ; i++) {
             if (A[i] & bitmask) { right[q] = A[i]; q++; }
             else                { left [p] = A[i]; p++; } }
-        radix_sort(left , p, bitmask >> 1 ) ;
-        radix_sort(right, q, bitmask >> 1 ) ;
+        radix_sort(left , p, bitmask >> 1) ;
+        radix_sort(right, q, bitmask >> 1) ;
         for (int i=0 ; i<p ; i++) { A[  i] = left [i] ; }
         for (int i=0 ; i<q ; i++) { A[p+i] = right[i] ; }
     }
@@ -103,24 +209,49 @@
        | 5 | 2 | 4 | 5 | 7 | 8 | 9 | 65 | 4 | 3 | 3 | 5 |  6 | 72 |
        | 2 | 3 | 3 | 4 | 4 | 5 | 5 |  5 | 6 | 7 | 8 | 9 | 65 | 72 |
 
-** Double-determinant
- A Slater determinant can be written as a Waller-Hartree double determinant
+** Double-determinant representation of $\Psi$
 
- \[ |I \rangle = \hat{I}\,  |\rangle =
-     -1^p \times \hat{I}_\uparrow \, \hat{I}_\downarrow\, |\rangle =
-     -1^p \times \hat{I}_\uparrow \, |\rangle \otimes\ \hat{I}_\downarrow \, |\rangle
-     \]
- 
  \[
  \mPsi = \sum_\mi \mci |\mi\rangle = \sum_{\mk=1}^{\Ndetup}
  \sum_{\mm=1}^{\Ndetdn} C_{\textcolor{darkgreen}{km}} \mDup \mDdn
  \]
 
  - If $\mDup$ and $\mDdn$ are represented as $\nmo$ -bit strings, this
-   transformation can be done in $\mathcal{O}(\Ndet \times \nmo)$.
+   transformation can be done in $\mathcal{O}(\Ndet \times \nmo)$ (sorting).
+ - Searching for same-spin excitations: looping over $\mk$ or $\mm$ :
+   $\mathcal{O}(\Ndetup) \sim \mathcal{O}(\sqrt{\Ndet})$
  
+** $\mathcal{H} | \Psi \rangle$
+   For all $\mi = \mDup \mDdn$ in $\mPsi$:
+   - Find indices $p$ of $\uparrow$ singles and $\uparrow \uparrow$ doubles 
+     \[ \langle \mi | \mathcal{H} | \Psi \rangle =
+     \sum_J \langle \mi | \mathcal{H}| J \rangle c_J =
+     \sum_p \langle \mDup \mDdn | \mathcal{H} | D_p^\uparrow \mDdn
+     \rangle C_{pm} \]
+   - Find indices $q$ of $\downarrow$ singles and $\downarrow \downarrow$ doubles 
+     \[ \langle \mi | \mathcal{H} | \Psi \rangle =
+     \sum_J \langle \mi | \mathcal{H}| J \rangle c_J =
+     \sum_q \langle \mDup \mDdn | \mathcal{H} | \mDup
+     D_q^\downarrow \rangle C_{kq} \]
+   - Find indices $pq$ of $\uparrow \downarrow$ doubles:
+     - Find indices $p$ of $\uparrow$ singles 
+     - Find indices $q$ of $\downarrow$ singles \\
+     \[ \langle \mi | \mathcal{H} | \Psi \rangle =
+     \sum_J \langle \mi | \mathcal{H}| J \rangle c_J =
+     \sum_{pq} \langle \mDup \mDdn | \mathcal{H} | D_p^\uparrow
+     D_q^\downarrow \rangle C_{pq} \]
+     
+** Scaling with $N_{\text{det}}$ 
+   #+ATTR_LATEX: :height 0.8\textheight
+   [[./scaling_davidson_ndet.pdf]]
 
-* Stochastic evaluation of the PT2 correction
+** Parallel efficiency
+   #+ATTR_LATEX: :height 0.8\textheight
+   [[./scaling_davidson.pdf]]
+
+
+
+* Stochastic evaluation of the PT2 correction and selection
   
 ** Epstein-Nesbet Second order correction
 
@@ -209,7 +340,7 @@
    \rangle \right)^2}{\Evar - \langle \ma_\mi | \mH | \ma_\mi \rangle}
    \]
 
-** From $\mathcal{O}(N_\text{det}^2)$ to $\mathcal{O}(N_\text{det})$
+** From $\mathcal{O}(N_\text{det}^2)$ to $\mathcal{O}(N_\text{det}^{3/2})$
  \[
  \mPsi = \sum_\mi \mci |\mi\rangle = \sum_{\mk=1}^{\Ndetup} \sum_{\mm=1}^{\Ndetdn} C_{\textcolor{red}{km}} \mDup \mDdn
  \]
@@ -224,17 +355,17 @@
 ** Finer filtering
 
    \[
-   \epsi  = \sum_{\ma \in \mA_\mi} \frac{ \langle \mPsi | \mH | \ma_\mi \rangle
-   \langle \ma_\mi | \mH | \mPsi' \rangle}{\Evar - \langle \ma_\mi | \mH | \ma_\mi
+   \epsi  = \sum_{\ma \in \mA_\mi} \frac{ \langle \mPsi'_I | \mH | \ma_\mi \rangle
+   \langle \ma_\mi | \mH | \mPsi'_I \rangle}{\Evar - \langle \ma_\mi | \mH | \ma_\mi
    \rangle}
    \]
 
    - We know that all the $| \ma_\mi \rangle$ are singles and doubles
      with respect to $| \mi \rangle$
-   - $|\mPsi'\rangle$ is the projection of $\mPsi$ on the subspace of
+   - $|\mPsi'_I\rangle$ is the projection of $|\mPsi \rangle$ on the subspace of
      determinants in $\mD$ which are no more than quadruply excited
      with respect to $| \mi \rangle$
-   - For a subset of excitations $\mi \rightarrow \ma$,
+   - For a subset of excitations $ij \rightarrow ab$,
      $|\mPsi'\rangle$ is filtered further with possible hole/particle constraints
 
 ** Monte Carlo sampling
@@ -456,9 +587,13 @@ def lazy_e(i):
    |         | ---                                 | $\sim$ 29 hr (estimated) |
    |---------+-------------------------------------+--------------------------|
 
+** Scaling with $\Ndet$ 
+   #+ATTR_LATEX: :height 0.8\textheight
+   [[./scaling_sel_det.pdf]]
+
 ** Parallel efficiency
    #+ATTR_LATEX: :height 0.8\textheight
-   [[./speedup_pt2.png]]
+   [[./scaling_sel_node.pdf]]