From a20ec081252b3c10d529b74944a375de98738e70 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sat, 11 Jan 2020 23:46:04 +0100 Subject: [PATCH] Implemented fast f12 in native code --- CI/CI.ml | 1 + CI/CIMatrixElementF12.ml | 147 ------ CI/F12CI.ml | 325 +++---------- MOBasis/HF12.ml | 993 +++++++++++++++++++++++++++++++-------- 4 files changed, 870 insertions(+), 596 deletions(-) delete mode 100644 CI/CIMatrixElementF12.ml diff --git a/CI/CI.ml b/CI/CI.ml index ab1ef95..2a87baa 100644 --- a/CI/CI.ml +++ b/CI/CI.ml @@ -290,6 +290,7 @@ let create_matrix_spin ?(nmax=2) f det_space = let singles = List.filter (fun (i,d,det_j) -> d < 2) doubles in + let triples = List.map (fun (i,_,det_j) -> (i,det_j)) triples in diff --git a/CI/CIMatrixElementF12.ml b/CI/CIMatrixElementF12.ml deleted file mode 100644 index 47019f1..0000000 --- a/CI/CIMatrixElementF12.ml +++ /dev/null @@ -1,147 +0,0 @@ -open Lacaml.D - -module De = Determinant -module Ex = Excitation -module Sp = Spindeterminant - -type t = float list - - -let non_zero integrals degree_a degree_b ki kj = - let kia = De.alfa ki and kib = De.beta ki - and kja = De.alfa kj and kjb = De.beta kj - in - let diag_element = - let mo_a = Sp.to_list kia - and mo_b = Sp.to_list kib - in - fun one_e two_e -> - let one = - (List.fold_left (fun accu i -> accu +. one_e i i Spin.Alfa) 0. mo_a) - +. - (List.fold_left (fun accu i -> accu +. one_e i i Spin.Beta) 0. mo_b) - in - let two = - let rec aux_same spin accu = function - | [] -> accu - | i :: rest -> - let new_accu = - List.fold_left (fun accu j -> accu +. two_e i j i j spin spin) accu rest - in - (aux_same [@tailcall]) spin new_accu rest - in - let rec aux_opposite accu other = function - | [] -> accu - | i :: rest -> - let new_accu = - List.fold_left (fun accu j -> accu +. two_e i j i j Spin.Alfa Spin.Beta) accu other - in - (aux_opposite [@tailcall]) new_accu other rest - in - (aux_same Spin.Alfa 0. mo_a) +. (aux_same Spin.Beta 0. mo_b) +. - (aux_opposite 0. mo_a mo_b) - in - one +. two - in - let result = - match degree_a, degree_b with - | 1, 1 -> (* alpha-beta double *) - begin - let ha, pa, phase_a = Ex.single_of_spindet kia kja in - let hb, pb, phase_b = Ex.single_of_spindet kib kjb in - let s1 = - match phase_a, phase_b with - | Phase.Pos, Phase.Pos - | Phase.Neg, Phase.Neg -> fun _ two_e -> two_e ha hb pa pb Spin.Alfa Spin.Beta - | Phase.Neg, Phase.Pos - | Phase.Pos, Phase.Neg -> fun _ two_e -> -. two_e ha hb pa pb Spin.Alfa Spin.Beta - in - let kk = Determinant.double_excitation Spin.Alfa ha pb Spin.Beta hb pa ki in - let kka = De.alfa kk and kkb = De.beta kk in - match Spindeterminant.(degree kia kka, degree kib kkb) with - | (1,1) -> - let s2 = - begin - let ha, pa, phase_a = Ex.single_of_spindet kia kka in - let hb, pb, phase_b = Ex.single_of_spindet kib kkb in - match phase_a, phase_b with - | Phase.Pos, Phase.Pos - | Phase.Neg, Phase.Neg -> fun _ two_e -> two_e ha hb pa pb Spin.Alfa Spin.Beta - | Phase.Neg, Phase.Pos - | Phase.Pos, Phase.Neg -> fun _ two_e -> -. two_e ha hb pa pb Spin.Alfa Spin.Beta - end - in fun x two_e -> 0.75 *. (s1 x two_e) +. 0.25 *. (s2 x two_e) - | _ -> fun x two_e -> 0. - end - - | 2, 0 -> (* alpha double *) - begin - let h1, p1, h2, p2, phase = Ex.double_of_spindet kia kja in - match phase with - | Phase.Pos -> fun _ two_e -> two_e h1 h2 p1 p2 Spin.Alfa Spin.Alfa - | Phase.Neg -> fun _ two_e -> -. two_e h1 h2 p1 p2 Spin.Alfa Spin.Alfa - end - - | 0, 2 -> (* beta double *) - begin - let h1, p1, h2, p2, phase = Ex.double_of_spindet kib kjb in - match phase with - | Phase.Pos -> fun _ two_e -> two_e h1 h2 p1 p2 Spin.Beta Spin.Beta - | Phase.Neg -> fun _ two_e -> -. two_e h1 h2 p1 p2 Spin.Beta Spin.Beta - end - - | 1, 0 (* alpha single *) - | 0, 1 (* beta single *) - -> fun _ _ -> 0. - | 0, 0 -> (* diagonal element *) - diag_element - - | _ -> assert false - in - List.map (fun (one_e, two_e) -> result one_e two_e) integrals - - -let make integrals ki kj = - let degree_a, degree_b = - De.degrees ki kj - in - if degree_a+degree_b > 2 then - List.map (fun _ -> 0.) integrals - else - non_zero integrals degree_a degree_b ki kj - - - - -let make_s2 ki kj = - let degree_a = De.degree_alfa ki kj in - let kia = De.alfa ki in - let kja = De.alfa kj in - if degree_a > 1 then 0. - else - let degree_b = De.degree_beta ki kj in - let kib = De.beta ki in - let kjb = De.beta kj in - match degree_a, degree_b with - | 1, 1 -> (* alpha-beta double *) - let ha, pa, phase_a = Ex.single_of_spindet kia kja in - let hb, pb, phase_b = Ex.single_of_spindet kib kjb in - if ha = pb && hb = pa then - begin - match phase_a, phase_b with - | Phase.Pos, Phase.Pos - | Phase.Neg, Phase.Neg -> -1. - | Phase.Neg, Phase.Pos - | Phase.Pos, Phase.Neg -> 1. - end - else 0. - | 0, 0 -> - let ba = Sp.bitstring kia and bb = Sp.bitstring kib in - let tmp = Bitstring.logxor ba bb in - let n_a = Bitstring.logand ba tmp |> Bitstring.popcount in - let n_b = Bitstring.logand bb tmp |> Bitstring.popcount in - let s_z = 0.5 *. float_of_int (n_a - n_b) in - float_of_int n_a +. s_z *. (s_z -. 1.) - | _ -> 0. - - diff --git a/CI/F12CI.ml b/CI/F12CI.ml index 5ee8586..2b241ca 100644 --- a/CI/F12CI.ml +++ b/CI/F12CI.ml @@ -4,7 +4,7 @@ module Ds = DeterminantSpace module De = Determinant module Sp = Spindeterminant -type t = +type t = { mo_basis : MOBasis.t ; det_space : DeterminantSpace.t ; @@ -20,245 +20,6 @@ let mo_class t = Ds.mo_class @@ det_space t let eigensystem t = Lazy.force t.eigensystem -(* - -let single_matrices hf12_integrals density = - let nocc = Mat.dim1 density in - let nvir = Mat.dim2 density in - let { HF12. - simulation ; aux_basis ; - hf12 ; hf12_anti ; - hf12_single ; hf12_single_anti } = hf12_integrals - in - let f12 = MOBasis.f12_ints aux_basis in - let eri = MOBasis.two_e_ints aux_basis in - - let d = Mat.as_vec density in - let v_s = Mat.create nocc nvir in - let v_o = Mat.create nocc nvir in - let h_o, h_s, f_o, f_s = - Mat.create nocc nvir , - Mat.create nocc nvir , - Mat.create nocc nvir , - Mat.create nocc nvir , - in - for a=1 to nvir do - for m=1 to occ do - for u=1 to nocc do - for t=1 to nocc do - let hmtau = ERI.get_phys eri m t a u - and hmtua = ERI.get_phys eri m t u a - in - v_o.{t,u} <- hmtau; - v_s.{t,u} <- hmtau -. hmtua; - done - done; - h_o.{m,a} <- dot d_o @@ Mat.as_vec v_o; - h_s.{m,a} <- dot d_s @@ Mat.as_vec v_s - for u=1 to nocc do - for t=1 to nocc do - let fmtau = ERI.get_phys f12 m t a u - and fmtua = ERI.get_phys f12 m t u a - in - v_o.{t,u} <- 0.375 *. fmtau +. 0.125 *. fmtua; - v_s.{t,u} <- 0.25 *, (fmtau -. fmtua); - done - done; - f_o.{m,a} <- dot d_o @@ Mat.as_vec v_o; - f_s.{m,a} <- dot d_s @@ Mat.as_vec v_s - done - done; -*) - -(*--- -let hf_ij_non_zero mo_basis hf12_integrals deg_a deg_b ki kj = - let integrals = [ - - let { HF12. - simulation ; aux_basis ; - hf12 ; hf12_anti ; - hf12_single ; hf12_single_anti } = hf12_integrals - in - - let kia = De.alfa ki and kib = De.beta ki in - let kja = De.alfa kj and kjb = De.beta kj in - - let mo_a = - Bitstring.logand (Sp.bitstring kia) (Sp.bitstring kja) - |> Bitstring.to_list - |> Array.of_list - and mo_b = - Bitstring.logand (Sp.bitstring kib) (Sp.bitstring kjb) - |> Bitstring.to_list - |> Array.of_list - in - - let aux_mos = - Util.list_range (MOBasis.size mo_basis) (MOBasis.size aux_basis) - |> Array.of_list - in - - let one_e _ _ _ = 0. in - - let h = MOBasis.ee_ints aux_basis in - let two_e_h i j k l s s' = - if s' <> s then - ERI.get_phys h i j k l - else - (ERI.get_phys h i j k l) -. (ERI.get_phys h i j l k) - in - - let f = MOBasis.f12_ints aux_basis in - let two_e_f i j k l s s' = - let fijkl = F12.get_phys f i j k l - and fijlk = F12.get_phys f i j l k - in - if s' <> s then - 0.325 *. fijkl +. 0.125 *. fijlk - else - 0.25 *. (fijkl -. fijlk) - in - - (* - let mo_of_s = function - | Spin.Alfa -> mo_a - | Spin.Beta -> mo_b - in - *) - - let two_e i j k l s s' = - ( - if s = s' then - hf12_anti.{i,j,k,l} -. ( - (Array.fold_left (fun accu m -> accu +. hf12_single_anti.{m,i,j,k,l}) 0. mo_a) +. - (Array.fold_left (fun accu m -> accu +. hf12_single_anti.{m,j,i,l,k}) 0. mo_b) ) - else - hf12.{i,j,k,l} -. ( - (Array.fold_left (fun accu m -> accu +. hf12_single.{m,i,j,k,l}) 0. mo_a) +. - (Array.fold_left (fun accu m -> accu +. hf12_single.{m,j,i,l,k}) 0. mo_b) ) - ) - (* - +. Array.fold_left ( fun accu a -> accu +. - (Array.fold_left ( fun accu m -> accu +. two_e_h m i m a s s) 0. (mo_of_s s) +. - Array.fold_left ( fun accu m -> accu +. two_e_h m i m a s s) 0. (mo_of_s @@ Spin.other s) ) - *. (two_e_f a j k l s s') ) 0. aux_mos - *) - in - - let three_e i j k l m n s s' s'' = - Array.fold_left (fun accu a -> accu +. two_e_h i j l a s s' *. two_e_f a k m n s' s'') 0. aux_mos - -. Array.fold_left (fun accu a -> accu +. two_e_h j i m a s' s *. two_e_f a k l n s s'') 0. aux_mos - +. Array.fold_left (fun accu a -> accu +. two_e_h j k m a s' s'' *. two_e_f a i n l s'' s ) 0. aux_mos - -. Array.fold_left (fun accu a -> accu +. two_e_h k j n a s'' s' *. two_e_f a i m l s' s ) 0. aux_mos - +. Array.fold_left (fun accu a -> accu +. two_e_h k i n a s'' s *. two_e_f a j l m s s' ) 0. aux_mos - -. Array.fold_left (fun accu a -> accu +. two_e_h i k l a s s'' *. two_e_f a j n m s'' s' ) 0. aux_mos - in - (one_e, two_e, Some three_e) - ] - in - CIMatrixElement.non_zero integrals deg_a deg_b ki kj - |> List.hd - - -let dressing_vector ~frozen_core hf12_integrals f12_amplitudes ci = - - if Parallel.master then - Printf.printf "Building matrix\n%!"; - let det_space = - ci.CI.det_space - in - let mo_basis = - Ds.mo_basis det_space - in - - let m_HF = - - let f = - match Ds.determinants det_space with - | Ds.Arbitrary _ -> CI.create_matrix_arbitrary - | Ds.Spin _ -> CI.create_matrix_spin_computed ~nmax:3 - in - f (fun deg_a deg_b ki kj -> - hf_ij_non_zero mo_basis hf12_integrals deg_a deg_b ki kj - ) det_space - - in - - Matrix.mm (Lazy.force m_HF) (Matrix.dense_of_mat f12_amplitudes) - ---- *) - - -let hf_ij_non_zero hf12_integrals deg_a deg_b ki kj = - let integrals = [ - - let { HF12. - simulation ; aux_basis ; - hf12 ; hf12_anti ; - hf12_single ; hf12_single_anti } = hf12_integrals - in - - let kia = De.alfa ki and kib = De.beta ki in - let kja = De.alfa kj and kjb = De.beta kj in - - let mo_a = - Bitstring.logand (Sp.bitstring kia) (Sp.bitstring kja) - |> Bitstring.to_list - and mo_b = - Bitstring.logand (Sp.bitstring kib) (Sp.bitstring kjb) - |> Bitstring.to_list - in - - let one_e _ _ _ = 0. in - - let two_e i j k l s s' = - if s = s' then - hf12_anti.{i,j,k,l} -. ( - (List.fold_left (fun accu m -> accu +. hf12_single_anti.{m,i,j,k,l}) 0. mo_a) +. - (List.fold_left (fun accu m -> accu +. hf12_single_anti.{m,j,i,l,k}) 0. mo_b) - ) - else - hf12.{i,j,k,l} -. ( - (List.fold_left (fun accu m -> accu +. hf12_single.{m,i,j,k,l}) 0. mo_a) +. - (List.fold_left (fun accu m -> accu +. hf12_single.{m,j,i,l,k}) 0. mo_b) - ) - in - - let h = MOBasis.ee_ints aux_basis in - let two_e_h i j k l s s' = - if s' <> s then - ERI.get_phys h l k j i - else - (ERI.get_phys h l k j i) -. (ERI.get_phys h k l j i) - in - - let f = MOBasis.f12_ints aux_basis in - let two_e_f i j k l s s' = - if s' <> s then - F12.get_phys f i j k l - else - (F12.get_phys f i j k l) -. (F12.get_phys f i j l k) - in - - let mo_of_s = function - | Spin.Alfa -> mo_a - | Spin.Beta -> mo_b - in - - let three_e i j k l m n s s' s'' = - List.fold_left (fun accu a -> accu +. two_e_h i j l a s s' *. two_e_f a k m n s' s'') 0. (mo_of_s s' ) - -. List.fold_left (fun accu a -> accu +. two_e_h j i m a s' s *. two_e_f a k l n s s'') 0. (mo_of_s s ) - +. List.fold_left (fun accu a -> accu +. two_e_h j k m a s' s'' *. two_e_f a i n l s'' s ) 0. (mo_of_s s'') - -. List.fold_left (fun accu a -> accu +. two_e_h k j n a s'' s' *. two_e_f a i m l s' s ) 0. (mo_of_s s' ) - +. List.fold_left (fun accu a -> accu +. two_e_h k i n a s'' s *. two_e_f a j l m s s' ) 0. (mo_of_s s ) - -. List.fold_left (fun accu a -> accu +. two_e_h i k l a s s'' *. two_e_f a j n m s'' s' ) 0. (mo_of_s s'') - in - - (one_e, two_e, Some three_e) - ] - in - CIMatrixElement.non_zero integrals deg_a deg_b ki kj - |> List.hd let dressing_vector ~frozen_core hf12_integrals f12_amplitudes ci = @@ -269,6 +30,11 @@ let dressing_vector ~frozen_core hf12_integrals f12_amplitudes ci = ci.CI.det_space in + let { HF12. + simulation ; aux_basis ; + f_0 ; f_1 ; f_2 ; f_3 } = hf12_integrals + in + let m_HF = let f = @@ -277,7 +43,12 @@ let dressing_vector ~frozen_core hf12_integrals f12_amplitudes ci = | Ds.Spin _ -> CI.create_matrix_spin_computed ~nmax:3 in f (fun deg_a deg_b ki kj -> - hf_ij_non_zero hf12_integrals deg_a deg_b ki kj + match deg_a + deg_b with + | 0 -> f_0 ki + | 1 -> f_1 ki kj + | 2 -> f_2 ki kj + | 3 -> f_3 ki kj + | _ -> assert false ) det_space in @@ -285,10 +56,59 @@ let dressing_vector ~frozen_core hf12_integrals f12_amplitudes ci = Matrix.mm (Lazy.force m_HF) (Matrix.dense_of_mat f12_amplitudes) +let sum l f = List.fold_left (fun accu i -> accu +. f i) 0. l + +let array_3_init d1 d2 d3 f = + let result = + Bigarray.(Array3.create Float64 fortran_layout) d1 d2 d3 + in + for k=1 to d3 do + for j=1 to d2 do + for i=1 to d1 do + result.{i,j,k} <- f i j k + done + done + done; + result + +let array_4_init d1 d2 d3 d4 f = + let result = + Bigarray.(Genarray.create Float64 fortran_layout) [| d1;d2;d3;d4 |] + in + for l=1 to d4 do + for k=1 to d3 do + for j=1 to d2 do + for i=1 to d1 do + result.{i,j,k,l} <- f i j k l + done + done + done + done; + result + +let array_5_init d1 d2 d3 d4 d5 f = + let result = + Bigarray.(Genarray.create Float64 fortran_layout) [| d1;d2;d3;d4;d5 |] + in + for m=1 to d5 do + for l=1 to d4 do + for k=1 to d3 do + for j=1 to d2 do + for i=1 to d1 do + result.{i,j,k,l,m} <- f i j k l m + done + done + done + done + done; + result + + + let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filename ?(state=1) () = - let det_space = - DeterminantSpace.fci_of_mo_basis mo_basis ~frozen_core + let det_space = + DeterminantSpace.fci_of_mo_basis mo_basis ~frozen_core in let ci = CI.make ~n_states:state det_space in @@ -302,20 +122,21 @@ let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filen Parallel.broadcast (lazy x) in + let eigensystem = lazy ( let m_H = Lazy.force ci.CI.m_H in - - let rec iteration ~state psi = + + let rec iteration ~state psi = (* Format.printf "%a@." DeterminantSpace.pp_det_space @@ CI.det_space ci; Format.printf "%a@." Matrix.pp_matrix @@ Matrix.dense_of_mat psi; *) let column_idx = iamax (Mat.to_col_vecs psi).(state-1) in - let delta = + let delta = (* delta_i = {% $\sum_j c_j H_{ij}$ %} *) dressing_vector ~frozen_core hf12_integrals psi ci |> Matrix.to_mat @@ -339,7 +160,7 @@ Format.printf "%a@." Matrix.pp_matrix @@ Matrix.dense_of_mat psi; delta.{column_idx,state} <- delta.{column_idx,state} -. delta_00; - let eigenvectors, eigenvalues = + let eigenvectors, eigenvalues = let delta = lacpy delta in Mat.scal f delta; @@ -357,8 +178,8 @@ Format.printf "%a@." Matrix.pp_matrix @@ Matrix.dense_of_mat psi; ) in - let matrix_prod c = - let w = + let matrix_prod c = + let w = Matrix.mm ~transa:`T m_H c |> Matrix.to_mat in @@ -382,7 +203,7 @@ Format.printf "%a@." Matrix.pp_matrix @@ Matrix.dense_of_mat psi; )) in - let eigenvectors = + let eigenvectors = Conventions.rephase @@ Util.remove_epsilons eigenvectors in @@ -396,11 +217,11 @@ Format.printf "%a@." Matrix.pp_matrix @@ Matrix.dense_of_mat psi; (Mat.to_col_vecs eigenvectors).(0) ) in if Parallel.master then - Printf.printf "F12 Convergence : %e %f\n" conv (eigenvalues.{state} + Printf.printf "F12 Convergence : %e %f\n" conv (eigenvalues.{state} +. Simulation.nuclear_repulsion simulation); (* - let cabs_singles = + let cabs_singles = let f = Fock.make_rhf ~density ~ao_basis:large_ao_basis in diff --git a/MOBasis/HF12.ml b/MOBasis/HF12.ml index 30794c3..2e9a6db 100644 --- a/MOBasis/HF12.ml +++ b/MOBasis/HF12.ml @@ -10,21 +10,61 @@ type q = (float, Bigarray.float64_elt, Bigarray.fortran_layout) Bigarray.Genarra type t = { simulation : Simulation.t ; aux_basis : MOBasis.t ; - hf12 : q ; (* hf12.{i,j,k,l} = \sum_{ab} - + \sum_{am} - + \sum_{bm} - *) - hf12_anti: q ; (* hf12_anti.{i,j,k,l} = - \sum_{ab} ( - ) ( - ) - + \sum_{am} ( - ) ( - ) - + \sum_{bm} ( - ) ( - ) - *) - hf12_single : q ; (* hf12.{m,i,j,k,l} = \sum_{a} *) - hf12_single_anti: q ; (* hf12_anti.{m,i,j,k,l} = - \sum_{ab} ( - ) ( - ) *) + f_0 : Determinant.t -> float ; + f_1 : Determinant.t -> Determinant.t -> float ; + f_2 : Determinant.t -> Determinant.t -> float ; + f_3 : Determinant.t -> Determinant.t -> float ; } +let sum l f = List.fold_left (fun accu i -> accu +. f i) 0. l + +let array_3_init d1 d2 d3 f = + let result = + Bigarray.(Array3.create Float64 fortran_layout) d1 d2 d3 + in + for k=1 to d3 do + for j=1 to d2 do + for i=1 to d1 do + result.{i,j,k} <- f i j k + done + done + done; + result + +let array_4_init d1 d2 d3 d4 f = + let result = + Bigarray.(Genarray.create Float64 fortran_layout) [| d1;d2;d3;d4 |] + in + for l=1 to d4 do + for k=1 to d3 do + for j=1 to d2 do + for i=1 to d1 do + result.{i,j,k,l} <- f i j k l + done + done + done + done; + result + +let array_5_init d1 d2 d3 d4 d5 f = + let result = + Bigarray.(Genarray.create Float64 fortran_layout) [| d1;d2;d3;d4;d5 |] + in + for m=1 to d5 do + for l=1 to d4 do + for k=1 to d3 do + for j=1 to d2 do + for i=1 to d1 do + result.{i,j,k,l,m} <- f i j k l m + done + done + done + done + done; + result + + let make ~simulation ~mo_basis ~aux_basis_filename () = let f12 = Util.of_some @@ Simulation.f12 simulation in @@ -53,209 +93,768 @@ let make ~simulation ~mo_basis ~aux_basis_filename () = let aux_num = MOBasis.size aux_basis in (* Fire calculation of F12 and ERI *) - let f12 = - MOBasis.f12_ints aux_basis - in - let eri = - MOBasis.two_e_ints aux_basis - in + ignore @@ MOBasis.f12_ints aux_basis ; + ignore @@ MOBasis.two_e_ints aux_basis ; (* Compute the integrals *) if Parallel.master then Printf.eprintf "Computing HF12 integrals\n%!"; - let hf12, hf12_anti, hf12_single, hf12_single_anti = - let create4 n = - Bigarray.Genarray.create Float64 Bigarray.fortran_layout [| n ; n ; n ; n |] - in - let create5 n = - Bigarray.Genarray.create Float64 Bigarray.fortran_layout [| n ; n ; n ; n ; n |] - in - create4 mo_num , - create4 mo_num , - create5 mo_num , - create5 mo_num + let mos_cabs = + Util.list_range (mo_num+1) aux_num in - - let h_s = Bigarray.Array3.create Float64 Bigarray.fortran_layout mo_num aux_num aux_num in - let f_s = Bigarray.Array3.create Float64 Bigarray.fortran_layout aux_num aux_num mo_num in - let h_o = Bigarray.Array3.create Float64 Bigarray.fortran_layout mo_num aux_num aux_num in - let f_o = Bigarray.Array3.create Float64 Bigarray.fortran_layout aux_num aux_num mo_num in - let hf_s = Mat.create mo_num mo_num in - let hf_o = Mat.create mo_num mo_num in - - let hfm_s = Bigarray.Array3.create Float64 Bigarray.fortran_layout mo_num mo_num mo_num in - let hfm_o = Bigarray.Array3.create Float64 Bigarray.fortran_layout mo_num mo_num mo_num in - - for a=1 to mo_num do - for b=1 to mo_num do - for i=1 to mo_num do - h_s.{i, a, b} <- 0. ; - h_o.{i, a, b} <- 0. - done - done - done; - - for k=1 to mo_num do - for b=1 to mo_num do - for a=1 to mo_num do - f_s.{a, b, k} <- 0. ; - f_o.{a, b, k} <- 0. - done - done - done; - - let task (j,l) = - - let h i a b = - let ijab = ERI.get_phys eri i j a b - and ijba = ERI.get_phys eri i j b a - in - h_s.{i, a, b} <- ijab -. ijba ; - h_o.{i, a, b} <- ijab - and f a b k = - let abkl = F12.get_phys f12 a b k l - and ablk = F12.get_phys f12 a b l k - in - f_s.{a, b, k} <- 0.25 *. (abkl -. ablk) ; - f_o.{a, b, k} <- 0.375 *. abkl +. 0.125 *. ablk - in - - for a=mo_num+1 to aux_num do - for b=mo_num+1 to aux_num do - for i=1 to mo_num do - h i a b - done - done - done; - - for k=1 to mo_num do - for b=mo_num+1 to aux_num do - for a=mo_num+1 to aux_num do - f a b k - done - done - done; - -(* - let h i a b = - h_s.{i, a, b} <- 0. ; - h_o.{i, a, b} <- 0. - and f a b k = - f_s.{a, b, k} <- 0. ; - f_o.{a, b, k} <- 0. - in -*) - - for m=1 to mo_num do - for a=mo_num+1 to aux_num do - for i=1 to mo_num do - h i a m ; - h i m a - done - done - done; - - for k=1 to mo_num do - for m=1 to mo_num do - for a=mo_num+1 to aux_num do - f a m k ; - f m a k - done - done - done; - - let h_o' = - Bigarray.(reshape (genarray_of_array3 h_o)) [| mo_num ; aux_num*aux_num |] - |> Bigarray.array2_of_genarray - in - let f_o' = - Bigarray.(reshape (genarray_of_array3 f_o)) [| aux_num*aux_num ; mo_num |] - |> Bigarray.array2_of_genarray - in - let h_s' = - Bigarray.(reshape (genarray_of_array3 h_s)) [| mo_num ; aux_num*aux_num |] - |> Bigarray.array2_of_genarray - in - let f_s' = - Bigarray.(reshape (genarray_of_array3 f_s)) [| aux_num*aux_num ; mo_num |] - |> Bigarray.array2_of_genarray - in - let hf_s = gemm ~c:hf_s h_s' f_s' in - let hf_o = gemm ~c:hf_o h_o' f_o' in - - let () = - for m=1 to mo_num do - let h_o' = - Mat.init_cols mo_num aux_num (fun i a -> h_o.{i,m,a}) - in - let f_o' = - Mat.init_cols aux_num mo_num (fun a k -> f_o.{m,a,k}) - in - let h_s' = - Mat.init_cols mo_num aux_num (fun i a -> h_s.{i,m,a}) - in - let f_s' = - Mat.init_cols aux_num mo_num (fun a k -> f_s.{m,a,k}) - in - let r_s, r_o = - gemm h_s' f_s' , - gemm h_o' f_o' - in - for k = 1 to mo_num do - for i = 1 to mo_num do - hfm_s.{m,i,k} <- r_s.{i,k}; - hfm_o.{m,i,k} <- r_o.{i,k} - done - done - done - in - hf_s, hf_o, hfm_s, hfm_o, j, l + let mos_in = + Util.list_range 1 mo_num in - let tasks = - let rec next accu = function - | _, 0 -> accu - | 0, l -> next accu (mo_num, l-1) - | j, l -> next ((j,l) :: accu) ((j-1), l) + let mos_a k = + Determinant.alfa k + |> Spindeterminant.to_list + in + + let mos_b k = + Determinant.beta k + |> Spindeterminant.to_list + in + + let h_one = + let h = + MOBasis.one_e_ints aux_basis + in fun i j _ -> h.{i,j} + in + + let h_two = + let two_e_ints = MOBasis.two_e_ints aux_basis in + let h2 i j k l (s:Spin.t) (s':Spin.t) = + if s' <> s then + ERI.get_phys two_e_ints i j k l + else + (ERI.get_phys two_e_ints i j k l) -. + (ERI.get_phys two_e_ints i j l k) in - next [] (mo_num, mo_num) - |> Stream.of_list + h2 in + let f_two = + let f12_ints = MOBasis.f12_ints aux_basis in + let f2 i j k l (s:Spin.t) (s':Spin.t) = + if s' <> s then + 0.375 *. F12.get_phys f12_ints i j k l +. + 0.125 *. F12.get_phys f12_ints i j l k + else + 0.25 *. ( + (F12.get_phys f12_ints i j k l) -. + (F12.get_phys f12_ints i j l k) ) + in + f2 + in - Farm.run ~f:task ~ordered:true tasks - |> Stream.iter (fun (hf_s, hf_o, hfm_s, hfm_o, j, l) -> - for k=1 to mo_num do - for i=1 to mo_num do - hf12.{i,j,k,l} <- hf_o.{i,k} ; - hf12_anti.{i,j,k,l} <- hf_s.{i,k} ; - for m=1 to mo_num do - hf12_single.{m,i,j,k,l} <- hfm_o.{m,i,k} ; - hf12_single_anti.{m,i,j,k,l} <- hfm_s.{m,i,k} - done - done - done ); + let f_one = fun _ _ _ -> 0. in - (* - for l=1 to mo_num do - for k=1 to mo_num do - for j=1 to mo_num do - for i=1 to mo_num do - Printf.printf "%d %d %d %d %e\n" i j k l result.{i,j,k,l} - done - done - done - done; - Printf.printf "%!"; - *) + (* Pre-compute dressed integrals *) + let m_0111_1H_1F = + Vec.init mo_num (fun i -> + sum mos_cabs (fun a -> + h_one i a Spin.Alfa *. f_one a i Spin.Alfa )) + in + + let m_0111_1H_2Fa, m_0111_2Ha_2Fa = + + let m_0122_Haa = + array_3_init mo_num mo_num mo_num (fun i j k -> + sum mos_cabs (fun a -> + h_two i j a k Spin.Alfa Spin.Alfa *. f_two a k i j Spin.Alfa Spin.Alfa + ) ) + in + + let m_0111_1H_2Fa = + Mat.init_cols mo_num mo_num (fun i j -> + sum mos_cabs (fun a -> + h_one i a Spin.Alfa *. f_two a j i j Spin.Alfa Spin.Alfa +. + h_two i j a j Spin.Alfa Spin.Alfa *. f_one a i Spin.Alfa + ) +. + if i < j then 0. else + begin + sum mos_cabs (fun a -> + sum mos_cabs (fun b -> if b >= a then 0. else + h_two i j a b Spin.Alfa Spin.Alfa *. f_two a b i j Spin.Alfa Spin.Alfa + ) + ) +. + sum mos_in (fun k -> m_0122_Haa.{i,j,k}) + end + ) + in + + let m_0111_2Ha_2Fa = + array_3_init mo_num mo_num mo_num (fun i j k -> + sum mos_cabs (fun a -> + h_two i j a j Spin.Alfa Spin.Alfa *. + f_two a k i k Spin.Alfa Spin.Alfa + ) -. if i < j then 0. else m_0122_Haa.{i,j,k} + ) + in m_0111_1H_2Fa, m_0111_2Ha_2Fa + in + + let m_0111_1H_2Fb, m_0111_2Hb_2Fb = + let m_0122_Hab = + array_3_init mo_num mo_num mo_num (fun i j k -> + sum mos_cabs (fun a -> + h_two i j a k Spin.Alfa Spin.Beta *. f_two a k i j Spin.Alfa Spin.Beta + ) ) + in + + let m_0111_1H_2Fb = + Mat.init_cols mo_num mo_num (fun i j -> + sum mos_cabs (fun a -> + h_one i a Spin.Alfa *. f_two a j i j Spin.Alfa Spin.Beta +. + h_two i j a j Spin.Alfa Spin.Beta *. f_one a i Spin.Alfa +. + h_one j a Spin.Alfa *. f_two a i j i Spin.Alfa Spin.Beta +. + h_two j i a i Spin.Alfa Spin.Beta *. f_one a j Spin.Alfa + ) +. + sum mos_in (fun k -> m_0122_Hab.{i,j,k} +. m_0122_Hab.{j,i,k} ) +. + sum mos_cabs (fun a -> + sum mos_cabs (fun b -> + h_two i j a b Spin.Alfa Spin.Beta *. f_two a b i j Spin.Alfa Spin.Beta + ) + ) + ) + in + + let m_0111_2Hb_2Fb = + array_3_init mo_num mo_num mo_num (fun i j k -> + sum mos_cabs (fun a -> + h_two k i a i Spin.Alfa Spin.Beta *. + f_two a j k j Spin.Alfa Spin.Beta +. + h_two i k a k Spin.Alfa Spin.Beta *. + f_two a j i j Spin.Alfa Spin.Alfa + ) -. m_0122_Hab.{k,i,j} + ) + in + m_0111_1H_2Fb, m_0111_2Hb_2Fb + in + + let m_0111_2Ha_2Fb = + array_3_init mo_num mo_num mo_num (fun i j k -> + sum mos_cabs (fun a -> + h_two i j a j Spin.Alfa Spin.Alfa *. + f_two a k i k Spin.Alfa Spin.Beta + ) + ) + in + + let f_0 ki = + + let mos_i, mos_i' = mos_a ki, mos_b ki in + let same = (mos_i = mos_i') in + + (* Alpha *) + let a = + sum mos_i (fun i -> m_0111_1H_1F.{i}) + in + + let b = + if same then a else + sum mos_i' (fun i -> m_0111_1H_1F.{i}) + in + + let aa = + sum mos_i (fun j -> + sum mos_i (fun i -> m_0111_1H_2Fa.{i,j} )) + in + + let bb = + if same then aa else + sum mos_i' (fun j -> + sum mos_i' (fun i -> m_0111_1H_2Fa.{i,j} )) + in + + let ab = + sum mos_i' (fun j -> + sum mos_i (fun i -> m_0111_1H_2Fb.{i,j} )) + in + + let aaa = + sum mos_i (fun k -> + sum mos_i (fun j -> + sum mos_i (fun i -> m_0111_2Ha_2Fa.{i,j,k} ))) + in + + let bbb = + if same then aaa else + sum mos_i' (fun k -> + sum mos_i' (fun j -> + sum mos_i' (fun i -> m_0111_2Ha_2Fa.{i,j,k} ))) + in + + let baa = + sum mos_i' (fun k -> + sum mos_i (fun j -> + sum mos_i (fun i -> + m_0111_2Ha_2Fb.{i,j,k} +. m_0111_2Hb_2Fb.{i,j,k} + ))) + in + + let bba = + sum mos_i (fun k -> + sum mos_i' (fun j -> + sum mos_i' (fun i -> + m_0111_2Ha_2Fb.{i,j,k} +. m_0111_2Hb_2Fb.{j,i,k} + ))) + in + + a +. b +. aa +. bb +. ab +. aaa +. baa +. bba +. bbb + in + + let m_1111_1H_1F = + Mat.init_cols mo_num mo_num (fun i k -> + sum mos_cabs (fun a -> h_one i a Spin.Alfa *. f_one a k Spin.Alfa )) + in + + let m_1111_2Ha_2Fa = + array_4_init mo_num mo_num mo_num mo_num (fun l j i k -> + if l=i then + sum mos_cabs (fun a -> + h_two j l a l Spin.Alfa Spin.Alfa *. + f_two a i j k Spin.Alfa Spin.Alfa + ) + else + sum mos_cabs (fun a -> + h_two i j a j Spin.Alfa Spin.Alfa *. + f_two a l k l Spin.Alfa Spin.Alfa +. + h_two j l a l Spin.Alfa Spin.Alfa *. + f_two a i j k Spin.Alfa Spin.Alfa ) + ) + in + + let m_1111_2Hb_2Fa = + array_4_init mo_num mo_num mo_num mo_num (fun l j i k -> + if l=i then + sum mos_cabs (fun a -> + h_two j l a l Spin.Alfa Spin.Beta *. + f_two a i j k Spin.Alfa Spin.Beta + ) + else + sum mos_cabs (fun a -> + h_two i j a j Spin.Alfa Spin.Beta *. + f_two a l k l Spin.Alfa Spin.Alfa +. + h_two j l a l Spin.Alfa Spin.Beta *. + f_two a i j k Spin.Alfa Spin.Beta + ) + ) + in + + let m_1111_2Ha_2Fb = + array_4_init mo_num mo_num mo_num mo_num (fun l j i k -> + sum mos_cabs (fun a -> + h_two i j a j Spin.Alfa Spin.Alfa *. + f_two a l k l Spin.Alfa Spin.Beta +. + h_two j l a l Spin.Alfa Spin.Beta *. + f_two a i j k Spin.Alfa Spin.Alfa + ) + ) + in + + let m_1111_2Hb_2Fb = + array_4_init mo_num mo_num mo_num mo_num (fun l j i k -> + sum mos_cabs (fun a -> + h_two i j a j Spin.Alfa Spin.Beta *. + f_two a l k l Spin.Alfa Spin.Beta +. + h_two j l a l Spin.Alfa Spin.Alfa *. + f_two a i j k Spin.Alfa Spin.Beta + ) + ) + in + + let m_1121_2Ha_2Fa = + array_4_init mo_num mo_num mo_num mo_num (fun l j i k -> + sum mos_cabs (fun a -> + h_two i j k a Spin.Alfa Spin.Alfa *. + f_two l a l j Spin.Alfa Spin.Alfa + ) + ) + in + + let m_1121_2Hb_2Fa = + array_4_init mo_num mo_num mo_num mo_num (fun l j i k -> + sum mos_cabs (fun a -> + h_two i j k a Spin.Alfa Spin.Beta *. + f_two l a l j Spin.Alfa Spin.Alfa + ) + ) + in + + let m_1121_2Ha_2Fb = + array_4_init mo_num mo_num mo_num mo_num (fun l j i k -> + sum mos_cabs (fun a -> + h_two i j k a Spin.Alfa Spin.Alfa *. + f_two l a l j Spin.Alfa Spin.Beta + ) + ) + in + + let m_1121_2Hb_2Fb = + array_4_init mo_num mo_num mo_num mo_num (fun l j i k -> + sum mos_cabs (fun a -> + h_two i j k a Spin.Alfa Spin.Beta *. + f_two l a l j Spin.Alfa Spin.Beta + ) + ) + in + + let m_1122_va = + array_4_init mo_num mo_num mo_num mo_num (fun l j i k -> + sum mos_cabs (fun a -> + h_two j i a l Spin.Alfa Spin.Alfa *. + f_two l a k j Spin.Alfa Spin.Alfa + ) + ) + in + + let m_1111_1H_2Fa = + array_3_init mo_num mo_num mo_num (fun j i k -> + sum mos_in (fun l -> m_1122_va.{l,j,i,k}) +. + sum mos_cabs (fun a -> + h_one i a Spin.Alfa *. f_two a j k j Spin.Alfa Spin.Alfa +. + h_two i j a j Spin.Alfa Spin.Alfa *. f_one a k Spin.Alfa +. + h_one j a Spin.Alfa *. f_two a i j k Spin.Alfa Spin.Alfa +. + h_two i j k a Spin.Alfa Spin.Alfa *. f_one a j Spin.Alfa +. + sum mos_cabs (fun b -> if b > a then 0. else + h_two i j a b Spin.Alfa Spin.Alfa *. + f_two a b k j Spin.Alfa Spin.Alfa ) + ) + ) + in + + let m_1122_v2 = + array_4_init mo_num mo_num mo_num mo_num (fun l j i k -> + sum mos_cabs (fun a -> + h_two i j l a Spin.Alfa Spin.Beta *. + f_two l a k j Spin.Alfa Spin.Beta + ) + ) + in + + let m_1122_v3 = + array_4_init mo_num mo_num mo_num mo_num (fun l j i k -> + sum mos_cabs (fun a -> + h_two i j a l Spin.Alfa Spin.Beta *. + f_two a l k j Spin.Alfa Spin.Beta + ) + ) + in + + let m_1111_1H_2Fb = + array_3_init mo_num mo_num mo_num (fun j i k -> + sum mos_in (fun l -> m_1122_v2.{l,j,i,k} +. m_1122_v3.{l,j,i,k}) +. + sum mos_cabs (fun a -> + h_one i a Spin.Alfa *. f_two a j k j Spin.Alfa Spin.Beta +. + h_two i j a j Spin.Alfa Spin.Beta *. f_one a k Spin.Alfa +. + h_one j a Spin.Beta *. f_two a i j k Spin.Alfa Spin.Beta +. + h_two i j k a Spin.Alfa Spin.Beta *. f_one a j Spin.Beta +. + sum mos_cabs (fun b -> + h_two i j a b Spin.Alfa Spin.Beta *. + f_two a b k j Spin.Alfa Spin.Beta + ) + ) + ) + in + + let m_1122_oa = + array_4_init mo_num mo_num mo_num mo_num (fun l j i k -> + if l > j then + sum mos_cabs (fun a -> + h_two j l a k Spin.Alfa Spin.Alfa *. + f_two i a l j Spin.Alfa Spin.Alfa + ) + else 0. + ) + in + + let m_1122_o = + array_4_init mo_num mo_num mo_num mo_num (fun l j i k -> + sum mos_cabs (fun a -> + h_two l j k a Spin.Alfa Spin.Beta *. + f_two i a l j Spin.Alfa Spin.Beta + ) + ) + in + + let f_1 ki kj = + let i, k, s, phase = + match Excitation.of_det ki kj with + | Excitation.(Single (phase, { hole ; particle ; spin })) -> + hole, particle, spin, phase + | _ -> assert false + in + + let mos_novirt, mos_novirt' = + let alfa = + let i = Spindeterminant.bitstring @@ Determinant.alfa ki in + let j = Spindeterminant.bitstring @@ Determinant.alfa kj in + Bitstring.to_list (Bitstring.logor i j) + in + let beta = + let i = Spindeterminant.bitstring @@ Determinant.beta ki in + let j = Spindeterminant.bitstring @@ Determinant.beta kj in + Bitstring.to_list (Bitstring.logor i j) + in + match s with + | Spin.Alfa -> alfa, beta + | Spin.Beta -> beta, alfa + in + + let mos_ij, mos_ij' = + let alfa = + let i = Spindeterminant.bitstring @@ Determinant.alfa ki in + let j = Spindeterminant.bitstring @@ Determinant.alfa kj in + Bitstring.to_list (Bitstring.logand i j) + in + let beta = + let i = Spindeterminant.bitstring @@ Determinant.beta ki in + let j = Spindeterminant.bitstring @@ Determinant.beta kj in + Bitstring.to_list (Bitstring.logand i j) + in + match s with + | Spin.Alfa -> alfa, beta + | Spin.Beta -> beta, alfa + in + + let mos_i, mos_i' = + match s with + | Spin.Alfa -> mos_a ki, mos_b ki + | Spin.Beta -> mos_b ki, mos_a ki + in + + let mos_j, mos_j' = + match s with + | Spin.Alfa -> mos_a kj, mos_b kj + | Spin.Beta -> mos_b kj, mos_a kj + in + + let result = + m_1111_1H_1F.{i,k} +. + sum mos_ij (fun j -> + m_1111_1H_2Fa.{j,i,k} + +. sum mos_i (fun l -> m_1111_2Ha_2Fa.{l,j,i,k}) + +. sum mos_i' (fun l -> m_1111_2Ha_2Fb.{l,j,i,k}) + +. sum mos_j (fun l -> m_1121_2Ha_2Fa.{l,j,i,k}) + +. sum mos_j' (fun l -> m_1121_2Ha_2Fb.{l,j,i,k}) + -. sum mos_novirt (fun l -> m_1122_va.{l,j,i,k}) + -. sum mos_ij (fun l -> m_1122_oa.{l,j,i,k} ) + ) +. + sum mos_ij' (fun j -> + m_1111_1H_2Fb.{j,i,k} + +. sum mos_i (fun l -> m_1111_2Hb_2Fa.{l,j,i,k}) + +. sum mos_i' (fun l -> m_1111_2Hb_2Fb.{l,j,i,k}) + +. sum mos_j (fun l -> m_1121_2Hb_2Fb.{l,j,i,k}) + +. sum mos_j' (fun l -> m_1121_2Hb_2Fa.{l,j,i,k}) + -. sum mos_novirt (fun l -> m_1122_v2.{l,j,i,k}) + -. sum mos_novirt' (fun l -> m_1122_v3.{l,j,i,k}) + -. sum mos_ij (fun l -> m_1122_o.{l,j,i,k}) + ) + in + match phase with + | Phase.Pos -> result + | Phase.Neg -> -. result + in + + let m_2112_1H_2Fa = + array_4_init mo_num mo_num mo_num mo_num (fun i j k l -> + sum mos_cabs (fun a -> + h_one i a Spin.Alfa *. f_two a j k l Spin.Alfa Spin.Alfa +. + h_one j a Spin.Alfa *. f_two i a k l Spin.Alfa Spin.Alfa +. + h_two i j a l Spin.Alfa Spin.Alfa *. f_one a k Spin.Alfa +. + h_two i j k a Spin.Alfa Spin.Alfa *. f_one a l Spin.Alfa +. + sum mos_in (fun m -> -. h_two i j a m Spin.Alfa Spin.Alfa *. + f_two m a k l Spin.Alfa Spin.Alfa) +. + sum mos_cabs (fun b -> if b >= a then 0. else + h_two i j a b Spin.Alfa Spin.Alfa *. f_two a b k l Spin.Alfa Spin.Alfa + ) + ) + ) + in + + let m_2112_1H_2Fb = + array_4_init mo_num mo_num mo_num mo_num (fun i j k l -> + sum mos_cabs (fun a -> + h_one i a Spin.Alfa *. f_two a j k l Spin.Alfa Spin.Beta +. + h_one j a Spin.Alfa *. f_two i a k l Spin.Alfa Spin.Beta +. + h_two i j a l Spin.Alfa Spin.Beta *. f_one a k Spin.Alfa +. + h_two i j k a Spin.Alfa Spin.Beta *. f_one a l Spin.Alfa +. + sum mos_in (fun m -> + h_two i j a m Spin.Alfa Spin.Beta *. f_two a m k l Spin.Alfa Spin.Beta +. + h_two i j m a Spin.Alfa Spin.Beta *. f_two m a k l Spin.Alfa Spin.Beta ) +. + sum mos_cabs (fun b -> + h_two i j a b Spin.Alfa Spin.Beta *. f_two a b k l Spin.Alfa Spin.Beta + ) + ) + ) + in + + let m_2112_2Ha_2Fa = + array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l -> + sum mos_cabs (fun a -> + h_two i n a n Spin.Alfa Spin.Alfa *. f_two a j k l Spin.Alfa Spin.Alfa +. + h_two j n a n Spin.Alfa Spin.Alfa *. f_two a i l k Spin.Alfa Spin.Alfa + ) + ) + in + + let m_2112_2Hb_2Fa = + array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l -> + sum mos_cabs (fun a -> + h_two i n a n Spin.Alfa Spin.Beta *. f_two a j k l Spin.Alfa Spin.Alfa +. + h_two j n a n Spin.Alfa Spin.Beta *. f_two a i l k Spin.Alfa Spin.Alfa + ) + ) + in + + let m_2112_2Ha_2Fb = + array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l -> + sum mos_cabs (fun a -> + h_two i n a n Spin.Alfa Spin.Alfa *. f_two a j k l Spin.Alfa Spin.Beta +. + h_two j n a n Spin.Alfa Spin.Beta *. f_two a i l k Spin.Alfa Spin.Beta ) + ) + in + + let m_2121_2Ha_2Fa = + array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l -> + sum mos_cabs (fun a -> + h_two i j a l Spin.Alfa Spin.Alfa *. f_two a n k n Spin.Alfa Spin.Alfa +. + h_two j i a k Spin.Alfa Spin.Alfa *. f_two a n l n Spin.Alfa Spin.Alfa + ) + ) + in + + let m_2121_2Hb_2Fa = + array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l -> + sum mos_cabs (fun a -> + h_two i j a l Spin.Alfa Spin.Beta *. f_two a n k n Spin.Alfa Spin.Alfa +. + h_two j i a k Spin.Alfa Spin.Beta *. f_two a n l n Spin.Alfa Spin.Beta + ) + ) + in + + let m_2121_2Ha_2Fb = + array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l -> + sum mos_cabs (fun a -> + h_two i j a l Spin.Alfa Spin.Alfa *. f_two a n k n Spin.Alfa Spin.Beta +. + h_two j i a k Spin.Alfa Spin.Alfa *. f_two a n l n Spin.Alfa Spin.Beta + ) + ) + in + + let m_2122_2Ha_2Fa_ij = + let s = Spin.Alfa in + array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l -> + sum mos_cabs (fun a -> + h_two i n a k s s *. f_two j a n l s s + +. h_two i n a l s s *. f_two j a k n s s + -. h_two j n a k s s *. f_two i a n l s s + -. h_two j n a l s s *. f_two i a k n s s + ) + ) + in + + let m_2122_2Hb_2Fb_ij = + let s, s' = Spin.(Alfa, Beta) in + array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l -> + sum mos_cabs (fun a -> + h_two n i a k s s *. f_two a j n l s s' + +. h_two n j a l s s' *. f_two i a k n s s + -. h_two n j k a s s' *. f_two i a n l s s' + ) + ) + in + + let m_2122_2Hb_2Fb_ij2 = + let s, s' = Spin.(Alfa, Beta) in + array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l -> + sum mos_cabs (fun a -> + -. h_two i n a l s s' *. f_two a j k n s s' +. + (if n < j then + h_two i n k a s s' *. f_two j a l n s' s' + +. h_two n j a l s' s' *. f_two i a k n s s' + else + -. h_two i n k a s s' *. f_two j a n l s' s' + -. h_two j n a l s' s' *. f_two i a k n s s' + ) + ) + ) + in + + let m_2122_2Ha_2Fa_ij2 = + let s, s' = Spin.(Alfa, Beta) in + array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l -> + sum mos_cabs (fun a -> + h_two i n k a s s' *. f_two j a l n s s' + +. h_two j n l a s s' *. f_two i a k n s s' + -. h_two i n l a s s' *. f_two j a k n s s' + -. h_two j n k a s s' *. f_two i a l n s s' + ) + ) + in + + let m_2122_2Ha_2Fa_nv = + let s = Spin.Alfa in + array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l -> + sum mos_cabs (fun a -> h_two i j a n s s *. f_two n a k l s s ) ) + in + + let m_2122_2Hb_2Fb_nv = + let s, s' = Spin.(Alfa, Beta) in + array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l -> + sum mos_cabs (fun a -> -. h_two i j a n s s' *. f_two a n k l s s' ) ) + in + + let m_2122_2Hb_2Fb_nv2 = + let s, s' = Spin.(Alfa, Beta) in + array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l -> + sum mos_cabs (fun a -> -. h_two i j n a s s' *. f_two n a k l s s' ) ) + in + + let f_2 ki kj = + let i, j, k, l, s, s', phase = + match Excitation.of_det ki kj with + | Excitation.(Double (phase, { hole ; particle ; spin }, {hole=hole' ; particle=particle' ; spin=spin' })) -> + hole, hole', particle, particle', spin, spin', phase + | _ -> assert false + in + let mos_i, mos_i' = + match s with + | Spin.Alfa -> mos_a ki, mos_b ki + | Spin.Beta -> mos_b ki, mos_a ki + in + + let mos_j, mos_j' = + match s with + | Spin.Alfa -> mos_a kj, mos_b kj + | Spin.Beta -> mos_b kj, mos_a kj + in + + let mos_ij, mos_ij' = + let alfa = + let i = Spindeterminant.bitstring @@ Determinant.alfa ki in + let j = Spindeterminant.bitstring @@ Determinant.alfa kj in + Bitstring.to_list (Bitstring.logand i j) + in + let beta = + let i = Spindeterminant.bitstring @@ Determinant.beta ki in + let j = Spindeterminant.bitstring @@ Determinant.beta kj in + Bitstring.to_list (Bitstring.logand i j) + in + match s with + | Spin.Alfa -> alfa, beta + | Spin.Beta -> beta, alfa + in + + let mos_novirt, mos_novirt' = + let alfa = + let i = Spindeterminant.bitstring @@ Determinant.alfa ki in + let j = Spindeterminant.bitstring @@ Determinant.alfa kj in + Bitstring.to_list (Bitstring.logor i j) + in + let beta = + let i = Spindeterminant.bitstring @@ Determinant.beta ki in + let j = Spindeterminant.bitstring @@ Determinant.beta kj in + Bitstring.to_list (Bitstring.logor i j) + in + match s with + | Spin.Alfa -> alfa, beta + | Spin.Beta -> beta, alfa + in + + let result = + if s = s' then + m_2112_1H_2Fa.{i,j,k,l} +. + sum mos_i (fun n -> m_2112_2Ha_2Fa.{n,i,j,k,l} ) +. + sum mos_i' (fun n -> m_2112_2Hb_2Fa.{n,i,j,k,l} ) +. + sum mos_j (fun n -> m_2121_2Ha_2Fa.{n,i,j,k,l} ) +. + sum mos_j' (fun n -> m_2121_2Ha_2Fb.{n,i,j,k,l} ) +. + sum mos_ij (fun n -> m_2122_2Ha_2Fa_ij.{n,i,j,k,l} ) +. + sum mos_ij' (fun n -> m_2122_2Ha_2Fa_ij2.{n,i,j,k,l} ) +. + sum mos_novirt (fun n -> m_2122_2Ha_2Fa_nv.{n,i,j,k,l} ) + else + m_2112_1H_2Fb.{i,j,k,l} +. + sum mos_i (fun n -> m_2112_2Ha_2Fb.{n,i,j,k,l} ) +. + sum mos_i' (fun n -> m_2112_2Ha_2Fb.{n,j,i,l,k} ) +. + sum mos_j (fun n -> m_2121_2Hb_2Fa.{n,i,j,k,l} ) +. + sum mos_j' (fun n -> m_2121_2Hb_2Fa.{n,j,i,l,k} ) +. + sum mos_novirt'(fun n -> m_2122_2Hb_2Fb_nv.{n,i,j,k,l} ) +. + sum mos_novirt (fun n -> m_2122_2Hb_2Fb_nv2.{n,i,j,k,l} )+. + sum mos_ij (fun n -> m_2122_2Hb_2Fb_ij.{n,i,j,k,l} ) +. + sum mos_ij' (fun n -> m_2122_2Hb_2Fb_ij2.{n,i,j,k,l} ) + in + + match phase with + | Phase.Pos -> result + | Phase.Neg -> -. result + + in + + let f_3 ki kj = + let i, j, m, k, l, n, s1, s2, s3, phase = + match Excitation.of_det ki kj with + | Excitation.(Triple (phase, + { hole=h1 ; particle=p1 ; spin=s1 }, + { hole=h2 ; particle=p2 ; spin=s2 }, + { hole=h3 ; particle=p3 ; spin=s3 }) ) -> + h1, h2, h3, p1, p2, p3, s1, s2, s3, phase + | _ -> assert false + in + + let result = + let open Spin in + match s1, s2, s3 with + | Alfa, Alfa, Alfa + | Beta, Beta, Beta -> + sum mos_cabs (fun a -> + h_two i j a k s1 s2 *. f_two m a l n s3 s3 + +. h_two i j a n s1 s2 *. f_two m a k l s3 s2 + +. h_two i m a l s1 s3 *. f_two j a k n s2 s3 + +. h_two j m a k s2 s3 *. f_two i a l n s1 s3 + +. h_two j m a n s2 s3 *. f_two i a k l s1 s2 + -. h_two i j a l s1 s2 *. f_two m a k n s3 s3 + -. h_two i m a k s1 s3 *. f_two j a l n s2 s3 + -. h_two i m a n s1 s3 *. f_two j a k l s2 s2 + -. h_two j m a l s2 s3 *. f_two i a k n s1 s3 ) + | Alfa, Alfa, Beta + | Beta, Beta, Alfa -> + sum mos_cabs (fun a -> + h_two i j a l s1 s2 *. f_two a m k n s1 s3 + +. h_two i m k a s1 s3 *. f_two j a l n s2 s3 + +. h_two j m a n s2 s3 *. f_two i a k l s1 s2 + +. h_two j m l a s2 s3 *. f_two i a k n s1 s3 + -. h_two i j a k s1 s2 *. f_two a m l n s1 s3 + -. h_two i m a n s1 s3 *. f_two j a k l s2 s2 + -. h_two i m l a s1 s3 *. f_two j a k n s2 s3 + -. h_two j m k a s2 s3 *. f_two i a l n s1 s3 + ) + | Alfa, Beta, Beta + | Beta, Alfa, Alfa -> + sum mos_cabs (fun a -> + h_two i j a l s1 s2 *. f_two a m k n s1 s3 + +. h_two i m a n s1 s3 *. f_two a j k l s1 s2 + +. h_two i m k a s1 s3 *. f_two j a l n s2 s3 + +. h_two j m a n s2 s3 *. f_two i a k l s1 s2 + -. h_two i j a n s1 s2 *. f_two a m k l s1 s2 + -. h_two i j k a s1 s2 *. f_two m a l n s2 s3 + -. h_two i m a l s1 s3 *. f_two a j k n s1 s3 + -. h_two j m a l s2 s3 *. f_two i a k n s1 s3 + ) + | Beta, Alfa, Beta + | Alfa, Beta, Alfa -> assert false (*TODO *) + in + match phase with + | Phase.Pos -> result + | Phase.Neg -> -. result + in let result = { simulation ; aux_basis ; - hf12 ; hf12_anti ; - hf12_single ; hf12_single_anti + f_0 ; f_1 ; f_2 ; f_3 } in Parallel.broadcast (lazy result)