feat: implemented up to V the integrals

src/bqn/hf.bqn

@@ -1,9 +1,12 @@
 Sin‿Cos‿Tan‿Erf ← •math
+F ← (π÷4)⊸÷×⟜Erf○√⊢
 
 STO3G ← {
   c‿e ← ⟨0.316894‿0.381531‿0.109991, 0.109818‿0.405771‿2.22766⟩
-  G ← ⋆⟜(3÷4)×·⋆-⊸×⟜(ט)
+  G ← ⋆⟜0.75×·⋆-⊸×⟜(ט)
   +´c×e⊸G𝕩
 }
 
-1
+S ← {a‿b𝕊𝕩: (1.5⋆˜π÷a+b) × ⋆-𝕩×a(×÷+)b}
+T ← {a‿b𝕊𝕩: f←a(×÷+)b ⋄ ×´⟨⋆-𝕩×f, 1.5⋆˜π÷a+b, (3×f)-2×𝕩×טf⟩}
+V ← {a‿b‿z𝕊r‿rz: ×´⟨-2×z×π÷a+b, F rz×a+b, ⋆-r×a(×÷+)b⟩}

src/hf.org

@@ -14,52 +14,52 @@ Four BQN system functions are required:
 
 #+begin_src bqn :results none :tangle ./bqn/hf.bqn
   Sin‿Cos‿Tan‿Erf ← •math
+  F ← (π÷4)⊸÷×⟜Erf○√⊢
 #+end_src
 
 ** STO-3G
 
 Basis sets are used to transform the PDEs into linear algebra problems. Physical intuition suggests that
-Slater type orbitals ([[https://en.wikipedia.org/wiki/Slater-type_orbital][STO]]s) should be a good choice, but the computation of the integrals is a lot easier
-if we approximate them with Gaussian functions[fn:3], three in our case. Since our molecule has only \(1s\)
+Slater type orbitals[fn:3] should be a good choice, but the computation of the integrals is a lot easier
+if we approximate them with Gaussian functions[fn:4], three in our case. Since our molecule has only \(1s\)
 orbitals:
 
 #+begin_src bqn :results none :tangle ./bqn/hf.bqn
   STO3G ← {
     c‿e ← ⟨0.316894‿0.381531‿0.109991, 0.109818‿0.405771‿2.22766⟩
-    G ← ⋆⟜(3÷4)×·⋆-⊸×⟜(ט)
+    G ← ⋆⟜0.75×·⋆-⊸×⟜(ט)
     +´c×e⊸G𝕩
   }
 #+end_src
 
 ** Electronic integrals
 
-Constructing the ERIs tensor is complicated[fn:4], and the reason of the poor scaling of electronic structure
-methods. Since we are restricting ourselves to \(1s\) orbitals, we can compute the integrals almost analytically.
+Constructing the integrals' tensor is complicated[fn:5] and is the main reason for the poor scaling
+of electronic structure methods. For \(1s\) orbitals, some integrals are analytical (S, T)
+while others lack a closed-form solution (V, ERI):
 
-#+begin_src bqn :tangle ./bqn/hf.bqn
-  1
+#+begin_src bqn :tangle ./bqn/hf.bqn 
+  S ← {a‿b𝕊𝕩: (1.5⋆˜π÷a+b) × ⋆-𝕩×a(×÷+)b}
+  T ← {a‿b𝕊𝕩: f←a(×÷+)b ⋄ ×´⟨⋆-𝕩×f, 1.5⋆˜π÷a+b, (3×f)-2×𝕩×טf⟩}
+  V ← {a‿b‿z𝕊r‿rz: ×´⟨-2×z×π÷a+b, F rz×a+b, ⋆-r×a(×÷+)b⟩}
 #+end_src
 
-#+RESULTS:
-: 1
-
 #+begin_export html
+<br/>
 <details>
-<summary>Further details</summary>
+<summary>Derivation strategy</summary>
 #+end_export
 
-We need to compute the overlap, kinetic energy, nuclear attraction, and four-center integrals.
+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, yields the first two set of integrals. The remaining
+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}
+  I(x) = \int_0^{\infty}{{{e^ {- a\,k^2 }\,\sin \left(k\,x\right)}\over{k}}\;dk} \sim \text{Erf}(x)
 \end{equation*}
 
-Solving the above integral is a fun exercise, with the solution expressed in terms of the error function.
-
 #+begin_export html
 </details>
 #+end_export
@@ -72,6 +72,7 @@ Compare the electronic energy with the one computed using the original [[./supp/
 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] These are called STO-nG: a non-linear least-squares fit of functions \(r^ne^{-\zeta r}\) as
-a weighted sum of n Gaussians.
-[fn:4] See for example [[https://arxiv.org/abs/2007.12057][arXiv:2007.12057]].
+[fn:3] 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:4] STO-nG: a non-linear least-squares fit of an STO as a weighted sum of n Gaussians.
+[fn:5] See for example [[https://arxiv.org/abs/2007.12057][arXiv:2007.12057]].