From b858e8a4b8ac4d3a417fe8fe4cee90870ddea145 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 19 Feb 2019 17:36:07 +0100 Subject: [PATCH] Added Anderson's algorithm --- CI/Spindeterminant.ml | 3 +++ CI/Spindeterminant.mli | 6 ++++++ CI/Spindeterminant_space.ml | 26 ++++++++++++++++++++++++++ MOBasis/MOBasis.ml | 4 ++++ MOBasis/MOBasis.mli | 4 +++- Parallel_mpi/Farm.ml | 5 +++-- Utils/Util.ml | 26 ++++++++++++++++++++++++++ Utils/Util.mli | 19 +++++++++++++++++++ 8 files changed, 90 insertions(+), 3 deletions(-) create mode 100644 CI/Spindeterminant_space.ml diff --git a/CI/Spindeterminant.ml b/CI/Spindeterminant.ml index cd13a88..dc75a9c 100644 --- a/CI/Spindeterminant.ml +++ b/CI/Spindeterminant.ml @@ -96,6 +96,9 @@ let negate_phase = function | None -> None +let of_bitstring ?(phase=Phase.Pos) bitstring = + Some { bitstring ; phase } + let of_list l = List.rev l |> List.fold_left (fun accu p -> creation p accu) vac diff --git a/CI/Spindeterminant.mli b/CI/Spindeterminant.mli index 17ee494..a11bf95 100644 --- a/CI/Spindeterminant.mli +++ b/CI/Spindeterminant.mli @@ -26,6 +26,7 @@ val is_none : t -> bool val negate_phase : t -> t (** Returns a spin-determinant with the phase reversed. *) + (** {1 Second quantization operators} *) val vac : t @@ -53,6 +54,11 @@ val particles_of : t -> t -> int list (** Returns the list of particles in the excitation from one determinant to another. *) +(** {1 Creation} *) + +val of_bitstring : ?phase:Phase.t -> Z.t -> t +(** Creates from a bitstring and an optional phase.*) + 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/CI/Spindeterminant_space.ml b/CI/Spindeterminant_space.ml new file mode 100644 index 0000000..5d6d9a6 --- /dev/null +++ b/CI/Spindeterminant_space.ml @@ -0,0 +1,26 @@ +type t = +{ + elec_num : int; + mo_basis : MOBasis.t; + spindets : Spindeterminant.t array; +} + +let fci_of_mo_basis mo_basis elec_num = + let mo_num = MOBasis.size mo_basis in + let spindets = + Util.bit_permtutations elec_num mo_num + |> List.map (fun b -> Spindeterminant.of_bitstring b) + |> Array.of_list + in + { elec_num ; mo_basis ; spindets } + + +let pp_spindet_space ppf t = + Format.fprintf ppf "@[["; + Array.iteri (fun i d -> Format.fprintf ppf "@[@[%d@]@;@[%a@]@]@," i Spindeterminant.pp_spindet d) t.spindets; + Format.fprintf ppf "]@]" + + + + + diff --git a/MOBasis/MOBasis.ml b/MOBasis/MOBasis.ml index 8f4d47e..026ceb3 100644 --- a/MOBasis/MOBasis.ml +++ b/MOBasis/MOBasis.ml @@ -32,6 +32,9 @@ type t = } +let size t = + Mat.dim2 t.mo_coef + let mo_matrix_of_ao_matrix ~mo_coef ao_matrix = xt_o_x ~x:mo_coef ~o:ao_matrix @@ -204,3 +207,4 @@ let of_hartree_fock ~frozen_core = function | HF.RHF hf -> of_rhf ~frozen_core hf | _ -> assert false + diff --git a/MOBasis/MOBasis.mli b/MOBasis/MOBasis.mli index 07e5f76..aa85b69 100644 --- a/MOBasis/MOBasis.mli +++ b/MOBasis/MOBasis.mli @@ -28,7 +28,6 @@ type t = private } - val make : ao_basis:AOBasis.t -> mo_type:mo_type -> mo_class:mo_class array -> @@ -41,3 +40,6 @@ val of_hartree_fock : frozen_core:bool -> HartreeFock_type.t -> t (** Build MOs from a Restricted Hartree-Fock calculation. *) +val size : t -> int +(** Number of molecular orbitals in the basis *) + diff --git a/Parallel_mpi/Farm.ml b/Parallel_mpi/Farm.ml index 2b619f3..fb89ba6 100644 --- a/Parallel_mpi/Farm.ml +++ b/Parallel_mpi/Farm.ml @@ -5,8 +5,9 @@ let run_sequential f stream = try let task = Stream.next stream in Some (f task) - with Stream.Failure -> None in - Stream.from next + with Stream.Failure -> None + in + Stream.from next (* Multi-process functions *) diff --git a/Utils/Util.ml b/Utils/Util.ml index fda3cf0..5d61f5c 100644 --- a/Utils/Util.ml +++ b/Utils/Util.ml @@ -96,6 +96,17 @@ let fact = function | i -> fact_memo.(i) +let binom n k = + (*TODO : slow function *) + assert (n > k); + let rec aux n k = + if k = 0 || k = n then + 1 + else + aux (n-1) (k-1) + aux (n-1) k + in aux n k + + let rec pow a = function | 0 -> 1. | 1 -> a @@ -225,6 +236,21 @@ let canonical_ortho ?thresh:(thresh=1.e-6) ~overlap c = +(** {2 Bitstring functions} *) +let bit_permtutations m n = + + let rec aux k u rest = + if k=1 then + List.rev (u :: rest) + else + let t = Z.(logor u (u-one)) in + let t' = Z.(t+one) in + let t'' = Z.(shift_right ((logand (lognot t) t') - one)) (Z.trailing_zeros u + 1) in + aux (k-1) (Z.logor t' t'') (u :: rest) + in + aux (binom n m) Z.(shift_left one m - one) [] + + (** {2 Printers} *) let pp_float_array_size ppf a = diff --git a/Utils/Util.mli b/Utils/Util.mli index f97e666..60762fb 100644 --- a/Utils/Util.mli +++ b/Utils/Util.mli @@ -23,6 +23,9 @@ val fact : int -> float @raise Invalid_argument for negative arguments or arguments >100. *) +val binom : int -> int -> int +(** Binomial coefficient. [binom n k] {% $= C_n^k$ %}. *) + val pow : float -> int -> float (** Fast implementation of the power function for small integer powers *) @@ -103,6 +106,22 @@ val sym_matrix_of_file : string -> Lacaml.D.mat (** Reads a symmetric matrix from a file with format "%d %d %f" corresponding to [i, j, A.{i,j}]. *) + +(** {2 Bitstring functions} *) + +val bit_permtutations : int -> int -> Z.t list +(** [bit_permtutations m n] generates the list of all possible [n]-bit + strings with [m] bits set to 1. + Algorithm adapted from + {{:https://graphics.stanford.edu/~seander/bithacks.html#NextBitPermutation} + Bit twiddling hacks}. + Example: +{[ + bit_permtutations 2 4 = [ 0011 ; 0101 ; 0110 ; 1001 ; 1010 ; 1100 ] +]} +*) + + (** {2 Printers} *) val pp_float_array_size : Format.formatter -> float array -> unit (** Example: