Accelerated vectorized RR

Basis/TwoElectronRRVectorized.ml

@@ -458,11 +458,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
           hrr_v angMom_a angMom_b angMom_c dm totAngMom_a totAngMom_b totAngMom_c (totAngMom_d-1)
         in h1 +. f *. h2
   in
-  hrr_v 
+  hrr_v angMom_a angMom_b angMom_c angMom_d
-      (angMom_a.(0),angMom_a.(1),angMom_a.(2))
-      (angMom_b.(0),angMom_b.(1),angMom_b.(2))
-      (angMom_c.(0),angMom_c.(1),angMom_c.(2))
-      (angMom_d.(0),angMom_d.(1),angMom_d.(2))
     totAngMom_a totAngMom_b totAngMom_c totAngMom_d 
 
 
@@ -497,6 +493,11 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
     Array.make (Array.length class_indices) 0.;
   in
 
+  let monocentric =
+      shell_p.ContractedShellPair.monocentric &&
+      shell_q.ContractedShellPair.monocentric
+  in
+
   (** Screening on the product of coefficients *)
   let coef_max_p =
     Array.fold_left (fun accu x ->
@@ -671,6 +672,26 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
         in
           (* Compute the integral class from the primitive shell quartet *)
             Array.iteri (fun i key ->
+                let (angMomA,angMomB,angMomC,angMomD) =
+                  match Zkey.to_int_tuple ~kind:Zkey.Kind_12 key with
+                  | Zkey.Twelve x -> x
+                  | _             -> assert false
+                in
+                try
+                    if monocentric then
+                      begin
+                        let ax,ay,az = angMomA
+                        and bx,by,bz = angMomB
+                        and cx,cy,cz = angMomC
+                        and dx,dy,dz = angMomD
+                        in
+                        if ( ((1 land ax+bx+cx+dx)=1) ||
+                             ((1 land ay+by+cy+dy)=1) ||
+                             ((1 land az+bz+cz+dz)=1)
+                           ) then
+                          raise NullQuartet
+                      end;
+  (*
                   let a = Zkey.to_int_array Zkey.Kind_12 key in
                   let (angMomA,angMomB,angMomC,angMomD) =
                     ( [| a.(0) ;  a.(1) ;  a.(2) |], 
@@ -678,6 +699,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
                       [| a.(6) ;  a.(7) ;  a.(8) |], 
                       [| a.(9) ; a.(10) ; a.(11) |] )
                   in
+  *)
                   let integral = 
                     hvrr_two_e_vector (angMomA, angMomB, angMomC, angMomD)
                       (Contracted_shell.totAngMom shell_a, Contracted_shell.totAngMom shell_b,
@@ -690,6 +712,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
                       map_1d map_2d np nq
                   in
                   contracted_class.(i) <- contracted_class.(i) +. integral *. norm.(i)
+                with NullQuartet -> ()
               ) class_indices
 
   end;