From 963441b434c156a5d0f7ed30f93dab66ba6c5ad1 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Tue, 19 Mar 2019 00:10:34 +0100
Subject: [PATCH] CAS+EN OK for 1 det

---
 CI/CI.ml   | 58 +++++++++++++++++++++++++++++++++---------------------
 run_cas.ml |  8 +++++++-
 2 files changed, 43 insertions(+), 23 deletions(-)

diff --git a/CI/CI.ml b/CI/CI.ml
index b6531d0..ac7d336 100644
--- a/CI/CI.ml
+++ b/CI/CI.ml
@@ -442,8 +442,32 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
         *)
 
 
+  let symmetric = i_o1_alfa == alfa_o2_i in
+
+  let is_internal alfa =
+    let rec aux = function
+      | -1 -> false
+      | j  -> if (Determinant.degree (fst psi0.(j)) alfa = 0) then true
+              else aux (j-1)
+    in
+    aux (Array.length psi0 - 1)
+  in
+
+
   let det_contribution i =
 
+    let already_generated alfa =
+      if is_internal alfa then
+        true
+      else
+        let rec aux = function
+          | -1 -> false
+          | j  -> if (Determinant.degree (fst psi0.(j)) alfa <= 2) then true
+                  else aux (j-1)
+        in
+        aux (i-1)
+    in
+
     let psi_filtered_idx =
       let rec aux accu =  function
         | j when j < i -> List.rev accu
@@ -458,8 +482,6 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
       List.map (fun i -> psi0.(i)) psi_filtered_idx
     in
 
-    let symmetric = i_o1_alfa == alfa_o2_i in
-
     let psi_h_alfa alfa =
       List.fold_left (fun accu (det, coef) ->
           accu +. coef *. (i_o1_alfa det alfa)) 0. psi_filtered
@@ -481,25 +503,6 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
         (psi_h_alfa alfa) *. (alfa_h_psi alfa)
     in
 
-    let is_internal alfa =
-      let rec aux = function
-        | -1 -> false
-        | j  -> Determinant.degree (fst psi0.(j)) alfa = 0 || aux (j-1)
-      in
-      aux (Array.length psi0 - 1)
-    in
-
-    let already_generated alfa =
-      if is_internal alfa then
-        true
-      else
-        let rec aux = function
-          | -1 -> false
-          | j  -> Determinant.degree (fst psi0.(j)) alfa <= 2 || aux (j-1)
-        in
-        aux (i-1)
-    in
-
     let det_i = fst psi0.(i) in
     let w_alfa = w_alfa det_i in
 
@@ -584,7 +587,7 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
     same_spin +. opposite_spin
   in
 
-  Array.mapi (fun i (_,c_i) -> c_i *. det_contribution i) psi0
+  Array.mapi (fun i (_,c_i) -> det_contribution i) psi0
   |> Array.fold_left (+.) 0.
 
 
@@ -627,3 +630,14 @@ let pt2_mp ci =
   second_order_sum  ci  i_o1_alfa  i_o1_alfa  w_alfa
 
 
+let variance ci =
+
+  let mo_basis  = Ds.mo_basis ci.det_space in
+
+  let i_o1_alfa = h_ij mo_basis in
+
+  let w_alfa _ _ = 1. in
+
+  second_order_sum  ci  i_o1_alfa  i_o1_alfa  w_alfa
+
+
diff --git a/run_cas.ml b/run_cas.ml
index 6c77ab9..75266ea 100644
--- a/run_cas.ml
+++ b/run_cas.ml
@@ -72,8 +72,14 @@ let () =
   let ci = CI.make space in
   Format.fprintf ppf "CAS-CI energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s);
 
+  let pt2 = CI.pt2_mp ci in
+  Format.fprintf ppf "CAS-MP2 energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s +. pt2);
+
   let pt2 = CI.pt2_en ci in
-  Format.fprintf ppf "CAS-EN2 energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s +. pt2)
+  Format.fprintf ppf "CAS-EN2 energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s +. pt2);
+
+  let variance = CI.variance ci in
+  Format.fprintf ppf "CAS variance   : %20.16f@." variance
   (*
   let s2 = Util.xt_o_x ~o:(CI.s2_matrix ci) ~x:(CI.eigenvectors ci) in
   Util.list_range 1 (DeterminantSpace.size space)