From de1791f9d67f092706b4a27fae3e2fa9c1512a56 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Fri, 25 May 2018 14:47:05 +0200
Subject: [PATCH] Find_pi_space updated

---
 ocaml/qp_find_pi_space.ml     | 40 +++++++++++++++++++++++++++++------
 scripts/qp_set_frozen_core.py | 14 +++++++++---
 2 files changed, 45 insertions(+), 9 deletions(-)

diff --git a/ocaml/qp_find_pi_space.ml b/ocaml/qp_find_pi_space.ml
index e23b4382..b3702d42 100644
--- a/ocaml/qp_find_pi_space.ml
+++ b/ocaml/qp_find_pi_space.ml
@@ -2,7 +2,7 @@ open Core
 open Qputils
 open Qptypes
 
-let run ?(sym="None") apply ezfio_filename =
+let run ?(sym="None") ?(apply="no") fc ezfio_filename =
   Ezfio.set_file ezfio_filename ;
 
   let aos = 
@@ -75,11 +75,11 @@ let run ?(sym="None") apply ezfio_filename =
   begin
     match sym with
     | "x" | "X"  | "y" | "Y"  | "z" | "Z" -> 
-        if apply then
+        if apply = "cas" then
           begin
             let ne = Ezfio.get_electrons_elec_alpha_num () in
             let command = 
-              (Printf.sprintf "qp_set_mo_class -act \"[") ^
+              "qp_set_mo_class -act \"[" ^
               (String.concat ~sep:"," @@ List.map ~f:string_of_int pi) ^
               "]\" -core \"[" ^
               (String.concat ~sep:"," @@ List.map ~f:string_of_int @@ List.filter ~f:(fun mo_i -> mo_i <= ne) sigma) ^
@@ -91,6 +91,32 @@ let run ?(sym="None") apply ezfio_filename =
             if Sys.command command <> 0 then
               failwith "Command failed"
           end
+        else if apply = "cassd" then
+          begin
+            let ne = Ezfio.get_electrons_elec_alpha_num () in
+            let nfc = 
+              if fc then
+                let ic =
+                  Unix.open_process_in ("qp_set_frozen_core.py -q "^ezfio_filename)
+                in
+                let result = int_of_string @@ input_line ic in
+                close_in ic; result
+              else 0
+            in
+            let command = 
+              "qp_set_mo_class" ^
+              ( if nfc > 0 then
+                Printf.sprintf " -core \"[1-%d]\"" nfc else "") ^
+              " -inact \"[" ^
+              (String.concat ~sep:"," @@ List.map ~f:string_of_int @@ List.filter ~f:(fun mo_i -> mo_i > nfc && mo_i <= ne) sigma) ^ "]\"" ^
+              " -act \"[" ^ (String.concat ~sep:"," @@ List.map ~f:string_of_int pi) ^ "]\"" ^
+              " -virt \"[" ^
+              (String.concat ~sep:"," @@ List.map ~f:string_of_int @@ List.filter ~f:(fun mo_i -> mo_i > ne) sigma) ^ "]\" " ^ ezfio_filename
+            in
+            print_endline command;
+            if Sys.command command <> 0 then
+              failwith "Command failed"
+          end
         else
           begin
             Printf.printf "Pi: [";
@@ -111,8 +137,10 @@ let spec =
   empty
   +> flag "sym" (optional string)
     ~doc:"{x,y,z} Axis perpendicular to the plane"
-  +> flag "apply" no_arg
-    ~doc:"Set the pi space as the active space, all other MOs frozen"
+  +> flag "apply" (optional string)
+    ~doc:"[cas|cassd] Set the pi space as the active space, all other MOs frozen. If [cassd] is chosen, set inactive and virtual mos."
+  +> flag "fc" no_arg
+    ~doc:"freeze core MOs with -apply cassd"
   +> anon ("ezfio_filename" %: string)
 
 
@@ -124,7 +152,7 @@ let command =
      "Find all the pi molecular orbitals to create a pi space.
      ")
     spec
-    (fun sym apply ezfio_filename () -> run ?sym apply ezfio_filename)
+    (fun sym apply fc ezfio_filename () -> run ?sym ?apply fc ezfio_filename)
 
 
 let () =
diff --git a/scripts/qp_set_frozen_core.py b/scripts/qp_set_frozen_core.py
index 49aaae29..aee32574 100755
--- a/scripts/qp_set_frozen_core.py
+++ b/scripts/qp_set_frozen_core.py
@@ -6,8 +6,16 @@ import sys
 sys.path = [ os.environ["QP_ROOT"]+"/install/EZFIO/Python" ] + sys.path
 from ezfio import ezfio
 
-ezfio.set_filename(sys.argv[1])
-if ezfio.get_pseudo_do_pseudo:
+
+filename = sys.argv[1]
+if filename == '-q': filename = sys.argv[2]
+
+ezfio.set_filename(filename)
+
+if ezfio.pseudo_do_pseudo:
+  if '-q' in sys.argv:
+    print 0
+    sys.exit(0)
   print "I will not set frozen core with pseudopotentials."
   sys.exit(0)
 
@@ -23,7 +31,7 @@ for charge in ezfio.nuclei_nucl_charge:
 mo_tot_num = ezfio.mo_basis_mo_tot_num
 
 if len(sys.argv)>2:
-  if sys.argv[2] == '-q':
+  if '-q' in sys.argv:
     print nb
     sys.exit(0)