Schwartz

Basis/ERI.ml

@@ -25,8 +25,13 @@ let zero_m ~maxm ~expo_pq_inv ~norm_pq_sq =
 let contracted_class  shell_a shell_b shell_c shell_d : float Zmap.t =
   TwoElectronRR.contracted_class ~zero_m shell_a shell_b shell_c shell_d 
 
+(** Compute all the integrals of a contracted class *)
+let contracted_class_shell_pairs ?schwartz_p ?schwartz_q shell_p shell_q : float Zmap.t =
+  TwoElectronRR.contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
+
 
 type n_cls = { n : int ; cls : Z.t array }
+exception NullIntegral
 
 (** Write all integrals to a file with the <ij|kl> convention *)
 let to_file ~filename basis =
@@ -38,14 +43,44 @@ let to_file ~filename basis =
   in
 
   let oc = open_out filename in
+
+  (* Pre-compute all shell pairs *)
+  let shell_pairs =
+    Array.mapi (fun i shell_a -> Array.map (fun shell_b -> 
+        Shell_pair.create_array shell_a shell_b) (Array.sub basis 0 (i+1)) ) basis
+  in
+
+  (* Pre-compute diagonal integrals for Schwartz *)
+  let schwartz = 
+    Array.map (fun pair_array -> Array.map (fun pair ->
+        let cls = 
+          contracted_class_shell_pairs pair pair
+        in
+        (cls, Zmap.fold (fun key value accu -> max (abs_float value) accu) cls 0. )
+      ) pair_array ) shell_pairs
+  in
+
   for i=0 to (Array.length basis) - 1 do
     print_int basis.(i).Contracted_shell.indice ; print_newline ();
     for j=0 to i do
+      let schwartz_p, schwartz_p_max = schwartz.(i).(j) in
+      try
+        if (schwartz_p_max < cutoff) then raise NullIntegral;
+        let
+          shell_p = shell_pairs.(i).(j)
+        in
         for k=0 to i do
           for l=0 to k do
+            let schwartz_q, schwartz_q_max = schwartz.(k).(l) in
+            try
+            if schwartz_p_max *. schwartz_q_max < cutoff *. cutoff then
+              raise NullIntegral;
+            let
+              shell_q = shell_pairs.(k).(l)
+            in
             (* Compute all the integrals of the class *)
             let cls =
-            contracted_class basis.(i) basis.(j) basis.(k) basis.(l)
+              contracted_class_shell_pairs ~schwartz_p ~schwartz_q shell_p shell_q 
             in
 
             (* Write the data in the output file *)
@@ -74,8 +109,10 @@ let to_file ~filename basis =
                       ) basis.(k).Contracted_shell.powers
                   ) basis.(j).Contracted_shell.powers
               ) basis.(i).Contracted_shell.powers;
+              with NullIntegral -> print_endline "Schwartz"; ()
           done;
         done;
+      with NullIntegral -> print_endline "Big Schwartz"; ()
     done;
   done;
   close_out oc

Basis/Shell_pair.ml

@@ -12,6 +12,8 @@ type t = {
   norm_fun : int array -> int array -> float;
   i        : int;
   j        : int;
+  shell_a  : Contracted_shell.t;
+  shell_b  : Contracted_shell.t;
 }
 
 exception Null_contribution
@@ -80,7 +82,7 @@ let create_array ?cutoff p_a p_b =
             let center_a =
               Coordinate.(center |- Contracted_shell.center p_a)
             in
-            Some { i ; j ; norm_fun ; norm ; coef ; expo ; expo_inv ; center ; center_a ; center_ab ; norm_sq }
+            Some { i ; j ; shell_a=p_a ; shell_b=p_b ; norm_fun ; norm ; coef ; expo ; expo_inv ; center ; center_a ; center_ab ; norm_sq }
           with
           | Null_contribution -> None
         )

Basis/TwoElectronRR.ml

@@ -168,13 +168,14 @@ let hvrr_two_e m (angMom_a, angMom_b, angMom_c, angMom_d)
 
 
 
+let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q : float Zmap.t =
 
-(** Computes all the two-electron integrals of the contracted shell quartet *)
+  let shell_a = shell_p.(0).Shell_pair.shell_a
-let contracted_class ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t =
+  and shell_b = shell_p.(0).Shell_pair.shell_b
-
+  and shell_c = shell_q.(0).Shell_pair.shell_a
-  let shell_p = Shell_pair.create_array shell_a shell_b
+  and shell_d = shell_q.(0).Shell_pair.shell_b
-  and shell_q = Shell_pair.create_array shell_c shell_d
+  in
-  and maxm = 
+  let maxm = 
     let open Angular_momentum in
     (to_int @@ Contracted_shell.totAngMom shell_a) + (to_int @@ Contracted_shell.totAngMom shell_b) 
     + (to_int @@ Contracted_shell.totAngMom shell_c) + (to_int @@ Contracted_shell.totAngMom shell_d)
@@ -248,32 +249,45 @@ let contracted_class ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t =
           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) =
               let a = Zkey.to_int_array Zkey.Kind_12 key in
+              let (angMomA,angMomB,angMomC,angMomD) =
                 ( [| a.(0) ;  a.(1) ;  a.(2) |], 
                   [| a.(3) ;  a.(4) ;  a.(5) |], 
                   [| a.(6) ;  a.(7) ;  a.(8) |], 
                   [| a.(9) ; a.(10) ; a.(11) |] )
               in
               try
-              (*
                 (* Schwartz screening *)
+                (*
                 let schwartz_p = 
-                  Zmap.find overlaps_p @@ Zkey.of_int_array ~kind:Zkey.Kind_6
+                  let key =
-                    [| angMomA.(0) ; angMomA.(1) ; angMomA.(2) ; 
+                  Zkey.of_int_array Zkey.Kind_12 
-                       angMomB.(0) ; angMomB.(1) ; angMomB.(2) |]
+                      [| a.(0) ;  a.(1) ;  a.(2) ;
-                  |> abs_float
+                         a.(3) ;  a.(4) ;  a.(5) ;
+                         a.(0) ;  a.(1) ;  a.(2) ;
+                         a.(3) ;  a.(4) ;  a.(5) |]
                   in
+                  match schwartz_p with
+                  | None -> 1.
+                  | Some schwartz_p -> Zmap.find schwartz_p key
+                in
+                if schwartz_p < cutoff then raise NullQuartet;
                 let schwartz_q = 
-                  Zmap.find overlaps_q @@ Zkey.of_int_array ~kind:Zkey.Kind_6
+                  let key =
-                    [| angMomC.(0) ; angMomC.(1) ; angMomC.(2) ; 
+                  Zkey.of_int_array Zkey.Kind_12 
-                       angMomD.(0) ; angMomD.(1) ; angMomD.(2) |]
+                      [| a.(6) ;  a.(7) ;  a.(8) ;
-                  |> abs_float
+                         a.(9) ; a.(10) ; a.(11) ;
+                         a.(6) ;  a.(7) ;  a.(8) ;
+                         a.(9) ; a.(10) ; a.(11) |]
                   in
-                if schwartz_p*.schwartz_q = 0. then
+                  match schwartz_q with
-                  () ; (*raise NullQuartet; *)
+                  | None -> 1.
+                  | Some schwartz_q -> Zmap.find schwartz_q key
+                in
+                if schwartz_p *. schwartz_q < cutoff2 then raise NullQuartet;
                 *)
 
+                
                 let norm = 
                   shell_p.(ab).Shell_pair.norm_fun angMomA angMomB *. shell_q.(cd).Shell_pair.norm_fun angMomC angMomD 
                 in
@@ -288,9 +302,6 @@ let contracted_class ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t =
                       map )
                 in
                 contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral
-                (*
-;if (schwartz_p*.schwartz_q < cutoff2) then Printf.printf "%e %e\n" (schwartz_p*.schwartz_q) integral;
-              *)
               with NullQuartet -> ()
             ) class_indices
       with NullQuartet -> ()
@@ -304,3 +315,12 @@ let contracted_class ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t =
 
 
 
+(** Computes all the two-electron integrals of the contracted shell quartet *)
+let contracted_class ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t =
+
+  let shell_p = Shell_pair.create_array shell_a shell_b
+  and shell_q = Shell_pair.create_array shell_c shell_d
+  in
+  contracted_class_shell_pairs ~zero_m shell_p shell_q 
+
+