From 2a53df423f2dade6bc7e17943341c2e206f6bb8d Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Thu, 6 Oct 2016 00:36:05 +0200
Subject: [PATCH] Added qp_create_guess

---
 ocaml/qp_create_guess.ml            | 141 ++++++++++++++++++++++++++++
 plugins/MRCC_Utils/mrcc_utils.irp.f |   1 +
 2 files changed, 142 insertions(+)
 create mode 100644 ocaml/qp_create_guess.ml

diff --git a/ocaml/qp_create_guess.ml b/ocaml/qp_create_guess.ml
new file mode 100644
index 00000000..62af57de
--- /dev/null
+++ b/ocaml/qp_create_guess.ml
@@ -0,0 +1,141 @@
+open Qputils
+open Qptypes
+open Core.Std
+
+let run ~multiplicity ezfio_file =
+  if (not (Sys.file_exists_exn ezfio_file)) then
+    failwith ("EZFIO directory "^ezfio_file^" not found");
+  Ezfio.set_file ezfio_file;
+  let d =
+    Input.Determinants_by_hand.read ()
+  in
+  let m =
+    Multiplicity.of_int multiplicity
+  in
+  let ne = 
+    Ezfio.get_electrons_elec_alpha_num () +
+    Ezfio.get_electrons_elec_beta_num () 
+    |> Elec_number.of_int
+  in
+  let alpha, beta = 
+    let (a,b) = 
+      Multiplicity.to_alpha_beta ne m 
+    in
+    (Elec_alpha_number.to_int a, Elec_beta_number.to_int b)
+  in
+  let n_open_shells =
+    alpha - beta
+  in
+  let mo_tot_num = 
+    Ezfio.get_mo_basis_mo_tot_num ()
+  in
+  let build_list_of_dets ne n_closed n_open =
+    let init = 
+      Array.create ~len:n_closed Bit.One 
+      |> Array.to_list
+    in
+    let rec set_electron accu = function
+    | 1 ->  [ Bit.One :: accu ]
+    | i -> 
+          assert (i>1);
+          let rest =
+            set_electron (Bit.Zero :: accu) (i-1)
+          in
+          (Bit.One::accu) :: rest
+    in
+    let rec extend accu = function 
+    | 0 -> List.rev accu
+    | i -> extend (Bit.Zero::accu) (i-1)
+    in
+    let rec set_n_electrons accu imax = function
+    | 0 -> []
+    | 1 -> set_electron accu imax
+    | i ->
+          assert (i>1);
+          let l =
+            set_electron accu (imax-1)
+          in
+          List.map ~f:(fun x -> set_n_electrons x (imax-1) (i-1)) l
+          |> List.concat
+    in
+    set_n_electrons init n_open ne
+          |> List.filter ~f:(fun x -> List.length x <= n_closed+n_open) 
+          |> List.map ~f:(fun x -> extend x (((mo_tot_num-1)/64+1)*64 - List.length x))
+  in
+
+  let alpha_new = 
+    (Elec_number.to_int ne + 1)/2
+  and beta_new = 
+    Elec_number.to_int ne/2
+  in
+  let l_alpha = 
+    build_list_of_dets ((alpha-beta+1)/2) beta n_open_shells
+  in 
+  let l_beta = 
+    if alpha_new = beta_new then
+      l_alpha
+    else
+      build_list_of_dets ((alpha-beta)/2)beta n_open_shells 
+  in
+
+  let n_int = 
+    Bitlist.n_int_of_mo_tot_num mo_tot_num
+  in
+  let determinants = 
+    List.map l_alpha ~f:(fun x -> List.map l_beta ~f:(fun y -> (x,y) ))
+    |> List.concat
+    |> List.map ~f:(fun pair -> Determinant.of_bitlist_couple ~n_int
+      ~alpha:(Elec_alpha_number.of_int alpha_new)
+      ~beta:(Elec_beta_number.of_int beta_new) pair )
+  in
+  let c = 
+    Array.create ~len:(List.length determinants) (Det_coef.of_float 1.)
+  in
+ 
+  determinants 
+  |> List.map ~f:(fun x -> Determinant.to_string ~mo_tot_num:(MO_number.of_int mo_tot_num) x) 
+  |> List.iter ~f:(fun x -> Printf.printf "%s\n\n%!" x);
+
+  let l =
+    List.length determinants
+  in
+  if l > 0 then
+    begin
+      let d = 
+        let s = (Float.of_int (alpha - beta)) *. 0.5 in
+        let open Input.Determinants_by_hand in
+        { d with n_int ; 
+          n_det = Det_number.of_int ~min:1 ~max:l  l;
+          expected_s2 = Positive_float.of_float (s *. (s +. 1.)) ;
+          psi_coef = c;
+          psi_det = Array.of_list determinants;
+        }
+      in
+      Input.Determinants_by_hand.write d;
+      Ezfio.set_determinants_read_wf true
+    end
+  else
+    Ezfio.set_determinants_read_wf false
+  
+
+
+let spec =
+  let open Command.Spec in
+  empty
+  +> flag "m" (required int)
+     ~doc:"int Spin multiplicity"
+  +> anon ("ezfio_file" %: string)
+
+let () = 
+  Command.basic
+  ~summary: "Quantum Package command"
+  ~readme:( fun () -> "
+Creates an open-shell multiplet initial guess\n\n" )
+  spec
+  (fun multiplicity ezfio_file () ->
+    run ~multiplicity ezfio_file
+  )
+  |> Command.run   ~version: Git.sha1   ~build_info: Git.message
+
+
+
diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f
index fcfb5dd1..7687b451 100644
--- a/plugins/MRCC_Utils/mrcc_utils.irp.f
+++ b/plugins/MRCC_Utils/mrcc_utils.irp.f
@@ -587,6 +587,7 @@ end subroutine
   N_ex_exists = 0
 
   n = 0
+  !TODO Openmp
   do i=1, N_det_ref
     do l=1, N_det_non_ref
       call get_excitation(psi_ref(1,1,i), psi_non_ref(1,1,l), exc, degree, phase, N_int)