Small optimizations

Basis/Contracted_shell.ml

@@ -45,14 +45,17 @@ let compute_norm_coef s =
     Angular_momentum.to_int s.totAngMom
   in
   Array.mapi (fun i alpha ->
+      let alpha_2 =
+        alpha +. alpha
+      in
       let c =
-        ((alpha +. alpha) *. pi_inv)**(1.5) *. (pow (4. *. alpha) atot)
+        (alpha_2 *. pi_inv)**(1.5) *. (pow (alpha_2 +. alpha_2) atot)
       in
       let result a = 
         let dfa = Array.map (fun j ->
-            fact (j+j) /. ( float_of_int (1 lsl j) *. fact j)
+            ( float_of_int (1 lsl j) *. fact j) /. fact (j+j) 
           ) a
-        in sqrt (c /.  (dfa.(0) *.dfa.(1) *. dfa.(2)))
+        in sqrt (c *.  dfa.(0) *.dfa.(1) *. dfa.(2))
       in
       result
     ) s.expo

Basis/TwoElectronRR.ml

@@ -22,11 +22,11 @@ let hvrr_two_e m (angMom_a, angMom_b, angMom_c, angMom_d)
 
   (** Vertical recurrence relations *)
   let rec vrr0 m angMom_a = function
-    | 0 -> zero_m_array.(m)
     | 1 -> let i = if angMom_a.(0) = 1 then 0 else if angMom_a.(1) = 1 then 1 else 2
       in expo_inv_p *.( (Coordinate.coord center_pq i) *. zero_m_array.(m+1)
                         -. expo_b *. (Coordinate.coord center_ab i) *. zero_m_array.(m) )
 
+    | 0 -> zero_m_array.(m)
     | totAngMom_a ->
       let key = Zkey.of_int_tuple (Zkey.Three
                                      (angMom_a.(0)+1, angMom_a.(1)+1, angMom_a.(2)+1) )
@@ -196,74 +196,40 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
 
   (* Compute all integrals in the shell for each pair of significant shell pairs *)
 
+  for ab=0 to (Array.length shell_p - 1) do
+    let cab = shell_p.(ab).Shell_pair.coef  in
+    let b = shell_p.(ab).Shell_pair.j in
+    for cd=0 to (Array.length shell_q - 1) do
+      let coef_prod =
+        cab *. shell_q.(cd).Shell_pair.coef
+      in
+      (** Screening on the product of coefficients *)
+      try
+        if (abs_float coef_prod) < 1.e-4*.cutoff then
+          raise NullQuartet;
+
+        let expo_pq_inv =
+          shell_p.(ab).Shell_pair.expo_inv +. shell_q.(cd).Shell_pair.expo_inv
+        in
+        let center_pq = 
+          Coordinate.(shell_p.(ab).Shell_pair.center |- shell_q.(cd).Shell_pair.center)
+        in
+        let norm_pq_sq =
+          Coordinate.dot center_pq center_pq
+        in
+
+        let zero_m_array = 
+          zero_m ~maxm ~expo_pq_inv ~norm_pq_sq
+        in
         begin
           match Contracted_shell.(totAngMom shell_a, totAngMom shell_b,
                                   totAngMom shell_c, totAngMom shell_d) with
           | Angular_momentum.(S,S,S,S) ->
-      begin
-        for ab=0 to (Array.length shell_p - 1) do
-          for cd=0 to (Array.length shell_q - 1) do
-            let coef_prod =
-              shell_p.(ab).Shell_pair.coef *. shell_q.(cd).Shell_pair.coef
-            in
-            (** Screening on the product of coefficients *)
-            try
-              if (abs_float coef_prod) < 1.e-4*.cutoff then
-                raise NullQuartet;
-
-              let expo_pq_inv =
-                shell_p.(ab).Shell_pair.expo_inv +. shell_q.(cd).Shell_pair.expo_inv
-              in
-              let center_pq = 
-                Coordinate.(shell_p.(ab).Shell_pair.center |- shell_q.(cd).Shell_pair.center)
-              in
-              let norm_pq_sq =
-                Coordinate.dot center_pq center_pq
-              in
-
-              let zero_m_array = 
-                zero_m ~maxm ~expo_pq_inv ~norm_pq_sq
-              in
-
-              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.(0) <- contracted_class.(0) +. coef_prod *. integral
-            with NullQuartet -> ()
-          done
-        done;
-      end
-
           | _ -> 
-      begin
-        for ab=0 to (Array.length shell_p - 1) do
-          let b = shell_p.(ab).Shell_pair.j in
-          for cd=0 to (Array.length shell_q - 1) do
-            try
-              let coef_prod =
-                shell_p.(ab).Shell_pair.coef *. shell_q.(cd).Shell_pair.coef
-              in
-              (** Screening on the product of coefficients *)
-              if (abs_float coef_prod) < 1.e-4*.cutoff then
-                raise NullQuartet;
-
-              let expo_pq_inv =
-                shell_p.(ab).Shell_pair.expo_inv +. shell_q.(cd).Shell_pair.expo_inv
-              in
-              let center_pq = 
-                Coordinate.(shell_p.(ab).Shell_pair.center |- shell_q.(cd).Shell_pair.center)
-              in
-              let norm_pq_sq =
-                Coordinate.dot center_pq center_pq
-              in
-
-              let zero_m_array = 
-                zero_m ~maxm ~expo_pq_inv ~norm_pq_sq
-              in
-
               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 *)
@@ -323,11 +289,10 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
                     contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral
                   with NullQuartet -> ()
                 ) class_indices
+        end
     with NullQuartet -> ()
     done
   done;
-      end
-  end;
 
   let result = 
     Zmap.create (Array.length contracted_class)

Basis/TwoElectronRRVectorized.ml

@@ -27,7 +27,6 @@ let hvrr_two_e_vector m (angMom_a, angMom_b, angMom_c, angMom_d)
 
   (** Vertical recurrence relations *)
   let rec vrr0_v m angMom_a = function
-    | 0 -> Array.mapi (fun k c -> c *. zero_m_array.(k).(m)) coef_prod
     | 1 -> let i = if angMom_a.(0) = 1 then 0 else if angMom_a.(1) = 1 then 1 else 2
       in
       let f = expo_b *. (Coordinate.coord center_ab i) in
@@ -35,6 +34,7 @@ let hvrr_two_e_vector m (angMom_a, angMom_b, angMom_c, angMom_d)
         ( (Coordinate.coord center_pq.(k) i) *. zero_m_array.(k).(m+1)
           -. f *. zero_m_array.(k).(m) ) ) coef_prod
 
+    | 0 -> Array.mapi (fun k c -> c *. zero_m_array.(k).(m)) coef_prod
     | totAngMom_a ->
       let key = Zkey.of_int_tuple (Zkey.Three
                                      (angMom_a.(0)+1, angMom_a.(1)+1, angMom_a.(2)+1) )

run_integrals.ml

@@ -11,6 +11,10 @@ let speclist = [
 ]
 
 let run ~out =
+  (*
+  let gc = Gc.get () in
+  Gc.set { gc with minor_heap_size=(262144 / 16) };
+  *)
   let out_file = 
     match out with 
     | None -> raise (Invalid_argument "Output file should be specified with -o")