More functional

Basis/AtomicShellPair.ml

@@ -28,12 +28,14 @@ let make ?(cutoff=Constants.epsilon) atomic_shell_a atomic_shell_b =
   let contracted_shell_pairs =
     List.map (fun s_a -> 
       List.map (fun s_b ->
-        Csp.make ~cutoff s_a s_b
+        if Cs.index s_b <= Cs.index s_a then 
+          Csp.make ~cutoff s_a s_b
+        else
+          None
       ) l_b
     ) l_a
     |> List.concat
-    |> List.filter (function None -> false | _ -> true)
+    |> list_some
-    |> List.map (function None -> assert false | Some x -> x)
   in
   match contracted_shell_pairs with
   | [] -> None

Basis/AtomicShellPairCouple.ml

@@ -1,3 +1,5 @@
+open Util
+
 type t = 
 {
   atomic_shell_pair_p: AtomicShellPair.t ;
@@ -32,8 +34,7 @@ let make ?(cutoff=Constants.epsilon) atomic_shell_pair_p atomic_shell_pair_q =
       ) (Asp.contracted_shell_pairs atomic_shell_pair_q)
     ) (Asp.contracted_shell_pairs atomic_shell_pair_p)
     |> List.concat
-    |> List.filter (function None -> false | _ -> true)
+    |> list_some
-    |> List.map (function None -> assert false | Some x -> x)
   in
   match contracted_shell_pair_couples with
   | [] -> None 

Basis/ContractedShell.ml

@@ -75,6 +75,8 @@ let size_of_shell x = Array.length x.norm_coef_scale
 
 let primitives x = x.prim
 
+let zkey_array x = Ps.zkey_array x.prim.(0)
+
 
 (** {2 Printers} *)
 

Basis/ContractedShell.mli

@@ -65,6 +65,9 @@ val norm_scales : t -> float array
 val size_of_shell : t -> int
 (** Number of contracted functions in the shell: length of {!norm_coef_scale}. *)
 
+val zkey_array : t -> Zkey.t array
+(** Returns the array of Zkeys associated with the contracted shell. *)
+
 
 (** {2 Printers} *)
 

Basis/ContractedShellPair.ml

@@ -120,11 +120,17 @@ let cmp a b =
     of contracted shells.
 *)
 let of_contracted_shell_array ?(cutoff=Constants.epsilon) basis =
-  Array.mapi (fun i shell_a ->
+  let rec loop accu = function
-    Array.mapi (fun j shell_b -> 
+  | [] -> List.rev accu
-      make ~cutoff shell_a shell_b)
+  | (s_a :: rest) as l ->
-    (Array.sub basis 0 (i+1)) 
+    let new_accu = 
-  ) basis
+      (List.map (fun s_b -> make ~cutoff s_a s_b) l) :: accu
+    in loop new_accu rest
+  in
+  loop [] (Array.to_list basis)
+  |> List.concat
+  |> list_some
+
 
 
 let equivalent x y =

Basis/ContractedShellPair.mli

@@ -22,18 +22,16 @@ val make : ?cutoff:float -> ContractedShell.t -> ContractedShell.t -> t option
     The contracted shell pair contains the only pairs of primitives for which
     the norm is greater than [cutoff]. 
 
-    If all the primitive shell pairs are not significant, the function returns
+    The function returns [None] if all the primitive shell pairs are not
-    [None].
+    significant, or if the index of the 1st primitive is smaller than the index
+    of the second primitive.
+
 *)
 
 
-val of_contracted_shell_array : ?cutoff:float -> ContractedShell.t array -> t option array array
+val of_contracted_shell_array : ?cutoff:float -> ContractedShell.t array -> t list
-(** Creates all possible contracted shell pairs from a list of contracted shells.
+(** Creates all possible contracted shell pairs from an array of contracted shells.
-    If a shell pair is not significant, sets the value to [None]:
+    Only significant shell pairs are kept.
-
-{[
-    (of_contracted_shell_array p).(i).(j) = create p.(i) p.(j)
-]}
 *)
 
 val shell_a : t -> ContractedShell.t

Basis/ContractedShellPairCouple.ml

@@ -1,3 +1,5 @@
+open Util
+
 type t = 
 {
   shell_pair_p: ContractedShellPair.t ;
@@ -35,8 +37,7 @@ let make ?(cutoff=Constants.epsilon) shell_pair_p shell_pair_q =
       ) (Csp.coefs_and_shell_pairs shell_pair_q)
     ) (Csp.coefs_and_shell_pairs shell_pair_p)
     |> List.concat
-    |> List.filter (function None -> false | _ -> true)
+    |> list_some
-    |> List.map (function None -> assert false | Some x -> x)
   in
   match coefs_and_shell_pair_couples with
   | [] -> None 

Basis/ERI.ml

@@ -7,6 +7,8 @@ open Bigarray
 type t = (float, float32_elt, fortran_layout) Bigarray.Genarray.t
 
 module Am  = AngularMomentum
+module As  = AtomicShell
+module Asp = AtomicShellPair
 module Bs  = Basis
 module Cs  = ContractedShell
 module Csp = ContractedShellPair
@@ -49,6 +51,15 @@ let contracted_class_shell_pairs_vec ?schwartz_p ?schwartz_q shell_p shell_q : f
   TwoElectronRRVectorized.contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
 
 
+let contracted_class_atomic_shell_pairs ?schwartz_p ?schwartz_q shell_p shell_q : float Zmap.t =
+  TwoElectronRR.contracted_class_atomic_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
+
+(*
+let contracted_class_atomic_shell_pairs_vec ?schwartz_p ?schwartz_q shell_p shell_q : float Zmap.t =
+  TwoElectronRRVectorized.contracted_class_atomic_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
+*)
+
+
 let cutoff2 = cutoff *. cutoff
 (*
 type n_cls = { n : int ; cls : Zkey.t array }
@@ -78,6 +89,9 @@ let of_basis basis =
 
   let n = Bs.size basis
   and shell = Bs.contracted_shells basis
+  (*TODO
+    and atomic_shells = Bs.atomic_shells basis
+  *)
   in
 
 
@@ -85,31 +99,31 @@ let of_basis basis =
   let shell_pairs =
     Csp.of_contracted_shell_array shell
   in
+  (*TODO
+    let atomic_shell_pairs =
+    Asp.of_atomic_shell_array ~cutoff atomic_shells
+    in
+  *)
 
   (* Pre-compute diagonal integrals for Schwartz *)
   let t0 = Unix.gettimeofday () in
 
   let schwartz = 
-    Array.map (fun pair_array -> Array.map (function 
+    List.map (fun pair -> 
-        | None -> (Zmap.create 0, 0.)
+        let cls = 
-        | Some pair ->
+          contracted_class_shell_pairs pair pair
-          let cls = 
+          (*TODO
-            contracted_class_shell_pairs pair pair
+            contracted_class_atomic_shell_pairs pair pair
-          in
+          *)
-          (cls, Zmap.fold (fun key value accu -> max (abs_float value) accu) cls 0. )
+        in
-      ) pair_array ) shell_pairs
+        (pair, cls, Zmap.fold (fun key value accu -> max (abs_float value) accu) cls 0. )
+      ) shell_pairs
+    |> List.filter (fun (_, _, schwartz_p_max) -> schwartz_p_max >= cutoff)
   in
 
-  let icount = ref 0 in
+  Printf.printf "%d shell pairs computed in %f seconds\n"
-  for i=0 to (Array.length shell) - 1 do
+    (List.length schwartz) (Unix.gettimeofday () -. t0);
-    print_int (Cs.index shell.(i)) ; print_newline ();
+
-    for j=0 to i do
-      let schwartz_p, schwartz_p_max = schwartz.(i).(j) in
-        if (schwartz_p_max >= cutoff) then
-          icount := !icount + 1;
-    done;
-  done;
-  Printf.printf "%d shell pairs computed in %f seconds\n" !icount (Unix.gettimeofday () -. t0);
 
   (* Group shell pairs by common pairs of atoms *)
 
@@ -125,81 +139,86 @@ let of_basis basis =
   let t0 = Unix.gettimeofday () in
   let inn = ref 0 and out = ref 0 in
 
-  for i=0 to (Array.length shell) - 1 do
+  (*TODO
-    print_int (Cs.index shell.(i)) ; print_newline ();
+    for i=0 to (Array.length atomic_shells) - 1 do
-    for j=0 to i do
+  *)
-      let schwartz_p, schwartz_p_max = schwartz.(i).(j) in
+  let ishell = ref 0 in
+  List.iter (fun (shell_p, schwartz_p, schwartz_p_max) ->
+      let () = 
+        if (Cs.index (Csp.shell_a shell_p) > !ishell) then
+          (ishell := Cs.index (Csp.shell_a shell_p) ; print_int !ishell ; print_newline ())
+      in
+
+      let sp =
+        Csp.shell_pairs shell_p
+        (*TODO
+          Asp.atomic_shell_pairs shell_p
+        *)
+      in
+
       try
-        if (schwartz_p_max < cutoff) then raise NullIntegral;
+        List.iter (fun (shell_q, schwartz_q, schwartz_q_max) ->
+            let () = 
+              if  Cs.index (Csp.shell_a shell_q) >
+                  Cs.index (Csp.shell_a shell_p) then
+                raise Exit
+            in
 
-        let shell_p = 
-          match shell_pairs.(i).(j) with
-          | None -> raise NullIntegral
-          | Some x -> x
-        in
-
-        let sp =
-          Csp.shell_pairs shell_p
-        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 < cutoff2 then
+              if schwartz_p_max *. schwartz_q_max < cutoff2 then
-              raise NullIntegral;
+                raise NullIntegral;
 
-            let shell_q =
+              let sq =
-              match shell_pairs.(k).(l) with
+                Csp.shell_pairs shell_q
-              | None -> raise NullIntegral
+                (*TODO
-              | Some x -> x
+                  Asp.atomic_shell_pairs shell_q
-            in
+                *)
+              in
 
-            let sq =
+              let swap = 
-              Csp.shell_pairs shell_q
-            in
-
-            let swap = 
                 Array.length sp > Array.length sq 
-            in
+              in
 
-            (* Compute all the integrals of the class *)
+              (* Compute all the integrals of the class *)
-            let cls =
+              let cls =
-              if swap then
+                if swap then
-                if (Array.length sp) + (Array.length sq)  < 4 then 
+                  (*TODO
-                  contracted_class_shell_pairs ~schwartz_p:schwartz_q ~schwartz_q:schwartz_p shell_q shell_p 
+                    contracted_class_atomic_shell_pairs ~schwartz_p:schwartz_q ~schwartz_q:schwartz_p shell_q shell_p 
+                  *)
+                  if (Array.length sp) + (Array.length sq)  < 4 then 
+                    contracted_class_shell_pairs ~schwartz_p:schwartz_q ~schwartz_q:schwartz_p shell_q shell_p 
+                  else
+                    contracted_class_shell_pairs_vec ~schwartz_p:schwartz_q ~schwartz_q:schwartz_p shell_q shell_p 
                 else
-                  contracted_class_shell_pairs_vec ~schwartz_p:schwartz_q ~schwartz_q:schwartz_p shell_q shell_p 
-              else
                 if (Array.length sp) + (Array.length sq)  < 4 then 
                   contracted_class_shell_pairs ~schwartz_p ~schwartz_q shell_p shell_q 
                 else
                   contracted_class_shell_pairs_vec ~schwartz_p ~schwartz_q shell_p shell_q 
-            in
+              in
 
 
-            (* Write the data in the output file *)
+              (* Write the data in the output file *)
-            Array.iteri (fun i_c powers_i ->
+              Array.iteri (fun i_c powers_i ->
-                let i_c = Cs.index shell.(i) + i_c + 1 in
+                  let i_c = Cs.index (Csp.shell_a shell_p) + i_c + 1 in
-                let xi = to_powers powers_i in
+                  let xi = to_powers powers_i in
-                Array.iteri (fun j_c powers_j ->
+                  Array.iteri (fun j_c powers_j ->
-                    let j_c = Cs.index shell.(j) + j_c + 1 in
+                      let j_c = Cs.index (Csp.shell_b shell_p) + j_c + 1 in
-                    let xj = to_powers powers_j in
+                      let xj = to_powers powers_j in
-                    Array.iteri (fun k_c powers_k ->
+                      Array.iteri (fun k_c powers_k ->
-                        let k_c = Cs.index shell.(k) + k_c + 1 in 
+                          let k_c = Cs.index (Csp.shell_a shell_q) + k_c + 1 in 
-                        let xk = to_powers powers_k in
+                          let xk = to_powers powers_k in
-                        Array.iteri (fun l_c powers_l ->
+                          Array.iteri (fun l_c powers_l ->
-                            let l_c = Cs.index shell.(l) + l_c + 1 in
+                              let l_c = Cs.index (Csp.shell_b shell_q) + l_c + 1 in
-                            let xl = to_powers powers_l in
+                              let xl = to_powers powers_l in
-                            let key = 
+                              let key = 
-                              if swap then
+                                if swap then
-                                Zkey.of_powers_twelve  xk xl xi xj 
+                                  Zkey.of_powers_twelve  xk xl xi xj 
-                              else
+                                else
-                                Zkey.of_powers_twelve  xi xj xk xl 
+                                  Zkey.of_powers_twelve  xi xj xk xl 
-                            in
+                              in
-                            let value = 
+                              let value = 
-                              Zmap.find cls key 
+                                Zmap.find cls key 
-                            in
+                              in
                               eri_array.{i_c,k_c,j_c,l_c} <- value;
                               eri_array.{j_c,k_c,i_c,l_c} <- value;
                               eri_array.{i_c,l_c,j_c,k_c} <- value;
@@ -208,21 +227,19 @@ let of_basis basis =
                               eri_array.{k_c,j_c,l_c,i_c} <- value;
                               eri_array.{l_c,i_c,k_c,j_c} <- value;
                               eri_array.{l_c,j_c,k_c,i_c} <- value;
-                            if (abs_float value > cutoff) then
+                              if (abs_float value > cutoff) then
-                              (inn := !inn + 1;
+                                (inn := !inn + 1;
-                              )
+                                )
-                            else
+                              else
-                              out := !out + 1;
+                                out := !out + 1;
-                          ) Am.(zkey_array (Singlet (Cs.ang_mom shell.(l))))
+                            ) (Cs.zkey_array (Csp.shell_b shell_q))
-                      ) Am.(zkey_array (Singlet (Cs.ang_mom shell.(k))))
+                        ) (Cs.zkey_array (Csp.shell_a shell_q))
-                  ) Am.(zkey_array (Singlet (Cs.ang_mom shell.(j))))
+                    ) (Cs.zkey_array (Csp.shell_b shell_p))
-              ) Am.(zkey_array (Singlet (Cs.ang_mom shell.(i))))
+                ) (Cs.zkey_array (Csp.shell_a shell_p))
             with NullIntegral ->  ()
-          done;
+          ) schwartz
-        done;
+      with Exit -> ()
-      with NullIntegral -> ()
+    ) schwartz;
-    done;
-  done;
   Printf.printf "In: %d  Out:%d\n" !inn !out ;
   Printf.printf "Computed ERIs in %f seconds\n%!" (Unix.gettimeofday () -. t0);
   eri_array

Basis/PrimitiveShell.ml

@@ -84,3 +84,5 @@ let norm_scales x = Lazy.force x.norm_scales
 
 let size_of_shell x = Array.length (norm_scales x)
 
+let zkey_array x = Am.(zkey_array (Singlet (x.ang_mom)))
+

Basis/PrimitiveShell.mli

@@ -61,3 +61,6 @@ val norm_scales : t -> float array
 val size_of_shell : t -> int
 (** Number of functions in the shell. *)
 
+val zkey_array : t -> Zkey.t array
+(** Returns the array of Zkeys associated with the primitive shell. *)
+

Basis/PrimitiveShellPair.ml

@@ -65,57 +65,57 @@ let create_make_of p_a p_b =
 
   function p_a -> 
 
-      let norm_coef_a =
+    let norm_coef_a =
-        Ps.normalization p_a
+      Ps.normalization p_a
+    in
+
+    let alpha_a = 
+      Co.( Ps.exponent p_a |. Ps.center p_a )
+    in
+
+    function p_b -> 
+
+      let normalization =
+        norm_coef_a *. Ps.normalization p_b
       in
 
-      let alpha_a = 
+      let exponent =
-        Co.( Ps.exponent p_a |. Ps.center p_a )
+        Ps.exponent p_a +. Ps.exponent p_b
       in
 
-      function p_b -> 
+      let exponent_inv = 1. /. exponent in
 
-          let normalization =
+      let normalization = 
-            norm_coef_a *. Ps.normalization p_b
+        let argexpo = 
+          Ps.exponent p_a *. Ps.exponent p_b *. a_minus_b_sq *. exponent_inv
+        in
+        normalization *. (pi *. exponent_inv)**1.5 *. exp (-. argexpo)
+      in
+
+      function cutoff ->
+
+        if abs_float normalization > cutoff then (
+
+          let beta_b = 
+            Co.( Ps.exponent p_b |. Ps.center p_b )
           in
 
-          let exponent =
+          let center =
-            Ps.exponent p_a +. Ps.exponent p_b
+            Co.(exponent_inv |. (alpha_a |+ beta_b))
           in
 
-          let exponent_inv = 1. /. exponent in
+          let center_minus_a = 
-
+            Co.(center |- Ps.center p_a)
-          let normalization = 
-            let argexpo = 
-              Ps.exponent p_a *. Ps.exponent p_b *. a_minus_b_sq *. exponent_inv
-            in
-            normalization *. (pi *. exponent_inv)**1.5 *. exp (-. argexpo)
           in
 
-          function cutoff ->
+          Some {
+            ang_mom ; 
+            exponent ; exponent_inv ; center ; center_minus_a ; a_minus_b ; 
+            a_minus_b_sq ; normalization ; norm_scales ; shell_a = p_a;
+            shell_b = p_b }
 
-              if abs_float normalization > cutoff then (
+        )
-
+        else None
-                  let beta_b = 
-                    Co.( Ps.exponent p_b |. Ps.center p_b )
-                  in
-
-                  let center =
-                    Co.(exponent_inv |. (alpha_a |+ beta_b))
-                  in
-
-                  let center_minus_a = 
-                    Co.(center |- Ps.center p_a)
-                  in
-
-                  Some {
-                    ang_mom ; 
-                    exponent ; exponent_inv ; center ; center_minus_a ; a_minus_b ; 
-                    a_minus_b_sq ; normalization ; norm_scales ; shell_a = p_a;
-                    shell_b = p_b }
-
-              )
-              else None
 
       
 let make p_a p_b =

Basis/PrimitiveShellPair.mli

@@ -48,8 +48,7 @@ val create_make_of : PrimitiveShell.t -> PrimitiveShell.t ->
     are pre-computed.
 
     The result is None if the normalization coefficient of the resulting
-    function is below the cutoff, given as a last argument.
+    function is below the cutoff given as a last argument.
-
 *)
 
 val ang_mom : t -> AngularMomentum.t

Utils/Util.ml

@@ -152,6 +152,15 @@ let boys_function ~maxm t =
     end
 
 
+(** {2 List functions} *)
+
+let list_some l =
+  List.filter (function None -> false | _ -> true) l
+  |> List.map (function Some x -> x | _ -> assert false)
+
+
+(** {2 Linear algebra} *)
+
 
 let array_sum a = 
    Array.fold_left ( +. ) 0. a

Utils/Util.mli

@@ -57,6 +57,10 @@ val array_product : float array -> float
 (** Returns the product of all the elements of the array *)
 
 
+(** {2 Extension of the List module} *)
+val list_some : 'a option list -> 'a list
+
+
 (** {2 Linear algebra } *)
 
 val diagonalize_symm : Lacaml.D.mat -> Lacaml.D.mat * Lacaml.D.vec