From 08e8b494ea1cc86c6b3e6b33983a032b802121d7 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 19 Mar 2019 19:07:55 +0100 Subject: [PATCH] Accelerated Haa --- CI/CI.ml | 78 +++++++++++++++++++++++++++++------------ CI/Determinant.mli | 1 - CI/DeterminantSpace.ml | 72 +++++++++++++++++++++++++++++++++++++ CI/DeterminantSpace.mli | 5 +++ CI/Spindeterminant.ml | 6 +++- CI/Spindeterminant.mli | 4 +++ 6 files changed, 141 insertions(+), 25 deletions(-) diff --git a/CI/CI.ml b/CI/CI.ml index 88ffd91..e9b0a7c 100644 --- a/CI/CI.ml +++ b/CI/CI.ml @@ -1,4 +1,5 @@ open Lacaml.D +open Util module Ds = DeterminantSpace module Sd = Spindeterminant @@ -418,6 +419,7 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states } | MOClass.Core _ -> false | _ -> true ) + |> Array.of_list in let psi0 = @@ -429,7 +431,7 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states } Array.init (Ds.size det_space) (fun i -> Stream.next stream, psi0.{i+1,1}) in - (* + let is_internal = let m l = List.fold_left (fun accu i -> @@ -455,20 +457,10 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states } in Z.logand neg_active_mask beta = occ_mask in - *) let symmetric = i_o1_alfa == alfa_o2_i in - let is_internal alfa = - let rec aux = function - | -1 -> false - | j -> Determinant.degree (fst psi0.(j)) alfa = 0 - || aux (j-1) - in - aux (Array.length psi0 - 1) - in - let det_contribution i = @@ -487,7 +479,9 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states } let psi_filtered_idx = let rec aux accu = function | j when j < i -> List.rev accu - | j -> if Determinant.degree (fst psi0.(i)) (fst psi0.(j)) < 4 then + | j -> + if Determinant.degree (fst psi0.(i)) (fst psi0.(j)) < 4 + then aux (j::accu) (j-1) else aux accu (j-1) @@ -522,13 +516,12 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states } let det_i = fst psi0.(i) in let w_alfa = w_alfa det_i in - let same_spin = List.fold_left (fun accu spin -> accu +. - List.fold_left (fun accu particle -> + Array.fold_left (fun accu particle -> accu +. - List.fold_left (fun accu hole -> + Array.fold_left (fun accu hole -> if hole = particle then accu else let alfa = Determinant.single_excitation spin hole particle det_i @@ -541,12 +534,12 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states } in let double = - List.fold_left (fun accu particle' -> + Array.fold_left (fun accu particle' -> if particle' > particle || particle' = hole then accu else accu +. - List.fold_left (fun accu hole' -> + Array.fold_left (fun accu hole' -> if hole' = particle' || hole' = particle || hole' < hole then accu else @@ -570,9 +563,9 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states } in let opposite_spin = - List.fold_left (fun accu particle -> + Array.fold_left (fun accu particle -> accu +. - List.fold_left (fun accu hole -> + Array.fold_left (fun accu hole -> if hole = particle then accu else let alfa = Determinant.single_excitation Spin.Alfa hole particle det_i @@ -580,9 +573,9 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states } if Determinant.is_none alfa then accu else let double = - List.fold_left (fun accu particle' -> + Array.fold_left (fun accu particle' -> accu +. - List.fold_left (fun accu hole' -> + Array.fold_left (fun accu hole' -> if hole' = particle' then accu else let alfa = Determinant.double_excitation @@ -616,9 +609,48 @@ let pt2_en ci = let i_o1_alfa = h_ij mo_basis in - let w_alfa _ = + let w_alfa det_i = + let one_e, two_e = h_integrals mo_basis in + let fock_diag_alfa, fock_diag_beta = + Ds.fock_diag ci.det_space det_i + in + let h_aa alfa = + match Excitation.of_det det_i alfa with + | Excitation.Double (_, + {hole = h ; particle = p ; spin = s }, + {hole = h'; particle = p'; spin = s'}) -> + let fock_diag1 = + if s = Spin.Alfa then + fock_diag_alfa + else + fock_diag_beta + in + let fock_diag2 = + if s' = Spin.Alfa then + fock_diag_alfa + else + fock_diag_beta + in + fock_diag1.(0) -. fock_diag1.(h) + +. (fock_diag1.(p ) -. two_e h p h p s s) + -. (fock_diag2.(h') -. two_e h h' h h' s s' +. two_e p h' p h' s s') + +. (fock_diag2.(p') -. two_e h p' h p' s s' + +. two_e p p' p p' s s' -. two_e h' p' h' p' s' s') + | Excitation.Single (_, + {hole = h ; particle = p ; spin = s }) -> + let fock_diag = + if s = Spin.Alfa then + fock_diag_alfa + else + fock_diag_beta + in + fock_diag.(0) -. fock_diag.(h) + +. (fock_diag.(p) -. two_e h p h p s s) + |> ignore; + h_ij mo_basis alfa alfa + | _ -> assert false + in let e0 = e0.{1} in - let h_aa alfa = h_ij mo_basis alfa alfa in fun alfa -> 1. /. (e0 -. h_aa alfa) in diff --git a/CI/Determinant.mli b/CI/Determinant.mli index 839f846..2c75204 100644 --- a/CI/Determinant.mli +++ b/CI/Determinant.mli @@ -58,7 +58,6 @@ val degree : t -> t -> int (** Returns degree of excitation between two determinants. *) - (** {1 Creators} *) val of_spindeterminants : Spindeterminant.t -> Spindeterminant.t -> t diff --git a/CI/DeterminantSpace.ml b/CI/DeterminantSpace.ml index 1fcf2e5..a0818a8 100644 --- a/CI/DeterminantSpace.ml +++ b/CI/DeterminantSpace.ml @@ -154,6 +154,78 @@ let determinant t i = Determinant.of_spindeterminants alfa beta + +let fock_diag det_space det = + + let alfa_list = + Determinant.alfa det + |> Spindeterminant.to_list + in + let beta_list = + Determinant.beta det + |> Spindeterminant.to_list + in + let mo_basis = mo_basis det_space in + let mo_num = MOBasis.size mo_basis in + let one_e_ints = MOBasis.one_e_ints mo_basis + and two_e_ints = MOBasis.two_e_ints mo_basis + in + let two_e i j k l = ERI.get_phys two_e_ints i j k l in + let build_array list1 list2 = + let result = Array.make (mo_num+1) 0. in + + (* Occupied *) + List.iter (fun i -> + let x = one_e_ints.{i,i} in + result.(i) <- result.(i) +. x; + result.(0) <- result.(0) +. x; + List.iter (fun j -> + if j <> i then + begin + let x = two_e i j i j -. two_e i j j i in + result.(i) <- result.(i) +. x; + result.(0) <- result.(0) +. 0.5 *. x; + end + ) list1; + List.iter (fun j -> + begin + let x = two_e i j i j in + result.(i) <- result.(i) +. x; + result.(0) <- result.(0) +. 0.5 *. x; + end + ) list2; + ) list1; + + (* Virtuals*) + List.iter (fun i -> + if result.(i) = 0. then + begin + let x = one_e_ints.{i,i} in + result.(i) <- result.(i) +. x; + List.iter (fun j -> + let x = two_e i j i j -. two_e i j j i in + result.(i) <- result.(i) +. x; + ) list1; + List.iter (fun j -> + begin + let x = two_e i j i j in + result.(i) <- result.(i) +. x; + end + ) list2; + end + ) (Util.list_range 1 mo_num); + result + in + let alfa, beta = + build_array alfa_list beta_list, + build_array beta_list alfa_list + in + let e = alfa.(0) +. beta.(0) in + alfa.(0) <- e; + beta.(0) <- e; + alfa, beta + + let cas_of_mo_basis mo_basis n m = let s = MOBasis.simulation mo_basis in let e = Simulation.electrons s in diff --git a/CI/DeterminantSpace.mli b/CI/DeterminantSpace.mli index 3280615..f3cedc4 100644 --- a/CI/DeterminantSpace.mli +++ b/CI/DeterminantSpace.mli @@ -44,6 +44,11 @@ val determinant_stream : t -> Determinant.t Stream.t val size : t -> int (** Size of the determinant space *) +val fock_diag : t -> Determinant.t -> float array * float array +(** Returns the diagonal of the {% $\alpha$ %} and {% $\beta$ %} Fock matrices. + The zero elements contain the energy of the determinant. +*) + val fci_of_mo_basis : ?frozen_core:bool -> MOBasis.t -> t (** Creates a space of all possible ways to put [n_alfa] electrons in the {% $\alpha$ %} diff --git a/CI/Spindeterminant.ml b/CI/Spindeterminant.ml index 429f333..0070119 100644 --- a/CI/Spindeterminant.ml +++ b/CI/Spindeterminant.ml @@ -113,7 +113,6 @@ let of_list l = |> List.fold_left (fun accu p -> creation p accu) vac - let rec to_list = function | None -> [] | Some spindet -> @@ -125,6 +124,11 @@ let rec to_list = function in aux [] spindet.bitstring +let n_electrons = function + | Some t -> Z.popcount t.bitstring + | None -> 0 + + let pp_spindet n ppf = function | None -> Format.fprintf ppf "@[None@]" | Some s -> diff --git a/CI/Spindeterminant.mli b/CI/Spindeterminant.mli index 7cd3630..aee066a 100644 --- a/CI/Spindeterminant.mli +++ b/CI/Spindeterminant.mli @@ -57,6 +57,10 @@ val holes_particles_of : t -> t -> (int*int) list (** Returns the list of pairs of holes/particles in the excitation from one determinant to another. *) +val n_electrons : t -> int +(** Returns the number of electrons in the determinant. *) + + (** {1 Creation} *) val of_bitstring : ?phase:Phase.t -> Z.t -> t