Double excitations

CI/Determinant.ml

@@ -37,6 +37,10 @@ let of_spindeterminants a 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 of_lists a b =
   let alfa = Spindeterminant.of_list a 
   and beta  = Spindeterminant.of_list b

CI/Determinant.mli

@@ -56,6 +56,10 @@ val of_spindeterminants : Spindeterminant.t -> Spindeterminant.t -> t
 val of_lists : int list -> int list -> t
 (** Creates a Slater determinant from a two lists of orbital indices.  *)
 
+val negate_phase : t -> t
+(** Returns the same determinant with the phase negated. *)
+
+
 (** {1 Printers} *)
 
 val pp_det : Format.formatter -> t -> unit

CI/Excitation.ml

@@ -3,14 +3,14 @@ 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
+| 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);
@@ -31,97 +31,143 @@ let single_of_spindet t t' =
   in
   let mask =
     let h = high-1 in
+    let l = low in
     let mask_up = Z.(shift_left one h - one)
-    and mask_dn = Z.(neg (shift_left one low) + one)
+    and mask_dn = Z.(lognot (shift_left one l) + 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
+    Phase.add (Phase.add (Spindeterminant.phase t) (Spindeterminant.phase t'))
+      (Phase.of_nperm (Z.(popcount @@ 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
-    { hole ; particle ; phase ; spin=Spin.Alfa }
+    Single (phase, { hole ; particle ; spin=Spin.Alfa })
   else
     let hole, particle, phase =
       single_of_spindet (Determinant.beta t) (Determinant.beta t') 
     in
-    { hole ; particle ; phase ; spin=Spin.Beta }
+    Single (phase, { hole ; particle ; 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'
+let multiple_of_spindet t t' =
+  let holes     = Spindeterminant.holes_of t t'
+  and particles = Spindeterminant.particles_of t 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
+  let t'' = 
+    List.fold_left (fun accu h -> Spindeterminant.annihilation h accu) t holes
   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
+  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 = 
-    Phase.(to_nperm d + to_nperm d') + Z.(popcount @@ logand d mask )
-    |> Phase.of_nperm
+    if Spindeterminant.phase t' = Spindeterminant.phase t'' then
+      Phase.Pos
+    else
+      Phase.Neg
   in
-  (hole, particle, phase)
-*)
+  (phase, List.map2 (fun hole particle -> (hole, particle)) holes (List.rev particles) )
 
-let pp_exc ppf t =
-  Format.fprintf ppf "@[%cT^{%s}_{%d->%d}@]"
-    (match t.phase with
-      | Phase.Pos -> '+'
-      | Phase.Neg -> '-' )
+
+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 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 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 
+    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" (t = { hole=2 ; particle=7 ; spin=Spin.Alfa ; phase=Phase.Neg} ) true;
+    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" (t = { hole=2 ; particle=7 ; spin=Spin.Beta ; phase=Phase.Pos} ) true;
+    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" (t = { hole=3 ; particle=256 ; spin=Spin.Beta ; phase=Phase.Pos} ) true;
+    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 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_single;
-    *)
+    "Double", `Quick, test_double;
   ]
 
 

CI/Phase.ml

@@ -10,15 +10,20 @@ let to_nperm = function
 | Pos -> 0
 | Neg -> 1
 
+let add t t' =
+  match t, t' with
+  | Pos, Pos 
+  | Neg, Neg -> Pos
+  | Pos, Neg 
+  | Neg, Pos -> Neg
+
+let neg = function
+  | Pos -> Neg
+  | Neg -> Pos
+
 let add_nperm phase = function
   | 0 -> phase
-  | nperm ->
-    begin
-      match (phase, of_nperm nperm) with
-      | (Pos,Pos) | (Neg,Neg) -> Pos
-      | _ -> Neg
-    end
-
+  | nperm -> add phase (of_nperm nperm)
 
 let pp_phase ppf = function
   | Pos -> Format.fprintf ppf "@[<h>+1@]"

CI/Phase.mli

@@ -8,9 +8,14 @@ val of_nperm : int -> t
 val to_nperm : t -> int
 (** Converts the phase to [1] or [0] permutations. *)
 
+val add : t -> t -> t
+(** Add a given phase to an existing phase. *)
+
 val add_nperm : t -> int -> t
 (** Add to an existing phase a given number of permutations. *)
 
+val neg : t -> t
+(** Negate the phase. *)
 
 (** {1 Printers} *)
 

CI/Spindeterminant.ml

@@ -91,6 +91,9 @@ let particles_of t t' =
   Z.logand (bitstring t') (Z.logxor (bitstring t) (bitstring t'))
   |> bits_to_list []
 
+let negate_phase = function
+  | Some t -> Some { t with phase = Phase.neg t.phase }
+  | None -> None
 
 
 let of_list l = 
@@ -110,14 +113,10 @@ let rec to_list = function
     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
+      Format.fprintf ppf "@[<h>%a %a@]" Phase.pp_phase s.phase Util.pp_bitstring s.bitstring
 
 
 

CI/Spindeterminant.mli

@@ -23,6 +23,8 @@ val bitstring : t -> Z.t
 val is_none : t -> bool
 (** Tests if a spin-determinant is [None]. *)
 
+val negate_phase : t -> t
+(** Returns a spin-determinant with the phase reversed. *)
 
 (** {1 Second quantization operators} *)
 
@@ -44,6 +46,13 @@ val double_excitation : hole -> particle -> hole -> particle -> t -> t
 val degree : t -> t -> int
 (** Returns degree of excitation between two determinants. *)
 
+val holes_of : t -> t -> int list
+(** Returns the list of holes in the excitation from one determinant to another. *)
+
+val particles_of : t -> t -> int list
+(** Returns the list of particles in the excitation from one determinant to another. *)
+
+
 val of_list : int list -> t
 (** Builds a spin-determinant from a list of orbital indices. If the creation of the
     spin-determinant is not possible because of Pauli's exclusion principle, a [None]

Utils/Util.ml

@@ -268,6 +268,11 @@ let pp_matrix ppf m =
   aux 1 cols
 
 
+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 string_of_matrix m = 

Utils/Util.mli

@@ -137,5 +137,8 @@ val pp_float_2darray : Format.formatter -> float array array -> unit
 ]}
 *)
 
+val pp_bitstring : Format.formatter -> Z.t -> unit
+(** Example: [ +++++------+-- ] *)
+
 val pp_matrix : Format.formatter -> Lacaml.D.mat -> unit