From 2da264acd113ce93cdbef58b7973c5dfdfb8b910 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 13 Sep 2017 20:41:12 +0200 Subject: [PATCH] Fixed OCaml --- ocaml/Input_ao_basis.ml | 387 +++++++++++++++++++++++----------- src/AO_Basis/EZFIO.cfg | 7 +- src/AO_Basis/ao_overlap.irp.f | 2 +- 3 files changed, 263 insertions(+), 133 deletions(-) diff --git a/ocaml/Input_ao_basis.ml b/ocaml/Input_ao_basis.ml index e01ca713..348ddd7f 100644 --- a/ocaml/Input_ao_basis.ml +++ b/ocaml/Input_ao_basis.ml @@ -1,38 +1,101 @@ -(* =~=~ *) -(* Init *) -(* =~=~ *) - open Qptypes;; open Qputils;; open Core;; module Ao_basis : sig -(* Generate type *) - type t = - { - threshold_overlap_ao_eigenvalues : Threshold.t; - } [@@deriving sexp] - ;; - val read : unit -> t option - val write : t-> unit + type t = + { ao_basis : AO_basis_name.t; + ao_num : AO_number.t ; + ao_prim_num : AO_prim_number.t array; + ao_prim_num_max : AO_prim_number.t; + ao_nucl : Nucl_number.t array; + ao_power : Symmetry.Xyz.t array; + ao_coef : AO_coef.t array; + ao_expo : AO_expo.t array; + ao_cartesian : bool; + } [@@deriving sexp] + ;; + val read : unit -> t option val to_string : t -> string + val to_basis : t -> Basis.t + val write : t -> unit + val to_md5 : t -> MD5.t val to_rst : t -> Rst_string.t - val of_rst : Rst_string.t -> t option end = struct -(* Generate type *) - type t = - { - threshold_overlap_ao_eigenvalues : Threshold.t; - } [@@deriving sexp] - ;; + type t = + { ao_basis : AO_basis_name.t; + ao_num : AO_number.t ; + ao_prim_num : AO_prim_number.t array; + ao_prim_num_max : AO_prim_number.t; + ao_nucl : Nucl_number.t array; + ao_power : Symmetry.Xyz.t array; + ao_coef : AO_coef.t array; + ao_expo : AO_expo.t array; + ao_cartesian : bool; + } [@@deriving sexp] + ;; let get_default = Qpackage.get_ezfio_default "ao_basis";; -(* =~=~=~=~=~=~==~=~=~=~=~=~ *) -(* Generate Special Function *) -(* =~=~=~==~=~~=~=~=~=~=~=~=~ *) + let read_ao_basis () = + Ezfio.get_ao_basis_ao_basis () + |> AO_basis_name.of_string + ;; + + let read_ao_num () = + Ezfio.get_ao_basis_ao_num () + |> AO_number.of_int + ;; + + let read_ao_prim_num () = + Ezfio.get_ao_basis_ao_prim_num () + |> Ezfio.flattened_ezfio + |> Array.map ~f:AO_prim_number.of_int + ;; + + let read_ao_prim_num_max () = + Ezfio.get_ao_basis_ao_prim_num () + |> Ezfio.flattened_ezfio + |> Array.fold ~f:(fun x y -> if x>y then x else y) ~init:0 + |> AO_prim_number.of_int + ;; + + let read_ao_nucl () = + let nmax = Nucl_number.get_max () in + Ezfio.get_ao_basis_ao_nucl () + |> Ezfio.flattened_ezfio + |> Array.map ~f:(fun x-> Nucl_number.of_int ~max:nmax x) + ;; + + let read_ao_power () = + let x = Ezfio.get_ao_basis_ao_power () in + let dim = x.Ezfio.dim.(0) in + let data = Ezfio.flattened_ezfio x in + let result = Array.init dim ~f:(fun x -> "") in + for i=1 to dim + do + if (data.(i-1) > 0) then + result.(i-1) <- result.(i-1)^"x"^(Int.to_string data.(i-1)); + if (data.(dim+i-1) > 0) then + result.(i-1) <- result.(i-1)^"y"^(Int.to_string data.(dim+i-1)); + if (data.(2*dim+i-1) > 0) then + result.(i-1) <- result.(i-1)^"z"^(Int.to_string data.(2*dim+i-1)); + done; + Array.map ~f:Symmetry.Xyz.of_string result + ;; + + let read_ao_coef () = + Ezfio.get_ao_basis_ao_coef () + |> Ezfio.flattened_ezfio + |> Array.map ~f:AO_coef.of_float + ;; + + let read_ao_expo () = + Ezfio.get_ao_basis_ao_expo () + |> Ezfio.flattened_ezfio + |> Array.map ~f:AO_expo.of_float + ;; -(* Read snippet for ao_cartesian *) let read_ao_cartesian () = if not (Ezfio.has_ao_basis_ao_cartesian ()) then get_default "ao_cartesian" @@ -41,119 +104,191 @@ end = struct ; Ezfio.get_ao_basis_ao_cartesian () ;; -(* Write snippet for ao_cartesian *) - let write_ao_cartesian = - Ezfio.set_ao_basis_ao_cartesian + + let to_long_basis b = + let ao_num = AO_number.to_int b.ao_num in + let gto_array = Array.init (AO_number.to_int b.ao_num) + ~f:(fun i -> + let s = Symmetry.Xyz.to_symmetry b.ao_power.(i) in + let ao_prim_num = AO_prim_number.to_int b.ao_prim_num.(i) in + let prims = List.init ao_prim_num ~f:(fun j -> + let prim = { GaussianPrimitive.sym = s ; + GaussianPrimitive.expo = b.ao_expo.(ao_num*j+i) + } + in + let coef = b.ao_coef.(ao_num*j+i) in + (prim,coef) + ) in + Gto.of_prim_coef_list prims + ) + in + let rec do_work accu sym gto nucl = + match (sym, gto, nucl) with + | (s::srest, g::grest, n::nrest) -> + do_work ((s,g,n)::accu) srest grest nrest + | ([],[],[]) -> List.rev accu + | _ -> assert false + in + do_work [] + (Array.to_list b.ao_power) + (Array.to_list gto_array) + (Array.to_list b.ao_nucl) + ;; + let to_basis b = + to_long_basis b + |> Long_basis.to_basis + ;; + + let to_md5 b = + let short_basis = to_basis b in + Basis.to_md5 short_basis + ;; + + + + let write_md5 b = + to_md5 b + |> MD5.to_string + |> Ezfio.set_ao_basis_ao_md5 ;; -(* Read snippet for ao_prim_num_max *) - let read_ao_prim_num_max () = - if not (Ezfio.has_ao_basis_ao_prim_num_max ()) then - get_default "ao_prim_num_max" - |> Int.of_string - |> Ezfio.set_ao_basis_ao_prim_num_max - ; - Ezfio.get_ao_basis_ao_prim_num_max () - ;; -(* Write snippet for ao_prim_num_max *) - let write_ao_prim_num_max = - Ezfio.set_ao_basis_ao_prim_num_max + let write_ao_basis name = + AO_basis_name.to_string name + |> Ezfio.set_ao_basis_ao_basis ;; -(* Read snippet for integral_kinetic *) - let read_integral_kinetic () = - if not (Ezfio.has_ao_basis_integral_kinetic ()) then - get_default "integral_kinetic" - |> Float.of_string - |> Ezfio.set_ao_basis_integral_kinetic - ; - Ezfio.get_ao_basis_integral_kinetic () - ;; -(* Write snippet for integral_kinetic *) - let write_integral_kinetic = - Ezfio.set_ao_basis_integral_kinetic + let write b = + let { ao_basis ; + ao_num ; + ao_prim_num ; + ao_prim_num_max ; + ao_nucl ; + ao_power ; + ao_coef ; + ao_expo ; + ao_cartesian ; + } = b + in + write_md5 b ; + write_ao_basis ao_basis; ;; -(* Read snippet for integral_nuclear *) - let read_integral_nuclear () = - if not (Ezfio.has_ao_basis_integral_nuclear ()) then - get_default "integral_nuclear" - |> Float.of_string - |> Ezfio.set_ao_basis_integral_nuclear - ; - Ezfio.get_ao_basis_integral_nuclear () + + let read () = + if (Ezfio.has_ao_basis_ao_basis ()) then + begin + let result = + { ao_basis = read_ao_basis (); + ao_num = read_ao_num () ; + ao_prim_num = read_ao_prim_num (); + ao_prim_num_max = read_ao_prim_num_max (); + ao_nucl = read_ao_nucl (); + ao_power = read_ao_power (); + ao_coef = read_ao_coef () ; + ao_expo = read_ao_expo () ; + ao_cartesian = read_ao_cartesian () ; + } + in + to_md5 result + |> MD5.to_string + |> Ezfio.set_ao_basis_ao_md5 ; + Some result + end + else + None ;; -(* Write snippet for integral_nuclear *) - let write_integral_nuclear = - Ezfio.set_ao_basis_integral_nuclear + + + let to_rst b = + let print_sym = + let l = List.init (Array.length b.ao_power) ~f:( + fun i -> ( (i+1),b.ao_nucl.(i),b.ao_power.(i) ) ) in + let rec do_work = function + | [] -> [] + | (i,n,x)::tail -> + (Printf.sprintf " %5d %6d %-8s\n" i (Nucl_number.to_int n) (Symmetry.Xyz.to_string x)):: + (do_work tail) + in do_work l + |> String.concat + in + + let short_basis = to_basis b in + Printf.sprintf " +Name of the AO basis :: + + ao_basis = %s + +Cartesian coordinates (6d,10f,...) :: + + ao_cartesian = %s + +Basis set (read-only) :: + +%s + + +======= ========= =========== + Basis Nucleus Symmetries +======= ========= =========== +%s +======= ========= =========== + +" (AO_basis_name.to_string b.ao_basis) + (Bool.to_string b.ao_cartesian) + (Basis.to_string short_basis + |> String.split ~on:'\n' + |> List.map ~f:(fun x-> " "^x) + |> String.concat ~sep:"\n" + ) print_sym + + |> Rst_string.of_string ;; -(* Read snippet for integral_overlap *) - let read_integral_overlap () = - if not (Ezfio.has_ao_basis_integral_overlap ()) then - get_default "integral_overlap" - |> Float.of_string - |> Ezfio.set_ao_basis_integral_overlap - ; - Ezfio.get_ao_basis_integral_overlap () - ;; -(* Write snippet for integral_overlap *) - let write_integral_overlap = - Ezfio.set_ao_basis_integral_overlap + let read_rst s = + let s = Rst_string.to_string s + |> String.split ~on:'\n' + in + let rec extract_basis = function + | [] -> failwith "Error in basis set" + | line :: tail -> + let line = String.strip line in + if line = "Basis set (read-only) ::" then + String.concat tail ~sep:"\n" + else + extract_basis tail + in + extract_basis s ;; -(* Read snippet for threshold_overlap_ao_eigenvalues *) - let read_threshold_overlap_ao_eigenvalues () = - if not (Ezfio.has_ao_basis_threshold_overlap_ao_eigenvalues ()) then - get_default "threshold_overlap_ao_eigenvalues" - |> Float.of_string - |> Ezfio.set_ao_basis_threshold_overlap_ao_eigenvalues - ; - Ezfio.get_ao_basis_threshold_overlap_ao_eigenvalues () - |> Threshold.of_float - ;; -(* Write snippet for threshold_overlap_ao_eigenvalues *) - let write_threshold_overlap_ao_eigenvalues var = - Threshold.to_float var - |> Ezfio.set_ao_basis_threshold_overlap_ao_eigenvalues + let to_string b = + Printf.sprintf " +ao_basis = %s +ao_num = %s +ao_prim_num = %s +ao_prim_num_max = %s +ao_nucl = %s +ao_power = %s +ao_coef = %s +ao_expo = %s +ao_cartesian = %s +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 + ~f:(AO_prim_number.to_string) |> String.concat ~sep:", " ) + (AO_prim_number.to_string b.ao_prim_num_max) + (b.ao_nucl |> Array.to_list |> List.map ~f:Nucl_number.to_string |> + String.concat ~sep:", ") + (b.ao_power |> Array.to_list |> List.map ~f:(fun x-> + "("^(Symmetry.Xyz.to_string x)^")" )|> String.concat ~sep:", ") + (b.ao_coef |> Array.to_list |> List.map ~f:AO_coef.to_string + |> String.concat ~sep:", ") + (b.ao_expo |> Array.to_list |> List.map ~f:AO_expo.to_string + |> String.concat ~sep:", ") + (b.ao_cartesian |> Bool.to_string) + (to_md5 b |> MD5.to_string ) + ;; +end -(* =~=~=~=~=~=~=~=~=~=~=~=~ *) -(* Generate Global Function *) -(* =~=~=~=~=~=~=~=~=~=~=~=~ *) - -(* Read all *) - let read() = - Some - { - threshold_overlap_ao_eigenvalues = read_threshold_overlap_ao_eigenvalues (); - } - ;; -(* Write all *) - let write{ - threshold_overlap_ao_eigenvalues; - } = - write_threshold_overlap_ao_eigenvalues threshold_overlap_ao_eigenvalues; - ;; -(* to_string*) - let to_string b = - Printf.sprintf " - threshold_overlap_ao_eigenvalues = %s - " - (Threshold.to_string b.threshold_overlap_ao_eigenvalues) - ;; -(* to_rst*) - let to_rst b = - Printf.sprintf " - Threshold on the magnitude of the smallest eigenvalues of the overlap matrix in the AO basis :: - - threshold_overlap_ao_eigenvalues = %s - - " - (Threshold.to_string b.threshold_overlap_ao_eigenvalues) - |> Rst_string.of_string - ;; - include Generic_input_of_rst;; - let of_rst = of_rst t_of_sexp;; - -end \ No newline at end of file diff --git a/src/AO_Basis/EZFIO.cfg b/src/AO_Basis/EZFIO.cfg index e02f791e..07712a4a 100644 --- a/src/AO_Basis/EZFIO.cfg +++ b/src/AO_Basis/EZFIO.cfg @@ -16,7 +16,7 @@ interface: ezfio, provider [ao_prim_num_max] type: integer -doc: number of primitive maximun +doc: maximum number of primitives default: =maxval(ao_basis.ao_prim_num) interface: ezfio @@ -76,9 +76,4 @@ size: (ao_basis.ao_num,ao_basis.ao_num) interface: ezfio default: false -[threshold_overlap_ao_eigenvalues] -type: Threshold -doc: Threshold on the magnitude of the smallest eigenvalues of the overlap matrix in the AO basis -interface: ezfio,provider,ocaml -default: 1.e-6 diff --git a/src/AO_Basis/ao_overlap.irp.f b/src/AO_Basis/ao_overlap.irp.f index 1d3b24e6..fd3fa8ec 100644 --- a/src/AO_Basis/ao_overlap.irp.f +++ b/src/AO_Basis/ao_overlap.irp.f @@ -143,6 +143,7 @@ BEGIN_PROVIDER [ double precision, S_half_inv, (AO_num,AO_num) ] integer :: LDA, LDC double precision, allocatable :: U(:,:),Vt(:,:), D(:) integer :: info, i, j, k + double precision, parameter :: threshold_overlap_AO_eigenvalues = 1.d-6 LDA = size(AO_overlap,1) LDC = size(S_half_inv,1) @@ -174,7 +175,6 @@ BEGIN_PROVIDER [ double precision, S_half_inv, (AO_num,AO_num) ] enddo enddo write(*,*) 'linear dependencies',num_linear_dependencies -! stop do k=1,AO_num if(D(k) /= 0.d0) then