Swapping of 2 electrons for HRR

Basis/TwoElectronRR.ml

@@ -8,13 +8,13 @@ let cutoff2 = cutoff *. cutoff
 exception NullQuartet
 
 (** Horizontal and Vertical Recurrence Relations (HVRR) *)
-let hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
+let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
     (totAngMom_a_in, totAngMom_b_in, totAngMom_c_in, totAngMom_d_in)
     (maxm, zero_m_array)
     (expo_b, expo_d)
     (expo_inv_p, expo_inv_q)
     (center_ab, center_cd, center_pq)
-    map_1d map_2d
+    map_1d map_2d map_1d' map_2d'
   =
 
   let maxsze = maxm+1 in
@@ -23,6 +23,17 @@ let hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
   and totAngMom_c = Angular_momentum.to_int totAngMom_c_in
   and totAngMom_d = Angular_momentum.to_int totAngMom_d_in
   in
+
+  (* Swap electrons 1 and 2 so that the max angular momentum is on 1 *)
+  if (totAngMom_a+totAngMom_b < totAngMom_c+totAngMom_d) then
+    hvrr_two_e (angMom_c, angMom_d, angMom_a, angMom_b)
+      (totAngMom_c_in, totAngMom_d_in, totAngMom_a_in, totAngMom_b_in)
+      (maxm, zero_m_array)
+      (expo_d, expo_b)
+      (expo_inv_q, expo_inv_p)
+      (center_cd, center_ab, (Coordinate.neg center_pq) )
+      map_1d' map_2d' map_1d map_2d
+  else
     let maxm = totAngMom_a + totAngMom_b + totAngMom_c + totAngMom_d in
     let empty = Array.make (maxm+1) 0. 
     in
@@ -362,6 +373,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
             let d = shell_q.(cd).ShellPair.j in
             let map_1d = Zmap.create (4*maxm) in
             let map_2d = Zmap.create (Array.length class_indices) in
+            let map_1d' = Zmap.create (4*maxm) in
+            let map_2d' = Zmap.create (Array.length class_indices) in
             let norm_coef_scale_q = shell_q.(cd).ShellPair.norm_coef_scale in
             let norm_coef_scale =
               Array.map (fun v1 ->
@@ -435,7 +448,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
                       (Contracted_shell.expo shell_b b, Contracted_shell.expo shell_d d)
                       (shell_p.(ab).ShellPair.expo_inv, shell_q.(cd).ShellPair.expo_inv)
                       (shell_p.(ab).ShellPair.center_ab, shell_q.(cd).ShellPair.center_ab, center_pq)
-                      map_1d map_2d 
+                      map_1d map_2d map_1d' map_2d'
                   in
                   contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral
                 with NullQuartet -> ()

Utils/Coordinate.ml

@@ -54,6 +54,8 @@ let (|-), (|+) =
   )
 
 
+let neg a = -1. |. a
+
 let rec dot p q =
   match (p,q) with 
     | Bohr (x,y,z), Bohr (x',y',z') ->  x*.x' +. y*.y' +. z*.z' 

Utils/Coordinate.mli

@@ -7,6 +7,7 @@ val of_3_floats : float -> float -> float -> [< `Angstrom | `Bohr ] -> t
 val ( |. ) : float -> t -> t
 val ( |- ) : t -> t -> t
 val ( |+ ) : t -> t -> t
+val neg : t -> t 
 val dot : t -> t -> float
 val norm : t -> float
 val to_string : t -> string