Single/Double excitation operators

CI/Determinant.ml

@@ -1,17 +1,27 @@
 type t =
 {
-  alpha : Spindeterminant.t ;
+  alfa : Spindeterminant.t ;
-  beta  : Spindeterminant.t ;
+  beta : Spindeterminant.t ;
 }
 
+type hole = int
+type particle = int
 
 
-let alpha t = t.alpha
+
+
+let alfa t = t.alfa
 
 let beta  t = t.beta 
 
+let vac =
+  {
+    alfa = Spindeterminant.vac ;
+    beta = Spindeterminant.vac ;
+  }
+
 let phase t = 
-  match Spindeterminant.(phase t.alpha, phase t.beta) with
+  match Spindeterminant.(phase t.alfa, phase t.beta) with
   | Phase.Pos, Phase.Pos 
   | Phase.Neg, Phase.Neg  -> Phase.Pos
   | _ -> Phase.Neg
@@ -19,32 +29,66 @@ let phase t =
 
 let of_spindeterminants a b =
   {
-    alpha = a ;
+    alfa = a ;
     beta  = b
   }
 
+
+let is_none t = Spindeterminant.(is_none t.alfa || is_none t.beta)
+
+
 let of_lists a b =
-  {
+  let alfa = Spindeterminant.of_list a 
-    alpha = Spindeterminant.of_list a ;
+  and beta  = Spindeterminant.of_list b
-    beta  = Spindeterminant.of_list 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 ppf t =
-  Format.fprintf ppf "@[<v> a: %a @; b: %a @]@."
+  Format.fprintf ppf "@[<v>@[phase:%a@]@;@[a:%a@];@[b:%a@]@]@."
-    Spindeterminant.pp_spindet t.alpha
+    Phase.pp_phase (phase t)
+    Spindeterminant.pp_spindet t.alfa
     Spindeterminant.pp_spindet 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 l_a l_b in
-      let z_a = alpha det
+      let z_a = alfa det
       and z_b = beta  det in
-      Alcotest.(check (list int )) "alpha" (Spindeterminant.to_list z_a) l_a;
+      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
@@ -67,9 +111,78 @@ let test_case () =
       let det = of_lists 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
+      in
+      Alcotest.(check bool) "creation 1" true
+        (det = of_lists [ 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 [ 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 [ 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 [ 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 [ 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 [ 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 [ 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;
   ]
 
 

CI/Determinant.mli

@@ -7,10 +7,12 @@ The {% $\alpha$ %} and {% $\beta$ %} determinants are of type [Spindeterminant.t
 
 
 type t
+type hole = int
+type particle = int
 
 (** {1 Accessors} *)
 
-val alpha : t -> Spindeterminant.t
+val alfa : t -> Spindeterminant.t
 (** Get the {% $\alpha$ %} spin-determinant. *)
 
 val beta : t -> Spindeterminant.t
@@ -21,6 +23,29 @@ val phase : t -> Phase.t
     spin-determinants.
 *)
 
+val is_none : t -> bool
+(** Tests if a Determinant is [None]. *)
+
+
+(** {1 Second quantization operators} *)
+
+val vac : t
+(** Vacuum state, [vac = Some ]{% $|\rangle$ %}  *)
+
+val creation : Spin.t -> particle -> t -> t
+(** [creation spin p] is the creation operator {% $a^\dagger_p$ %}. *)
+
+val annihilation : Spin.t -> hole -> t -> t
+(** [annihilation spin h] is the annihilation operator {% $a_h$ %}. *)
+
+val single_excitation : Spin.t -> hole -> particle -> t -> t
+(** Single excitation operator {% $T_h^p = a^\dagger_p a_h$ %}. *)
+
+val double_excitation : Spin.t -> hole -> particle -> Spin.t -> hole -> particle -> t -> t
+(** Double excitation operator {% $T_{hh'}^{pp'} = a^\dagger_p a^\dagger_{p'} a_{h'} a_h$ %}. *)
+
+
+
 (** {1 Creators} *)
 
 val of_spindeterminants : Spindeterminant.t -> Spindeterminant.t -> t

CI/Spindeterminant.ml

@@ -10,7 +10,7 @@ type particle = int
 
 let phase = function 
   | Some s -> s.phase
-  | None -> invalid_arg "Spindeterminant is None"
+  | None -> Phase.Pos
 
 
 let is_none = function

CI/Spindeterminant.mli

@@ -38,6 +38,9 @@ val annihilation : hole -> t -> t
 val single_excitation : hole -> particle -> t -> t
 (** Single excitation operator {% $T_h^p = a^\dagger_p a_h$ %}. *)
 
+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 of_list : int list -> t
 (** Builds a spin-determinant from a list of orbital indices. If the creation of the

Utils/Spin.ml

@@ -0,0 +1,6 @@
+(** Electron spin *)
+
+type t = Alfa | Beta
+
+
+

Utils/Spin.mli

@@ -0,0 +1,12 @@
+(** Electron spin.
+
+Note :
+[Alfa] if written with an 'f' instead of 'ph' because it has the same number of
+letters as [Beta], so the alignment of the code is nicer.
+
+*)
+
+type t = Alfa | Beta
+
+
+