Single excitations

CI/Excitation.ml

@@ -0,0 +1,127 @@
+type single_exc =
+  {
+    hole     : int    ;
+    particle : int    ;
+    spin     : Spin.t ;
+    phase    : Phase.t;
+  }
+
+type t =
+| Identity
+| Single    of single_exc
+| Double    of single_exc * single_exc
+| Multiple  of 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 = Z.logxor d d' in
+  let hole_z     = Z.logand (Spindeterminant.bitstring t ) tmp
+  and particle_z = Z.logand (Spindeterminant.bitstring t') tmp
+  in
+  let hole     = 1 + Z.trailing_zeros hole_z
+  and particle = 1 + Z.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 mask_up = Z.(shift_left one h - one)
+    and mask_dn = Z.(neg (shift_left one low) + one)
+    in Z.logand mask_up mask_dn
+  in
+  let phase = 
+    Phase.to_nperm (Spindeterminant.phase t ) +
+    Phase.to_nperm (Spindeterminant.phase t') + 
+    Z.(popcount @@ logand d mask )
+    |> Phase.of_nperm
+  in
+  (hole, particle, phase)
+
+
+let single_of_det t t' =
+  if Determinant.(beta t = beta t') then
+    let hole, particle, phase =
+      single_of_spindet (Determinant.alfa t) (Determinant.alfa t') 
+    in
+    { hole ; particle ; phase ; spin=Spin.Alfa }
+  else
+    let hole, particle, phase =
+      single_of_spindet (Determinant.beta t) (Determinant.beta t') 
+    in
+    { hole ; particle ; phase ; spin=Spin.Beta }
+
+(*
+
+let double_of_spindet t t' =
+  assert (Spindeterminant.degree t t' = 2);
+  let d  = Spindeterminant.bitstring t
+  and d' = Spindeterminant.bitstring t'
+  in
+  let tmp = Z.logxor d d' in
+  let hole_z     = Z.logand (Spindeterminant.bitstring t ) tmp
+  and particle_z = Z.logand (Spindeterminant.bitstring t') tmp
+  in
+  let hole     = 1 + Z.trailing_zeros hole_z
+  and particle = 1 + Z.trailing_zeros particle_z
+  in
+  (* Phase calculation *)
+  let low, high = 
+    if particle > hole then hole, particle
+    else particle, hole
+  in
+  let mask =
+    let mask_up = Z.(shift_left one (high-1) - one) 
+    and mask_dn = Z.(neg (shift_left one low) + one)
+    in Z.logand mask_up mask_dn
+  in
+  let phase = 
+    Phase.(to_nperm d + to_nperm d') + Z.(popcount @@ logand d mask )
+    |> Phase.of_nperm
+  in
+  (hole, particle, phase)
+*)
+
+let pp_exc ppf t =
+  Format.fprintf ppf "@[%cT^{%s}_{%d->%d}@]"
+    (match t.phase with
+      | Phase.Pos -> '+'
+      | Phase.Neg -> '-' )
+    (match t.spin with
+      | Spin.Alfa -> "\\alpha"
+      | Spin.Beta -> "\\beta " )
+    t.hole t.particle
+
+
+let test_case () =
+
+  let test_single () =
+    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.single_excitation Spin.Alfa 3 7 det1 in 
+    let t = single_of_det det1 det2 in
+    Alcotest.(check bool) "single 1" (t = { hole=3 ; particle=7 ; spin=Spin.Alfa ; phase=Phase.Neg} ) true;
+    let det2 = Determinant.single_excitation Spin.Alfa 2 7 det1 in 
+    let t = single_of_det det1 det2 in
+    Alcotest.(check bool) "single 2" (t = { hole=2 ; particle=7 ; spin=Spin.Alfa ; phase=Phase.Neg} ) true;
+    let det2 = Determinant.single_excitation Spin.Beta 2 7 det1 in 
+    let t = single_of_det det1 det2 in
+    Alcotest.(check bool) "single 3" (t = { hole=2 ; particle=7 ; spin=Spin.Beta ; phase=Phase.Pos} ) true;
+    let det2 = Determinant.single_excitation Spin.Beta 3 256 det1 in 
+    let t = single_of_det det1 det2 in
+    Alcotest.(check bool) "single 4" (t = { hole=3 ; particle=256 ; spin=Spin.Beta ; phase=Phase.Pos} ) true;
+  in
+  [
+    "Single", `Quick, test_single;
+    (*
+    "Double", `Quick, test_single;
+    *)
+  ]
+
+

CI/Phase.ml

@@ -6,6 +6,10 @@ let of_nperm nperm =
   if (nperm land 1) = 1 then Neg
   else Pos
 
+let to_nperm = function
+| Pos -> 0
+| Neg -> 1
+
 let add_nperm phase = function
   | 0 -> phase
   | nperm ->

CI/Phase.mli

@@ -5,6 +5,9 @@ type t =
 val of_nperm : int -> t
 (** Returns the phase obtained by a given number of permuations. *)
 
+val to_nperm : t -> int
+(** Converts the phase to [1] or [0] permutations. *)
+
 val add_nperm : t -> int -> t
 (** Add to an existing phase a given number of permutations. *)
 

CI/Spindeterminant.ml

@@ -73,6 +73,25 @@ let double_excitation_reference h' p' h p spindet =
 let double_excitation h' p' h p =
   double_excitation_reference h' p' h p
 
+let rec bits_to_list accu = function
+  | t when (t = Z.zero) -> List.rev accu
+  | t -> let newlist =
+           (Z.trailing_zeros t + 1)::accu
+         in
+         bits_to_list newlist Z.(logand t (t-one))
+
+let degree t t' =
+  Z.hamdist (bitstring t) (bitstring t') / 2
+
+let holes_of t t' =
+  Z.logand (bitstring t) (Z.logxor (bitstring t) (bitstring t'))
+  |> bits_to_list []
+
+let particles_of t t' =
+  Z.logand (bitstring t') (Z.logxor (bitstring t) (bitstring t'))
+  |> bits_to_list []
+
+
 
 let of_list l = 
   List.rev l
@@ -182,10 +201,26 @@ let test_case () =
     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 l_a in
+    let l_b = [ 1 ; 7 ; 3 ; 5 ] in
+    let det2 = of_list 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 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;
   ]
 
 

CI/Spindeterminant.mli

@@ -41,6 +41,8 @@ val single_excitation : hole -> particle -> t -> t
 val double_excitation : hole -> particle -> hole -> particle -> t -> t
 (** Double excitation operator {% $T_{hh'}^{pp'} = a^\dagger_p a^\dagger_{p'} a_{h'} a_h$ %}. *)
 
+val degree : t -> t -> int
+(** Returns degree of excitation between two determinants. *)
 
 val of_list : int list -> t
 (** Builds a spin-determinant from a list of orbital indices. If the creation of the

run_tests.ml

@@ -17,6 +17,7 @@ let test_water_dz () =
     "Guess", Guess.test_case ao_basis;
     "Spindeterminant", Spindeterminant.test_case ();
     "Determinant", Determinant.test_case ();
+    "Excitation", Excitation.test_case ();
   ]