mu_erf propagated from input

Basis/ERI.ml

@@ -36,17 +36,17 @@ module T = struct
     aux f 0 maxm;
     result
 
-  let class_of_contracted_shell_pair_couple shell_pair_couple = 
+  let class_of_contracted_shell_pair_couple ~basis shell_pair_couple = 
     let shell_p = Cspc.shell_pair_p shell_pair_couple
     and shell_q = Cspc.shell_pair_q shell_pair_couple                                         
     in
     if Array.length (Csp.shell_pairs shell_p) +
        (Array.length (Csp.shell_pairs shell_q))  < 4 then
       TwoElectronRR.contracted_class_shell_pair_couple
-        ~zero_m shell_pair_couple
+        ~basis ~zero_m shell_pair_couple
     else
       TwoElectronRRVectorized.contracted_class_shell_pairs
-        ~zero_m shell_p shell_q
+        ~basis ~zero_m shell_p shell_q
 
 end
 

Basis/ERI_lr.ml

@@ -16,9 +16,9 @@ module T = struct
 
   let zero_m z =
     let mu_erf = 
-      match z.range_separation with
+      match Basis.range_separation z.basis with
       | Some x -> x
-      | None   -> 0.5  (* TODO invalid_arg "range_separation is None" *)
+      | None   -> invalid_arg "range_separation is None" 
     in
     let m = mu_erf *. mu_erf in
     let expo_pq_inv = z.expo_p_inv +. z.expo_q_inv in
@@ -45,17 +45,17 @@ module T = struct
     aux f 0 maxm;
     result
 
-  let class_of_contracted_shell_pair_couple shell_pair_couple = 
+  let class_of_contracted_shell_pair_couple ~basis shell_pair_couple = 
     let shell_p = Cspc.shell_pair_p shell_pair_couple
     and shell_q = Cspc.shell_pair_q shell_pair_couple                                         
     in
     if Array.length (Csp.shell_pairs shell_p) +
        (Array.length (Csp.shell_pairs shell_q))  < 4 then
       TwoElectronRR.contracted_class_shell_pair_couple
-        ~zero_m shell_pair_couple
+        ~basis ~zero_m shell_pair_couple
     else
       TwoElectronRRVectorized.contracted_class_shell_pairs
-        ~zero_m shell_p shell_q
+        ~basis ~zero_m shell_p shell_q
 
 end
 

Basis/F12.ml

@@ -11,9 +11,9 @@ module T = struct
 
   let f12_factor = F12factor.gaussian_geminal 1.0
 
-  let class_of_contracted_shell_pair_couple shell_pair_couple = 
+  let class_of_contracted_shell_pair_couple ~basis shell_pair_couple = 
     let g = f12_factor.F12factor.gaussian in
-    F12RR.contracted_class_shell_pair_couple 
+    F12RR.contracted_class_shell_pair_couple ~basis 
       g.GaussianOperator.expo_g_inv
       g.GaussianOperator.coef_g
       shell_pair_couple

Basis/F12RR.ml

@@ -171,7 +171,7 @@ let rec hvrr angMom_a angMom_b angMom_c angMom_d
     aux 0. (Array.length x - 1)
 
 
-let contracted_class_shell_pair_couple expo_g_inv coef_g shell_pair_couple : float Zmap.t =
+let contracted_class_shell_pair_couple ~basis expo_g_inv coef_g shell_pair_couple : float Zmap.t =
 
     (* Pre-computation of integral class indices *)
     let class_indices = Cspc.zkey_array shell_pair_couple in

Basis/ScreenedERI.ml

@@ -76,17 +76,17 @@ module T = struct
       ) tmp_result) ;
     result
 
-  let class_of_contracted_shell_pair_couple shell_pair_couple = 
+  let class_of_contracted_shell_pair_couple ~basis shell_pair_couple = 
     let shell_p = Cspc.shell_pair_p shell_pair_couple
     and shell_q = Cspc.shell_pair_q shell_pair_couple                                         
     in
     if Array.length (Csp.shell_pairs shell_p) +
        (Array.length (Csp.shell_pairs shell_q))  < 4 then
       TwoElectronRR.contracted_class_shell_pair_couple
-        ~zero_m shell_pair_couple
+        ~basis ~zero_m shell_pair_couple
     else
       TwoElectronRRVectorized.contracted_class_shell_pairs
-        ~zero_m shell_p shell_q
+        ~basis ~zero_m shell_p shell_q
 
 end
 

Basis/TwoElectronIntegrals.ml

@@ -13,7 +13,8 @@ module Fis = FourIdxStorage
 module type TwoEI_structure =
   sig
     val name : string
-    val class_of_contracted_shell_pair_couple : ContractedShellPairCouple.t -> float Zmap.t
+    val class_of_contracted_shell_pair_couple :
+      basis:Basis.t -> ContractedShellPairCouple.t -> float Zmap.t
   end
 
 
@@ -23,12 +24,12 @@ module Make(T : TwoEI_structure) = struct
 
   let class_of_contracted_shell_pair_couple = T.class_of_contracted_shell_pair_couple 
 
-  let filter_contracted_shell_pairs ?(cutoff=integrals_cutoff) shell_pairs =
+  let filter_contracted_shell_pairs ?(cutoff=integrals_cutoff) ~basis shell_pairs =
     List.rev_map (fun pair ->
 	match Cspc.make ~cutoff pair pair with
 	| Some cspc ->
-	  let cls = class_of_contracted_shell_pair_couple cspc in
+	  let cls = class_of_contracted_shell_pair_couple ~basis cspc in
-	  (pair, Zmap.fold (fun key value accu -> max (abs_float value) accu) cls 0. )
+	  (pair, Zmap.fold (fun _key value accu -> max (abs_float value) accu) cls 0. )
 	(* TODO \sum_k |coef_k * integral_k| *)
 	| None -> (pair, -1.)
       ) shell_pairs
@@ -105,7 +106,7 @@ module Make(T : TwoEI_structure) = struct
 
     let shell_pairs =
       Csp.of_contracted_shell_array shell
-      |> filter_contracted_shell_pairs ~cutoff
+      |> filter_contracted_shell_pairs ~basis ~cutoff
     in
 
     Printf.printf "%d significant shell pairs computed in %f seconds\n"
@@ -143,7 +144,7 @@ module Make(T : TwoEI_structure) = struct
               match cspc with
               | Some cspc ->
                 let cls =
-                  class_of_contracted_shell_pair_couple cspc
+                  class_of_contracted_shell_pair_couple ~basis cspc
                 in
                 store_class ~cutoff eri_array cspc cls
               | None -> ()
@@ -211,7 +212,7 @@ module Make(T : TwoEI_structure) = struct
 
     let shell_pairs =
       Csp.of_contracted_shell_array shell
-      |> filter_contracted_shell_pairs ~cutoff
+      |> filter_contracted_shell_pairs ~basis ~cutoff
     in
 
     if Parallel.master then
@@ -254,7 +255,7 @@ module Make(T : TwoEI_structure) = struct
               match cspc with
               | Some cspc ->
                 let cls =
-                  class_of_contracted_shell_pair_couple cspc
+                  class_of_contracted_shell_pair_couple ~basis cspc
                 in
                 result := (store_class_parallel ~cutoff cspc cls) :: !result;
               | None -> ()

Basis/TwoElectronIntegrals.mli

@@ -11,7 +11,8 @@ module type TwoEI_structure =
 val name : string
 (** Name of the kind of integrals, for printing purposes. *)
 
-val class_of_contracted_shell_pair_couple : ContractedShellPairCouple.t -> float Zmap.t                                              
+val class_of_contracted_shell_pair_couple :
+  basis:Basis.t -> ContractedShellPairCouple.t -> float Zmap.t                                              
 (** Returns an integral class from a couple of contracted shells.
     The results is stored in a Zmap.
 *)
@@ -24,7 +25,7 @@ module Make : functor (T : TwoEI_structure) ->
       include module type of FourIdxStorage
 
       val filter_contracted_shell_pairs :
-        ?cutoff:float ->
+        ?cutoff:float -> basis:Basis.t ->
         ContractedShellPair.t list -> ContractedShellPair.t list
       (** Uses Schwartz screening on contracted shell pairs. *)
 

Basis/TwoElectronRR.ml

@@ -301,7 +301,7 @@ let rec hvrr_two_e
 
 
 
-let contracted_class_shell_pair_couple ~zero_m ?range_separation shell_pair_couple : float Zmap.t =
+let contracted_class_shell_pair_couple ~basis ~zero_m shell_pair_couple : float Zmap.t =
 
     let maxm = Am.to_int (Cspc.ang_mom shell_pair_couple) in
 
@@ -340,7 +340,7 @@ let contracted_class_shell_pair_couple ~zero_m ?range_separation shell_pair_coup
         let norm_pq_sq = Co.dot center_pq center_pq in
         let expo_p_inv = Psp.exponent_inv sp_ab in
         let expo_q_inv = Psp.exponent_inv sp_cd in
-        let zero = Zp.zero ?range_separation zero_m in
+        let zero = Zp.zero basis zero_m in
         let zero_m_array = zero_m
           { zero with
             maxm ; expo_p_inv ; expo_q_inv ; norm_pq_sq ;
@@ -411,7 +411,7 @@ let contracted_class_shell_pair_couple ~zero_m ?range_separation shell_pair_coup
 
 
 
-let contracted_class_atomic_shell_pair_couple ~zero_m ?range_separation atomic_shell_pair_couple : float Zmap.t =
+let contracted_class_atomic_shell_pair_couple ~basis ~zero_m atomic_shell_pair_couple : float Zmap.t =
 
     let maxm = Am.to_int (Aspc.ang_mom atomic_shell_pair_couple) in
 
@@ -454,7 +454,7 @@ let contracted_class_atomic_shell_pair_couple ~zero_m ?range_separation atomic_s
         let norm_pq_sq = Co.dot center_pq center_pq in
         let expo_q_inv = Psp.exponent_inv sp_cd in
 
-        let zero = Zp.zero ?range_separation zero_m in
+        let zero = Zp.zero basis zero_m in
         let zero_m_array = zero_m
           { zero with
             maxm ; expo_p_inv ; expo_q_inv ; norm_pq_sq ;

Basis/TwoElectronRRVectorized.ml

@@ -577,7 +577,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
 
 
 
-let contracted_class_shell_pairs ~zero_m ?range_separation ?schwartz_p ?schwartz_q shell_p shell_q : float Zmap.t =
+let contracted_class_shell_pairs ~basis ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q : float Zmap.t =
 
   let sp = Csp.shell_pairs shell_p
   and sq = Csp.shell_pairs shell_q
@@ -652,7 +652,7 @@ let contracted_class_shell_pairs ~zero_m ?range_separation ?schwartz_p ?schwartz
                         Co.dot center_pq center_pq
                       in
 
-                      let zero = Zp.zero ?range_separation zero_m in
+                      let zero = Zp.zero basis zero_m in
                       let zero_m_array = 
                         zero_m
                         {zero with
@@ -783,7 +783,7 @@ let contracted_class_shell_pairs ~zero_m ?range_separation ?schwartz_p ?schwartz
                         let norm_pq_sq =
                           x *. x +. y *. y +. z *. z
                         in
-                        let zero = Zp.zero ?range_separation zero_m in
+                        let zero = Zp.zero basis zero_m in
                         zero_m {zero with
                                 maxm ; expo_p_inv ; expo_q_inv ; norm_pq_sq ;
                                 center_pq = Coordinate.make Coordinate.{x ; y ; z} ;

Basis/Zero_m_parameters.ml

@@ -7,16 +7,14 @@ type t =
   center_pq   : Coordinate.t ;
   center_pa   : Coordinate.t ;
   center_qc   : Coordinate.t ;
-  f12_factor  : F12factor.t option;
-  range_separation : float option;
   zero_m_func : t -> float array ;
+  basis       : Basis.t ;
 }
 
-let zero ?f12_factor ?range_separation zero_m_func = 
+let zero basis zero_m_func = 
 {
   zero_m_func ;
-  f12_factor ;
+  basis ;
-  range_separation ;
   maxm=0 ; 
   expo_p_inv = 0.;
   expo_q_inv = 0.;

MOBasis/HF12.ml

@@ -208,8 +208,8 @@ let make ~frozen_core ~simulation ~mo_basis ~aux_basis_filename () =
     GeneralBasis.combine [
       general_basis ; GeneralBasis.read aux_basis_filename
     ]
-    |> Basis.of_nuclei_and_general_basis nuclei
+    |> Basis.of_nuclei_and_general_basis ~f12 nuclei
-    |> Simulation.make ~f12 ~charge ~multiplicity ~nuclei 
+    |> Simulation.make ~charge ~multiplicity ~nuclei 
   in
 
   let aux_basis = 

Simulation.ml

@@ -4,8 +4,6 @@ type t = {
   nuclei            : Nuclei.t;
   basis             : Basis.t;
   ao_basis          : AOBasis.t;
-  f12               : F12factor.t option;
-  range_separation  : float option;
   nuclear_repulsion : float;
 }
 
@@ -15,14 +13,12 @@ let electrons         t = t.electrons
 let basis             t = t.basis
 let ao_basis          t = t.ao_basis
 let nuclear_repulsion t = t.nuclear_repulsion
-let range_separation  t = t.range_separation
+let range_separation  t = Basis.range_separation t.basis 
-let f12               t = t.f12
+let f12               t = Basis.f12 t.basis
 
 let make ?cartesian:(cartesian=false)
       ?multiplicity:(multiplicity=1)
       ?charge:(charge=0)
-      ?f12
-      ?range_separation
       ~nuclei
       basis
   =
@@ -40,7 +36,7 @@ let make ?cartesian:(cartesian=false)
   in
 
   let ao_basis =
-    AOBasis.make ?f12 ~basis ~cartesian nuclei
+    AOBasis.make ~basis ~cartesian nuclei
   in
 
   let nuclear_repulsion =
@@ -48,7 +44,7 @@ let make ?cartesian:(cartesian=false)
   in
 
   {
-    charge ; basis ; nuclei ; electrons ; ao_basis ; f12 ; range_separation ;
+    charge ; basis ; nuclei ; electrons ; ao_basis ; 
     nuclear_repulsion ; 
   }
 
@@ -58,9 +54,9 @@ let of_filenames ?(cartesian=false) ?(multiplicity=1) ?(charge=0) ?f12 ?range_se
     Nuclei.of_filename nuclei
   in
   let basis =
-    Basis.of_nuclei_and_basis_filenames ~nuclei (basis_filename :: aux_basis_filenames)
+    Basis.of_nuclei_and_basis_filenames ?f12 ?range_separation ~nuclei (basis_filename :: aux_basis_filenames)
   in
 
-  lazy (make ?range_separation ?f12 ~cartesian ~charge ~multiplicity ~nuclei basis)
+  lazy (make ~cartesian ~charge ~multiplicity ~nuclei basis)
   |> Parallel.broadcast
 

Simulation.mli

@@ -28,11 +28,18 @@ val range_separation : t -> float option
 
 val make :
   ?cartesian:bool ->
-  ?multiplicity:int -> ?charge:int -> ?f12:F12factor.t -> 
+  ?multiplicity:int ->
-  ?range_separation:float -> nuclei:Nuclei.t -> Basis.t -> t
+  ?charge:int -> 
+  nuclei:Nuclei.t ->
+  Basis.t -> t
 
 val of_filenames :
   ?cartesian:bool ->
-  ?multiplicity:int -> ?charge:int -> ?f12:F12factor.t -> 
+  ?multiplicity:int ->
-  ?range_separation: float -> nuclei:string -> ?aux_basis_filenames:string list -> string -> t
+  ?charge:int ->
+  ?f12:F12factor.t -> 
+  ?range_separation: float ->
+  nuclei:string ->
+  ?aux_basis_filenames:string list ->
+  string -> t