110 lines
3.1 KiB
OCaml
110 lines
3.1 KiB
OCaml
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+1,
|
|
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_spec
|
|
~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
|
|
|