From bf3ffb652d43a918df34f03568414d27a8c77493 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 18 Mar 2019 19:17:15 +0100 Subject: [PATCH] PT2 is not working yet --- CI/CI.ml | 163 +++++++++++++++++++++++++++--------- CI/DeterminantSpace.ml | 40 +++++++++ CI/DeterminantSpace.mli | 4 + CI/SpindeterminantSpace.ml | 26 ++++++ CI/SpindeterminantSpace.mli | 10 ++- 5 files changed, 203 insertions(+), 40 deletions(-) diff --git a/CI/CI.ml b/CI/CI.ml index f77450a..bc6ebe6 100644 --- a/CI/CI.ml +++ b/CI/CI.ml @@ -385,12 +385,14 @@ let make ?(n_states=1) det_space = -(* let pt2 { det_space ; m_H ; m_S2 ; eigensystem ; n_states } = + let mo_basis = Ds.mo_basis det_space in + + let mo_class = DeterminantSpace.mo_class det_space in let mo_indices = - let cls = MOClass.mo_class_array (DeterminantSpace.mo_class det_space) in - Util.list_range 1 (Ds.mo_basis det_space |> MOBasis.size) + let cls = MOClass.mo_class_array mo_class in + Util.list_range 1 (MOBasis.size mo_basis) |> List.filter (fun i -> match cls.(i) with | MOClass.Deleted _ | MOClass.Core _ -> false @@ -398,6 +400,7 @@ let pt2 { det_space ; m_H ; m_S2 ; eigensystem ; n_states } = ) in +(* Only the gournd state is computed here *) let psi0, e0 = Lazy.force eigensystem in let psi0 = @@ -405,44 +408,128 @@ let pt2 { det_space ; m_H ; m_S2 ; eigensystem ; n_states } = Ds.determinant_stream det_space in Array.init (Ds.size det_space) (fun i -> - (Stream.next stream, psi0.{i,1}) + Stream.next stream, psi0.{i+1,1}) in let e0 = e0.{1} in - let det_contribution i = - let psi_filtered = - List.filter (fun (det, _) -> - Determinant.degree det i < 4) psi0 + (* + let is_internal = + let m l = + List.fold_left (fun accu i -> + let j = i-1 in Z.(logor accu (shift_left one j)) + ) Z.zero l in - List.fold_left (fun accu spin1 -> - accu +. - List.fold_left (fun accu particle -> - accu +. - List.fold_left (fun accu hole -> - let alfa = - Determinant.single_excitation spin1 hole particle i - in - if Determinant.is_none alfa then - accu - else - let psi_h_alfa = - Array.fold_left (fun (det, coef) -> - coef *. (h_ij det alfa) - ) 0. psi_filtered - in - let h_aa = h_ij alfa alfa in - accu +. psi_h_alfa *. psi_h_alfa / (e0 - h_aa) - ) 0. mo_indices - ) 0. mo_indices - ) 0. [ Spin.Alfa ; Spin.Beta ] - + let active_mask = m (MOClass.active_mos mo_class) in + let occ_mask = m (MOClass.core_mos mo_class) in + let inactive_mask = m (MOClass.inactive_mos mo_class) in + let occ_mask = Z.logor occ_mask inactive_mask in + let neg_active_mask = Z.lognot active_mask in + fun a -> + let alfa = + Determinant.alfa a + |> Spindeterminant.bitstring + in + if Z.logand neg_active_mask alfa <> occ_mask then + false + else + let beta = + Determinant.beta a + |> Spindeterminant.bitstring + in + Z.logand neg_active_mask beta = occ_mask in + *) + + + let det_contribution i = + + 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 + aux (j::accu) (j-1) + else + aux accu (j-1) + in aux [] (Array.length psi0 - 1) + in + + let psi_filtered = + List.map (fun i -> psi0.(i)) psi_filtered_idx + in + + let psi_h alfa = + let hij = h_ij mo_basis alfa in + List.fold_left (fun accu (det, coef) -> + accu +. coef *. (hij det)) 0. psi_filtered + 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 already_generated alfa = + if is_internal alfa then + true + else + let rec aux = function + | -1 -> false + | j -> Determinant.degree (fst psi0.(j)) alfa <= 2 || aux (j-1) + in + aux (i-1) + in + + let det_i = fst psi0.(i) in + + List.fold_left (fun accu particle -> + accu +. + List.fold_left (fun accu hole -> + if hole = particle then accu else + accu +. + List.fold_left (fun accu spin -> + let alfa = + Determinant.single_excitation spin hole particle det_i + in + if Determinant.is_none alfa then accu else + + let single = + if already_generated alfa then 0. else + let h_aa = h_ij mo_basis alfa alfa in + let psi_h_alfa = psi_h alfa in + psi_h_alfa *. psi_h_alfa /. (e0 -. h_aa) + in + + let double = + List.fold_left (fun accu particle' -> + accu +. + List.fold_left (fun accu hole' -> + accu +. + List.fold_left (fun accu spin' -> + let alfa = + Determinant.double_excitation + spin hole particle + spin' hole' particle' det_i + in + if Determinant.is_none alfa || + already_generated alfa then + accu + else + let h_aa = h_ij mo_basis alfa alfa in + let psi_h_alfa = psi_h alfa in + accu +. psi_h_alfa *. psi_h_alfa /. (e0 -. h_aa) + ) 0. [ Spin.Alfa ; Spin.Beta ] + ) 0. mo_indices + ) 0. mo_indices + in + accu +. single +. double + ) 0. [ Spin.Alfa ; Spin.Beta ] + ) 0. mo_indices + ) 0. mo_indices + in + + Array.mapi (fun i (_,c_i) -> c_i *. det_contribution i) psi0 + |> Array.fold_left (+.) 0. - match det_space with - | Ds.Arbitrary -> assert false - | Ds.Spin alfa_dets beta_dets -> - Array.map ( fun alfa -> - Array.map ( fun beta -> - ) beta_dets - ) alfa_dets - *) diff --git a/CI/DeterminantSpace.ml b/CI/DeterminantSpace.ml index d329404..1fcf2e5 100644 --- a/CI/DeterminantSpace.ml +++ b/CI/DeterminantSpace.ml @@ -154,6 +154,46 @@ let determinant t i = Determinant.of_spindeterminants alfa beta +let cas_of_mo_basis mo_basis n m = + let s = MOBasis.simulation mo_basis in + let e = Simulation.electrons s in + let n_alfa = Electrons.n_alfa e + and n_beta = Electrons.n_beta e in + let det_a = + Ss.cas_of_mo_basis mo_basis n_alfa n m + and det_b = + Ss.cas_of_mo_basis mo_basis n_beta n m + in + let mo_class = Ss.mo_class det_a in + let determinants = + let det_alfa = Ss.spin_determinants det_a + and det_beta = Ss.spin_determinants det_b + in + let n_det_beta = Array.length det_beta in + let n_det_alfa = Array.length det_alfa in + + let ndet = n_det_alfa * n_det_beta in + Format.printf "Number of determinants : %d %d %d\n%!" + n_det_alfa n_det_beta ndet; + Spin (det_alfa, det_beta) + + (* + let det = Array.make n_det_alfa + (Array.init n_det_beta (fun i -> i)) + in + let index_start = Array.init (n_det_alfa+1) (fun i -> i*n_det_beta) in + let ndet = (index_start.(n_det_alfa)) in + + Format.printf "Number of determinants : %d %d %d\n%!" + n_det_alfa n_det_beta ndet; + Arbitrary { + det_alfa ; det_beta ; det ; index_start + } + *) + in + { n_alfa ; n_beta ; mo_class ; mo_basis ; determinants } + + let fci_of_mo_basis ?(frozen_core=true) mo_basis = let s = MOBasis.simulation mo_basis in let e = Simulation.electrons s in diff --git a/CI/DeterminantSpace.mli b/CI/DeterminantSpace.mli index e8e6280..3280615 100644 --- a/CI/DeterminantSpace.mli +++ b/CI/DeterminantSpace.mli @@ -51,6 +51,10 @@ val fci_of_mo_basis : ?frozen_core:bool -> MOBasis.t -> t All other MOs are untouched. *) +val cas_of_mo_basis : MOBasis.t -> int -> int -> t +(** Creates a CAS(n,m) space of determinants. +*) + (** {2 Printing} *) val pp_det_space : Format.formatter -> t -> unit diff --git a/CI/SpindeterminantSpace.ml b/CI/SpindeterminantSpace.ml index 701631f..a88fc65 100644 --- a/CI/SpindeterminantSpace.ml +++ b/CI/SpindeterminantSpace.ml @@ -24,6 +24,8 @@ let fci_of_mo_basis ?(frozen_core=true) mo_basis elec_num = and active_mask = m (MOClass.active_mos mo_class) in let neg_active_mask = Z.lognot active_mask in + (* Here we generate the FCI space and filter out unwanted determinants + with excitations involving the core electrons. This should be improved. *) let spin_determinants = Util.bit_permtutations elec_num mo_num |> List.filter (fun b -> Z.logand neg_active_mask b = occ_mask) @@ -33,6 +35,30 @@ let fci_of_mo_basis ?(frozen_core=true) mo_basis elec_num = { elec_num ; mo_basis ; mo_class ; spin_determinants } +let cas_of_mo_basis mo_basis elec_num n m = + let mo_num = MOBasis.size mo_basis in + let mo_class = MOClass.cas_sd mo_basis n m in + let m l = + List.fold_left (fun accu i -> let j = i-1 in Z.(logor accu (shift_left one j)) + ) Z.zero l + in + let active_mask = m (MOClass.active_mos mo_class) in + let occ_mask = m (MOClass.core_mos mo_class) in + let inactive_mask = m (MOClass.inactive_mos mo_class) in + let occ_mask = Z.logor occ_mask inactive_mask in + let neg_active_mask = Z.lognot active_mask in + (* Here we generate the FCI space and filter out all the unwanted determinants. + This should be improved. *) + let spin_determinants = + Util.bit_permtutations elec_num mo_num + |> List.filter (fun b -> Z.logand neg_active_mask b = occ_mask) + |> List.map (fun b -> Spindeterminant.of_bitstring b) + |> Array.of_list + in + { elec_num ; mo_basis ; mo_class ; spin_determinants } + + + let pp_spindet_space ppf t = Format.fprintf ppf "@[[ "; Array.iteri (fun i d -> Format.fprintf ppf "@[@[%8d@] @[%a@]@]@;" i diff --git a/CI/SpindeterminantSpace.mli b/CI/SpindeterminantSpace.mli index 1f6cc08..e1c2653 100644 --- a/CI/SpindeterminantSpace.mli +++ b/CI/SpindeterminantSpace.mli @@ -25,8 +25,14 @@ val mo_basis : t -> MOBasis.t (** {1 Creation} *) val fci_of_mo_basis : ?frozen_core:bool -> MOBasis.t -> int -> t -(** Create a space of all possible ways to put [n_elec] electrons in the [Active] MOs. - All other MOs are untouched. +(** Create a space of all possible ways to put [n_elec-ncore] electrons in the + [Active] MOs. All other MOs are untouched. +*) + +val cas_of_mo_basis : MOBasis.t -> int -> int -> int -> t +(** [cas_of_mo_basis mo_basis n_elec n m] creates a CAS(n,m) space of + [Active] MOs. The unoccupied MOs are [Virtual], and the occupied MOs + are [Core] and [Inactive]. *)