diff --git a/ocaml/Generic_input_of_rst.ml b/ocaml/Generic_input_of_rst.ml index ec7de8c9..fec7c0da 100644 --- a/ocaml/Generic_input_of_rst.ml +++ b/ocaml/Generic_input_of_rst.ml @@ -1,6 +1,7 @@ open Sexplib open Sexplib.Std open Qptypes +open Qputils let fail_msg str (ex,range) = @@ -25,7 +26,7 @@ let fail_msg str (ex,range) = in let str = String_ext.tr str ~target:'(' ~replacement:' ' |> String_ext.split ~on:')' - |> List.map String_ext.strip + |> list_map String_ext.strip |> List.filter (fun x -> match String_ext.substr_index ~pos:0 ~pattern:"##" x with | None -> false @@ -48,7 +49,7 @@ let of_rst t_of_sexp s = Rst_string.to_string s |> String_ext.split ~on:'\n' |> List.filter (fun line -> String.contains line '=') - |> List.map (fun line -> + |> list_map (fun line -> "("^( String_ext.tr ~target:'=' ~replacement:' ' line )^")" ) diff --git a/ocaml/Input_ao_basis.ml b/ocaml/Input_ao_basis.ml index f0b3e92c..a73d641c 100644 --- a/ocaml/Input_ao_basis.ml +++ b/ocaml/Input_ao_basis.ml @@ -202,14 +202,14 @@ end = struct in let ao_prim_num = Array.to_list ao_prim_num - |> List.map AO_prim_number.to_int + |> list_map AO_prim_number.to_int in Ezfio.set_ao_basis_ao_prim_num (Ezfio.ezfio_array_of_list ~rank:1 ~dim:[| ao_num |] ~data:ao_prim_num) ; let ao_nucl = Array.to_list ao_nucl - |> List.map Nucl_number.to_int + |> list_map Nucl_number.to_int in Ezfio.set_ao_basis_ao_nucl(Ezfio.ezfio_array_of_list ~rank:1 ~dim:[| ao_num |] ~data:ao_nucl) ; @@ -217,9 +217,9 @@ end = struct let ao_power = let l = Array.to_list ao_power in List.concat [ - (List.map (fun a -> Positive_int.to_int a.Symmetry.Xyz.x) l) ; - (List.map (fun a -> Positive_int.to_int a.Symmetry.Xyz.y) l) ; - (List.map (fun a -> Positive_int.to_int a.Symmetry.Xyz.z) l) ] + (list_map (fun a -> Positive_int.to_int a.Symmetry.Xyz.x) l) ; + (list_map (fun a -> Positive_int.to_int a.Symmetry.Xyz.y) l) ; + (list_map (fun a -> Positive_int.to_int a.Symmetry.Xyz.z) l) ] in Ezfio.set_ao_basis_ao_power(Ezfio.ezfio_array_of_list ~rank:2 ~dim:[| ao_num ; 3 |] ~data:ao_power) ; @@ -230,14 +230,14 @@ end = struct let ao_coef = Array.to_list ao_coef - |> List.map AO_coef.to_float + |> list_map AO_coef.to_float in Ezfio.set_ao_basis_ao_coef(Ezfio.ezfio_array_of_list ~rank:2 ~dim:[| ao_num ; ao_prim_num_max |] ~data:ao_coef) ; let ao_expo = Array.to_list ao_expo - |> List.map AO_expo.to_float + |> list_map AO_expo.to_float in Ezfio.set_ao_basis_ao_expo(Ezfio.ezfio_array_of_list ~rank:2 ~dim:[| ao_num ; ao_prim_num_max |] ~data:ao_expo) ; @@ -299,8 +299,8 @@ end = struct | Some (s', g', n') -> if s <> s' || n <> n' then find2 (s,g,n) a (i+1) else - let lc = List.map (fun (prim, _) -> prim) g.Gto.lc - and lc' = List.map (fun (prim, _) -> prim) g'.Gto.lc + let lc = list_map (fun (prim, _) -> prim) g.Gto.lc + and lc' = list_map (fun (prim, _) -> prim) g'.Gto.lc in if lc <> lc' then find2 (s,g,n) a (i+1) else (a.(i) <- None ; i) in @@ -314,14 +314,14 @@ end = struct let of_long_basis long_basis name ao_cartesian = let ao_num = List.length long_basis |> AO_number.of_int in let ao_prim_num = - List.map (fun (_,g,_) -> List.length g.Gto.lc + list_map (fun (_,g,_) -> List.length g.Gto.lc |> AO_prim_number.of_int ) long_basis |> Array.of_list and ao_nucl = - List.map (fun (_,_,n) -> n) long_basis + list_map (fun (_,_,n) -> n) long_basis |> Array.of_list and ao_power = - List.map (fun (x,_,_) -> x) long_basis + list_map (fun (x,_,_) -> x) long_basis |> Array.of_list in let ao_prim_num_max = Array.fold_left (fun s x -> @@ -331,15 +331,15 @@ end = struct in let gtos = - List.map (fun (_,x,_) -> x) long_basis + list_map (fun (_,x,_) -> x) long_basis in let create_expo_coef ec = let coefs = begin match ec with - | `Coefs -> List.map (fun x-> - List.map (fun (_,coef) -> AO_coef.to_float coef) x.Gto.lc ) gtos - | `Expos -> List.map (fun x-> - List.map (fun (prim,_) -> AO_expo.to_float + | `Coefs -> list_map (fun x-> + list_map (fun (_,coef) -> AO_coef.to_float coef) x.Gto.lc ) gtos + | `Expos -> list_map (fun x-> + list_map (fun (prim,_) -> AO_expo.to_float prim.GaussianPrimitive.expo) x.Gto.lc ) gtos end in @@ -450,7 +450,7 @@ Basis set (read-only) :: (string_of_bool b.ao_normalized) (Basis.to_string short_basis |> String_ext.split ~on:'\n' - |> List.map (fun x-> " "^x) + |> list_map (fun x-> " "^x) |> String.concat "\n" ) print_sym @@ -490,16 +490,16 @@ md5 = %s " (AO_basis_name.to_string b.ao_basis) (AO_number.to_string b.ao_num) - (b.ao_prim_num |> Array.to_list |> List.map + (b.ao_prim_num |> Array.to_list |> list_map (AO_prim_number.to_string) |> String.concat ", " ) (AO_prim_number.to_string b.ao_prim_num_max) - (b.ao_nucl |> Array.to_list |> List.map Nucl_number.to_string |> + (b.ao_nucl |> Array.to_list |> list_map Nucl_number.to_string |> String.concat ", ") - (b.ao_power |> Array.to_list |> List.map (fun x-> + (b.ao_power |> Array.to_list |> list_map (fun x-> "("^(Symmetry.Xyz.to_string x)^")" )|> String.concat ", ") - (b.ao_coef |> Array.to_list |> List.map AO_coef.to_string + (b.ao_coef |> Array.to_list |> list_map AO_coef.to_string |> String.concat ", ") - (b.ao_expo |> Array.to_list |> List.map AO_expo.to_string + (b.ao_expo |> Array.to_list |> list_map AO_expo.to_string |> String.concat ", ") (b.ao_cartesian |> string_of_bool) (b.ao_normalized |> string_of_bool) diff --git a/ocaml/Input_determinants_by_hand.ml b/ocaml/Input_determinants_by_hand.ml index dee338e3..fb0aef7f 100644 --- a/ocaml/Input_determinants_by_hand.ml +++ b/ocaml/Input_determinants_by_hand.ml @@ -377,7 +377,7 @@ end = struct (coefs_string i) (Determinant.to_string ~mo_num:mo_num b.psi_det.(i) |> String_ext.split ~on:'\n' - |> List.map (fun x -> " "^x) + |> list_map (fun x -> " "^x) |> String.concat "\n" ) ) @@ -427,7 +427,7 @@ psi_det = %s (b.n_det |> Det_number.to_string) (b.n_states |> States_number.to_string) (b.expected_s2 |> Positive_float.to_string) - (b.state_average_weight |> Array.to_list |> List.map Positive_float.to_string |> String.concat ",") + (b.state_average_weight |> Array.to_list |> list_map Positive_float.to_string |> String.concat ",") (b.psi_coef |> Array.map Det_coef.to_string |> Array.to_list |> String.concat ", ") (b.psi_det |> Array.map (Determinant.to_string ~mo_num) |> Array.to_list @@ -457,7 +457,7 @@ psi_det = %s else ( (String.contains line '=') && (line.[0] = ' ') ) ) - |> List.map (fun line -> + |> list_map (fun line -> "("^( String_ext.tr line ~target:'=' ~replacement:' ' |> String.trim @@ -468,7 +468,7 @@ psi_det = %s (* Handle determinant coefs *) let dets = match ( dets |> String_ext.split ~on:'\n' - |> List.map String.trim + |> list_map String.trim ) with | _::lines -> lines | _ -> failwith "Error in determinants" @@ -481,7 +481,7 @@ psi_det = %s | ""::c::tail -> let c = String_ext.split ~on:'\t' c - |> List.map (fun x -> Det_coef.of_float (Float.of_string x)) + |> list_map (fun x -> Det_coef.of_float (Float.of_string x)) |> Array.of_list in read_coefs (c::accu) tail @@ -499,7 +499,7 @@ psi_det = %s let i = i-1 in - List.map (fun x -> Det_coef.to_string x.(i)) buffer + list_map (fun x -> Det_coef.to_string x.(i)) buffer |> String.concat " " in let rec build_result = function diff --git a/ocaml/Input_mo_basis.ml b/ocaml/Input_mo_basis.ml index 94435349..a4e6176a 100644 --- a/ocaml/Input_mo_basis.ml +++ b/ocaml/Input_mo_basis.ml @@ -257,9 +257,9 @@ mo_coef = %s " (MO_label.to_string b.mo_label) (MO_number.to_string b.mo_num) - (b.mo_class |> Array.to_list |> List.map + (b.mo_class |> Array.to_list |> list_map (MO_class.to_string) |> String.concat ", " ) - (b.mo_occ |> Array.to_list |> List.map + (b.mo_occ |> Array.to_list |> list_map (MO_occ.to_string) |> String.concat ", " ) (b.mo_coef |> Array.map (fun x-> Array.map MO_coef.to_string x |> diff --git a/ocaml/Input_nuclei_by_hand.ml b/ocaml/Input_nuclei_by_hand.ml index f195a2de..44747f20 100644 --- a/ocaml/Input_nuclei_by_hand.ml +++ b/ocaml/Input_nuclei_by_hand.ml @@ -50,7 +50,7 @@ end = struct in let labels = Array.to_list labels - |> List.map Element.to_string + |> list_map Element.to_string in Ezfio.ezfio_array_of_list ~rank:1 ~dim:[| nucl_num |] ~data:labels @@ -70,7 +70,7 @@ end = struct in let charges = Array.to_list charges - |> List.map Charge.to_float + |> list_map Charge.to_float in Ezfio.ezfio_array_of_list ~rank:1 ~dim:[| nucl_num |] ~data:charges @@ -101,9 +101,9 @@ end = struct in let coord = Array.to_list coord in let coord = - (List.map (fun x-> x.Point3d.x) coord) @ - (List.map (fun x-> x.Point3d.y) coord) @ - (List.map (fun x-> x.Point3d.z) coord) + (list_map (fun x-> x.Point3d.x) coord) @ + (list_map (fun x-> x.Point3d.y) coord) @ + (list_map (fun x-> x.Point3d.z) coord) in Ezfio.ezfio_array_of_list ~rank:2 ~dim:[| nucl_num ; 3 |] ~data:coord @@ -159,11 +159,11 @@ nucl_charge = %s nucl_coord = %s " (Nucl_number.to_string b.nucl_num) - (b.nucl_label |> Array.to_list |> List.map + (b.nucl_label |> Array.to_list |> list_map (Element.to_string) |> String.concat ", " ) - (b.nucl_charge |> Array.to_list |> List.map + (b.nucl_charge |> Array.to_list |> list_map (Charge.to_string) |> String.concat ", " ) - (b.nucl_coord |> Array.to_list |> List.map + (b.nucl_coord |> Array.to_list |> list_map (Point3d.to_string ~units:Units.Bohr) |> String.concat "\n" ) ;; @@ -226,11 +226,11 @@ Nuclear coordinates in xyz format (Angstroms) :: let result = { nucl_num = List.length atom_list |> Nucl_number.of_int ~max:nmax; - nucl_label = List.map (fun x -> + nucl_label = list_map (fun x -> x.Atom.element) atom_list |> Array.of_list ; - nucl_charge = List.map (fun x -> + nucl_charge = list_map (fun x -> x.Atom.charge ) atom_list |> Array.of_list ; - nucl_coord = List.map (fun x -> + nucl_coord = list_map (fun x -> x.Atom.coord ) atom_list |> Array.of_list ; } in Some result diff --git a/ocaml/Long_basis.ml b/ocaml/Long_basis.ml index dd5af64a..6c88954a 100644 --- a/ocaml/Long_basis.ml +++ b/ocaml/Long_basis.ml @@ -1,4 +1,5 @@ open Qptypes +open Qputils open Sexplib.Std type t = (Symmetry.Xyz.t * Gto.t * Nucl_number.t ) list [@@deriving sexp] @@ -39,7 +40,7 @@ let to_basis b = let to_string b = - let middle = List.map (fun (x,y,z) -> + let middle = list_map (fun (x,y,z) -> "( "^((string_of_int (Nucl_number.to_int z)))^", "^ (Symmetry.Xyz.to_string x)^", "^(Gto.to_string y) ^" )" diff --git a/ocaml/MO_class.ml b/ocaml/MO_class.ml index f4a2d7d2..090e3adf 100644 --- a/ocaml/MO_class.ml +++ b/ocaml/MO_class.ml @@ -1,4 +1,5 @@ open Qptypes +open Qputils open Sexplib.Std @@ -13,7 +14,7 @@ type t = let to_string x = let print_list l = - let s = List.map (fun x-> MO_number.to_int x |> string_of_int ) l + let s = list_map (fun x-> MO_number.to_int x |> string_of_int ) l |> (String.concat ", ") in "("^s^")" @@ -43,7 +44,7 @@ let of_string s = let _mo_number_list_of_range range = - Range.of_string range |> List.map MO_number.of_int + Range.of_string range |> list_map MO_number.of_int let create_core range = Core (_mo_number_list_of_range range) diff --git a/ocaml/Message.ml b/ocaml/Message.ml index 2ea1d38c..b7d77430 100644 --- a/ocaml/Message.ml +++ b/ocaml/Message.ml @@ -1,5 +1,6 @@ open Sexplib.Std open Qptypes +open Qputils (** New job : Request to create a new multi-tasked job *) @@ -193,12 +194,12 @@ end = struct } let create ~state ~task_ids = { state = State.of_string state ; - task_ids = List.map Id.Task.of_int task_ids + task_ids = list_map Id.Task.of_int task_ids } let to_string x = Printf.sprintf "del_task %s %s" (State.to_string x.state) - (String.concat "|" @@ List.map Id.Task.to_string x.task_ids) + (String.concat "|" @@ list_map Id.Task.to_string x.task_ids) end @@ -219,7 +220,7 @@ end = struct else "done" in Printf.sprintf "del_task_reply %s %s" - more (String.concat "|" @@ List.map Id.Task.to_string x.task_ids) + more (String.concat "|" @@ list_map Id.Task.to_string x.task_ids) end @@ -303,7 +304,7 @@ end = struct "get_tasks_reply ok" let to_string_list x = "get_tasks_reply ok" :: ( - List.map (fun (task_id, task) -> + list_map (fun (task_id, task) -> match task_id with | Some task_id -> Printf.sprintf "%d %s" (Id.Task.to_int task_id) task | None -> Printf.sprintf "0 terminate" @@ -408,14 +409,14 @@ end = struct let create ~state ~client_id ~task_ids = { client_id = Id.Client.of_int client_id ; state = State.of_string state ; - task_ids = List.map Id.Task.of_int task_ids; + task_ids = list_map Id.Task.of_int task_ids; } let to_string x = Printf.sprintf "task_done %s %d %s" (State.to_string x.state) (Id.Client.to_int x.client_id) - (String.concat "|" @@ List.map Id.Task.to_string x.task_ids) + (String.concat "|" @@ list_map Id.Task.to_string x.task_ids) end (** Terminate *) diff --git a/ocaml/Molecule.ml b/ocaml/Molecule.ml index 78ceff0c..186f4d01 100644 --- a/ocaml/Molecule.ml +++ b/ocaml/Molecule.ml @@ -1,4 +1,5 @@ open Qptypes +open Qputils open Sexplib.Std exception MultiplicityError of string @@ -96,7 +97,7 @@ let to_string_general ~f m = let title = name m in - [ string_of_int n ; title ] @ (List.map f nuclei) + [ string_of_int n ; title ] @ (list_map f nuclei) |> String.concat "\n" let to_string = @@ -112,7 +113,7 @@ let of_xyz_string s = let l = String_ext.split s ~on:'\n' |> List.filter (fun x -> x <> "") - |> List.map (fun x -> Atom.of_string units x) + |> list_map (fun x -> Atom.of_string units x) in let ne = ( get_charge { nuclei=l ; @@ -186,7 +187,7 @@ let of_file let distance_matrix molecule = let coord = molecule.nuclei - |> List.map (fun x -> x.Atom.coord) + |> list_map (fun x -> x.Atom.coord) |> Array.of_list in let n = diff --git a/ocaml/Point3d.ml b/ocaml/Point3d.ml index ab086dee..57b02bfe 100644 --- a/ocaml/Point3d.ml +++ b/ocaml/Point3d.ml @@ -1,4 +1,5 @@ open Qptypes +open Qputils open Sexplib.Std type t = { @@ -23,7 +24,7 @@ let of_string ~units s = let l = s |> String_ext.split ~on:' ' |> List.filter (fun x -> x <> "") - |> List.map float_of_string + |> list_map float_of_string |> Array.of_list in { x = l.(0) *. f ; diff --git a/ocaml/Pseudo.ml b/ocaml/Pseudo.ml index 0fd2c263..9bddca02 100644 --- a/ocaml/Pseudo.ml +++ b/ocaml/Pseudo.ml @@ -1,4 +1,5 @@ open Sexplib.Std +open Qputils open Qptypes @@ -81,7 +82,7 @@ let to_string_local = function | t -> "Local component:" :: ( Printf.sprintf "%20s %8s %20s" "Coeff." "r^n" "Exp." ) :: - ( List.map (fun (l,c) -> Printf.sprintf "%20f %8d %20f" + ( list_map (fun (l,c) -> Printf.sprintf "%20f %8d %20f" (AO_coef.to_float c) (R_power.to_int l.GaussianPrimitive_local.r_power) (AO_expo.to_float l.GaussianPrimitive_local.expo) @@ -95,7 +96,7 @@ let to_string_non_local = function | t -> "Non-local component:" :: ( Printf.sprintf "%20s %8s %20s %8s" "Coeff." "r^n" "Exp." "Proj") :: - ( List.map (fun (l,c) -> + ( list_map (fun (l,c) -> let p = Positive_int.to_int l.GaussianPrimitive_non_local.proj in diff --git a/ocaml/Qpackage.ml b/ocaml/Qpackage.ml index 5099a231..d78d3b73 100644 --- a/ocaml/Qpackage.ml +++ b/ocaml/Qpackage.ml @@ -30,7 +30,7 @@ let bit_kind_size = lazy ( in begin match (String_ext.rsplit2 ~on:':' line) with | Some (_ ,buffer) -> - begin match (String_ext.split ~on:'=' buffer |> List.map String.trim) with + begin match (String_ext.split ~on:'=' buffer |> list_map String.trim) with | ["bit_kind_size"; x] -> int_of_string x |> Bit_kind_size.of_int | _ -> get_data tail @@ -58,7 +58,7 @@ let executables = lazy ( result in lines - |> List.map (fun x -> + |> list_map (fun x -> let e = String_ext.split ~on:' ' x |> List.filter (fun x -> x <> "") in diff --git a/ocaml/Qputils.ml b/ocaml/Qputils.ml index 392a6764..270e069f 100644 --- a/ocaml/Qputils.ml +++ b/ocaml/Qputils.ml @@ -53,3 +53,6 @@ let input_lines ic = let string_of_string s = s +let list_map f l = + List.rev_map f l + |> List.rev diff --git a/ocaml/qp_create_ezfio.ml b/ocaml/qp_create_ezfio.ml index 8d09d55b..f4c7f3b6 100644 --- a/ocaml/qp_create_ezfio.ml +++ b/ocaml/qp_create_ezfio.ml @@ -38,7 +38,7 @@ let dummy_centers ~threshold ~molecule ~nuclei = | _ -> assert false in aux [] (n-1,n-1) - |> List.map (fun (i,x,j,y,r) -> + |> list_map (fun (i,x,j,y,r) -> let f = x /. (x +. y) in @@ -270,7 +270,7 @@ let run ?o b au c d m p cart xyz_file = (* Write Pseudo *) let pseudo = - List.map (fun x -> + list_map (fun x -> match pseudo_channel x.Atom.element with | Some channel -> Pseudo.read_element channel x.Atom.element | None -> Pseudo.empty x.Atom.element @@ -292,7 +292,7 @@ let run ?o b au c d m p cart xyz_file = |> Elec_beta_number.of_int; Molecule.nuclei = let charges = - List.map (fun x -> Positive_int.to_int x.Pseudo.n_elec + list_map (fun x -> Positive_int.to_int x.Pseudo.n_elec |> Float.of_int) pseudo |> Array.of_list in @@ -315,13 +315,13 @@ let run ?o b au c d m p cart xyz_file = (* Write Nuclei *) let labels = - List.map (fun x->Element.to_string x.Atom.element) nuclei + list_map (fun x->Element.to_string x.Atom.element) nuclei and charges = - List.map (fun x-> Atom.(Charge.to_float x.charge)) nuclei + list_map (fun x-> Atom.(Charge.to_float x.charge)) nuclei and coords = - (List.map (fun x-> x.Atom.coord.Point3d.x) nuclei) @ - (List.map (fun x-> x.Atom.coord.Point3d.y) nuclei) @ - (List.map (fun x-> x.Atom.coord.Point3d.z) nuclei) in + (list_map (fun x-> x.Atom.coord.Point3d.x) nuclei) @ + (list_map (fun x-> x.Atom.coord.Point3d.y) nuclei) @ + (list_map (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 @@ -365,7 +365,7 @@ let run ?o b au c d m p cart xyz_file = let kmax = Array.init (lmax+1) (fun i-> - List.map (fun x -> + list_map (fun x -> List.filter (fun (y,_) -> (Positive_int.to_int y.Pseudo.GaussianPrimitive_non_local.proj) = i) x.Pseudo.non_local @@ -478,7 +478,7 @@ let run ?o b au c d m p cart xyz_file = in let result = do_work [] 1 nuclei |> List.rev - |> List.map (fun (x,i) -> + |> list_map (fun (x,i) -> try let e = match x.Atom.element with @@ -512,30 +512,30 @@ let run ?o b au c d m p cart xyz_file = 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 (fun (_,g,_) -> List.length g.Gto.lc) long_basis - and ao_nucl = List.map (fun (_,_,n) -> Nucl_number.to_int n) long_basis + let ao_prim_num = list_map (fun (_,g,_) -> List.length g.Gto.lc) long_basis + and ao_nucl = list_map (fun (_,_,n) -> Nucl_number.to_int n) long_basis and ao_power= - let l = List.map (fun (x,_,_) -> x) long_basis in - (List.map (fun t -> Positive_int.to_int Symmetry.Xyz.(t.x)) l)@ - (List.map (fun t -> Positive_int.to_int Symmetry.Xyz.(t.y)) l)@ - (List.map (fun t -> Positive_int.to_int Symmetry.Xyz.(t.z)) l) + let l = list_map (fun (x,_,_) -> x) long_basis in + (list_map (fun t -> Positive_int.to_int Symmetry.Xyz.(t.x)) l)@ + (list_map (fun t -> Positive_int.to_int Symmetry.Xyz.(t.y)) l)@ + (list_map (fun t -> Positive_int.to_int Symmetry.Xyz.(t.z)) l) in let ao_prim_num_max = List.fold_left (fun s x -> if x > s then x else s) 0 ao_prim_num in let gtos = - List.map (fun (_,x,_) -> x) long_basis + list_map (fun (_,x,_) -> x) long_basis in let create_expo_coef ec = let coefs = begin match ec with - | `Coefs -> List.map (fun x-> - List.map (fun (_,coef) -> + | `Coefs -> list_map (fun x-> + list_map (fun (_,coef) -> AO_coef.to_float coef) x.Gto.lc) gtos - | `Expos -> List.map (fun x-> - List.map (fun (prim,_) -> AO_expo.to_float + | `Expos -> list_map (fun x-> + list_map (fun (prim,_) -> AO_expo.to_float prim.GaussianPrimitive.expo) x.Gto.lc) gtos end in diff --git a/src/ao_one_e_ints/ao_ortho_canonical.irp.f b/src/ao_one_e_ints/ao_ortho_canonical.irp.f index 21deed41..45275a06 100644 --- a/src/ao_one_e_ints/ao_ortho_canonical.irp.f +++ b/src/ao_one_e_ints/ao_ortho_canonical.irp.f @@ -79,7 +79,7 @@ BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_inv, (ao_cart_to_sphe_num,ao_ call get_pseudo_inverse(ao_cart_to_sphe_coef,size(ao_cart_to_sphe_coef,1),& ao_num,ao_cart_to_sphe_num, & - ao_cart_to_sphe_inv, size(ao_cart_to_sphe_inv,1)) + ao_cart_to_sphe_inv, size(ao_cart_to_sphe_inv,1), lin_dep_cutoff) END_PROVIDER @@ -107,16 +107,13 @@ END_PROVIDER ao_ortho_canonical_coef(i,i) = 1.d0 enddo -!call ortho_lowdin(ao_overlap,size(ao_overlap,1),ao_num,ao_ortho_canonical_coef,size(ao_ortho_canonical_coef,1),ao_num) -!ao_ortho_canonical_num=ao_num -!return - + call write_double(6, lin_dep_cutoff, "Linear dependencies cut-off") if (ao_cartesian) then ao_ortho_canonical_num = ao_num call ortho_canonical(ao_overlap,size(ao_overlap,1), & ao_num,ao_ortho_canonical_coef,size(ao_ortho_canonical_coef,1), & - ao_ortho_canonical_num) + ao_ortho_canonical_num, lin_dep_cutoff) else @@ -131,7 +128,7 @@ END_PROVIDER ao_ortho_canonical_num = ao_cart_to_sphe_num call ortho_canonical(ao_cart_to_sphe_overlap, size(ao_cart_to_sphe_overlap,1), & - ao_cart_to_sphe_num, S, size(S,1), ao_ortho_canonical_num) + ao_cart_to_sphe_num, S, size(S,1), ao_ortho_canonical_num, lin_dep_cutoff) call dgemm('N','N', ao_num, ao_ortho_canonical_num, ao_cart_to_sphe_num, 1.d0, & ao_cart_to_sphe_coef, size(ao_cart_to_sphe_coef,1), & diff --git a/src/ao_one_e_ints/ao_overlap.irp.f b/src/ao_one_e_ints/ao_overlap.irp.f index d7300936..11c95e42 100644 --- a/src/ao_one_e_ints/ao_overlap.irp.f +++ b/src/ao_one_e_ints/ao_overlap.irp.f @@ -162,7 +162,8 @@ BEGIN_PROVIDER [ double precision, S_inv,(ao_num,ao_num) ] BEGIN_DOC ! Inverse of the overlap matrix END_DOC - call get_pseudo_inverse(ao_overlap,size(ao_overlap,1),ao_num,ao_num,S_inv,size(S_inv,1)) + call get_pseudo_inverse(ao_overlap,size(ao_overlap,1),ao_num,ao_num,S_inv, & + size(S_inv,1),lin_dep_cutoff) END_PROVIDER BEGIN_PROVIDER [ complex*16, S_inv_complex,(ao_num,ao_num) ] @@ -170,8 +171,8 @@ BEGIN_PROVIDER [ complex*16, S_inv_complex,(ao_num,ao_num) ] BEGIN_DOC ! Inverse of the overlap matrix END_DOC - call get_pseudo_inverse_complex(ao_overlap_complex, & - size(ao_overlap_complex,1),ao_num,ao_num,S_inv_complex,size(S_inv_complex,1)) + call get_pseudo_inverse_complex(ao_overlap_complex, size(ao_overlap_complex,1),& + ao_num,ao_num,S_inv_complex,size(S_inv_complex,1),lin_dep_cutoff) END_PROVIDER BEGIN_PROVIDER [ double precision, S_half_inv, (AO_num,AO_num) ] diff --git a/src/mo_guess/mo_ortho_lowdin.irp.f b/src/mo_guess/mo_ortho_lowdin.irp.f index 47a1d24c..8edf7c48 100644 --- a/src/mo_guess/mo_ortho_lowdin.irp.f +++ b/src/mo_guess/mo_ortho_lowdin.irp.f @@ -13,7 +13,7 @@ BEGIN_PROVIDER [double precision, ao_ortho_lowdin_coef, (ao_num,ao_num)] do j=1, ao_num tmp_matrix(j,j) = 1.d0 enddo - call ortho_lowdin(ao_overlap,ao_num,ao_num,tmp_matrix,ao_num,ao_num) + call ortho_lowdin(ao_overlap,ao_num,ao_num,tmp_matrix,ao_num,ao_num,lin_dep_cutoff) do i=1, ao_num do j=1, ao_num ao_ortho_lowdin_coef(j,i) = tmp_matrix(i,j) diff --git a/src/mo_one_e_ints/EZFIO.cfg b/src/mo_one_e_ints/EZFIO.cfg index 0f31b16a..06c91024 100644 --- a/src/mo_one_e_ints/EZFIO.cfg +++ b/src/mo_one_e_ints/EZFIO.cfg @@ -48,3 +48,8 @@ doc: Read/Write |MO| one-electron integrals from/to disk [ Write | Read | None ] interface: ezfio,provider,ocaml default: None +[lin_dep_cutoff] +type: Threshold +doc: Remove linear dependencies when the eigenvalues of the overlap matrix are below this value +interface: ezfio,provider,ocaml +default: 1.e-6 diff --git a/src/mo_one_e_ints/orthonormalize.irp.f b/src/mo_one_e_ints/orthonormalize.irp.f index cffc1993..3a5d5488 100644 --- a/src/mo_one_e_ints/orthonormalize.irp.f +++ b/src/mo_one_e_ints/orthonormalize.irp.f @@ -3,7 +3,7 @@ subroutine orthonormalize_mos integer :: m,p,s m = size(mo_coef,1) p = size(mo_overlap,1) - call ortho_lowdin(mo_overlap,p,mo_num,mo_coef,m,ao_num) + call ortho_lowdin(mo_overlap,p,mo_num,mo_coef,m,ao_num,lin_dep_cutoff) mo_label = 'Orthonormalized' SOFT_TOUCH mo_coef mo_label end diff --git a/src/scf_utils/huckel.irp.f b/src/scf_utils/huckel.irp.f index ac104a72..2d110e69 100644 --- a/src/scf_utils/huckel.irp.f +++ b/src/scf_utils/huckel.irp.f @@ -23,8 +23,6 @@ subroutine huckel_guess Fock_matrix_ao_alpha(1:ao_num,1:ao_num) = A(1:ao_num,1:ao_num) Fock_matrix_ao_beta (1:ao_num,1:ao_num) = A(1:ao_num,1:ao_num) -! TOUCH mo_coef - TOUCH Fock_matrix_ao_alpha Fock_matrix_ao_beta mo_coef = eigenvectors_fock_matrix_mo SOFT_TOUCH mo_coef diff --git a/src/utils/extrapolation.irp.f b/src/utils/extrapolation.irp.f index feb550bb..6af4c0b7 100644 --- a/src/utils/extrapolation.irp.f +++ b/src/utils/extrapolation.irp.f @@ -23,7 +23,7 @@ subroutine extrapolate_data(N_data, data, pt2, output) x(i,2) = pt2_rev(i) enddo do ifit=2,N_data - call get_pseudo_inverse(x,size(x,1),ifit,2,x_inv,size(x_inv,1)) + call get_pseudo_inverse(x,size(x,1),ifit,2,x_inv,size(x_inv,1),1.d-10) ab = matmul(x_inv(1:2,1:ifit),y(1:ifit)) output(ifit) = ab(1) enddo diff --git a/src/utils/linear_algebra.irp.f b/src/utils/linear_algebra.irp.f index 93b367aa..39ffc873 100644 --- a/src/utils/linear_algebra.irp.f +++ b/src/utils/linear_algebra.irp.f @@ -47,14 +47,14 @@ subroutine svd_complex(A,LDA,U,LDU,D,Vt,LDVt,m,n) implicit none BEGIN_DOC ! Compute A = U.D.Vt - ! + ! ! LDx : leftmost dimension of x ! ! Dimension of A is m x n ! A,U,Vt are complex*16 ! D is double precision END_DOC - + integer, intent(in) :: LDA, LDU, LDVt, m, n complex*16, intent(in) :: A(LDA,n) complex*16, intent(out) :: U(LDU,m) @@ -63,12 +63,12 @@ subroutine svd_complex(A,LDA,U,LDU,D,Vt,LDVt,m,n) complex*16,allocatable :: work(:) double precision,allocatable :: rwork(:) integer :: info, lwork, i, j, k, lrwork - + complex*16,allocatable :: A_tmp(:,:) - allocate (A_tmp(LDA,n)) + allocate (A_tmp(LDA,n)) A_tmp = A lrwork = 5*min(m,n) - + ! Find optimal size for temp arrays allocate(work(1),rwork(lrwork)) lwork = -1 @@ -76,25 +76,25 @@ subroutine svd_complex(A,LDA,U,LDU,D,Vt,LDVt,m,n) D, U, LDU, Vt, LDVt, work, lwork, rwork, info) lwork = int(work(1)) deallocate(work) - + allocate(work(lwork)) call zgesvd('A','A', m, n, A_tmp, LDA, & D, U, LDU, Vt, LDVt, work, lwork, rwork, info) deallocate(work,rwork,A_tmp) - + if (info /= 0) then print *, info, ': SVD failed' stop endif - + end - -subroutine ortho_canonical_complex(overlap,LDA,N,C,LDC,m) + +subroutine ortho_canonical_complex(overlap,LDA,N,C,LDC,m,cutoff) implicit none BEGIN_DOC ! Compute C_new=C_old.U.s^-1/2 canonical orthogonalization. ! - ! overlap : overlap matrix + ! overlap : overlap matrix ! ! LDA : leftmost dimension of overlap array ! @@ -108,10 +108,11 @@ subroutine ortho_canonical_complex(overlap,LDA,N,C,LDC,m) ! m : Coefficients matrix is MxN, ( array is (LDC,N) ) ! END_DOC - + integer, intent(in) :: lda, ldc, n integer, intent(out) :: m complex*16, intent(in) :: overlap(lda,n) + double precision, intent(in) :: cutoff complex*16, intent(inout) :: C(ldc,n) complex*16, allocatable :: U(:,:) complex*16, allocatable :: Vt(:,:) @@ -119,19 +120,19 @@ subroutine ortho_canonical_complex(overlap,LDA,N,C,LDC,m) complex*16, allocatable :: S(:,:) !DIR$ ATTRIBUTES ALIGN : 64 :: U, Vt, D integer :: info, i, j - + if (n < 2) then return endif - + allocate (U(ldc,n), Vt(lda,n), D(n), S(lda,n)) - + call svd_complex(overlap,lda,U,ldc,D,Vt,lda,n,n) - + D(:) = dsqrt(D(:)) m=n do i=1,n - if ( D(i) >= 1.d-6 ) then + if ( D(i) >= cutoff ) then D(i) = 1.d0/D(i) else m = i-1 @@ -139,39 +140,39 @@ subroutine ortho_canonical_complex(overlap,LDA,N,C,LDC,m) exit endif enddo - do i=m+1,n + do i=m+1,n D(i) = 0.d0 - enddo - + enddo + do i=1,m if ( D(i) >= 1.d5 ) then print *, 'Warning: Basis set may have linear dependence problems' endif enddo - + do j=1,n do i=1,n S(i,j) = U(i,j)*D(j) enddo enddo - + do j=1,n do i=1,n U(i,j) = C(i,j) enddo enddo - + call zgemm('N','N',n,n,n,(1.d0,0.d0),U,size(U,1),S,size(S,1),(0.d0,0.d0),C,size(C,1)) - deallocate (U, Vt, D, S) - -end - - + deallocate (U, Vt, D, S) + +end + + subroutine ortho_qr_complex(A,LDA,m,n) implicit none BEGIN_DOC ! Orthogonalization using Q.R factorization - ! + ! ! A : matrix to orthogonalize ! ! LDA : leftmost dimension of A @@ -183,7 +184,7 @@ subroutine ortho_qr_complex(A,LDA,m,n) END_DOC integer, intent(in) :: m,n, LDA complex*16, intent(inout) :: A(LDA,n) - + integer :: lwork, info integer, allocatable :: jpvt(:) complex*16, allocatable :: tau(:), work(:) @@ -215,7 +216,7 @@ subroutine ortho_qr_unblocked_complex(A,LDA,m,n) END_DOC integer, intent(in) :: m,n, LDA double precision, intent(inout) :: A(LDA,n) - + integer :: info integer, allocatable :: jpvt(:) double precision, allocatable :: tau(:), work(:) @@ -228,13 +229,13 @@ subroutine ortho_qr_unblocked_complex(A,LDA,m,n) ! call dorg2r(m, n, n, A, LDA, tau, WORK, INFO) ! deallocate(WORK,jpvt,tau) end - -subroutine ortho_lowdin_complex(overlap,LDA,N,C,LDC,m) + +subroutine ortho_lowdin_complex(overlap,LDA,N,C,LDC,m,cutoff) implicit none BEGIN_DOC ! Compute C_new=C_old.S^-1/2 orthogonalization. ! - ! overlap : overlap matrix + ! overlap : overlap matrix ! ! LDA : leftmost dimension of overlap array ! @@ -248,7 +249,7 @@ subroutine ortho_lowdin_complex(overlap,LDA,N,C,LDC,m) ! M : Coefficients matrix is MxN, ( array is (LDC,N) ) ! END_DOC - + integer, intent(in) :: LDA, ldc, n, m complex*16, intent(in) :: overlap(lda,n) complex*16, intent(inout) :: C(ldc,n) @@ -256,8 +257,9 @@ subroutine ortho_lowdin_complex(overlap,LDA,N,C,LDC,m) complex*16, allocatable :: Vt(:,:) double precision, allocatable :: D(:) complex*16, allocatable :: S(:,:) + double precision, intent(in) :: cutoff integer :: info, i, j, k - + if (n < 2) then return endif @@ -267,12 +269,13 @@ subroutine ortho_lowdin_complex(overlap,LDA,N,C,LDC,m) call svd_complex(overlap,lda,U,ldc,D,Vt,lda,n,n) !$OMP PARALLEL DEFAULT(NONE) & - !$OMP SHARED(S,U,D,Vt,n,C,m) & + !$OMP SHARED(S,U,D,Vt,n,C,m,cutoff) & !$OMP PRIVATE(i,j,k) !$OMP DO do i=1,n - if ( D(i) < 1.d-6 ) then + if ( D(i) < cutoff) then + print *, 'Removed Linear dependencies :', 1.d0/D(i) D(i) = 0.d0 else D(i) = 1.d0/dsqrt(D(i)) @@ -294,7 +297,7 @@ subroutine ortho_lowdin_complex(overlap,LDA,N,C,LDC,m) !$OMP END DO NOWAIT endif enddo - + !$OMP BARRIER !$OMP DO do j=1,n @@ -303,11 +306,11 @@ subroutine ortho_lowdin_complex(overlap,LDA,N,C,LDC,m) enddo enddo !$OMP END DO - + !$OMP END PARALLEL call zgemm('N','N',m,n,n,(1.d0,0.d0),U,size(U,1),S,size(S,1),(0.d0,0.d0),C,size(C,1)) - + deallocate(U,Vt,S,D) end @@ -340,15 +343,16 @@ subroutine get_inverse_complex(A,LDA,m,C,LDC) end -subroutine get_pseudo_inverse_complex(A,LDA,m,n,C,LDC) +subroutine get_pseudo_inverse_complex(A,LDA,m,n,C,LDC,cutoff) implicit none BEGIN_DOC ! Find C = A^-1 END_DOC integer, intent(in) :: m,n, LDA, LDC complex*16, intent(in) :: A(LDA,n) + double precision, intent(in) :: cutoff complex*16, intent(out) :: C(LDC,m) - + double precision, allocatable :: D(:), rwork(:) complex*16, allocatable :: U(:,:), Vt(:,:), work(:), A_tmp(:,:) integer :: info, lwork @@ -373,15 +377,15 @@ subroutine get_pseudo_inverse_complex(A,LDA,m,n,C,LDC) print *, info, ':: SVD failed' stop 1 endif - + do i=1,n - if (D(i)/D(1) > 1.d-10) then + if (D(i)/D(1) > cutoff) then D(i) = 1.d0/D(i) else D(i) = 0.d0 endif enddo - + C = (0.d0,0.d0) do i=1,m do j=1,n @@ -390,9 +394,9 @@ subroutine get_pseudo_inverse_complex(A,LDA,m,n,C,LDC) enddo enddo enddo - + deallocate(U,D,Vt,work,A_tmp,rwork) - + end subroutine lapack_diagd_diag_in_place_complex(eigvalues,eigvectors,nmax,n) @@ -475,7 +479,7 @@ subroutine lapack_diagd_diag_in_place_complex(eigvalues,eigvectors,nmax,n) end subroutine lapack_diagd_diag_complex(eigvalues,eigvectors,H,nmax,n) - implicit none + implicit none BEGIN_DOC ! Diagonalize matrix H(complex) ! @@ -617,7 +621,7 @@ subroutine lapack_diagd_complex(eigvalues,eigvectors,H,nmax,n) allocate (work(lwork),iwork(liwork),rwork(lrwork)) call ZHEEVD( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, & rwork, lrwork, iwork, liwork, info ) - deallocate(work,iwork,rwork) + deallocate(work,iwork,rwork) if (info < 0) then @@ -640,7 +644,7 @@ subroutine lapack_diagd_complex(eigvalues,eigvectors,H,nmax,n) end subroutine lapack_diag_complex(eigvalues,eigvectors,H,nmax,n) - implicit none + implicit none BEGIN_DOC ! Diagonalize matrix H (complex) ! @@ -695,10 +699,10 @@ subroutine lapack_diag_complex(eigvalues,eigvectors,H,nmax,n) do j=1,n print *, H(i,j) enddo - enddo + enddo stop 1 end if - + eigvectors = (0.d0,0.d0) eigvalues = 0.d0 do j = 1, n @@ -708,12 +712,12 @@ subroutine lapack_diag_complex(eigvalues,eigvectors,H,nmax,n) enddo enddo deallocate(A,eigenvalues) -end - +end + subroutine matrix_vector_product_complex(u0,u1,matrix,sze,lda) implicit none BEGIN_DOC -! performs u1 += u0 * matrix +! performs u1 += u0 * matrix END_DOC integer, intent(in) :: sze,lda complex*16, intent(in) :: u0(sze) @@ -727,7 +731,7 @@ subroutine matrix_vector_product_complex(u0,u1,matrix,sze,lda) call zhemv('U', sze, (1.d0,0.d0), matrix, lda, u0, incx, (1.d0,0.d0), u1, incy) end -subroutine ortho_canonical(overlap,LDA,N,C,LDC,m) +subroutine ortho_canonical(overlap,LDA,N,C,LDC,m,cutoff) implicit none BEGIN_DOC ! Compute C_new=C_old.U.s^-1/2 canonical orthogonalization. @@ -750,6 +754,7 @@ subroutine ortho_canonical(overlap,LDA,N,C,LDC,m) integer, intent(in) :: lda, ldc, n integer, intent(out) :: m double precision, intent(in) :: overlap(lda,n) + double precision, intent(in) :: cutoff double precision, intent(inout) :: C(ldc,n) double precision, allocatable :: U(:,:) double precision, allocatable :: Vt(:,:) @@ -769,7 +774,7 @@ subroutine ortho_canonical(overlap,LDA,N,C,LDC,m) D(:) = dsqrt(D(:)) m=n do i=1,n - if ( D(i) >= 1.d-6 ) then + if ( D(i) >= cutoff ) then D(i) = 1.d0/D(i) else m = i-1 @@ -840,7 +845,7 @@ subroutine ortho_qr(A,LDA,m,n) call dorgqr(m, n, n, A, LDA, TAU, WORK, LWORK, INFO) ! /!\ int(WORK(1)) becomes negative when WORK(1) > 2147483648 LWORK=max(n,int(WORK(1))) - + deallocate(WORK) allocate(WORK(LWORK)) call dorgqr(m, n, n, A, LDA, TAU, WORK, LWORK, INFO) @@ -874,7 +879,7 @@ subroutine ortho_qr_unblocked(A,LDA,m,n) deallocate(WORK,TAU) end -subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m) +subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m,cutoff) implicit none BEGIN_DOC ! Compute C_new=C_old.S^-1/2 orthogonalization. @@ -896,6 +901,7 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m) integer, intent(in) :: LDA, ldc, n, m double precision, intent(in) :: overlap(lda,n) + double precision, intent(in) :: cutoff double precision, intent(inout) :: C(ldc,n) double precision, allocatable :: U(:,:) double precision, allocatable :: Vt(:,:) @@ -912,12 +918,13 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m) call svd(overlap,lda,U,ldc,D,Vt,lda,n,n) !$OMP PARALLEL DEFAULT(NONE) & - !$OMP SHARED(S,U,D,Vt,n,C,m) & + !$OMP SHARED(S,U,D,Vt,n,C,m,cutoff) & !$OMP PRIVATE(i,j,k) !$OMP DO do i=1,n - if ( D(i) < 1.d-6 ) then + if ( D(i) < cutoff ) then + print *, 'Removed Linear dependencies :', 1.d0/D(i) D(i) = 0.d0 else D(i) = 1.d0/dsqrt(D(i)) @@ -986,13 +993,14 @@ subroutine get_inverse(A,LDA,m,C,LDC) deallocate(ipiv,work) end -subroutine get_pseudo_inverse(A,LDA,m,n,C,LDC) +subroutine get_pseudo_inverse(A,LDA,m,n,C,LDC,cutoff) implicit none BEGIN_DOC ! Find C = A^-1 END_DOC integer, intent(in) :: m,n, LDA, LDC double precision, intent(in) :: A(LDA,n) + double precision, intent(in) :: cutoff double precision, intent(out) :: C(LDC,m) double precision, allocatable :: U(:,:), D(:), Vt(:,:), work(:), A_tmp(:,:) @@ -1020,7 +1028,7 @@ subroutine get_pseudo_inverse(A,LDA,m,n,C,LDC) endif do i=1,n - if (D(i)/D(1) > 1.d-10) then + if (D(i)/D(1) > cutoff) then D(i) = 1.d0/D(i) else D(i) = 0.d0 @@ -1053,7 +1061,7 @@ subroutine find_rotation(A,LDA,B,m,C,n) double precision, allocatable :: A_inv(:,:) allocate(A_inv(LDA,n)) - call get_pseudo_inverse(A,LDA,m,n,A_inv,LDA) + call get_pseudo_inverse(A,LDA,m,n,A_inv,LDA,1.d-10) integer :: i,j,k call dgemm('N','N',n,n,m,1.d0,A_inv,n,B,LDA,0.d0,C,n)