From cc000f6bba15a2dbf63be7d43f5f741df0e4caed Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 18 Feb 2019 19:45:41 +0100 Subject: [PATCH] Double excitations --- CI/Determinant.ml | 4 ++ CI/Determinant.mli | 4 ++ CI/Excitation.ml | 144 +++++++++++++++++++++++++++-------------- CI/Phase.ml | 19 ++++-- CI/Phase.mli | 5 ++ CI/Spindeterminant.ml | 9 ++- CI/Spindeterminant.mli | 9 +++ Utils/Util.ml | 5 ++ Utils/Util.mli | 3 + 9 files changed, 141 insertions(+), 61 deletions(-) diff --git a/CI/Determinant.ml b/CI/Determinant.ml index 61ab53e..349df7f 100644 --- a/CI/Determinant.ml +++ b/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 diff --git a/CI/Determinant.mli b/CI/Determinant.mli index f63271b..6bbf6b8 100644 --- a/CI/Determinant.mli +++ 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 diff --git a/CI/Excitation.ml b/CI/Excitation.ml index fb34e84..ad14f53 100644 --- a/CI/Excitation.ml +++ b/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; ] diff --git a/CI/Phase.ml b/CI/Phase.ml index 93cdf4d..cfa9e07 100644 --- a/CI/Phase.ml +++ b/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 "@[+1@]" diff --git a/CI/Phase.mli b/CI/Phase.mli index dc1b3e2..f186370 100644 --- a/CI/Phase.mli +++ b/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} *) diff --git a/CI/Spindeterminant.ml b/CI/Spindeterminant.ml index d0f3504..cd13a88 100644 --- a/CI/Spindeterminant.ml +++ b/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 "@[%s@]" - let pp_spindet ppf = function | None -> Format.fprintf ppf "@[None@]" | Some s -> - Format.fprintf ppf "@[%a %a@]" Phase.pp_phase s.phase pp_bitstring s.bitstring + Format.fprintf ppf "@[%a %a@]" Phase.pp_phase s.phase Util.pp_bitstring s.bitstring diff --git a/CI/Spindeterminant.mli b/CI/Spindeterminant.mli index ded02f5..17ee494 100644 --- a/CI/Spindeterminant.mli +++ b/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] diff --git a/Utils/Util.ml b/Utils/Util.ml index 25c11ec..fda3cf0 100644 --- a/Utils/Util.ml +++ b/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 "@[%s@]" + + let string_of_matrix m = diff --git a/Utils/Util.mli b/Utils/Util.mli index 61c9fac..f97e666 100644 --- a/Utils/Util.mli +++ b/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