10
1
mirror of https://gitlab.com/scemama/QCaml.git synced 2024-10-31 19:23:40 +01:00
QCaml/CI/Excitation.ml

189 lines
5.8 KiB
OCaml

type single_exc =
{
hole : int ;
particle : int ;
spin : Spin.t ;
}
type t =
| Identity of Phase.t
| Single of Phase.t * single_exc
| Double of Phase.t * single_exc * single_exc
| Multiple of Phase.t * single_exc list
let single_of_spindet t t' =
assert (Spindeterminant.degree t t' = 1);
let d = Spindeterminant.bitstring t
and d' = Spindeterminant.bitstring t'
in
let tmp = Bitstring.logxor d d' in
let shift_left_one = Bitstring.(shift_left_one (numbits tmp)) in
let hole_z = Bitstring.logand (Spindeterminant.bitstring t ) tmp
and particle_z = Bitstring.logand (Spindeterminant.bitstring t') tmp
in
let hole = 1 + Bitstring.trailing_zeros hole_z
and particle = 1 + Bitstring.trailing_zeros particle_z
in
(* Phase calculation *)
let low, high =
if particle > hole then hole, particle
else particle, hole
in
let mask =
let h = high-1 in
let l = low in
let mask_up = shift_left_one h |> Bitstring.minus_one
and mask_dn = Bitstring.plus_one @@ Bitstring.lognot (shift_left_one l)
in Bitstring.logand mask_up mask_dn
in
let phase =
Phase.add (Phase.add (Spindeterminant.phase t) (Spindeterminant.phase t'))
(Phase.of_nperm (Bitstring.popcount @@ Bitstring.logand d mask ))
in
(hole, particle, phase)
let single_of_det t t' =
assert Determinant.(beta t = beta t' || alfa t = alfa t');
if Determinant.(beta t = beta t') then
let hole, particle, phase =
single_of_spindet (Determinant.alfa t) (Determinant.alfa t')
in
Single (phase, { hole ; particle ; spin=Spin.Alfa })
else
let hole, particle, phase =
single_of_spindet (Determinant.beta t) (Determinant.beta t')
in
Single (phase, { hole ; particle ; spin=Spin.Beta })
let multiple_of_spindet t t' =
let holes = Spindeterminant.holes_of t t'
and particles = Spindeterminant.particles_of t t'
in
let t'' =
List.fold_left (fun accu h -> Spindeterminant.annihilation h accu) t holes
in
let t'' =
List.fold_left (fun accu p -> Spindeterminant.creation p accu) t'' particles
in
assert (t' = t'' || t' = Spindeterminant.negate_phase t'');
let phase =
if Spindeterminant.phase t' = Spindeterminant.phase t'' then
Phase.Pos
else
Phase.Neg
in
(phase, List.map2 (fun hole particle -> (hole, particle)) holes (List.rev particles) )
let double_of_spindet t t' =
match multiple_of_spindet t t' with
| (phase, (h1,p1)::(h2,p2)::[]) -> (h1, p1, h2, p2, phase)
| _ -> invalid_arg "t and t' are not doubly excited"
let multiple_of_det t t' =
let pa, a =
multiple_of_spindet (Determinant.alfa t) (Determinant.alfa t')
and pb, b =
multiple_of_spindet (Determinant.beta t) (Determinant.beta t')
in
let phase = Phase.add pa pb in
Multiple (phase, List.concat [
List.map (fun (hole, particle) -> { hole ; particle ; spin=Spin.Alfa }) a ;
List.map (fun (hole, particle) -> { hole ; particle ; spin=Spin.Beta }) b ])
let double_of_det t t' =
match multiple_of_det t t' with
| Multiple (phase, [e1 ; e2]) -> Double (phase, e1, e2)
| _ -> assert false
let of_det t t' =
match Determinant.degree t t' with
| 0 -> if Determinant.phase t = Determinant.phase t' then
Identity Phase.Pos
else
Identity Phase.Neg
| 1 -> single_of_det t t'
| 2 -> double_of_det t t'
| _ -> multiple_of_det t t'
let pp_s_exc ppf t =
Format.fprintf ppf "@[T^{%s}_{%d->%d}@]"
(match t.spin with
| Spin.Alfa -> "\\alpha"
| Spin.Beta -> "\\beta " )
t.hole t.particle
let pp_exc ppf t =
let phase, l =
match t with
| Identity p -> p, []
| Single (p,x) -> p, x::[]
| Double (p,x,y) -> p, x::y::[]
| Multiple (p,l) -> p, l
in
Format.fprintf ppf "@[%c"
(match phase with
| Phase.Pos -> '+'
| Phase.Neg -> '-' );
List.iter (fun x -> Format.fprintf ppf "@[T^{%s}_{%d->%d}@]"
(match x.spin with
| Spin.Alfa -> "\\alpha"
| Spin.Beta -> "\\beta " )
x.hole x.particle) l;
Format.fprintf ppf "@]"
let test_case () =
(*
let test_id () =
let l_a = [ 1 ; 2 ; 3 ; 5 ; 64 ]
and l_b = [ 2 ; 3 ; 5 ; 65 ] in
let det1 = Determinant.of_lists l_a l_b in
let det2 = Determinant.negate_phase det1 in
in
*)
let test_single () =
let l_a = [ 1 ; 2 ; 3 ; 5 ; 64 ]
and l_b = [ 2 ; 3 ; 5 ; 65 ] in
let det1 = Determinant.of_lists 66 l_a l_b in
let det2 = Determinant.single_excitation Spin.Alfa 3 7 det1 in
let t = single_of_det det1 det2 in
Alcotest.(check bool) "single 1" true (t = Single (Phase.Pos, { hole=3 ; particle=7 ; spin=Spin.Alfa}) );
let det2 =
Determinant.single_excitation Spin.Alfa 2 7 det1
|> Determinant.negate_phase
in
let t = single_of_det det1 det2 in
Alcotest.(check bool) "single 2" true (t = Single (Phase.Neg, { hole=2 ; particle=7 ; spin=Spin.Alfa }) );
let det2 = Determinant.single_excitation Spin.Beta 2 7 det1 in
let t = single_of_det det1 det2 in
Alcotest.(check bool) "single 3" true (t = Single (Phase.Pos, { hole=2 ; particle=7 ; spin=Spin.Beta}) );
let det2 = Determinant.single_excitation Spin.Beta 3 256 det1 in
let t = single_of_det det1 det2 in
Alcotest.(check bool) "single 4" true (t = Single (Phase.Pos, { hole=3 ; particle=256 ; spin=Spin.Beta}) );
in
let test_double () =
let l_a = [ 1 ; 2 ; 3 ; 5 ; 64 ]
and l_b = [ 2 ; 3 ; 5 ; 65 ] in
let det1 = Determinant.of_lists 66 l_a l_b in
let det2 = Determinant.double_excitation Spin.Alfa 3 7 Spin.Alfa 2 6 det1 in
let t = double_of_det det1 det2 in
Alcotest.(check bool) "double 1" true
(t = Double (Phase.Neg,
{ hole=2 ; particle=7 ; spin=Spin.Alfa},
{ hole=3 ; particle=6 ; spin=Spin.Alfa}));
in
[
"Single", `Quick, test_single;
"Double", `Quick, test_double;
]