mirror of
https://github.com/LCPQ/quantum_package
synced 2024-12-23 04:43:50 +01:00
Corrected bug when writing pseudo to EZFIO
This commit is contained in:
parent
ccc061fbeb
commit
d94138cfed
@ -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
|
||||
|
||||
|
||||
|
||||
|
||||
@ -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 ()
|
||||
|
||||
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user