mirror of
https://gitlab.com/scemama/QCaml.git
synced 2024-08-30 16:23:41 +02:00
148 lines
4.8 KiB
OCaml
148 lines
4.8 KiB
OCaml
open Lacaml.D
|
|
|
|
module De = Determinant
|
|
module Ex = Excitation
|
|
module Sp = Spindeterminant
|
|
|
|
type t = float list
|
|
|
|
|
|
let non_zero integrals degree_a degree_b ki kj =
|
|
let kia = De.alfa ki and kib = De.beta ki
|
|
and kja = De.alfa kj and kjb = De.beta kj
|
|
in
|
|
let diag_element =
|
|
let mo_a = Sp.to_list kia
|
|
and mo_b = Sp.to_list kib
|
|
in
|
|
fun one_e two_e ->
|
|
let one =
|
|
(List.fold_left (fun accu i -> accu +. one_e i i Spin.Alfa) 0. mo_a)
|
|
+.
|
|
(List.fold_left (fun accu i -> accu +. one_e i i Spin.Beta) 0. mo_b)
|
|
in
|
|
let two =
|
|
let rec aux_same spin accu = function
|
|
| [] -> accu
|
|
| i :: rest ->
|
|
let new_accu =
|
|
List.fold_left (fun accu j -> accu +. two_e i j i j spin spin) accu rest
|
|
in
|
|
(aux_same [@tailcall]) spin new_accu rest
|
|
in
|
|
let rec aux_opposite accu other = function
|
|
| [] -> accu
|
|
| i :: rest ->
|
|
let new_accu =
|
|
List.fold_left (fun accu j -> accu +. two_e i j i j Spin.Alfa Spin.Beta) accu other
|
|
in
|
|
(aux_opposite [@tailcall]) new_accu other rest
|
|
in
|
|
(aux_same Spin.Alfa 0. mo_a) +. (aux_same Spin.Beta 0. mo_b) +.
|
|
(aux_opposite 0. mo_a mo_b)
|
|
in
|
|
one +. two
|
|
in
|
|
let result =
|
|
match degree_a, degree_b with
|
|
| 1, 1 -> (* alpha-beta double *)
|
|
begin
|
|
let ha, pa, phase_a = Ex.single_of_spindet kia kja in
|
|
let hb, pb, phase_b = Ex.single_of_spindet kib kjb in
|
|
let s1 =
|
|
match phase_a, phase_b with
|
|
| Phase.Pos, Phase.Pos
|
|
| Phase.Neg, Phase.Neg -> fun _ two_e -> two_e ha hb pa pb Spin.Alfa Spin.Beta
|
|
| Phase.Neg, Phase.Pos
|
|
| Phase.Pos, Phase.Neg -> fun _ two_e -> -. two_e ha hb pa pb Spin.Alfa Spin.Beta
|
|
in
|
|
let kk = Determinant.double_excitation Spin.Alfa ha pb Spin.Beta hb pa ki in
|
|
let kka = De.alfa kk and kkb = De.beta kk in
|
|
match Spindeterminant.(degree kia kka, degree kib kkb) with
|
|
| (1,1) ->
|
|
let s2 =
|
|
begin
|
|
let ha, pa, phase_a = Ex.single_of_spindet kia kka in
|
|
let hb, pb, phase_b = Ex.single_of_spindet kib kkb in
|
|
match phase_a, phase_b with
|
|
| Phase.Pos, Phase.Pos
|
|
| Phase.Neg, Phase.Neg -> fun _ two_e -> two_e ha hb pa pb Spin.Alfa Spin.Beta
|
|
| Phase.Neg, Phase.Pos
|
|
| Phase.Pos, Phase.Neg -> fun _ two_e -> -. two_e ha hb pa pb Spin.Alfa Spin.Beta
|
|
end
|
|
in fun x two_e -> 0.75 *. (s1 x two_e) +. 0.25 *. (s2 x two_e)
|
|
| _ -> fun x two_e -> 0.
|
|
end
|
|
|
|
| 2, 0 -> (* alpha double *)
|
|
begin
|
|
let h1, p1, h2, p2, phase = Ex.double_of_spindet kia kja in
|
|
match phase with
|
|
| Phase.Pos -> fun _ two_e -> two_e h1 h2 p1 p2 Spin.Alfa Spin.Alfa
|
|
| Phase.Neg -> fun _ two_e -> -. two_e h1 h2 p1 p2 Spin.Alfa Spin.Alfa
|
|
end
|
|
|
|
| 0, 2 -> (* beta double *)
|
|
begin
|
|
let h1, p1, h2, p2, phase = Ex.double_of_spindet kib kjb in
|
|
match phase with
|
|
| Phase.Pos -> fun _ two_e -> two_e h1 h2 p1 p2 Spin.Beta Spin.Beta
|
|
| Phase.Neg -> fun _ two_e -> -. two_e h1 h2 p1 p2 Spin.Beta Spin.Beta
|
|
end
|
|
|
|
| 1, 0 (* alpha single *)
|
|
| 0, 1 (* beta single *)
|
|
-> fun _ _ -> 0.
|
|
| 0, 0 -> (* diagonal element *)
|
|
diag_element
|
|
|
|
| _ -> assert false
|
|
in
|
|
List.map (fun (one_e, two_e) -> result one_e two_e) integrals
|
|
|
|
|
|
let make integrals ki kj =
|
|
let degree_a, degree_b =
|
|
De.degrees ki kj
|
|
in
|
|
if degree_a+degree_b > 2 then
|
|
List.map (fun _ -> 0.) integrals
|
|
else
|
|
non_zero integrals degree_a degree_b ki kj
|
|
|
|
|
|
|
|
|
|
let make_s2 ki kj =
|
|
let degree_a = De.degree_alfa ki kj in
|
|
let kia = De.alfa ki in
|
|
let kja = De.alfa kj in
|
|
if degree_a > 1 then 0.
|
|
else
|
|
let degree_b = De.degree_beta ki kj in
|
|
let kib = De.beta ki in
|
|
let kjb = De.beta kj in
|
|
match degree_a, degree_b with
|
|
| 1, 1 -> (* alpha-beta double *)
|
|
let ha, pa, phase_a = Ex.single_of_spindet kia kja in
|
|
let hb, pb, phase_b = Ex.single_of_spindet kib kjb in
|
|
if ha = pb && hb = pa then
|
|
begin
|
|
match phase_a, phase_b with
|
|
| Phase.Pos, Phase.Pos
|
|
| Phase.Neg, Phase.Neg -> -1.
|
|
| Phase.Neg, Phase.Pos
|
|
| Phase.Pos, Phase.Neg -> 1.
|
|
end
|
|
else 0.
|
|
| 0, 0 ->
|
|
let ba = Sp.bitstring kia and bb = Sp.bitstring kib in
|
|
let tmp = Bitstring.logxor ba bb in
|
|
let n_a = Bitstring.logand ba tmp |> Bitstring.popcount in
|
|
let n_b = Bitstring.logand bb tmp |> Bitstring.popcount in
|
|
let s_z = 0.5 *. float_of_int (n_a - n_b) in
|
|
float_of_int n_a +. s_z *. (s_z -. 1.)
|
|
| _ -> 0.
|
|
|
|
|