From 4c5bade8fb52421188ebd326f2733828c6487d29 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 20 Mar 2019 19:18:36 +0100 Subject: [PATCH] CAS+EN2 in the presence of auxiliary basis OK --- Basis/Basis.ml | 7 +++ Basis/Basis.mli | 5 ++ Basis/GeneralBasis.ml | 12 ++++ Basis/GeneralBasis.mli | 7 +++ CI/CI.ml | 9 ++- CI/DeterminantSpace.ml | 118 +++++++++++++++++++++++++++------------- CI/DeterminantSpace.mli | 11 +++- MOBasis/MOClass.ml | 49 ++++++++++++----- MOBasis/MOClass.mli | 11 ++-- Simulation.ml | 4 +- Simulation.mli | 4 +- run_fci_f12.ml | 89 ++++++++++++++++++++++++++++++ 12 files changed, 262 insertions(+), 64 deletions(-) create mode 100644 run_fci_f12.ml diff --git a/Basis/Basis.ml b/Basis/Basis.ml index 6da5954..dec6c92 100644 --- a/Basis/Basis.ml +++ b/Basis/Basis.ml @@ -91,3 +91,10 @@ let of_nuclei_and_basis_filename ~nuclei filename = of_nuclei_and_general_basis nuclei general_basis +let of_nuclei_and_basis_filenames ~nuclei filenames = + let general_basis = + GeneralBasis.read_many filenames + in + of_nuclei_and_general_basis nuclei general_basis + + diff --git a/Basis/Basis.mli b/Basis/Basis.mli index ac2b055..b5e43d2 100644 --- a/Basis/Basis.mli +++ b/Basis/Basis.mli @@ -18,6 +18,11 @@ val of_nuclei_and_basis_filename : nuclei:Nuclei.t -> string -> t from a file. *) +val of_nuclei_and_basis_filenames : nuclei:Nuclei.t -> string list -> t +(** Same as {!of_nuclei_and_general_basis}, but taking the {!GeneralBasis.t} + from multiple files. + *) + val size : t -> int (** Number of contracted basis functions. *) diff --git a/Basis/GeneralBasis.ml b/Basis/GeneralBasis.ml index 73abc9e..bf2c5d5 100644 --- a/Basis/GeneralBasis.ml +++ b/Basis/GeneralBasis.ml @@ -78,6 +78,18 @@ let read filename = loop () +let read_many filenames = + let h = Hashtbl.create 63 in + let l = + List.map read filenames + |> List.concat + in + List.iter (fun (k,v) -> + let l = Hashtbl.find_all h k in + Hashtbl.replace h k (Array.concat (v::l) ) + ) l; + Hashtbl.fold (fun k v accu -> (k, v)::accu) h [] + let string_of_primitive ?id prim = match id with diff --git a/Basis/GeneralBasis.mli b/Basis/GeneralBasis.mli index 40aad8f..45421c6 100644 --- a/Basis/GeneralBasis.mli +++ b/Basis/GeneralBasis.mli @@ -51,6 +51,13 @@ val read : string -> t *) +val read_many : string list -> t +(** Reads multiple basis set files and return an association list where + the key is an {!Element.t} and the value is the parsed basis set. + This function is usually used to add an auxiliary basis set. +*) + + val read_element : in_channel -> element_basis option (** Reads an element from the input channel [ic]. For example, {[ diff --git a/CI/CI.ml b/CI/CI.ml index e9b0a7c..3a97d7c 100644 --- a/CI/CI.ml +++ b/CI/CI.ml @@ -394,7 +394,7 @@ let make ?(n_states=1) det_space = Matrix.mm ~transa:`T m_H psi in let eigenvectors, eigenvalues = - Davidson.make ~threshold:1.e-12 ~n_states diagonal matrix_prod + Davidson.make ~threshold:1.e-6 ~n_states diagonal matrix_prod in let eigenvalues = Vec.map (fun x -> x +. e_shift) eigenvalues in eigenvectors, eigenvalues @@ -415,9 +415,12 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states } let cls = MOClass.mo_class_array mo_class in Util.list_range 1 (MOBasis.size mo_basis) |> List.filter (fun i -> match cls.(i) with - | MOClass.Deleted _ + | MOClass.Inactive _ + | MOClass.Active _ + | MOClass.Virtual _ -> true + | MOClass.Auxiliary _ + | MOClass.Deleted _ | MOClass.Core _ -> false - | _ -> true ) |> Array.of_list in diff --git a/CI/DeterminantSpace.ml b/CI/DeterminantSpace.ml index a0818a8..e9a195b 100644 --- a/CI/DeterminantSpace.ml +++ b/CI/DeterminantSpace.ml @@ -226,15 +226,15 @@ let fock_diag det_space det = alfa, beta -let cas_of_mo_basis mo_basis n m = + +let spin_of_mo_basis mo_basis f = + 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.cas_of_mo_basis mo_basis n_alfa n m - and det_b = - Ss.cas_of_mo_basis mo_basis n_beta n m + let det_a = f n_alfa + and det_b = f n_beta in let mo_class = Ss.mo_class det_a in let determinants = @@ -248,8 +248,31 @@ let cas_of_mo_basis mo_basis n m = Format.printf "Number of determinants : %d %d %d\n%!" n_det_alfa n_det_beta ndet; Spin (det_alfa, det_beta) + in + { n_alfa ; n_beta ; mo_class ; mo_basis ; determinants } + + +let arbitrary_of_mo_basis mo_basis f = + + 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 = f n_alfa + and det_b = f n_beta + in + let mo_class = Ss.mo_class det_a in + let determinants = + let det_alfa = Ss.spin_determinants det_a + and det_beta = Ss.spin_determinants det_b + in + let n_det_beta = Array.length det_beta in + let n_det_alfa = Array.length det_alfa in + + let ndet = n_det_alfa * n_det_beta in + Format.printf "Number of determinants : %d %d %d\n%!" + n_det_alfa n_det_beta ndet; - (* let det = Array.make n_det_alfa (Array.init n_det_beta (fun i -> i)) in @@ -261,49 +284,70 @@ let cas_of_mo_basis mo_basis n m = Arbitrary { det_alfa ; det_beta ; det ; index_start } - *) in { n_alfa ; n_beta ; mo_class ; mo_basis ; determinants } + +let cas_of_mo_basis mo_basis n m = + let f n_alfa = + Ss.cas_of_mo_basis mo_basis n_alfa n m + in + spin_of_mo_basis mo_basis f + + let fci_of_mo_basis ?(frozen_core=true) mo_basis = + let f n_alfa = + Ss.fci_of_mo_basis ~frozen_core mo_basis n_alfa + in + spin_of_mo_basis mo_basis f + + +let fci_f12_of_mo_basis ?(frozen_core=true) mo_basis mo_num = 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 + let n_core = + if frozen_core then + (Nuclei.small_core @@ Simulation.nuclei @@ MOBasis.simulation mo_basis) / 2 + else 0 in - let mo_class = Ss.mo_class det_a in - let determinants = - let det_alfa = Ss.spin_determinants det_a - and det_beta = Ss.spin_determinants det_b - in - let n_det_beta = Array.length det_beta in - let n_det_alfa = Array.length det_alfa in + let n, m = + (n_alfa + n_beta - n_core), + (mo_num - n_core) + in + let f n_alfa = + Ss.cas_of_mo_basis mo_basis n_alfa n m + in + let r = + spin_of_mo_basis mo_basis f + in + { r with mo_class = + MOClass.to_list r.mo_class + |> List.map (fun i -> + match i with + | MOClass.Virtual i when i > mo_num -> MOClass.Auxiliary i + | i -> i) + |> MOClass.of_list } + - let ndet = n_det_alfa * n_det_beta in - Format.printf "Number of determinants : %d %d %d\n%!" - n_det_alfa n_det_beta ndet; - Spin (det_alfa, det_beta) - - (* - let det = Array.make n_det_alfa - (Array.init n_det_beta (fun i -> i)) - in - let index_start = Array.init (n_det_alfa+1) (fun i -> i*n_det_beta) in - let ndet = (index_start.(n_det_alfa)) in - - Format.printf "Number of determinants : %d %d %d\n%!" - n_det_alfa n_det_beta ndet; - Arbitrary { - det_alfa ; det_beta ; det ; index_start +let cas_f12_of_mo_basis mo_basis n m mo_num = + let f n_alfa = + Ss.cas_of_mo_basis mo_basis n_alfa n m + in + let r = + spin_of_mo_basis mo_basis f + in + { r with mo_class = + MOClass.to_list r.mo_class + |> List.map (fun i -> + match i with + | MOClass.Virtual i when i > mo_num -> MOClass.Auxiliary i + | i -> i) + |> MOClass.of_list } - *) - in - { n_alfa ; n_beta ; mo_class ; mo_basis ; determinants } + let pp_det_space ppf t = diff --git a/CI/DeterminantSpace.mli b/CI/DeterminantSpace.mli index f3cedc4..e8aa4ee 100644 --- a/CI/DeterminantSpace.mli +++ b/CI/DeterminantSpace.mli @@ -57,7 +57,16 @@ val fci_of_mo_basis : ?frozen_core:bool -> MOBasis.t -> t *) val cas_of_mo_basis : MOBasis.t -> int -> int -> t -(** Creates a CAS(n,m) space of determinants. +(** Creates a CAS(n,m) space of determinants. *) + +val fci_f12_of_mo_basis : ?frozen_core:bool -> MOBasis.t -> int -> t +(** Creates the active space to perform a FCI-F12 with an + auxiliary basis set. *) + +val cas_f12_of_mo_basis : MOBasis.t -> int -> int -> int -> t +(** [cas_of_mo_basis mo_basis m n mo_num] Creates a CAS(n,m) space + of determinants with an auxiliary basis set defined as the MOs from + [mo_num+1] to [MOBasis.size mo_basis]. *) (** {2 Printing} *) diff --git a/MOBasis/MOClass.ml b/MOBasis/MOClass.ml index 682708a..b22c17b 100644 --- a/MOBasis/MOClass.ml +++ b/MOBasis/MOClass.ml @@ -1,19 +1,21 @@ type mo_class = - | Core of int (* Always doubly occupied *) - | Inactive of int (* With 0,1 or 2 holes *) - | Active of int (* With 0,1 or 2 holes or particles *) - | Virtual of int (* With 0,1 or 2 particles *) - | Deleted of int (* Always unoccupied *) + | Core of int (* Always doubly occupied *) + | Inactive of int (* With 0,1 or 2 holes *) + | Active of int (* With 0,1 or 2 holes or particles *) + | Virtual of int (* With 0,1 or 2 particles *) + | Deleted of int (* Always unoccupied *) + | Auxiliary of int (* Auxiliary basis function *) type t = mo_class list let pp_mo_occ ppf = function - | Core i -> Format.fprintf ppf "@[Core %d@]" i - | Inactive i -> Format.fprintf ppf "@[Inactive %d@]" i - | Active i -> Format.fprintf ppf "@[Active %d@]" i - | Virtual i -> Format.fprintf ppf "@[Virtual %d@]" i - | Deleted i -> Format.fprintf ppf "@[Deleted %d@]" i + | Core i -> Format.fprintf ppf "@[Core %d@]" i + | Inactive i -> Format.fprintf ppf "@[Inactive %d@]" i + | Active i -> Format.fprintf ppf "@[Active %d@]" i + | Virtual i -> Format.fprintf ppf "@[Virtual %d@]" i + | Deleted i -> Format.fprintf ppf "@[Deleted %d@]" i + | Auxiliary i -> Format.fprintf ppf "@[Auxiliary %d@]" i @@ -62,6 +64,14 @@ let deleted_mos t = |> Util.list_some +let virtual_mos t = + List.map (fun x -> + match x with + | Auxiliary i -> Some i + | _ -> None ) t + |> Util.list_some + + let mo_class_array t = let sze = List.length t + 1 in let result = Array.make sze (Deleted 0) in @@ -72,6 +82,7 @@ let mo_class_array t = | Active i -> result.(i) <- Active i | Virtual i -> result.(i) <- Virtual i | Deleted i -> result.(i) <- Deleted i + | Auxiliary i -> result.(i) <- Auxiliary i ) t; result @@ -96,15 +107,22 @@ let fci ?(frozen_core=true) mo_basis = let cas_sd mo_basis n m = let mo_num = MOBasis.size mo_basis in let n_alfa = MOBasis.simulation mo_basis |> Simulation.electrons |> Electrons.n_alfa in - let n_beta = MOBasis.simulation mo_basis |> Simulation.electrons |> Electrons.n_alfa in + let n_beta = MOBasis.simulation mo_basis |> Simulation.electrons |> Electrons.n_beta in let n_unpaired = n_alfa - n_beta in - let last_inactive = n_alfa - (n + n_unpaired)/2 in + let n_alfa_in_cas = (n - n_unpaired)/2 in + let last_inactive = n_alfa - n_alfa_in_cas in let last_active = last_inactive + m in - let ncore = (Nuclei.small_core @@ Simulation.nuclei @@ MOBasis.simulation mo_basis) / 2 in + let ncore = + (Nuclei.small_core @@ Simulation.nuclei @@ MOBasis.simulation mo_basis) / 2 + |> min last_inactive + in of_list ( List.concat [ - Util.list_range 1 ncore - |> List.map (fun i -> Core i) ; + if ncore > 0 then + Util.list_range 1 ncore + |> List.map (fun i -> Core i) + else + [] ; Util.list_range (ncore+1) last_inactive |> List.map (fun i -> Inactive i) ; Util.list_range (last_inactive+1) last_active @@ -114,3 +132,4 @@ let cas_sd mo_basis n m = ] ) + diff --git a/MOBasis/MOClass.mli b/MOBasis/MOClass.mli index c1c43cf..6abd6f2 100644 --- a/MOBasis/MOClass.mli +++ b/MOBasis/MOClass.mli @@ -1,11 +1,12 @@ (** CI Classes of MOs : active, inactive, etc *) type mo_class = - | Core of int (* Always doubly occupied *) - | Inactive of int (* With 0,1 or 2 holes *) - | Active of int (* With 0,1 or 2 holes or particles *) - | Virtual of int (* With 0,1 or 2 particles *) - | Deleted of int (* Always unoccupied *) + | Core of int (* Always doubly occupied *) + | Inactive of int (* With 0,1 or 2 holes *) + | Active of int (* With 0,1 or 2 holes or particles *) + | Virtual of int (* With 0,1 or 2 particles *) + | Deleted of int (* Always unoccupied *) + | Auxiliary of int (* Function of the auxiliary basis set *) type t diff --git a/Simulation.ml b/Simulation.ml index 5265284..6a3f91f 100644 --- a/Simulation.ml +++ b/Simulation.ml @@ -47,12 +47,12 @@ let make ?cartesian:(cartesian=false) } -let of_filenames ?cartesian:(cartesian=false) ?multiplicity:(multiplicity=1) ?charge:(charge=0) ~nuclei basis_filename = +let of_filenames ?(cartesian=false) ?(multiplicity=1) ?(charge=0) ~nuclei ?(aux_basis_filenames=[]) basis_filename = let nuclei = Nuclei.of_filename nuclei in let basis = - Basis.of_nuclei_and_basis_filename ~nuclei basis_filename + Basis.of_nuclei_and_basis_filenames ~nuclei (basis_filename :: aux_basis_filenames) in lazy (make ~cartesian ~charge ~multiplicity ~nuclei basis) |> Parallel.broadcast diff --git a/Simulation.mli b/Simulation.mli index 7a4ce83..947c4ba 100644 --- a/Simulation.mli +++ b/Simulation.mli @@ -26,4 +26,6 @@ val make : val of_filenames : ?cartesian:bool -> - ?multiplicity:int -> ?charge:int -> nuclei:string -> string -> t + ?multiplicity:int -> ?charge:int -> nuclei:string -> + ?aux_basis_filenames:string list -> string -> t + diff --git a/run_fci_f12.ml b/run_fci_f12.ml new file mode 100644 index 0000000..eae08d3 --- /dev/null +++ b/run_fci_f12.ml @@ -0,0 +1,89 @@ +let () = + let open Command_line in + begin + set_header_doc (Sys.argv.(0) ^ " - QCaml command"); + set_description_doc "Runs a Hartree-Fock calculation"; + set_specs + [ { short='b' ; long="basis" ; opt=Mandatory; + arg=With_arg ""; + doc="Name of the file containing the basis set"; } ; + + { short='a' ; long="aux_basis" ; opt=Mandatory; + arg=With_arg ""; + doc="Name of the file containing the auxiliary basis set"; } ; + + { short='x' ; long="xyz" ; opt=Mandatory; + arg=With_arg ""; + doc="Name of the file containing the nuclear coordinates in xyz format"; } ; + + { short='m' ; long="multiplicity" ; opt=Optional; + arg=With_arg ""; + doc="Spin multiplicity (2S+1). Default is singlet"; } ; + + { short='c' ; long="charge" ; opt=Optional; + arg=With_arg ""; + doc="Total charge of the molecule. Default is 0"; } ; + ] + end; + + (* Handle options *) + let basis_file = Util.of_some @@ Command_line.get "basis" in + + let aux_basis_file = Util.of_some @@ Command_line.get "aux_basis" in + + let nuclei_file = Util.of_some @@ Command_line.get "xyz" in + + let charge = + match Command_line.get "charge" with + | Some x -> int_of_string x + | None -> 0 + in + + let multiplicity = + match Command_line.get "multiplicity" with + | Some x -> int_of_string x + | None -> 1 + in + + let ppf = + if Parallel.master then Format.std_formatter + else Printing.ppf_dev_null + in + + let s' = + Simulation.of_filenames ~charge ~multiplicity ~nuclei:nuclei_file basis_file + in + + let hf = HartreeFock.make s' in + Format.fprintf ppf "@[%a@]@." HartreeFock.pp_hf hf; + + let mos' = + MOBasis.of_hartree_fock hf + in + + let mo_num = + MOBasis.size mos' + in + + + let s = + Simulation.of_filenames ~charge ~multiplicity ~nuclei:nuclei_file + ~aux_basis_filenames:[aux_basis_file] basis_file + in + + let mos = + MOBasis.of_mo_basis s mos' + in + + let space = + DeterminantSpace.fci_f12_of_mo_basis mos ~frozen_core:false mo_num + in + + let ci = CI.make space in + Format.fprintf ppf "FCI energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s); + (* + let s2 = Util.xt_o_x ~o:(CI.s2_matrix ci) ~x:(CI.eigenvectors ci) in + Util.list_range 1 (DeterminantSpace.size space) + |> List.iter (fun i -> Format.printf "@[%f@]@;" s2.{i,i}); + *) +