Schwartz in VectorRR

Basis/TwoElectronRR.ml

@@ -350,7 +350,9 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
 
   let monocentric =
       shell_p.ContractedShellPair.monocentric &&
-      shell_q.ContractedShellPair.monocentric 
+      shell_q.ContractedShellPair.monocentric &&
+      Contracted_shell.center shell_p.ContractedShellPair.shell_a = 
+      Contracted_shell.center shell_q.ContractedShellPair.shell_a
   in
 
   (* Compute all integrals in the shell for each pair of significant shell pairs *)
@@ -369,10 +371,6 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
         if (abs_float coef_prod) < 1.e-3*.cutoff then
           raise NullQuartet;
 
-        let expo_pq_inv =
-          shell_p.ContractedShellPair.expo_inv.(ab) +.
-          shell_q.ContractedShellPair.expo_inv.(cd) 
-        in
         let center_pq = 
           Coordinate.(sp.(ab).ShellPair.center |- sq.(cd).ShellPair.center)
         in
@@ -380,6 +378,11 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
           Coordinate.dot center_pq center_pq
         in
 
+        let expo_pq_inv =
+          shell_p.ContractedShellPair.expo_inv.(ab) +.
+          shell_q.ContractedShellPair.expo_inv.(cd) 
+        in
+
         let zero_m_array = 
           zero_m ~maxm ~expo_pq_inv ~norm_pq_sq
         in
@@ -403,14 +406,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
                   Array.map (fun v2 -> v1 *. v2) norm_coef_scale_q) 
               |> Array.concat
             in
-            let () = 
+
-              if debug then (
-              if monocentric then
-                Printf.printf "Mono-centric\n"
-              else
-                Printf.printf "Multi-centric\n"
-              )
-            in
             (* Compute the integral class from the primitive shell quartet *)
             class_indices
             |> Array.iteri (fun i key ->
@@ -433,15 +429,13 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
                           ) then
                         raise NullQuartet
                     end;
-                  (* Schwartz screening *)
                   (*
+                  (* Schwartz screening *)
+                  if (maxm > 2) then
+                  (
                     let schwartz_p = 
-                    let key =
+                      let key = Zkey.of_int_tuple (Zkey.Twelve
-                    Zkey.of_int_array Zkey.Kind_12 
+                      (angMomA, angMomB, angMomA, angMomB) )
-                        [| a.(0) ;  a.(1) ;  a.(2) ;
-                          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.
@@ -449,20 +443,18 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
                     in
                     if schwartz_p < cutoff then raise NullQuartet;
                     let schwartz_q = 
-                    let key =
+                      let key = Zkey.of_int_tuple (Zkey.Twelve
-                    Zkey.of_int_array Zkey.Kind_12 
+                      (angMomC, angMomD, angMomC, angMomD) ) 
-                        [| a.(6) ;  a.(7) ;  a.(8) ;
-                          a.(9) ; a.(10) ; a.(11) ;
-                          a.(6) ;  a.(7) ;  a.(8) ;
-                          a.(9) ; a.(10) ; a.(11) |]
                       in
                       match schwartz_q with
                       | None -> 1.
                       | Some schwartz_q -> Zmap.find schwartz_q key
                     in
                     if schwartz_p *. schwartz_q < cutoff2 then raise NullQuartet;
+                    );
                   *)
 
+
                   let norm =  norm_coef_scale.(i) in
                   let coef_prod = coef_prod *. norm in
                   let integral = 

Basis/TwoElectronRRVectorized.ml

@@ -495,7 +495,9 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
 
   let monocentric =
       shell_p.ContractedShellPair.monocentric &&
-      shell_q.ContractedShellPair.monocentric
+      shell_q.ContractedShellPair.monocentric &&
+      Contracted_shell.center shell_p.ContractedShellPair.shell_a =
+      Contracted_shell.center shell_q.ContractedShellPair.shell_a
   in
 
   (** Screening on the product of coefficients *)
@@ -566,6 +568,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
                   let expo_pq_inv =
                     expo_inv_p.{i} +. expo_inv_q.{j}
                   in
+
                   let center_pq =
                     Coordinate.(sp.(i-1).ShellPair.center |- sq.(j-1).ShellPair.center)
                   in
@@ -691,15 +694,30 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
                            ) then
                           raise NullQuartet
                       end;
-  (*
+
-                  let a = Zkey.to_int_array Zkey.Kind_12 key in
+                   (* Schwartz screening *)
-                  let (angMomA,angMomB,angMomC,angMomD) =
+                  if (np+nq> 24) then
-                    ( [| a.(0) ;  a.(1) ;  a.(2) |], 
+                  (
-                      [| a.(3) ;  a.(4) ;  a.(5) |], 
+                    let schwartz_p = 
-                      [| a.(6) ;  a.(7) ;  a.(8) |], 
+                      let key = Zkey.of_int_tuple (Zkey.Twelve
-                      [| a.(9) ; a.(10) ; a.(11) |] )
+                      (angMomA, angMomB, angMomA, angMomB) )
                       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 = 
+                      let key = Zkey.of_int_tuple (Zkey.Twelve
+                      (angMomC, angMomD, angMomC, angMomD) ) 
+                      in
+                      match schwartz_q with
+                      | None -> 1.
+                      | Some schwartz_q -> Zmap.find schwartz_q key
+                    in
+                    if schwartz_p *. schwartz_q < cutoff2 then raise NullQuartet;
+                    );
+ 
                   let integral = 
                     hvrr_two_e_vector (angMomA, angMomB, angMomC, angMomD)
                       (Contracted_shell.totAngMom shell_a, Contracted_shell.totAngMom shell_b,

Utils/Zkey.ml

@@ -187,6 +187,8 @@ let to_int_tuple ~kind a =
   | Kind_1  -> One    ( Z.to_int a )
 
 
+let zero = Z.of_int 0
+
 include Z
 
 let to_string ~kind a =