mirror of
https://gitlab.com/scemama/QCaml.git
synced 2024-09-01 17:23:41 +02:00
2310 lines
66 KiB
OCaml
2310 lines
66 KiB
OCaml
|
|
let png_image = print_endline ;;
|
|
|
|
|
|
open Lacaml.D
|
|
|
|
|
|
let basis_filename = "/home/scemama/qp2/data/basis/6-31g"
|
|
let aux_basis_filename = "/home/scemama/qp2/data/basis/cc-pvdz"
|
|
let nuclei = Nuclei.of_zmt_string "li"
|
|
let frozen_core = false
|
|
let multiplicity = 2
|
|
let state = 1
|
|
|
|
let f12 = F12factor.gaussian_geminal 1.0
|
|
let basis = Basis.of_nuclei_and_basis_filenames ~f12 ~nuclei [basis_filename]
|
|
let charge = 0
|
|
|
|
|
|
let simulation =
|
|
Simulation.make
|
|
~charge ~multiplicity ~nuclei
|
|
~cartesian:true
|
|
basis
|
|
|
|
|
|
let n_elec_alfa, n_elec_beta, n_elec =
|
|
let e = Simulation.electrons simulation in
|
|
Electrons.(n_alfa e, n_beta e, n_elec e)
|
|
|
|
let hf = HartreeFock.make ~guess:`Hcore ~max_scf:1 simulation ;;
|
|
|
|
let mo_basis = MOBasis.of_hartree_fock hf
|
|
|
|
let f12 = Util.of_some @@ Simulation.f12 simulation
|
|
|
|
let mo_num = MOBasis.size mo_basis
|
|
|
|
let pp_spindet = Spindeterminant.pp mo_num
|
|
|
|
let pp_det = Determinant.pp mo_num
|
|
|
|
;;
|
|
|
|
|
|
|
|
let simulation_aux =
|
|
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 ~f12 nuclei
|
|
|> Simulation.make ~charge ~multiplicity ~nuclei
|
|
|
|
|
|
let aux_basis =
|
|
MOBasis.of_mo_basis simulation_aux mo_basis
|
|
|
|
let aux_num =
|
|
MOBasis.size aux_basis
|
|
|
|
|
|
let () = ignore @@ MOBasis.f12_ints aux_basis
|
|
let () = ignore @@ MOBasis.two_e_ints aux_basis
|
|
|
|
let cancel_singles = false
|
|
|
|
let mos_cabs =
|
|
Util.list_range (mo_num+1) aux_num
|
|
|
|
let mos_abs =
|
|
Util.list_range 1 aux_num
|
|
|
|
let mos_in =
|
|
Util.list_range 1 mo_num
|
|
|
|
let mos_a k =
|
|
Determinant.alfa k
|
|
|> Spindeterminant.to_list
|
|
|
|
let mos_b k =
|
|
Determinant.beta k
|
|
|> Spindeterminant.to_list
|
|
|
|
|
|
let h_one =
|
|
let h =
|
|
MOBasis.one_e_ints aux_basis
|
|
in fun i j _ -> h.{i,j}
|
|
|
|
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
|
|
|
|
|
|
|
|
|
|
let f_two =
|
|
let two_e_ints = MOBasis.f12_ints aux_basis in
|
|
let f2 i j k l (s:Spin.t) (s':Spin.t) =
|
|
if s' <> s then
|
|
0.5 *. F12.get_phys two_e_ints i j k l
|
|
else
|
|
0.5 *. (
|
|
(F12.get_phys two_e_ints i j k l) -.
|
|
(F12.get_phys two_e_ints i j l k) )
|
|
in
|
|
let f3 i j k l (s:Spin.t) (s':Spin.t) =
|
|
if (i=k && j<>l) || (j=l && i<>k) then
|
|
0.
|
|
else
|
|
f2 i j k l s s'
|
|
in
|
|
if cancel_singles then f3 else f2
|
|
|
|
let f_one = fun _ _ _ -> 0.
|
|
|
|
(*
|
|
let f_two = h_two
|
|
|
|
let f_one = h_one
|
|
*)
|
|
|
|
let f12_integrals mo_basis =
|
|
( f_one, f_two, None )
|
|
|
|
let h_ij mo_basis ki kj =
|
|
let integrals =
|
|
List.map (fun f -> f mo_basis)
|
|
[ CI.h_integrals ]
|
|
in
|
|
CIMatrixElement.make integrals ki kj
|
|
|> List.hd
|
|
|
|
|
|
let f_ij mo_basis ki kj =
|
|
let integrals =
|
|
List.map (fun f -> f mo_basis)
|
|
[ f12_integrals ]
|
|
in
|
|
CIMatrixElement.make integrals ki kj
|
|
|> List.hd
|
|
|
|
|
|
let hf_ij mo_basis ki kj =
|
|
let integrals =
|
|
List.map (fun f -> f mo_basis)
|
|
[ CI.h_integrals ; f12_integrals ]
|
|
in
|
|
CIMatrixElement.make integrals ki kj
|
|
|
|
|
|
|
|
|
|
let is_a_double det_space =
|
|
let mo_class = DeterminantSpace.mo_class det_space in
|
|
let mo_num = Array.length @@ MOClass.mo_class_array mo_class in
|
|
let m l =
|
|
List.fold_left (fun accu i ->
|
|
let j = i-1 in Bitstring.logor accu (Bitstring.shift_left_one mo_num j)
|
|
) (Bitstring.zero mo_num) l
|
|
in
|
|
let aux_mask = m (MOClass.auxiliary_mos mo_class) in
|
|
fun k ->
|
|
let alfa =
|
|
Determinant.alfa k
|
|
|> Spindeterminant.bitstring
|
|
in
|
|
let beta =
|
|
Determinant.beta k
|
|
|> Spindeterminant.bitstring
|
|
in
|
|
let a = Bitstring.logand aux_mask alfa
|
|
and b = Bitstring.logand aux_mask beta
|
|
in
|
|
match Bitstring.popcount a + Bitstring.popcount b with
|
|
| 2 | 1 -> true
|
|
| 0 | _ -> false
|
|
|
|
|
|
|
|
let in_space =
|
|
DeterminantSpace.fci_f12_of_mo_basis aux_basis ~frozen_core mo_num
|
|
|
|
let aux_space =
|
|
DeterminantSpace.fci_of_mo_basis aux_basis ~frozen_core
|
|
|
|
let det_space_in () =
|
|
DeterminantSpace.determinant_stream in_space
|
|
|
|
let det_space_out () =
|
|
let s =
|
|
DeterminantSpace.determinant_stream aux_space
|
|
in
|
|
Stream.from (fun _ ->
|
|
try
|
|
let is_a_double = is_a_double in_space in
|
|
let rec result () =
|
|
let ki = Stream.next s in
|
|
if is_a_double ki then
|
|
Some (ki,ki)
|
|
else
|
|
result ()
|
|
in
|
|
result ()
|
|
with Stream.Failure -> None
|
|
)
|
|
|
|
|
|
let ci = CI.make ~n_states:state in_space
|
|
|
|
let ci_coef, ci_energy = Lazy.force ci.eigensystem
|
|
|
|
let _ = print_newline ()
|
|
|
|
let p12 det_space =
|
|
let mo_class = DeterminantSpace.mo_class det_space in
|
|
let mo_num = Array.length @@ MOClass.mo_class_array mo_class in
|
|
let m l =
|
|
List.fold_left (fun accu i ->
|
|
let j = i-1 in Bitstring.logor accu (Bitstring.shift_left_one mo_num j)
|
|
) (Bitstring.zero mo_num) l
|
|
in
|
|
let aux_mask = m (MOClass.auxiliary_mos mo_class) in
|
|
let not_aux_mask =
|
|
Bitstring.(shift_left_one mo_num (mo_num-1) |> minus_one |> logxor aux_mask)
|
|
in
|
|
fun k ->
|
|
let alfa =
|
|
Determinant.alfa k
|
|
|> Spindeterminant.bitstring
|
|
in
|
|
let beta =
|
|
Determinant.beta k
|
|
|> Spindeterminant.bitstring
|
|
in
|
|
let a = Bitstring.logand aux_mask alfa
|
|
and b = Bitstring.logand aux_mask beta
|
|
in
|
|
match Bitstring.popcount a, Bitstring.popcount b with
|
|
| 2, 0
|
|
| 0, 2 -> Some (Determinant.negate_phase k)
|
|
| 1, 1 -> Some (Determinant.of_spindeterminants
|
|
(Spindeterminant.of_bitstring @@
|
|
Bitstring.(logor b (logand not_aux_mask alfa)) )
|
|
(Spindeterminant.of_bitstring @@
|
|
Bitstring.(logor a (logand not_aux_mask beta))
|
|
) )
|
|
| 1, 0
|
|
| 0, 1 -> Some k
|
|
| _ -> None
|
|
|
|
|
|
|
|
let out_list =
|
|
Util.stream_to_list (det_space_out ())
|
|
|
|
let in_list =
|
|
Util.stream_to_list (det_space_in ())
|
|
|
|
let det_a = Array.of_list out_list
|
|
|> Array.map (fun (i,_) -> i)
|
|
|
|
let det_I =
|
|
let n = 123456789 in
|
|
in_list
|
|
|> List.map (fun k -> (Random.int n, k))
|
|
|> List.sort compare
|
|
|> List.map (fun (_,k) -> k)
|
|
|> Array.of_list
|
|
|
|
|
|
let generate_alphas ki kj exc cabs l r =
|
|
(* Check input excitation degree *)
|
|
let d = Determinant.degree ki kj in
|
|
if d <> exc then
|
|
Printf.sprintf "Invalid excitation degree. Expected %d and found %d." exc d
|
|
|> failwith;
|
|
|
|
(* Generate single excitations *)
|
|
let all_singles ki =
|
|
let mos_a, mos_b = Determinant.to_lists ki in
|
|
[ List.map (fun hole ->
|
|
List.map (fun particle ->
|
|
if hole = particle then None else
|
|
Some (Determinant.single_excitation Spin.Alfa hole particle ki)
|
|
) mos_abs
|
|
) mos_a
|
|
;
|
|
List.map (fun hole ->
|
|
List.map (fun particle ->
|
|
if hole = particle then None else
|
|
Some (Determinant.single_excitation Spin.Beta hole particle ki)
|
|
) mos_abs
|
|
) mos_b
|
|
]
|
|
|> List.concat
|
|
|> List.concat
|
|
|> Util.list_some
|
|
|> List.filter (fun x -> Determinant.is_none x = false)
|
|
|> List.map (fun x -> Determinant.set_phase Phase.Pos x)
|
|
in
|
|
|
|
(* Generate double excitations *)
|
|
let all_doubles ki =
|
|
let mos_a, mos_b = Determinant.to_lists ki in
|
|
[ List.map (fun hole -> (* Alpha-Alpha *)
|
|
List.map (fun particle ->
|
|
List.map (fun hole' ->
|
|
List.map (fun particle' ->
|
|
if hole = particle || hole' = particle' then None else
|
|
if hole' > hole && particle' < particle then
|
|
Some (Determinant.double_excitation Spin.Alfa hole particle Spin.Alfa hole' particle' ki)
|
|
else None
|
|
) mos_abs
|
|
) mos_a
|
|
) mos_abs
|
|
) mos_a
|
|
;
|
|
List.map (fun hole -> (* Beta-Beta *)
|
|
List.map (fun particle ->
|
|
List.map (fun hole' ->
|
|
List.map (fun particle' ->
|
|
if hole' > hole && particle' < particle && hole <> particle && hole' <> particle' then
|
|
Some (Determinant.double_excitation Spin.Beta hole particle Spin.Beta hole' particle' ki)
|
|
else None
|
|
) mos_abs
|
|
) mos_b
|
|
) mos_abs
|
|
) mos_b
|
|
;
|
|
List.map (fun hole -> (* Alpha-Beta *)
|
|
List.map (fun particle ->
|
|
List.map (fun hole' ->
|
|
List.map (fun particle' ->
|
|
if hole = particle || hole' = particle' then None else
|
|
Some (Determinant.double_excitation Spin.Alfa hole particle Spin.Beta hole' particle' ki)
|
|
) mos_abs
|
|
) mos_b
|
|
) mos_abs
|
|
) mos_a
|
|
]
|
|
|> List.concat
|
|
|> List.concat
|
|
|> List.concat
|
|
|> List.concat
|
|
|> Util.list_some
|
|
|> List.filter (fun x -> Determinant.is_none x = false)
|
|
|> List.map (fun x -> Determinant.set_phase Phase.Pos x)
|
|
in
|
|
|
|
(* Generate left and right excitations *)
|
|
let al =
|
|
match l with
|
|
| 1 -> all_singles ki
|
|
| 2 -> all_doubles ki
|
|
| _ -> assert false
|
|
in
|
|
let ar =
|
|
match r with
|
|
| 1 -> all_singles kj
|
|
| 2 -> all_doubles kj
|
|
| _ -> assert false
|
|
in
|
|
|
|
let mo_class = DeterminantSpace.mo_class in_space in
|
|
let m l =
|
|
List.fold_left (fun accu i ->
|
|
let j = i-1 in Bitstring.logor accu (Bitstring.shift_left_one mo_num j)
|
|
) (Bitstring.zero mo_num) l
|
|
in
|
|
let aux_mask = m (MOClass.auxiliary_mos mo_class) in
|
|
let good_cabs k =
|
|
let alfa =
|
|
Determinant.alfa k
|
|
|> Spindeterminant.bitstring
|
|
in
|
|
let beta =
|
|
Determinant.beta k
|
|
|> Spindeterminant.bitstring
|
|
in
|
|
let a = Bitstring.logand aux_mask alfa
|
|
and b = Bitstring.logand aux_mask beta
|
|
in
|
|
Bitstring.popcount a + Bitstring.popcount b = cabs
|
|
in
|
|
let good_lr k =
|
|
Determinant.degree ki k = l &&
|
|
Determinant.degree k kj = r
|
|
in
|
|
|
|
|
|
|
|
(* Merge lists in a set of unique determinants *)
|
|
List.concat [ al; ar ]
|
|
|> List.sort_uniq compare
|
|
|
|
(* Filter out all determinants with incorrect numbers of electrons in the CABS *)
|
|
|> List.filter good_cabs
|
|
|
|
(* Filter out all determinants with incorrect excitation numbers *)
|
|
|> List.filter good_lr
|
|
|
|
|
|
let compute_HaaF ki alphas kj =
|
|
List.fold_left (fun accu alpha -> accu
|
|
+. h_ij aux_basis ki alpha
|
|
*. f_ij aux_basis alpha kj
|
|
) 0. alphas
|
|
|
|
|
|
let check n integral_value exc cabs lexc rexc =
|
|
let cpudet, cpuint = ref 0., ref 0. in
|
|
let det_list =
|
|
match n with
|
|
| 0 -> det_I
|
|
| n -> Array.sub det_I 0 n
|
|
in
|
|
let result =
|
|
if Parallel.master then Printf.printf "Checking ... \n%!";
|
|
let percent = ref 0 in
|
|
let task (i,ki) =
|
|
(let p = (10 * (i+1))/(Array.length det_list) in
|
|
if p <> !percent then
|
|
( percent := p ; if Parallel.master then Printf.printf " - %3d %%\n%!" (p*10) ));
|
|
Array.mapi (fun j kj ->
|
|
if i > j || Determinant.degree ki kj <> exc then
|
|
(0,0,0.,0.)
|
|
else
|
|
begin
|
|
let alphas = generate_alphas ki kj exc cabs lexc rexc in
|
|
let det_value =
|
|
let t0 = Unix.gettimeofday () in
|
|
let result = compute_HaaF ki alphas kj in
|
|
cpudet := !cpudet +. Unix.gettimeofday () -. t0;
|
|
result
|
|
in
|
|
let int_value =
|
|
let t0 = Unix.gettimeofday () in
|
|
let result = integral_value ki kj in
|
|
cpuint := !cpuint +. Unix.gettimeofday () -. t0;
|
|
result
|
|
in
|
|
(* Printf.printf "%d %d %e %e\n%!" i j det_value int_value; *)
|
|
(i,j,det_value,int_value)
|
|
end
|
|
) det_list
|
|
in
|
|
det_list
|
|
|> Array.mapi (fun i ki -> (i,ki))
|
|
|> Array.to_list
|
|
|> Stream.of_list
|
|
|> Farm.run ~f:task
|
|
|> Util.stream_to_list
|
|
|> Array.concat
|
|
in
|
|
let i,j,d,v =
|
|
let rec aux k imax jmax emax dmax vmax =
|
|
if k = -1 then
|
|
imax, jmax, dmax, vmax
|
|
else
|
|
let i, j, d, v = result.(k) in
|
|
let e = abs_float (d -. v) in
|
|
if e >= emax then
|
|
aux (k-1) i j e d v
|
|
else
|
|
aux (k-1) imax jmax emax dmax vmax
|
|
in
|
|
aux (Array.length result - 1) 0 0 0. 0. 0.
|
|
in
|
|
let error = abs_float (d -. v) in
|
|
if error < epsilon_float then
|
|
(*
|
|
Printf.printf "OK: %e\n%!" error
|
|
*)
|
|
if Parallel.master then Printf.printf "OK: (%d, %d) | %e %e | %e | cpu : %f %f\n%!" i j d v error !cpudet !cpuint
|
|
else
|
|
if Parallel.master then Printf.printf "Failed: (%d, %d) | %e %e | %e | cpu : %f %f\n%!" i j d v error !cpudet !cpuint
|
|
|
|
|
|
let sum l f = List.fold_left (fun accu i -> accu +. f i) 0. l
|
|
|
|
let array_3_init d1 d2 d3 f =
|
|
let result =
|
|
Bigarray.(Array3.create Float64 fortran_layout) d1 d2 d3
|
|
in
|
|
for k=1 to d3 do
|
|
for j=1 to d2 do
|
|
for i=1 to d1 do
|
|
result.{i,j,k} <- f i j k
|
|
done
|
|
done
|
|
done;
|
|
result
|
|
|
|
|
|
let array_4_init d1 d2 d3 d4 f =
|
|
let result =
|
|
Bigarray.(Genarray.create Float64 fortran_layout) [| d1;d2;d3;d4 |]
|
|
in
|
|
for l=1 to d4 do
|
|
for k=1 to d3 do
|
|
for j=1 to d2 do
|
|
for i=1 to d1 do
|
|
result.{i,j,k,l} <- f i j k l
|
|
done
|
|
done
|
|
done
|
|
done;
|
|
result
|
|
|
|
|
|
let array_5_init d1 d2 d3 d4 d5 f =
|
|
let result =
|
|
Bigarray.(Genarray.create Float64 fortran_layout) [| d1;d2;d3;d4;d5 |]
|
|
in
|
|
for m=1 to d5 do
|
|
for l=1 to d4 do
|
|
for k=1 to d3 do
|
|
for j=1 to d2 do
|
|
for i=1 to d1 do
|
|
result.{i,j,k,l,m} <- f i j k l m
|
|
done
|
|
done
|
|
done
|
|
done
|
|
done;
|
|
result
|
|
|
|
(* ----- *)
|
|
|
|
let _ = png_image "0_1_11.png"
|
|
|
|
let integral_value ki kj =
|
|
(* Alpha-Beta *)
|
|
let s, s' = Spin.(Alfa, Beta) in
|
|
let mos_a, mos_b = mos_a ki, mos_b ki in
|
|
|
|
sum mos_cabs (fun a ->
|
|
sum mos_a (fun i ->
|
|
(h_one i a s +.
|
|
sum mos_a (fun j -> h_two i j a j s s ) +.
|
|
sum mos_b (fun j -> h_two i j a j s s') ) *.
|
|
(f_one a i s +.
|
|
sum mos_a (fun j -> f_two a j i j s s ) +.
|
|
sum mos_b (fun j -> f_two a j i j s s') )
|
|
)
|
|
+.
|
|
sum mos_b (fun i ->
|
|
(h_one i a s +.
|
|
sum mos_b (fun j -> h_two i j a j s' s') +.
|
|
sum mos_a (fun j -> h_two i j a j s' s ) ) *.
|
|
(f_one a i s +.
|
|
sum mos_a (fun j -> f_two a j i j s' s ) +.
|
|
sum mos_b (fun j -> f_two a j i j s' s') )
|
|
)
|
|
)
|
|
|
|
let _ =
|
|
check 0 integral_value 0 1 1 1
|
|
|
|
|
|
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 ))
|
|
|
|
|
|
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 ))
|
|
|
|
|
|
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 ))
|
|
|
|
|
|
let m_0111_2Ha_1F =
|
|
Mat.init_cols mo_num mo_num (fun i j -> sum mos_cabs (fun a ->
|
|
h_two i j a j Spin.Alfa Spin.Alfa *. f_one a j Spin.Alfa ))
|
|
|
|
|
|
let m_0111_2Hb_1F =
|
|
Mat.init_cols mo_num mo_num (fun i j -> sum mos_cabs (fun a ->
|
|
h_two i j a j Spin.Alfa Spin.Beta *. f_one a j Spin.Alfa ))
|
|
|
|
|
|
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
|
|
)
|
|
)
|
|
|
|
|
|
let m_0111_2Hb_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.Beta *.
|
|
f_two a k i k Spin.Alfa Spin.Alfa
|
|
)
|
|
)
|
|
|
|
|
|
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
|
|
)
|
|
)
|
|
|
|
|
|
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 i j a j Spin.Alfa Spin.Beta *.
|
|
f_two a k i k Spin.Alfa Spin.Beta
|
|
)
|
|
)
|
|
|
|
|
|
|
|
let integral_value ki kj =
|
|
|
|
let mos_a, mos_b = mos_a ki, mos_b ki in
|
|
|
|
sum mos_a (fun i ->
|
|
m_0111_1H_1F.{i} +.
|
|
sum mos_a (fun j ->
|
|
m_0111_1H_2Fa.{i,j} +. m_0111_2Ha_1F.{i,j} +.
|
|
sum mos_a (fun k -> m_0111_2Ha_2Fa.{i,j,k}) +.
|
|
sum mos_b (fun k -> m_0111_2Ha_2Fb.{i,j,k})) +.
|
|
sum mos_b (fun j ->
|
|
m_0111_1H_2Fb.{i,j} +. m_0111_2Hb_1F.{i,j} +.
|
|
sum mos_a (fun k -> m_0111_2Hb_2Fa.{i,j,k}) +.
|
|
sum mos_b (fun k -> m_0111_2Hb_2Fb.{i,j,k}))
|
|
)
|
|
+.
|
|
sum mos_b (fun i ->
|
|
m_0111_1H_1F.{i} +.
|
|
sum mos_b (fun j ->
|
|
m_0111_1H_2Fa.{i,j} +. m_0111_2Ha_1F.{i,j} +.
|
|
sum mos_b (fun k -> m_0111_2Ha_2Fa.{i,j,k}) +.
|
|
sum mos_a (fun k -> m_0111_2Ha_2Fb.{i,j,k})) +.
|
|
sum mos_a (fun j ->
|
|
m_0111_1H_2Fb.{i,j} +. m_0111_2Hb_1F.{i,j} +.
|
|
sum mos_b (fun k -> m_0111_2Hb_2Fa.{i,j,k}) +.
|
|
sum mos_a (fun k -> m_0111_2Hb_2Fb.{i,j,k}))
|
|
)
|
|
|
|
let _ =
|
|
check 0 integral_value 0 1 1 1
|
|
|
|
|
|
let _ = png_image "0_1_22.png"
|
|
|
|
let integral_value ki kj =
|
|
(* mos unoccupied in both I and J *)
|
|
let mos_virt_a, mos_virt_b =
|
|
Array.init mo_num (fun i -> Some (i+1)) ,
|
|
Array.init mo_num (fun i -> Some (i+1))
|
|
in
|
|
List.iter (fun i -> mos_virt_a.(i-1) <- None) (mos_a ki);
|
|
List.iter (fun i -> mos_virt_a.(i-1) <- None) (mos_a kj);
|
|
List.iter (fun i -> mos_virt_b.(i-1) <- None) (mos_b ki);
|
|
List.iter (fun i -> mos_virt_b.(i-1) <- None) (mos_b kj);
|
|
|
|
let mos_virt_a, mos_virt_b =
|
|
Array.to_list mos_virt_a |> Util.list_some,
|
|
Array.to_list mos_virt_b |> Util.list_some
|
|
in
|
|
|
|
(* Alpha-Beta *)
|
|
let s, s' = Spin.(Alfa, Beta) in
|
|
sum mos_cabs (fun a ->
|
|
sum (mos_b ki) (fun j ->
|
|
sum (mos_a ki) (fun i ->
|
|
sum mos_virt_b (fun k ->
|
|
h_two i j a k s s' *. f_two a k i j s s'
|
|
)))
|
|
+.
|
|
sum (mos_a ki) (fun j ->
|
|
sum (mos_b ki) (fun i ->
|
|
sum mos_virt_a (fun k ->
|
|
h_two i j a k s s' *. f_two a k i j s s'
|
|
))))
|
|
+.
|
|
(* Alpha-Alpha / Beta-Beta *)
|
|
List.fold_left (fun accu (mos_virt,mos,s) ->
|
|
let s' = s in accu +.
|
|
sum mos (fun j ->
|
|
sum mos (fun i -> if i < j then 0. else
|
|
sum mos_virt (fun k ->
|
|
sum mos_cabs (fun a ->
|
|
h_two i j a k s s' *. f_two a k i j s s'
|
|
))))
|
|
) 0. [ (mos_virt_a,mos_a ki,Spin.Alfa) ; (mos_virt_b,mos_b ki,Spin.Beta)]
|
|
|
|
let _ =
|
|
check 0 integral_value 0 1 2 2
|
|
|
|
|
|
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
|
|
) )
|
|
|
|
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
|
|
) )
|
|
|
|
|
|
let integral_value ki kj =
|
|
(* mos unoccupied in both I and J *)
|
|
let mos_virt_a, mos_virt_b =
|
|
Array.init mo_num (fun i -> Some (i+1)) ,
|
|
Array.init mo_num (fun i -> Some (i+1))
|
|
in
|
|
List.iter (fun i -> mos_virt_a.(i-1) <- None) (mos_a ki);
|
|
List.iter (fun i -> mos_virt_a.(i-1) <- None) (mos_a kj);
|
|
List.iter (fun i -> mos_virt_b.(i-1) <- None) (mos_b ki);
|
|
List.iter (fun i -> mos_virt_b.(i-1) <- None) (mos_b kj);
|
|
|
|
let mos_virt_a, mos_virt_b =
|
|
Array.to_list mos_virt_a |> Util.list_some,
|
|
Array.to_list mos_virt_b |> Util.list_some
|
|
in
|
|
|
|
let result =
|
|
sum (mos_a ki) (fun i ->
|
|
sum (mos_b ki) (fun j ->
|
|
sum mos_virt_b (fun k -> m_0122_Hab.{i,j,k} )))
|
|
+.
|
|
sum (mos_b ki) (fun i ->
|
|
sum (mos_a ki) (fun j ->
|
|
sum mos_virt_a (fun k -> m_0122_Hab.{i,j,k} )))
|
|
+.
|
|
(* Alpha-Alpha / Beta-Beta *)
|
|
List.fold_left (fun accu (mos_virt,mos) ->
|
|
accu +.
|
|
sum mos (fun j ->
|
|
sum mos (fun i -> if i <= j then 0. else
|
|
sum mos_virt (fun k -> m_0122_Haa.{i,j,k} )))
|
|
) 0. [ (mos_virt_a,mos_a ki) ; (mos_virt_b,mos_b ki)]
|
|
in
|
|
result
|
|
|
|
let _ =
|
|
check 0 integral_value 0 1 2 2
|
|
|
|
|
|
let _ = png_image "0_2_22.png"
|
|
|
|
let integral_value ki kj =
|
|
(* Alpha-Beta *)
|
|
let s, s' = Spin.(Alfa, Beta) in
|
|
sum mos_cabs (fun a ->
|
|
sum mos_cabs (fun b ->
|
|
sum (mos_a ki) (fun i ->
|
|
sum (mos_b ki) (fun j ->
|
|
h_two i j a b s s' *. f_two a b i j s s'
|
|
))))
|
|
+.
|
|
(* Alpha-Alpha / Beta-Beta *)
|
|
List.fold_left (fun accu (mos,s) ->
|
|
let s' = s in accu +.
|
|
sum mos_cabs (fun b ->
|
|
sum mos_cabs (fun a -> if b > a then 0. else
|
|
sum mos (fun j ->
|
|
sum mos (fun i -> if i < j then 0. else
|
|
h_two i j a b s s' *. f_two a b i j s s'
|
|
))))
|
|
) 0. [ (mos_a ki,Spin.Alfa) ; (mos_b ki,Spin.Beta)]
|
|
|
|
let _ =
|
|
check 0 integral_value 0 2 2 2
|
|
|
|
|
|
let m_0222_Haa =
|
|
Mat.init_cols mo_num mo_num (fun i j ->
|
|
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
|
|
) )
|
|
)
|
|
|
|
let m_0222_Hab =
|
|
Mat.init_cols mo_num mo_num (fun i j ->
|
|
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
|
|
) )
|
|
)
|
|
|
|
let integral_value ki kj =
|
|
(* Alpha-Beta *)
|
|
sum (mos_a ki) (fun i ->
|
|
sum (mos_b ki) (fun j ->
|
|
m_0222_Hab.{i,j}
|
|
))
|
|
+.
|
|
(* Alpha-Alpha / Beta-Beta *)
|
|
List.fold_left (fun accu (mos) ->
|
|
accu +.
|
|
sum mos (fun j ->
|
|
sum mos (fun i -> if i < j then 0. else
|
|
m_0222_Haa.{i,j}
|
|
))
|
|
) 0. [ (mos_a ki) ; (mos_b ki)]
|
|
|
|
let _ =
|
|
check 0 integral_value 0 2 2 2
|
|
|
|
|
|
let _ = png_image "1_1_11.png"
|
|
|
|
let integral_value 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_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 =
|
|
let s' = Spin.other s in
|
|
|
|
sum mos_cabs (fun a ->
|
|
(h_one i a s +.
|
|
sum mos_i (fun j -> h_two i j a j s s ) +.
|
|
sum mos_i' (fun j -> h_two i j a j s s') ) *.
|
|
(f_one a k s +.
|
|
sum mos_j (fun j -> f_two a j k j s s ) +.
|
|
sum mos_j' (fun j -> f_two a j k j s s') )
|
|
)
|
|
in
|
|
match phase with
|
|
| Phase.Pos -> result
|
|
| Phase.Neg -> -. result
|
|
|
|
|
|
let _ = check 0 integral_value 1 1 1 1
|
|
|
|
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 ))
|
|
|
|
|
|
let m_1111_1H_2Fa =
|
|
array_3_init mo_num mo_num mo_num (fun i j k ->
|
|
sum mos_cabs (fun a ->
|
|
h_one i a Spin.Alfa *. f_two a j k j Spin.Alfa Spin.Alfa ))
|
|
|
|
|
|
let m_1111_1H_2Fb =
|
|
array_3_init mo_num mo_num mo_num (fun i j k ->
|
|
sum mos_cabs (fun a ->
|
|
h_one i a Spin.Alfa *. f_two a j k j Spin.Alfa Spin.Beta ))
|
|
|
|
|
|
let m_1111_2Ha_1F =
|
|
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_one a k Spin.Alfa ))
|
|
|
|
|
|
let m_1111_2Hb_1F =
|
|
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.Beta *. f_one a k Spin.Alfa ))
|
|
|
|
|
|
let m_1111_2Ha_2Fa =
|
|
array_4_init mo_num mo_num mo_num mo_num (fun i j k l ->
|
|
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
|
|
)
|
|
)
|
|
|
|
|
|
let m_1111_2Hb_2Fa =
|
|
array_4_init mo_num mo_num mo_num mo_num (fun i j k l ->
|
|
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
|
|
)
|
|
)
|
|
|
|
|
|
let m_1111_2Ha_2Fb =
|
|
array_4_init mo_num mo_num mo_num mo_num (fun i j k l ->
|
|
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
|
|
)
|
|
)
|
|
|
|
|
|
let m_1111_2Hb_2Fb =
|
|
array_4_init mo_num mo_num mo_num mo_num (fun i j k l ->
|
|
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
|
|
)
|
|
)
|
|
|
|
|
|
let integral_value 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_a, mos_b =
|
|
match s with
|
|
| Spin.Alfa -> mos_a ki, mos_b ki
|
|
| Spin.Beta -> mos_b ki, mos_a ki
|
|
in
|
|
|
|
let result =
|
|
m_1111_1H_1F.{i,k} +.
|
|
sum mos_a (fun j -> if j = i then 0. else
|
|
m_1111_1H_2Fa.{i,j,k} +. m_1111_2Ha_1F.{i,j,k} +.
|
|
sum mos_a (fun l -> if l = i then 0. else m_1111_2Ha_2Fa.{i,j,k,l}) +.
|
|
sum mos_b (fun l -> m_1111_2Ha_2Fb.{i,j,k,l})) +.
|
|
sum mos_b (fun j ->
|
|
m_1111_1H_2Fb.{i,j,k} +. m_1111_2Hb_1F.{i,j,k} +.
|
|
sum mos_a (fun l -> if l = i then 0. else m_1111_2Hb_2Fa.{i,j,k,l}) +.
|
|
sum mos_b (fun l -> m_1111_2Hb_2Fb.{i,j,k,l}))
|
|
in
|
|
|
|
match phase with
|
|
| Phase.Pos -> result
|
|
| Phase.Neg -> -. result
|
|
|
|
let _ =
|
|
check 0 integral_value 1 1 1 1
|
|
|
|
|
|
let _ = png_image "1_1_12.png"
|
|
|
|
let integral_value 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_j, mos_j' =
|
|
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 result =
|
|
let s' = Spin.other s in
|
|
sum mos_cabs (fun a ->
|
|
sum mos_j (fun j ->
|
|
(h_one j a s +.
|
|
sum mos_i (fun m -> h_two j m a m s s ) +.
|
|
sum mos_i' (fun m -> h_two j m a m s s') ) *.
|
|
(f_two a i j k s s )
|
|
)
|
|
+.
|
|
sum mos_j' (fun j ->
|
|
(h_one j a s +.
|
|
sum mos_i (fun m -> h_two j m a m s' s ) +.
|
|
sum mos_i' (fun m -> h_two j m a m s' s') ) *.
|
|
(f_two a i j k s' s )
|
|
)
|
|
)
|
|
in
|
|
match phase with
|
|
| Phase.Pos -> result
|
|
| Phase.Neg -> -. result
|
|
|
|
|
|
let _ = check 0 integral_value 1 1 1 2
|
|
|
|
let m_1112_1H_2Fa =
|
|
array_3_init mo_num mo_num mo_num (fun i j k ->
|
|
sum mos_cabs (fun a ->
|
|
h_one j a Spin.Alfa *. f_two a i j k Spin.Alfa Spin.Alfa ))
|
|
|
|
|
|
let m_1112_1H_2Fb =
|
|
array_3_init mo_num mo_num mo_num (fun i j k ->
|
|
sum mos_cabs (fun a ->
|
|
h_one j a Spin.Beta *. f_two a i j k Spin.Alfa Spin.Beta ))
|
|
|
|
|
|
let m_1112_2Ha_2Fa =
|
|
array_4_init mo_num mo_num mo_num mo_num (fun i j k m ->
|
|
sum mos_cabs (fun a ->
|
|
h_two j m a m Spin.Alfa Spin.Alfa *.
|
|
f_two a i j k Spin.Alfa Spin.Alfa
|
|
)
|
|
)
|
|
|
|
|
|
let m_1112_2Hb_2Fa =
|
|
array_4_init mo_num mo_num mo_num mo_num (fun i j k m ->
|
|
sum mos_cabs (fun a ->
|
|
h_two j m a m Spin.Alfa Spin.Beta *.
|
|
f_two a i j k Spin.Alfa Spin.Alfa
|
|
)
|
|
)
|
|
|
|
let m_1112_2Ha_2Fb =
|
|
array_4_init mo_num mo_num mo_num mo_num (fun i j k m ->
|
|
sum mos_cabs (fun a ->
|
|
h_two j m a m Spin.Alfa Spin.Alfa *.
|
|
f_two a i j k Spin.Alfa Spin.Beta
|
|
)
|
|
)
|
|
|
|
|
|
let m_1112_2Hb_2Fb =
|
|
array_4_init mo_num mo_num mo_num mo_num (fun i j k m ->
|
|
sum mos_cabs (fun a ->
|
|
h_two j m a m Spin.Alfa Spin.Beta *.
|
|
f_two a i j k Spin.Alfa Spin.Beta
|
|
)
|
|
)
|
|
|
|
|
|
let integral_value 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_j, mos_j' =
|
|
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 result =
|
|
sum mos_j (fun j ->
|
|
m_1112_1H_2Fa.{i,j,k} +.
|
|
sum mos_i (fun m -> m_1112_2Ha_2Fa.{i,j,k,m}) +.
|
|
sum mos_i' (fun m -> m_1112_2Hb_2Fa.{i,j,k,m})) +.
|
|
sum mos_j' (fun j ->
|
|
m_1112_1H_2Fb.{i,j,k} +.
|
|
sum mos_i (fun m -> m_1112_2Hb_2Fb.{i,j,k,m}) +.
|
|
sum mos_i' (fun m -> m_1112_2Ha_2Fb.{i,j,k,m}))
|
|
in
|
|
|
|
match phase with
|
|
| Phase.Pos -> result
|
|
| Phase.Neg -> -. result
|
|
|
|
|
|
|
|
let _ = check 0 integral_value 1 1 1 2
|
|
|
|
let _ = png_image "1_1_21.png"
|
|
|
|
let integral_value 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_j, mos_j' =
|
|
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 kj, mos_b kj
|
|
| Spin.Beta -> mos_b kj, mos_a kj
|
|
in
|
|
|
|
let result =
|
|
let s' = Spin.other s in
|
|
sum mos_cabs (fun a ->
|
|
sum mos_j (fun j ->
|
|
(h_two i j k a s s ) *.
|
|
(f_one a j s +.
|
|
sum mos_i (fun m -> f_two m a m j s s) +.
|
|
sum mos_i' (fun m -> f_two m a m j s' s) )
|
|
)
|
|
+.
|
|
sum mos_j' (fun j ->
|
|
(h_two i j k a s s') *.
|
|
(f_one a j s' +.
|
|
sum mos_i (fun m -> f_two m a m j s s') +.
|
|
sum mos_i' (fun m -> f_two m a m j s' s') )
|
|
))
|
|
in
|
|
match phase with
|
|
| Phase.Pos -> result
|
|
| Phase.Neg -> -. result
|
|
|
|
|
|
let _ = check 0 integral_value 1 1 2 1
|
|
|
|
let m_1121_2Ha_1F =
|
|
array_3_init mo_num mo_num mo_num (fun i j k ->
|
|
sum mos_cabs (fun a ->
|
|
h_two i j k a Spin.Alfa Spin.Alfa *. f_one a j Spin.Alfa))
|
|
|
|
|
|
let m_1121_2Hb_1F =
|
|
array_3_init mo_num mo_num mo_num (fun i j k ->
|
|
sum mos_cabs (fun a ->
|
|
h_two i j k a Spin.Alfa Spin.Beta *. f_one a j Spin.Beta))
|
|
|
|
|
|
let m_1121_2Ha_2Fa =
|
|
array_4_init mo_num mo_num mo_num mo_num (fun i j k m ->
|
|
sum mos_cabs (fun a ->
|
|
h_two i j k a Spin.Alfa Spin.Alfa *.
|
|
f_two m a m j Spin.Alfa Spin.Alfa
|
|
)
|
|
)
|
|
|
|
|
|
let m_1121_2Hb_2Fa =
|
|
array_4_init mo_num mo_num mo_num mo_num (fun i j k m ->
|
|
sum mos_cabs (fun a ->
|
|
h_two i j k a Spin.Alfa Spin.Beta *.
|
|
f_two m a m j Spin.Alfa Spin.Alfa
|
|
)
|
|
)
|
|
|
|
let m_1121_2Ha_2Fb =
|
|
array_4_init mo_num mo_num mo_num mo_num (fun i j k m ->
|
|
sum mos_cabs (fun a ->
|
|
h_two i j k a Spin.Alfa Spin.Alfa *.
|
|
f_two m a m j Spin.Alfa Spin.Beta
|
|
)
|
|
)
|
|
|
|
|
|
let m_1121_2Hb_2Fb =
|
|
array_4_init mo_num mo_num mo_num mo_num (fun i j k m ->
|
|
sum mos_cabs (fun a ->
|
|
h_two i j k a Spin.Alfa Spin.Beta *.
|
|
f_two m a m j Spin.Alfa Spin.Beta
|
|
)
|
|
)
|
|
|
|
|
|
let integral_value 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_j, mos_j' =
|
|
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 kj, mos_b kj
|
|
| Spin.Beta -> mos_b kj, mos_a kj
|
|
in
|
|
|
|
let result =
|
|
sum mos_j (fun j ->
|
|
m_1121_2Ha_1F.{i,j,k} +.
|
|
sum mos_i (fun m -> m_1121_2Ha_2Fa.{i,j,k,m}) +.
|
|
sum mos_i' (fun m -> m_1121_2Ha_2Fb.{i,j,k,m})) +.
|
|
sum mos_j' (fun j ->
|
|
m_1121_2Hb_1F.{i,j,k} +.
|
|
sum mos_i (fun m -> m_1121_2Hb_2Fb.{i,j,k,m}) +.
|
|
sum mos_i' (fun m -> m_1121_2Hb_2Fa.{i,j,k,m}))
|
|
in
|
|
|
|
match phase with
|
|
| Phase.Pos -> result
|
|
| Phase.Neg -> -. result
|
|
|
|
|
|
|
|
let _ = check 0 integral_value 1 1 2 1
|
|
|
|
let _ = png_image "1_1_22.png"
|
|
|
|
(*
|
|
let ki = det_I.(6)
|
|
let kj = det_I.(58)
|
|
|
|
let _ =
|
|
Format.printf "|I> -> |J> : %a |\n@." Excitation.pp (Excitation.of_det ki kj) ;
|
|
generate_alphas ki kj 1 1 2 2
|
|
|> Array.of_list
|
|
|> Array.mapi (fun kk alpha ->
|
|
|
|
(*
|
|
let _ = Determinant.to_lists ki
|
|
let _ = Determinant.to_lists alpha
|
|
let _ = Determinant.to_lists kj
|
|
|
|
|
|
let _ =
|
|
Format.printf "|I> -> |J> : %a\n@." Excitation.pp (Excitation.of_det ki kj);
|
|
Format.printf "|I> -> |a> : %a\n@." Excitation.pp (Excitation.of_det ki alpha);
|
|
Format.printf "|a> -> |J> : %a\n@." Excitation.pp (Excitation.of_det alpha kj)
|
|
*)
|
|
|
|
let integral_value ki kj =
|
|
|
|
let exc0 = Array.init (aux_num+1) (fun _ -> [|"-";"-"|]) in
|
|
let exc1 = Array.init (aux_num+1) (fun _ -> [|"-";"-"|]) in
|
|
let exc2 = Array.init (aux_num+1) (fun _ -> [|"-";"-"|]) in
|
|
|
|
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 spin = function
|
|
| Spin.Alfa -> 0
|
|
| _ -> 1
|
|
in
|
|
exc0.(i).(spin s ) <- "i"
|
|
; exc0.(k).(spin s ) <- "k"
|
|
;
|
|
let s0 = s in
|
|
|
|
let i, j, k, l, s, s', p1 =
|
|
match Excitation.of_det ki alpha with
|
|
| Excitation.(Double (phase, { hole ; particle ; spin }, {hole=hole' ; particle=particle' ; spin=spin' })) ->
|
|
hole, hole', particle, particle', spin, spin', phase
|
|
| _ -> assert false
|
|
in
|
|
|
|
if exc0.(i).(spin s ) = "-" then exc0.(i).(spin s ) <- "n";
|
|
if exc0.(j).(spin s') = "-" then exc0.(j).(spin s') <- "n";
|
|
if exc0.(k).(spin s ) = "-" then exc0.(k).(spin s ) <- if k > mo_num then "a" else "m";
|
|
if exc0.(l).(spin s') = "-" then exc0.(l).(spin s') <- if l > mo_num then "a" else "m";
|
|
|
|
let string_h =
|
|
Printf.sprintf "h_two %s %s %s %s %s %s *. "
|
|
exc0.(i).(spin s )
|
|
exc0.(j).(spin s')
|
|
exc0.(k).(spin s )
|
|
exc0.(l).(spin s')
|
|
(if s = s0 then "s " else "s'")
|
|
(if s' = s0 then "s " else "s'")
|
|
(*
|
|
(if exc0.(i).(spin s ) = "i" || exc0.(k).(spin s ) = "k" then "s " else
|
|
if exc0.(i).(spin s ) = "j" || exc0.(k).(spin s ) = "l" then "s' " else "s''")
|
|
(if exc0.(j).(spin s') = "i" || exc0.(l).(spin s') = "k" then "s " else
|
|
if exc0.(j).(spin s') = "j" || exc0.(l).(spin s') = "l" then "s' " else "s''")
|
|
*)
|
|
in
|
|
|
|
let i, j, k, l, s, s', p2 =
|
|
match Excitation.of_det alpha 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 string_f =
|
|
Printf.sprintf "f_two %s %s %s %s %s %s"
|
|
exc0.(i).(spin s )
|
|
exc0.(j).(spin s')
|
|
exc0.(k).(spin s )
|
|
exc0.(l).(spin s')
|
|
(*
|
|
(if exc0.(i).(spin s ) = "i" || exc0.(k).(spin s ) = "k" then "s " else
|
|
if exc0.(i).(spin s ) = "j" || exc0.(k).(spin s ) = "l" then "s' " else
|
|
if exc0.(i).(spin s ) = "n" || exc0.(k).(spin s ) = "n" then "s''" else
|
|
if s = s0 then "s" else "s'")
|
|
(if exc0.(j).(spin s') = "i" || exc0.(l).(spin s') = "k" then "s " else
|
|
if exc0.(j).(spin s') = "j" || exc0.(l).(spin s') = "l" then "s' " else
|
|
if exc0.(j).(spin s') = "n" || exc0.(l).(spin s') = "n" then "s''" else
|
|
if s' = s0 then "s" else "s'")
|
|
*)
|
|
(if s = s0 then "s " else "s'")
|
|
(if s' = s0 then "s " else "s'")
|
|
in
|
|
|
|
Format.printf "|I> -> |a> : %a | %s\n@." Excitation.pp (Excitation.of_det ki alpha) string_h ;
|
|
Format.printf "|a> -> |J> : %a | %s\n@." Excitation.pp (Excitation.of_det alpha kj) string_f ;
|
|
|
|
let sign =
|
|
if Phase.add p1 p2 = phase then "+. " else "-. "
|
|
in
|
|
sign ^ string_h ^ string_f
|
|
(*
|
|
|
|
let mos, mos' =
|
|
match s with
|
|
| Spin.Alfa -> mos_a ki, mos_b ki
|
|
| Spin.Beta -> mos_b ki, mos_a ki
|
|
in
|
|
|
|
let result =
|
|
let k=4 and l=5 in
|
|
let i=1 and j=2 in
|
|
let n=2 and a=7 in
|
|
let s = Spin.Alfa
|
|
and s' = Spin.Alfa
|
|
and s'' = Spin.Beta
|
|
in
|
|
(*
|
|
h_two j n l a s' s'' *. f_two a i n k s s''
|
|
*)
|
|
h_two i n k a s s'' *. f_two j a l n s s''
|
|
in
|
|
|
|
match phase with
|
|
| Phase.Pos -> result
|
|
| Phase.Neg -> -. result
|
|
*)
|
|
|
|
in
|
|
integral_value ki kj
|
|
(*
|
|
let a = (compute_HaaF ki [alpha] kj)
|
|
and b = (integral_value ki kj)
|
|
in
|
|
if kk = 31 then
|
|
Format.printf "%6d %e %e@.@." (kk) a b
|
|
*)
|
|
)
|
|
|> Array.to_list
|
|
|> List.sort_uniq compare
|
|
|
|
let _ = compute_HaaF ki (generate_alphas ki kj 1 1 2 2) kj
|
|
|
|
let _ = integral_value ki kj
|
|
|
|
*)
|
|
|
|
let integral_value 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
|
|
|
|
(* MOs unoccupied in both I and J *)
|
|
let mos_virt_a, mos_virt_b =
|
|
Array.init mo_num (fun i -> Some (i+1)) ,
|
|
Array.init mo_num (fun i -> Some (i+1))
|
|
in
|
|
List.iter (fun i -> mos_virt_a.(i-1) <- None) (mos_a ki);
|
|
List.iter (fun i -> mos_virt_a.(i-1) <- None) (mos_a kj);
|
|
List.iter (fun i -> mos_virt_b.(i-1) <- None) (mos_b ki);
|
|
List.iter (fun i -> mos_virt_b.(i-1) <- None) (mos_b kj);
|
|
|
|
let mos_virt_a, mos_virt_b =
|
|
Array.to_list mos_virt_a |> Util.list_some,
|
|
Array.to_list mos_virt_b |> Util.list_some
|
|
in
|
|
|
|
let mos_virt, mos_virt' =
|
|
match s with
|
|
| Spin.Alfa -> mos_virt_a, mos_virt_b
|
|
| Spin.Beta -> mos_virt_b, mos_virt_a
|
|
in
|
|
|
|
let mos_j, mos_j' =
|
|
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 result =
|
|
let s' = Spin.other s in
|
|
sum mos_cabs (fun a ->
|
|
sum mos_j (fun n ->
|
|
sum mos_virt (fun m ->
|
|
h_two n i a m s s *. f_two m a k n s s
|
|
))) +.
|
|
sum mos_cabs (fun a' ->
|
|
sum mos_j' (fun n' ->
|
|
sum mos_virt (fun m ->
|
|
h_two i n' m a' s s' *. f_two m a' k n' s s'
|
|
))) +.
|
|
sum mos_cabs (fun a ->
|
|
sum mos_j' (fun n' ->
|
|
sum mos_virt' (fun m' ->
|
|
h_two i n' a m' s s' *. f_two a m' k n' s s'
|
|
))) -.
|
|
sum mos_cabs (fun a ->
|
|
sum mos_j (fun n ->
|
|
sum mos_j (fun l -> if l <=n then 0. else
|
|
h_two n l a k s s *. f_two i a l n s s
|
|
))) -.
|
|
sum mos_cabs (fun a' ->
|
|
sum mos_j' (fun n' ->
|
|
sum mos_j (fun l ->
|
|
h_two l n' k a' s s' *. f_two i a' l n' s s'
|
|
)))
|
|
in
|
|
match phase with
|
|
| Phase.Pos -> result
|
|
| Phase.Neg -> -. result
|
|
|
|
let _ = check 0 integral_value 1 1 2 2
|
|
|
|
|
|
|
|
let m_1122_va =
|
|
array_4_init mo_num mo_num mo_num mo_num (fun i j k m ->
|
|
sum mos_cabs (fun a ->
|
|
h_two j i a m Spin.Alfa Spin.Alfa *.
|
|
f_two m a k j Spin.Alfa Spin.Alfa
|
|
)
|
|
)
|
|
|
|
let m_1122_v2 =
|
|
array_4_init mo_num mo_num mo_num mo_num (fun i j k m ->
|
|
sum mos_cabs (fun a ->
|
|
h_two i j m a Spin.Alfa Spin.Beta *.
|
|
f_two m a k j Spin.Alfa Spin.Beta
|
|
)
|
|
)
|
|
|
|
let m_1122_v3 =
|
|
array_4_init mo_num mo_num mo_num mo_num (fun i j k m ->
|
|
sum mos_cabs (fun a ->
|
|
h_two i j a m Spin.Alfa Spin.Beta *.
|
|
f_two a m k j Spin.Alfa Spin.Beta
|
|
)
|
|
)
|
|
|
|
let m_1122_oa =
|
|
array_4_init mo_num mo_num mo_num mo_num (fun i j k l ->
|
|
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
|
|
)
|
|
)
|
|
|
|
let m_1122_o =
|
|
array_4_init mo_num mo_num mo_num mo_num (fun i j k l ->
|
|
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
|
|
)
|
|
)
|
|
|
|
|
|
let integral_value 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
|
|
|
|
(* MOs unoccupied in both I and J *)
|
|
let mos_virt_a, mos_virt_b =
|
|
Array.init mo_num (fun i -> Some (i+1)) ,
|
|
Array.init mo_num (fun i -> Some (i+1))
|
|
in
|
|
List.iter (fun i -> mos_virt_a.(i-1) <- None) (mos_a ki);
|
|
List.iter (fun i -> mos_virt_a.(i-1) <- None) (mos_a kj);
|
|
List.iter (fun i -> mos_virt_b.(i-1) <- None) (mos_b ki);
|
|
List.iter (fun i -> mos_virt_b.(i-1) <- None) (mos_b kj);
|
|
|
|
let mos_virt_a, mos_virt_b =
|
|
Array.to_list mos_virt_a |> Util.list_some,
|
|
Array.to_list mos_virt_b |> Util.list_some
|
|
in
|
|
|
|
let mos_virt, mos_virt' =
|
|
match s with
|
|
| Spin.Alfa -> mos_virt_a, mos_virt_b
|
|
| Spin.Beta -> mos_virt_b, mos_virt_a
|
|
in
|
|
|
|
let mos_j, mos_j' =
|
|
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 result =
|
|
sum mos_j (fun j ->
|
|
sum mos_virt (fun m ->
|
|
m_1122_va.{i,j,k,m}
|
|
)) +.
|
|
sum mos_j' (fun j ->
|
|
sum mos_virt (fun m ->
|
|
m_1122_v2.{i,j,k,m}
|
|
)) +.
|
|
sum mos_j' (fun j ->
|
|
sum mos_virt' (fun m ->
|
|
m_1122_v3.{i,j,k,m}
|
|
)) -.
|
|
sum mos_j (fun j ->
|
|
sum mos_j (fun l -> if l <= j then 0. else
|
|
m_1122_oa.{i,j,k,l}
|
|
)) -.
|
|
sum mos_j' (fun j ->
|
|
sum mos_j (fun l ->
|
|
m_1122_o.{i,j,k,l}
|
|
))
|
|
in
|
|
match phase with
|
|
| Phase.Pos -> result
|
|
| Phase.Neg -> -. result
|
|
|
|
|
|
let _ = check 0 integral_value 1 1 2 2
|
|
|
|
|
|
|
|
|
|
let _ = png_image "1_2_22.png"
|
|
|
|
let integral_value ki kj =
|
|
let h, p, s, phase =
|
|
match Excitation.of_det ki kj with
|
|
| Excitation.(Single (phase, { hole ; particle ; spin })) ->
|
|
hole, particle, spin, phase
|
|
| _ -> assert false
|
|
in
|
|
|
|
let mos, mos' =
|
|
match s with
|
|
| Spin.Alfa -> mos_a ki, mos_b ki
|
|
| Spin.Beta -> mos_b ki, mos_a ki
|
|
in
|
|
|
|
|
|
|
|
let result =
|
|
(* Alpha-Beta *)
|
|
let s' = Spin.other s in
|
|
sum mos_cabs (fun b ->
|
|
sum mos_cabs (fun a ->
|
|
sum mos' (fun j -> h_two h j a b s s' *. f_two a b p j s s')
|
|
))
|
|
+.
|
|
(* Alpha-Alpha / Beta-Beta *)
|
|
sum mos_cabs (fun b ->
|
|
sum mos_cabs (fun a -> if b >= a then 0. else
|
|
sum mos (fun j -> h_two h j a b s s *. f_two a b p j s s)
|
|
))
|
|
in
|
|
match phase with
|
|
| Phase.Pos -> result
|
|
| Phase.Neg -> -. result
|
|
|
|
|
|
let _ = check 0 integral_value 1 2 2 2
|
|
|
|
let m_1222a =
|
|
array_3_init mo_num mo_num mo_num (fun i j k ->
|
|
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 k j Spin.Alfa Spin.Alfa
|
|
)
|
|
)
|
|
)
|
|
|
|
let m_1222 =
|
|
array_3_init mo_num mo_num mo_num (fun i j 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 k j Spin.Alfa Spin.Beta
|
|
)
|
|
)
|
|
)
|
|
|
|
|
|
let integral_value 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
|
|
|
|
(* MOs unoccupied in both I and J *)
|
|
let mos_j, mos_j' =
|
|
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 result =
|
|
sum mos_j (fun j -> m_1222a.{i,j,k} ) +.
|
|
sum mos_j' (fun j -> m_1222.{i,j,k} )
|
|
in
|
|
match phase with
|
|
| Phase.Pos -> result
|
|
| Phase.Neg -> -. result
|
|
|
|
|
|
let _ = check 0 integral_value 1 2 2 2
|
|
|
|
|
|
|
|
|
|
let _ = png_image "2_1_12.png"
|
|
|
|
let integral_value 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, mos' =
|
|
match s with
|
|
| Spin.Alfa -> mos_a ki, mos_b ki
|
|
| Spin.Beta -> mos_b ki, mos_a ki
|
|
in
|
|
|
|
let mos2, mos2' =
|
|
match s' with
|
|
| Spin.Alfa -> mos_a ki, mos_b ki
|
|
| Spin.Beta -> mos_b ki, mos_a ki
|
|
in
|
|
|
|
let result =
|
|
sum mos_cabs (fun a ->
|
|
let s'' = Spin.other s in
|
|
( h_one i a s +.
|
|
sum mos (fun n -> h_two i n a n s s ) +.
|
|
sum mos' (fun n -> h_two i n a n s s'')
|
|
) *. f_two a j k l s s'
|
|
) +.
|
|
sum mos_cabs (fun a ->
|
|
let s'' = Spin.other s' in
|
|
( h_one j a s' +.
|
|
sum mos2 (fun n -> h_two j n a n s' s' ) +.
|
|
sum mos2' (fun n -> h_two j n a n s' s'')
|
|
) *. f_two i a k l s s'
|
|
)
|
|
in
|
|
|
|
match phase with
|
|
| Phase.Pos -> result
|
|
| Phase.Neg -> -. result
|
|
|
|
|
|
let _ = check 0 integral_value 2 1 1 2
|
|
|
|
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 )
|
|
)
|
|
|
|
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)
|
|
)
|
|
|
|
let m_2112_2Ha_2Fa =
|
|
array_5_init mo_num mo_num mo_num mo_num mo_num (fun i j k l n ->
|
|
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 )
|
|
)
|
|
|
|
let m_2112_2Hb_2Fa =
|
|
array_5_init mo_num mo_num mo_num mo_num mo_num (fun i j k l n ->
|
|
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 )
|
|
)
|
|
|
|
let m_2112_2Ha_2Fb =
|
|
array_5_init mo_num mo_num mo_num mo_num mo_num (fun i j k l n ->
|
|
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 )
|
|
)
|
|
|
|
let m_2112_2Hb_2Fb =
|
|
array_5_init mo_num mo_num mo_num mo_num mo_num (fun i j k l n ->
|
|
sum mos_cabs (fun a ->
|
|
h_two i n a n Spin.Alfa Spin.Beta *.
|
|
f_two a j k l Spin.Alfa Spin.Beta )
|
|
)
|
|
|
|
let integral_value 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, mos' =
|
|
match s with
|
|
| Spin.Alfa -> mos_a ki, mos_b ki
|
|
| Spin.Beta -> mos_b ki, mos_a ki
|
|
in
|
|
|
|
let result =
|
|
if s = s' then
|
|
m_2112_1H_2Fa.{i,j,k,l} +.
|
|
sum mos (fun n ->
|
|
m_2112_2Ha_2Fa.{i,j,k,l,n} +. m_2112_2Ha_2Fa.{j,i,l,k,n}
|
|
) +.
|
|
sum mos' (fun n ->
|
|
m_2112_2Hb_2Fa.{i,j,k,l,n} +. m_2112_2Hb_2Fa.{j,i,l,k,n}
|
|
)
|
|
else
|
|
m_2112_1H_2Fb.{i,j,k,l} +.
|
|
sum mos (fun n ->
|
|
m_2112_2Ha_2Fb.{i,j,k,l,n} +. m_2112_2Hb_2Fb.{j,i,l,k,n}
|
|
) +.
|
|
sum mos' (fun n ->
|
|
m_2112_2Hb_2Fb.{i,j,k,l,n} +. m_2112_2Ha_2Fb.{j,i,l,k,n}
|
|
)
|
|
in
|
|
|
|
match phase with
|
|
| Phase.Pos -> result
|
|
| Phase.Neg -> -. result
|
|
|
|
|
|
let _ = check 0 integral_value 2 1 1 2
|
|
|
|
|
|
|
|
|
|
let _ = png_image "2_1_21.png"
|
|
|
|
let integral_value 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 result =
|
|
|
|
let mos, mos', s'' =
|
|
match s with
|
|
| Spin.Alfa -> mos_a kj, mos_b kj, Spin.Beta
|
|
| Spin.Beta -> mos_b kj, mos_a kj, Spin.Alfa
|
|
in
|
|
sum mos_cabs (fun a ->
|
|
h_two i j a l s s' *.
|
|
( f_one a k s +.
|
|
sum mos (fun n -> f_two a n k n s s) +.
|
|
sum mos' (fun n -> f_two a n k n s s'')
|
|
) ) +.
|
|
|
|
let mos, mos', s'' =
|
|
match s' with
|
|
| Spin.Alfa -> mos_a kj, mos_b kj, Spin.Beta
|
|
| Spin.Beta -> mos_b kj, mos_a kj, Spin.Alfa
|
|
in
|
|
sum mos_cabs (fun a ->
|
|
h_two j i a k s' s *.
|
|
( f_one a l s' +.
|
|
sum mos (fun n -> f_two a n l n s' s') +.
|
|
sum mos' (fun n -> f_two a n l n s' s'')
|
|
)
|
|
)
|
|
in
|
|
|
|
match phase with
|
|
| Phase.Pos -> result
|
|
| Phase.Neg -> -. result
|
|
|
|
|
|
let _ = check 0 integral_value 2 1 2 1
|
|
|
|
let m_2121_2Ha_1F =
|
|
array_4_init mo_num mo_num mo_num mo_num (fun i j k l ->
|
|
sum mos_cabs (fun a ->
|
|
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)
|
|
)
|
|
|
|
let m_2121_2Hb_1F =
|
|
array_4_init mo_num mo_num mo_num mo_num (fun i j k l ->
|
|
sum mos_cabs (fun a ->
|
|
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)
|
|
)
|
|
|
|
let m_2121_2Ha_2Fa =
|
|
array_5_init mo_num mo_num mo_num mo_num mo_num (fun i j k l n ->
|
|
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 )
|
|
)
|
|
|
|
let m_2121_2Hb_2Fa =
|
|
array_5_init mo_num mo_num mo_num mo_num mo_num (fun i j k l n ->
|
|
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 )
|
|
)
|
|
|
|
let m_2121_2Ha_2Fb =
|
|
array_5_init mo_num mo_num mo_num mo_num mo_num (fun i j k l n ->
|
|
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 )
|
|
)
|
|
|
|
let m_2121_2Hb_2Fb =
|
|
array_5_init mo_num mo_num mo_num mo_num mo_num (fun i j k l n ->
|
|
sum mos_cabs (fun a ->
|
|
h_two i j a l Spin.Alfa Spin.Beta *.
|
|
f_two a n k n Spin.Alfa Spin.Beta )
|
|
)
|
|
|
|
let integral_value 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, mos' =
|
|
match s with
|
|
| Spin.Alfa -> mos_a kj, mos_b kj
|
|
| Spin.Beta -> mos_b kj, mos_a kj
|
|
in
|
|
|
|
let result =
|
|
if s = s' then
|
|
m_2121_2Ha_1F.{i,j,k,l} +.
|
|
sum mos (fun n ->
|
|
m_2121_2Ha_2Fa.{i,j,k,l,n} +. m_2121_2Ha_2Fa.{j,i,l,k,n}
|
|
) +.
|
|
sum mos' (fun n ->
|
|
m_2121_2Ha_2Fb.{i,j,k,l,n} +. m_2121_2Ha_2Fb.{j,i,l,k,n}
|
|
)
|
|
else
|
|
m_2121_2Hb_1F.{i,j,k,l} +.
|
|
sum mos (fun n ->
|
|
m_2121_2Hb_2Fa.{i,j,k,l,n} +. m_2121_2Hb_2Fb.{j,i,l,k,n}
|
|
) +.
|
|
sum mos' (fun n ->
|
|
m_2121_2Hb_2Fb.{i,j,k,l,n} +. m_2121_2Hb_2Fa.{j,i,l,k,n}
|
|
)
|
|
in
|
|
|
|
match phase with
|
|
| Phase.Pos -> result
|
|
| Phase.Neg -> -. result
|
|
|
|
|
|
let _ = check 0 integral_value 2 1 2 1
|
|
|
|
|
|
|
|
|
|
let _ = png_image "2_1_22.png"
|
|
|
|
let integral_value 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_virt_a, mos_virt_b =
|
|
Array.init mo_num (fun i -> Some (i+1)) ,
|
|
Array.init mo_num (fun i -> Some (i+1))
|
|
in
|
|
List.iter (fun i -> mos_virt_a.(i-1) <- None) (mos_a ki);
|
|
List.iter (fun i -> mos_virt_a.(i-1) <- None) (mos_a kj);
|
|
List.iter (fun i -> mos_virt_b.(i-1) <- None) (mos_b ki);
|
|
List.iter (fun i -> mos_virt_b.(i-1) <- None) (mos_b kj);
|
|
|
|
let mos_virt_a, mos_virt_b =
|
|
Array.to_list mos_virt_a |> Util.list_some,
|
|
Array.to_list mos_virt_b |> Util.list_some
|
|
in
|
|
|
|
let result =
|
|
|
|
let mos_virt, mos_virt' =
|
|
match s with
|
|
| Spin.Alfa -> mos_virt_a, mos_virt_b
|
|
| Spin.Beta -> mos_virt_b, mos_virt_a
|
|
in
|
|
|
|
let mos, mos' =
|
|
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
|
|
|
|
if s = s' then
|
|
let s'' = Spin.other s' in
|
|
sum mos_cabs (fun a ->
|
|
sum mos (fun n ->
|
|
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
|
|
)
|
|
+. sum mos_virt (fun m ->
|
|
-. h_two i j a m s s *. f_two m a k l s s )
|
|
+. sum mos' (fun n ->
|
|
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''
|
|
)
|
|
)
|
|
else
|
|
sum mos_cabs (fun a ->
|
|
sum mos_virt' (fun m ->
|
|
h_two i j a m s s' *. f_two a m k l s s' ) +.
|
|
sum mos_virt (fun m ->
|
|
h_two i j m a s s' *. f_two m a k l s s' ) +.
|
|
sum mos (fun n ->
|
|
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'
|
|
) +.
|
|
sum mos' (fun n -> if n >= j then 0. else
|
|
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' ) +.
|
|
sum mos' (fun n -> if n <= j then 0. 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' ) +.
|
|
sum mos' (fun n ->
|
|
-. h_two i n a l s s' *. f_two a j k n s s'
|
|
)
|
|
)
|
|
in
|
|
|
|
match phase with
|
|
| Phase.Pos -> result
|
|
| Phase.Neg -> -. result
|
|
|
|
|
|
let _ = check 0 integral_value 2 1 2 2
|
|
|
|
let _ = png_image "2_2_22.png"
|
|
|
|
let integral_value ki kj =
|
|
let h, h', p, p', s, s', phase =
|
|
match Excitation.of_det ki kj with
|
|
| Excitation.(Double (phase,
|
|
{ hole=h ; particle=p ; spin=s },
|
|
{ hole=h'; particle=p'; spin=s'}) ) -> h, h', p, p', s, s', phase
|
|
| _ -> assert false
|
|
in
|
|
|
|
let result =
|
|
if s <> s' then (* Alpha-Beta *)
|
|
sum mos_cabs (fun b ->
|
|
sum mos_cabs (fun a ->
|
|
h_two h h' a b s s' *. f_two a b p p' s s'
|
|
))
|
|
else (* Alpha-Alpha / Beta-Beta *)
|
|
sum mos_cabs (fun b ->
|
|
sum mos_cabs (fun a -> if b >= a then 0. else
|
|
h_two h h' a b s s' *. f_two a b p p' s s'
|
|
))
|
|
in
|
|
match phase with
|
|
| Phase.Pos -> result
|
|
| Phase.Neg -> -. result
|
|
|
|
|
|
let _ = check 0 integral_value 2 2 2 2
|
|
|
|
let _ = png_image "3_1_22.png"
|
|
|
|
let integral_value 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 -> 0. (*TODO *)
|
|
in
|
|
match phase with
|
|
| Phase.Pos -> result
|
|
| Phase.Neg -> -. result
|
|
|
|
|
|
let _ = check 0 integral_value 3 1 2 2
|
|
|
|
let ki = det_I.(129)
|
|
let kj = det_I.(349)
|
|
|
|
let alpha_to_string alpha =
|
|
|
|
let exc0 = Array.init (aux_num+1) (fun _ -> [|"-";"-"|]) in
|
|
|
|
let i, j, m, k, l, n, s, s', s'', phase =
|
|
match Excitation.of_det ki kj with
|
|
| Excitation.(Triple (phase,
|
|
{ hole ; particle ; spin },
|
|
{hole=hole' ; particle=particle' ; spin=spin' },
|
|
{hole=hole''; particle=particle''; spin=spin''} )) ->
|
|
hole, hole', hole'', particle, particle', particle'', spin, spin', spin'', phase
|
|
| _ -> assert false
|
|
in
|
|
let spin = function
|
|
| Spin.Alfa -> 0
|
|
| _ -> 1
|
|
in
|
|
exc0.(i).(spin s ) <- "i"
|
|
; exc0.(j).(spin s' ) <- "j"
|
|
; exc0.(k).(spin s ) <- "k"
|
|
; exc0.(l).(spin s' ) <- "l"
|
|
; exc0.(m).(spin s'') <- "m"
|
|
; exc0.(n).(spin s'') <- "n"
|
|
;
|
|
let s0, s0', s0'' = s, s', s'' in
|
|
|
|
let i, j, k, l, s, s', p1 =
|
|
match Excitation.of_det ki alpha with
|
|
| Excitation.(Double (phase, { hole ; particle ; spin }, {hole=hole' ; particle=particle' ; spin=spin' })) ->
|
|
hole, hole', particle, particle', spin, spin', phase
|
|
| _ -> assert false
|
|
in
|
|
|
|
if exc0.(i).(spin s ) = "-" then exc0.(i).(spin s ) <- "p";
|
|
if exc0.(j).(spin s') = "-" then exc0.(j).(spin s') <- "p";
|
|
if exc0.(k).(spin s ) = "-" then exc0.(k).(spin s ) <- if k > mo_num then "a" else "q";
|
|
if exc0.(l).(spin s') = "-" then exc0.(l).(spin s') <- if l > mo_num then "a" else "q";
|
|
|
|
let string_h =
|
|
Printf.sprintf "h_two %s %s %s %s %s %s *. "
|
|
exc0.(i).(spin s )
|
|
exc0.(j).(spin s')
|
|
exc0.(k).(spin s )
|
|
exc0.(l).(spin s')
|
|
(if s = s0 then "s " else if s = s0' then "s'" else "s''")
|
|
(if s' = s0' then "s'" else if s = s0'' then "s''" else "s" )
|
|
in
|
|
|
|
let i, j, k, l, s, s', p2 =
|
|
match Excitation.of_det alpha 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 string_f =
|
|
Printf.sprintf "f_two %s %s %s %s %s %s"
|
|
exc0.(i).(spin s )
|
|
exc0.(j).(spin s')
|
|
exc0.(k).(spin s )
|
|
exc0.(l).(spin s')
|
|
(if s = s0 then "s " else if s = s0' then "s'" else "s''")
|
|
(if s' = s0' then "s'" else if s = s0'' then "s''" else "s" )
|
|
in
|
|
(*
|
|
Format.printf "|I> -> |a> : %a | %s\n@." Excitation.pp (Excitation.of_det ki alpha) string_h ;
|
|
Format.printf "|a> -> |J> : %a | %s\n@." Excitation.pp (Excitation.of_det alpha kj) string_f ;
|
|
*)
|
|
|
|
|
|
let sign =
|
|
if Phase.add p1 p2 = phase then "+." else "-."
|
|
in
|
|
sign ^ string_h ^ string_f
|
|
|
|
|
|
let alpha_debug alpha =
|
|
|
|
let i, j, k, l, s, s', p1 =
|
|
match Excitation.of_det ki alpha with
|
|
| Excitation.(Double (phase, { hole ; particle ; spin }, {hole=hole' ; particle=particle' ; spin=spin' })) ->
|
|
hole, hole', particle, particle', spin, spin', phase
|
|
| _ -> assert false
|
|
in
|
|
Printf.printf "%d %d %d %d " i j k l;
|
|
|
|
let i, j, k, l, s, s', p2 =
|
|
match Excitation.of_det alpha kj with
|
|
| Excitation.(Double (phase, { hole ; particle ; spin }, {hole=hole' ; particle=particle' ; spin=spin' })) ->
|
|
hole, hole', particle, particle', spin, spin', phase
|
|
| _ -> assert false
|
|
in
|
|
(*
|
|
Format.printf "|I> -> |a> : %a | \n@." Excitation.pp (Excitation.of_det ki alpha) ;
|
|
Format.printf "|a> -> |J> : %a | \n@." Excitation.pp (Excitation.of_det alpha kj) ;
|
|
*)
|
|
Printf.printf "%d %d %d %d \n%!" i j k l
|
|
|
|
|
|
let strings =
|
|
Format.printf "|I> -> |J> : %a |\n@." Excitation.pp (Excitation.of_det ki kj) ;
|
|
generate_alphas ki kj 3 1 2 2
|
|
|> Array.of_list
|
|
|> Array.mapi (fun kk alpha -> alpha_to_string alpha)
|
|
|> Array.to_list
|
|
|> List.sort_uniq compare
|
|
|> Array.of_list
|
|
|
|
let _ = Array.iteri (fun i x -> Printf.printf "%d %s \n%!" i x) strings
|
|
|
|
let _ =
|
|
let v =
|
|
let alphas =
|
|
generate_alphas ki kj 3 1 2 2
|
|
(*
|
|
|> List.filter (fun alpha ->
|
|
let x = alpha_to_string alpha in
|
|
x = strings.(6)
|
|
)
|
|
*)
|
|
in
|
|
(*
|
|
List.iter alpha_debug alphas ;
|
|
Printf.printf "\n%!";
|
|
*)
|
|
compute_HaaF ki alphas kj
|
|
in
|
|
let x = (integral_value ki kj) in
|
|
Printf.printf "%20.8e %20.8e %20.8e\n%!" x v (v-. x)
|
|
|
|
|
|
|