Fix 3e contributions

CI/F12CI.ml

@@ -70,7 +70,7 @@ let single_matrices hf12_integrals density =
     done;
 *)
   
-let hf_ij_non_zero hf12_integrals deg_a deg_b ki kj =
+let hf_ij_non_zero mo_basis hf12_integrals deg_a deg_b ki kj =
   let integrals = [
 
     let { HF12.
@@ -85,26 +85,20 @@ let hf_ij_non_zero hf12_integrals deg_a deg_b ki kj =
     let mo_a =
       Bitstring.logand (Sp.bitstring kia) (Sp.bitstring kja)
       |> Bitstring.to_list
+      |> Array.of_list
     and mo_b =
       Bitstring.logand (Sp.bitstring kib) (Sp.bitstring kjb)
       |> Bitstring.to_list
+      |> Array.of_list
+    in
+
+    let aux_mos = 
+      Util.list_range (MOBasis.size mo_basis) (MOBasis.size aux_basis)
+      |> Array.of_list
     in
 
     let one_e _ _ _ = 0. in
 
-    let two_e i j k l s s' =
-      if s = s' then
-	hf12_anti.{i,j,k,l} -. (
-	(List.fold_left (fun accu m -> accu +. hf12_single_anti.{m,i,j,k,l}) 0. mo_a) +.
-	(List.fold_left (fun accu m -> accu +. hf12_single_anti.{m,j,i,l,k}) 0. mo_b) 
-        )
-      else
-	hf12.{i,j,k,l}  -. (
-	(List.fold_left (fun accu m -> accu +. hf12_single.{m,i,j,k,l}) 0. mo_a) +.
-	(List.fold_left (fun accu m -> accu +. hf12_single.{m,j,i,l,k}) 0. mo_b) 
-        )
-    in
-
     let h = MOBasis.ee_ints aux_basis in
     let two_e_h i j k l s s' =
       if s' <> s then
@@ -115,24 +109,48 @@ let hf_ij_non_zero hf12_integrals deg_a deg_b ki kj =
 
     let f = MOBasis.f12_ints aux_basis in
     let two_e_f i j k l s s' =
+      let fijkl = F12.get_phys f i j k l
+      and fijlk = F12.get_phys f i j l k
+      in
       if s' <> s then
-        F12.get_phys f i j k l
+        0.325 *. fijkl +. 0.125 *. fijlk
       else
-        (F12.get_phys f i j k l) -.  (F12.get_phys f i j l k)
+        0.25 *. (fijkl -. fijlk)
     in
 
+    (*
     let mo_of_s = function
     | Spin.Alfa -> mo_a
     | Spin.Beta -> mo_b
     in
+    *)
+
+    let two_e i j k l s s' =
+      (
+      if s = s' then
+	hf12_anti.{i,j,k,l} -. (
+	(Array.fold_left (fun accu m -> accu +. hf12_single_anti.{m,i,j,k,l}) 0. mo_a) +.
+	(Array.fold_left (fun accu m -> accu +. hf12_single_anti.{m,j,i,l,k}) 0. mo_b) )
+      else
+	hf12.{i,j,k,l}  -. (
+	(Array.fold_left (fun accu m -> accu +. hf12_single.{m,i,j,k,l}) 0. mo_a) +.
+	(Array.fold_left (fun accu m -> accu +. hf12_single.{m,j,i,l,k}) 0. mo_b) )
+      )
+      (*
+      +. Array.fold_left ( fun accu a -> accu +.
+         (Array.fold_left ( fun accu m -> accu +. two_e_h m i m a s s) 0. (mo_of_s s) +.
+          Array.fold_left ( fun accu m -> accu +. two_e_h m i m a s s) 0. (mo_of_s @@ Spin.other s) )
+         *. (two_e_f a j k l s s') ) 0. aux_mos 
+         *)
+    in
 
     let three_e i j k l m n s s' s'' =
-         List.fold_left (fun accu a -> accu +. two_e_h i j l a s   s'  *. two_e_f a k m n s'  s'') 0. (mo_of_s s' )
+         Array.fold_left (fun accu a -> accu +. two_e_h i j l a s   s'  *. two_e_f a k m n s'  s'') 0. aux_mos
-      -. List.fold_left (fun accu a -> accu +. two_e_h j i m a s'  s   *. two_e_f a k l n s   s'') 0. (mo_of_s s  )
+      -. Array.fold_left (fun accu a -> accu +. two_e_h j i m a s'  s   *. two_e_f a k l n s   s'') 0. aux_mos
-      +. List.fold_left (fun accu a -> accu +. two_e_h j k m a s'  s'' *. two_e_f a i n l s'' s  ) 0. (mo_of_s s'')
+      +. Array.fold_left (fun accu a -> accu +. two_e_h j k m a s'  s'' *. two_e_f a i n l s'' s  ) 0. aux_mos
-      -. List.fold_left (fun accu a -> accu +. two_e_h k j n a s'' s'  *. two_e_f a i m l s'  s  ) 0. (mo_of_s s' )
+      -. Array.fold_left (fun accu a -> accu +. two_e_h k j n a s'' s'  *. two_e_f a i m l s'  s  ) 0. aux_mos
-      +. List.fold_left (fun accu a -> accu +. two_e_h k i n a s'' s   *. two_e_f a j l m s   s' ) 0. (mo_of_s s  )
+      +. Array.fold_left (fun accu a -> accu +. two_e_h k i n a s'' s   *. two_e_f a j l m s   s' ) 0. aux_mos
-      -. List.fold_left (fun accu a -> accu +. two_e_h i k l a s   s'' *. two_e_f a j n m s'' s' ) 0. (mo_of_s s'')
+      -. Array.fold_left (fun accu a -> accu +. two_e_h i k l a s   s'' *. two_e_f a j n m s'' s' ) 0. aux_mos
     in
     (one_e, two_e, Some three_e)
   ]
@@ -148,6 +166,9 @@ let dressing_vector ~frozen_core hf12_integrals f12_amplitudes ci =
   let det_space = 
     ci.CI.det_space 
   in
+  let mo_basis =
+    Ds.mo_basis det_space
+  in
 
   let m_HF = 
     
@@ -157,7 +178,7 @@ let dressing_vector ~frozen_core hf12_integrals f12_amplitudes ci =
       | Ds.Spin      _ -> CI.create_matrix_spin_computed ~nmax:3
     in
     f (fun deg_a deg_b ki kj ->
-         hf_ij_non_zero hf12_integrals deg_a deg_b ki kj
+         hf_ij_non_zero mo_basis hf12_integrals deg_a deg_b ki kj
       ) det_space 
 
   in