diff --git a/.merlin b/.merlin index c5ed8ad..9f51dcc 100644 --- a/.merlin +++ b/.merlin @@ -1,4 +1,4 @@ -PKG str unix bigarray lacaml +PKG str unix bigarray lacaml zarith mpi alcotest S . S Test S Nuclei diff --git a/CI/Determinant.ml b/CI/Determinant.ml index c0fa584..6d7c342 100644 --- a/CI/Determinant.ml +++ b/CI/Determinant.ml @@ -1,14 +1,18 @@ -type sign_type = +type phase_type = | Plus | Minus -type spin_determinant_type = Z.t +type spin_determinant_type = +{ + bitstring : Z.t ; + phase : phase_type; +} type determinant_type = { - alpha : spin_determinant_type ; - beta : spin_determinant_type ; - sign : sign_type; + alpha : Z.t ; + beta : Z.t ; + phase : phase_type; } @@ -20,11 +24,11 @@ let spindet_of_list l = | [] -> accu, if nperm mod 2 = 0 then Plus else Minus | i :: rest -> let i = pred i in - let x = Z.(one lsl i) in - let accu = Z.(x lor accu) in + let x = Z.(shift_left one i) in + let accu = Z.logor accu x in let nperm = let mask = Z.(x-one) in - let r = Z.(accu land mask) in + let r = Z.logand accu mask in if r = Z.zero then nperm else @@ -48,25 +52,25 @@ let rec to_list spindet = let of_lists a b = - let alpha, sign_a = + let alpha, phase_a = spindet_of_list a in - let beta, sign_b = + let beta, phase_b = spindet_of_list b in - let sign = - match sign_a, sign_b with + let phase = + match phase_a, phase_b with | Plus , Plus -> Plus | Minus, Minus-> Plus | _ -> Minus in - { alpha ; beta ; sign } + { alpha ; beta ; phase } let alpha det = det.alpha let beta det = det.beta let sgn det = - match det.sign with + match det.phase with | Plus -> 1. | Minus -> -1. diff --git a/CI/Phase.ml b/CI/Phase.ml new file mode 100644 index 0000000..d4febb2 --- /dev/null +++ b/CI/Phase.ml @@ -0,0 +1,12 @@ +type t = +| Pos +| Neg + +let of_nperm nperm = + if (nperm land 1) = 1 then Neg + else Pos + +let pp_phase ppf = function + | Pos -> Format.fprintf ppf "@[+1@]" + | Neg -> Format.fprintf ppf "@[-1@]" + diff --git a/CI/Phase.mli b/CI/Phase.mli new file mode 100644 index 0000000..5aaeec7 --- /dev/null +++ b/CI/Phase.mli @@ -0,0 +1,10 @@ +type t = +| Pos +| Neg + +val of_nperm : int -> t +(** Returns the phase obtained by a given number of permuations. *) + + +(** Formatters *) +val pp_phase : Format.formatter -> t -> unit diff --git a/CI/Spindeterminant.ml b/CI/Spindeterminant.ml new file mode 100644 index 0000000..02b6c66 --- /dev/null +++ b/CI/Spindeterminant.ml @@ -0,0 +1,64 @@ +type t = + { + bitstring : Z.t ; + phase : Phase.t ; + } + +let of_list l = + let rec aux accu nperm = function + | [] -> { bitstring = accu; + phase = Phase.of_nperm nperm; } + | i :: rest -> + let i = pred i in + let x = Z.(shift_left one i) in + let accu = Z.logor accu x in + let nperm = + let mask = Z.(x-one) in + let r = Z.logand accu mask in + if r = Z.zero then + nperm + else + nperm + (Z.popcount r) + in + aux accu nperm rest + in + List.rev l + |> aux Z.zero 0 + + +(* TODO Phase *) +let rec to_list spindet = + let rec aux accu z = + if z <> Z.zero then + let element = (Z.(trailing_zeros z)+1) in + aux (element::accu) Z.(z land (pred z)) + else List.rev accu + in aux [] spindet + + +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 s = + Format.fprintf ppf "@[%a %a@]" Phase.pp_phase s.phase pp_bitstring s.bitstring + + + +let test_case () = + + let test_creation () = + let l_a = [ 1 ; 2 ; 3 ; 5 ] in + let det = of_list l_a in + Alcotest.(check (list int )) "bitstring" (to_list det.bitstring) l_a; + Alcotest.(check bool) "phase" (det.phase = Phase.Pos) true; + let l_b = [ 1 ; 3 ; 2 ; 5 ] in + let det = of_list l_b in + Alcotest.(check (list int )) "bitstring" (to_list det.bitstring) l_a; + Alcotest.(check bool) "phase" (det.phase = Phase.Neg) true; + in + [ + "Creation", `Quick, test_creation; + ] + + diff --git a/run_tests.ml b/run_tests.ml index 6e66624..1b884ab 100644 --- a/run_tests.ml +++ b/run_tests.ml @@ -15,6 +15,7 @@ let test_water_dz () = Alcotest.run "Water, cc-pVDZ" [ "AO_Basis", AOBasis.test_case ao_basis; "Guess", Guess.test_case ao_basis; + "Spindeterminant", Spindeterminant.test_case (); "Determinant", Determinant.test_case (); ]