F12 seems ok

Basis/F12.ml

@@ -21,7 +21,10 @@ type gaussian_geminal =
   coef_g : 
   *)
 (** Exponent of the geminal *)
-let expo_s = 1.
+let expo_s = 0.6
+
+
+(** Slater geminal *)
 
 (** Coefficients and exponents of the Gaussian fit of the Slater Geminal*) 
 let coef_g =
@@ -32,6 +35,42 @@ let expo_sg_inv =
     [| 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*) 
+let coef_g =
+    [| 0.2454 ; 0.2938 ; 0.1815 ; 0.11281 ; 0.07502 ; 0.05280 |]
+        
+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 |]
+*)
+
+(*
+(** Coefficients and exponents of the Gaussian fit of the Slater Geminal*) 
+let coef_g =
+[|
+ -3.4793465193721626604883567779324948787689208984375 ;
+ -0.00571703486454788484955047422886309504974633455276489257812 ;
+ 4.14878218728681513738365538301877677440643310546875 ;
+ 0.202874298181392742623785352407139725983142852783203125 ;
+ 0.0819187742387294803858566183407674543559551239013671875 ;
+ 0.04225945671351955673644695821167260874062776565551757812 ;
+|]
+        
+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 ;
+|]
+*)
+
 
 (*
 
@@ -319,6 +358,7 @@ let of_basis_parallel basis =
     else
       Fis.create ~size:0 `Dense
   in
+
   let s =
     Overlap.of_basis basis
     |> Overlap.matrix

CI/F12CI.ml

@@ -48,16 +48,13 @@ let hf_ij mo_basis ki kj =
   CIMatrixElement.make integrals ki kj 
 
 
-let f_ij gamma mo_basis ki kj =
+let f_ij mo_basis ki kj =
   let integrals =
     List.map (fun f -> f mo_basis)
       [ f12_integrals ]
   in
-  let integral = 
-    CIMatrixElement.make integrals ki kj 
-    |> List.hd
-  in
-  gamma *. integral 
+  CIMatrixElement.make integrals ki kj 
+  |> List.hd
   (*
   match Determinant.degrees ki kj with
   | (2,0) 
@@ -138,7 +135,7 @@ let dressing_vector gamma aux_basis f12_amplitudes ci =
       let h, f = 
         List.map (fun kj ->
         match hf_ij aux_basis ki kj with
-        | [ a ; b ] -> a, gamma *. b
+        | [ a ; b ] -> a, b
         | _ -> assert false ) in_dets
         |> List.split
       in
@@ -166,7 +163,7 @@ let dressing_vector gamma aux_basis f12_amplitudes ci =
   let m_F_aux = 
     List.map (fun ki ->
       List.map (fun kj ->
-        f_ij gamma aux_basis ki kj
+        f_ij aux_basis ki kj
       ) out_dets
       |> Array.of_list
     ) in_dets
@@ -192,7 +189,7 @@ let dressing_vector gamma aux_basis f12_amplitudes ci =
 
 let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filename () =
 
-  let gamma = 0.5 in
+  let gamma = 1.0 in
 
   let mo_num = MOBasis.size mo_basis in
 
@@ -236,9 +233,9 @@ let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filen
 
     let f = fun ki kj ->
       if ki <> kj then
-        (f_ij gamma aux_basis ki kj)
+        (f_ij aux_basis ki kj)
       else
-        1. +. (f_ij gamma aux_basis ki kj)
+        1./. gamma +. (f_ij aux_basis ki kj)
     in
     let m_F = 
       CI.create_matrix_spin f det_space

run_fci_f12.ml

@@ -62,7 +62,7 @@ let () =
   in
 
   let fcif12 =
-    F12CI.make ~simulation ~frozen_core:false ~mo_basis ~aux_basis_filename ()
+    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);