From d94138cfedeede81283819a3f65de71ed9bd41ba Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 3 Feb 2016 00:37:03 +0100 Subject: [PATCH] Corrected bug when writing pseudo to EZFIO --- ocaml/Qputils.ml | 20 +- ocaml/qp_create_ezfio_from_xyz.ml | 553 ++++++++++++++++-------------- 2 files changed, 308 insertions(+), 265 deletions(-) diff --git a/ocaml/Qputils.ml b/ocaml/Qputils.ml index 7f840fde..36dc1102 100644 --- a/ocaml/Qputils.ml +++ b/ocaml/Qputils.ml @@ -1,4 +1,4 @@ -open Core.Std;; +open Core.Std (* let rec transpose = function @@ -24,5 +24,21 @@ let input_to_sexp s = print_endline ("("^result^")"); "("^result^")" |> Sexp.of_string -;; + +let rmdir dirname = + let rec remove_one dir = + Sys.chdir dir; + Sys.readdir "." + |> Array.iter ~f:(fun x -> + match (Sys.is_directory x, Sys.is_file x) with + | (`Yes, _) -> remove_one x + | (_, `Yes) -> Sys.remove x + | _ -> failwith ("Unable to remove file "^x^".") + ); + Sys.chdir ".."; + Unix.rmdir dir + in + remove_one dirname + + diff --git a/ocaml/qp_create_ezfio_from_xyz.ml b/ocaml/qp_create_ezfio_from_xyz.ml index bb7fd847..a6820d88 100644 --- a/ocaml/qp_create_ezfio_from_xyz.ml +++ b/ocaml/qp_create_ezfio_from_xyz.ml @@ -1,6 +1,6 @@ -open Qputils;; -open Qptypes;; -open Core.Std;; +open Qputils +open Qptypes +open Core.Std let spec = let open Command.Spec in @@ -316,289 +316,316 @@ let run ?o b c d m p cart xyz_file = if Sys.file_exists_exn ezfio_file then failwith (ezfio_file^" already exists"); - (* Create EZFIO *) - Ezfio.set_file ezfio_file; + let write_file () = + (* Create EZFIO *) + Ezfio.set_file ezfio_file; - (* Write Pseudo *) - let pseudo = - List.map nuclei ~f:(fun x -> - match pseudo_channel x.Atom.element with - | Some channel -> Pseudo.read_element channel x.Atom.element - | None -> Pseudo.empty x.Atom.element - ) - in + (* Write Pseudo *) + let pseudo = + List.map nuclei ~f:(fun x -> + match pseudo_channel x.Atom.element with + | Some channel -> Pseudo.read_element channel x.Atom.element + | None -> Pseudo.empty x.Atom.element + ) + in - let molecule = - let n_elec_to_remove = - List.fold pseudo ~init:0 ~f:(fun accu x -> - accu + (Positive_int.to_int x.Pseudo.n_elec)) - in - { Molecule.elec_alpha = - (Elec_alpha_number.to_int molecule.Molecule.elec_alpha) - - n_elec_to_remove/2 - |> Elec_alpha_number.of_int; - Molecule.elec_beta = - (Elec_beta_number.to_int molecule.Molecule.elec_beta) - - (n_elec_to_remove - n_elec_to_remove/2) - |> Elec_beta_number.of_int; - Molecule.nuclei = - let charges = - List.map pseudo ~f:(fun x -> Positive_int.to_int x.Pseudo.n_elec - |> Float.of_int) - |> Array.of_list + let molecule = + let n_elec_to_remove = + List.fold pseudo ~init:0 ~f:(fun accu x -> + accu + (Positive_int.to_int x.Pseudo.n_elec)) in - List.mapi molecule.Molecule.nuclei ~f:(fun i x -> - { x with Atom.charge = (Charge.to_float x.Atom.charge) -. charges.(i) - |> Charge.of_float } - ) - } - in - let nuclei = - molecule.Molecule.nuclei @ dummy - in - - - (* Write Electrons *) - Ezfio.set_electrons_elec_alpha_num ( Elec_alpha_number.to_int - molecule.Molecule.elec_alpha ) ; - Ezfio.set_electrons_elec_beta_num ( Elec_beta_number.to_int - molecule.Molecule.elec_beta ) ; - - (* Write Nuclei *) - let labels = - List.map ~f:(fun x->Element.to_string x.Atom.element) nuclei - and charges = - List.map ~f:(fun x-> Atom.(Charge.to_float x.charge)) nuclei - and coords = - (List.map ~f:(fun x-> x.Atom.coord.Point3d.x) nuclei) @ - (List.map ~f:(fun x-> x.Atom.coord.Point3d.y) nuclei) @ - (List.map ~f:(fun x-> x.Atom.coord.Point3d.z) nuclei) in - let nucl_num = (List.length labels) in - Ezfio.set_nuclei_nucl_num nucl_num ; - Ezfio.set_nuclei_nucl_label (Ezfio.ezfio_array_of_list - ~rank:1 ~dim:[| nucl_num |] ~data:labels); - Ezfio.set_nuclei_nucl_charge (Ezfio.ezfio_array_of_list - ~rank:1 ~dim:[| nucl_num |] ~data:charges); - Ezfio.set_nuclei_nucl_coord (Ezfio.ezfio_array_of_list - ~rank:2 ~dim:[| nucl_num ; 3 |] ~data:coords); - - - (* Write pseudopotential *) - let () = - match p with - | None -> Ezfio.set_pseudo_do_pseudo false - | _ -> Ezfio.set_pseudo_do_pseudo true - in - - let klocmax = - List.fold pseudo ~init:0 ~f:(fun accu x -> - let x = - List.length x.Pseudo.local + { Molecule.elec_alpha = + (Elec_alpha_number.to_int molecule.Molecule.elec_alpha) + - n_elec_to_remove/2 + |> Elec_alpha_number.of_int; + Molecule.elec_beta = + (Elec_beta_number.to_int molecule.Molecule.elec_beta) + - (n_elec_to_remove - n_elec_to_remove/2) + |> Elec_beta_number.of_int; + Molecule.nuclei = + let charges = + List.map pseudo ~f:(fun x -> Positive_int.to_int x.Pseudo.n_elec + |> Float.of_int) + |> Array.of_list + in + List.mapi molecule.Molecule.nuclei ~f:(fun i x -> + { x with Atom.charge = (Charge.to_float x.Atom.charge) -. charges.(i) + |> Charge.of_float } + ) + } in - if (x > accu) then x - else accu - ) - and kmax = - List.fold pseudo ~init:0 ~f:(fun accu x -> - let x = - List.length x.Pseudo.non_local + let nuclei = + molecule.Molecule.nuclei @ dummy in - if (x > accu) then x - else accu - ) - and lmax = - List.fold pseudo ~init:0 ~f:(fun accu x -> - let x = - List.fold x.Pseudo.non_local ~init:0 ~f:(fun accu (x,_) -> + + + (* Write Electrons *) + Ezfio.set_electrons_elec_alpha_num ( Elec_alpha_number.to_int + molecule.Molecule.elec_alpha ) ; + Ezfio.set_electrons_elec_beta_num ( Elec_beta_number.to_int + molecule.Molecule.elec_beta ) ; + + (* Write Nuclei *) + let labels = + List.map ~f:(fun x->Element.to_string x.Atom.element) nuclei + and charges = + List.map ~f:(fun x-> Atom.(Charge.to_float x.charge)) nuclei + and coords = + (List.map ~f:(fun x-> x.Atom.coord.Point3d.x) nuclei) @ + (List.map ~f:(fun x-> x.Atom.coord.Point3d.y) nuclei) @ + (List.map ~f:(fun x-> x.Atom.coord.Point3d.z) nuclei) in + let nucl_num = (List.length labels) in + Ezfio.set_nuclei_nucl_num nucl_num ; + Ezfio.set_nuclei_nucl_label (Ezfio.ezfio_array_of_list + ~rank:1 ~dim:[| nucl_num |] ~data:labels); + Ezfio.set_nuclei_nucl_charge (Ezfio.ezfio_array_of_list + ~rank:1 ~dim:[| nucl_num |] ~data:charges); + Ezfio.set_nuclei_nucl_coord (Ezfio.ezfio_array_of_list + ~rank:2 ~dim:[| nucl_num ; 3 |] ~data:coords); + + + (* Write pseudopotential *) + let () = + match p with + | None -> Ezfio.set_pseudo_do_pseudo false + | _ -> Ezfio.set_pseudo_do_pseudo true + in + + let klocmax = + List.fold pseudo ~init:0 ~f:(fun accu x -> let x = - Positive_int.to_int x.Pseudo.Primitive_non_local.proj + List.length x.Pseudo.local + in + if (x > accu) then x + else accu + ) + and lmax = + List.fold pseudo ~init:0 ~f:(fun accu x -> + let x = + List.fold x.Pseudo.non_local ~init:0 ~f:(fun accu (x,_) -> + let x = + Positive_int.to_int x.Pseudo.Primitive_non_local.proj + in + if (x > accu) then x + else accu + ) in if (x > accu) then x else accu ) in - if (x > accu) then x - else accu - ) - in - - let () = - Ezfio.set_pseudo_pseudo_klocmax klocmax; - Ezfio.set_pseudo_pseudo_kmax kmax; - Ezfio.set_pseudo_pseudo_lmax lmax; - let tmp_array_v_k, tmp_array_dz_k, tmp_array_n_k = - Array.make_matrix ~dimx:klocmax ~dimy:nucl_num 0. , - Array.make_matrix ~dimx:klocmax ~dimy:nucl_num 0. , - Array.make_matrix ~dimx:klocmax ~dimy:nucl_num 0 - in - List.iteri pseudo ~f:(fun j x -> - List.iteri x.Pseudo.local ~f:(fun i (y,c) -> - tmp_array_v_k.(i).(j) <- AO_coef.to_float c; - let y, z = - AO_expo.to_float y.Pseudo.Primitive_local.expo, - R_power.to_int y.Pseudo.Primitive_local.r_power + + let kmax = + Array.init (lmax+1) ~f:(fun i-> + List.map pseudo ~f:(fun x -> + List.filter x.Pseudo.non_local ~f:(fun (y,_) -> + (Positive_int.to_int y.Pseudo.Primitive_non_local.proj) = i) + |> List.length ) + |> List.fold ~init:0 ~f:(fun accu x -> + if accu > x then accu else x) + ) + |> Array.fold ~init:0 ~f:(fun accu i -> + if i > accu then i else accu) + in + + + let () = + Ezfio.set_pseudo_pseudo_klocmax klocmax; + Ezfio.set_pseudo_pseudo_kmax kmax; + Ezfio.set_pseudo_pseudo_lmax lmax; + let tmp_array_v_k, tmp_array_dz_k, tmp_array_n_k = + Array.make_matrix ~dimx:klocmax ~dimy:nucl_num 0. , + Array.make_matrix ~dimx:klocmax ~dimy:nucl_num 0. , + Array.make_matrix ~dimx:klocmax ~dimy:nucl_num 0 in - tmp_array_dz_k.(i).(j) <- y; - tmp_array_n_k.(i).(j) <- z; - ) - ); - let concat_2d tmp_array = - let data = - Array.map tmp_array ~f:Array.to_list - |> Array.to_list - |> List.concat - in - Ezfio.ezfio_array_of_list ~rank:2 ~dim:[|nucl_num ; klocmax|] ~data - in - concat_2d tmp_array_v_k - |> Ezfio.set_pseudo_pseudo_v_k ; - concat_2d tmp_array_dz_k - |> Ezfio.set_pseudo_pseudo_dz_k; - concat_2d tmp_array_n_k - |> Ezfio.set_pseudo_pseudo_n_k; + List.iteri pseudo ~f:(fun j x -> + List.iteri x.Pseudo.local ~f:(fun i (y,c) -> + tmp_array_v_k.(i).(j) <- AO_coef.to_float c; + let y, z = + AO_expo.to_float y.Pseudo.Primitive_local.expo, + R_power.to_int y.Pseudo.Primitive_local.r_power + in + tmp_array_dz_k.(i).(j) <- y; + tmp_array_n_k.(i).(j) <- z; + ) + ); + let concat_2d tmp_array = + let data = + Array.map tmp_array ~f:Array.to_list + |> Array.to_list + |> List.concat + in + Ezfio.ezfio_array_of_list ~rank:2 ~dim:[|nucl_num ; klocmax|] ~data + in + concat_2d tmp_array_v_k + |> Ezfio.set_pseudo_pseudo_v_k ; + concat_2d tmp_array_dz_k + |> Ezfio.set_pseudo_pseudo_dz_k; + concat_2d tmp_array_n_k + |> Ezfio.set_pseudo_pseudo_n_k; - let tmp_array_v_kl, tmp_array_dz_kl, tmp_array_n_kl = - Array.create ~len:(lmax+1) - (Array.make_matrix ~dimx:kmax ~dimy:nucl_num 0. ), - Array.create ~len:(lmax+1) - (Array.make_matrix ~dimx:kmax ~dimy:nucl_num 0. ), - Array.create ~len:(lmax+1) - (Array.make_matrix ~dimx:kmax ~dimy:nucl_num 0 ) - in - List.iteri pseudo ~f:(fun j x -> - List.iteri x.Pseudo.non_local ~f:(fun i (y,c) -> - let k, y, z = - Positive_int.to_int y.Pseudo.Primitive_non_local.proj, - AO_expo.to_float y.Pseudo.Primitive_non_local.expo, - R_power.to_int y.Pseudo.Primitive_non_local.r_power + let tmp_array_v_kl, tmp_array_dz_kl, tmp_array_n_kl = + Array.init (lmax+1) ~f:(fun _ -> + (Array.make_matrix ~dimx:kmax ~dimy:nucl_num 0. )), + Array.init (lmax+1) ~f:(fun _ -> + (Array.make_matrix ~dimx:kmax ~dimy:nucl_num 0. )), + Array.init (lmax+1) ~f:(fun _ -> + (Array.make_matrix ~dimx:kmax ~dimy:nucl_num 0 )) in - tmp_array_v_kl.(k).(i).(j) <- AO_coef.to_float c; - tmp_array_dz_kl.(k).(i).(j) <- y; - tmp_array_n_kl.(k).(i).(j) <- z; - ) - ); - let concat_3d tmp_array = - let data = - Array.map tmp_array ~f:(fun x -> - Array.map x ~f:Array.to_list - |> Array.to_list - |> List.concat) - |> Array.to_list + List.iteri pseudo ~f:(fun j x -> + let last_idx = + Array.create ~len:(lmax+1) 0 + in + List.iter x.Pseudo.non_local ~f:(fun (y,c) -> + let k, y, z = + Positive_int.to_int y.Pseudo.Primitive_non_local.proj, + AO_expo.to_float y.Pseudo.Primitive_non_local.expo, + R_power.to_int y.Pseudo.Primitive_non_local.r_power + in + let i = + last_idx.(k) + in + tmp_array_v_kl.(k).(i).(j) <- AO_coef.to_float c; + tmp_array_dz_kl.(k).(i).(j) <- y; + tmp_array_n_kl.(k).(i).(j) <- z; + last_idx.(k) <- i+1; + ) + ); + let concat_3d tmp_array = + let data = + Array.map tmp_array ~f:(fun x -> + Array.map x ~f:Array.to_list + |> Array.to_list + |> List.concat) + |> Array.to_list + |> List.concat + in + Ezfio.ezfio_array_of_list ~rank:3 ~dim:[|nucl_num ; kmax ; lmax+1|] ~data + in + concat_3d tmp_array_v_kl + |> Ezfio.set_pseudo_pseudo_v_kl ; + concat_3d tmp_array_dz_kl + |> Ezfio.set_pseudo_pseudo_dz_kl ; + concat_3d tmp_array_n_kl + |> Ezfio.set_pseudo_pseudo_n_kl ; + in + + + + + (* Write Basis set *) + let basis = + + let nmax = + Nucl_number.get_max () + in + let rec do_work (accu:(Atom.t*Nucl_number.t) list) (n:int) = function + | [] -> accu + | e::tail -> + let new_accu = + (e,(Nucl_number.of_int ~max:nmax n))::accu + in + do_work new_accu (n+1) tail + in + let result = do_work [] 1 nuclei + |> List.rev + |> List.map ~f:(fun (x,i) -> + try + let e = + match x.Atom.element with + | Element.X -> Element.H + | e -> e + in + Basis.read_element (basis_channel x.Atom.element) i e + with + | End_of_file -> failwith + ("Element "^(Element.to_string x.Atom.element)^" not found in basis set.") + ) |> List.concat + in + (* close all in_channels *) + result in - Ezfio.ezfio_array_of_list ~rank:3 ~dim:[|nucl_num ; kmax ; lmax+1|] ~data - in - concat_3d tmp_array_v_kl - |> Ezfio.set_pseudo_pseudo_v_kl ; - concat_3d tmp_array_dz_kl - |> Ezfio.set_pseudo_pseudo_dz_kl ; - concat_3d tmp_array_n_kl - |> Ezfio.set_pseudo_pseudo_n_kl ; - in - - - - - (* Write Basis set *) - let basis = - - let nmax = - Nucl_number.get_max () - in - let rec do_work (accu:(Atom.t*Nucl_number.t) list) (n:int) = function - | [] -> accu - | e::tail -> - let new_accu = - (e,(Nucl_number.of_int ~max:nmax n))::accu + let long_basis = Long_basis.of_basis basis in + let ao_num = List.length long_basis in + Ezfio.set_ao_basis_ao_num ao_num; + Ezfio.set_ao_basis_ao_basis b; + let ao_prim_num = List.map long_basis ~f:(fun (_,g,_) -> List.length g.Gto.lc) + and ao_nucl = List.map long_basis ~f:(fun (_,_,n) -> Nucl_number.to_int n) + and ao_power= + let l = List.map long_basis ~f:(fun (x,_,_) -> x) in + (List.map l ~f:(fun t -> Positive_int.to_int Symmetry.Xyz.(t.x)) )@ + (List.map l ~f:(fun t -> Positive_int.to_int Symmetry.Xyz.(t.y)) )@ + (List.map l ~f:(fun t -> Positive_int.to_int Symmetry.Xyz.(t.z)) ) in - do_work new_accu (n+1) tail - in - let result = do_work [] 1 nuclei - |> List.rev - |> List.map ~f:(fun (x,i) -> - try - let e = - match x.Atom.element with - | Element.X -> Element.H - | e -> e - in - Basis.read_element (basis_channel x.Atom.element) i e - with - | End_of_file -> failwith - ("Element "^(Element.to_string x.Atom.element)^" not found in basis set.") - ) - |> List.concat - in - (* close all in_channels *) - result - in - let long_basis = Long_basis.of_basis basis in - let ao_num = List.length long_basis in - Ezfio.set_ao_basis_ao_num ao_num; - Ezfio.set_ao_basis_ao_basis b; - let ao_prim_num = List.map long_basis ~f:(fun (_,g,_) -> List.length g.Gto.lc) - and ao_nucl = List.map long_basis ~f:(fun (_,_,n) -> Nucl_number.to_int n) - and ao_power= - let l = List.map long_basis ~f:(fun (x,_,_) -> x) in - (List.map l ~f:(fun t -> Positive_int.to_int Symmetry.Xyz.(t.x)) )@ - (List.map l ~f:(fun t -> Positive_int.to_int Symmetry.Xyz.(t.y)) )@ - (List.map l ~f:(fun t -> Positive_int.to_int Symmetry.Xyz.(t.z)) ) - in - let ao_prim_num_max = List.fold ~init:0 ~f:(fun s x -> - if x > s then x - else s) ao_prim_num - in - let gtos = - List.map long_basis ~f:(fun (_,x,_) -> x) - in - - let create_expo_coef ec = - let coefs = - begin match ec with - | `Coefs -> List.map gtos ~f:(fun x-> - List.map x.Gto.lc ~f:(fun (_,coef) -> AO_coef.to_float coef) ) - | `Expos -> List.map gtos ~f:(fun x-> - List.map x.Gto.lc ~f:(fun (prim,_) -> AO_expo.to_float - prim.Primitive.expo) ) - end + let ao_prim_num_max = List.fold ~init:0 ~f:(fun s x -> + if x > s then x + else s) ao_prim_num in - let rec get_n n accu = function - | [] -> List.rev accu - | h::tail -> - let y = - begin match List.nth h n with - | Some x -> x - | None -> 0. + let gtos = + List.map long_basis ~f:(fun (_,x,_) -> x) + in + + let create_expo_coef ec = + let coefs = + begin match ec with + | `Coefs -> List.map gtos ~f:(fun x-> + List.map x.Gto.lc ~f:(fun (_,coef) -> AO_coef.to_float coef) ) + | `Expos -> List.map gtos ~f:(fun x-> + List.map x.Gto.lc ~f:(fun (prim,_) -> AO_expo.to_float + prim.Primitive.expo) ) end - in - get_n n (y::accu) tail + in + let rec get_n n accu = function + | [] -> List.rev accu + | h::tail -> + let y = + begin match List.nth h n with + | Some x -> x + | None -> 0. + end + in + get_n n (y::accu) tail + in + let rec build accu = function + | n when n=ao_prim_num_max -> accu + | n -> build ( accu @ (get_n n [] coefs) ) (n+1) + in + build [] 0 in - let rec build accu = function - | n when n=ao_prim_num_max -> accu - | n -> build ( accu @ (get_n n [] coefs) ) (n+1) - in - build [] 0 - in - let ao_coef = create_expo_coef `Coefs - and ao_expo = create_expo_coef `Expos + let ao_coef = create_expo_coef `Coefs + and ao_expo = create_expo_coef `Expos + in + let () = + Ezfio.set_ao_basis_ao_prim_num (Ezfio.ezfio_array_of_list + ~rank:1 ~dim:[| ao_num |] ~data:ao_prim_num) ; + Ezfio.set_ao_basis_ao_nucl(Ezfio.ezfio_array_of_list + ~rank:1 ~dim:[| ao_num |] ~data:ao_nucl) ; + Ezfio.set_ao_basis_ao_power(Ezfio.ezfio_array_of_list + ~rank:2 ~dim:[| ao_num ; 3 |] ~data:ao_power) ; + Ezfio.set_ao_basis_ao_coef(Ezfio.ezfio_array_of_list + ~rank:2 ~dim:[| ao_num ; ao_prim_num_max |] ~data:ao_coef) ; + Ezfio.set_ao_basis_ao_expo(Ezfio.ezfio_array_of_list + ~rank:2 ~dim:[| ao_num ; ao_prim_num_max |] ~data:ao_expo) ; + Ezfio.set_ao_basis_ao_cartesian(cart); + in + match Input.Ao_basis.read () with + | None -> failwith "Error in basis" + | Some x -> Input.Ao_basis.write x in let () = - Ezfio.set_ao_basis_ao_prim_num (Ezfio.ezfio_array_of_list - ~rank:1 ~dim:[| ao_num |] ~data:ao_prim_num) ; - Ezfio.set_ao_basis_ao_nucl(Ezfio.ezfio_array_of_list - ~rank:1 ~dim:[| ao_num |] ~data:ao_nucl) ; - Ezfio.set_ao_basis_ao_power(Ezfio.ezfio_array_of_list - ~rank:2 ~dim:[| ao_num ; 3 |] ~data:ao_power) ; - Ezfio.set_ao_basis_ao_coef(Ezfio.ezfio_array_of_list - ~rank:2 ~dim:[| ao_num ; ao_prim_num_max |] ~data:ao_coef) ; - Ezfio.set_ao_basis_ao_expo(Ezfio.ezfio_array_of_list - ~rank:2 ~dim:[| ao_num ; ao_prim_num_max |] ~data:ao_expo) ; - Ezfio.set_ao_basis_ao_cartesian(cart); - in - match Input.Ao_basis.read () with - | None -> failwith "Error in basis" - | Some x -> Input.Ao_basis.write x + try write_file () with + | ex -> + begin + begin + match Sys.is_directory ezfio_file with + | `Yes -> rmdir ezfio_file + | _ -> () + end; + raise ex; + end + in ()