Added qp_find_pi_space.ml

ocaml/Input_mo_basis.ml

@@ -2,7 +2,9 @@ open Qptypes
 open Qputils
 open Core
 
-type t_mo = 
+
+module Mo_basis : sig
+  type t = 
       { mo_tot_num      : MO_number.t ;
         mo_label        : MO_label.t;
         mo_class        : MO_class.t array;
@@ -10,14 +12,18 @@ type t_mo =
         mo_coef         : (MO_coef.t array) array;
         ao_md5          : MD5.t;
       } [@@deriving sexp]
-
-module Mo_basis : sig
-  type t = t_mo
   val read : unit -> t option
   val to_string : t -> string
   val to_rst : t -> Rst_string.t
 end = struct
-  type t = t_mo
+  type t = 
+      { mo_tot_num      : MO_number.t ;
+        mo_label        : MO_label.t;
+        mo_class        : MO_class.t array;
+        mo_occ          : MO_occ.t array;
+        mo_coef         : (MO_coef.t array) array;
+        ao_md5          : MD5.t;
+      } [@@deriving sexp]
   let get_default = Qpackage.get_ezfio_default "mo_basis"
 
   let read_mo_label () = 

ocaml/qp_find_pi_space.ml

@@ -0,0 +1,109 @@
+open Core
+open Qputils
+open Qptypes
+
+let run ?(sym="None") ezfio_filename =
+  Ezfio.set_file ezfio_filename ;
+
+  let aos = 
+    match Input.Ao_basis.read () with
+    | Some aos -> aos
+    | None -> failwith "Unable to read AOs"
+  and mos = 
+    match Input.Mo_basis.read () with
+    | Some mos -> mos
+    | None -> failwith "Unable to read MOs"
+  in
+  let rec find_power symmetry accu = function
+    | -1 -> accu
+    | i  -> let new_accu = 
+              if aos.Input.Ao_basis.ao_power.(i) = symmetry then (i::accu) 
+              else accu
+           in
+           find_power symmetry new_accu (i-1)
+  and n =
+    (AO_number.to_int aos.Input.Ao_basis.ao_num) 
+  and m =
+    (MO_number.to_int mos.Input.Mo_basis.mo_tot_num) 
+  and one  = Positive_int.of_int 1
+  and zero =  Positive_int.of_int 0
+  in
+
+  (* Indices of the px, py and pz AOs *)
+  let x_indices =
+    find_power Symmetry.Xyz.{x=one  ; y=zero ; z=zero} [] (n-1)
+  and y_indices =
+    find_power Symmetry.Xyz.{x=zero ; y=one  ; z=zero} [] (n-1)
+  and z_indices =
+    find_power Symmetry.Xyz.{x=zero ; y=zero ; z=one } [] (n-1)
+  in
+
+  (* Compute the relative weight of each MO on the px, py, pz spaces *)
+  let compute_weight mo_i list_aos =
+    let num =
+       List.fold_left ~f:(fun s i -> s +. (MO_coef.to_float @@ mos.Input.Mo_basis.mo_coef.(mo_i).(i)) ** 2.) ~init:0. list_aos
+    and denom =
+       Array.fold ~f:(fun s x -> s +. (MO_coef.to_float x) ** 2.) ~init:0. mos.Input.Mo_basis.mo_coef.(mo_i) 
+    in
+    num /. denom
+  in
+  let result = 
+    Array.init ~f:(fun mo_i -> 
+      (mo_i,
+      compute_weight mo_i x_indices,
+      compute_weight mo_i y_indices,
+      compute_weight mo_i z_indices) ) m
+    |> Array.to_list
+  in
+
+  let pi, sigma =
+    let rec aux test_xyz (accu_pi: int list) (accu_sigma: int list) = function
+    | [] -> (List.rev accu_pi, List.rev accu_sigma)
+    | (mo_i,x,y,z)::rest ->
+        if test_xyz (x,y,z) then
+           aux test_xyz (mo_i::accu_pi) accu_sigma rest
+         else
+           aux test_xyz accu_pi (mo_i::accu_sigma) rest
+    in
+    match sym with
+    | "x" | "X" -> aux (fun (x,y,z) -> (x>y && x>z)) [] [] result
+    | "y" | "Y" -> aux (fun (x,y,z) -> (y>x && y>z)) [] [] result
+    | "z" | "Z" -> aux (fun (x,y,z) -> (z>x && z>y)) [] [] result
+    | _         -> ([],[])
+  in
+
+  match sym with
+  | "x" | "X"  | "y" | "Y"  | "z" | "Z" -> 
+      begin
+        Printf.printf "Pi: [";
+        List.iter ~f:(fun mo_i -> Printf.printf "%d," mo_i) pi;
+        Printf.printf "\b]\n\nSigma: [";
+        List.iter ~f:(fun mo_i -> Printf.printf "%d," mo_i) sigma;
+        Printf.printf "\b]\n"
+      end
+  | _ -> List.iter ~f:(fun (mo_i,x,y,z) -> Printf.printf "%d: (%f,%f,%f)\n" mo_i x y z) result
+  
+
+
+
+
+let spec =
+  let open Command.Spec in
+  empty
+  +> flag "sym" (optional string) ~doc:"{x,y,z} Axis perpendicular to the plane"
+  +> anon ("ezfio_filename" %: string)
+
+
+let command =
+    Command.basic
+    ~summary: "Quantum Package command"
+    ~readme:(fun () ->
+     "Find all the pi molecular orbitals to create a pi space.
+     ")
+    spec
+    (fun sym ezfio_filename () -> run ?sym ezfio_filename )
+
+
+let () =
+    Command.run command
+

ocaml/qp_set_mo_class.ml

@@ -210,7 +210,7 @@ let get () =
   in
 
   let mo_tot_num =
-    MO_number.to_int data.Input_mo_basis.mo_tot_num
+    MO_number.to_int data.Input.Mo_basis.mo_tot_num
   in
 
   let n_int =
@@ -244,7 +244,7 @@ let get () =
     | (MO_class.Deleted  _) :: rest ->
         work ~del:(Printf.sprintf "%s,%d" del   i) ~inact ~act  ~virt ~core (i+1) rest
   in
-  work 1 (Array.to_list data.Input_mo_basis.mo_class)
+  work 1 (Array.to_list data.Input.Mo_basis.mo_class)
 
 
 

scripts/ezfio_interface/qp_convert_output_to_ezfio.py

@@ -3,7 +3,7 @@
 convert output of gamess/GAU$$IAN to ezfio
 
 Usage:
-    qp_convert_output_to_ezfio.py <file.out> [--ezfio=<ezfio_directory>]
+    qp_convert_output_to_ezfio.py <file.out> [-o <ezfio_directory>]
 
 Option:
     file.out is the file to check (like gamess.out)
@@ -429,8 +429,8 @@ if __name__ == '__main__':
 
     file_ = get_full_path(arguments['<file.out>'])
 
-    if arguments["--ezfio"]:
+    if arguments["-o"]:
-        ezfio_file = get_full_path(arguments["--ezfio"])
+        ezfio_file = get_full_path(arguments["<ezfio_directory>"])
     else:
         ezfio_file = "{0}.ezfio".format(file_)