diff --git a/ocaml/Input_mo_basis.ml b/ocaml/Input_mo_basis.ml index 78ab08ee..c51a7f4a 100644 --- a/ocaml/Input_mo_basis.ml +++ b/ocaml/Input_mo_basis.ml @@ -2,22 +2,28 @@ open Qptypes open Qputils open Core -type t_mo = - { 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] module Mo_basis : sig - 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] 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 () = diff --git a/ocaml/qp_find_pi_space.ml b/ocaml/qp_find_pi_space.ml new file mode 100644 index 00000000..86c1401f --- /dev/null +++ b/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 + diff --git a/ocaml/qp_set_mo_class.ml b/ocaml/qp_set_mo_class.ml index 2db69ed1..e82610d4 100644 --- a/ocaml/qp_set_mo_class.ml +++ b/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) diff --git a/scripts/ezfio_interface/qp_convert_output_to_ezfio.py b/scripts/ezfio_interface/qp_convert_output_to_ezfio.py index 6de221f3..7f4f30be 100755 --- a/scripts/ezfio_interface/qp_convert_output_to_ezfio.py +++ b/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 [--ezfio=] + qp_convert_output_to_ezfio.py [-o ] Option: file.out is the file to check (like gamess.out) @@ -429,8 +429,8 @@ if __name__ == '__main__': file_ = get_full_path(arguments['']) - if arguments["--ezfio"]: - ezfio_file = get_full_path(arguments["--ezfio"]) + if arguments["-o"]: + ezfio_file = get_full_path(arguments[""]) else: ezfio_file = "{0}.ezfio".format(file_)