LCPQ/QCaml
10
1
Fork 0
mirror of https://gitlab.com/scemama/QCaml.git synced 2024-11-07 14:43:41 +01:00
Code Releases Activity
3991f7e212
QCaml/MOBasis/HF12.ml

911 lines
25 KiB
OCaml
Raw Blame History

(** %{ $ \langle ij | H F | kl \rangle $ %} integrals. *)
open Lacaml.D
module Fis = FourIdxStorage
type t = {
simulation : Simulation.t ;
aux_basis : MOBasis.t ;
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 fx =
let f k =
let x =
Bigarray.(Array2.create Float64 fortran_layout d1 d2)
in
for j=1 to d2 do
for i=1 to d1 do
x.{i,j} <- fx i j k
done
done;
(k,x)
in
let result =
SharedMemory.create Bigarray.Float64 [| d1 ; d2 ; d3 |]
|> Bigarray.array3_of_genarray
in
Util.list_range 1 d3
|> Stream.of_list
|> Farm.run ~f
|> Stream.iter (fun (k,x) ->
for j=1 to d2 do
for i=1 to d1 do
result.{i,j,k} <- x.{i,j}
done
done)
;
Parallel.broadcast (lazy result)
let array_4_init d1 d2 d3 d4 fx =
let f l =
let x =
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
x.{i,j,k} <- fx i j k l
done
done
done;
(l,x)
in
let result =
SharedMemory.create Bigarray.Float64 [| d1;d2;d3;d4 |]
in
Util.list_range 1 d4
|> Stream.of_list
|> Farm.run ~f
|> Stream.iter (fun (l,x) ->
for k=1 to d3 do
for j=1 to d2 do
for i=1 to d1 do
result.{i,j,k,l} <- x.{i,j,k}
done
done
done)
;
Parallel.broadcast (lazy result)
let array_5_init d1 d2 d3 d4 d5 fx =
let f m =
let x =
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
x.{i,j,k,l} <- fx i j k l m
done
done
done
done;
(m,x)
in
let result =
SharedMemory.create Bigarray.Float64 [| d1;d2;d3;d4;d5 |]
in
Util.list_range 1 d5
|> Stream.of_list
|> Farm.run ~f
|> Stream.iter (fun (m,x) ->
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} <- x.{i,j,k,l}
done
done
done
done)
;
Parallel.broadcast (lazy result)
let make ~simulation ~mo_basis ~aux_basis_filename () =
let f12 = Util.of_some @@ Simulation.f12 simulation in
let mo_num = MOBasis.size mo_basis in
(* Add auxiliary basis set *)
let simulation =
let charge = Charge.to_int @@ Simulation.charge simulation
and multiplicity = Electrons.multiplicity @@ Simulation.electrons simulation
and nuclei = Simulation.nuclei simulation
in
let general_basis =
Basis.general_basis @@ Simulation.basis simulation
in
GeneralBasis.combine [
general_basis ; GeneralBasis.read aux_basis_filename
]
|> Basis.of_nuclei_and_general_basis nuclei
|> Simulation.make ~f12 ~charge ~multiplicity ~nuclei
in
let aux_basis =
MOBasis.of_mo_basis simulation mo_basis
in
let aux_num = MOBasis.size aux_basis in
(* Fire calculation of F12 and ERI *)
ignore @@ MOBasis.f12_ints aux_basis ;
ignore @@ MOBasis.two_e_ints aux_basis ;
(* Compute the <ij|QHF|kl> integrals *)
if Parallel.master then Printf.eprintf "Computing HF12 integrals\n%!";
let mos_cabs =
Util.list_range (mo_num+1) aux_num
in
let mos_in =
Util.list_range 1 mo_num
in
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
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
let f_one = fun _ _ _ -> 0. in
(* 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 ;
f_0 ; f_1 ; f_2 ; f_3
}
in
Parallel.broadcast (lazy result)
View Git Blame Copy Permalink