mirror of https://gitlab.com/scemama/QCaml.git
250 lines
7.5 KiB
OCaml
250 lines
7.5 KiB
OCaml
type s =
|
|
{
|
|
bitstring : Bitstring.t;
|
|
phase : Phase.t ;
|
|
}
|
|
|
|
type t = s option
|
|
type hole = int
|
|
type particle = int
|
|
|
|
let phase = function
|
|
| Some s -> s.phase
|
|
| None -> Phase.Pos
|
|
|
|
|
|
let is_none = function
|
|
| None -> true
|
|
| _ -> false
|
|
|
|
|
|
let bitstring = function
|
|
| Some s -> s.bitstring
|
|
| None -> invalid_arg "Spindeterminant is None"
|
|
|
|
|
|
let vac n =
|
|
Some { bitstring = Bitstring.zero n;
|
|
phase = Phase.Pos; }
|
|
|
|
|
|
let creation p = function
|
|
| None -> None
|
|
| Some spindet ->
|
|
let i = pred p in
|
|
if Bitstring.testbit spindet.bitstring i then
|
|
None
|
|
else
|
|
begin
|
|
let numbits = Bitstring.numbits spindet.bitstring in
|
|
let x = Bitstring.shift_left_one numbits i in
|
|
let bitstring = Bitstring.logor spindet.bitstring x in
|
|
let mask = Bitstring.minus_one x in
|
|
let r = Bitstring.logand bitstring mask in
|
|
let phase = Phase.add_nperm spindet.phase (Bitstring.popcount r) in
|
|
Some { bitstring ; phase }
|
|
end
|
|
|
|
|
|
let annihilation h = function
|
|
| None -> None
|
|
| Some spindet ->
|
|
let i = pred h in
|
|
if not (Bitstring.testbit spindet.bitstring i) then
|
|
None
|
|
else
|
|
begin
|
|
let numbits = Bitstring.numbits spindet.bitstring in
|
|
let x = Bitstring.shift_left_one numbits i in
|
|
let mask = Bitstring.minus_one x in
|
|
let r = Bitstring.logand spindet.bitstring mask in
|
|
let phase = Phase.add_nperm spindet.phase (Bitstring.popcount r) in
|
|
let bitstring = Bitstring.logand spindet.bitstring (Bitstring.lognot x) in
|
|
Some { bitstring ; phase }
|
|
end
|
|
|
|
let single_excitation_reference h p spindet =
|
|
creation p @@ annihilation h @@ spindet
|
|
|
|
let single_excitation h p =
|
|
single_excitation_reference h p
|
|
|
|
|
|
let double_excitation_reference h' p' h p spindet =
|
|
creation p' @@ creation p @@ annihilation h @@ annihilation h' @@ spindet
|
|
|
|
let double_excitation h' p' h p =
|
|
double_excitation_reference h' p' h p
|
|
|
|
|
|
let degree t t' =
|
|
Bitstring.hamdist (bitstring t) (bitstring t') / 2
|
|
|
|
let holes_of t t' =
|
|
Bitstring.logand (bitstring t) (Bitstring.logxor (bitstring t) (bitstring t'))
|
|
|> Bitstring.to_list
|
|
|
|
let particles_of t t' =
|
|
Bitstring.logand (bitstring t') (Bitstring.logxor (bitstring t) (bitstring t'))
|
|
|> Bitstring.to_list
|
|
|
|
let holes_particles_of t t' =
|
|
let x = Bitstring.logxor (bitstring t) (bitstring t') in
|
|
let holes = Bitstring.logand (bitstring t) x |> Bitstring.to_list
|
|
and particles = Bitstring.logand (bitstring t') x |> Bitstring.to_list
|
|
in
|
|
List.rev_map2 (fun h p -> (h,p)) holes particles
|
|
|> List.rev
|
|
|
|
|
|
let set_phase p = function
|
|
| Some t -> Some { t with phase = p }
|
|
| None -> None
|
|
|
|
let negate_phase = function
|
|
| Some t -> Some { t with phase = Phase.neg t.phase }
|
|
| None -> None
|
|
|
|
|
|
let of_bitstring ?(phase=Phase.Pos) bitstring =
|
|
Some { bitstring ; phase }
|
|
|
|
let of_list n l =
|
|
List.rev l
|
|
|> List.fold_left (fun accu p -> creation p accu) (vac n)
|
|
|
|
|
|
let to_list = function
|
|
| None -> []
|
|
| Some spindet ->
|
|
let rec aux accu z =
|
|
if not (Bitstring.is_zero z) then
|
|
let element = ((Bitstring.trailing_zeros z)+1) in
|
|
(aux [@tailcall]) (element::accu) (Bitstring.logand z (Bitstring.minus_one z) )
|
|
else List.rev accu
|
|
in aux [] spindet.bitstring
|
|
|
|
let to_array t =
|
|
to_list t
|
|
|> Array.of_list
|
|
|
|
let n_electrons = function
|
|
| Some t -> Bitstring.popcount t.bitstring
|
|
| None -> 0
|
|
|
|
|
|
let pp n ppf = function
|
|
| None -> Format.fprintf ppf "@[<h>None@]"
|
|
| Some s ->
|
|
Format.fprintf ppf "@[<h>%a %a@]" Phase.pp s.phase Bitstring.pp
|
|
s.bitstring
|
|
|
|
|
|
|
|
|
|
(*-----------------------------------------------------------------------------------*)
|
|
|
|
|
|
let test_case () =
|
|
|
|
let test_creation () =
|
|
let l_a = [ 1 ; 2 ; 3 ; 5 ] in
|
|
let det = of_list 10 l_a in
|
|
Alcotest.(check (list int )) "bitstring 1" l_a (to_list det);
|
|
Alcotest.(check bool) "phase 2" true (phase det = Phase.Pos);
|
|
let l_b = [ 1 ; 3 ; 2 ; 5 ] in
|
|
let det = of_list 10 l_b in
|
|
Alcotest.(check (list int )) "bitstring 2" l_a (to_list det);
|
|
Alcotest.(check bool) "phase 2" true (phase det = Phase.Neg);
|
|
in
|
|
|
|
let test_a_operators () =
|
|
let det =
|
|
creation 5 @@ creation 2 @@ creation 2 @@ creation 1 @@ (vac 10)
|
|
in
|
|
Alcotest.(check bool) "none 1" true (is_none det);
|
|
|
|
let det =
|
|
creation 5 @@ creation 3 @@ creation 2 @@ creation 1 @@ (vac 10)
|
|
in
|
|
let l_a = [ 1 ; 2 ; 3 ; 5 ] in
|
|
Alcotest.(check (list int )) "bitstring 1" l_a (to_list det);
|
|
Alcotest.(check bool) "phase 1" true (phase det = Phase.Pos);
|
|
|
|
let det =
|
|
creation 1 @@ creation 3 @@ creation 2 @@ creation 5 @@ (vac 10)
|
|
in
|
|
Alcotest.(check (list int )) "bitstring 2" l_a (to_list det);
|
|
Alcotest.(check bool) "phase 2" true (phase det = Phase.Neg);
|
|
|
|
let l_b = [ 1 ; 3 ; 2 ; 5 ] in
|
|
let det = of_list 10 l_b in
|
|
Alcotest.(check (list int )) "bitstring 3" l_a (to_list det);
|
|
Alcotest.(check bool) "phase 3" true (phase det = Phase.Neg);
|
|
|
|
Alcotest.(check bool) "none 1" true (annihilation 4 det |> is_none);
|
|
|
|
let det =
|
|
annihilation 1 det
|
|
in
|
|
Alcotest.(check (list int )) "bitstring 4" (List.tl l_a) (to_list det);
|
|
Alcotest.(check bool) "phase 4" true (phase det = Phase.Neg);
|
|
|
|
let det =
|
|
annihilation 3 det
|
|
in
|
|
Alcotest.(check (list int )) "bitstring 5" [ 2 ; 5 ] (to_list det);
|
|
Alcotest.(check bool) "phase 5" true (phase det = Phase.Pos);
|
|
|
|
let det =
|
|
annihilation 5 @@ annihilation 2 det
|
|
in
|
|
Alcotest.(check (list int )) "bitstring 6" [] (to_list det);
|
|
Alcotest.(check bool) "phase 6" true (phase det = Phase.Pos);
|
|
|
|
in
|
|
|
|
let test_exc_operators () =
|
|
let l_a = [ 1 ; 2 ; 3 ; 5 ] in
|
|
let det = of_list 10 l_a in
|
|
let l_b = [ 1 ; 7 ; 3 ; 5 ] in
|
|
let det2 = of_list 10 l_b in
|
|
Format.printf "%a@." (pp 7) det;
|
|
Format.printf "%a@." (pp 7) det2;
|
|
Format.printf "%a@." (pp 7) (single_excitation_reference 2 7 det);
|
|
Alcotest.(check bool) "single 1" true (single_excitation_reference 2 7 det = det2);
|
|
Alcotest.(check bool) "single 2" true (single_excitation 2 7 det = single_excitation_reference 2 7 det);
|
|
Alcotest.(check bool) "single 3" true (single_excitation_reference 4 7 det |> is_none);
|
|
Alcotest.(check bool) "single 4" true (single_excitation 4 7 det |> is_none);
|
|
|
|
let l_c = [ 1 ; 7 ; 6 ; 5 ] in
|
|
let det3 = of_list 10 l_c in
|
|
Alcotest.(check bool) "double 1" true (double_excitation_reference 2 7 3 6 det = det3);
|
|
Alcotest.(check bool) "double 2" true (double_excitation 2 7 3 6 det = double_excitation_reference 2 7 3 6 det);
|
|
Alcotest.(check bool) "double 3" true (double_excitation_reference 4 7 3 6 det |> is_none);
|
|
Alcotest.(check bool) "double 4" true (double_excitation 4 7 3 6 det |> is_none);
|
|
in
|
|
|
|
let test_exc_spindet () =
|
|
let l_a = [ 1 ; 2 ; 3 ; 5 ] in
|
|
let det = of_list 10 l_a in
|
|
let l_b = [ 1 ; 7 ; 3 ; 5 ] in
|
|
let det2 = of_list 10 l_b in
|
|
Alcotest.(check int) "single" 1 (degree det det2);
|
|
Alcotest.(check (list int)) "holes" [2] (holes_of det det2);
|
|
Alcotest.(check (list int)) "particles" [7] (particles_of det det2);
|
|
let l_b = [ 1 ; 7 ; 3 ; 6 ] in
|
|
let det2 = of_list 10 l_b in
|
|
Alcotest.(check int) "double" 2 (degree det det2);
|
|
Alcotest.(check (list int)) "holes" [2 ; 5] (holes_of det det2);
|
|
Alcotest.(check (list int)) "particles" [6 ; 7] (particles_of det det2);
|
|
in
|
|
[
|
|
"Creation", `Quick, test_creation;
|
|
"Creation/Annihilation Operators", `Quick, test_a_operators;
|
|
"Excitation Operators", `Quick, test_exc_operators;
|
|
"Excitation of spindet", `Quick, test_exc_spindet;
|
|
]
|
|
|
|
|