From 7f76f8e2c6b4e587f06223bb80d081960e2037c7 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 20 Oct 2014 12:19:12 +0200 Subject: [PATCH] Added Input_determinants.ml --- ocaml/Bitlist.ml | 12 ++ ocaml/Input.ml | 1 + ocaml/Input_determinants.ml | 262 ++++++++++++++++++++++++++++++++++++ ocaml/Qputils.ml | 1 - ocaml/qptypes_generator.ml | 44 +++++- ocaml/test_input.ml | 9 +- src/Dets/s2.irp.f | 13 +- 7 files changed, 336 insertions(+), 6 deletions(-) create mode 100644 ocaml/Input_determinants.ml diff --git a/ocaml/Bitlist.ml b/ocaml/Bitlist.ml index e5c8e4b7..d197c42b 100644 --- a/ocaml/Bitlist.ml +++ b/ocaml/Bitlist.ml @@ -152,6 +152,15 @@ let xor_operator a b = logical_operator2 Bit.xor_operator a b;; let or_operator a b = logical_operator2 Bit.or_operator a b;; let not_operator b = logical_operator1 Bit.not_operator b ;; + +let popcnt b = + let rec popcnt accu = function + | [] -> accu + | Bit.One::rest -> popcnt (accu+1) rest + | Bit.Zero::rest -> popcnt (accu) rest + in popcnt 0 b +;; + let test_module () = let test = of_int64_list ([-1231L;255L]) in print_string (to_string test); @@ -174,5 +183,8 @@ let test_module () = print_string (to_string (or_operator a b)); print_newline (); print_string (to_string (xor_operator a b)); + print_string (to_string a); + print_int (popcnt a); end ;; + diff --git a/ocaml/Input.ml b/ocaml/Input.ml index ccdc3efa..d434b016 100644 --- a/ocaml/Input.ml +++ b/ocaml/Input.ml @@ -7,5 +7,6 @@ include Input_ao_basis;; include Input_bi_integrals;; include Input_bitmasks;; include Input_cis;; +include Input_determinants;; diff --git a/ocaml/Input_determinants.ml b/ocaml/Input_determinants.ml new file mode 100644 index 00000000..a3f21ca7 --- /dev/null +++ b/ocaml/Input_determinants.ml @@ -0,0 +1,262 @@ +open Qptypes;; +open Qputils;; +open Core.Std;; + +module Determinants : sig + type t + val read : unit -> t + val to_string : t -> string +end = struct + type t = + { n_int : N_int_number.t; + bit_kind : Bit_kind.t; + mo_label : Non_empty_string.t; + n_det : Det_number.t; + n_states : States_number.t; + n_states_diag : States_number.t; + n_det_max_jacobi : Det_number.t; + threshold_generators : Positive_float.t; + threshold_selectors : Positive_float.t; + read_wf : bool; + expected_s2 : Positive_float.t; + s2_eig : bool; + psi_coef : Det_coef.t array; + psi_det : Determinant.t array; + } + ;; + + let get_default = Qpackage.get_ezfio_default "determinants";; + + let read_n_int () = + if not (Ezfio.has_determinants_n_int()) then + Ezfio.get_mo_basis_mo_tot_num () + |> Bitlist.n_int_of_mo_tot_num + |> N_int_number.to_int + |> Ezfio.set_determinants_n_int + ; + Ezfio.get_determinants_n_int () + |> N_int_number.of_int + ;; + + let read_bit_kind () = + if not (Ezfio.has_determinants_bit_kind ()) then + Lazy.force Qpackage.bit_kind + |> Bit_kind.to_int + |> Ezfio.set_determinants_bit_kind + ; + Ezfio.get_determinants_bit_kind () + |> Bit_kind.of_int + ;; + + let read_mo_label () = + if (not (Ezfio.has_determinants_mo_label ())) then + Ezfio.get_mo_basis_mo_label () + |> Ezfio.set_determinants_mo_label + ; + Ezfio.get_determinants_mo_label () + |> Non_empty_string.of_string + ;; + + let read_n_det () = + if not (Ezfio.has_determinants_n_det ()) then + Ezfio.set_determinants_n_det 1 + ; + Ezfio.get_determinants_n_det () + |> Det_number.of_int + ;; + + let read_n_states () = + if not (Ezfio.has_determinants_n_states ()) then + Ezfio.set_determinants_n_states 1 + ; + Ezfio.get_determinants_n_states () + |> States_number.of_int + ;; + + let read_n_states_diag () = + if not (Ezfio.has_determinants_n_states_diag ()) then + read_n_states () + |> States_number.to_int + |> Ezfio.set_determinants_n_states_diag + ; + Ezfio.get_determinants_n_states_diag () + |> States_number.of_int + ;; + + let read_n_det_max_jacobi () = + if not (Ezfio.has_determinants_n_det_max_jacobi ()) then + get_default "n_det_max_jacobi" + |> Int.of_string + |> Ezfio.set_determinants_n_det_max_jacobi + ; + Ezfio.get_determinants_n_det_max_jacobi () + |> Det_number.of_int + ;; + + let read_threshold_generators () = + if not (Ezfio.has_determinants_threshold_generators ()) then + get_default "threshold_generators" + |> Float.of_string + |> Ezfio.set_determinants_threshold_generators + ; + Ezfio.get_determinants_threshold_generators () + |> Positive_float.of_float + ;; + + let read_threshold_selectors () = + if not (Ezfio.has_determinants_threshold_selectors ()) then + get_default "threshold_selectors" + |> Float.of_string + |> Ezfio.set_determinants_threshold_selectors + ; + Ezfio.get_determinants_threshold_selectors () + |> Positive_float.of_float + ;; + + let read_read_wf () = + if not (Ezfio.has_determinants_read_wf ()) then + get_default "read_wf" + |> Bool.of_string + |> Ezfio.set_determinants_read_wf + ; + Ezfio.get_determinants_read_wf () + ;; + + let read_expected_s2 () = + if not (Ezfio.has_determinants_expected_s2 ()) then + begin + let na = Ezfio.get_electrons_elec_alpha_num () + and nb = Ezfio.get_electrons_elec_beta_num () + in + let s = 0.5 *. (Float.of_int (na - nb)) + in + Ezfio.set_determinants_expected_s2 ( s *. (s +. 1.) ) + end + ; + Ezfio.get_determinants_expected_s2 () + |> Positive_float.of_float + ;; + + let read_s2_eig () = + if not (Ezfio.has_determinants_s2_eig ()) then + get_default "s2_eig" + |> Bool.of_string + |> Ezfio.set_determinants_s2_eig + ; + Ezfio.get_determinants_s2_eig () + ;; + + let read_psi_coef () = + if not (Ezfio.has_determinants_psi_coef ()) then + Ezfio.ezfio_array_of_list ~rank:1 ~dim:[| 1 |] ~data:[1.] + |> Ezfio.set_determinants_psi_coef + ; + (Ezfio.get_determinants_psi_coef ()).Ezfio.data + |> Ezfio.flattened_ezfio_data + |> Array.map ~f:Det_coef.of_float + ;; + + let read_psi_det () = + let n_int = read_n_int () in + if not (Ezfio.has_determinants_psi_det ()) then + begin + let rec build_data accu = function + | 0 -> accu + | n -> build_data ((MO_number.of_int n)::accu) (n-1) + in + let det_a = build_data [] (Ezfio.get_electrons_elec_alpha_num ()) + |> Bitlist.of_mo_number_list n_int + and det_b = build_data [] (Ezfio.get_electrons_elec_beta_num ()) + |> Bitlist.of_mo_number_list n_int + in + let data = ( (Bitlist.to_int64_list det_a) @ + (Bitlist.to_int64_list det_b) ) + in + Ezfio.ezfio_array_of_list ~rank:3 ~dim:[| N_int_number.to_int n_int ; 2 ; 1 |] ~data:data + |> Ezfio.set_determinants_psi_det + end ; + let n_int = N_int_number.to_int n_int in + let rec transform accu1 accu2 n_rest = function + | [] -> + let accu1 = List.rev accu1 + |> Array.of_list + |> Determinant.of_int64_array + in + List.rev (accu1::accu2) |> Array.of_list + | i::rest -> + if (n_rest > 0) then + transform (i::accu1) accu2 (n_rest-1) rest + else + let accu1 = List.rev accu1 + |> Array.of_list + |> Determinant.of_int64_array + in + transform [] (accu1::accu2) (2*n_int) rest + in + (Ezfio.get_determinants_psi_det ()).Ezfio.data + |> Ezfio.flattened_ezfio_data + |> Array.to_list + |> transform [] [] (2*n_int) + ;; + + let read () = + { n_int = read_n_int () ; + bit_kind = read_bit_kind () ; + mo_label = read_mo_label () ; + n_det = read_n_det () ; + n_states = read_n_states () ; + n_states_diag = read_n_states_diag () ; + n_det_max_jacobi = read_n_det_max_jacobi () ; + threshold_generators = read_threshold_generators () ; + threshold_selectors = read_threshold_selectors () ; + read_wf = read_read_wf () ; + expected_s2 = read_expected_s2 () ; + s2_eig = read_s2_eig () ; + psi_coef = read_psi_coef () ; + psi_det = read_psi_det () ; + } + ;; + + let to_string b = + Printf.sprintf " +n_int = %s +bit_kind = %s +mo_label = \"%s\" +n_det = %s +n_states = %s +n_states_diag = %s +n_det_max_jacobi = %s +threshold_generators = %s +threshold_selectors = %s +read_wf = %s +expected_s2 = %s +s2_eig = %s +psi_coef = %s +psi_det = %s +" + (b.n_int |> N_int_number.to_string) + (b.bit_kind |> Bit_kind.to_string) + (b.mo_label |> Non_empty_string.to_string) + (b.n_det |> Det_number.to_string) + (b.n_states |> States_number.to_string) + (b.n_states_diag |> States_number.to_string) + (b.n_det_max_jacobi |> Det_number.to_string) + (b.threshold_generators |> Positive_float.to_string) + (b.threshold_selectors |> Positive_float.to_string) + (b.read_wf |> Bool.to_string) + (b.expected_s2 |> Positive_float.to_string) + (b.s2_eig |> Bool.to_string) + (b.psi_coef |> Array.to_list |> List.map ~f:Det_coef.to_string + |> String.concat ~sep:", ") + (b.psi_det |> Array.map ~f:(fun x -> Determinant.to_int64_array x + |> Array.map ~f:(fun x-> + print_endline (Int64.to_string x) ; + Int64.to_string x )|> Array.to_list |> + String.concat ~sep:", ") |> Array.to_list + |> String.concat ~sep:" | ") + ; +;; + +end + + diff --git a/ocaml/Qputils.ml b/ocaml/Qputils.ml index 305bb0ef..fd40547f 100644 --- a/ocaml/Qputils.ml +++ b/ocaml/Qputils.ml @@ -9,4 +9,3 @@ let rec transpose = function ;; - diff --git a/ocaml/qptypes_generator.ml b/ocaml/qptypes_generator.ml index 7a0ed9c1..24b08b25 100644 --- a/ocaml/qptypes_generator.ml +++ b/ocaml/qptypes_generator.ml @@ -22,6 +22,14 @@ let input_data = " * Negative_int : int assert (x <= 0) ; +* Det_coef : float + assert (x >= -1.) ; + assert (x <= 1.) ; + +* Normalized_float : float + assert (x <= 1.) ; + assert (x >= 0.) ; + * Strictly_negative_int : int assert (x < 0) ; @@ -70,6 +78,15 @@ let input_data = " if (Ezfio.has_determinants_det_num ()) then assert (x <= (Ezfio.get_determinants_det_num ())); +* States_number : int + assert (x > 0) ; + if (x > 100) then + warning \"More than 100 states\"; + if (Ezfio.has_determinants_n_states_diag ()) then + assert (x <= (Ezfio.get_determinants_n_states_diag ())) + else if (Ezfio.has_determinants_n_states ()) then + assert (x <= (Ezfio.get_determinants_n_states ())); + * Bit_kind_size : int begin match x with | 8 | 16 | 32 | 64 -> () @@ -89,6 +106,28 @@ let input_data = " " ;; +let untouched = " +module Determinant : sig + type t + val to_int64_array : t -> int64 array + val of_int64_array : int64 array -> t + val to_string : t -> string +end = struct + type t = int64 array + let to_int64_array x = x + let of_int64_array x = + if (Ezfio.has_determinants_n_int ()) then + begin + let n_int = Ezfio.get_determinants_n_int () in + assert ((Array.length x) = n_int*2) + end + ; x + let to_string x = Array.to_list x + |> List.map ~f:Int64.to_string + |> String.concat ~sep:\", \" +end + +" let template = format_of_string " module %s : sig @@ -129,6 +168,9 @@ let parse_input input= |> print_string ;; -let () = parse_input input_data;; +let () = + parse_input input_data ; + print_endline untouched +;; diff --git a/ocaml/test_input.ml b/ocaml/test_input.ml index 6deefb40..3e6fd458 100644 --- a/ocaml/test_input.ml +++ b/ocaml/test_input.ml @@ -26,4 +26,11 @@ let test_cis () = print_endline (Input.Cis_dressed.to_string b); ;; -test_cis ();; +let test_dets () = + Ezfio.set_file "F2.ezfio" ; + let b = Input.Determinants.read () + in + print_endline (Input.Determinants.to_string b); +;; + +test_dets();; diff --git a/src/Dets/s2.irp.f b/src/Dets/s2.irp.f index 18c3dfaf..588ab4cc 100644 --- a/src/Dets/s2.irp.f +++ b/src/Dets/s2.irp.f @@ -36,6 +36,9 @@ end BEGIN_PROVIDER [ double precision, S_z ] &BEGIN_PROVIDER [ double precision, S_z2_Sz ] implicit none + BEGIN_DOC +! z component of the Spin + END_DOC S_z = 0.5d0*dble(elec_alpha_num-elec_beta_num) S_z2_Sz = S_z*(S_z-1.d0) @@ -44,15 +47,19 @@ END_PROVIDER BEGIN_PROVIDER [ double precision, expected_s2] implicit none - PROVIDE ezfio_filename + BEGIN_DOC +! Expected value of S2 : S*(S+1) + END_DOC logical :: has_expected_s2 call ezfio_has_determinants_expected_s2(has_expected_s2) if (has_expected_s2) then call ezfio_get_determinants_expected_s2(expected_s2) else - expected_s2 = elec_alpha_num - elec_beta_num + 0.5d0 * ((elec_alpha_num - elec_beta_num)**2*0.5d0 - (elec_alpha_num-elec_beta_num)) -! call ezfio_set_determinants_expected_s2(expected_s2) + double precision :: S + S = (elec_alpha_num-elec_beta_num)*0.5d0 + expected_s2 = S * (S+1.d0) +! expected_s2 = elec_alpha_num - elec_beta_num + 0.5d0 * ((elec_alpha_num - elec_beta_num)**2*0.5d0 - (elec_alpha_num-elec_beta_num)) endif END_PROVIDER