From 23389acfc88531b175e56d22f9fa8f5de00a6eeb Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 22 Feb 2018 18:20:45 +0100 Subject: [PATCH] Hartree-Fock OK --- Basis/TwoElectronRRVectorized.ml | 231 +++++++++++++++---------------- Simulation.ml | 19 ++- Utils/Util.ml | 6 + run_integrals.ml | 2 +- 4 files changed, 134 insertions(+), 124 deletions(-) diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml index d66dfb8..fed0ff3 100644 --- a/Basis/TwoElectronRRVectorized.ml +++ b/Basis/TwoElectronRRVectorized.ml @@ -329,7 +329,6 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) result -(* and trr_v angMom_a angMom_c = match (angMom_a.tot, angMom_c.tot) with @@ -340,123 +339,120 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) try Zmap.find map_2d.(0) key with | Not_found -> - let result = - let xyz = get_xyz angMom_c in - let axyz = Powers.get xyz angMom_a in - let cm = Powers.decr xyz angMom_c in - let cmxyz = Powers.get xyz cm in - let expo_inv_q_over_p = - Array.mapi (fun l expo_inv_p_l -> - let expo_p_l = 1./.expo_inv_p_l in - Array.mapi (fun k expo_inv_q_k -> - expo_inv_q_k *. expo_p_l) expo_inv_q ) expo_inv_p - in - let result = None in + let xyz = get_xyz angMom_c in + let axyz = Powers.get xyz angMom_a in + let cm = Powers.decr xyz angMom_c in + let cmxyz = Powers.get xyz cm in + let expo_inv_q_over_p = + Array.mapi (fun l expo_inv_p_l -> + let expo_p_l = 1./.expo_inv_p_l in + Array.mapi (fun k expo_inv_q_k -> + expo_inv_q_k *. expo_p_l) expo_inv_q ) expo_inv_p + in + let result = None in - let result = - if cmxyz < 1 then result else - let f = 0.5 *. (float_of_int cmxyz) in - if abs_float f < cutoff then result else - begin - let cmm = Powers.decr xyz cm in - match result, trr_v angMom_a cmm with - | None, None -> None - | Some result, None -> Some result - | None, Some v3 -> - Some (Array.init np (fun l -> - let v3_l = v3.(l) in - Array.mapi (fun k v3_lk -> - expo_inv_q.(k) *. f *. v3_lk) v3_l - ) ) - | Some result, Some v3 -> - (Array.iteri (fun l v3_l -> - let result_l = result.(l) in - Array.iteri (fun k v3_lk -> - result_l.(k) <- result_l.(k) +. - expo_inv_q.(k) *. f *. v3_lk) v3_l - ) v3 ; Some result) - end - in - let result = - match result, trr_v angMom_a cm with + let result = + if cmxyz < 1 then result else + begin + let f = 0.5 *. (float_of_int cmxyz) in + let cmm = Powers.decr xyz cm in + match result, trr_v angMom_a cmm with | None, None -> None + | None, Some v3 -> + Some (Array.init np (fun l -> + let v3_l = v3.(l) in + Array.mapi (fun k v3_lk -> + expo_inv_q.(k) *. f *. v3_lk) v3_l + ) ) | Some result, None -> Some result - | None, Some v1 -> - Some (Array.init np (fun l -> - let v1_l = v1.(l) - and cpa = (center_pa xyz).(l) - and expo_inv_q_over_p_l = expo_inv_q_over_p.(l) - in - Array.mapi (fun k v1_lk -> - let cqc = (center_qc xyz).(k) in - (cqc +. expo_inv_q_over_p_l.(k) *. cpa) *. v1_lk - ) v1_l - ) ) - | Some result, Some v1 -> - (Array.iteri (fun l v1_l -> - let cpa = (center_pa xyz).(l) - and result_l = result.(l) - and expo_inv_q_over_p_l = expo_inv_q_over_p.(l) - in - Array.iteri (fun k v1_lk -> - let cqc = (center_qc xyz).(k) in - result_l.(k) <- result_l.(k) +. - (cqc +. expo_inv_q_over_p_l.(k) *. cpa) *. v1_lk - ) v1_l - ) v1 ; Some result) - in - let result = - if cmxyz < 0 then result else - begin - let ap = Powers.incr xyz angMom_a in - match result, trr_v ap cm with - | None, None -> None - | Some result, None -> Some result - | None, Some v4 -> - Some (Array.init np (fun l -> - let v4_l = v4.(l) in - let expo_inv_q_over_p_l = expo_inv_q_over_p.(l) in - Array.mapi (fun k v4_lk -> - -. expo_inv_q_over_p_l.(k) *. v4_lk) v4_l - ) ) - | Some result, Some v4 -> - (Array.iteri (fun l v4_l -> - let result_l = result.(l) in - Array.iteri (fun k v4_lk -> - let expo_inv_q_over_p_l = expo_inv_q_over_p.(l) in - result_l.(k) <- result_l.(k) - -. expo_inv_q_over_p_l.(k) *. v4_lk) v4_l - ) v4 ; Some result) - end - in - let result = - if axyz < 1 then result else - let f = 0.5 *. (float_of_int axyz) in - if abs_float f < cutoff then result else - begin - let am = Powers.decr xyz angMom_a in - match result, trr_v am cm with - | None, None -> None - | Some result, None -> Some result - | None, Some v2 -> - Some (Array.init np (fun l -> - let v2_l = v2.(l) in - Array.mapi (fun k v2_lk -> - expo_inv_q.(k) *. f *. v2_lk) v2_l - ) ) - | Some result, Some v2 -> - (Array.iteri (fun l v2_l -> - let result_l = result.(l) in - Array.iteri (fun k v2_lk -> - result_l.(k) <- result_l.(k) +. - expo_inv_q.(k) *. f *. v2_lk) v2_l - ) v2; Some result) - end - in result - in - Zmap.add map_2d.(0) key result; - result -*) + | Some result, Some v3 -> + (Array.iteri (fun l v3_l -> + let result_l = result.(l) in + Array.iteri (fun k v3_lk -> + result_l.(k) <- result_l.(k) +. + expo_inv_q.(k) *. f *. v3_lk) v3_l + ) v3 ; Some result) + end + in + let result = + begin + match result, trr_v angMom_a cm with + | Some result, None -> Some result + | Some result, Some v1 -> + (Array.iteri (fun l v1_l -> + let cpa = (center_pa xyz).(l) + and result_l = result.(l) + and expo_inv_q_over_p_l = expo_inv_q_over_p.(l) + in + Array.iteri (fun k v1_lk -> + let cqc = (center_qc xyz).(k) in + result_l.(k) <- result_l.(k) +. + (cqc +. expo_inv_q_over_p_l.(k) *. cpa) *. v1_lk + ) v1_l + ) v1 ; Some result) + | None, None -> None + | None, Some v1 -> + Some (Array.init np (fun l -> + let v1_l = v1.(l) + and cpa = (center_pa xyz).(l) + and expo_inv_q_over_p_l = expo_inv_q_over_p.(l) + in + Array.mapi (fun k v1_lk -> + let cqc = (center_qc xyz).(k) in + (cqc +. expo_inv_q_over_p_l.(k) *. cpa) *. v1_lk + ) v1_l + ) ) + end + in + let result = + if cmxyz < 0 then result else + begin + let ap = Powers.incr xyz angMom_a in + match result, trr_v ap cm with + | Some result, None -> Some result + | Some result, Some v4 -> + (Array.iteri (fun l v4_l -> + let result_l = result.(l) in + Array.iteri (fun k v4_lk -> + let expo_inv_q_over_p_l = expo_inv_q_over_p.(l) in + result_l.(k) <- result_l.(k) + -. expo_inv_q_over_p_l.(k) *. v4_lk) v4_l + ) v4 ; Some result) + | None, None -> None + | None, Some v4 -> + Some (Array.init np (fun l -> + let v4_l = v4.(l) in + let expo_inv_q_over_p_l = expo_inv_q_over_p.(l) in + Array.mapi (fun k v4_lk -> + -. expo_inv_q_over_p_l.(k) *. v4_lk) v4_l + ) ) + end + in + let result = + if axyz < 1 then result else + begin + let f = 0.5 *. (float_of_int axyz) in + let am = Powers.decr xyz angMom_a in + match result, trr_v am cm with + | Some result, None -> Some result + | Some result, Some v2 -> + (Array.iteri (fun l v2_l -> + let result_l = result.(l) in + Array.iteri (fun k v2_lk -> + result_l.(k) <- result_l.(k) +. + expo_inv_q.(k) *. f *. v2_lk) v2_l + ) v2; Some result) + | None, None -> None + | None, Some v2 -> + Some (Array.init np (fun l -> + let v2_l = v2.(l) in + Array.mapi (fun k v2_lk -> + expo_inv_q.(k) *. f *. v2_lk) v2_l + ) ) + end + in + Zmap.add map_2d.(0) key result; + result in let sum matrix = @@ -465,11 +461,8 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) let vrr_v a b = let v = - vrr_v 0 a b - (* - if Array.length zero_m_array < 10 then vrr_v 0 a b + if Array.length zero_m_array < 6 then vrr_v 0 a b else trr_v a b - *) in match v with | None -> 0. diff --git a/Simulation.ml b/Simulation.ml index 9809c80..2988ba3 100644 --- a/Simulation.ml +++ b/Simulation.ml @@ -1,4 +1,6 @@ type t = { + charge : Charge.t; + electrons : Electrons.t ; basis : Basis.t; nuclei : Nuclei.t; overlap : Overlap.t lazy_t; @@ -9,12 +11,21 @@ type t = { nuclear_repulsion : float; } -let make ~nuclei ~basis = +let make ?multiplicity:(multiplicity=1) ?charge:(charge=0) ~nuclei basis = + let electrons = + Electrons.make ~multiplicity ~charge nuclei + in + let charge = + Array.fold_left (fun accu (e, _) -> accu + Charge.to_int (Element.to_charge e) ) + 0 nuclei - Electrons.(electrons.n_alpha + electrons.n_beta) + |> Charge.of_int + in let overlap = lazy (Overlap.of_basis basis) in { - basis ; nuclei ; overlap ; + charge ; + basis ; nuclei ; electrons ; overlap ; overlap_ortho = lazy (Orthonormalization.make (Lazy.force overlap)); eN_ints = lazy (NucInt.of_basis_nuclei basis nuclei); kin_ints = lazy (KinInt.of_basis basis); @@ -23,12 +34,12 @@ let make ~nuclei ~basis = } -let of_filenames ~basis ~nuclei = +let of_filenames ?multiplicity:(multiplicity=1) ?charge:(charge=0) ~nuclei basis = let nuclei = Nuclei.of_filename ~filename:nuclei in let basis = Basis.of_nuclei_and_basis_filename ~nuclei ~filename:basis in - make ~nuclei ~basis + make ~charge ~multiplicity ~nuclei basis diff --git a/Utils/Util.ml b/Utils/Util.ml index 5c67741..7a66a05 100644 --- a/Utils/Util.ml +++ b/Utils/Util.ml @@ -161,3 +161,9 @@ let diagonalize_symm h = syevd ~vectors:true ~w v in v, result + +let xt_o_x o x = + gemm o x + |> gemm ~transa:`T x + + diff --git a/run_integrals.ml b/run_integrals.ml index 9c405c3..d5515c9 100644 --- a/run_integrals.ml +++ b/run_integrals.ml @@ -30,7 +30,7 @@ let run ~out = in let s = - Simulation.of_filenames ~nuclei:nuclei_file ~basis:basis_file + Simulation.of_filenames ~nuclei:nuclei_file basis_file in print_endline @@ Nuclei.to_string s.Simulation.nuclei;