Moved 4-idx into FourIndex.mli

Basis/Basis.ml

@@ -1,8 +1,9 @@
 type t = 
   {
-    size : int;
+    size              : int;
     contracted_shells : ContractedShell.t array ;
     atomic_shells     : AtomicShell.t array lazy_t;
+    general_basis     : GeneralBasis.t ;
   }
 
 module As = AtomicShell
@@ -10,6 +11,7 @@ module Cs = ContractedShell
 module Gb = GeneralBasis
 module Ps = PrimitiveShell
 
+let general_basis t = t.general_basis
 
 (** Returns an array of the basis set per atom *)
 let of_nuclei_and_general_basis nucl bas =
@@ -49,7 +51,8 @@ let of_nuclei_and_general_basis nucl bas =
     |> List.sort (fun x y -> compare (As.index x) (As.index y))
     |> Array.of_list
   ) in
-  { contracted_shells ; atomic_shells ; size = !index_ }
+  { contracted_shells ; atomic_shells ; size = !index_;
+    general_basis = bas }
 
 
 let size x = x.size

Basis/Basis.mli

@@ -33,6 +33,8 @@ val atomic_shells : t -> AtomicShell.t array
 val contracted_shells : t -> ContractedShell.t array
 (** Returns all the contracted basis functions. *)
 
+val general_basis : t -> GeneralBasis.t
+(** Returns the [!GeneralBasis] that was used to build the current basis. *)
 
 (** {2 Printers} *)
 

Basis/GeneralBasis.ml

@@ -78,19 +78,21 @@ let read filename =
   loop ()
 
 
-let read_many filenames =
+let combine basis_list =
   let h = Hashtbl.create 63 in
-  let l = 
-    List.map read filenames
-    |> List.concat
-  in
-  List.iter (fun (k,v) ->
+  List.concat basis_list
+  |> 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 read_many filenames =
+  List.map read filenames
+  |> combine
+
+  
 let string_of_primitive ?id prim =
   match id with
   | None -> (string_of_float prim.exponent)^" "^(string_of_float prim.coefficient)

Basis/GeneralBasis.mli

@@ -51,6 +51,12 @@ val read : string -> t
 *)
 
 
+val combine : t list -> t
+(** Combines the basis functions of each element among different basis sets.
+    Used to complement a basis with an auxiliary basis set.
+*)
+
+
 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.
@@ -58,6 +64,7 @@ val read_many : string list -> t
 *)
 
 
+
 val read_element : in_channel -> element_basis option
 (** Reads an element from the input channel [ic]. For example,
 {[

CI/CI.ml

@@ -592,12 +592,9 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
     same_spin +. opposite_spin
   in
 
-  let result = 
-    Util.stream_range 0 (Array.length psi0 - 1)
-    |> Farm.run ~ordered:false ~f:det_contribution
-    |> Util.stream_fold (+.) 0.
-  in
-  Parallel.broadcast (lazy result)
+  Util.stream_range 0 (Array.length psi0 - 1)
+  |> Farm.run ~ordered:true ~f:det_contribution
+  |> Util.stream_to_list
 
 
 let pt2_en ci =
@@ -662,6 +659,7 @@ let pt2_en ci =
 
   second_order_sum  ci  list_holes  list_particles
     i_o1_alfa  i_o1_alfa  w_alfa
+  |> List.fold_left (+.) 0.
 
 
 
@@ -692,6 +690,7 @@ let pt2_mp ci =
 
   second_order_sum  ci  list_holes  list_particles
     i_o1_alfa  i_o1_alfa  w_alfa
+  |> List.fold_left (+.) 0.
 
 
 let variance ci =
@@ -711,5 +710,6 @@ let variance ci =
 
   second_order_sum  ci  list_holes  list_particles
     i_o1_alfa  i_o1_alfa  w_alfa
+  |> List.fold_left (+.) 0.
 
 

CI/F12CI.ml

@@ -0,0 +1,102 @@
+open Lacaml.D
+
+type t = 
+  {
+    mo_basis  : MOBasis.t ;
+    aux_basis : MOBasis.t ;
+    det_space : DeterminantSpace.t ;
+    ci        : CI.t ;
+  }
+
+let ci t        =  t.ci
+let mo_basis t  =  t.mo_basis
+let det_space t =  t.det_space
+let mo_class t  =  DeterminantSpace.mo_class @@ det_space t
+
+let f12_integrals mo_basis =
+  let two_e_ints  = MOBasis.two_e_ints mo_basis in
+  ( (fun i j _ -> 0.),
+    (fun i j k l s s' ->
+       if s' = Spin.other s then
+         ERI.get_phys two_e_ints i j k l
+       else
+         (ERI.get_phys two_e_ints i j k l) -.
+         (ERI.get_phys two_e_ints i j l k)
+    ) )
+
+
+
+let h_ij mo_basis ki kj =
+  let integrals =
+    List.map (fun f -> f mo_basis)
+      [ CI.h_integrals ]
+  in
+  CIMatrixElement.make integrals ki kj 
+  |> List.hd
+
+
+let f_ij mo_basis ki kj =
+  let integrals =
+    List.map (fun f -> f mo_basis)
+      [ f12_integrals ]
+  in
+  CIMatrixElement.make integrals ki kj 
+  |> List.hd
+
+
+let dressing_vector ci =
+
+  let mo_basis  = DeterminantSpace.mo_basis ci.CI.det_space in
+
+  let i_o1_alfa = h_ij mo_basis in
+  let alfa_o2_i = f_ij mo_basis in
+
+  let w_alfa _ _ = 1. in
+
+  let mo_class = CI.mo_class ci in
+  let list_holes = List.concat
+    [ MOClass.inactive_mos mo_class ; MOClass.active_mos mo_class ]
+  and list_particles = MOClass.auxiliary_mos mo_class
+  in
+
+  CI.second_order_sum  ci  list_holes  list_particles
+    i_o1_alfa  alfa_o2_i  w_alfa
+  |> Vec.of_list
+
+
+
+
+
+
+let make ~simulation ?(frozen_core=true) ~mo_basis ~aux_basis_filename () =
+
+  let mo_num = MOBasis.size mo_basis in
+
+  (* Add auxiliary basis set *)
+  let s =
+    let charge        = Charge.to_int @@ Simulation.charge simulation 
+    and multiplicity  = Electrons.multiplicity @@ Simulation.electrons simulation
+    and nuclei        = Simulation.nuclei simulation
+    in
+    let general_basis = 
+       Basis.general_basis @@ Simulation.basis simulation
+    in
+    GeneralBasis.combine [
+         general_basis ; GeneralBasis.read aux_basis_filename
+         ]
+    |> Basis.of_nuclei_and_general_basis nuclei
+    |> Simulation.make ~charge ~multiplicity ~nuclei 
+  in
+
+  let aux_basis = 
+    MOBasis.of_mo_basis s mo_basis
+  in
+
+  let det_space = 
+     DeterminantSpace.fci_f12_of_mo_basis  aux_basis  ~frozen_core  mo_num
+  in
+
+  let ci = CI.make det_space in
+  { mo_basis ; aux_basis ; det_space ; ci }
+
+

MOBasis/MOBasis.ml

@@ -20,6 +20,7 @@ type t =
     mo_coef         : Mat.t;             (* Matrix of the MO coefficients in the AO basis *)
     eN_ints         : NucInt.t lazy_t;   (* Electron-nucleus potential integrals *)
     ee_ints         : ERI.t lazy_t;      (* Electron-electron potential integrals *)
+    f12_ints        : F12.t lazy_t;      (* F12 integrals *)
     kin_ints        : KinInt.t lazy_t;   (* Kinetic energy integrals *)
     one_e_ints      : Mat.t lazy_t;      (* Kinetic energy integrals *)
   }
@@ -37,6 +38,7 @@ let eN_ints       t = Lazy.force t.eN_ints
 let ee_ints       t = Lazy.force t.ee_ints
 let kin_ints      t = Lazy.force t.kin_ints
 let two_e_ints    t = Lazy.force t.ee_ints
+let f12_ints      t = Lazy.force t.f12_ints
 let one_e_ints    t = Lazy.force t.one_e_ints
 
 let mo_energies t = 
@@ -64,99 +66,6 @@ let ao_matrix_of_mo_matrix ~mo_coef ~ao_overlap mo_matrix =
   x_o_xt ~x:sc ~o:mo_matrix
 
 
-let four_index_transform ~mo_coef eri_ao =
-
-  let ao_num = Mat.dim1 mo_coef in
-  let mo_num = Mat.dim2 mo_coef in
-  let eri_mo =
-      ERI.create ~size:mo_num `Dense
-  in
-
-  let mo_num_2  = mo_num * mo_num in
-  let ao_num_2  = ao_num * ao_num in
-  let ao_mo_num = ao_num * mo_num in
-
-  let range_mo = list_range 1 mo_num in
-  let range_ao = list_range 1 ao_num in
-
-  let u = Mat.create mo_num_2 mo_num
-  and o = Mat.create ao_num ao_num_2 
-  and p = Mat.create ao_num_2 mo_num
-  and q = Mat.create ao_mo_num mo_num
-  in
-
-  if Parallel.master then Printf.eprintf "4-idx transformation \n%!";
-
-  let task delta =
-      Mat.fill u 0.;
-
-      List.iter (fun l ->
-          if abs_float mo_coef.{l,delta} > epsilon then
-            begin
-              let jk = ref 0 in
-              List.iter (fun k ->
-                  List.iter (fun j ->
-                      incr jk;
-                      ERI.get_chem_all_i eri_ao ~j ~k ~l
-                      |> Array.iteri (fun i x -> o.{i+1,!jk} <- x)
-                    ) range_ao
-                ) range_ao;
-              (* o_i_jk *)
-
-              let p =
-                gemm ~transa:`T ~c:p o mo_coef
-                (* p_jk_alpha  = \sum_i o_i_jk c_i_alpha *)
-              in
-              let p' = 
-                Bigarray.reshape_2 (Bigarray.genarray_of_array2 p) ao_num ao_mo_num
-                (* p_j_kalpha *)
-              in
-
-              let q =
-                gemm ~transa:`T ~c:q p' mo_coef
-                (* q_kalpha_beta  = \sum_j p_j_kalpha c_j_beta *)
-              in
-              let q' =
-                Bigarray.reshape_2 (Bigarray.genarray_of_array2 q) ao_num mo_num_2
-                (* q_k_alphabeta  = \sum_j p_j_kalpha c_j_beta *)
-              in
-
-              ignore @@
-              gemm ~transa:`T ~beta:1. ~alpha:mo_coef.{l,delta} ~c:u q' mo_coef ;
-              (* u_alphabeta_gamma = \sum_k q_k_alphabeta c_k_gamma *)
-            end
-        ) range_ao;
-      let u =
-        Bigarray.reshape 
-          (Bigarray.genarray_of_array2 u)
-          [| mo_num ; mo_num ; mo_num |]
-        |> Bigarray.array3_of_genarray 
-      in
-      let result = ref [] in
-      List.iter (fun gamma ->
-          List.iter (fun beta ->
-              List.iter (fun alpha ->
-                  let x = u.{alpha,beta,gamma} in
-                  if x <> 0. then
-                    result := (alpha, beta, gamma, delta, x) :: !result;
-                ) (list_range 1 beta) 
-            ) range_mo
-        ) (list_range 1 delta);
-      Array.of_list !result
-  in
-
-  let n = ref 0 in
-  Stream.of_list range_mo
-  |> Farm.run ~f:task ~ordered:false
-  |> Stream.iter (fun l ->
-     if Parallel.master then (Printf.eprintf "\r%d / %d%!" !n mo_num; incr n);
-     Array.iter (fun (alpha, beta, gamma, delta, x) ->
-      ERI.set_chem eri_mo alpha beta gamma delta x) l);
-
-  if Parallel.master then Printf.eprintf "\n";
-  Parallel.broadcast (lazy eri_mo) 
-
-
 let make ~simulation ~mo_type ~mo_occupation ~mo_coef ()  =
   let ao_basis =
     Si.ao_basis simulation
@@ -175,7 +84,11 @@ let make ~simulation ~mo_type ~mo_occupation ~mo_coef ()  =
   )
   and ee_ints = lazy (
     AOBasis.ee_ints ao_basis
-    |> four_index_transform ~mo_coef   
+    |> ERI.four_index_transform mo_coef   
+  )
+  and f12_ints = lazy (
+    AOBasis.f12_ints ao_basis
+    |> F12.four_index_transform mo_coef   
   )
   in
   let one_e_ints = lazy (
@@ -183,7 +96,8 @@ let make ~simulation ~mo_type ~mo_occupation ~mo_coef ()  =
         (KinInt.matrix @@ Lazy.force kin_ints) )
   in
   { simulation ; mo_type ; mo_occupation ; mo_coef ;
-    eN_ints ; ee_ints ; kin_ints ; one_e_ints }
+    eN_ints ; ee_ints ; kin_ints ; one_e_ints ;
+    f12_ints }
 
 
 

MOBasis/MOBasis.mli

@@ -46,6 +46,9 @@ val one_e_ints : t -> Mat.t
 val two_e_ints : t -> ERI.t
 (** Electron-electron repulsion integrals *)
 
+val f12_ints : t -> F12.t
+(** F12 integrals *)
+
 val size : t -> int
 (** Number of molecular orbitals in the basis *)
 

Utils/FourIdxStorage.ml

@@ -1,4 +1,6 @@
 open Util
+open Lacaml.D
+open Constants
 
 let max_index = 1 lsl 14
 
@@ -260,3 +262,101 @@ let to_list data =
   in
   append [] 
 
+
+
+
+
+let four_index_transform coef source =
+
+  let ao_num = Mat.dim1 coef in
+  let mo_num = Mat.dim2 coef in
+  let destination =
+      match source.four_index with
+      | Dense _  -> create ~size:mo_num `Dense
+      | Sparse _ -> create ~size:mo_num `Sparse
+  in
+
+  let mo_num_2  = mo_num * mo_num in
+  let ao_num_2  = ao_num * ao_num in
+  let ao_mo_num = ao_num * mo_num in
+
+  let range_mo = list_range 1 mo_num in
+  let range_ao = list_range 1 ao_num in
+
+  let u = Mat.create mo_num_2 mo_num
+  and o = Mat.create ao_num ao_num_2 
+  and p = Mat.create ao_num_2 mo_num
+  and q = Mat.create ao_mo_num mo_num
+  in
+
+  if Parallel.master then Printf.eprintf "4-idx transformation \n%!";
+
+  let task delta =
+      Mat.fill u 0.;
+
+      List.iter (fun l ->
+          if abs_float coef.{l,delta} > epsilon then
+            begin
+              let jk = ref 0 in
+              List.iter (fun k ->
+                  List.iter (fun j ->
+                      incr jk;
+                      get_chem_all_i source ~j ~k ~l
+                      |> Array.iteri (fun i x -> o.{i+1,!jk} <- x)
+                    ) range_ao
+                ) range_ao;
+              (* o_i_jk *)
+
+              let p =
+                gemm ~transa:`T ~c:p o coef
+                (* p_jk_alpha  = \sum_i o_i_jk c_i_alpha *)
+              in
+              let p' = 
+                Bigarray.reshape_2 (Bigarray.genarray_of_array2 p) ao_num ao_mo_num
+                (* p_j_kalpha *)
+              in
+
+              let q =
+                gemm ~transa:`T ~c:q p' coef
+                (* q_kalpha_beta  = \sum_j p_j_kalpha c_j_beta *)
+              in
+              let q' =
+                Bigarray.reshape_2 (Bigarray.genarray_of_array2 q) ao_num mo_num_2
+                (* q_k_alphabeta  = \sum_j p_j_kalpha c_j_beta *)
+              in
+
+              ignore @@
+              gemm ~transa:`T ~beta:1. ~alpha:coef.{l,delta} ~c:u q' coef ;
+              (* u_alphabeta_gamma = \sum_k q_k_alphabeta c_k_gamma *)
+            end
+        ) range_ao;
+      let u =
+        Bigarray.reshape 
+          (Bigarray.genarray_of_array2 u)
+          [| mo_num ; mo_num ; mo_num |]
+        |> Bigarray.array3_of_genarray 
+      in
+      let result = ref [] in
+      List.iter (fun gamma ->
+          List.iter (fun beta ->
+              List.iter (fun alpha ->
+                  let x = u.{alpha,beta,gamma} in
+                  if x <> 0. then
+                    result := (alpha, beta, gamma, delta, x) :: !result;
+                ) (list_range 1 beta) 
+            ) range_mo
+        ) (list_range 1 delta);
+      Array.of_list !result
+  in
+
+  let n = ref 0 in
+  Stream.of_list range_mo
+  |> Farm.run ~f:task ~ordered:false
+  |> Stream.iter (fun l ->
+     if Parallel.master then (Printf.eprintf "\r%d / %d%!" !n mo_num; incr n);
+     Array.iter (fun (alpha, beta, gamma, delta, x) ->
+      set_chem destination alpha beta gamma delta x) l);
+
+  if Parallel.master then Printf.eprintf "\n";
+  Parallel.broadcast (lazy destination) 
+

Utils/FourIdxStorage.mli

@@ -7,6 +7,8 @@ There are two kinds of ordering of indices:
 
 *)
 
+open Lacaml.D
+
 type t
 
 type element = (** Element for the stream *)
@@ -54,6 +56,8 @@ val to_stream : t -> element Stream.t
 val to_list : t -> element list
 (** Retrun the data structure as a list. *)
 
+val four_index_transform : Mat.t -> t -> t
+(** Perform a four-index transformation *)
 
 (** {2 I/O} *)
 

run_fci_f12.ml

@@ -2,7 +2,7 @@ let () =
   let open Command_line in
   begin
     set_header_doc (Sys.argv.(0) ^ " - QCaml command");
-    set_description_doc "Runs a Hartree-Fock calculation";
+    set_description_doc "Runs a F12-Full CI calculation";
     set_specs
       [ { short='b' ; long="basis" ; opt=Mandatory;
           arg=With_arg "<string>";
@@ -29,7 +29,7 @@ let () =
   (* 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 aux_basis_filename  = Util.of_some @@ Command_line.get "aux_basis" in
 
   let nuclei_file = Util.of_some @@ Command_line.get "xyz" in
 
@@ -50,37 +50,22 @@ let () =
     else Printing.ppf_dev_null
   in
 
-  let s' =
+  let simulation =
     Simulation.of_filenames ~charge ~multiplicity ~nuclei:nuclei_file basis_file
   in
 
-  let hf = HartreeFock.make s' in
+  let hf = HartreeFock.make simulation in
   Format.fprintf ppf "@[%a@]@." HartreeFock.pp_hf hf;
 
-  let mos' =
+  let mo_basis = 
     MOBasis.of_hartree_fock hf
   in
 
-  let mo_num =
-    MOBasis.size mos'
+  let fcif12 =
+    F12CI.make ~simulation ~frozen_core:false ~mo_basis ~aux_basis_filename ()
   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 ci = F12CI.ci fcif12 in
+  Format.fprintf ppf "FCI energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion simulation);
   (*
   let s2 = Util.xt_o_x ~o:(CI.s2_matrix ci) ~x:(CI.eigenvectors ci) in
   Util.list_range 1 (DeterminantSpace.size space)