QCaml/CI/Spindeterminant.ml

type s =
  {
    bitstring : Z.t ;
    phase     : Phase.t ;
  }

type t = s option
type hole = int
type particle = int

let phase = function 
  | Some s -> s.phase
  | None -> invalid_arg "Spindeterminant is None"


let is_none = function
  | None -> true
  | _ -> false


let bitstring = function 
  | Some s -> s.bitstring
  | None -> invalid_arg "Spindeterminant is None"


let vac = 
  Some  { bitstring = Z.zero;
          phase     = Phase.Pos; }

let creation p = function
  | None -> None
  | Some spindet ->
    let i = pred p in
    if Z.testbit spindet.bitstring i then
      None
    else
      begin
        let x = Z.(shift_left one i) in
        let bitstring = Z.logor spindet.bitstring x in
        let mask = Z.(x-one) in
        let r = Z.logand bitstring mask in
        let phase = Phase.add_nperm spindet.phase (Z.popcount r) in
        Some { bitstring ; phase }
      end


let annihilation h = function
  | None -> None
  | Some spindet ->
    let i = pred h in
    if not (Z.testbit spindet.bitstring i) then
      None
    else
      begin
        let x = Z.(shift_left one i) in
        let mask = Z.(x-one) in
        let r = Z.logand spindet.bitstring mask in
        let phase = Phase.add_nperm spindet.phase (Z.popcount r) in
        let bitstring = Z.logand spindet.bitstring (Z.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 of_list l = 
  List.rev l
  |> List.fold_left (fun accu p -> creation p accu) vac



let rec to_list = function
  | None -> []
  | Some spindet ->
    let rec aux accu z =
      if z <> Z.zero then
        let element = (Z.(trailing_zeros z)+1) in
        aux (element::accu) Z.(z land (pred z))
      else List.rev accu
    in aux [] spindet.bitstring


let pp_bitstring ppf bs =
  String.init (Z.numbits bs) (fun i -> if (Z.testbit bs i) then '+' else '-')
  |> Format.fprintf ppf "@[<h>%s@]"

let pp_spindet ppf = function
  | None -> Format.fprintf ppf "@[<h>None@]"
  | Some s -> 
      Format.fprintf ppf "@[<h>%a %a@]" Phase.pp_phase s.phase pp_bitstring s.bitstring




(*-----------------------------------------------------------------------------------*)


let test_case () =

  let test_creation () =
    let l_a = [ 1 ; 2 ; 3 ; 5 ] in
    let det = of_list 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 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
    in
    Alcotest.(check bool) "none 1" true (is_none det);

    let det = 
      creation 5 @@ creation 3 @@ creation 2 @@ creation 1 @@ vac
    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
    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 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 l_a in
    let l_b = [ 1 ; 7 ; 3 ; 5 ] in
    let det2 = of_list l_b in
    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 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
  [
    "Creation", `Quick, test_creation;
    "Creation/Annihilation Operators", `Quick, test_a_operators;
    "Excitation Operators", `Quick, test_exc_operators;
  ]