Minimal screening

Basis/TwoElectronRR.ml

@@ -202,8 +202,12 @@ let contracted_class ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t =
 
     for cd=0 to (Array.length shell_q - 1)
     do
-
+      let coef_prod =
-      let d = shell_q.(cd).Shell_pair.j in
+         shell_p.(ab).Shell_pair.coef *. shell_q.(cd).Shell_pair.coef
+      in
+      (** Screening on thr product of coefficients *)
+      if (abs_float coef_prod) > cutoff then
+        begin
 
             let expo_pq_inv =
               shell_p.(ab).Shell_pair.expo_inv +. shell_q.(cd).Shell_pair.expo_inv
@@ -218,7 +222,21 @@ let contracted_class ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t =
             let zero_m_array = 
               zero_m ~maxm ~expo_pq_inv ~norm_pq_sq
             in
-      let map = Array.init maxm (fun _ -> Zmap.create 129) in
+
+              match Contracted_shell.(totAngMom shell_a, totAngMom shell_b,
+                                      totAngMom shell_c, totAngMom shell_d) with
+              | Angular_momentum.(S,S,S,S) -> Array.iteri (fun i key ->
+                  let coef_prod =
+                    shell_p.(ab).Shell_pair.coef *. shell_q.(cd).Shell_pair.coef 
+                  in
+                  let integral =
+                    zero_m_array.(0)
+                  in
+                  contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral
+                ) class_indices
+              | _ ->
+                  let d = shell_q.(cd).Shell_pair.j in
+                  let map = Array.init maxm (fun _ -> Zmap.create (Array.length class_indices)) in
                   (* Compute the integral class from the primitive shell quartet *)
                   Array.iteri (fun i key ->
                       let (angMomA,angMomB,angMomC,angMomD) =
@@ -228,7 +246,10 @@ let contracted_class ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t =
                           [| a.(6) ;  a.(7) ;  a.(8) |], 
                           [| a.(9) ; a.(10) ; a.(11) |] )
                       in
-          let integral = 
+                      let norm = 
+                        shell_p.(ab).Shell_pair.norm_fun angMomA angMomB *. shell_q.(cd).Shell_pair.norm_fun angMomC angMomD 
+                      in
+                      let integral = chop norm (fun () ->
                         ghvrr 0 (angMomA, angMomB, angMomC, angMomD)
                           (Contracted_shell.totAngMom shell_a, Contracted_shell.totAngMom shell_b,
                           Contracted_shell.totAngMom shell_c, Contracted_shell.totAngMom shell_d)
@@ -236,17 +257,11 @@ let contracted_class ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t =
                           (Contracted_shell.expo shell_b b, Contracted_shell.expo shell_d d)
                           (shell_p.(ab).Shell_pair.expo_inv, shell_q.(cd).Shell_pair.expo_inv)
                           (shell_p.(ab).Shell_pair.center_ab, shell_q.(cd).Shell_pair.center_ab, center_pq)
-              map 
+                          map )
-          in
-          let norm = 
-            shell_p.(ab).Shell_pair.norm_fun angMomA angMomB *. shell_q.(cd).Shell_pair.norm_fun angMomC angMomD 
-          in
-          let coef_prod =
-            shell_p.(ab).Shell_pair.coef *. shell_q.(cd).Shell_pair.coef *. norm
                       in
                       contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral
                 ) class_indices
-
+        end
     done
   done;
   let result = 

Utils/Constants.ml

@@ -1,4 +1,4 @@
-let cutoff = 1.e-16
+let cutoff = 1.e-15
 
 (** Constants *)
 let pi = acos (-1.)

Utils/Zmap.ml

@@ -2,3 +2,4 @@
 
 module Zmap = Hashtbl.Make(Zkey)
 include Zmap
+