CAS+EN2 in the presence of auxiliary basis OK

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
+
+

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. *)

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

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,
 {[

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.Inactive _
+        | MOClass.Active _
+        | MOClass.Virtual _ -> true
+        | MOClass.Auxiliary _
         | MOClass.Deleted _
         | MOClass.Core     _ -> false
-        | _ -> true
       )
     |> Array.of_list
   in

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 =
+  let det_a = f n_alfa
-    Ss.cas_of_mo_basis mo_basis n_alfa n m
+  and det_b = f n_beta
-  and det_b =
-    Ss.cas_of_mo_basis mo_basis n_beta n m
   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 =
+  let n_core = 
-    Ss.fci_of_mo_basis ~frozen_core mo_basis n_alfa
+    if frozen_core then
-  and det_b =
+      (Nuclei.small_core @@ Simulation.nuclei @@ MOBasis.simulation mo_basis) / 2
-    Ss.fci_of_mo_basis ~frozen_core mo_basis n_beta
+    else 0
   in
-  let mo_class = Ss.mo_class det_a in
+  let n, m = 
-  let determinants =
+      (n_alfa + n_beta - n_core),
-    let det_alfa = Ss.spin_determinants det_a 
+      (mo_num - n_core)
-    and det_beta = Ss.spin_determinants det_b
   in
-    let n_det_beta = Array.length det_beta in
+  let f n_alfa =
-    let n_det_alfa = Array.length det_alfa in
+    Ss.cas_of_mo_basis mo_basis n_alfa n m
-
-    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 r = 
-    let ndet = (index_start.(n_det_alfa)) in
+    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 }
     
-    Format.printf "Number of determinants : %d %d %d\n%!" 
+
-      n_det_alfa n_det_beta ndet;
+let cas_f12_of_mo_basis mo_basis n m mo_num =
-    Arbitrary {
+  let f n_alfa =
-      det_alfa ; det_beta ; det ; index_start
+    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 =

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} *)

MOBasis/MOClass.ml

@@ -4,6 +4,7 @@ type mo_class =
   | 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
 
@@ -14,6 +15,7 @@ let pp_mo_occ ppf = function
   | 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 [
+        if ncore > 0 then
           Util.list_range 1 ncore
-        |> List.map (fun i -> Core i) ;
+          |> 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 =
       ]
   )
 
+

MOBasis/MOClass.mli

@@ -6,6 +6,7 @@ type mo_class =
   | 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 
 

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

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
+

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 "<string>";
+          doc="Name of the file containing the basis set"; } ;
+
+        { short='a' ; long="aux_basis" ; opt=Mandatory;
+          arg=With_arg "<string>";
+          doc="Name of the file containing the auxiliary basis set"; } ;
+
+        { short='x' ; long="xyz" ; opt=Mandatory;
+          arg=With_arg "<string>";
+          doc="Name of the file containing the nuclear coordinates in xyz format"; } ;
+
+        { short='m' ; long="multiplicity" ; opt=Optional;
+          arg=With_arg "<int>";
+          doc="Spin multiplicity (2S+1). Default is singlet"; } ;
+
+        { short='c' ; long="charge" ; opt=Optional;
+          arg=With_arg "<int>";
+          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});
+  *)
+