fix: actually add the posts.

src/posts/hf.org

@@ -0,0 +1,297 @@
+#+TITLE: Helonium's Hartree-Fock program
+#+INCLUDE: "../html-head.org"
+
+** Exordium
+
+We will implement the Hartree-Fock[fn:1] program from the classic [[https://store.doverpublications.com/products/9780486691862][Szabo-Ostlund]] text,
+a staple in quantum chemistry. If you have any experience in the field, chances are you know it well.
+If you don't, the BQN implementation will take you only a few cognitive units: the coarse mathematical
+description involves basis sets, the calculation of electronic integrals, and the self-consistent
+optimization of the Fock matrix. Using this program, we will compute the energy of the HeH\(^+\) molecule[fn:2].
+
+First, we import the required BQN system values and utility functions:
+
+#+begin_src bqn :results none :tangle ../bqn/hf.bqn
+  Sin‿Cos‿ATan‿Erf ← •math
+  Setplot‿Plot ← •Import "../bqn-utils/plots.bqn"
+#+end_src
+
+Some auxiliary functions are needed for the computation of the Boys function, the matrix and vector product,
+the spectrum of a matrix with shape 2‿2, and Löwdin's canonical orthogonalization. Additionally, we define
+a 1-modifier to compute the fixed point of a function:
+
+#+begin_src bqn :results none :tangle ../bqn/hf.bqn
+  E ← 1e¯6⊸<◶(1-÷⟜3)‿((π÷4)⊸÷×⟜Erf○√⊢)
+  P ← +˝∘×⎉1‿∞
+  D ← (⊢⋈·P´<∘⊢∾⋈)⟜([Cos‿Sin⋄Sin‿(-Cos)]{𝕎𝕩}¨·<2÷˜·ATan 2×0‿1⊸⊑÷·-´0‿0⊸⍉)
+  L ← (-⌾(1‿1⊸⊑)·÷∘√≢⥊2˙)⊸P(0¨⌾(0‿0⍉⌽˘)·÷∘√1+1‿¯1×⌽˘)
+  _fp ← {𝕊∘⊢⍟≢⟜𝔽𝕩}
+#+end_src
+
+Then, we define a function returning a namespace with the physical constants of the system[fn:3],
+as is customary in /ab-initio/ methods. We have the flexibility of providing different molecular
+geometries:
+
+#+begin_src bqn :results none :tangle ../bqn/hf.bqn
+  System ← {𝕊𝕩:
+    e1‿e2 ⇐ 2.0925‿1.24
+    z1‿z2 ⇐ 2‿1
+    r ⇐ 𝕩
+  }
+#+end_src
+
+** Basis set
+
+Basis sets are used to transform the PDEs into linear algebra problems. Physical intuition suggests that
+Slater type orbitals[fn:4] are a good choice for our Hamiltonian. However, the computation of the integrals
+is a lot easier if we approximate them using Gaussian functions[fn:5]. Specifically, the STO-3G approach defines
+the approximating function as follows:
+
+#+begin_src bqn :results none :tangle ../bqn/hf.bqn
+  STO ← {
+    e ← 0.109818‿0.405771‿2.22766 ×ט 𝕩
+    c ← 0.444635‿0.535328‿0.154329
+    e⋈c×(2×e÷π)⋆3÷4
+  }
+#+end_src
+
+** Electronic integrals
+
+Constructing the integrals' tensor is complicated[fn:6] and is the main reason for the poor scaling
+of electronic structure methods. The \(1s\) orbitals are the simplest case, and here two types of integrals
+are analytical (S, T) while the rest already lacks a closed-form solution (V, ERI):
+
+#+begin_src bqn :results none :tangle ../bqn/hf.bqn 
+  S ← {a‿b𝕊𝕩: (1.5⋆˜π÷a+b) × ⋆-𝕩×a(×÷+)b}
+  T ← {a‿b𝕊𝕩: f ← a(×÷+)b ⋄ ×´⟨1.5⋆˜π÷a+b, (3×f)-2×𝕩×טf, ⋆-𝕩×f⟩}
+  V ← {a‿b‿z𝕊r‿s: ×´⟨-2×z×π÷a+b, E s×a+b, ⋆-r×a(×÷+)b⟩}
+  ERI ← {a‿b‿c‿d𝕊r1‿r2‿r3‿r4:
+    r5 ← -´⟨a‿b ⋄ c‿d⟩ (+´∘×÷+´∘⊣)¨ ⟨r1‿r2 ⋄ r3‿r4⟩
+    f1‿f2‿f3 ← a‿b ({(√∘+××)⋈(×÷+)}○(+´)∾<∘⋈○((×÷+)´)) c‿d
+    ×´⟨f1÷˜2×π⋆5÷2, E f2×טr5, ⋆-+´f3×ט-´¨⟨r1‿r2, r3‿r4⟩⟩
+  }
+#+end_src
+
+#+begin_export html
+<br/>
+<details>
+<summary>Derivation strategy</summary>
+#+end_export
+
+We need to compute the overlap (S), kinetic energy (T), nuclear attraction (V), and four-center (ERI) integrals.
+Crucially, the product of two Gaussians at different centers is proportional to a Gaussian at a scaled center.
+This property, combined with the Laplacian of a Gaussian, readily yields S and T. The remaining
+two sets are more complex: we combine the Gaussians as before, then transform to reciprocal space where
+the delta distribution arises and simplifies the problem to this integration by reduction:
+
+\begin{equation*}
+  I(x) = \int_0^{\infty}{{{e^ {- a\,k^2 }\,\sin \left(k\,x\right)}\over{k}}\;dk} \sim \text{Erf}(x)
+\end{equation*}
+
+#+begin_export html
+</details>
+#+end_export
+
+The following namespace exports the corresponding integral arrays. Extending the code to an arbitrary number
+of atoms implies mapping over an array of coordinates, as opposed to fusing them in the implementation.
+
+#+begin_src bqn :tangle ../bqn/hf.bqn :results none
+  I ← {𝕊e1‿e2‿z1‿z2‿r:
+    bs‿na‿nb ← (<∾·≢⊏∘>)⍉>STO¨ e1‿e2
+    M ← {∾‿×({2: {nb(⋈˜/∘⋈˜)⊸⊔𝕎⌜˜𝕩}; 4: {𝕎⌜⍟3˜𝕩}}𝕩)¨○⊢<∘∾˘bs}
+
+    sm‿hcore ⇐ {e𝕊c:
+      mst ← ⌽⊸≍∾⟜0טr
+      r1‿r2 ← <˘⍉⁼> (r⊸-⊸⋈˜×⟜r÷+)⌜´ ⊏bs
+      mv ← ט∘{[0‿2,3‿1]⊏({0‿𝕨¨𝕩}⟜𝕩¨𝕨)∾⋈⟜⍉r⋈¨𝕩}´¨⟨0‿r, r1⟩‿⟨r‿0, r2⟩
+      (⊑⋈·+´1⊸↓)+´∘⥊¨¨ c<⊸× ({e𝕏¨¨mst}¨S‿T) ∾ z1‿z2{e∾⟜𝕨⊸V¨¨𝕩}¨mv
+    }´ M 2
+
+    erim ⇐ {e𝕊c:
+      meri ← (c⊸×⊣ERI¨⊢/˜·<¨≢∘⊣÷≢∘⊢)⟜{0‿r⊏˜⚇1↕na¨↕=𝕩} e
+      =⊸{+˝∘⥊⎉𝕨 (2×↕𝕨)⍉⁼(na‿nb⥊˜na×𝕨)⥊𝕩} meri
+    }´ M 4
+  }
+#+end_src
+
+#+begin_export html
+<br/>
+<details>
+<summary>Performance</summary>
+#+end_export
+
+The computation of the ERIs is expected to be the primary bottleneck, as there are =N⋆4= of them—in our case, 16.
+The required tensors have a shape of =6¨↕4=. As shown in the profile below, using an array-based strategy
+for the ERIs significantly improved their computational efficiency compared to the two-center integrals.
+For the latter, I increased the depth by grouping the tables (block matrices). The resulting code was significantly
+slower than replicating the elements to match each axis' length, like I do for the ERIs.
+
+#+begin_src bqn :exports both :tangle no :results raw :wrap example
+  )profile {𝕊: I∘System 1.4632}¨↕1e4
+#+end_src
+
+#+RESULTS:
+#+begin_example
+Got 38006 samples
+(self-hosted runtime1): 1067 samples
+(REPL): 36939 samples:
+    72│I ← {𝕊e1‿e2‿z1‿z2‿r:
+    68│  bs‿na‿nb ← (<∾·≢⊏∘>)⍉>STO¨ e1‿e2
+  2053│  M ← {∾‿×({2: {nb(⋈˜/∘⋈˜)⊸⊔𝕎⌜˜𝕩}; 4: {𝕎⌜⍟3˜𝕩}}𝕩)¨○⊢<∘∾˘bs}
+      │
+   245│  sm‿hcore ⇐ {e𝕊c:
+    75│    mst ← ⌽⊸≍∾⟜0טr
+  4181│    r1‿r2 ← <˘⍉⁼> (r⊸-⊸⋈˜×⟜r÷+)⌜´ ⊏bs
+ 16277│    mv ← ט∘{[0‿2,3‿1]⊏({0‿𝕨¨𝕩}⟜𝕩¨𝕨)∾⋈⟜⍉r⋈¨𝕩}´¨⟨0‿r, r1⟩‿⟨r‿0, r2⟩
+  8830│    (⊑⋈·+´1⊸↓)+´∘⥊¨¨ c<⊸× ({e𝕏¨¨mst}¨S‿T) ∾ z1‿z2{e∾⟜𝕨⊸V¨¨𝕩}¨mv
+  3708│  }´ M 2
+      │
+     8│  erim ⇐ {e𝕊c:
+  1100│    meri ← (c⊸×⊣ERI¨⊢/˜·<¨≢∘⊣÷≢∘⊢)⟜{0‿r⊏˜⚇1↕na¨↕=𝕩} e
+   318│    =⊸{+˝∘⥊⎉𝕨 (2×↕𝕨)⍉⁼(na‿nb⥊˜na×𝕨)⥊𝕩} meri
+     4│  }´ M 4
+      │}
+#+end_example
+
+Morals: Never underestimate the power of vectorization and reshaping operations are often computationally trivial.
+
+#+begin_export html
+</details>
+#+end_export
+
+** Fock matrix
+
+The following function constructs the Fock matrix, our approximation to the true Hamiltonian of the system:
+
+#+begin_src bqn :results none :tangle ../bqn/hf.bqn
+  F ← {𝕩.hcore + 𝕨 (+˝∘⥊⎉2⊣×⎉2⊢-2÷˜0‿3‿2⊸⍉⁼) 𝕩.erim}
+#+end_src
+
+#+begin_export html
+<br/>
+<details>
+<summary>Physical context</summary>
+#+end_export
+
+The Fock operator is an effective one-electron operator that arises after constrained
+minimization of the energy functional. The form of the functional is a consequence of
+the use of [[https://en.wikipedia.org/wiki/Slater_determinant][Slater determinants]] as wave functions.
+
+\begin{equation*}
+  \tilde{\mathcal{F}} \left[ \{\psi_i\} \right] = \sum_i h_i +
+  \frac{1}{2} \sum_{i,j} (J_{ij} - K_{ij}) - \sum_{i,j} \lambda_{ij}
+  \left( \langle \psi_i | \psi_j \rangle - \delta_{ij} \right)
+\end{equation*}
+
+where \(h_i\) is the core Hamiltonian matrix, \(J_{ij}, K_{ij}\) are the Coulomb and
+exchange components of the ERI matrix, and \(\lambda_{ij}\) are Lagrange multipliers.
+To fully understand the derivation, consider the variational derivative of this
+functional with respect to the complex conjugate of the one-particle wave function \(\psi_i^*\):
+
+\begin{align*}
+  \lim_{\epsilon \to 0} \frac{\tilde{\mathcal{F}} \left[ \psi_k^* + \epsilon \delta
+   \psi_k^* \right] - \tilde{\mathcal{F}} \left[ \psi_k^* \right]}{\epsilon} 
+  &= \langle \delta \psi_k | \hat{h} | \psi_k \rangle + \sum_j \left( \langle \delta
+   \psi_k \psi_j | \frac{1}{r} | \psi_k \psi_j \rangle - \langle \delta
+   \psi_k \psi_j | \frac{1}{r} | \psi_j \psi_k \rangle \right)
+   - \sum_j \lambda_{kj} \langle \delta \psi_k | \psi_j \rangle \\
+  &= \int \left[ \hat{h} \psi_k(x) + \sum_j
+     \left( \psi_k(x) \int \frac{|\psi_j(x')|^2}{|r - r'|} dx'
+   - \psi_j(x) \int \frac{\psi_j^*(x') \psi_k(x')}{|r - r'|} dx' \right) \right.
+   \left. - \sum_j \lambda_{kj} \psi_j(x) \right] \delta \psi_k^*(x) \, dx.
+\end{align*}
+
+As discussed earlier, basis sets are used to discretize the Hartree-Fock problem.
+This process results in the  [[https://en.wikipedia.org/wiki/Roothaan_equations][Roothaan equations]], which are implemented in the code below.
+
+#+begin_export html
+</details>
+#+end_export
+
+** Self-consistent field
+
+The final stage of the computation involves solving the pseudo-eigenvalue problem using a fixed-point iteration.
+This process is commonly known as the self-consistent field method, a term coined by D. R. Hartree.
+
+#+begin_src bqn :tangle ../bqn/hf.bqn :results none
+  SCF ← {
+    ints ← I𝕩
+    pm ← {2××⌜˜⊏⍉ xm(⊢P⟜⊑·D·P´⍉∘⊢<⊸∾⋈)˜𝕩F ints}_fp 0¨xm ← L ints.sm
+    2÷˜+´∘⥊ pm (⊣⍉⊸×{𝕩.hcore + 𝕨F𝕩}) ints
+  }
+#+end_src
+
+If you are receptive and humble, mathematics will lead you by the hand[fn:7]:
+
+#+begin_src bqn :tangle ../bqn/hf.bqn :exports both
+  SCF∘System 1.4632
+#+end_src
+
+#+RESULTS:
+: ¯4.22752930421725
+
+Compare the electronic energy with the one computed using the original [[./supp/hf_so/hf_so.html][Fortran]] program.
+In terms of performance, the CBQN implementation runs in 5 ms, which is of the same order
+of magnitude as the original program. Notably, the BQN version consists of just 45 lines of code,
+compared to 541 lines in the Fortran version.
+
+#+begin_export html
+<details>
+<summary>Potential Energy Surface</summary>
+#+end_export
+
+Here we compute the system's [[https://en.wikipedia.org/wiki/Potential_energy_surface][PES]]. To do this, we need to add to the electronic energy above
+the nuclear repulsion energy. We also catch the error of non-converged calculations, instead
+of fiddling with convergence thresholds and different starting points:
+
+#+begin_src bqn :results none :tangle ../bqn/hf.bqn
+  PES ← 2⊸÷+SCF⎊∞∘System
+#+end_src
+
+#+NAME: attr_wrap
+#+BEGIN_SRC sh :var data="" :results output :exports none :tangle no
+  echo "<br/>"
+  echo '<div style="display: flex; justify-content: center; width: 100%;">'
+  echo '<div style="width: 40%;">'
+  echo "$data"
+  echo "</div>"
+  echo "</div>"
+#+END_SRC
+
+Then we leverage my modified version of the =•Plot= namespace:
+
+#+begin_src bqn :results html :exports both :tangle ../bqn/hf.bqn :post attr_wrap(data=*this*)
+  )r Setplot "line" ⋄ •Out¨ Plot´ (⊢/¨˜·<∞>⊢´)(⊢⋈PES¨) ↕∘⌈⌾((0.5+1e¯2×⊢)⁼)3
+#+end_src
+
+#+RESULTS:
+#+begin_export html
+<br/>
+<div style="display: flex; justify-content: center; width: 100%;">
+<div style="width: 40%;">
+<svg viewBox='-10 -10 404 201.112'>
+  <g font-family='BQN,monospace' font-size='18px'>
+    <rect class='code' style='fill:none;stroke:black' stroke-width='1' rx='5' x='-5' y='-5' width='394' height='191.112'/>
+    <path class='code' style='fill:none;stroke:#267CB9' stroke-width='3' d='M0 0L1.542 8.542L3.084 16.664L6.169 31.743L7.711 38.744L10.795 51.771L18.506 79.417L20.048 84.216L21.59 88.799L23.133 93.178L24.675 97.362L26.217 101.361L27.759 105.183L29.301 108.836L30.843 112.33L35.47 121.919L37.012 124.841L41.639 132.864L43.181 135.309L46.265 139.883L49.349 144.065L50.892 146.02L57.06 153.015L58.602 154.574L63.229 158.843L64.771 160.14L66.313 161.377L67.855 162.557L70.94 164.756L72.482 165.779L78.651 169.403L80.193 170.202L81.735 170.961L83.277 171.681L87.904 173.628L89.446 174.209L90.988 174.759L94.072 175.769L95.614 176.231L97.157 176.667L98.699 177.076L100.241 177.46L101.783 177.82L106.41 178.763L107.952 179.035L109.494 179.287L114.12 179.929L117.205 180.269L118.747 180.415L121.831 180.661L124.916 180.849L126.458 180.923L128 180.984L129.542 181.032L131.084 181.069L132.627 181.094L134.169 181.108L135.711 181.112L137.253 181.105L141.88 181.026L143.422 180.982L146.506 180.868L148.048 180.799L151.133 180.638L152.675 180.547L155.759 180.345L158.843 180.119L160.386 179.997L163.47 179.736L168.096 179.308L169.639 179.156L172.723 178.84L174.265 178.676L177.349 178.337L178.892 178.162L180.434 177.984L181.976 177.804L183.518 177.62L185.06 177.433L186.602 177.244L188.145 177.053L191.229 176.663L194.313 176.265L197.398 175.859L198.94 175.654L206.651 174.608L208.193 174.395L209.735 174.182L211.277 173.967L212.819 173.752L214.361 173.537L215.904 173.32L218.988 172.886L220.53 172.669L223.614 172.233L225.157 172.015L226.699 171.796L228.241 171.578L231.325 171.141L232.867 170.923L234.41 170.705L235.952 170.488L237.494 170.27L239.036 170.053L242.12 169.62L243.663 169.404L245.205 169.189L249.831 168.547L252.916 168.121L254.458 167.91L256 167.699L257.542 167.489L259.084 167.28L262.169 166.864L265.253 166.452L269.88 165.841L271.422 165.639L272.964 165.438L276.048 165.04L277.59 164.842L279.133 164.645L280.675 164.45L282.217 164.256L283.759 164.063L285.301 163.871L286.843 163.68L288.386 163.491L289.928 163.302L291.47 163.115L293.012 162.929L294.554 162.744L296.096 162.561L302.265 161.84L303.807 161.663L305.349 161.487L306.892 161.313L314.602 160.46L316.145 160.294L317.687 160.129L319.229 159.965L322.313 159.641L323.855 159.482L325.398 159.323L326.94 159.166L328.482 159.011L330.024 158.857L331.566 158.704L333.108 158.552L334.651 158.402L336.193 158.253L339.277 157.96L340.819 157.816L342.361 157.672L343.904 157.53L345.446 157.39L346.988 157.251L350.072 156.977L353.157 156.708L363.952 155.81L365.494 155.688L368.578 155.446L374.747 154.978L377.831 154.752L379.373 154.641L382.458 154.422L384 154.315'/>
+  </g>
+</svg>
+</div>
+</div>
+#+end_export
+
+#+begin_export html
+</details>
+#+end_export
+
+[fn:1] Recasting of the [[https://en.wikipedia.org/wiki/Schr%C3%B6dinger_equation#Time-independent_equation][TISE]] into a set of coupled integro-differential equations. Derived by optimizing
+the expectation value of the energy subject to normalization constraints, then discretizing it using a suitable
+basis set.
+[fn:2] It may not look like much, but helonium was the very [[https://www.scientificamerican.com/article/the-first-molecule-in-the-universe/][first molecule]] formed in the universe.
+[fn:3] This program can compute the Hartree-Fock energy of any two-electron diatomic molecule.
+[fn:4] STO: functions of the form \(r^le^{-\zeta r}Y_l^m(\theta, \phi)\). For \(1s\) orbitals the
+spherical harmonics integrate out to 1.
+[fn:5] STO-nG: a non-linear least-squares fit of an STO as a weighted sum of n Gaussians.
+[fn:6] See for example [[https://arxiv.org/abs/2007.12057][arXiv:2007.12057]].
+[fn:7] Paul A.M. Dirac, 27 November, 1975
+
+#+INCLUDE: "../html-foot.org"