Generalized F12

Basis/AOBasis.ml

@@ -24,7 +24,7 @@ let cartesian t = t.cartesian
 
 
 
-let make ~cartesian ~basis nuclei =
+let make ~cartesian ~basis ?f12 nuclei =
 
   let overlap = lazy (
     Overlap.of_basis basis
@@ -47,8 +47,14 @@ let make ~cartesian ~basis nuclei =
   ) in
 
   let f12_ints = lazy (
-    F12.of_basis basis
+    let f12 = 
-  ) in
+      match f12 with
+      | Some f12 -> f12
+      | None -> failwith "Missing f12 factor"
+    in
+    F12.of_basis f12 basis 
+  )
+  in
 
   { basis ; overlap ; ortho ; eN_ints ; kin_ints ; ee_ints ; f12_ints ;
     cartesian ;

Basis/AOBasis.mli

@@ -31,7 +31,7 @@ val cartesian : t -> bool
 
 (** {1 Creators} *)
 
-val make : cartesian:bool -> basis:Basis.t -> Nuclei.t -> t
+val make : cartesian:bool -> basis:Basis.t ->  ?f12:F12.f12_factor ->  Nuclei.t -> t
 (** Creates the data structure for atomic orbitals from a {Basis.t} and the
     molecular geometry {Nuclei.t} *)
 

Basis/F12.ml

@@ -14,89 +14,86 @@ module Fis = FourIdxStorage
 
 include FourIdxStorage
 
-(*
+type f12_factor =
-type gaussian_geminal =
 {
   expo_s  : float ;
-  coef_g : 
+  coef_g  : float array;
-  *)
+  expo_sg : float array;
-(** Exponent of the geminal *)
+  expo_sg_inv : float array;
-let expo_s = 1.0
+}
+
+let make_gaussian_corr_factor expo_s coef_g expo_sg = 
+  let expo_sg_inv = 
+    Array.map (fun x -> 1. /. (x *. expo_s *. expo_s)) expo_sg 
+  in
+  { 
+    expo_s ; coef_g ; expo_sg ; expo_sg_inv 
+  }
 
 
-(** Slater geminal *)
+(* exp (-expo_s r) *)
-
+let gaussian_geminal expo_s =
-(** Coefficients and exponents of the Gaussian fit of the Slater Geminal*) 
+  let coef_g =
-let coef_g =
+    [| 0.3144 ; 0.3037 ; 0.1681 ; 0.09811 ; 0.06024 ; 0.03726 |];
-    [| 0.3144 ; 0.3037 ; 0.1681 ; 0.09811 ; 0.06024 ; 0.03726 |]
+  and expo_sg =
-        
+    [| 0.2209 ; 1.004 ;  3.622 ; 12.16 ;   45.87 ;  254.4 |];
-let expo_sg_inv =
+  in
-  Array.map (fun x -> 1. /. (x *. expo_s *. expo_s))
+  make_gaussian_corr_factor expo_s coef_g expo_sg
-    [| 0.2209 ; 1.004 ;  3.622 ; 12.16 ;   45.87 ;  254.4 |]
 
 
-(*
-(** r12 * Slater *)
 
-(** Coefficients and exponents of the Gaussian fit of the Slater Geminal*) 
+(** r12 * exp ( -expo_s * r) *)
-let coef_g =
+let gaussian_geminal_times_r12 expo_s =
+  let coef_g =
     [| 0.2454 ; 0.2938 ; 0.1815 ; 0.11281 ; 0.07502 ; 0.05280 |]
-        
+  and expo_sg =
-let expo_sg_inv =
-  Array.map (fun x -> 1. /. (x *. expo_s *. expo_s))
     [| 0.1824 ; 0.7118;  2.252 ; 6.474 ;   19.66 ;  77.92 |]
-*)
+  in make_gaussian_corr_factor expo_s coef_g expo_sg
 
-(*
+
-(** Coefficients and exponents of the Gaussian fit of the Slater Geminal*) 
+(* exp (-expo_s r) *)
-let coef_g =
+let gaussian_geminal' expo_s =
-[|
+  let coef_g =
+    [|
     -3.4793465193721626604883567779324948787689208984375 ;
     -0.00571703486454788484955047422886309504974633455276489257812 ;
     4.14878218728681513738365538301877677440643310546875 ;
     0.202874298181392742623785352407139725983142852783203125 ;
     0.0819187742387294803858566183407674543559551239013671875 ;
     0.04225945671351955673644695821167260874062776565551757812 ;
-|]
+    |]
-        
+  and expo_sg =
-let expo_sg_inv =
+    [|
-  Array.map (fun x -> 1. /. (x *. expo_s *. expo_s))
-[|
     0.63172472556807146570889699432882480323314666748046875; 
     26.3759196683467962429858744144439697265625; 
     0.63172102793029016876147352377302013337612152099609375; 
     7.08429025944207335641067402320913970470428466796875; 
     42.4442841447001910637482069432735443115234375;
     391.44036073596890901171718724071979522705078125 ;
-|]
+    |]
-*)
+  in make_gaussian_corr_factor expo_s coef_g expo_sg
 
 
-(*
 
-Fit of 1/r:
 
-let coef_g = [|
+let one_over_r =
+  
+  let coef_g = [|
     841.88478132     ;  70.590185207     ;  18.3616020768    ;  7.2608642093    ;
     3.57483416444    ;  2.01376031082    ;  1.24216542801    ;  0.81754348620   ;
     0.564546514023   ;  0.404228610699   ;  0.297458536575   ;  0.223321219537  ;
     0.169933732064   ;  0.130190978230   ;  0.099652303426   ;  0.075428246546  ;
     0.0555635614051  ;  0.0386791283055  ;  0.0237550435652  ;  0.010006278387  ;
     |]
-
+  and expo_sg_inv =
-let expo_sg_inv =
-  Array.map (fun x -> 1. /. (x *. expo_s *. expo_s))
     [| 84135.654509    ;  2971.58727634    ;  474.716025959    ;  130.676724560    ;
         47.3938388887   ;  20.2078651631    ;  9.5411021938     ;  4.8109546955     ;
         2.52795733067   ;  1.35894103210    ;  0.73586710268    ;  0.39557629706    ;
         0.20785895177   ;  0.104809693858   ;  0.049485682527   ;  0.021099788990   ;
         0.007652472186  ;  0.0021065225215  ;  0.0003365204879  ;  0.0000118855674  |]
-*)
+  in make_gaussian_corr_factor 1.0 coef_g expo_sg_inv
 
 
-let class_of_contracted_shell_pair_couple shell_pair_couple = 
-    F12RR.contracted_class_shell_pair_couple 
-      expo_sg_inv coef_g shell_pair_couple
 
 
 
@@ -105,13 +102,16 @@ module Zero_m = struct
   let name = "F12"
 end 
 
+let class_of_contracted_shell_pair_couple f12 shell_pair_couple = 
+    F12RR.contracted_class_shell_pair_couple 
+      f12.expo_sg_inv f12.coef_g shell_pair_couple
 
 
-let filter_contracted_shell_pairs ?(cutoff=integrals_cutoff) shell_pairs =
+let filter_contracted_shell_pairs f12 ?(cutoff=integrals_cutoff) shell_pairs =
   List.map (fun pair ->
       match Cspc.make ~cutoff pair pair with
       | Some cspc ->
-        let cls = class_of_contracted_shell_pair_couple cspc in                               
+        let cls = class_of_contracted_shell_pair_couple f12 cspc in                               
         (pair, Zmap.fold (fun key value accu -> max (abs_float value) accu) cls 0. )
       (* TODO \sum_k |coef_k * integral_k| *)
       | None -> (pair, -1.)
@@ -152,7 +152,7 @@ let store_class basis ?(cutoff=integrals_cutoff) data contracted_shell_pair_coup
     Overlap.of_basis basis
     |> Overlap.matrix
   in
-  let lambda_inv = 1. /. expo_s in
+  let lambda_inv = 1. /. f12.expo_s in
   *)
   Array.iteri (fun i_c powers_i ->
       let i_c = Cs.index (Csp.shell_a shell_p) + i_c + 1 in
@@ -183,7 +183,7 @@ let store_class basis ?(cutoff=integrals_cutoff) data contracted_shell_pair_coup
 
 
 
-let of_basis_serial basis = 
+let of_basis_serial f12 basis = 
 
   let n = Bs.size basis
   and shell = Bs.contracted_shells basis
@@ -200,7 +200,7 @@ let of_basis_serial basis =
 
   let shell_pairs =
     Csp.of_contracted_shell_array shell
-    |> filter_contracted_shell_pairs ~cutoff
+    |> filter_contracted_shell_pairs f12 ~cutoff
   in
 
   Printf.printf "%d significant shell pairs computed in %f seconds\n"
@@ -238,7 +238,7 @@ let of_basis_serial basis =
             match cspc with
             | Some cspc ->
               let cls =
-                class_of_contracted_shell_pair_couple cspc
+                class_of_contracted_shell_pair_couple f12 cspc
               in
               store_class basis ~cutoff eri_array cspc cls
             | None -> ()
@@ -258,7 +258,7 @@ let of_basis_serial basis =
 
 
 
-let of_basis_parallel basis = 
+let of_basis_parallel f12 basis = 
 
   let n = Bs.size basis
   and shell = Bs.contracted_shells basis
@@ -306,7 +306,7 @@ let of_basis_parallel basis =
 
   let shell_pairs =
     Csp.of_contracted_shell_array shell
-    |> filter_contracted_shell_pairs ~cutoff
+    |> filter_contracted_shell_pairs f12 ~cutoff
   in
 
   if Parallel.master then
@@ -348,7 +348,7 @@ let of_basis_parallel basis =
           match cspc with
           | Some cspc ->
             let cls =
-              class_of_contracted_shell_pair_couple cspc
+              class_of_contracted_shell_pair_couple f12 cspc
             in
             result := (store_class_parallel ~cutoff cspc cls) :: !result;
           | None -> ()
@@ -369,7 +369,7 @@ let of_basis_parallel basis =
     Overlap.of_basis basis
     |> Overlap.matrix
   in
-  let lambda_inv = 1. /. expo_s in
+  let lambda_inv = 1. /. f12.expo_s in
   *)
   Farm.run ~ordered:false ~f input_stream
   |> Stream.iter (fun l ->

CI/F12CI.ml

@@ -2,7 +2,6 @@ open Lacaml.D
 
 type t = 
   {
-    gamma     : float;
     mo_basis  : MOBasis.t ;
     aux_basis : MOBasis.t ;
     det_space : DeterminantSpace.t ;
@@ -55,12 +54,6 @@ let f_ij mo_basis ki kj =
   in
   CIMatrixElement.make integrals ki kj 
   |> List.hd
-  (*
-  match Determinant.degrees ki kj with
-  | (2,0) 
-  | (0,2) -> 0.5 *. gamma *. integral
-  | _     -> gamma *. integral 
-  *)
 
 
 let is_internal det_space =
@@ -97,7 +90,7 @@ let is_internal det_space =
         Bitstring.logand aux_mask beta |> Bitstring.is_zero 
 
 
-let dressing_vector gamma ~frozen_core aux_basis f12_amplitudes ci =
+let dressing_vector ~frozen_core aux_basis f12_amplitudes ci =
 
 (*
   let i_o1_alfa = h_ij aux_basis in
@@ -215,8 +208,7 @@ let dressing_vector gamma ~frozen_core aux_basis f12_amplitudes ci =
 
 let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filename () =
 
-  let gamma = -0.5 /. F12.expo_s  in
+  let f12 = Util.of_some @@ Simulation.f12 simulation in
-
   let mo_num = MOBasis.size mo_basis in
 
   (* Add auxiliary basis set *)
@@ -232,7 +224,7 @@ let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filen
       general_basis ; GeneralBasis.read aux_basis_filename
     ]
     |> Basis.of_nuclei_and_general_basis nuclei
-    |> Simulation.make ~charge ~multiplicity ~nuclei 
+    |> Simulation.make ~f12 ~charge ~multiplicity ~nuclei 
   in
 
   let aux_basis = 
@@ -257,11 +249,13 @@ let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filen
     *)
     ignore @@ MOBasis.f12_ints aux_basis;
 
+    let two_gamma_inv = -2. *. f12.F12.expo_s  in
+
     let f = fun ki kj ->
       if ki <> kj then
         (f_ij aux_basis ki kj)
       else
-        1./. gamma +. (f_ij aux_basis ki kj)
+        two_gamma_inv +. (f_ij aux_basis ki kj)
     in
     let m_F = 
       CI.create_matrix_spin f det_space
@@ -287,7 +281,7 @@ let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filen
 
     let rec iteration ?(state=1) psi = 
       let delta = 
-         dressing_vector gamma ~frozen_core aux_basis (f12_amplitudes psi) ci
+         dressing_vector ~frozen_core aux_basis (f12_amplitudes psi) ci
       in
 
       let f = 1.0 /. psi.{1,1} in
@@ -357,6 +351,6 @@ TODO SINGLE STATE HERE
 
   )
   in
-  { mo_basis ; aux_basis ; det_space ; ci ; eigensystem ; gamma }
+  { mo_basis ; aux_basis ; det_space ; ci ; eigensystem }
 
 

Simulation.ml

@@ -4,6 +4,7 @@ type t = {
   nuclei            : Nuclei.t;
   basis             : Basis.t;
   ao_basis          : AOBasis.t;
+  f12               : F12.f12_factor option;
   nuclear_repulsion : float;
 }
 
@@ -13,10 +14,12 @@ let electrons         t = t.electrons
 let basis             t = t.basis
 let ao_basis          t = t.ao_basis
 let nuclear_repulsion t = t.nuclear_repulsion
+let f12               t = t.f12
 
 let make ?cartesian:(cartesian=false)
          ?multiplicity:(multiplicity=1)
          ?charge:(charge=0)
+         ?f12
          ~nuclei
          basis
   =
@@ -34,7 +37,7 @@ let make ?cartesian:(cartesian=false)
   in
 
   let ao_basis = 
-    AOBasis.make ~basis ~cartesian nuclei
+    AOBasis.make ?f12 ~basis ~cartesian nuclei
   in
 
   let nuclear_repulsion = 
@@ -42,18 +45,18 @@ let make ?cartesian:(cartesian=false)
   in
 
   {
-    charge ; basis ; nuclei ; electrons ; ao_basis ;
+    charge ; basis ; nuclei ; electrons ; ao_basis ; f12 ;
     nuclear_repulsion ; 
   }
 
 
-let of_filenames ?(cartesian=false) ?(multiplicity=1) ?(charge=0) ~nuclei ?(aux_basis_filenames=[]) basis_filename =
+let of_filenames ?(cartesian=false) ?(multiplicity=1) ?(charge=0) ?f12 ~nuclei ?(aux_basis_filenames=[]) basis_filename =
   let nuclei =
     Nuclei.of_filename nuclei
   in
   let basis =
     Basis.of_nuclei_and_basis_filenames ~nuclei (basis_filename :: aux_basis_filenames)
   in
-  lazy (make ~cartesian ~charge ~multiplicity ~nuclei basis)
+  lazy (make ?f12 ~cartesian ~charge ~multiplicity ~nuclei basis)
   |> Parallel.broadcast
 

Simulation.mli

@@ -18,14 +18,18 @@ val ao_basis : t -> AOBasis.t
 val nuclear_repulsion : t -> float
 (** Nuclear repulsion energy *)
 
+val f12 : t -> F12.f12_factor option
+(** f12 correlation factor *)
+
 (** {1 Creation} *)
 
 val make :
   ?cartesian:bool ->
-  ?multiplicity:int -> ?charge:int -> nuclei:Nuclei.t -> Basis.t -> t
+  ?multiplicity:int -> ?charge:int -> ?f12:F12.f12_factor -> 
+  nuclei:Nuclei.t -> Basis.t -> t
 
 val of_filenames :
   ?cartesian:bool ->
-  ?multiplicity:int -> ?charge:int -> nuclei:string ->
+  ?multiplicity:int -> ?charge:int -> ?f12:F12.f12_factor -> 
-  ?aux_basis_filenames:string list -> string -> t
+  nuclei:string -> ?aux_basis_filenames:string list -> string -> t
 

Utils/FourIdxStorage.ml

@@ -358,7 +358,7 @@ let four_index_transform coef source =
   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);
+     if Parallel.master then (incr n ; Printf.eprintf "\r%d / %d%!" !n mo_num);
      Array.iter (fun (alpha, beta, gamma, delta, x) ->
       set_chem destination alpha beta gamma delta x) l);
 

run_fci_f12.ml

@@ -23,6 +23,10 @@ let () =
         { short='c' ; long="charge" ; opt=Optional;
           arg=With_arg "<int>";
           doc="Total charge of the molecule. Default is 0"; } ;
+
+        { short='e' ; long="expo" ; opt=Optional;
+          arg=With_arg "<float>";
+          doc="Exponent of the Gaussian geminal"; } ;
       ]
   end;
 
@@ -39,6 +43,12 @@ let () =
     | None   -> 0
   in
 
+  let expo = 
+    match Command_line.get "expo" with
+    | Some x -> float_of_string x
+    | None   -> 1.0
+  in
+
   let multiplicity =
     match Command_line.get "multiplicity" with
     | Some x -> int_of_string x
@@ -50,8 +60,12 @@ let () =
     else Printing.ppf_dev_null
   in
 
+  let f12 =
+    F12.gaussian_geminal expo
+  in
+
   let simulation =
-    Simulation.of_filenames ~charge ~multiplicity ~nuclei:nuclei_file basis_file
+    Simulation.of_filenames ~f12 ~charge ~multiplicity ~nuclei:nuclei_file basis_file
   in
 
   let hf = HartreeFock.make simulation in
@@ -64,7 +78,9 @@ let () =
   let fcif12 =
     F12CI.make ~simulation ~frozen_core:true ~mo_basis ~aux_basis_filename ()
   in
+
   let ci = F12CI.ci fcif12 in
   Format.fprintf ppf "FCI energy     : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion simulation);
+
   let _, e_cif12 = F12CI.eigensystem fcif12 in
   Format.fprintf ppf "FCI-F12 energy : %20.16f@." (e_cif12.{1} +. Simulation.nuclear_repulsion simulation);