From 8e7107cfbaef0bf38fd7dbb251dba6ae76a931a5 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Thu, 21 Mar 2019 11:02:58 +0100
Subject: [PATCH] Working on F12

---
 CI/F12CI.ml      | 55 ++++++++++++++++++++++++++++++++++++++++++------
 Utils/Matrix.ml  | 17 ++++++++++++++-
 Utils/Matrix.mli |  3 +++
 3 files changed, 67 insertions(+), 8 deletions(-)

diff --git a/CI/F12CI.ml b/CI/F12CI.ml
index d02a831..b2240b4 100644
--- a/CI/F12CI.ml
+++ b/CI/F12CI.ml
@@ -6,6 +6,7 @@ type t =
     aux_basis : MOBasis.t ;
     det_space : DeterminantSpace.t ;
     ci        : CI.t ;
+    f12_amplitudes : Mat.t;
   }
 
 let ci t        =  t.ci
@@ -13,6 +14,7 @@ let mo_basis t  =  t.mo_basis
 let det_space t =  t.det_space
 let mo_class t  =  DeterminantSpace.mo_class @@ det_space t
 
+
 let f12_integrals mo_basis =
   let two_e_ints  = MOBasis.two_e_ints mo_basis in
   ( (fun i j _ -> 0.),
@@ -26,6 +28,7 @@ let f12_integrals mo_basis =
 
 
 
+
 let h_ij mo_basis ki kj =
   let integrals =
     List.map (fun f -> f mo_basis)
@@ -54,15 +57,29 @@ let dressing_vector ci =
   let w_alfa _ _ = 1. in
 
   let mo_class = CI.mo_class ci in
-  let list_holes = List.concat
-    [ MOClass.inactive_mos mo_class ; MOClass.active_mos mo_class ]
-  and list_particles = MOClass.auxiliary_mos mo_class
+
+  let internal_term =
+    let list_holes = List.concat
+      [ MOClass.inactive_mos mo_class ; MOClass.active_mos mo_class ]
+    and list_particles = List.concat
+      [ MOClass.active_mos mo_class ; MOClass.virtual_mos mo_class ]
+    in
+    CI.second_order_sum  ci  list_holes  list_particles
+      i_o1_alfa  alfa_o2_i  w_alfa
+    |> Vec.of_list
   in
 
-  CI.second_order_sum  ci  list_holes  list_particles
-    i_o1_alfa  alfa_o2_i  w_alfa
-  |> Vec.of_list
+  let external_term = 
+    let list_holes = List.concat
+      [ MOClass.inactive_mos mo_class ; MOClass.active_mos mo_class ]
+    and list_particles = MOClass.auxiliary_mos mo_class
+    in
+    CI.second_order_sum  ci  list_holes  list_particles
+      i_o1_alfa  alfa_o2_i  w_alfa
+    |> Vec.of_list
+  in
 
+  Vec.sub external_term internal_term
 
 
 
@@ -97,6 +114,30 @@ let make ~simulation ?(frozen_core=true) ~mo_basis ~aux_basis_filename () =
   in
 
   let ci = CI.make det_space in
-  { mo_basis ; aux_basis ; det_space ; ci }
+
+  let ci_coef, ci_energy = Parallel.broadcast ci.eigensystem in
+
+  let f12_amplitudes = 
+    (* While in a sequential region, initiate the parallel
+       4-idx transformation to avoid nested parallel jobs
+    *)
+    ignore @@ MOBasis.f12_ints mo_basis;
+
+    let f = fun ki kj ->
+         if ki <> kj then
+            f_ij mo_basis ki kj
+         else
+            f_ij mo_basis ki kj +. 1.
+    in
+    let m_F = 
+      CI.create_matrix_spin f det_space
+      |> Lazy.force
+    in
+    Matrix.ax_eq_b (Matrix.dense_of_sparse m_F) (Matrix.dense_of_mat ci_coef)
+    |> Matrix.to_mat
+  in
+
+
+  { mo_basis ; aux_basis ; det_space ; ci ; f12_amplitudes }
 
 
diff --git a/Utils/Matrix.ml b/Utils/Matrix.ml
index 690405f..4370796 100644
--- a/Utils/Matrix.ml
+++ b/Utils/Matrix.ml
@@ -306,8 +306,23 @@ let mv ?(sparse=false) ?(trans=`N) ?(threshold=epsilon) a b =
     Vector.sparse_of_vec dense_result
   else
     Vector.dense_of_vec  dense_result
-      
 
+let iterative_ax_eq_b ~trans a b = 
+  failwith "Not implemented"
+
+let rec ax_eq_b ?(trans=`N) a b =
+  match a, b with 
+  | (Dense  a), (Dense  b) ->
+      let x = lacpy a in
+      (getrs ~trans x b; Dense x)
+  | (Dense  _), (Sparse _) -> 
+      let b = dense_of_sparse b in
+      ax_eq_b ~trans a b
+  | _ -> iterative_ax_eq_b ~trans a b
+
+
+
+(* ------------ Printers ------------ *)
 
 let rec pp_matrix ppf = function
   | Dense  m -> Util.pp_matrix ppf m
diff --git a/Utils/Matrix.mli b/Utils/Matrix.mli
index 0d37233..e743a9f 100644
--- a/Utils/Matrix.mli
+++ b/Utils/Matrix.mli
@@ -60,6 +60,9 @@ val mm : ?transa:trans3 -> ?transb:trans3 -> ?threshold:float -> t -> t -> t
 val mv : ?sparse:bool -> ?trans:trans3 -> ?threshold:float -> t -> Vector.t -> Vector.t
 (** Matrix Vector product *)
 
+val ax_eq_b : ?trans:trans3 -> t -> t -> t
+(** Solves A.X = B or A'.X = B *)
+
 
 (** {1 Printers } *)