Added ScreenedERI

Basis/AOBasis.ml

@@ -10,6 +10,7 @@ type t =
   kin_ints          : KinInt.t lazy_t;
   ee_ints           : ERI.t lazy_t;
   f12_ints          : F12.t lazy_t;
+  f12_over_r12_ints : ScreenedERI.t lazy_t;
   cartesian     : bool;
 }
 
@@ -20,6 +21,7 @@ let eN_ints   t = Lazy.force t.eN_ints
 let kin_ints  t = Lazy.force t.kin_ints
 let ee_ints   t = Lazy.force t.ee_ints
 let f12_ints  t = Lazy.force t.f12_ints
+let f12_over_r12_ints t = Lazy.force t.f12_over_r12_ints
 let cartesian t = t.cartesian
 
 
@@ -50,7 +52,11 @@ let make ~cartesian ~basis ?f12 nuclei =
     F12.of_basis basis 
   ) in
 
-  { basis ; overlap ; ortho ; eN_ints ; kin_ints ; ee_ints ; f12_ints ;
+  let f12_over_r12_ints = lazy (
+    ScreenedERI.of_basis basis
+  ) in
+
+  { basis ; overlap ; ortho ; eN_ints ; kin_ints ; ee_ints ; f12_ints ; f12_over_r12_ints ;
     cartesian ;
   }
 

Basis/F12.ml

@@ -14,7 +14,7 @@ module T = struct
   let class_of_contracted_shell_pair_couple shell_pair_couple = 
     let g = f12_factor.F12factor.gaussian in
     F12RR.contracted_class_shell_pair_couple 
-      g.GaussianOperator.expo_sg_inv
+      g.GaussianOperator.expo_g_inv
       g.GaussianOperator.coef_g
       shell_pair_couple
 

Basis/F12RR.ml

@@ -12,7 +12,6 @@ module Po   = Powers
 module Psp  = PrimitiveShellPair
 module Pspc = PrimitiveShellPairCouple
 module Ps   = PrimitiveShell
-module Zp   = Zero_m_parameters
 
 let cutoff = Constants.integrals_cutoff 
 let cutoff2 = cutoff *. cutoff
@@ -83,7 +82,7 @@ let rec hvrr angMom_a angMom_b angMom_c angMom_d
                 let v1 = 
                   vrr (angMom_a-1) 0
                 in
-                let a = float_of_int (angMom_a-1) in
+                let a = float_of_int_fast (angMom_a-1) in
                 let v2 = 
                   vrr (angMom_a-2) 0
                 in
@@ -92,7 +91,7 @@ let rec hvrr angMom_a angMom_b angMom_c angMom_d
                 let v1 = 
                   vrr angMom_a 0
                 in
-                let a = float_of_int angMom_a in
+                let a = float_of_int_fast angMom_a in
                 let v2 = 
                   vrr (angMom_a-1) 0
                 in
@@ -102,7 +101,7 @@ let rec hvrr angMom_a angMom_b angMom_c angMom_d
                 let v1 = 
                   vrr 0 (angMom_c-1)
                 in
-                let b = float_of_int (angMom_c-1) in
+                let b = float_of_int_fast (angMom_c-1) in
                 let v3 = 
                   vrr 0 (angMom_c-2)
                 in
@@ -112,8 +111,8 @@ let rec hvrr angMom_a angMom_b angMom_c angMom_d
                 let v1 = 
                   vrr angMom_a (angMom_c-1)
                 in
-                let a = float_of_int angMom_a in
+                let a = float_of_int_fast angMom_a in
-                let b = float_of_int (angMom_c-1) in
+                let b = float_of_int_fast (angMom_c-1) in
                 let v2 = 
                   vrr (angMom_a-1) (angMom_c-1)
                 in

Basis/GaussianOperator.ml

@@ -5,15 +5,15 @@ open Constants
 type t =
 {
   coef_g     : float array;
-  expo_sg : float array;
+  expo_g     : float array;
-  expo_sg_inv : float array;
+  expo_g_inv : float array;
 }
 
-let make coef_g expo_sg =
+let make coef_g expo_g =
-  let expo_sg_inv = 
+  let expo_g_inv = 
-      Array.map (fun x -> 1. /. x ) expo_sg 
+      Array.map (fun x -> 1. /. x ) expo_g 
   in
-  { coef_g ; expo_sg ; expo_sg_inv }
+  { coef_g ; expo_g ; expo_g_inv }
 
 
 let one_over_r =
@@ -25,13 +25,13 @@ let one_over_r =
     0.169933732064   ;  0.130190978230   ;  0.099652303426   ;  0.075428246546  ;
     0.0555635614051  ;  0.0386791283055  ;  0.0237550435652  ;  0.010006278387  ;
     |]
-  and expo_sg =
+  and expo_g =
     [| 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 coef_g expo_sg
+  in make coef_g expo_g
 
 
 

Basis/KinInt.ml

@@ -90,8 +90,8 @@ let contracted_class shell_a shell_b : float Zmap.t =
                   and s2 = ov (Po.get xyz angMomA + 1) (Po.get xyz angMomB - 1) k
                   and s3 = ov (Po.get xyz angMomA - 1) (Po.get xyz angMomB + 1) k
                   and s4 = ov (Po.get xyz angMomA + 1) (Po.get xyz angMomB + 1) k
-                  and a = float_of_int (Po.get xyz angMomA)
+                  and a = float_of_int_fast (Po.get xyz angMomA)
-                  and b = float_of_int (Po.get xyz angMomB) 
+                  and b = float_of_int_fast (Po.get xyz angMomB) 
                   in
                   0.5 *. a *. b *. s1 -. expo_a *. b *. s2 -. expo_b *. a *. s3 +. 
                     2.0 *. expo_a *. expo_b *. s4

Basis/OneElectronRR.ml

@@ -58,7 +58,7 @@ let hvrr_one_e  angMom_a  angMom_b
                   vrr amm
                 in
                 let v1 = vrr am in
-                let f3 = float_of_int amxyz *. expo_inv_p *. 0.5 in
+                let f3 = float_of_int_fast amxyz *. expo_inv_p *. 0.5 in
                 Array.init (maxsze - angMom_a.Po.tot)
                   (fun m -> f1 *. v1.(m) -.  f2 *. v1.(m+1) +.
                             f3 *. (v3.(m) -. expo_inv_p *. v3.(m+1))

Basis/Overlap_primitives.ml

@@ -10,7 +10,7 @@ let hvrr (angMom_a, angMom_b) (expo_inv_p) (center_ab, center_pa)
 
     if angMom_a > 0 then
       chop center_pa (fun () -> vrr (angMom_a-1))
-      +. chop (0.5 *. float_of_int (angMom_a-1) *. expo_inv_p) 
+      +. chop (0.5 *. float_of_int_fast (angMom_a-1) *. expo_inv_p) 
         (fun () -> vrr (angMom_a-2))
     else if angMom_a = 0 then 1.
     else 0.

Basis/ScreenedERI.ml

@@ -0,0 +1,95 @@
+(** Screened Electron-electron repulsion integrals (Yukawa potential). Required for F12/r12. *)
+
+open Constants
+open Util
+
+module Csp = ContractedShellPair
+module Cspc = ContractedShellPairCouple
+
+module T = struct
+
+  let name = "Screened electron repulsion integrals"
+
+  let f12_factor = F12factor.gaussian_geminal 1.0
+
+  open Zero_m_parameters
+
+  let zero_m z =
+    let expo_pq_inv = z.expo_p_inv +. z.expo_q_inv in
+    assert (expo_pq_inv <> 0.);
+    (*
+    let f12_factor = 
+      match z.f12_factor with
+      | Some f -> f
+      | None -> assert false
+    in
+    *)
+
+    let expo_G_inv, coef_G =
+      f12_factor.F12factor.gaussian.GaussianOperator.expo_g_inv,
+      f12_factor.F12factor.gaussian.GaussianOperator.coef_g
+    in
+
+    let expo_pq = 1. /. expo_pq_inv in
+    let maxm = z.maxm in
+
+    let result = Array.init (maxm+1) (fun _ -> 0.) in
+    array_range 0 (Array.length coef_G)
+    |> Array.iter (fun i -> 
+      let fG_inv = expo_pq_inv +. expo_G_inv.(i) in
+      let fG = 1. /. fG_inv in
+      let t =
+        if z.norm_pq_sq > integrals_cutoff then
+          z.norm_pq_sq *. (expo_pq -. fG)
+          else 0.
+      in
+      let fm = boys_function ~maxm t in
+
+      let tmp_array =
+        let result = Array.init (maxm+1) (fun k -> 1.) in
+        array_range 1 maxm
+        |> Array.iter (fun k -> result.(k) <- result.(k-1) *. expo_pq *. expo_G_inv.(i));
+        result
+      in
+
+      let tmp_result = Array.init (maxm+1) (fun _ -> 0.) in
+      let rec aux accu m = function
+        | 0 -> tmp_result.(m) <- fm.(m) *. accu
+        | j ->
+          begin
+            let f =
+              array_range 0 m
+              |> Array.fold_left (fun v k -> 
+                v +. (binom_float m k) *. tmp_array.(k) *. fm.(k)) 0. 
+            in
+            tmp_result.(m) <- accu *. f;
+            let new_accu = -. accu *. expo_pq in
+            (aux [@tailcall]) new_accu (m+1) (j-1)
+          end
+      in
+      let f =
+        two_over_sq_pi *. (sqrt expo_pq) *. fG *. expo_G_inv.(i) *.  exp (-.fG *. z.norm_pq_sq)
+      in
+      aux f 0 maxm;
+      Array.iteri (fun k v ->
+        result.(k) <- result.(k) +. coef_G.(i) *. v
+      ) tmp_result) ;
+    result
+
+  let class_of_contracted_shell_pair_couple shell_pair_couple = 
+    let shell_p = Cspc.shell_pair_p shell_pair_couple
+    and shell_q = Cspc.shell_pair_q shell_pair_couple                                         
+    in
+    if Array.length (Csp.shell_pairs shell_p) +
+       (Array.length (Csp.shell_pairs shell_q))  < 4 then
+      TwoElectronRR.contracted_class_shell_pair_couple
+        ~zero_m shell_pair_couple
+    else
+      TwoElectronRRVectorized.contracted_class_shell_pairs
+        ~zero_m shell_p shell_q
+
+end
+
+module M = TwoElectronIntegrals.Make(T) 
+include M
+

Basis/TwoElectronRR.ml

@@ -16,6 +16,7 @@ module Zp   = Zero_m_parameters
 let cutoff = Constants.integrals_cutoff 
 let cutoff2 = cutoff *. cutoff
   
+
 exception NullQuartet
 
 type four_idx_intermediates =
@@ -105,7 +106,7 @@ let rec hvrr_two_e
                   let amm = Po.decr xyz am in
                   let v3 = vrr0 amm in
                   let v1 = vrr0 am in
-                  let f3 = (float_of_int amxyz) *. expo_p_inv *. 0.5 in
+                  let f3 = (float_of_int_fast amxyz) *. expo_p_inv *. 0.5 in
                   Array.iteri (fun m _ ->
                     result.(m) <- f1 *. v1.(m+1) -. f2 *. v1.(m) 
                       +. f3 *. (v3.(m) +. expo_p_inv *. v3.(m+1)) ) result
@@ -142,7 +143,7 @@ let rec hvrr_two_e
                 begin
                   let am = Po.decr xyz angMom_a in
                   let f5 = 
-                    (float_of_int axyz) *. expo_p_inv *. expo_q_inv *. 0.5 
+                    (float_of_int_fast axyz) *. expo_p_inv *. expo_q_inv *. 0.5 
                   in
                   if (abs_float f5 > cutoff) then 
                     let v5 =
@@ -154,7 +155,7 @@ let rec hvrr_two_e
             if cmxyz > 0 then 
                 begin
                   let f3 =
-                    (float_of_int cmxyz) *. expo_q_inv *. 0.5 
+                    (float_of_int_fast cmxyz) *. expo_q_inv *. 0.5 
                   in
                   if (abs_float f3 > cutoff) ||
                       (abs_float (f3 *. expo_q_inv) > cutoff) then 
@@ -206,7 +207,7 @@ let rec hvrr_two_e
 
             let result = 
               if cmxyz < 1 then result else
-                  let f = 0.5 *. (float_of_int cmxyz) *. expo_q_inv in
+                  let f = 0.5 *. (float_of_int_fast cmxyz) *. expo_q_inv in
                   if abs_float f < cutoff then 0.  else
                       let cmm = Po.decr xyz cm in
                       let v3 = trr angMom_a cmm in
@@ -226,7 +227,7 @@ let rec hvrr_two_e
             in
             let result = 
               if axyz < 1 then result else
-                  let f = 0.5 *. (float_of_int axyz) *. expo_q_inv in
+                  let f = 0.5 *. (float_of_int_fast axyz) *. expo_q_inv in
                   if abs_float f < cutoff then result else
                       let am = Po.decr xyz angMom_a in
                       let v2 = trr am cm in
@@ -333,17 +334,16 @@ let contracted_class_shell_pair_couple ~zero_m shell_pair_couple : float Zmap.t
         and sp_cd = Pspc.shell_pair_q spc
         in
 
-        let expo_p_inv = Psp.exponent_inv sp_ab
-        in
-
         let center_pq = Co.(Psp.center sp_ab |- Psp.center sp_cd) in
         let center_pa = Psp.center_minus_a sp_ab in
         let center_qc = Psp.center_minus_a sp_cd in
         let norm_pq_sq = Co.dot center_pq center_pq in
+        let expo_p_inv = Psp.exponent_inv sp_ab in
         let expo_q_inv = Psp.exponent_inv sp_cd in
         let normalization = Psp.normalization sp_ab *. Psp.normalization sp_cd in
-        let zero_m_array = 
+        let zero = Zp.zero zero_m in
-          zero_m Zp.{
+        let zero_m_array = zero_m
+          { zero with
             maxm ; expo_p_inv ; expo_q_inv ; norm_pq_sq ;
             center_pq ; center_pa ; center_qc ; normalization ;
           }
@@ -456,8 +456,9 @@ let contracted_class_atomic_shell_pair_couple ~zero_m atomic_shell_pair_couple :
         let expo_q_inv = Psp.exponent_inv sp_cd in
         let normalization = Psp.normalization sp_ab *. Psp.normalization sp_cd in
 
-        let zero_m_array = 
+        let zero = Zp.zero zero_m in
-           zero_m Zp.{
+        let zero_m_array = zero_m
+          { zero with
             maxm ; expo_p_inv ; expo_q_inv ; norm_pq_sq ;
             center_pq ; center_pa ; center_qc ; normalization ;
           }

Basis/TwoElectronRRVectorized.ml

@@ -112,7 +112,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
           else
               begin
                 let amm = Po.decr xyz am in
-                let amxyz = float_of_int amxyz  in
+                let amxyz = float_of_int_fast amxyz  in
                 let v_amm = vrr0_v amm in
                 Array.iteri (fun l expo_inv_p_l ->
                     let f = amxyz *. expo_p_inv.(l) *. 0.5 
@@ -231,7 +231,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
             let p2 = 
               if cxyz < 2 then p1 else
                 let cmm = Po.decr xyz cm in
-                let fcm = (float_of_int (cxyz-1)) *. 0.5 in
+                let fcm = (float_of_int_fast (cxyz-1)) *. 0.5 in
                 let f1 =
                   Array.init nq (fun k ->
                       let x = fcm *. expo_q_inv.(k) in
@@ -356,7 +356,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
                     | Some p2 -> p2
                   in
                   for l=0 to np-1 do
-                    let fa = (float_of_int axyz) *. expo_p_inv.(l) *. 0.5 in
+                    let fa = (float_of_int_fast axyz) *. expo_p_inv.(l) *. 0.5 in
                     let p2_l = p2.(l)
                     and v_l = v.(l)
                     in
@@ -396,7 +396,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
           let result =
             if cmxyz < 1 then result else
               begin
-                let f = 0.5 *. (float_of_int cmxyz)  in
+                let f = 0.5 *. (float_of_int_fast cmxyz)  in
                 let cmm = Po.decr xyz cm in
                 match result, trr_v angMom_a cmm with
                 | None, None -> None
@@ -473,7 +473,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
           let result =
             if axyz < 1 then result else
               begin
-                let f = 0.5 *. (float_of_int axyz)  in
+                let f = 0.5 *. (float_of_int_fast axyz)  in
                 let am = Po.decr xyz angMom_a in
                 match result, trr_v am cm with
                 | Some result, None -> Some result
@@ -655,9 +655,11 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
                         Psp.normalization sq.(i-1)
                       in
 
+                      let zero = Zp.zero zero_m in
                       let zero_m_array = 
-                        zero_m Zp.{
+                        zero_m
-                          maxm=0 ; expo_p_inv ; expo_q_inv ; norm_pq_sq ;
+                        {zero with
+                          expo_p_inv ; expo_q_inv ; norm_pq_sq ;
                           center_pq ; center_pa ; center_qc ; normalization ;
                         }
                       in
@@ -787,7 +789,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
                         let normalization = Psp.normalization shell_ab *.
                           Psp.normalization shell_cd
                         in
-                        zero_m Zp.{
+                        let zero = Zp.zero zero_m in
+                        zero_m {zero with
                                 maxm ; expo_p_inv ; expo_q_inv ; norm_pq_sq ;
                                 center_pq = Coordinate.make Coordinate.{x ; y ; z} ;
                                 center_pa ; center_qc = center_qc_tmp.(cd) ;

Basis/Zero_m_parameters.ml

@@ -8,10 +8,14 @@ type t =
   center_pq   : Coordinate.t ;
   center_pa   : Coordinate.t ;
   center_qc   : Coordinate.t ;
+  f12_factor  : F12factor.t option;
+  zero_m_func : t -> float array ;
 }
 
-let zero = 
+let zero ?f12_factor zero_m_func = 
 {
+  zero_m_func ;
+  f12_factor ;
   maxm=0 ; 
   expo_p_inv = 0.;
   expo_q_inv = 0.;

Utils/Util.ml

@@ -45,12 +45,22 @@ let () =
   ;;
 
 
+let memo_float_of_int =
+  Array.init 64 float_of_int
+
+let float_of_int_fast i =
+    if Int.logand i 63 = i then
+      memo_float_of_int.(i)
+    else
+      float_of_int i
+
 
 let factmax = 150
 
-(* Incomplete gamma function : Int_0^x   exp(-t) t^(a-1) dt  
+(* Incomplete gamma function :
-  p:  1 / Gamma(a) * Int_0^x   exp(-t) t^(a-1) dt  
+  {% $\gamma(\alpha,x) = \int_0^x e^{-t}  t^{\alpha-1} dt$ %}
-  q:  1 / Gamma(a) * Int_x^inf exp(-t) t^(a-1) dt  
+  {% $p: \frac{1}{\Gamma(\alpha)} \int_0^x e^{-t} t^{\alpha-1} dt$ %} 
+  {% $q: \frac{1}{\Gamma(\alpha)} \int_x^\infty e^{-t} t^{\alpha-1} dt$ %} 
 
   reference - Haruhiko Okumura: C-gengo niyoru saishin algorithm jiten
   (New Algorithm handbook in C language) (Gijyutsu hyouron
@@ -123,7 +133,7 @@ let fact = function
   | i -> fact_memo.(i)
 
 
-let binom n k =
+let binom_func n k =
   (*TODO : slow function *)
   assert (n >= k);
   let rec aux n k =
@@ -133,6 +143,23 @@ let binom n k =
       aux (n-1) (k-1) + aux (n-1) k
   in aux n k
 
+let binom =
+  let memo = 
+    Array.init 64 (fun n -> Array.init 64 (fun k ->
+      if n >= k then binom_func n k else 0) )
+  in
+  fun n k ->
+    if n < 64 then 
+      begin
+        assert (n >= k);
+        memo.(n).(k)
+      end
+    else
+      binom_func n k
+
+let binom_float n k =
+  binom n k
+  |> float_of_int_fast
 
 let rec pow a = function
   | 0 -> 1.
@@ -157,6 +184,14 @@ let chop f g =
 
 (** Generalized Boys function. 
     maxm : Maximum total angular momentum
+    {% $F_m(x) = \frac{\gamma(m+1/2,x)}{2x^{m+1/2}}$ %}
+    where %{ $\gamma$ %} is the incomplete gamma function.
+
+    {% $F_0(0.) = 1$ %}
+    {% $F_0(t) = \frac{\sqrt{\pi}}{2\sqrt{t}} \text{erf} ( \sqrt{t} )$ %}
+    {% $F_m(0.) = \frac{1}{2m+1}$ %}
+    {% $F_m(t)  = \frac{\gamma{m+1/2,t}}{2t^{m+1/2}}
+    {% $F_m(t)  = \frac{ 2t\, F_{m+1}(t) + e^{-t} }{2m+1}$ %}
 *)
 let boys_function ~maxm t =
   match maxm with
@@ -182,10 +217,9 @@ let boys_function ~maxm t =
             let dm = 0.5 +. n in
             let f  = (pow t_inv (maxm+maxm+1) ) in
             match classify_float f with
+              | FP_normal -> (incomplete_gamma dm t) *. 0.5 *. f
               | FP_zero
-              | FP_subnormal
+              | FP_subnormal -> 0.
-              | FP_normal ->
-                (incomplete_gamma dm t) *. 0.5 *. f
               | _ -> invalid_arg "zero_m overflow"
           in
           let emt = exp (-. t) in

Utils/Util.mli

@@ -49,6 +49,9 @@ val fact : int -> float
 val binom : int -> int -> int
 (** Binomial coefficient. [binom n k] {% $= C_n^k$ %}. *)
 
+val binom_float : int -> int -> float
+(** Binomial coefficient. [binom n k] {% $= C_n^k$ %}. *)
+
 val pow : float -> int -> float
 (** Fast implementation of the power function for small integer powers *)
 
@@ -57,8 +60,10 @@ val chop : float -> (unit -> float) -> float
   than {!Constants.epsilon}, and return [a *. f ()].
   *)
 
-val of_some : 'a option -> 'a
+val float_of_int_fast : int -> float
+(* Faster implementation of float_of_int for small positive ints *)
 
+val of_some : 'a option -> 'a
 
 (** {2 Functions related to the Boys function} *)
 
@@ -77,6 +82,9 @@ val boys_function : maxm:int -> float -> float array
 
 (** {2 Extension of the Array module} *)
 
+val array_range : int -> int -> int array
+(** [array_range first last] returns an array [| first; first+1 ; ... ; last-1 ; last |]. *)
+
 val array_sum : float array -> float
 (** Returns the sum of all the elements of the array *)