From f94f54e4d9ae4bd24491ad71b63611bb3135d55c Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 18 Feb 2019 15:39:26 +0100 Subject: [PATCH] Single excitations --- CI/Excitation.ml | 127 +++++++++++++++++++++++++++++++++++++++++ CI/Phase.ml | 4 ++ CI/Phase.mli | 3 + CI/Spindeterminant.ml | 35 ++++++++++++ CI/Spindeterminant.mli | 2 + run_tests.ml | 1 + 6 files changed, 172 insertions(+) create mode 100644 CI/Excitation.ml diff --git a/CI/Excitation.ml b/CI/Excitation.ml new file mode 100644 index 0000000..fb34e84 --- /dev/null +++ b/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; + *) + ] + + diff --git a/CI/Phase.ml b/CI/Phase.ml index 4861c57..93cdf4d 100644 --- a/CI/Phase.ml +++ b/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 -> diff --git a/CI/Phase.mli b/CI/Phase.mli index 4c04ea9..dc1b3e2 100644 --- a/CI/Phase.mli +++ b/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. *) diff --git a/CI/Spindeterminant.ml b/CI/Spindeterminant.ml index be7efe7..d0f3504 100644 --- a/CI/Spindeterminant.ml +++ b/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; ] diff --git a/CI/Spindeterminant.mli b/CI/Spindeterminant.mli index dcebabc..ded02f5 100644 --- a/CI/Spindeterminant.mli +++ b/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 diff --git a/run_tests.ml b/run_tests.ml index 1b884ab..7b8ad30 100644 --- a/run_tests.ml +++ b/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 (); ]