From 61ba79ee59a8e2cd681a75f17394e7c83ddcecc2 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 4 Nov 2019 23:51:47 +0100 Subject: [PATCH] Added three-e contrbutions --- CI/CI.ml | 314 +++++++++++++++++++++++++++++----------- CI/CIMatrixElement.ml | 63 +++++++- CI/Excitation.ml | 17 +++ CI/F12CI.ml | 114 +++++++++++++-- MOBasis/HF12.ml | 14 +- Utils/FourIdxStorage.ml | 1 + 6 files changed, 419 insertions(+), 104 deletions(-) diff --git a/CI/CI.ml b/CI/CI.ml index 65d648f..171fa6d 100644 --- a/CI/CI.ml +++ b/CI/CI.ml @@ -33,17 +33,17 @@ let eigenvalues t = let (_,x) = eigensystem t in x let h_integrals mo_basis = - let one_e_ints = MOBasis.one_e_ints mo_basis - and two_e_ints = MOBasis.two_e_ints mo_basis + let one_e_ints = MOBasis.one_e_ints mo_basis + and two_e_ints = MOBasis.two_e_ints mo_basis in ( (fun i j _ -> one_e_ints.{i,j}), (fun i j k l s s' -> - if s' = Spin.other s then + 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) - ) ) + ), None ) @@ -231,7 +231,7 @@ let create_matrix_arbitrary f det_space = (* Create a matrix using the fact that the determinant space is made of the outer product of spindeterminants. *) -let create_matrix_spin f det_space = +let create_matrix_spin ?(nmax=2) f det_space = lazy ( let ndet = Ds.size det_space in let a, b = @@ -242,36 +242,73 @@ let create_matrix_spin f det_space = let n_beta = Array.length b in - (** Array of (list of singles, list of doubles) in the beta spin *) + (** Array of (list of singles, list of doubles, list of triples) in the beta spin *) let degree_bb = - Array.map (fun det_i -> - let deg = Spindeterminant.degree det_i in - let doubles = - Array.mapi (fun i det_j -> - let d = deg det_j in - if d < 3 then - Some (i,d,det_j) - else - None - ) b - |> Array.to_list - |> Util.list_some - in - let singles = - List.filter (fun (i,d,det_j) -> d < 2) doubles - |> List.map (fun (i,_,det_j) -> (i,det_j)) - in - let doubles = - List.map (fun (i,_,det_j) -> (i,det_j)) doubles - in - (singles, doubles) - ) b + match nmax with + | 2 -> begin + Array.map (fun det_i -> + let deg = Spindeterminant.degree det_i in + let triples = [] in + let doubles = + Array.mapi (fun i det_j -> + let d = deg det_j in + if d < 3 then + Some (i,d,det_j) + else + None + ) b + |> Array.to_list + |> Util.list_some + in + let singles = + List.filter (fun (i,d,det_j) -> d < 2) doubles + |> List.map (fun (i,_,det_j) -> (i,det_j)) + in + let doubles = + List.map (fun (i,_,det_j) -> (i,det_j)) doubles + in + (singles, doubles, triples) + ) b + end + | 3 -> begin + Array.map (fun det_i -> + let deg = Spindeterminant.degree det_i in + let triples = + Array.mapi (fun i det_j -> + let d = deg det_j in + if d < 4 then + Some (i,d,det_j) + else + None + ) b + |> Array.to_list + |> Util.list_some + in + let doubles = + List.filter (fun (i,d,det_j) -> d < 3) triples + in + 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 + let doubles = + List.map (fun (i,_,det_j) -> (i,det_j)) doubles + in + let singles = + List.map (fun (i,_,det_j) -> (i,det_j)) singles + in + (singles, doubles, triples) + ) b + end + | _ -> assert false in let a = Array.to_list a and b = Array.to_list b in - let task (i,i_alfa) = + let task = let result = Array.init n_beta (fun _ -> []) in @@ -283,61 +320,124 @@ let create_matrix_spin f det_space = result.(i) <- (j, x) :: result.(i) ; in - let j = ref 1 in - let deg_a = Spindeterminant.degree i_alfa in - List.iter (fun j_alfa -> - let degree_a = deg_a j_alfa in - begin - match degree_a with - | 2 -> - let i' = ref 0 in - List.iteri (fun ib i_beta -> - let ki = Determinant.of_spindeterminants i_alfa i_beta in - let kj = Determinant.of_spindeterminants j_alfa i_beta in - update !i' (ib + !j) 2 0 ki kj; - incr i'; - ) b; - | 1 -> - let i' = ref 0 in - List.iteri (fun ib i_beta -> - let ki = Determinant.of_spindeterminants i_alfa i_beta in - let singles, _ = degree_bb.(ib) in - List.iter (fun (j', j_beta) -> - let kj = Determinant.of_spindeterminants j_alfa j_beta in - update !i' (j' + !j) 1 (Determinant.degree_beta ki kj) ki kj - ) singles; - incr i'; - ) b; - | 0 -> - let i' = ref 0 in - List.iteri (fun ib i_beta -> - let ki = Determinant.of_spindeterminants i_alfa i_beta in - let _singles, doubles = degree_bb.(ib) in - List.iter (fun (j', j_beta) -> - let kj = Determinant.of_spindeterminants j_alfa j_beta in - update !i' (j' + !j) 0 (Determinant.degree_beta ki kj) ki kj - ) doubles; - incr i'; - ) b; - | _ -> (); + fun (i,i_alfa) -> + begin match nmax with + | 2 -> begin + let j = ref 1 in + + let deg_a = Spindeterminant.degree i_alfa in + List.iter (fun j_alfa -> + let degree_a = deg_a j_alfa in + begin + match degree_a with + | 2 -> + let i' = ref 0 in + List.iteri (fun ib i_beta -> + let ki = Determinant.of_spindeterminants i_alfa i_beta in + let kj = Determinant.of_spindeterminants j_alfa i_beta in + update !i' (ib + !j) 2 0 ki kj; + incr i'; + ) b; + | 1 -> + let i' = ref 0 in + List.iteri (fun ib i_beta -> + let ki = Determinant.of_spindeterminants i_alfa i_beta in + let singles, _doubles, _triples = degree_bb.(ib) in + List.iter (fun (j', j_beta) -> + let kj = Determinant.of_spindeterminants j_alfa j_beta in + update !i' (j' + !j) 1 (Determinant.degree_beta ki kj) ki kj + ) singles; + incr i'; + ) b; + | 0 -> + let i' = ref 0 in + List.iteri (fun ib i_beta -> + let ki = Determinant.of_spindeterminants i_alfa i_beta in + let _singles, doubles, _triples = degree_bb.(ib) in + List.iter (fun (j', j_beta) -> + let kj = Determinant.of_spindeterminants j_alfa j_beta in + update !i' (j' + !j) 0 (Determinant.degree_beta ki kj) ki kj + ) doubles; + incr i'; + ) b; + | _ -> (); + end; + j := !j + n_beta + ) a + end + | 3 -> begin + let j = ref 1 in + + let deg_a = Spindeterminant.degree i_alfa in + List.iter (fun j_alfa -> + let degree_a = deg_a j_alfa in + begin + match degree_a with + | 3 -> + let i' = ref 0 in + List.iteri (fun ib i_beta -> + let ki = Determinant.of_spindeterminants i_alfa i_beta in + let kj = Determinant.of_spindeterminants j_alfa i_beta in + update !i' (ib + !j) 3 0 ki kj; + incr i'; + ) b; + | 2 -> + let i' = ref 0 in + List.iteri (fun ib i_beta -> + let ki = Determinant.of_spindeterminants i_alfa i_beta in + let singles, _doubles, _triples = degree_bb.(ib) in + List.iter (fun (j', j_beta) -> + let kj = Determinant.of_spindeterminants j_alfa j_beta in + update !i' (j' + !j) 2 (Determinant.degree_beta ki kj) ki kj + ) singles; + incr i'; + ) b; + | 1 -> + let i' = ref 0 in + List.iteri (fun ib i_beta -> + let ki = Determinant.of_spindeterminants i_alfa i_beta in + let _singles, doubles, _triples = degree_bb.(ib) in + List.iter (fun (j', j_beta) -> + let kj = Determinant.of_spindeterminants j_alfa j_beta in + update !i' (j' + !j) 1 (Determinant.degree_beta ki kj) ki kj + ) doubles; + incr i'; + ) b; + | 0 -> + let i' = ref 0 in + List.iteri (fun ib i_beta -> + let ki = Determinant.of_spindeterminants i_alfa i_beta in + let _singles, _doubles, triples = degree_bb.(ib) in + List.iter (fun (j', j_beta) -> + let kj = Determinant.of_spindeterminants j_alfa j_beta in + update !i' (j' + !j) 0 (Determinant.degree_beta ki kj) ki kj + ) triples; + incr i'; + ) b; + | _ -> (); + end; + j := !j + n_beta + ) a + end + | _ -> assert false end; - j := !j + n_beta - ) a; - let r = - Array.map (fun l -> - List.rev l - |> Vector.sparse_of_assoc_list ndet - ) result - in (i,r) + let r = + Array.map (fun l -> + List.rev l + |> Vector.sparse_of_assoc_list ndet + ) result + in (i,r) in + let result = Array.init ndet (fun _ -> Vector.sparse_of_assoc_list ndet []) in + List.mapi (fun i i_alfa -> i*n_beta, i_alfa) a |> Stream.of_list - |> Farm.run ~ordered:false ~f:task + |> Farm.run ~ordered:false ~f:task |> Stream.iter (fun (k, r) -> Array.iteri (fun j r_j -> result.(k+j) <- r_j) r; Printf.eprintf "%8d / %8d\r%!" (k+Array.length r) ndet; @@ -349,7 +449,7 @@ let create_matrix_spin f det_space = (* Create a matrix using the fact that the determinant space is made of the outer product of spindeterminants. *) -let create_matrix_spin_computed f det_space = +let create_matrix_spin_computed ?(nmax=2) f det_space = lazy ( let ndet = Ds.size det_space in let a, b = @@ -359,7 +459,7 @@ let create_matrix_spin_computed f det_space = in let n_beta = Array.length b in - let h i_alfa j_alfa = + let h2 i_alfa j_alfa = let deg_a = Spindeterminant.degree a.(i_alfa) a.(j_alfa) in match deg_a with | 2 -> @@ -397,11 +497,63 @@ let create_matrix_spin_computed f det_space = | _ -> (fun _ _ -> 0.) in + let h3 i_alfa j_alfa = + let deg_a = Spindeterminant.degree a.(i_alfa) a.(j_alfa) in + match deg_a with + | 3 -> + let ai, aj = a.(i_alfa), a.(j_alfa) in + (fun i_beta j_beta -> + if i_beta <> j_beta then 0. else + let deg_b = 0 in + let ki = Determinant.of_spindeterminants ai b.(i_beta) in + let kj = Determinant.of_spindeterminants aj b.(j_beta) in + f deg_a deg_b ki kj + ) + | 2 -> + let ai, aj = a.(i_alfa), a.(j_alfa) in + (fun i_beta j_beta -> + let deg_b = Spindeterminant.degree b.(i_beta) b.(j_beta) in + match deg_b with + | 0 | 1 -> + let ki = Determinant.of_spindeterminants ai b.(i_beta) in + let kj = Determinant.of_spindeterminants aj b.(j_beta) in + f deg_a deg_b ki kj + | _ -> 0. + ) + | 1 -> + let ai, aj = a.(i_alfa), a.(j_alfa) in + (fun i_beta j_beta -> + let deg_b = Spindeterminant.degree b.(i_beta) b.(j_beta) in + match deg_b with + | 0 | 1 | 2 -> + let deg_b = Spindeterminant.degree b.(i_beta) b.(j_beta) in + let ki = Determinant.of_spindeterminants ai b.(i_beta) in + let kj = Determinant.of_spindeterminants aj b.(j_beta) in + f deg_a deg_b ki kj + | _ -> 0. + ) + | 0 -> + let ai, aj = a.(i_alfa), a.(j_alfa) in + (fun i_beta j_beta -> + let deg_b = Spindeterminant.degree b.(i_beta) b.(j_beta) in + match deg_b with + | 0 | 1 | 2 | 3 -> + let deg_b = Spindeterminant.degree b.(i_beta) b.(j_beta) in + let ki = Determinant.of_spindeterminants ai b.(i_beta) in + let kj = Determinant.of_spindeterminants aj b.(j_beta) in + f deg_a deg_b ki kj + | _ -> 0. + ) + | _ -> (fun _ _ -> 0.) + in + let i_prev = ref (-10) and result = ref (fun _ -> 0.) in let h123 = ref (fun _ -> 0.) in + let h = if nmax = 2 then h2 else h3 in + let g i = (* i : index of the i-th determinant. 1 <= i <= ndet @@ -478,9 +630,9 @@ let make ?(n_states=1) ?(algo=`Direct) det_space = | Ds.Arbitrary _ -> create_matrix_arbitrary | Ds.Spin _ -> if algo = `Direct then - create_matrix_spin_computed + create_matrix_spin_computed ~nmax:2 else - create_matrix_spin + create_matrix_spin ~nmax:2 in f (fun deg_a deg_b ki kj -> if deg_a + deg_b > 0 then @@ -494,7 +646,7 @@ let make ?(n_states=1) ?(algo=`Direct) det_space = let f = match Ds.determinants det_space with | Ds.Arbitrary _ -> create_matrix_arbitrary - | Ds.Spin _ -> create_matrix_spin + | Ds.Spin _ -> create_matrix_spin ~nmax:2 in f (fun _deg_a _deg_b ki kj -> CIMatrixElement.make_s2 ki kj) det_space (*TODO*) @@ -827,7 +979,7 @@ let _pt2_en ci = let i_o1_alfa = h_ij mo_basis in let w_alfa det_i = - let one_e, two_e = h_integrals mo_basis in + 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 diff --git a/CI/CIMatrixElement.ml b/CI/CIMatrixElement.ml index cac190c..0822a66 100644 --- a/CI/CIMatrixElement.ml +++ b/CI/CIMatrixElement.ml @@ -59,7 +59,7 @@ let non_zero integrals degree_a degree_b ki kj = one +. two in - let result = + let result_2e = lazy ( match degree_a, degree_b with | 1, 1 -> (* alpha-beta double *) begin @@ -108,8 +108,65 @@ let non_zero integrals degree_a degree_b ki kj = diag_element | _ -> assert false - in - List.map (fun (one_e, two_e) -> result one_e two_e) integrals + ) in + + let result_3e = lazy ( + 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 + 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 + + | 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 *) + begin + let h, p, phase = Ex.single_of_spindet kia kja in + match phase with + | Phase.Pos -> fun one_e two_e _ -> single h p Spin.Alfa kia kib one_e two_e + | Phase.Neg -> fun one_e two_e _ -> -. single h p Spin.Alfa kia kib one_e two_e + end + + | 0, 1 -> (* beta single *) + begin + let h, p, phase = Ex.single_of_spindet kib kjb in + match phase with + | Phase.Pos -> fun one_e two_e _ -> single h p Spin.Beta kib kia one_e two_e + | Phase.Neg -> fun one_e two_e _ -> -. single h p Spin.Beta kib kia one_e two_e + end + + | 0, 0 -> (* diagonal element *) + fun one_e two_e _ -> diag_element one_e two_e + + | _ -> fun _ _ _ -> 0. + ) in + + List.map (fun (one_e, two_e, x) -> + match x with + | None -> (Lazy.force result_2e) one_e two_e + | Some three_e -> (Lazy.force result_3e) one_e two_e three_e + ) integrals + let make integrals ki kj = diff --git a/CI/Excitation.ml b/CI/Excitation.ml index 6b154da..09b40da 100644 --- a/CI/Excitation.ml +++ b/CI/Excitation.ml @@ -9,6 +9,7 @@ type t = | Identity of Phase.t | Single of Phase.t * single_exc | Double of Phase.t * single_exc * single_exc +| Triple of Phase.t * single_exc * single_exc * single_exc | Multiple of Phase.t * single_exc list @@ -77,11 +78,19 @@ let multiple_of_spindet t t' = in (phase, List.map2 (fun hole particle -> (hole, particle)) holes (List.rev particles) ) + let double_of_spindet t t' = match multiple_of_spindet t t' with | (phase, (h1,p1)::(h2,p2)::[]) -> (h1, p1, h2, p2, phase) | _ -> invalid_arg "t and t' are not doubly excited" + +let triple_of_spindet t t' = + match multiple_of_spindet t t' with + | (phase, (h1,p1)::(h2,p2)::(h3,p3)::[]) -> (h1, p1, h2, p2, h3, p3, phase) + | _ -> invalid_arg "t and t' are not doubly excited" + + let multiple_of_det t t' = let pa, a = multiple_of_spindet (Determinant.alfa t) (Determinant.alfa t') @@ -100,6 +109,12 @@ let double_of_det t t' = | _ -> assert false +let triple_of_det t t' = + match multiple_of_det t t' with + | Multiple (phase, [e1 ; e2 ; e3]) -> Triple (phase, e1, e2, e3) + | _ -> assert false + + let of_det t t' = match Determinant.degree t t' with | 0 -> if Determinant.phase t = Determinant.phase t' then @@ -108,6 +123,7 @@ let of_det t t' = Identity Phase.Neg | 1 -> single_of_det t t' | 2 -> double_of_det t t' + | 3 -> triple_of_det t t' | _ -> multiple_of_det t t' let pp_s_exc ppf t = @@ -123,6 +139,7 @@ let pp_exc ppf t = | Identity p -> p, [] | Single (p,x) -> p, x::[] | Double (p,x,y) -> p, x::y::[] + | Triple (p,x,y,z) -> p, x::y::z::[] | Multiple (p,l) -> p, l in Format.fprintf ppf "@[%c" diff --git a/CI/F12CI.ml b/CI/F12CI.ml index d3c8430..701ee6d 100644 --- a/CI/F12CI.ml +++ b/CI/F12CI.ml @@ -20,18 +20,68 @@ let mo_class t = Ds.mo_class @@ det_space t let eigensystem t = Lazy.force t.eigensystem +(* -let hf_ij_non_zero hf12_integrals deg_a deg_b ki kj = - let integrals = [ - let one_e _ _ _ = 0. in +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 kia = De.alfa ki and kib = De.beta ki - and kja = De.alfa kj and kjb = De.beta kj + 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 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 @@ -39,27 +89,59 @@ let hf_ij_non_zero hf12_integrals deg_a deg_b ki kj = 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' = Spin.other s then - hf12.{i,j,k,l} -. 1. *. ( - (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) - ) - else - hf12_anti.{i,j,k,l} -. 1. *. ( + 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 - (one_e, two_e) + + 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' = + 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 = if Parallel.master then @@ -73,7 +155,7 @@ let dressing_vector ~frozen_core hf12_integrals f12_amplitudes ci = let f = match Ds.determinants det_space with | Ds.Arbitrary _ -> CI.create_matrix_arbitrary - | Ds.Spin _ -> CI.create_matrix_spin_computed + | 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 diff --git a/MOBasis/HF12.ml b/MOBasis/HF12.ml index df8d77a..30794c3 100644 --- a/MOBasis/HF12.ml +++ b/MOBasis/HF12.ml @@ -108,11 +108,17 @@ let make ~simulation ~mo_basis ~aux_basis_filename () = let task (j,l) = let h i a b = - h_s.{i, a, b} <- ERI.get_phys eri i j a b -. ERI.get_phys eri i j b a ; - h_o.{i, a, b} <- ERI.get_phys eri i j 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 = - f_s.{a, b, k} <- 0.25 *. (F12.get_phys f12 a b k l -. F12.get_phys f12 a b l k) ; - f_o.{a, b, k} <- 0.375 *. F12.get_phys f12 a b k l +. 0.125 *. F12.get_phys f12 b a k l + 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 diff --git a/Utils/FourIdxStorage.ml b/Utils/FourIdxStorage.ml index ebf2f4c..d224c98 100644 --- a/Utils/FourIdxStorage.ml +++ b/Utils/FourIdxStorage.ml @@ -504,6 +504,7 @@ let four_index_transform_dense_sparse ds coef source = Array.of_list !result in + let n = ref 0 in Stream.of_list range_mo |> Farm.run ~f:task ~ordered:true ~comm:Parallel.Node.comm