LCPQ/QCaml
10
1
Fork 0
mirror of https://gitlab.com/scemama/QCaml.git synced 2025-01-03 01:55:40 +01:00
Code Releases Activity

Added three-e contrbutions

Browse Source
This commit is contained in:
Anthony Scemama 2019-11-04 23:51:47 +01:00
parent d6dd38d706
commit 61ba79ee59
6 changed files with 419 additions and 104 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

314
CI/CI.ml
View File

@ -33,17 +33,17 @@ let eigenvalues t =
let (_,x) = eigensystem t in x let (_,x) = eigensystem t in x
let h_integrals mo_basis = let h_integrals mo_basis =
let one_e_ints = MOBasis.one_e_ints mo_basis let one_e_ints = MOBasis.one_e_ints mo_basis
and two_e_ints = MOBasis.two_e_ints mo_basis and two_e_ints = MOBasis.two_e_ints mo_basis
in in
( (fun i j _ -> one_e_ints.{i,j}), ( (fun i j _ -> one_e_ints.{i,j}),
(fun i j k l s s' -> (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 ERI.get_phys two_e_ints i j k l
else else
(ERI.get_phys two_e_ints i j k l) -. (ERI.get_phys two_e_ints i j k l) -.
(ERI.get_phys two_e_ints i j l k) (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 (* Create a matrix using the fact that the determinant space is made of
the outer product of spindeterminants. *) the outer product of spindeterminants. *)
let create_matrix_spin f det_space = let create_matrix_spin ?(nmax=2) f det_space =
lazy ( lazy (
let ndet = Ds.size det_space in let ndet = Ds.size det_space in
let a, b = let a, b =
@ -242,36 +242,73 @@ let create_matrix_spin f det_space =
let n_beta = Array.length b in 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 = let degree_bb =
Array.map (fun det_i -> match nmax with
let deg = Spindeterminant.degree det_i in | 2 -> begin
let doubles = Array.map (fun det_i ->
Array.mapi (fun i det_j -> let deg = Spindeterminant.degree det_i in
let d = deg det_j in let triples = [] in
if d < 3 then let doubles =
Some (i,d,det_j) Array.mapi (fun i det_j ->
else let d = deg det_j in
None if d < 3 then
) b Some (i,d,det_j)
|> Array.to_list else
|> Util.list_some None
in ) b
let singles = |> Array.to_list
List.filter (fun (i,d,det_j) -> d < 2) doubles |> Util.list_some
|> List.map (fun (i,_,det_j) -> (i,det_j)) in
in let singles =
let doubles = List.filter (fun (i,d,det_j) -> d < 2) doubles
List.map (fun (i,_,det_j) -> (i,det_j)) doubles |> List.map (fun (i,_,det_j) -> (i,det_j))
in in
(singles, doubles) let doubles =
) b 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 in
let a = Array.to_list a let a = Array.to_list a
and b = Array.to_list b and b = Array.to_list b
in in
let task (i,i_alfa) = let task =
let result = let result =
Array.init n_beta (fun _ -> []) Array.init n_beta (fun _ -> [])
in in
@ -283,61 +320,124 @@ let create_matrix_spin f det_space =
result.(i) <- (j, x) :: result.(i) ; result.(i) <- (j, x) :: result.(i) ;
in in
let j = ref 1 in fun (i,i_alfa) ->
let deg_a = Spindeterminant.degree i_alfa in begin match nmax with
List.iter (fun j_alfa -> | 2 -> begin
let degree_a = deg_a j_alfa in let j = ref 1 in
begin
match degree_a with let deg_a = Spindeterminant.degree i_alfa in
| 2 -> List.iter (fun j_alfa ->
let i' = ref 0 in let degree_a = deg_a j_alfa in
List.iteri (fun ib i_beta -> begin
let ki = Determinant.of_spindeterminants i_alfa i_beta in match degree_a with
let kj = Determinant.of_spindeterminants j_alfa i_beta in | 2 ->
update !i' (ib + !j) 2 0 ki kj; let i' = ref 0 in
incr i'; List.iteri (fun ib i_beta ->
) b; let ki = Determinant.of_spindeterminants i_alfa i_beta in
| 1 -> let kj = Determinant.of_spindeterminants j_alfa i_beta in
let i' = ref 0 in update !i' (ib + !j) 2 0 ki kj;
List.iteri (fun ib i_beta -> incr i';
let ki = Determinant.of_spindeterminants i_alfa i_beta in ) b;
let singles, _ = degree_bb.(ib) in | 1 ->
List.iter (fun (j', j_beta) -> let i' = ref 0 in
let kj = Determinant.of_spindeterminants j_alfa j_beta in List.iteri (fun ib i_beta ->
update !i' (j' + !j) 1 (Determinant.degree_beta ki kj) ki kj let ki = Determinant.of_spindeterminants i_alfa i_beta in
) singles; let singles, _doubles, _triples = degree_bb.(ib) in
incr i'; List.iter (fun (j', j_beta) ->
) b; let kj = Determinant.of_spindeterminants j_alfa j_beta in
| 0 -> update !i' (j' + !j) 1 (Determinant.degree_beta ki kj) ki kj
let i' = ref 0 in ) singles;
List.iteri (fun ib i_beta -> incr i';
let ki = Determinant.of_spindeterminants i_alfa i_beta in ) b;
let _singles, doubles = degree_bb.(ib) in | 0 ->
List.iter (fun (j', j_beta) -> let i' = ref 0 in
let kj = Determinant.of_spindeterminants j_alfa j_beta in List.iteri (fun ib i_beta ->
update !i' (j' + !j) 0 (Determinant.degree_beta ki kj) ki kj let ki = Determinant.of_spindeterminants i_alfa i_beta in
) doubles; let _singles, doubles, _triples = degree_bb.(ib) in
incr i'; List.iter (fun (j', j_beta) ->
) b; 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; end;
j := !j + n_beta let r =
) a; Array.map (fun l ->
let r = List.rev l
Array.map (fun l -> |> Vector.sparse_of_assoc_list ndet
List.rev l ) result
|> Vector.sparse_of_assoc_list ndet in (i,r)
) result
in (i,r)
in in
let result = let result =
Array.init ndet (fun _ -> Vector.sparse_of_assoc_list ndet []) Array.init ndet (fun _ -> Vector.sparse_of_assoc_list ndet [])
in in
List.mapi (fun i i_alfa -> i*n_beta, i_alfa) a List.mapi (fun i i_alfa -> i*n_beta, i_alfa) a
|> Stream.of_list |> Stream.of_list
|> Farm.run ~ordered:false ~f:task |> Farm.run ~ordered:false ~f:task
|> Stream.iter (fun (k, r) -> |> Stream.iter (fun (k, r) ->
Array.iteri (fun j r_j -> result.(k+j) <- r_j) r; Array.iteri (fun j r_j -> result.(k+j) <- r_j) r;
Printf.eprintf "%8d / %8d\r%!" (k+Array.length r) ndet; 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 (* Create a matrix using the fact that the determinant space is made of
the outer product of spindeterminants. *) the outer product of spindeterminants. *)
let create_matrix_spin_computed f det_space = let create_matrix_spin_computed ?(nmax=2) f det_space =
lazy ( lazy (
let ndet = Ds.size det_space in let ndet = Ds.size det_space in
let a, b = let a, b =
@ -359,7 +459,7 @@ let create_matrix_spin_computed f det_space =
in in
let n_beta = Array.length b 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 let deg_a = Spindeterminant.degree a.(i_alfa) a.(j_alfa) in
match deg_a with match deg_a with
| 2 -> | 2 ->
@ -397,11 +497,63 @@ let create_matrix_spin_computed f det_space =
| _ -> (fun _ _ -> 0.) | _ -> (fun _ _ -> 0.)
in 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) let i_prev = ref (-10)
and result = ref (fun _ -> 0.) and result = ref (fun _ -> 0.)
in in
let h123 = ref (fun _ -> 0.) in let h123 = ref (fun _ -> 0.) in
let h = if nmax = 2 then h2 else h3 in
let g i = let g i =
(* (*
i : index of the i-th determinant. 1 <= i <= ndet 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.Arbitrary _ -> create_matrix_arbitrary
| Ds.Spin _ -> | Ds.Spin _ ->
if algo = `Direct then if algo = `Direct then
create_matrix_spin_computed create_matrix_spin_computed ~nmax:2
else else
create_matrix_spin create_matrix_spin ~nmax:2
in in
f (fun deg_a deg_b ki kj -> f (fun deg_a deg_b ki kj ->
if deg_a + deg_b > 0 then if deg_a + deg_b > 0 then
@ -494,7 +646,7 @@ let make ?(n_states=1) ?(algo=`Direct) det_space =
let f = let f =
match Ds.determinants det_space with match Ds.determinants det_space with
| Ds.Arbitrary _ -> create_matrix_arbitrary | Ds.Arbitrary _ -> create_matrix_arbitrary
| Ds.Spin _ -> create_matrix_spin | Ds.Spin _ -> create_matrix_spin ~nmax:2
in in
f (fun _deg_a _deg_b ki kj -> CIMatrixElement.make_s2 ki kj) det_space f (fun _deg_a _deg_b ki kj -> CIMatrixElement.make_s2 ki kj) det_space
(*TODO*) (*TODO*)
@ -827,7 +979,7 @@ let _pt2_en ci =
let i_o1_alfa = h_ij mo_basis in let i_o1_alfa = h_ij mo_basis in
let w_alfa det_i = 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 = let fock_diag_alfa, fock_diag_beta =
Ds.fock_diag ci.det_space det_i Ds.fock_diag ci.det_space det_i
in in

63
CI/CIMatrixElement.ml
View File

@ -59,7 +59,7 @@ let non_zero integrals degree_a degree_b ki kj =
one +. two one +. two
in in
let result = let result_2e = lazy (
match degree_a, degree_b with match degree_a, degree_b with
| 1, 1 -> (* alpha-beta double *) | 1, 1 -> (* alpha-beta double *)
begin begin
@ -108,8 +108,65 @@ let non_zero integrals degree_a degree_b ki kj =
diag_element diag_element
| _ -> assert false | _ -> assert false
in ) in
List.map (fun (one_e, two_e) -> result one_e two_e) integrals
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 = let make integrals ki kj =

17
CI/Excitation.ml
View File

@ -9,6 +9,7 @@ type t =
| Identity of Phase.t | Identity of Phase.t
| Single of Phase.t * single_exc | Single of Phase.t * single_exc
| Double of Phase.t * single_exc * 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 | Multiple of Phase.t * single_exc list
@ -77,11 +78,19 @@ let multiple_of_spindet t t' =
in in
(phase, List.map2 (fun hole particle -> (hole, particle)) holes (List.rev particles) ) (phase, List.map2 (fun hole particle -> (hole, particle)) holes (List.rev particles) )
let double_of_spindet t t' = let double_of_spindet t t' =
match multiple_of_spindet t t' with match multiple_of_spindet t t' with
| (phase, (h1,p1)::(h2,p2)::[]) -> (h1, p1, h2, p2, phase) | (phase, (h1,p1)::(h2,p2)::[]) -> (h1, p1, h2, p2, phase)
| _ -> invalid_arg "t and t' are not doubly excited" | _ -> 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 multiple_of_det t t' =
let pa, a = let pa, a =
multiple_of_spindet (Determinant.alfa t) (Determinant.alfa t') multiple_of_spindet (Determinant.alfa t) (Determinant.alfa t')
@ -100,6 +109,12 @@ let double_of_det t t' =
| _ -> assert false | _ -> 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' = let of_det t t' =
match Determinant.degree t t' with match Determinant.degree t t' with
| 0 -> if Determinant.phase t = Determinant.phase t' then | 0 -> if Determinant.phase t = Determinant.phase t' then
@ -108,6 +123,7 @@ let of_det t t' =
Identity Phase.Neg Identity Phase.Neg
| 1 -> single_of_det t t' | 1 -> single_of_det t t'
| 2 -> double_of_det t t' | 2 -> double_of_det t t'
| 3 -> triple_of_det t t'
| _ -> multiple_of_det t t' | _ -> multiple_of_det t t'
let pp_s_exc ppf t = let pp_s_exc ppf t =
@ -123,6 +139,7 @@ let pp_exc ppf t =
| Identity p -> p, [] | Identity p -> p, []
| Single (p,x) -> p, x::[] | Single (p,x) -> p, x::[]
| Double (p,x,y) -> p, x::y::[] | Double (p,x,y) -> p, x::y::[]
| Triple (p,x,y,z) -> p, x::y::z::[]
| Multiple (p,l) -> p, l | Multiple (p,l) -> p, l
in in
Format.fprintf ppf "@[%c" Format.fprintf ppf "@[%c"

114
CI/F12CI.ml
View File

@ -20,18 +20,68 @@ let mo_class t = Ds.mo_class @@ det_space t
let eigensystem t = Lazy.force t.eigensystem let eigensystem t = Lazy.force t.eigensystem
(*
let hf_ij_non_zero hf12_integrals deg_a deg_b ki kj = let single_matrices hf12_integrals density =
let integrals = [ let nocc = Mat.dim1 density in
let one_e _ _ _ = 0. in let nvir = Mat.dim2 density in
let { HF12. let { HF12.
simulation ; aux_basis ; simulation ; aux_basis ;
hf12 ; hf12_anti ; hf12 ; hf12_anti ;
hf12_single ; hf12_single_anti } = hf12_integrals hf12_single ; hf12_single_anti } = hf12_integrals
in in
let kia = De.alfa ki and kib = De.beta ki let f12 = MOBasis.f12_ints aux_basis in
and kja = De.alfa kj and kjb = De.beta kj 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 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 = let mo_a =
Bitstring.logand (Sp.bitstring kia) (Sp.bitstring kja) Bitstring.logand (Sp.bitstring kia) (Sp.bitstring kja)
|> Bitstring.to_list |> 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.logand (Sp.bitstring kib) (Sp.bitstring kjb)
|> Bitstring.to_list |> Bitstring.to_list
in in
let one_e _ _ _ = 0. in
let two_e i j k l s s' = let two_e i j k l s s' =
if s' = Spin.other s then if s = s' then
hf12.{i,j,k,l} -. 1. *. ( hf12_anti.{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)
)
else
hf12_anti.{i,j,k,l} -. 1. *. (
(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,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) (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 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 in
CIMatrixElement.non_zero integrals deg_a deg_b ki kj CIMatrixElement.non_zero integrals deg_a deg_b ki kj
|> List.hd |> List.hd
let dressing_vector ~frozen_core hf12_integrals f12_amplitudes ci = let dressing_vector ~frozen_core hf12_integrals f12_amplitudes ci =
if Parallel.master then if Parallel.master then
@ -73,7 +155,7 @@ let dressing_vector ~frozen_core hf12_integrals f12_amplitudes ci =
let f = let f =
match Ds.determinants det_space with match Ds.determinants det_space with
| Ds.Arbitrary _ -> CI.create_matrix_arbitrary | Ds.Arbitrary _ -> CI.create_matrix_arbitrary
| Ds.Spin _ -> CI.create_matrix_spin_computed | Ds.Spin _ -> CI.create_matrix_spin_computed ~nmax:3
in in
f (fun deg_a deg_b ki kj -> f (fun deg_a deg_b ki kj ->
hf_ij_non_zero hf12_integrals deg_a deg_b ki kj hf_ij_non_zero hf12_integrals deg_a deg_b ki kj

14
MOBasis/HF12.ml
View File

@ -108,11 +108,17 @@ let make ~simulation ~mo_basis ~aux_basis_filename () =
let task (j,l) = let task (j,l) =
let h i a b = 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 ; let ijab = ERI.get_phys eri i j a b
h_o.{i, a, b} <- 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 = 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) ; let abkl = F12.get_phys f12 a b k l
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 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 in
for a=mo_num+1 to aux_num do for a=mo_num+1 to aux_num do

1
Utils/FourIdxStorage.ml
View File

@ -504,6 +504,7 @@ let four_index_transform_dense_sparse ds coef source =
Array.of_list !result Array.of_list !result
in in
let n = ref 0 in let n = ref 0 in
Stream.of_list range_mo Stream.of_list range_mo
|> Farm.run ~f:task ~ordered:true ~comm:Parallel.Node.comm |> Farm.run ~f:task ~ordered:true ~comm:Parallel.Node.comm