Missing file

CI/CIMatrixElementF12.ml

@@ -0,0 +1,147 @@
+open Lacaml.D
+
+module De = Determinant
+module Ex = Excitation
+module Sp = Spindeterminant
+
+type t = float list
+
+
+let non_zero integrals degree_a degree_b ki kj =
+  let kia = De.alfa ki and kib = De.beta ki
+  and kja = De.alfa kj and kjb = De.beta kj
+  in
+  let diag_element =
+    let mo_a = Sp.to_list kia
+    and mo_b = Sp.to_list kib
+    in
+    fun one_e two_e ->
+      let one =
+        (List.fold_left (fun accu i -> accu +. one_e i i Spin.Alfa) 0. mo_a)
+        +.
+        (List.fold_left (fun accu i -> accu +. one_e i i Spin.Beta) 0. mo_b)
+      in
+      let two =
+        let rec aux_same spin accu = function
+          | [] -> accu
+          | i :: rest -> 
+            let new_accu = 
+              List.fold_left (fun accu j -> accu +. two_e i j i j spin spin) accu rest
+            in
+            aux_same spin new_accu rest
+        in
+        let rec aux_opposite accu other = function
+          | [] -> accu
+          | i :: rest ->
+            let new_accu =
+              List.fold_left (fun accu j -> accu +. two_e i j i j Spin.Alfa Spin.Beta) accu other
+            in
+            aux_opposite new_accu other rest
+        in
+        (aux_same Spin.Alfa 0. mo_a)  +.  (aux_same Spin.Beta 0. mo_b)  +.
+        (aux_opposite 0. mo_a mo_b)
+      in
+      one +. two
+  in
+  let result =
+    match degree_a, degree_b with
+    | 1, 1 -> (* alpha-beta double *)
+      begin
+        let ha, pa, phase_a = Ex.single_of_spindet kia kja in
+        let hb, pb, phase_b = Ex.single_of_spindet kib kjb in
+        let s1 = 
+          match phase_a, phase_b with
+          | Phase.Pos, Phase.Pos
+          | Phase.Neg, Phase.Neg -> fun _ two_e ->    two_e ha hb pa pb Spin.Alfa Spin.Beta
+          | Phase.Neg, Phase.Pos                       
+          | Phase.Pos, Phase.Neg -> fun _ two_e -> -. two_e ha hb pa pb Spin.Alfa Spin.Beta
+        in
+        let kk = Determinant.double_excitation Spin.Alfa ha pb Spin.Beta hb pa ki in
+        let kka = De.alfa kk and kkb = De.beta kk in
+        match Spindeterminant.(degree kia kka, degree kib kkb) with
+        | (1,1) ->
+          let s2 = 
+            begin
+              let ha, pa, phase_a = Ex.single_of_spindet kia kka in
+              let hb, pb, phase_b = Ex.single_of_spindet kib kkb in
+              match phase_a, phase_b with
+              | Phase.Pos, Phase.Pos
+              | Phase.Neg, Phase.Neg -> fun _ two_e ->    two_e ha hb pa pb Spin.Alfa Spin.Beta
+              | Phase.Neg, Phase.Pos                       
+              | Phase.Pos, Phase.Neg -> fun _ two_e -> -. two_e ha hb pa pb Spin.Alfa Spin.Beta
+            end
+          in fun x two_e -> 0.75 *. (s1 x two_e) +. 0.25 *. (s2 x two_e)  
+        | _ -> fun x two_e -> 0.
+      end
+
+    | 2, 0 -> (* alpha double *)
+      begin
+        let h1, p1, h2, p2, phase = Ex.double_of_spindet kia kja in
+        match phase with 
+        | Phase.Pos -> fun _ two_e ->    two_e h1 h2 p1 p2 Spin.Alfa Spin.Alfa
+        | Phase.Neg -> fun _ two_e -> -. two_e h1 h2 p1 p2 Spin.Alfa Spin.Alfa
+      end
+
+    | 0, 2 -> (* beta double *)
+      begin
+        let h1, p1, h2, p2, phase = Ex.double_of_spindet kib kjb in
+        match phase with 
+        | Phase.Pos -> fun _ two_e ->    two_e h1 h2 p1 p2 Spin.Beta Spin.Beta
+        | Phase.Neg -> fun _ two_e -> -. two_e h1 h2 p1 p2 Spin.Beta Spin.Beta
+      end
+
+    | 1, 0    (* alpha single *)
+    | 0, 1    (* beta single *)
+           -> fun _ _ -> 0.
+    | 0, 0 -> (* diagonal element *)
+      diag_element
+
+    | _ -> assert false
+  in
+  List.map (fun (one_e, two_e) -> result one_e two_e) integrals
+
+
+let make integrals ki kj =
+  let degree_a, degree_b =
+    De.degrees ki kj 
+  in
+  if degree_a+degree_b > 2 then
+    List.map (fun _ -> 0.) integrals
+  else
+    non_zero integrals degree_a degree_b ki kj
+    
+
+
+
+let make_s2 ki kj =
+  let degree_a = De.degree_alfa ki kj in
+  let kia = De.alfa ki in
+  let kja = De.alfa kj in
+  if degree_a > 1 then 0.
+  else
+    let degree_b = De.degree_beta ki kj in
+    let kib = De.beta ki in
+    let kjb = De.beta kj in
+    match degree_a, degree_b with
+    | 1, 1 -> (* alpha-beta double *)
+      let ha, pa, phase_a = Ex.single_of_spindet kia kja in
+      let hb, pb, phase_b = Ex.single_of_spindet kib kjb in
+      if ha = pb && hb = pa then
+        begin
+          match phase_a, phase_b with
+          | Phase.Pos, Phase.Pos
+          | Phase.Neg, Phase.Neg -> -1.
+          | Phase.Neg, Phase.Pos                       
+          | Phase.Pos, Phase.Neg -> 1.
+        end
+      else 0.
+    | 0, 0 -> 
+      let ba = Sp.bitstring kia and bb = Sp.bitstring kib in
+      let tmp = Bitstring.logxor ba bb in
+      let n_a = Bitstring.logand ba tmp |> Bitstring.popcount in
+      let n_b = Bitstring.logand bb tmp |> Bitstring.popcount in
+      let s_z = 0.5 *. float_of_int (n_a - n_b) in
+      float_of_int n_a +. s_z *. (s_z -. 1.)
+    | _ -> 0.
+
+