From 47691d0a0a69e4057e6b7f7fccd1fac00e69587d Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 20 Feb 2019 19:43:16 +0100 Subject: [PATCH] Determinant space OK --- Basis/PrimitiveShellPair.ml | 20 +++++++-------- CI/Determinant.ml | 8 +++--- CI/Determinant.mli | 3 ++- CI/Determinant_space.ml | 48 ++++++++++++++++++++++++++++++++++++ CI/Determinant_space.mli | 32 ++++++++++++++++++++++++ CI/Spindeterminant.ml | 4 +-- CI/Spindeterminant.mli | 3 ++- CI/Spindeterminant_space.ml | 44 ++++++++++++++------------------- CI/Spindeterminant_space.mli | 31 +++++++++++++++++------ MOBasis/MOClass.ml | 16 ++++++++++++ MOBasis/MOClass.mli | 5 ++++ SCF/HartreeFock.ml | 2 +- SCF/RHF.ml | 2 +- Utils/Electrons.ml | 17 +++++++------ Utils/Electrons.mli | 16 +++++++----- Utils/Util.ml | 4 +-- Utils/Util.mli | 4 +-- 17 files changed, 190 insertions(+), 69 deletions(-) create mode 100644 CI/Determinant_space.ml create mode 100644 CI/Determinant_space.mli diff --git a/Basis/PrimitiveShellPair.ml b/Basis/PrimitiveShellPair.ml index 3646e1d..3c07783 100644 --- a/Basis/PrimitiveShellPair.ml +++ b/Basis/PrimitiveShellPair.ml @@ -3,15 +3,15 @@ open Constants type t = { - exponent : float; (* alpha + beta *) - exponent_inv : float; (* 1/(alpha + beta) *) - center : Coordinate.t; (* P = (alpha * A + beta B)/(alpha+beta) *) - center_minus_a : Coordinate.t; (* P - A *) - a_minus_b : Coordinate.t; (* A - B *) - a_minus_b_sq : float; (* |A-B|^2 *) + exponent : float; (* {% $\alpha + \beta$ %} *) + exponent_inv : float; (* {% $1/(\alpha + \beta)$ %} *) + center : Coordinate.t; (* {% $P = (\alpha A + \beta B)/(\alpha+\beta)$ %} *) + center_minus_a : Coordinate.t; (* {% $P - A$ %} *) + a_minus_b : Coordinate.t; (* {% $A - B$ %} *) + a_minus_b_sq : float; (* {% $|A-B|^2$ %} *) norm_scales : float array lazy_t; - normalization : float; (* norm_coef_a * norm_coef_b * g, with - g = (pi/(alpha+beta))^(3/2) exp (-|A-B|^2 * alpha*beta/(alpha+beta)) *) + normalization : float; (* [norm_coef_a * norm_coef_b * g], with + {% $g = (\pi/(\alpha+\beta))^(3/2) \exp (-|A-B|^2 \alpha\beta/(\alpha+\beta))$ %} *) ang_mom : AngularMomentum.t; shell_a : PrimitiveShell.t; shell_b : PrimitiveShell.t; @@ -69,7 +69,7 @@ let create_make_of p_a p_b = Ps.normalization p_a in - let alpha_a = + let alfa_a = Co.( Ps.exponent p_a |. Ps.center p_a ) in @@ -101,7 +101,7 @@ let create_make_of p_a p_b = in let center = - Co.(exponent_inv |. (alpha_a |+ beta_b)) + Co.(exponent_inv |. (alfa_a |+ beta_b)) in let center_minus_a = diff --git a/CI/Determinant.ml b/CI/Determinant.ml index 349df7f..ee5fba3 100644 --- a/CI/Determinant.ml +++ b/CI/Determinant.ml @@ -76,11 +76,11 @@ let double_excitation spin h p spin' h' p' t = -let pp_det ppf t = - Format.fprintf ppf "@[@[phase:%a@]@;@[a:%a@];@[b:%a@]@]@." +let pp_det n ppf t = + Format.fprintf ppf "@[@[phase:%a@]@;@[a:%a@]@;@[b:%a@]@]@." Phase.pp_phase (phase t) - Spindeterminant.pp_spindet t.alfa - Spindeterminant.pp_spindet t.beta + (Spindeterminant.pp_spindet n) t.alfa + (Spindeterminant.pp_spindet n) t.beta diff --git a/CI/Determinant.mli b/CI/Determinant.mli index 6bbf6b8..84b4a87 100644 --- a/CI/Determinant.mli +++ b/CI/Determinant.mli @@ -62,7 +62,8 @@ val negate_phase : t -> t (** {1 Printers} *) -val pp_det : Format.formatter -> t -> unit +val pp_det : int -> Format.formatter -> t -> unit +(** First [int] is the number of MOs to print. *) (** {1 Unit tests} *) diff --git a/CI/Determinant_space.ml b/CI/Determinant_space.ml new file mode 100644 index 0000000..6b2dbe7 --- /dev/null +++ b/CI/Determinant_space.ml @@ -0,0 +1,48 @@ +type t = +{ + n_alfa : int ; + n_beta : int ; + mo_class : MOClass.t ; + mo_basis : MOBasis.t ; + determinants : Determinant.t array; +} + +module Ss = Spindeterminant_space + +let n_alfa t = t.n_alfa +let n_beta t = t.n_beta +let mo_class t = t.mo_class +let mo_basis t = t.mo_basis +let determinants t = t.determinants + +let fci_of_mo_basis ?(frozen_core=true) mo_basis = + let s = MOBasis.simulation mo_basis in + let e = Simulation.electrons s in + let n_alfa = Electrons.n_alfa e + and n_beta = Electrons.n_beta e in + let det_a = + Ss.fci_of_mo_basis ~frozen_core mo_basis n_alfa + and det_b = + Ss.fci_of_mo_basis ~frozen_core mo_basis n_beta + in + let mo_class = Ss.mo_class det_a in + let determinants = + let a = Ss.spin_determinants det_a + and b = Ss.spin_determinants det_b + in + let sb = Ss.size det_b in + Array.to_list a + |> List.map (fun a_i -> + Array.init sb (fun j -> + Determinant.of_spindeterminants a_i b.(j)) ) + |> Array.concat + in + { n_alfa ; n_beta ; mo_class ; mo_basis ; determinants } + + +let pp_det_space ppf t = + Format.fprintf ppf "@[[ "; + Array.iteri (fun i d -> Format.fprintf ppf "@[@[%8d@]@;@[%a@]@]@;" i + (Determinant.pp_det (MOBasis.size (mo_basis t))) d) t.determinants ; + Format.fprintf ppf "]@]" + diff --git a/CI/Determinant_space.mli b/CI/Determinant_space.mli new file mode 100644 index 0000000..7fd8d22 --- /dev/null +++ b/CI/Determinant_space.mli @@ -0,0 +1,32 @@ +(** +The determinant space in which we solve the Schrodinger equation. +*) + +type t + +(** {1 Accessors} *) + +val n_alfa : t -> int +(** Number of {% $\alpha$ %} electrons in the {% $\alpha$ %} MOs. *) + +val n_beta : t -> int +(** Number of {% $\beta$ %} electrons in the {% $\beta$ %} MOs. *) + +val mo_class : t -> MOClass.t +(** The MO classes used to generate the space. *) + +val mo_basis : t -> MOBasis.t +(** The MO basis on which the determinants are expanded. *) + +val determinants : t -> Determinant.t array +(** All the determinants belonging to the space. *) + +val fci_of_mo_basis : ?frozen_core:bool -> MOBasis.t -> t +(** Creates a space of all possible ways to put [n_alfa] electrons in the {% $\alpha$ %} + [Active] MOs and [n_beta] electrons in the {% $\beta$ %} [Active] MOs. + All other MOs are untouched. +*) + +(** {2 Printing} *) + +val pp_det_space : Format.formatter -> t -> unit diff --git a/CI/Spindeterminant.ml b/CI/Spindeterminant.ml index dc75a9c..f56a9f4 100644 --- a/CI/Spindeterminant.ml +++ b/CI/Spindeterminant.ml @@ -116,10 +116,10 @@ let rec to_list = function in aux [] spindet.bitstring -let pp_spindet ppf = function +let pp_spindet n ppf = function | None -> Format.fprintf ppf "@[None@]" | Some s -> - Format.fprintf ppf "@[%a %a@]" Phase.pp_phase s.phase Util.pp_bitstring s.bitstring + Format.fprintf ppf "@[%a %a@]" Phase.pp_phase s.phase (Util.pp_bitstring n) s.bitstring diff --git a/CI/Spindeterminant.mli b/CI/Spindeterminant.mli index a11bf95..b90d848 100644 --- a/CI/Spindeterminant.mli +++ b/CI/Spindeterminant.mli @@ -69,7 +69,8 @@ val to_list : t -> int list (** {1 Printers}. *) -val pp_spindet : Format.formatter -> t -> unit +val pp_spindet : int -> Format.formatter -> t -> unit +(** First [int] is the number of MOs to print *) (** {1 Unit testing} *) diff --git a/CI/Spindeterminant_space.ml b/CI/Spindeterminant_space.ml index 38ca23e..701631f 100644 --- a/CI/Spindeterminant_space.ml +++ b/CI/Spindeterminant_space.ml @@ -1,48 +1,42 @@ type t = -{ - elec_num : int; - mo_basis : MOBasis.t; - mo_class : MOClass.t; - spindets : Spindeterminant.t array; -} + { + elec_num : int; + mo_basis : MOBasis.t; + mo_class : MOClass.t; + spin_determinants : Spindeterminant.t array; + } +let spin_determinants t = t.spin_determinants +let elec_num t = t.elec_num +let mo_basis t = t.mo_basis +let mo_class t = t.mo_class +let size t = Array.length t.spin_determinants + let fci_of_mo_basis ?(frozen_core=true) mo_basis elec_num = let mo_num = MOBasis.size mo_basis in - let ncore = (Nuclei.small_core @@ Simulation.nuclei @@ MOBasis.simulation mo_basis) / 2 in - let mo_class = - MOClass.of_list @@ - if frozen_core then - List.concat [ - Util.list_range 1 ncore - |> List.map (fun i -> MOClass.Core i) ; - Util.list_range (ncore+1) mo_num - |> List.map (fun i -> MOClass.Active i) - ] - else - Util.list_range 1 mo_num - |> List.map (fun i -> MOClass.Active i) - in + let mo_class = MOClass.fci ~frozen_core mo_basis in let m l = List.fold_left (fun accu i -> let j = i-1 in Z.(logor accu (shift_left one j)) - ) Z.zero l + ) Z.zero l in let occ_mask = m (MOClass.core_mos mo_class) and active_mask = m (MOClass.active_mos mo_class) in let neg_active_mask = Z.lognot active_mask in - let spindets = + let spin_determinants = Util.bit_permtutations elec_num mo_num |> List.filter (fun b -> Z.logand neg_active_mask b = occ_mask) |> List.map (fun b -> Spindeterminant.of_bitstring b) |> Array.of_list in - { elec_num ; mo_basis ; mo_class ; spindets } + { elec_num ; mo_basis ; mo_class ; spin_determinants } 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 "@[[ "; + Array.iteri (fun i d -> Format.fprintf ppf "@[@[%8d@] @[%a@]@]@;" i + (Spindeterminant.pp_spindet (MOBasis.size (mo_basis t))) d) (spin_determinants t) ; Format.fprintf ppf "]@]" diff --git a/CI/Spindeterminant_space.mli b/CI/Spindeterminant_space.mli index e743f41..1f6cc08 100644 --- a/CI/Spindeterminant_space.mli +++ b/CI/Spindeterminant_space.mli @@ -1,11 +1,28 @@ -type t = -{ - elec_num : int ; (** Number of electrons *) - mo_basis : MOBasis.t ; (** MO basis on which the space is built *) - mo_class : MOClass.t ; (** CI Classes of the MOs *) - spindets : Spindeterminant.t array ; (** Spin-determinants belonging to the space *) -} +(** +The space built with determinants made with same-spin spinorbitals. +*) +type t + +(** {1 Accessors} *) + +val size : t -> int +(** Number of determinants in the space. *) + +val spin_determinants : t -> Spindeterminant.t array +(** All the spin-determinants belonging to the space. *) + +val elec_num : t -> int +(** Number of (same-spin) electrons occupying the MOs. *) + +val mo_class : t -> MOClass.t +(** The MO classes used to generate the space. *) + +val mo_basis : t -> MOBasis.t +(** The MO basis on which the determinants are expanded. *) + + +(** {1 Creation} *) val fci_of_mo_basis : ?frozen_core:bool -> MOBasis.t -> int -> t (** Create a space of all possible ways to put [n_elec] electrons in the [Active] MOs. diff --git a/MOBasis/MOClass.ml b/MOBasis/MOClass.ml index d39594c..7f23a98 100644 --- a/MOBasis/MOClass.ml +++ b/MOBasis/MOClass.ml @@ -55,3 +55,19 @@ let deleted_mos t = |> Util.list_some +let fci ?(frozen_core=true) mo_basis = + let mo_num = MOBasis.size mo_basis in + let ncore = (Nuclei.small_core @@ Simulation.nuclei @@ MOBasis.simulation mo_basis) / 2 in + of_list ( + if frozen_core then + List.concat [ + Util.list_range 1 ncore + |> List.map (fun i -> Core i) ; + Util.list_range (ncore+1) mo_num + |> List.map (fun i -> Active i) + ] + else + Util.list_range 1 mo_num + |> List.map (fun i -> Active i) + ) + diff --git a/MOBasis/MOClass.mli b/MOBasis/MOClass.mli index 495ea47..b43ef33 100644 --- a/MOBasis/MOClass.mli +++ b/MOBasis/MOClass.mli @@ -12,6 +12,11 @@ type t (** Creation *) val of_list : mo_class list -> t +val fci : ?frozen_core:bool -> MOBasis.t -> t +(** Creates the MO classes for FCI calculations : all [Active]. The + [n] lowest MOs are [Core] if [frozen_core = true]. +*) + val core_mos : t -> int list (** Returns a list containing the indices of the core MOs. *) diff --git a/SCF/HartreeFock.ml b/SCF/HartreeFock.ml index 396efc3..5b8305f 100644 --- a/SCF/HartreeFock.ml +++ b/SCF/HartreeFock.ml @@ -1,7 +1,7 @@ open Util let make ?guess simulation = - if (Simulation.electrons simulation).Electrons.multiplicity = 1 then + if Electrons.multiplicity @@ Simulation.electrons simulation = 1 then RHF.make ?guess simulation else invalid_arg "UHF or ROHF not implemented" diff --git a/SCF/RHF.ml b/SCF/RHF.ml index b029a8f..62b705e 100644 --- a/SCF/RHF.ml +++ b/SCF/RHF.ml @@ -12,7 +12,7 @@ let make ?guess:(guess=`Huckel) ?max_scf:(max_scf=64) ?level_shift:(level_shift= (* Number of occupied MOs *) let nocc = - (Si.electrons simulation).El.n_alpha + El.n_alfa @@ Si.electrons simulation in let nuclear_repulsion = diff --git a/Utils/Electrons.ml b/Utils/Electrons.ml index b8b41f8..642ae93 100644 --- a/Utils/Electrons.ml +++ b/Utils/Electrons.ml @@ -1,8 +1,8 @@ -(** Number of alpha and beta electrons *) +(** Number of {% $\alpha$ %} and {% $\beta$ %} electrons *) type t = { - n_alpha : int ; - n_beta : int ; + n_alfa : int ; + n_beta : int ; multiplicity : int; } @@ -15,13 +15,16 @@ let make ?multiplicity:(multiplicity=1) ?charge:(charge=0) nuclei = let negative_charges = charge - positive_charges in let n_elec = - negative_charges in let n_beta = ((n_elec - multiplicity)+1)/2 in - let n_alpha = n_elec - n_beta in - let result = { n_alpha ; n_beta ; multiplicity } in - if multiplicity <> (n_alpha - n_beta)+1 then + let n_alfa = n_elec - n_beta in + let result = { n_alfa ; n_beta ; multiplicity } in + if multiplicity <> (n_alfa - n_beta)+1 then invalid_arg (__FILE__^": make"); result let charge e = - - (e.n_alpha + e.n_beta) + - (e.n_alfa + e.n_beta) |> Charge.of_int +let n_alfa t = t.n_alfa +let n_beta t = t.n_beta +let multiplicity t = t.multiplicity diff --git a/Utils/Electrons.mli b/Utils/Electrons.mli index 68c6cca..caf0122 100644 --- a/Utils/Electrons.mli +++ b/Utils/Electrons.mli @@ -1,11 +1,6 @@ (** Information related to electrons. *) -type t = { - n_alpha : int ; (** Number of alpha electrons *) - n_beta : int ; (** Number of beta electrons *) - multiplicity : int ; (** Spin multiplicity: {% $2S+1$ %} *) -} - +type t val make : ?multiplicity:int -> ?charge:int -> Nuclei.t -> t (** Creates the data relative to electrons for a molecular system @@ -19,3 +14,12 @@ val make : ?multiplicity:int -> ?charge:int -> Nuclei.t -> t val charge : t -> Charge.t (** Sum of the charges of the electrons. *) +val n_alfa : t -> int +(** Number of alpha electrons *) + +val n_beta : t -> int +(** Number of beta electrons *) + +val multiplicity : t -> int +(** Spin multiplicity: {% $2S+1$ %} *) + diff --git a/Utils/Util.ml b/Utils/Util.ml index fc69f93..ff3ec7c 100644 --- a/Utils/Util.ml +++ b/Utils/Util.ml @@ -298,8 +298,8 @@ 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 '-') +let pp_bitstring n ppf bs = + String.init n (fun i -> if (Z.testbit bs i) then '+' else '-') |> Format.fprintf ppf "@[%s@]" diff --git a/Utils/Util.mli b/Utils/Util.mli index f667cf7..01a6d4e 100644 --- a/Utils/Util.mli +++ b/Utils/Util.mli @@ -156,8 +156,8 @@ val pp_float_2darray : Format.formatter -> float array array -> unit ]} *) -val pp_bitstring : Format.formatter -> Z.t -> unit -(** Example: [ +++++------+-- ] *) +val pp_bitstring : int -> Format.formatter -> Z.t -> unit +(** Example: [ pp_bitstring 14 ppf z -> +++++------+-- ] *) val pp_matrix : Format.formatter -> Lacaml.D.mat -> unit