Introduced TRR

Basis/TwoElectronRR.ml

@@ -15,19 +15,23 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
     (expo_b, expo_d)
     (expo_inv_p, expo_inv_q)
     (center_ab, center_cd, center_pq)
-    map_1d map_2d map_1d' map_2d'
+    (center_pa, center_qc)
+    map_1d map_2d 
   =
 
   (* Swap electrons 1 and 2 so that the max angular momentum is on 1 *)
+  (*
   if  angMom_a.tot + angMom_b.tot < angMom_c.tot + angMom_d.tot then
     hvrr_two_e (angMom_c, angMom_d, angMom_a, angMom_b)
       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
+      (center_qc, center_pa)
+      map_1d map_2d
 
   else
+  *)
 
     let maxm = angMom_a.tot + angMom_b.tot + angMom_c.tot + angMom_d.tot in
     let maxsze = maxm+1 in
@@ -98,9 +102,9 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
         try Zmap.find map_2d key with
         | Not_found -> 
           let result = 
-            (* Invariant : angMom_c.tot > 0  so cm.tot >= 0 *)
+            (* angMom_c.tot > 0  so cm.tot >= 0 *)
-            let xyz = get_xyz angMom_c in
+            let xyz   = get_xyz angMom_c in
-            let cm = Powers.decr xyz angMom_c in
+            let cm    = Powers.decr xyz angMom_c in
             let cmxyz = Powers.get xyz cm in
             let axyz  = Powers.get xyz angMom_a in
 
@@ -162,19 +166,65 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
     and trr angMom_a angMom_c =
 
       match (angMom_a.tot, angMom_c.tot) with
-      | (i,0) -> if (i>0) then vrr0 angMom_a
+      | (i,0) -> if (i>0) then (vrr0 angMom_a).(0)
-                 else zero_m_array 
+                 else zero_m_array.(0)
       | (_,_) ->
+        (*
         let key = Zkey.of_powers (Zkey.Six (angMom_a, angMom_c)) in
 
         try Zmap.find map_2d key with
         | Not_found -> 
+        *)
           let result = 
-            let xyz = get_xyz angMom_c in
+            let xyz   = get_xyz angMom_c in
-            let cm = Powers.decr xyz angMom_c in
+            let axyz  = Powers.get  xyz angMom_a in
-            let cmxyz = Powers.get xyz cm in
+            let cm    = Powers.decr xyz angMom_c in
-            let axyz  = Powers.get xyz angMom_a in
+            let cmxyz = Powers.get  xyz cm       in
-            *)
+
+            let expo_inv_q_over_p = expo_inv_q /. expo_inv_p in
+            let f =
+              Coordinate.get xyz center_qc +. expo_inv_q_over_p *.
+              (Coordinate.get xyz center_pa)
+            in
+            let result = 
+              if abs_float f < cutoff then 0.  else
+                  let v1 = trr angMom_a cm in
+                  f *. v1
+            in
+            let result = 
+              if axyz < 1 then result else
+                  let f = 0.5 *. (float_of_int axyz) *. expo_inv_q in
+                  if abs_float f < cutoff then result else
+                      let am = Powers.decr xyz angMom_a in
+                      let v2 = trr am cm in
+                      result +. f *. v2
+            in
+            let result = 
+              if cmxyz < 1 then result else
+                  let f = 0.5 *. (float_of_int cmxyz) *. expo_inv_q in
+                  if abs_float f < cutoff then 0.  else
+                      let cmm = Powers.decr xyz cm in
+                      let v3 = trr angMom_a cmm in
+                      result +. f *. v3
+            in
+            let result = 
+              if cmxyz < 0 then result else
+                  let f = -. expo_inv_q_over_p in
+                  let ap = Powers.incr xyz angMom_a in
+                  let v4 = trr ap cm in
+                  result +. v4 *. f
+            in
+            result
+          in
+            (*
+          Zmap.add map_2d key result;
+          *)
+          result
+          *)
+
+
+
+
 
 
 
@@ -183,7 +233,6 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
     and hrr0 angMom_a angMom_b angMom_c =
 
       match angMom_b.tot with
-      | 0 -> (vrr angMom_a angMom_c).(0)
       | 1 ->
         let xyz  = get_xyz angMom_b in
         let ap = Powers.incr xyz angMom_a in
@@ -192,6 +241,10 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
         if (abs_float f2 < cutoff) then v1.(0) else
             let v2 = vrr angMom_a angMom_c in
             v1.(0) +. f2 *. v2.(0)
+      | 0 -> (vrr angMom_a angMom_c).(0)
+      (*
+      | 0 -> trr angMom_a angMom_c
+      *)
       | _ ->
         let xyz  = get_xyz angMom_b in
         let bxyz = Powers.get xyz angMom_b in
@@ -212,6 +265,9 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
       | (_,0) ->
         if (angMom_b.tot = 0) then 
             (vrr angMom_a angMom_c).(0)
+            (*
+            trr angMom_a angMom_c
+            *)
         else
             hrr0 angMom_a angMom_b angMom_c
       | (_,_) ->
@@ -306,8 +362,6 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
             let d = shell_q.ContractedShellPair.shell_pairs.(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 =
               Array.to_list norm_coef_scale_p
               |> List.map (fun v1 ->
@@ -361,6 +415,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
 
                   let norm =  norm_coef_scale.(i) in
                   let coef_prod = coef_prod *. norm in
+
                   let integral = 
                     hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
                       zero_m_array
@@ -368,7 +423,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
                       (shell_p.ContractedShellPair.expo_inv.(ab), 
                        shell_q.ContractedShellPair.expo_inv.(cd) )
                       (sp.(ab).ShellPair.center_ab, sq.(cd).ShellPair.center_ab, center_pq)
-                      map_1d map_2d map_1d' map_2d'
+                      (sp.(ab).ShellPair.center_a , sq.(cd).ShellPair.center_a)
+                      map_1d map_2d 
                   in
                   contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral
                 with NullQuartet -> ()