mirror of
https://gitlab.com/scemama/QCaml.git
synced 2024-09-02 01:33:40 +02:00
1171 lines
41 KiB
OCaml
1171 lines
41 KiB
OCaml
(** %{ $ \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 = Array.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 ~ordered:false
|
|
|> 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)
|
|
;
|
|
if Parallel.master then Printf.printf "Broadcast d3\n%!" ;
|
|
try Parallel.InterNode.broadcast (lazy result)
|
|
with Invalid_argument _ ->
|
|
begin
|
|
Printf.eprintf "Array too large.\n%!";
|
|
failwith "Array too large"
|
|
end
|
|
|
|
let array_4_init d1 d2 d3 d4 fx =
|
|
let f (k,l) =
|
|
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 l
|
|
done
|
|
done;
|
|
(k,l,x)
|
|
in
|
|
let result =
|
|
SharedMemory.create Bigarray.Float64 [| d1;d2;d3;d4 |]
|
|
in
|
|
Util.list_range 1 d4
|
|
|> List.rev_map (fun l ->
|
|
Util.list_range 1 d3
|
|
|> List.rev_map (fun k -> (k,l)) )
|
|
|> List.concat
|
|
|> List.rev
|
|
|> Stream.of_list
|
|
|> Farm.run ~f ~ordered:false
|
|
|> Stream.iter (fun (k,l,x) ->
|
|
for j=1 to d2 do
|
|
for i=1 to d1 do
|
|
result.{i,j,k,l} <- x.{i,j}
|
|
done
|
|
done)
|
|
;
|
|
if Parallel.master then Printf.printf "Broadcast d4\n%!" ;
|
|
try Parallel.InterNode.broadcast (lazy result)
|
|
with Invalid_argument _ ->
|
|
begin
|
|
Printf.eprintf "Array too large... splitting.\n%!";
|
|
let x =
|
|
Bigarray.(Array3.create Float64 fortran_layout d1 d2 d3)
|
|
in
|
|
for l=1 to d4 do
|
|
if Parallel.master then
|
|
begin
|
|
for k=1 to d3 do
|
|
for j=1 to d2 do
|
|
for i=1 to d1 do
|
|
x.{i,j,k} <- result.{i,j,k,l}
|
|
done
|
|
done
|
|
done;
|
|
ignore @@ Parallel.InterNode.broadcast (lazy x)
|
|
end
|
|
else
|
|
begin
|
|
ignore @@ Parallel.InterNode.broadcast (lazy 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
|
|
end
|
|
done;
|
|
result
|
|
end
|
|
|
|
let array_5_init d1 d2 d3 d4 d5 fx =
|
|
let f (l,m) =
|
|
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 m
|
|
done
|
|
done
|
|
done;
|
|
(l,m,x)
|
|
in
|
|
let result =
|
|
SharedMemory.create Bigarray.Float64 [| d1;d2;d3;d4;d5 |]
|
|
in
|
|
Util.list_range 1 d5
|
|
|> List.rev_map (fun m ->
|
|
Util.list_range 1 d4
|
|
|> List.rev_map (fun l -> (l,m)) )
|
|
|> List.concat
|
|
|> List.rev
|
|
|> Stream.of_list
|
|
|> Farm.run ~f ~ordered:false
|
|
|> Stream.iter (fun (l,m,x) ->
|
|
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}
|
|
done
|
|
done
|
|
done)
|
|
;
|
|
if Parallel.master then Printf.printf "Broadcast d5\n%!" ;
|
|
try
|
|
Parallel.InterNode.broadcast (lazy result)
|
|
with Invalid_argument _ ->
|
|
begin
|
|
Printf.eprintf "Array too large... splitting.\n%!";
|
|
let x =
|
|
Bigarray.(Genarray.create Float64 fortran_layout [| d1; d2; d3; d4 |])
|
|
in
|
|
for m=1 to d5 do
|
|
if Parallel.master then
|
|
begin
|
|
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} <- result.{i,j,k,l,m}
|
|
done
|
|
done
|
|
done
|
|
done;
|
|
ignore @@ Parallel.InterNode.broadcast (lazy x)
|
|
end
|
|
else
|
|
begin
|
|
ignore @@ Parallel.InterNode.broadcast (lazy 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
|
|
end
|
|
done;
|
|
result
|
|
end
|
|
|
|
|
|
let make ~frozen_core ~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
|
|
|
|
(* While in a sequential region, initiate the parallel
|
|
4-idx transformation to avoid nested parallel jobs
|
|
*)
|
|
ignore @@ MOBasis.two_e_ints aux_basis;
|
|
ignore @@ MOBasis.f12_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
|
|
|> Array.of_list
|
|
in
|
|
|
|
let n_core =
|
|
if frozen_core then
|
|
(Nuclei.small_core @@ Simulation.nuclei @@ MOBasis.simulation mo_basis) / 2
|
|
else 0
|
|
in
|
|
|
|
let mos_in =
|
|
Util.list_range (n_core+1) mo_num
|
|
|> Array.of_list
|
|
in
|
|
|
|
let mos_a k =
|
|
Determinant.alfa k
|
|
|> Spindeterminant.to_list
|
|
|> List.filter (fun i -> i > n_core)
|
|
|> Array.of_list
|
|
in
|
|
|
|
let mos_b k =
|
|
Determinant.beta k
|
|
|> Spindeterminant.to_list
|
|
|> List.filter (fun i -> i > n_core)
|
|
|> Array.of_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
|
|
|
|
|
|
(*
|
|
let s, s' = Spin.(Alfa, Beta) in
|
|
let m_3_Hschwarz_aa =
|
|
array_3_init mo_num mo_num mo_num (fun i j k ->
|
|
sum mos_cabs (fun a -> let x = h_two a i j k s s in x *. x)
|
|
)
|
|
in
|
|
let m_3_Hschwarz_ab =
|
|
array_3_init mo_num mo_num mo_num (fun i j k ->
|
|
sum mos_cabs (fun a -> let x = h_two a i j k s s' in x *. x)
|
|
)
|
|
in
|
|
let m_3_Fschwarz_aa =
|
|
array_3_init mo_num mo_num mo_num (fun i j k ->
|
|
sum mos_cabs (fun a -> let x = f_two a i j k s s in x *. x)
|
|
)
|
|
in
|
|
let m_3_Fschwarz_ab =
|
|
array_3_init mo_num mo_num mo_num (fun i j k ->
|
|
sum mos_cabs (fun a -> let x = f_two a i j k s s' in x *. x)
|
|
)
|
|
in
|
|
|
|
let eps = 1.e-15 in
|
|
|
|
let sum_schwarz i j k l m n s1 s2 s3 s4 =
|
|
if (s1 = s2 && m_3_Hschwarz_aa.{i,j,k} < eps ) ||
|
|
(s1 <> s2 && m_3_Hschwarz_ab.{i,j,k} < eps ) ||
|
|
(s3 = s4 && m_3_Fschwarz_aa.{l,m,n} < eps ) ||
|
|
(s3 <> s4 && m_3_Fschwarz_ab.{l,m,n} < eps )
|
|
then
|
|
0.
|
|
else
|
|
Array.fold_left (fun accu a ->
|
|
accu +. h_two a i j k s1 s2 *. f_two a l m n s3 s4
|
|
) 0. mos_cabs
|
|
in
|
|
*)
|
|
|
|
(* Pre-compute dressed integrals *)
|
|
if Parallel.master then Printf.printf "Computing m_0111_1H_1F\n%!" ;
|
|
let m_0111_1H_1F =
|
|
Vec.init mo_num (fun i ->
|
|
sum mos_cabs (fun a ->
|
|
h_one a i Spin.Alfa *. f_one a i Spin.Alfa ))
|
|
in
|
|
|
|
if Parallel.master then Printf.printf "Computing m_0111_1H_2Fa\n%!" ;
|
|
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 a k i j 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 a i Spin.Alfa *. f_two a j i j Spin.Alfa Spin.Alfa +.
|
|
h_two a j i 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 a j i 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
|
|
|
|
if Parallel.master then Printf.printf "Computing m_0111_1H_2Fb\n%!" ;
|
|
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 a k i j 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 a i Spin.Alfa *. f_two a j i j Spin.Alfa Spin.Beta +.
|
|
h_two a j i j Spin.Alfa Spin.Beta *. f_one a i Spin.Alfa +.
|
|
h_one a j Spin.Alfa *. f_two a i j i Spin.Alfa Spin.Beta +.
|
|
h_two a i j 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 b ->
|
|
sum mos_cabs (fun a ->
|
|
h_two a b i j 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 a i k i Spin.Alfa Spin.Beta *.
|
|
f_two a j k j Spin.Alfa Spin.Beta +.
|
|
h_two a k i 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
|
|
|
|
if Parallel.master then Printf.printf "Computing m_0111_2Ha_2Fb\n%!" ;
|
|
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 a j i 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
|
|
|
|
if Parallel.master then Printf.printf "Computing m_1111_1H_1F\n%!" ;
|
|
let m_1111_1H_1F =
|
|
Mat.init_cols mo_num mo_num (fun i k ->
|
|
sum mos_cabs (fun a -> h_one a i Spin.Alfa *. f_one a k Spin.Alfa ))
|
|
in
|
|
|
|
if Parallel.master then Printf.printf "Computing m_1111_2Ha_2Fa\n%!" ;
|
|
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 a l j l Spin.Alfa Spin.Alfa *.
|
|
f_two a i j k Spin.Alfa Spin.Alfa )
|
|
else
|
|
sum mos_cabs (fun a ->
|
|
h_two a j i j Spin.Alfa Spin.Alfa *.
|
|
f_two a l k l Spin.Alfa Spin.Alfa +.
|
|
h_two a l j l Spin.Alfa Spin.Alfa *.
|
|
f_two a i j k Spin.Alfa Spin.Alfa )
|
|
)
|
|
in
|
|
|
|
if Parallel.master then Printf.printf "Computing m_1111_2Hb_2Fa\n%!" ;
|
|
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 a l j l Spin.Alfa Spin.Beta *.
|
|
f_two a i j k Spin.Alfa Spin.Beta )
|
|
else
|
|
sum mos_cabs (fun a ->
|
|
h_two a j i j Spin.Alfa Spin.Beta *.
|
|
f_two a l k l Spin.Alfa Spin.Alfa +.
|
|
h_two a l j l Spin.Alfa Spin.Beta *.
|
|
f_two a i j k Spin.Alfa Spin.Beta )
|
|
)
|
|
in
|
|
|
|
if Parallel.master then Printf.printf "Computing m_1111_2Ha_2Fb\n%!" ;
|
|
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 a j i j Spin.Alfa Spin.Alfa *.
|
|
f_two a l k l Spin.Alfa Spin.Beta +.
|
|
h_two a l j l Spin.Alfa Spin.Beta *.
|
|
f_two a i j k Spin.Alfa Spin.Alfa
|
|
)
|
|
)
|
|
in
|
|
|
|
if Parallel.master then Printf.printf "Computing m_1111_2Hb_2Fb\n%!" ;
|
|
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 a j i j Spin.Alfa Spin.Beta *.
|
|
f_two a l k l Spin.Alfa Spin.Beta +.
|
|
h_two a l j l Spin.Alfa Spin.Alfa *.
|
|
f_two a i j k Spin.Alfa Spin.Beta
|
|
)
|
|
)
|
|
in
|
|
|
|
if Parallel.master then Printf.printf "Computing m_1121_2Ha_2Fa\n%!" ;
|
|
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 a k j i Spin.Alfa Spin.Alfa *.
|
|
f_two a l j l Spin.Alfa Spin.Alfa
|
|
)
|
|
)
|
|
in
|
|
|
|
if Parallel.master then Printf.printf "Computing m_1121_2Hb_2Fa\n%!" ;
|
|
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 a k j i Spin.Alfa Spin.Beta *.
|
|
f_two a l j l Spin.Alfa Spin.Alfa
|
|
)
|
|
)
|
|
in
|
|
|
|
if Parallel.master then Printf.printf "Computing m_1121_2Ha_2Fb\n%!" ;
|
|
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 a k j i Spin.Alfa Spin.Alfa *.
|
|
f_two a l j l Spin.Alfa Spin.Beta
|
|
)
|
|
)
|
|
in
|
|
|
|
if Parallel.master then Printf.printf "Computing m_1121_2Hb_2Fb\n%!" ;
|
|
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 a k j i Spin.Alfa Spin.Beta *.
|
|
f_two a l j l Spin.Alfa Spin.Beta
|
|
)
|
|
)
|
|
in
|
|
|
|
if Parallel.master then Printf.printf "Computing m_1122_va\n%!" ;
|
|
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 a l j i Spin.Alfa Spin.Alfa *.
|
|
f_two a l j k Spin.Alfa Spin.Alfa
|
|
)
|
|
)
|
|
in
|
|
|
|
if Parallel.master then Printf.printf "Computing m_1111_1H_2Fa\n%!" ;
|
|
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 a i Spin.Alfa *. f_two a j k j Spin.Alfa Spin.Alfa +.
|
|
h_two a j i j Spin.Alfa Spin.Alfa *. f_one a k Spin.Alfa +.
|
|
h_one a j Spin.Alfa *. f_two a i j k Spin.Alfa Spin.Alfa +.
|
|
h_two a k j i Spin.Alfa Spin.Alfa *. f_one a j Spin.Alfa +.
|
|
sum mos_cabs (fun b -> if b > a then 0. else
|
|
h_two b a j i Spin.Alfa Spin.Alfa *.
|
|
f_two b a j k Spin.Alfa Spin.Alfa )
|
|
)
|
|
)
|
|
in
|
|
|
|
if Parallel.master then Printf.printf "Computing m_1122_v2\n%!" ;
|
|
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 a l j i Spin.Alfa Spin.Beta *.
|
|
f_two a l j k Spin.Alfa Spin.Beta
|
|
)
|
|
)
|
|
in
|
|
|
|
if Parallel.master then Printf.printf "Computing m_1122_v3\n%!" ;
|
|
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 a l i j Spin.Alfa Spin.Beta *.
|
|
f_two a l k j Spin.Alfa Spin.Beta
|
|
)
|
|
)
|
|
in
|
|
|
|
if Parallel.master then Printf.printf "Computing m_1111_1H_2Fb\n%!" ;
|
|
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 a i Spin.Alfa *. f_two a j k j Spin.Alfa Spin.Beta +.
|
|
h_two a j i j Spin.Alfa Spin.Beta *. f_one a k Spin.Alfa +.
|
|
h_one a j Spin.Beta *. f_two a i j k Spin.Alfa Spin.Beta +.
|
|
h_two a k j i Spin.Alfa Spin.Beta *. f_one a j Spin.Beta +.
|
|
sum mos_cabs (fun b ->
|
|
h_two b a j i Spin.Alfa Spin.Beta *.
|
|
f_two b a j k Spin.Alfa Spin.Beta
|
|
)
|
|
)
|
|
)
|
|
in
|
|
|
|
if Parallel.master then Printf.printf "Computing m_1122_oa\n%!" ;
|
|
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 a k j l Spin.Alfa Spin.Alfa *.
|
|
f_two a i j l Spin.Alfa Spin.Alfa
|
|
)
|
|
else 0.
|
|
)
|
|
in
|
|
|
|
if Parallel.master then Printf.printf "Computing m_1122_o\n%!" ;
|
|
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 a k j l Spin.Alfa Spin.Beta *.
|
|
f_two a i j l 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)
|
|
|> Array.of_list
|
|
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)
|
|
|> Array.of_list
|
|
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)
|
|
|> Array.of_list
|
|
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)
|
|
|> Array.of_list
|
|
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
|
|
|
|
if Parallel.master then Printf.printf "Computing m_2112_1H_2Fa\n%!" ;
|
|
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 a i Spin.Alfa *. f_two a j k l Spin.Alfa Spin.Alfa +.
|
|
h_one a j Spin.Alfa *. f_two i a k l Spin.Alfa Spin.Alfa +.
|
|
h_two a l i j Spin.Alfa Spin.Alfa *. f_one a k Spin.Alfa +.
|
|
h_two a k j i Spin.Alfa Spin.Alfa *. f_one a l Spin.Alfa ) +.
|
|
sum mos_cabs (fun a ->
|
|
sum mos_in (fun m -> -. h_two m a j i Spin.Alfa Spin.Alfa *.
|
|
f_two m a k l Spin.Alfa Spin.Alfa) ) +.
|
|
sum mos_cabs (fun a ->
|
|
sum mos_cabs (fun b -> if b >= a then 0. else
|
|
h_two b a j i Spin.Alfa Spin.Alfa *. f_two b a l k Spin.Alfa Spin.Alfa
|
|
)
|
|
)
|
|
)
|
|
in
|
|
|
|
if Parallel.master then Printf.printf "Computing m_2112_1H_2Fb\n%!" ;
|
|
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 a i Spin.Alfa *. f_two a j k l Spin.Alfa Spin.Beta +.
|
|
h_one a j Spin.Alfa *. f_two a i l k Spin.Alfa Spin.Beta +.
|
|
h_two a l i j Spin.Alfa Spin.Beta *. f_one a k Spin.Alfa +.
|
|
h_two a k j i Spin.Alfa Spin.Beta *. f_one a l Spin.Alfa ) +.
|
|
sum mos_cabs (fun a ->
|
|
sum mos_in (fun m ->
|
|
h_two m a j i Spin.Alfa Spin.Beta *. f_two m a l k Spin.Alfa Spin.Beta +.
|
|
h_two m a i j Spin.Alfa Spin.Beta *. f_two m a k l Spin.Alfa Spin.Beta ) ) +.
|
|
sum mos_cabs (fun a ->
|
|
sum mos_cabs (fun b ->
|
|
h_two b a j i Spin.Alfa Spin.Beta *. f_two b a l k Spin.Alfa Spin.Beta
|
|
)
|
|
)
|
|
)
|
|
in
|
|
|
|
if Parallel.master then Printf.printf "Computing m_2112_2Ha_2Fa\n%!" ;
|
|
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 a n i n Spin.Alfa Spin.Alfa *. f_two a j k l Spin.Alfa Spin.Alfa +.
|
|
h_two a n j n Spin.Alfa Spin.Alfa *. f_two a i l k Spin.Alfa Spin.Alfa
|
|
)
|
|
)
|
|
in
|
|
|
|
if Parallel.master then Printf.printf "Computing m_2112_2Hb_2Fa\n%!" ;
|
|
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 a n i n Spin.Alfa Spin.Beta *. f_two a j k l Spin.Alfa Spin.Alfa +.
|
|
h_two a n j n Spin.Alfa Spin.Beta *. f_two a i l k Spin.Alfa Spin.Alfa
|
|
)
|
|
)
|
|
in
|
|
|
|
if Parallel.master then Printf.printf "Computing m_2112_2Ha_2Fb\n%!" ;
|
|
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 a n i n Spin.Alfa Spin.Alfa *. f_two a j k l Spin.Alfa Spin.Beta +.
|
|
h_two a n j n Spin.Alfa Spin.Beta *. f_two a i l k Spin.Alfa Spin.Beta
|
|
)
|
|
)
|
|
in
|
|
|
|
if Parallel.master then Printf.printf "Computing m_2121_2Ha_2Fa\n%!" ;
|
|
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 a l i j Spin.Alfa Spin.Alfa *. f_two a n k n Spin.Alfa Spin.Alfa +.
|
|
h_two a k j i Spin.Alfa Spin.Alfa *. f_two a n l n Spin.Alfa Spin.Alfa
|
|
)
|
|
)
|
|
in
|
|
|
|
if Parallel.master then Printf.printf "Computing m_2121_2Hb_2Fa\n%!" ;
|
|
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 a l i j Spin.Alfa Spin.Beta *. f_two a n k n Spin.Alfa Spin.Alfa +.
|
|
h_two a k j i Spin.Alfa Spin.Beta *. f_two a n l n Spin.Alfa Spin.Beta
|
|
)
|
|
)
|
|
in
|
|
|
|
if Parallel.master then Printf.printf "Computing m_2121_2Ha_2Fb\n%!" ;
|
|
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 a l i j Spin.Alfa Spin.Alfa *. f_two a n k n Spin.Alfa Spin.Beta +.
|
|
h_two a k j i Spin.Alfa Spin.Alfa *. f_two a n l n Spin.Alfa Spin.Beta
|
|
)
|
|
)
|
|
in
|
|
|
|
if Parallel.master then Printf.printf "Computing m_2122_2Ha_2Fa_ij\n%!" ;
|
|
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 a k i n s s *. f_two a j l n s s
|
|
+. h_two a l i n s s *. f_two a j n k s s
|
|
-. h_two a k j n s s *. f_two a i l n s s
|
|
-. h_two a l j n s s *. f_two a i n k s s
|
|
)
|
|
)
|
|
in
|
|
|
|
if Parallel.master then Printf.printf "Computing m_2122_2Hb_2Fb_ij\n%!" ;
|
|
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 a k n i s s *. f_two a j n l s s'
|
|
+. h_two a l n j s s' *. f_two a i n k s s
|
|
-. h_two a k j n s s' *. f_two a i l n s s'
|
|
)
|
|
)
|
|
in
|
|
|
|
if Parallel.master then Printf.printf "Computing m_2122_2Hb_2Fb_ij2\n%!" ;
|
|
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 a l i n s s' *. f_two a j k n s s' +.
|
|
(if n < j then
|
|
h_two a k n i s s' *. f_two a j n l s' s'
|
|
+. h_two a l n j s' s' *. f_two a i n k s s'
|
|
else
|
|
-. h_two a k n i s s' *. f_two j a n l s' s'
|
|
-. h_two a l j n s' s' *. f_two i a k n s s'
|
|
)
|
|
)
|
|
)
|
|
in
|
|
|
|
if Parallel.master then Printf.printf "Computing m_2122_2Ha_2Fa_ij2\n%!" ;
|
|
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 a k n i s s' *. f_two a j n l s s') +.
|
|
sum mos_cabs (fun a -> h_two a l n j s s' *. f_two a i n k s s') -.
|
|
sum mos_cabs (fun a -> h_two a l n i s s' *. f_two a j n k s s') -.
|
|
sum mos_cabs (fun a -> h_two a k n j s s' *. f_two a i n l s s')
|
|
)
|
|
in
|
|
|
|
if Parallel.master then Printf.printf "Computing m_2122_2Ha_2Fa_nv\n%!" ;
|
|
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 a n i j s s *. f_two a n l k s s ) )
|
|
in
|
|
|
|
if Parallel.master then Printf.printf "Computing m_2122_2Hb_2Fb_nv\n%!" ;
|
|
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 a n i j s s' *. f_two a n k l s s' ) )
|
|
in
|
|
|
|
if Parallel.master then Printf.printf "Computing m_2122_2Hb_2Fb_nv2\n%!" ;
|
|
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 a n j i s s' *. f_two a n l k 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)
|
|
|> Array.of_list
|
|
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)
|
|
|> Array.of_list
|
|
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)
|
|
|> Array.of_list
|
|
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)
|
|
|> Array.of_list
|
|
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 ncabs = Array.length mos_cabs in
|
|
let tmp_h = Mat.make0 ncabs 9 in
|
|
let tmp_f = Mat.make0 ncabs 9 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 a k i j s1 s2 *. f_two a m n l s3 s3)
|
|
+. sum mos_cabs (fun a -> h_two a n i j s1 s2 *. f_two a m l k s2 s3)
|
|
-. sum mos_cabs (fun a -> h_two a l i j s1 s2 *. f_two a m n k s3 s3)
|
|
+. sum mos_cabs (fun a -> h_two a l i m s1 s3 *. f_two a j n k s3 s2)
|
|
-. sum mos_cabs (fun a -> h_two a n i m s1 s3 *. f_two a j l k s2 s2)
|
|
-. sum mos_cabs (fun a -> h_two a k i m s1 s3 *. f_two a j n l s3 s2)
|
|
+. sum mos_cabs (fun a -> h_two a n j m s2 s3 *. f_two a i l k s2 s1)
|
|
+. sum mos_cabs (fun a -> h_two a k j m s2 s3 *. f_two a i n l s3 s1)
|
|
-. sum mos_cabs (fun a -> h_two a l j m s2 s3 *. f_two a i n k s3 s1)
|
|
|
|
sum_schwarz k i j m n l s1 s2 s3 s3
|
|
+. sum_schwarz n i j m l k s1 s2 s2 s3
|
|
-. sum_schwarz l i j m n k s1 s2 s3 s3
|
|
+. sum_schwarz l i m j n k s1 s3 s3 s2
|
|
-. sum_schwarz n i m j l k s1 s3 s2 s2
|
|
-. sum_schwarz k i m j n l s1 s3 s3 s2
|
|
+. sum_schwarz n j m i l k s2 s3 s2 s1
|
|
+. sum_schwarz k j m i n l s2 s3 s3 s1
|
|
-. sum_schwarz l j m i n k s2 s3 s3 s1
|
|
*)
|
|
Array.iteri (fun idx a -> tmp_f.{idx+1,7} <- f_two a i l k s2 s1) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_f.{idx+1,5} <- f_two a j l k s2 s2) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_f.{idx+1,2} <- f_two a m l k s2 s3) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_f.{idx+1,9} <- f_two a i n k s3 s1) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_f.{idx+1,4} <- f_two a j n k s3 s2) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_f.{idx+1,3} <- f_two a m n k s3 s3) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_f.{idx+1,8} <- f_two a i n l s3 s1) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_f.{idx+1,6} <- f_two a j n l s3 s2) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_f.{idx+1,1} <- f_two a m n l s3 s3) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_h.{idx+1,3} <- if tmp_f.{idx+1,3} = 0. then 0. else -. h_two a l i j s1 s2) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_h.{idx+1,1} <- if tmp_f.{idx+1,1} = 0. then 0. else h_two a k i j s1 s2) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_h.{idx+1,2} <- if tmp_f.{idx+1,2} = 0. then 0. else h_two a n i j s1 s2) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_h.{idx+1,5} <- if tmp_f.{idx+1,5} = 0. then 0. else -. h_two a n i m s1 s3) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_h.{idx+1,6} <- if tmp_f.{idx+1,6} = 0. then 0. else -. h_two a k i m s1 s3) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_h.{idx+1,4} <- if tmp_f.{idx+1,4} = 0. then 0. else h_two a l i m s1 s3) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_h.{idx+1,9} <- if tmp_f.{idx+1,9} = 0. then 0. else -. h_two a l j m s2 s3) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_h.{idx+1,7} <- if tmp_f.{idx+1,7} = 0. then 0. else h_two a n j m s2 s3) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_h.{idx+1,8} <- if tmp_f.{idx+1,8} = 0. then 0. else h_two a k j m s2 s3) mos_cabs;
|
|
dot (Mat.as_vec tmp_h) (Mat.as_vec tmp_f)
|
|
| Alfa, Alfa, Beta
|
|
| Beta, Beta, Alfa ->
|
|
(*
|
|
sum mos_cabs (fun a -> h_two a l i j s1 s2 *. f_two a m k n s1 s3)
|
|
-. sum mos_cabs (fun a -> h_two a k i j s1 s2 *. f_two a m l n s1 s3)
|
|
+. sum mos_cabs (fun a -> h_two a l m j s3 s2 *. f_two a i n k s3 s1)
|
|
-. sum mos_cabs (fun a -> h_two a k m j s3 s2 *. f_two a i n l s3 s1)
|
|
+. sum mos_cabs (fun a -> h_two a k m i s3 s1 *. f_two a j n l s3 s2)
|
|
-. sum mos_cabs (fun a -> h_two a l m i s3 s1 *. f_two a j n k s3 s2)
|
|
+. sum mos_cabs (fun a -> h_two a n j m s2 s3 *. f_two a i l k s2 s1)
|
|
-. sum mos_cabs (fun a -> h_two a n i m s1 s3 *. f_two a j l k s2 s2)
|
|
|
|
sum_schwarz l i j m k n s1 s2 s1 s3
|
|
-. sum_schwarz k i j m l n s1 s2 s1 s3
|
|
+. sum_schwarz l m j i n k s3 s2 s3 s1
|
|
-. sum_schwarz k m j i n l s3 s2 s3 s1
|
|
+. sum_schwarz k m i j n l s3 s1 s3 s2
|
|
-. sum_schwarz l m i j n k s3 s1 s3 s2
|
|
+. sum_schwarz n j m i l k s2 s3 s2 s1
|
|
-. sum_schwarz n i m j l k s1 s3 s2 s2
|
|
*)
|
|
Array.iteri (fun idx a -> tmp_f.{idx+1,1} <- f_two a m k n s1 s3) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_f.{idx+1,2} <- f_two a m l n s1 s3) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_f.{idx+1,4} <- f_two a i n l s3 s1) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_f.{idx+1,5} <- f_two a j n l s3 s2) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_f.{idx+1,3} <- f_two a i n k s3 s1) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_f.{idx+1,6} <- f_two a j n k s3 s2) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_f.{idx+1,7} <- f_two a i l k s2 s1) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_f.{idx+1,8} <- f_two a j l k s2 s2) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_h.{idx+1,1} <- if tmp_f.{idx+1,1} = 0. then 0. else h_two a l i j s1 s2) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_h.{idx+1,2} <- if tmp_f.{idx+1,2} = 0. then 0. else -. h_two a k i j s1 s2) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_h.{idx+1,3} <- if tmp_f.{idx+1,3} = 0. then 0. else +. h_two a l m j s3 s2) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_h.{idx+1,4} <- if tmp_f.{idx+1,4} = 0. then 0. else -. h_two a k m j s3 s2) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_h.{idx+1,5} <- if tmp_f.{idx+1,5} = 0. then 0. else +. h_two a k m i s3 s1) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_h.{idx+1,6} <- if tmp_f.{idx+1,6} = 0. then 0. else -. h_two a l m i s3 s1) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_h.{idx+1,7} <- if tmp_f.{idx+1,7} = 0. then 0. else +. h_two a n j m s2 s3) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_h.{idx+1,8} <- if tmp_f.{idx+1,8} = 0. then 0. else -. h_two a n i m s1 s3) mos_cabs;
|
|
dot ~n:(8*ncabs) (Mat.as_vec tmp_h) (Mat.as_vec tmp_f)
|
|
| Alfa, Beta, Beta
|
|
| Beta, Alfa, Alfa ->
|
|
(*
|
|
sum mos_cabs (fun a -> h_two a l i j s1 s2 *. f_two a m k n s1 s3)
|
|
-. sum mos_cabs (fun a -> h_two a n i j s1 s2 *. f_two a m k l s1 s2)
|
|
+. sum mos_cabs (fun a -> h_two a n i m s1 s3 *. f_two a j k l s1 s2)
|
|
+. sum mos_cabs (fun a -> h_two a n j m s2 s3 *. f_two a i l k s2 s1)
|
|
-. sum mos_cabs (fun a -> h_two a l i m s1 s3 *. f_two a j k n s1 s3)
|
|
-. sum mos_cabs (fun a -> h_two a l j m s2 s3 *. f_two a i n k s3 s1)
|
|
+. sum mos_cabs (fun a -> h_two a k m i s3 s1 *. f_two a j n l s3 s2)
|
|
-. sum mos_cabs (fun a -> h_two a k j i s2 s1 *. f_two a m n l s3 s2)
|
|
|
|
sum_schwarz l i j m k n s1 s2 s1 s3
|
|
-. sum_schwarz n i j m k l s1 s2 s1 s2
|
|
+. sum_schwarz n i m j k l s1 s3 s1 s2
|
|
+. sum_schwarz n j m i l k s2 s3 s2 s1
|
|
-. sum_schwarz l i m j k n s1 s3 s1 s3
|
|
-. sum_schwarz l j m i n k s2 s3 s3 s1
|
|
+. sum_schwarz k m i j n l s3 s1 s3 s2
|
|
-. sum_schwarz k j i m n l s2 s1 s3 s2
|
|
*)
|
|
Array.iteri (fun idx a -> tmp_f.{idx+1,6} <- f_two a i n k s3 s1) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_f.{idx+1,4} <- f_two a i l k s2 s1) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_f.{idx+1,5} <- f_two a j k n s1 s3) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_f.{idx+1,1} <- f_two a m k n s1 s3) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_f.{idx+1,2} <- f_two a m k l s1 s2) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_f.{idx+1,3} <- f_two a j k l s1 s2) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_f.{idx+1,7} <- f_two a j n l s3 s2) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_f.{idx+1,8} <- f_two a m n l s3 s2) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_h.{idx+1,1} <- if tmp_f.{idx+1,1} = 0. then 0. else h_two a l i j s1 s2) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_h.{idx+1,2} <- if tmp_f.{idx+1,2} = 0. then 0. else -. h_two a n i j s1 s2) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_h.{idx+1,3} <- if tmp_f.{idx+1,3} = 0. then 0. else +. h_two a n i m s1 s3) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_h.{idx+1,4} <- if tmp_f.{idx+1,4} = 0. then 0. else +. h_two a n j m s2 s3) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_h.{idx+1,5} <- if tmp_f.{idx+1,5} = 0. then 0. else -. h_two a l i m s1 s3) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_h.{idx+1,6} <- if tmp_f.{idx+1,6} = 0. then 0. else -. h_two a l j m s2 s3) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_h.{idx+1,7} <- if tmp_f.{idx+1,7} = 0. then 0. else +. h_two a k m i s3 s1) mos_cabs;
|
|
Array.iteri (fun idx a -> tmp_h.{idx+1,8} <- if tmp_f.{idx+1,8} = 0. then 0. else -. h_two a k j i s2 s1) mos_cabs;
|
|
dot ~n:(8*ncabs) (Mat.as_vec tmp_h) (Mat.as_vec tmp_f)
|
|
| Beta, Alfa, Beta
|
|
| Alfa, Beta, Alfa -> assert false
|
|
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
|
|
result
|
|
|
|
|
|
|