QCaml/CI/Determinant.ml

type t =
{
  alfa : Spindeterminant.t ;
  beta : Spindeterminant.t ;
}

type hole = int
type particle = int




let alfa t = t.alfa

let beta  t = t.beta 

let vac n =
  {
    alfa = Spindeterminant.vac n;
    beta = Spindeterminant.vac n;
  }

let phase t = 
  match Spindeterminant.(phase t.alfa, phase t.beta) with
  | Phase.Pos, Phase.Pos 
  | Phase.Neg, Phase.Neg  -> Phase.Pos
  | _ -> Phase.Neg


let of_spindeterminants a b =
  {
    alfa = a ;
    beta  = b
  }


let is_none t = Spindeterminant.(is_none t.alfa || is_none t.beta)


let negate_phase t =
  { t with alfa = Spindeterminant.negate_phase t.alfa }


let degree_alfa t t' =
  Spindeterminant.degree t.alfa t'.alfa

let degree_beta t t' =
  Spindeterminant.degree t.beta t'.beta

let degrees t t' =
  degree_alfa t t', degree_beta t t' 


let degree t t' =
  (degree_alfa t t') + (degree_beta t t')


let of_lists n a b =
  let alfa = Spindeterminant.of_list n a 
  and beta  = Spindeterminant.of_list n b
  in of_spindeterminants alfa beta


let creation spin p t =
  match spin with
  | Spin.Alfa -> { t with alfa = Spindeterminant.creation p t.alfa }
  | Spin.Beta -> { t with beta = Spindeterminant.creation p t.beta }


let annihilation spin h t =
  match spin with
  | Spin.Alfa -> { t with alfa = Spindeterminant.annihilation h t.alfa }
  | Spin.Beta -> { t with beta = Spindeterminant.annihilation h t.beta }


let single_excitation spin h p t =
  match spin with
  | Spin.Alfa -> { t with alfa = Spindeterminant.single_excitation h p t.alfa }
  | Spin.Beta -> { t with beta = Spindeterminant.single_excitation h p t.beta }


let double_excitation spin h p spin' h' p' t = 
  match spin, spin' with
  | Spin.(Alfa, Beta)  ->        { alfa = Spindeterminant.single_excitation h  p      t.alfa ;
                                   beta = Spindeterminant.single_excitation h' p'     t.beta }
  | Spin.(Beta, Alfa)  ->        { beta = Spindeterminant.single_excitation h  p      t.beta ;
                                   alfa = Spindeterminant.single_excitation h' p'     t.alfa }
  | Spin.(Alfa, Alfa)  -> { t with alfa = Spindeterminant.double_excitation h p h' p' t.alfa }
  | Spin.(Beta, Beta)  -> { t with beta = Spindeterminant.double_excitation h p h' p' t.beta }



let pp_det n ppf t =
  Format.fprintf ppf "@[<v>@[phase:%a@]@;@[a:%a@]@;@[b:%a@]@]@."
    Phase.pp_phase (phase t)
    (Spindeterminant.pp_spindet n) t.alfa
    (Spindeterminant.pp_spindet n) t.beta



let test_case () =

  let test_creation () =
      let l_a = [ 1 ; 2 ; 3 ; 5 ; 64 ]
      and l_b = [ 2 ; 3 ; 5 ; 65 ] in
      let det = of_lists 66 l_a l_b in
      let z_a = alfa det
      and z_b = beta  det in
      Alcotest.(check (list int )) "alfa" (Spindeterminant.to_list z_a) l_a;
      Alcotest.(check (list int )) "beta"  (Spindeterminant.to_list z_b) l_b;
      Alcotest.(check bool) "phase" (phase det = Phase.Pos)  true;
  in

  let test_phase () =
      let l_a = [ 1 ; 2 ; 3 ; 64 ; 5 ]
      and l_b = [ 2 ; 3 ; 5 ; 65 ] in
      let det = of_lists 66 l_a l_b in
      Alcotest.(check bool) "phase" (phase det = Phase.Neg)  true;
      let l_a = [ 1 ; 2 ; 3 ; 64 ; 5 ]
      and l_b = [ 3 ; 2 ; 5 ; 65 ] in
      let det = of_lists 66 l_a l_b in
      Alcotest.(check bool) "phase" (phase det = Phase.Pos)  true;
      let l_a = [ 1 ; 3 ; 2 ; 64 ; 5 ]
      and l_b = [ 3 ; 2 ; 5 ; 65 ] in
      let det = of_lists 66 l_a l_b in
      Alcotest.(check bool) "phase" (phase det = Phase.Neg)  true;
      let l_a = [ 1 ; 3 ; 2 ; 64 ; 5 ]
      and l_b = [ 3 ; 2 ; 65 ; 5 ] in
      let det = of_lists 66 l_a l_b in
      Alcotest.(check bool) "phase" (phase det = Phase.Pos)  true;
  in

  let test_operators () =
      let det =
        let open Spin in
        creation Alfa 1 @@ creation Alfa 3 @@ creation Alfa 2 @@ creation Alfa 5 @@ 
        creation Beta 1 @@ creation Beta 3 @@ creation Beta 4 @@ creation Beta 5 @@ vac 10
      in
      Alcotest.(check bool) "creation 1" true
        (det = of_lists 10 [ 1 ; 3 ; 2 ; 5 ] [1 ; 3 ; 4 ; 5 ] );

      let det' =
        single_excitation Spin.Alfa 3 6 det
      in
      Alcotest.(check bool) "single_exc 1" true
        (det' = of_lists 10 [ 1 ; 6 ; 2 ; 5 ] [1 ; 3 ; 4 ; 5 ] );

      let det' =
        single_excitation Spin.Beta  3 6 det
      in
      Alcotest.(check bool) "single_exc 2" true
        (det' = of_lists 10 [ 1 ; 3 ; 2 ; 5 ] [1 ; 6 ; 4 ; 5 ] );

      let det' =
        single_excitation Spin.Alfa 4 6 det
      in
      Alcotest.(check bool) "single_exc 3" true (is_none det');

      let det' =
        single_excitation Spin.Beta  1 5 det
      in
      Alcotest.(check bool) "single_exc 4" true (is_none det');

      let det' =
        double_excitation Spin.Alfa 3 6 Spin.Alfa 2 7 det
      in
      let det'' = of_lists 10 [ 1 ; 6 ; 7 ; 5 ] [1 ; 3 ; 4 ; 5 ] in
      Alcotest.(check bool) "double_exc 1" true (det' = det'');

      let det' =
        double_excitation Spin.Beta  3 6 Spin.Beta 5 7  det
      in
      Alcotest.(check bool) "double_exc 2" true
        (det' = of_lists 10 [ 1 ; 3 ; 2 ; 5 ] [1 ; 6 ; 4 ; 7 ] );

      let det' =
        double_excitation Spin.Alfa 3 6 Spin.Beta  5 7 det
      in
      Alcotest.(check bool) "double_exc 3" true
        (det' = of_lists 10 [ 1 ; 6 ; 2 ; 5 ] [1 ; 3 ; 4 ; 7 ] );

      let det' =
        double_excitation Spin.Beta 5 7 Spin.Alfa 3 6 det
      in
      Alcotest.(check bool) "double_exc 4" true
        (det' = of_lists 10 [ 1 ; 6 ; 2 ; 5 ] [1 ; 3 ; 4 ; 7 ] );

      let det' =
        double_excitation Spin.Alfa 4 6 Spin.Alfa 2 7 det
      in
      Alcotest.(check bool) "double_exc 5" true (is_none det');

      let det' =
        double_excitation Spin.Beta  1 5 Spin.Alfa 2 7  det
      in
      Alcotest.(check bool) "double_exc 6" true (is_none det');

  in

  [
    "Creation", `Quick, test_creation;
    "Phase",    `Quick, test_phase;
    "Operators",`Quick, test_operators;
  ]