Merge branch 'master' of gitlab.com:scemama/QCaml

Basis/AOBasis.ml

@@ -3,26 +3,28 @@ open Util
 
 type t =
 {
-  basis             : Basis.t ;
+  basis              :  Basis.t                     ;
-  overlap           : Overlap.t lazy_t;
+  overlap            :  Overlap.t             lazy_t;
-  ortho             : Orthonormalization.t lazy_t;
+  ortho              :  Orthonormalization.t  lazy_t;
-  eN_ints           : NucInt.t lazy_t;
+  eN_ints            :  NucInt.t              lazy_t;
-  kin_ints          : KinInt.t lazy_t;
+  kin_ints           :  KinInt.t              lazy_t;
-  ee_ints           : ERI.t lazy_t;
+  ee_ints            :  ERI.t                 lazy_t;
-  f12_ints          : F12.t lazy_t;
+  ee_lr_ints         :  ERI_lr.t              lazy_t;
-  f12_over_r12_ints : ScreenedERI.t lazy_t;
+  f12_ints           :  F12.t                 lazy_t;
-  cartesian     : bool;
+  f12_over_r12_ints  :  ScreenedERI.t         lazy_t;
+  cartesian          :  bool                        ;
 }
 
-let basis     t = t.basis
+let basis        t = t.basis
-let overlap   t = Lazy.force t.overlap
+let overlap      t = Lazy.force t.overlap
-let ortho     t = Lazy.force t.ortho
+let ortho        t = Lazy.force t.ortho
-let eN_ints   t = Lazy.force t.eN_ints
+let eN_ints      t = Lazy.force t.eN_ints
-let kin_ints  t = Lazy.force t.kin_ints
+let kin_ints     t = Lazy.force t.kin_ints
-let ee_ints   t = Lazy.force t.ee_ints
+let ee_ints      t = Lazy.force t.ee_ints
-let f12_ints  t = Lazy.force t.f12_ints
+let ee_lr_ints   t = Lazy.force t.ee_lr_ints
+let f12_ints     t = Lazy.force t.f12_ints
 let f12_over_r12_ints t = Lazy.force t.f12_over_r12_ints
-let cartesian t = t.cartesian
+let cartesian    t = t.cartesian
 
 
 
@@ -48,45 +50,58 @@ let make ~cartesian ~basis ?f12 nuclei =
     ERI.of_basis basis
   ) in
 
+  let ee_lr_ints  = lazy (
+    ERI_lr.of_basis basis
+  ) in
+
   let f12_ints = lazy (
-    F12.of_basis basis 
+    F12.of_basis basis
   ) in
 
   let f12_over_r12_ints = lazy (
     ScreenedERI.of_basis basis
   ) in
 
-  { basis ; overlap ; ortho ; eN_ints ; kin_ints ; ee_ints ; f12_ints ; f12_over_r12_ints ;
+  { basis ; overlap ; ortho ; eN_ints ; kin_ints ; ee_ints ;
-    cartesian ;
+    ee_lr_ints ; f12_ints ; f12_over_r12_ints ; cartesian ;
   }
 
 
 
-let test_case t =
+let test_case name t =
 
   let check_matrix title a r =
     let a = Mat.to_array a in
-    Array.iteri (fun i x ->
+    Mat.to_array r 
-      let message =
+    |> Array.iteri (fun i x ->
-        Printf.sprintf "%s line %d" title i
+          let message =
-      in
+            Printf.sprintf "%s line %d" title i
-      Alcotest.(check (array (float 1.e-10))) message a.(i) x
+          in
-    ) (Mat.to_array r)
+          Alcotest.(check (array (float 1.e-10))) message a.(i) x
+        ) 
   in
 
   let check_eri title a r =
     let f { ERI.i_r1 ; j_r2 ; k_r1 ; l_r2 ; value } =
-      (i_r1, (j_r2, (k_r1, (l_r2, value)))) 
+      (i_r1, (j_r2, (k_r1, (l_r2, value))))
     in
-    let a = ERI.to_list a |> List.map f
+    let a = ERI.to_list a |> List.rev_map f |> List.rev in
-    and r = ERI.to_list r |> List.map f
+    let r = ERI.to_list r |> List.rev_map f |> List.rev in
+    Alcotest.(check (list (pair int (pair int (pair int (pair int (float 1.e-10))))))) "ERI" a r
+  in
+
+  let check_eri_lr title a r =
+    let f { ERI_lr.i_r1 ; j_r2 ; k_r1 ; l_r2 ; value } =
+      (i_r1, (j_r2, (k_r1, (l_r2, value))))
     in
-    Alcotest.(check (list (pair int (pair int (pair int (pair int (float 1.e-12))))))) "ERI" a r
+    let a = ERI_lr.to_list a |> List.rev_map f |> List.rev in
+    let r = ERI_lr.to_list r |> List.rev_map f |> List.rev in
+    Alcotest.(check (list (pair int (pair int (pair int (pair int (float 1.e-10))))))) "ERI_lr" a r
   in
 
   let test_overlap () =
     let reference =
-      sym_matrix_of_file "test_files/ao_overlap.ref"
+      sym_matrix_of_file ("test_files/"^name^"_overlap.ref")
     in
     let overlap =
       Lazy.force t.overlap |> Overlap.matrix
@@ -96,7 +111,7 @@ let test_case t =
 
   let test_eN_ints () =
     let reference =
-      sym_matrix_of_file "test_files/ao_nuc.ref"
+      sym_matrix_of_file ("test_files/"^name^"_nuc.ref")
     in
     let eN_ints =
       Lazy.force t.eN_ints |> NucInt.matrix
@@ -106,7 +121,7 @@ let test_case t =
 
   let test_kin_ints () =
     let reference =
-      sym_matrix_of_file "test_files/ao_kin.ref"
+      sym_matrix_of_file ("test_files/"^name^"_kin.ref")
     in
     let kin_ints =
       Lazy.force t.kin_ints |> KinInt.matrix
@@ -116,19 +131,32 @@ let test_case t =
 
   let test_ee_ints () =
     let reference =
-      ERI.of_file "test_files/ao_eri.ref" ~sparsity:`Dense ~size:(Basis.size t.basis)
+      ERI.of_file ("test_files/"^name^"_eri.ref") ~sparsity:`Dense ~size:(Basis.size t.basis)
     in
     let ee_ints =
-      Lazy.force t.ee_ints 
+      Lazy.force t.ee_ints
     in
     check_eri "ee_ints" ee_ints reference
+    ;
+  in
+
+  let test_ee_lr_ints () =
+    let reference =
+      ERI_lr.of_file ("test_files/"^name^"_eri_lr.ref") ~sparsity:`Dense
+          ~size:(Basis.size t.basis)
+    in
+    let ee_lr_ints =
+      Lazy.force t.ee_lr_ints
+    in
+    check_eri_lr "ee_lr_ints" ee_lr_ints reference
   in
 
   [
-    "Overlap", `Quick, test_overlap;
+    "Overlap",     `Quick,  test_overlap;
-    "eN_ints", `Quick, test_eN_ints;
+    "eN_ints",     `Quick,  test_eN_ints;
-    "kin_ints", `Quick, test_kin_ints;
+    "kin_ints",    `Quick,  test_kin_ints;
-    "ee_ints", `Quick, test_ee_ints;
+    "ee_ints",     `Quick,  test_ee_ints;
+    "ee_lr_ints",  `Quick,  test_ee_lr_ints;
   ]
 
 

Basis/AOBasis.mli

@@ -1,31 +1,34 @@
 (** Data structure for Atomic Orbitals. *)
 
-type t 
+type t
 
 (** {1 Accessors} *)
 
-val basis : t -> Basis.t 
+val basis : t -> Basis.t
 (** One-electron basis set *)
 
-val overlap : t -> Overlap.t 
+val overlap : t -> Overlap.t
 (** Overlap matrix *)
 
-val ortho : t -> Orthonormalization.t 
+val ortho : t -> Orthonormalization.t
 (** Orthonormalization matrix of the overlap *)
 
-val eN_ints : t -> NucInt.t 
+val eN_ints : t -> NucInt.t
 (** Electron-nucleus potential integrals *)
 
-val ee_ints : t -> ERI.t 
+val ee_ints : t -> ERI.t
 (** Electron-electron potential integrals *)
 
-val f12_ints : t -> F12.t 
+val ee_lr_ints : t -> ERI_lr.t
+(** Electron-electron long-range potential integrals *)
+
+val f12_ints : t -> F12.t
 (** Electron-electron potential integrals *)
 
-val kin_ints : t -> KinInt.t 
+val kin_ints : t -> KinInt.t
 (** Kinetic energy integrals *)
 
-val cartesian : t -> bool 
+val cartesian : t -> bool
 (** If true, use cartesian Gaussians (6d, 10f, ...) *)
 
 
@@ -38,4 +41,4 @@ val make : cartesian:bool -> basis:Basis.t ->  ?f12:F12factor.t ->  Nuclei.t ->
 
 (** {2 Tests} *)
 
-val test_case : t -> unit Alcotest.test_case list
+val test_case : string -> t -> unit Alcotest.test_case list

Basis/AtomicShellPair.ml

@@ -26,7 +26,7 @@ let make ?(cutoff=Constants.epsilon) atomic_shell_a atomic_shell_b =
   in
 
   let contracted_shell_pairs =
-    List.map (fun s_a -> 
+    List.concat_map (fun s_a -> 
       List.map (fun s_b ->
         if Cs.index s_b <= Cs.index s_a then 
           Csp.make ~cutoff s_a s_b
@@ -34,7 +34,6 @@ let make ?(cutoff=Constants.epsilon) atomic_shell_a atomic_shell_b =
           None
       ) l_b
     ) l_a
-    |> List.concat
     |> list_some
   in
   match contracted_shell_pairs with

Basis/AtomicShellPairCouple.ml

@@ -28,12 +28,11 @@ let make ?(cutoff=Constants.epsilon) atomic_shell_pair_p atomic_shell_pair_q =
   and atomic_shell_d = Asp.atomic_shell_b atomic_shell_pair_q
   in
   let contracted_shell_pair_couples =
-    List.map (fun ap_ab -> 
+    List.concat_map (fun ap_ab -> 
       List.map (fun ap_cd -> 
         ContractedShellPairCouple.make ~cutoff ap_ab ap_cd
       ) (Asp.contracted_shell_pairs atomic_shell_pair_q)
     ) (Asp.contracted_shell_pairs atomic_shell_pair_p)
-    |> List.concat
     |> list_some
   in
   match contracted_shell_pair_couples with

Basis/ContractedShellPairCouple.ml

@@ -29,14 +29,13 @@ let make ?(cutoff=Constants.epsilon) shell_pair_p shell_pair_q =
   in
   let cutoff = 1.e-3 *. cutoff in
   let coefs_and_shell_pair_couples =
-    List.map (fun (c_ab, sp_ab) -> 
+    List.concat_map (fun (c_ab, sp_ab) -> 
       List.map (fun (c_cd, sp_cd) -> 
         let coef_prod = c_ab *. c_cd in
         if abs_float coef_prod < cutoff then None
         else Some (coef_prod, Pspc.make sp_ab sp_cd)
       ) (Csp.coefs_and_shell_pairs shell_pair_q)
     ) (Csp.coefs_and_shell_pairs shell_pair_p)
-    |> List.concat
     |> list_some
   in
   match coefs_and_shell_pair_couples with

Basis/ERI.ml

@@ -15,10 +15,10 @@ module T = struct
   let zero_m z =
     let expo_pq_inv = z.expo_p_inv +. z.expo_q_inv in
     assert (expo_pq_inv <> 0.);
-    let exp_pq = 1. /. expo_pq_inv in
+    let expo_pq = 1. /. expo_pq_inv in
     let t =
       if z.norm_pq_sq > integrals_cutoff then
-        z.norm_pq_sq *. exp_pq
+        z.norm_pq_sq *. expo_pq
         else 0.
     in
     let maxm = z.maxm in
@@ -28,11 +28,11 @@ module T = struct
       | l -> 
         begin
           result.(k) <- result.(k) *. accu;
-          let new_accu = -. accu *. exp_pq in
+          let new_accu = -. accu *. expo_pq in
           (aux [@tailcall]) new_accu (k+1) (l-1)
         end
     in
-    let f = two_over_sq_pi *. (sqrt exp_pq) in
+    let f = two_over_sq_pi *. (sqrt expo_pq) in
     aux f 0 maxm;
     result
 

Basis/TwoElectronIntegrals.ml

@@ -24,7 +24,7 @@ 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 =
-    List.map (fun pair ->
+    List.rev_map (fun pair ->
 	match Cspc.make ~cutoff pair pair with
 	| Some cspc ->
 	  let cls = class_of_contracted_shell_pair_couple cspc in
@@ -33,20 +33,20 @@ module Make(T : TwoEI_structure) = struct
 	| None -> (pair, -1.)
       ) shell_pairs
     |> List.filter (fun (_, schwartz_p_max) -> schwartz_p_max >= cutoff)
-    |> List.map fst
+    |> List.rev_map fst
 
 
 (*  TODO                                                                                      
     let filter_contracted_shell_pair_couples
     ?(cutoff=integrals_cutoff) shell_pair_couples = 
-    List.map (fun pair -> 
+    List.rev_map (fun pair -> 
         let cls = 
           class_of_contracted_shell_pairs pair pair
         in
         (pair, 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)
-    |> List.map fst
+    |> List.rev_map fst
 
 *)
 

Basis/TwoElectronRR.ml

@@ -301,7 +301,7 @@ let rec hvrr_two_e
 
 
 
-let contracted_class_shell_pair_couple ~zero_m shell_pair_couple : float Zmap.t =
+let contracted_class_shell_pair_couple ~zero_m ?mu_erf shell_pair_couple : float Zmap.t =
 
     let maxm = Am.to_int (Cspc.ang_mom shell_pair_couple) in
 
@@ -340,12 +340,11 @@ let contracted_class_shell_pair_couple ~zero_m shell_pair_couple : float Zmap.t
         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 normalization = Psp.normalization sp_ab *. Psp.normalization sp_cd in
+        let zero = Zp.zero ?mu_erf zero_m in
-        let zero = Zp.zero zero_m in
         let zero_m_array = zero_m
           { zero with
             maxm ; expo_p_inv ; expo_q_inv ; norm_pq_sq ;
-            center_pq ; center_pa ; center_qc ; normalization ;
+            center_pq ; center_pa ; center_qc ; 
           }
         in
 
@@ -412,7 +411,7 @@ let contracted_class_shell_pair_couple ~zero_m shell_pair_couple : float Zmap.t
 
 
 
-let contracted_class_atomic_shell_pair_couple ~zero_m atomic_shell_pair_couple : float Zmap.t =
+let contracted_class_atomic_shell_pair_couple ~zero_m ?mu_erf atomic_shell_pair_couple : float Zmap.t =
 
     let maxm = Am.to_int (Aspc.ang_mom atomic_shell_pair_couple) in
 
@@ -454,13 +453,12 @@ let contracted_class_atomic_shell_pair_couple ~zero_m atomic_shell_pair_couple :
         let center_pa = Psp.center_minus_a sp_ab in
         let norm_pq_sq = Co.dot center_pq center_pq in
         let expo_q_inv = Psp.exponent_inv sp_cd in
-        let normalization = Psp.normalization sp_ab *. Psp.normalization sp_cd in
 
-        let zero = Zp.zero zero_m in
+        let zero = Zp.zero ?mu_erf zero_m in
         let zero_m_array = zero_m
           { zero with
             maxm ; expo_p_inv ; expo_q_inv ; norm_pq_sq ;
-            center_pq ; center_pa ; center_qc ; normalization ;
+            center_pq ; center_pa ; center_qc ; 
           }
         in
 

Basis/TwoElectronRRVectorized.ml

@@ -24,7 +24,7 @@ let at_least_one_valid arr =
     Array.iter (fun x -> if (abs_float x > cutoff) then raise Found) arr ; false
   with Found -> true
 
-type fou_idx_intermediate =
+type four_idx_intermediate =
 {
   expo_b : float array;
   expo_d : float array;
@@ -577,7 +577,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
 
 
 
-let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q : float Zmap.t =
+let contracted_class_shell_pairs ~zero_m ?mu_erf ?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
@@ -651,16 +651,13 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
                       let norm_pq_sq =
                         Co.dot center_pq center_pq
                       in
-                      let normalization = Psp.normalization sp.(i-1) *.
-                        Psp.normalization sq.(i-1)
-                      in
 
-                      let zero = Zp.zero zero_m in
+                      let zero = Zp.zero ?mu_erf zero_m in
                       let zero_m_array = 
                         zero_m
                         {zero with
                           expo_p_inv ; expo_q_inv ; norm_pq_sq ;
-                          center_pq ; center_pa ; center_qc ; normalization ;
+                          center_pq ; center_pa ; center_qc ; 
                         }
                       in
                       zero_m_array.(0)
@@ -786,15 +783,11 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
                         let norm_pq_sq =
                           x *. x +. y *. y +. z *. z
                         in
-                        let normalization = Psp.normalization shell_ab *.
+                        let zero = Zp.zero ?mu_erf zero_m in
-                          Psp.normalization shell_cd
-                        in
-                        let zero = Zp.zero 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} ;
                                 center_pa ; center_qc = center_qc_tmp.(cd) ;
-                                normalization ;
                               }
                     ) sq
                 in

Basis/Zero_m_parameters.ml

@@ -3,25 +3,27 @@ type t =
   expo_p_inv  : float ;
   expo_q_inv  : float ;
   norm_pq_sq  : float ;
-  normalization : float ;
   maxm        : int   ;
   center_pq   : Coordinate.t ;
   center_pa   : Coordinate.t ;
   center_qc   : Coordinate.t ;
   f12_factor  : F12factor.t option;
+  mu_erf      : float option;
   zero_m_func : t -> float array ;
 }
 
-let zero ?f12_factor zero_m_func = 
+let zero ?f12_factor ?mu_erf zero_m_func = 
 {
   zero_m_func ;
   f12_factor ;
+  mu_erf;
   maxm=0 ; 
   expo_p_inv = 0.;
   expo_q_inv = 0.;
   norm_pq_sq = 0.;
-  normalization = 0.;
   center_pq = Coordinate.zero ;
   center_pa = Coordinate.zero ;
   center_qc = Coordinate.zero ;
 }
+
+

CI/CI.ml

@@ -97,10 +97,12 @@ let create_matrix_arbitrary f det_space =
           in
           let singles =
             List.filter (fun (i,d,det_j) -> d < 2) doubles
-            |> List.map (fun (i,_,det_j) -> (i,det_j))
+            |> List.rev_map (fun (i,_,det_j) -> (i,det_j))
+            |> List.rev
           in
           let doubles =
-            List.map (fun (i,_,det_j) -> (i,det_j)) doubles
+            List.rev_map (fun (i,_,det_j) -> (i,det_j)) doubles
+            |> List.rev
           in
           (singles, doubles)
         ) det_beta
@@ -262,10 +264,12 @@ let create_matrix_spin ?(nmax=2) f det_space =
                     in
                     let singles =
                       List.filter (fun (i,d,det_j) -> d < 2) doubles
-                      |> List.map (fun (i,_,det_j) -> (i,det_j))
+                      |> List.rev_map (fun (i,_,det_j) -> (i,det_j))
+                      |> List.rev
                     in
                     let doubles =
-                      List.map (fun (i,_,det_j) -> (i,det_j)) doubles
+                      List.rev_map (fun (i,_,det_j) -> (i,det_j)) doubles
+                      |> List.rev
                     in
                     (singles, doubles, triples)
                   ) b
@@ -292,13 +296,16 @@ let create_matrix_spin ?(nmax=2) f det_space =
                     in
 
                     let triples =
-                      List.map (fun (i,_,det_j) -> (i,det_j)) triples
+                      List.rev_map (fun (i,_,det_j) -> (i,det_j)) triples
+                      |> List.rev
                     in
                     let doubles =
-                      List.map (fun (i,_,det_j) -> (i,det_j)) doubles
+                      List.rev_map (fun (i,_,det_j) -> (i,det_j)) doubles
+                      |> List.rev
                     in
                     let singles =
-                      List.map (fun (i,_,det_j) -> (i,det_j)) singles
+                      List.rev_map (fun (i,_,det_j) -> (i,det_j)) singles
+                      |> List.rev
                     in
                     (singles, doubles, triples)
                   ) b
@@ -769,7 +776,8 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
     in
 
     let psi_filtered =
-      List.map (fun i -> psi0.(i)) psi_filtered_idx
+      List.rev_map (fun i -> psi0.(i)) psi_filtered_idx
+      |> List.rev
     in
 
     let psi_h_alfa alfa =
@@ -896,39 +904,34 @@ let second_order_sum2 { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
 
     Ds.determinants_array det_space
     |> Array.to_list
-    |> List.map (fun det_i ->
+    |> List.concat_map (fun det_i ->
         [ Spin.Alfa ; Spin.Beta ]
-        |>  List.map (fun spin ->
+        |>  List.concat_map (fun spin ->
-            List.map (fun particle ->
+            List.concat_map (fun particle ->
               List.map (fun hole ->
                   [ [ Determinant.single_excitation spin hole particle det_i ] ;
-                    List.map (fun particle' ->
+                    List.concat_map (fun particle' ->
                         List.map (fun hole' ->
                           Determinant.double_excitation
                             spin  hole particle
                             spin  hole' particle' det_i
                         ) list_holes
                     ) list_particles
-                    |> List.concat
                     ;
-                    List.map (fun particle' ->
+                    List.concat_map (fun particle' ->
                         List.map (fun hole' ->
                           Determinant.double_excitation
                             spin  hole particle
                             (Spin.other spin)  hole' particle' det_i
                         ) list_holes
                     ) list_particles
-                    |> List.concat
                   ]
                   |> List.concat
                 ) list_holes
             ) list_particles
-            |> List.concat
           )
-          |> List.concat
         )
         |> List.concat
-        |> List.concat
     |> List.filter (fun alfa -> not (Determinant.is_none alfa))
     |> List.sort_uniq compare
   in
@@ -988,16 +991,16 @@ let is_external det_space =
         Determinant.alfa a
         |> Spindeterminant.bitstring
       in
-      let n_a = 
+      let n_a =
-	Bitstring.(popcount @@ logand inactive_mask alfa) 
+	Bitstring.(popcount @@ logand inactive_mask alfa)
       in
       match n_a with
-      | 0 | 1 | 2 -> 
+      | 0 | 1 | 2 ->
 	let beta =
 	  Determinant.beta a
 	  |> Spindeterminant.bitstring
 	in
-	n_a + Bitstring.(popcount @@ logand inactive_mask beta) < 3 
+	n_a + Bitstring.(popcount @@ logand inactive_mask beta) < 3
       | _ ->   false
 
 

CI/DeterminantSpace.ml

@@ -323,10 +323,11 @@ let fci_f12_of_mo_basis mo_basis ~frozen_core mo_num =
   in
   { r with mo_class =
     MOClass.to_list r.mo_class
-    |> List.map (fun i ->
+    |> List.rev_map (fun i ->
       match i with
       | MOClass.Virtual i when i > mo_num -> MOClass.Auxiliary i 
       | i -> i)
+    |> List.rev
     |> MOClass.of_list }
     
 
@@ -339,10 +340,11 @@ let cas_f12_of_mo_basis mo_basis ~frozen_core n m mo_num =
   in
   { r with mo_class =
     MOClass.to_list r.mo_class
-    |> List.map (fun i ->
+    |> List.rev_map (fun i ->
       match i with
       | MOClass.Virtual i when i > mo_num -> MOClass.Auxiliary i 
       | i -> i)
+    |> List.rev
     |> MOClass.of_list
     }
 

CI/Excitation.ml

@@ -76,7 +76,7 @@ let multiple_of_spindet t t' =
     else
       Phase.Neg
   in
-  (phase, List.map2 (fun hole particle -> (hole, particle)) holes (List.rev particles) )
+  (phase, List.rev @@ List.rev_map2 (fun hole particle -> (hole, particle)) holes (List.rev particles) )
 
 
 let double_of_spindet t t' =
@@ -99,8 +99,8 @@ let multiple_of_det t t' =
   in
   let phase = Phase.add pa pb in
   Multiple (phase, List.concat [
-    List.map (fun (hole, particle) -> { hole ; particle ; spin=Spin.Alfa }) a ;
+    List.rev @@ List.rev_map (fun (hole, particle) -> { hole ; particle ; spin=Spin.Alfa }) a ;
-    List.map (fun (hole, particle) -> { hole ; particle ; spin=Spin.Beta }) b ])
+    List.rev @@ List.rev_map (fun (hole, particle) -> { hole ; particle ; spin=Spin.Beta }) b ])
   
 
 let double_of_det t t' =

CI/Spindeterminant.ml

@@ -93,7 +93,8 @@ let holes_particles_of t t' =
   let holes = Bitstring.logand (bitstring t) x |> Bitstring.to_list 
   and particles = Bitstring.logand (bitstring t') x |> Bitstring.to_list 
   in
-  List.map2 (fun h p -> (h,p)) holes particles
+  List.rev_map2 (fun h p -> (h,p)) holes particles
+  |> List.rev
 
 
 let set_phase p = function

CI/SpindeterminantSpace.ml

@@ -30,8 +30,8 @@ let fci_of_mo_basis ~frozen_core mo_basis elec_num =
   let spin_determinants = 
     Bitstring.permtutations elec_num mo_num 
     |> List.filter (fun b -> Bitstring.logand neg_active_mask b = occ_mask)
-    |> List.map (fun b -> Spindeterminant.of_bitstring b)
     |> Array.of_list
+    |> Array.map (fun b -> Spindeterminant.of_bitstring b)
   in
   { elec_num ; mo_basis ; mo_class ; spin_determinants }
 
@@ -54,8 +54,8 @@ let cas_of_mo_basis mo_basis ~frozen_core elec_num n m =
   let spin_determinants = 
     Bitstring.permtutations elec_num mo_num 
     |> List.filter (fun b -> Bitstring.logand neg_active_mask b = occ_mask)
-    |> List.map (fun b -> Spindeterminant.of_bitstring b)
     |> Array.of_list
+    |> Array.map (fun b -> Spindeterminant.of_bitstring b)
   in
   { elec_num ; mo_basis ; mo_class ; spin_determinants }
 

INSTALL.md

@@ -1,7 +1,7 @@
 # Generic installation
 
 ```bash
-opam install lacaml mpi getopt alcotest zarith
+opam install ocamlbuild lacaml mpi getopt alcotest zarith
 ```
 
 

MOBasis/HF12.ml

@@ -67,10 +67,11 @@ let array_4_init d1 d2 d3 d4 fx =
       SharedMemory.create Bigarray.Float64 [| d1;d2;d3;d4 |] 
   in
   Util.list_range 1 d4
-  |> List.map (fun l -> 
+  |> List.rev_map (fun l -> 
     Util.list_range 1 d3
-    |> List.map (fun k -> (k,l)) )
+    |> List.rev_map (fun k -> (k,l)) )
   |> List.concat
+  |> List.rev
   |> Stream.of_list 
   |> Farm.run ~f ~ordered:false
   |> Stream.iter (fun (k,l,x) ->
@@ -133,10 +134,11 @@ let array_5_init d1 d2 d3 d4 d5 fx =
       SharedMemory.create Bigarray.Float64 [| d1;d2;d3;d4;d5 |]
   in
   Util.list_range 1 d5
-  |> List.map (fun m -> 
+  |> List.rev_map (fun m -> 
     Util.list_range 1 d4
-    |> List.map (fun l -> (l,m)) )
+    |> List.rev_map (fun l -> (l,m)) )
   |> List.concat
+  |> List.rev
   |> Stream.of_list 
   |> Farm.run ~f ~ordered:false
   |> Stream.iter (fun (l,m,x) ->

MOBasis/MOClass.ml

@@ -34,51 +34,45 @@ let to_list t = t
 
 
 let core_mos t =
-  List.map (fun x ->
+  List.filter_map (fun x ->
       match x with 
       | Core i -> Some i
       | _ -> None) t
-  |> Util.list_some
 
 
 let inactive_mos t =
-  List.map (fun x ->
+  List.filter_map (fun x ->
       match x with 
       | Inactive i -> Some i
       | _ -> None ) t
-  |> Util.list_some
 
 
 let active_mos t =
-  List.map (fun x ->
+  List.filter_map (fun x ->
       match x with 
       | Active i -> Some i
       | _ -> None ) t
-  |> Util.list_some
 
 
 let virtual_mos t =
-  List.map (fun x ->
+  List.filter_map (fun x ->
       match x with 
       | Virtual i -> Some i
       | _ -> None ) t
-  |> Util.list_some
 
 
 let deleted_mos t =
-  List.map (fun x ->
+  List.filter_map (fun x ->
       match x with 
       | Deleted i -> Some i
       | _ -> None ) t
-  |> Util.list_some
 
 
 let auxiliary_mos t =
-  List.map (fun x ->
+  List.filter_map (fun x ->
       match x with 
       | Auxiliary i -> Some i
       | _ -> None ) t
-  |> Util.list_some
 
 
 let mo_class_array t =

Notebooks/F12_matrix.ipynb

@@ -9724,7 +9724,7 @@
  "metadata": {
   "celltoolbar": "Raw Cell Format",
   "kernelspec": {
-   "display_name": "OCaml 4.07.1",
+   "display_name": "OCaml default",
    "language": "OCaml",
    "name": "ocaml-jupyter"
   },

Notebooks/SpVec.ipynb

@@ -1746,7 +1746,7 @@
  "metadata": {
   "celltoolbar": "Raw Cell Format",
   "kernelspec": {
-   "display_name": "OCaml 4.07.1",
+   "display_name": "OCaml default",
    "language": "OCaml",
    "name": "ocaml-jupyter"
   },

Notebooks/TemplateQCaml.ipynb

@@ -0,0 +1,159 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "val png_image : string -> unit = <fun>\n"
+      ]
+     },
+     "execution_count": 4,
+     "metadata": {},
+     "output_type": "execute_result"
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "- : unit = ()\n",
+      "Findlib has been successfully loaded. Additional directives:\n",
+      "  #require \"package\";;      to load a package\n",
+      "  #list;;                   to list the available packages\n",
+      "  #camlp4o;;                to load camlp4 (standard syntax)\n",
+      "  #camlp4r;;                to load camlp4 (revised syntax)\n",
+      "  #predicates \"p,q,...\";;   to set these predicates\n",
+      "  Topfind.reset();;         to force that packages will be reloaded\n",
+      "  #thread;;                 to enable threads\n",
+      "\n",
+      "- : unit = ()\n"
+     ]
+    },
+    {
+     "data": {
+      "text/plain": [
+       "val png_image : string -> Jupyter_notebook.display_id = <fun>\n"
+      ]
+     },
+     "execution_count": 4,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "let png_image = print_endline ;;\n",
+    "\n",
+    "(* --------- *)\n",
+    "\n",
+    "#cd \"/home/scemama/QCaml\";;\n",
+    "#use \"topfind\";;\n",
+    "#require \"jupyter.notebook\";;\n",
+    "\n",
+    "let png_image name = \n",
+    "    Jupyter_notebook.display_file ~base64:true \"image/png\" (\"Notebooks/images/\"^name)\n",
+    ";;\n",
+    "\n",
+    "#require \"lacaml.top\";;\n",
+    "#require \"alcotest\";;\n",
+    "#require \"str\";;\n",
+    "#require \"bigarray\";;\n",
+    "#require \"zarith\";;\n",
+    "#require \"getopt\";;\n",
+    "#directory \"_build\";;\n",
+    "#directory \"_build/Basis\";;\n",
+    "#directory \"_build/CI\";;\n",
+    "#directory \"_build/MOBasis\";;\n",
+    "#directory \"_build/Nuclei\";;\n",
+    "#directory \"_build/Parallel\";;\n",
+    "#directory \"_build/Perturbation\";;\n",
+    "#directory \"_build/SCF\";;\n",
+    "#directory \"_build/Utils\";;\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "File Constants.cmo is not a bytecode object file.\n",
+      "File Util.cma is not a bytecode object file.\n",
+      "File Matrix.cmo is not a bytecode object file.\n",
+      "File Simulation.cmo is not a bytecode object file.\n",
+      "File Spindeterminant.cmo is not a bytecode object file.\n",
+      "File Determinant.cmo is not a bytecode object file.\n",
+      "File HartreeFock.cmo is not a bytecode object file.\n",
+      "File MOBasis.cmo is not a bytecode object file.\n",
+      "File F12CI.cmo is not a bytecode object file.\n"
+     ]
+    }
+   ],
+   "source": [
+    "#load \"Constants.cmo\";;\n",
+    "#load_rec \"Util.cma\";;\n",
+    "#load_rec \"Matrix.cmo\";;\n",
+    "#load_rec \"Simulation.cmo\";;\n",
+    "#load_rec \"Spindeterminant.cmo\";;\n",
+    "#load_rec \"Determinant.cmo\";;\n",
+    "#load_rec \"HartreeFock.cmo\";;\n",
+    "#load_rec \"MOBasis.cmo\";;\n",
+    "#load_rec \"F12CI.cmo\";;"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#install_printer AngularMomentum.pp_string ;;\n",
+    "#install_printer Basis.pp ;;\n",
+    "#install_printer Charge.pp ;;\n",
+    "#install_printer Coordinate.pp ;;\n",
+    "#install_printer Vector.pp;;\n",
+    "#install_printer Matrix.pp;;\n",
+    "#install_printer Util.pp_float_2darray;;\n",
+    "#install_printer Util.pp_float_array;;\n",
+    "#install_printer Util.pp_matrix;;\n",
+    "#install_printer HartreeFock.pp ;;\n",
+    "#install_printer Fock.pp ;;\n",
+    "#install_printer MOClass.pp ;;\n",
+    "let pp_mo ppf t = MOBasis.pp ~start:1 ~finish:0 ppf t ;;\n",
+    "#install_printer pp_mo;;\n",
+    "(*\n",
+    "#install_printer DeterminantSpace.pp;;\n",
+    "*)\n",
+    "#install_printer SpindeterminantSpace.pp;;\n",
+    "#install_printer Bitstring.pp;;\n",
+    "\n",
+    "(* --------- *)\n",
+    "\n",
+    "open Lacaml.D\n"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "OCaml default",
+   "language": "OCaml",
+   "name": "ocaml-jupyter"
+  },
+  "language_info": {
+   "codemirror_mode": "text/x-ocaml",
+   "file_extension": ".ml",
+   "mimetype": "text/x-ocaml",
+   "name": "OCaml",
+   "nbconverter_exporter": null,
+   "pygments_lexer": "OCaml",
+   "version": "4.07.1"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}

SCF/Guess.ml

@@ -66,7 +66,6 @@ let test_case ao_basis =
         Lacaml.D.Mat.add 
           (AOBasis.eN_ints  ao_basis |> NucInt.matrix)
           (AOBasis.kin_ints ao_basis |> KinInt.matrix)
-        |> Conventions.rephase
         |> Lacaml.D.Mat.to_array
       in
       Array.iteri (fun i x ->

Simulation.ml

@@ -5,6 +5,7 @@ type t = {
   basis             : Basis.t;
   ao_basis          : AOBasis.t;
   f12               : F12factor.t option;
+  mu_erf            : float option;
   nuclear_repulsion : float;
 }
 
@@ -14,12 +15,14 @@ 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 mu_erf            t = t.mu_erf
 let f12               t = t.f12
 
 let make ?cartesian:(cartesian=false)
          ?multiplicity:(multiplicity=1)
          ?charge:(charge=0)
          ?f12
+         ?mu_erf
          ~nuclei
          basis
   =
@@ -45,18 +48,19 @@ let make ?cartesian:(cartesian=false)
   in
 
   {
-    charge ; basis ; nuclei ; electrons ; ao_basis ; f12 ;
+    charge ; basis ; nuclei ; electrons ; ao_basis ; f12 ; mu_erf ;
     nuclear_repulsion ; 
   }
 
 
-let of_filenames ?(cartesian=false) ?(multiplicity=1) ?(charge=0) ?f12 ~nuclei ?(aux_basis_filenames=[]) basis_filename =
+let of_filenames ?(cartesian=false) ?(multiplicity=1) ?(charge=0) ?f12 ?mu_erf ~nuclei ?(aux_basis_filenames=[]) basis_filename =
   let nuclei =
     Nuclei.of_filename nuclei
   in
   let basis =
     Basis.of_nuclei_and_basis_filenames ~nuclei (basis_filename :: aux_basis_filenames)
   in
-  lazy (make ?f12 ~cartesian ~charge ~multiplicity ~nuclei basis)
+
+  lazy (make ?mu_erf ?f12 ~cartesian ~charge ~multiplicity ~nuclei basis)
   |> Parallel.broadcast
 

Simulation.mli

@@ -21,15 +21,18 @@ val nuclear_repulsion : t -> float
 val f12 : t -> F12factor.t option
 (** f12 correlation factor *)
 
+val mu_erf : t -> float option
+(** Range-separation parameter of the electron repulsion integrals potential *)
+
 (** {1 Creation} *)
 
 val make :
   ?cartesian:bool ->
   ?multiplicity:int -> ?charge:int -> ?f12:F12factor.t -> 
-  nuclei:Nuclei.t -> Basis.t -> t
+  ?mu_erf:float -> nuclei:Nuclei.t -> Basis.t -> t
 
 val of_filenames :
   ?cartesian:bool ->
   ?multiplicity:int -> ?charge:int -> ?f12:F12factor.t -> 
-  nuclei:string -> ?aux_basis_filenames:string list -> string -> t
+  ?mu_erf: float -> nuclei:string -> ?aux_basis_filenames:string list -> string -> t
 

Utils/AngularMomentum.ml

@@ -107,19 +107,19 @@ let zkey_array a =
       begin
         match a with
         | Singlet l1 -> 
-          List.map (fun x -> Zkey.of_powers_three x) (keys_1d @@ to_int l1)
+          List.rev_map (fun x -> Zkey.of_powers_three x) (keys_1d @@ to_int l1)
 
         | Doublet (l1, l2) -> 
-          List.map (fun a ->
+          List.rev_map (fun a ->
-              List.map (fun b -> Zkey.of_powers_six a b) (keys_1d @@ to_int l2)
+              List.rev_map (fun b -> Zkey.of_powers_six a b) (keys_1d @@ to_int l2)
             ) (keys_1d @@ to_int l1)
           |> List.concat
 
         | Triplet (l1, l2, l3) ->
 
-          List.map (fun a ->
+          List.rev_map (fun a ->
-              List.map (fun b ->
+              List.rev_map (fun b ->
-                  List.map (fun c ->
+                  List.rev_map (fun c ->
                       Zkey.of_powers_nine a b c) (keys_1d @@ to_int l3)
                 ) (keys_1d @@ to_int l2)
               |> List.concat
@@ -128,10 +128,10 @@ let zkey_array a =
 
         | Quartet (l1, l2, l3, l4) ->
 
-          List.map (fun a ->
+          List.rev_map (fun a ->
-              List.map (fun b ->
+              List.rev_map (fun b ->
-                  List.map (fun c ->
+                  List.rev_map (fun c ->
-                      List.map (fun d ->
+                      List.rev_map (fun d ->
                           Zkey.of_powers_twelve a b c d) (keys_1d @@ to_int l4)
                     ) (keys_1d @@ to_int l3)
                   |> List.concat
@@ -140,6 +140,7 @@ let zkey_array a =
             ) (keys_1d @@ to_int l1)
           |> List.concat
       end
+      |> List.rev
       |> Array.of_list
     in
     Hashtbl.add zkey_array_memo a result;

Utils/Davidson.ml

@@ -117,7 +117,7 @@ let make
     in
 
 
-    let residual_norms = List.map nrm2 u_proposed in
+    let residual_norms = List.rev @@ List.rev_map nrm2 u_proposed in
     let residual_norm  =
       List.fold_left (fun accu i -> accu +. i *. i) 0. residual_norms
       |> sqrt

Utils/Matrix.ml

@@ -80,11 +80,12 @@ let sparse_of_computed ?(threshold=epsilon) = function
   | Computed {m ; n ; f} ->
     Sparse { m ; n ; v=Array.init n (fun j ->
         Util.list_range 1 m
-        |> List.map (fun i ->
+        |> List.rev_map (fun i ->
             let x = f i (j+1) in
             if abs_float x > threshold then Some (i, x)
             else None)
         |> Util.list_some
+        |> List.rev
         |> Vector.sparse_of_assoc_list m
       ) }
   | _ -> invalid_arg "Expected a computed matrix"
@@ -176,7 +177,7 @@ let outer_product ?(threshold=epsilon) v1 v2 =
     in
     let v = 
       Array.init (Vector.dim v2) (fun j -> 
-          List.map (fun (i, x) -> 
+          List.rev_map (fun (i, x) -> 
               let z = x *. v'.{j+1} in
               if abs_float z < threshold then
                 None
@@ -184,6 +185,7 @@ let outer_product ?(threshold=epsilon) v1 v2 =
                 Some (i, z)
             ) v
           |> Util.list_some
+          |> List.rev
           |> Vector.sparse_of_assoc_list (Vector.dim v1)
         )
     in
@@ -500,22 +502,23 @@ let split_cols nrows = function
       Mat.to_col_vecs a
       |> Array.to_list
       |> Util.list_pack nrows
-      |> List.map (fun l ->
+      |> List.rev_map (fun l ->
             Dense (Mat.of_col_vecs @@ Array.of_list l) )
+      |> List.rev
     end
   | Sparse a -> 
     begin
       Array.to_list a.v 
       |> Util.list_pack nrows
-      |> List.map Array.of_list
+      |> List.rev_map Array.of_list
-      |> List.map (fun v -> Sparse { m=a.m ; n= Array.length v ; v })
+      |> List.rev_map (fun v -> Sparse { m=a.m ; n= Array.length v ; v })
     end
   | Computed a -> 
     begin
       Util.list_range 0 (a.n-1)
       |> Util.list_pack nrows
-      |> List.map Array.of_list
+      |> List.rev_map Array.of_list
-      |> List.map (fun v -> Computed { m=a.m ; n= Array.length v ; f = (fun i j -> a.f i (j+v.(0)) ) })
+      |> List.rev_map (fun v -> Computed { m=a.m ; n= Array.length v ; f = (fun i j -> a.f i (j+v.(0)) ) })
     end
   
 
@@ -534,7 +537,9 @@ let join_cols l =
     | [] -> Sparse { m=0 ; n=0 ; v=[| |] }
     | (Dense    a) :: rest -> aux_dense  [] ((Dense  a) :: rest)
     | (Sparse   a) :: rest -> aux_sparse 0 0 [] ((Sparse a) :: rest)
-    | (Computed a) :: rest -> aux_sparse 0 0 [] (List.map sparse_of_computed ( (Computed a) :: rest ))
+    | (Computed a) :: rest -> aux_sparse 0 0 []
+        (List.rev_map sparse_of_computed ( (Computed a) :: rest )
+         |> List.rev)
 
   in aux (List.rev l)
 

Utils/Util.ml

@@ -133,34 +133,37 @@ let fact = function
   | i -> fact_memo.(i)
 
 
-let binom_func n k =
-  (*TODO : slow function *)
-  assert (n >= k);
-  let rec aux n k =
-    if k = 0 || k = n then
-      1
-    else
-      aux (n-1) (k-1) + aux (n-1) k
-  in aux n k
-
 let binom =
   let memo = 
-    Array.init 64 (fun n -> Array.init 64 (fun k ->
+    let m = 
-      if n >= k then binom_func n k else 0) )
+      Array.make_matrix 64 64 0 
+    in
+    for n=0 to Array.length m - 1 do
+      m.(n).(0) <- 1;
+      m.(n).(n) <- 1;
+      for k=1 to (n - 1) do
+        m.(n).(k) <- m.(n-1).(k-1) + m.(n-1).(k)
+      done
+    done;
+    m
   in
-  fun n k ->
+  let rec f n k =
-    if n < 64 then 
+    assert (k >= 0);
-      begin
+    assert (n >= k);
-        assert (n >= k);
+    if k = 0 || k = n then
+      1
+    else if n < 64 then 
         memo.(n).(k)
-      end
     else
-      binom_func n k
+      f (n-1) (k-1) + f (n-1) k
+  in f
+
 
 let binom_float n k =
   binom n k
   |> float_of_int_fast
 
+
 let rec pow a = function
   | 0 -> 1.
   | 1 -> a
@@ -206,29 +209,26 @@ let boys_function ~maxm t =
       let result =
         Array.init (maxm+1) (fun m -> 1. /. float_of_int (2*m+1))
       in
-      (*
+      let power_t_inv = (maxm+maxm+1)  in
-      assert (abs_float t > 1.e-10);
+      try  
-      *)
-      if t <> 0. then
-        begin
           let fmax = 
             let t_inv = sqrt (1. /. t) in
             let n = float_of_int maxm in
             let dm = 0.5 +. n in
-            let f  = (pow t_inv (maxm+maxm+1) ) in
+            let f  = (pow t_inv power_t_inv ) in
             match classify_float f with
               | FP_normal -> (incomplete_gamma dm t) *. 0.5 *. f
               | FP_zero
               | FP_subnormal -> 0.
-              | _ -> invalid_arg "zero_m overflow"
+              | _ -> raise Exit
           in
           let emt = exp (-. t) in
           result.(maxm) <- fmax;
           for n=maxm-1 downto 0 do
             result.(n) <- ( (t+.t) *. result.(n+1) +. emt) *. result.(n)
-          done
+          done;
-        end;
+          result
-      result
+      with Exit -> result
     end
 
 
@@ -241,7 +241,8 @@ let of_some = function
 
 let list_some l =
   List.filter (function None -> false | _ -> true) l
-  |> List.map (function Some x -> x | _ -> assert false)
+  |> List.rev_map (function Some x -> x | _ -> assert false)
+  |> List.rev
 
 
 let list_range first last = 

Utils/Vector.ml

@@ -91,8 +91,8 @@ let sparse_of_vec ?(threshold=epsilon) v =
 
 let sparse_of_assoc_list n v =
   Sparse { n ;
-    v = List.map (fun (index, value) -> {index ; value}) v
+    v = Array.of_list v
-        |> Array.of_list
+        |> Array.map (fun (index, value) -> {index ; value})
   }
 
 

run_integrals.ml

@@ -1,6 +1,7 @@
 let out_file    : string option ref = ref None
 let basis_file  : string option ref = ref None
 let nuclei_file : string option ref = ref None
+let mu_erf      : float  option ref = ref None
 
 
 let speclist = [
@@ -8,6 +9,8 @@ let speclist = [
     "File containing the atomic basis set") ;
   ( "-x" , Arg.String (fun x -> nuclei_file := Some x),
     "File containing the nuclear coordinates") ;
+  ( "-m" , Arg.Float  (fun x -> mu_erf := Some x),
+    "Value of mu, the range separation factor.") ;
 ]
 
 let run () =
@@ -19,10 +22,11 @@ let run () =
     match !nuclei_file with
     | None -> raise (Invalid_argument "Coordinate file should be specified with -x")
     | Some x -> x
+  and mu_erf = !mu_erf
   in
 
   let s =
-    Simulation.of_filenames ~nuclei:nuclei_file basis_file 
+    Simulation.of_filenames ?mu_erf ~nuclei:nuclei_file basis_file 
   in
 
   print_endline @@ Nuclei.to_string @@ Simulation.nuclei s;
@@ -38,8 +42,11 @@ let run () =
   Overlap.to_file ~filename:("Ov.dat") overlap;
   NucInt.to_file ~filename:("Nuc.dat") eN_ints;
   KinInt.to_file ~filename:("Kin.dat") kin_ints;
-  ERI.to_file ~filename:("ERI.dat") ee_ints;
+  ERI.to_file    ~filename:("ERI.dat") ee_ints;
-  ()
+  match mu_erf with
+  | Some mu_erf ->
+      ERI_lr.to_file ~filename:("ERI_lr.dat") (AOBasis.ee_lr_ints ao_basis)
+  | None -> ()
 
 
 let () =

run_tests.ml

@@ -1,32 +1,54 @@
   
 let () =
 
+
+  (* H2O / cc-pVDZ *)
+  let test_integrals_water, test_guess_water = 
+    let basis_file   = "test_files/cc-pvdz"
+    and nuclei_file  = "test_files/water.xyz"
+    in
+
+    let simulation_closed_shell =
+      Simulation.of_filenames ~charge:0 ~multiplicity:1 ~nuclei:nuclei_file basis_file
+    in
+
+    let ao_basis = 
+      Simulation.ao_basis simulation_closed_shell
+    in
+    (fun () -> AOBasis.test_case "water" ao_basis),
+    (fun () -> Guess.test_case ao_basis);
+  in
+
+
+  (* LiH / cc-pVTZ *)
+  let test_integrals_lih =
+    let basis_file   = "test_files/cc-pvdz"
+    and nuclei_file  = "test_files/lih.xyz"
+    in
+
+    let simulation_closed_shell =
+      Simulation.of_filenames ~mu_erf:0.5 ~charge:0 ~cartesian:true ~multiplicity:1 ~nuclei:nuclei_file basis_file
+    in
+
+    let ao_basis = 
+      Simulation.ao_basis simulation_closed_shell
+    in
+    (fun () -> AOBasis.test_case "lih" ao_basis)
+  in
+
   Alcotest.run "Unit tests" [
-    "Util", Util.test_case ();
+      "Util", Util.test_case ();
-    "Bitstring", Bitstring.test_case ();
+      "Bitstring", Bitstring.test_case ();
-    "Spindeterminant", Spindeterminant.test_case ();
+      "Spindeterminant", Spindeterminant.test_case ();
-    "Determinant", Determinant.test_case ();
+      "Determinant", Determinant.test_case ();
-    "Excitation", Excitation.test_case ();
+      "Excitation", Excitation.test_case ();
-    "Sparse vectors", Vector.test_case ();
+      "Sparse vectors", Vector.test_case ();
-    "Sparse matrices", Matrix.test_case ();
+      "Sparse matrices", Matrix.test_case ();
-    "Cholesky", Cholesky.test_case ();
+      "Cholesky", Cholesky.test_case ();
-  ];
+      "Guess water", test_guess_water ();
+      "AO Integrals water", test_integrals_water ();
+      "AO Integrals LiH", test_integrals_lih ();
+    ]
 
 
-  let basis_file   = "test_files/cc-pvdz"
-  and nuclei_file  = "test_files/h2o.xyz"
-  in
-
-  let simulation_closed_shell =
-    Simulation.of_filenames ~charge:0 ~multiplicity:1 ~nuclei:nuclei_file basis_file
-  in
-
-  let ao_basis = 
-    Simulation.ao_basis simulation_closed_shell
-  in
-  Alcotest.run "Water, cc-pVDZ" [
-    "AO_Basis", AOBasis.test_case ao_basis;
-    "Guess", Guess.test_case ao_basis;
-  ]
-
 

test_files/lih.xyz

@@ -0,0 +1,5 @@
+2
+LiH
+Li       0.000000      0.000000      0.000000
+H        0.000000      0.000000      1.639920
+

test_files/lih_kin.ref

@@ -0,0 +1,69 @@
+   1    1   3.611833164750e+00
+   1    2  -7.043972403440e-01
+   2    2   1.714942506008e-01
+   1    3   1.334243764489e-02
+   2    3   2.523613418040e-02
+   3    3   4.207500000000e-02
+   4    4   8.984286842706e-02
+   5    5   8.984286842706e-02
+   6    6   8.984286842706e-02
+   4    7   3.465216246844e-02
+   7    7   6.007500000000e-02
+   5    8   3.465216246844e-02
+   8    8   6.007500000000e-02
+   6    9   3.465216246844e-02
+   9    9   6.007500000000e-02
+   1   10  -7.615906158701e-02
+   2   10   5.555021308733e-02
+   3   10   3.749351798583e-02
+  10   10   2.684500000000e-01
+  11   11   4.336500000000e-01
+  12   12   4.336500000000e-01
+   1   13  -7.615906158701e-02
+   2   13   5.555021308733e-02
+   3   13   3.749351798583e-02
+  10   13  -2.065000000000e-02
+  13   13   2.684500000000e-01
+  14   14   4.336500000000e-01
+   1   15  -7.615906158701e-02
+   2   15   5.555021308733e-02
+   3   15   3.749351798583e-02
+  10   15  -2.065000000000e-02
+  13   15  -2.065000000000e-02
+  15   15   2.684500000000e-01
+   1   16  -1.642492953508e-02
+   2   16   1.617265204482e-02
+   3   16   8.786521484803e-03
+   6   16   3.703320786321e-02
+   9   16   1.186053439203e-02
+  10   16  -2.459048298758e-02
+  13   16  -2.459048298758e-02
+  15   16   1.143468816632e-01
+  16   16   3.347671997727e-01
+   1   17   1.504210795985e-02
+   2   17   2.915477195499e-02
+   3   17   3.231574311466e-02
+   6   17   7.703233575788e-02
+   9   17   3.915758653648e-02
+  10   17  -2.973530349695e-03
+  13   17  -2.973530349695e-03
+  15   17   1.311020407537e-01
+  16   17   1.182520739417e-01
+  17   17   1.830000000000e-01
+   4   18   1.982822516835e-02
+   7   18   6.226173253833e-03
+  12   18   1.482007009374e-01
+  18   18   1.817500000000e+00
+   5   19   1.982822516835e-02
+   8   19   6.226173253833e-03
+  14   19   1.482007009374e-01
+  19   19   1.817500000000e+00
+   1   20   4.721148005029e-02
+   2   20  -2.863428424260e-02
+   3   20  -4.264822214538e-03
+   6   20  -4.834234856015e-02
+   9   20   2.222472845464e-03
+  10   20   2.574462749493e-02
+  13   20   2.574462749493e-02
+  15   20  -4.716114478488e-02
+  20   20   1.817500000000e+00

test_files/lih_nuc.ref

@@ -0,0 +1,87 @@
+     1     1   -8.381828186445269
+     1     2    1.179001296158989
+     2     2   -0.549381198537747
+     1     3   -0.691991161919778
+     2     3   -0.455903587224674
+     3     3   -1.027906943724010
+     4     4   -0.465989084986239
+     5     5   -0.465989084986239
+     1     6   -0.009311304353181
+     2     6   -0.037430525449801
+     3     6   -0.059079445909559
+     6     6   -0.500424104643515
+     4     7   -0.359176298697607
+     7     7   -0.645815634254628
+     5     8   -0.359176298697607
+     8     8   -0.645815634254628
+     1     9   -0.002980968207840
+     2     9   -0.029486033852596
+     3     9   -0.061347174703540
+     6     9   -0.386436952594028
+     9     9   -0.678822272157784
+     1    10   -0.147278589160974
+     2    10   -0.441648276259208
+     3    10   -0.713237808727169
+     6    10   -0.022721765213106
+     9    10   -0.019544779518061
+    10    10   -1.125890026980850
+    11    11   -1.125890026980850
+     4    12   -0.039355251786751
+     7    12   -0.033852551148014
+    12    12   -1.172194150656758
+     1    13   -0.147278589160974
+     2    13   -0.441648276259208
+     3    13   -0.713237808727169
+     6    13   -0.022721765213106
+     9    13   -0.019544779518061
+    10    13   -0.375296675660283
+    13    13   -1.125890026980850
+     5    14   -0.039355251786751
+     8    14   -0.033852551148014
+    14    14   -1.172194150656758
+     1    15   -0.151198126172260
+     2    15   -0.482700544116671
+     3    15   -0.774599673502485
+     6    15   -0.114136808458971
+     9    15   -0.094187507956906
+    10    15   -0.390731383552253
+    13    15   -0.390731383552253
+    15    15   -1.290917179075705
+     1    16   -0.070723878218244
+     2    16   -0.169030847542230
+     3    16   -0.249803229886425
+     6    16   -0.296525237835227
+     9    16   -0.201280582486981
+    10    16   -0.058887783283701
+    13    16   -0.058887783283701
+    15    16   -0.606125224490307
+    16    16   -0.810069851165359
+     1    17   -0.650096934511009
+     2    17   -0.396410154179119
+     3    17   -0.785967649140876
+     6    17   -0.569284102460052
+     9    17   -0.497476954555701
+    10    17   -0.431045748470551
+    13    17   -0.431045748470551
+    15    17   -1.044365316472985
+    16    17   -0.771046294920814
+    17    17   -1.496005797103620
+     4    18   -0.139242826009133
+     7    18   -0.092708609257056
+    12    18   -0.467812962335489
+    18    18   -1.840472494075094
+     5    19   -0.139242826009133
+     8    19   -0.092708609257056
+    14    19   -0.467812962335489
+    19    19   -1.840472494075094
+     1    20    0.138618867744086
+     2    20    0.115004494251478
+     3    20    0.126444897306323
+     6    20    0.189479103637283
+     9    20   -0.005399573459954
+    10    20    0.055409942558440
+    13    20    0.055409942558440
+    15    20    0.183195153579426
+    16    20    0.115478769718145
+    17    20    0.184168640350032
+    20    20   -1.944456654208564

test_files/lih_overlap.ref

@@ -0,0 +1,69 @@
+   1    1   1.001041982979e+00
+   1    2  -8.270896233604e-02
+   2    2   2.844647710602e-01
+   1    3   1.656106292542e-01
+   2    3   4.322930475687e-01
+   3    3   1.000000000000e+00
+   4    4   3.594196016605e-01
+   5    5   3.594196016605e-01
+   6    6   3.594196016605e-01
+   4    7   3.756686814419e-01
+   7    7   1.000000000000e+00
+   5    8   3.756686814419e-01
+   8    8   1.000000000000e+00
+   6    9   3.756686814419e-01
+   9    9   1.000000000000e+00
+   1   10   6.111087257716e-02
+   2   10   3.738979197853e-01
+   3   10   6.435570260814e-01
+  10   10   1.000000000000e+00
+  11   11   1.000000000000e+00
+  12   12   1.000000000000e+00
+   1   13   6.111087257716e-02
+   2   13   3.738979197853e-01
+   3   13   6.435570260814e-01
+  10   13   3.333333333333e-01
+  13   13   1.000000000000e+00
+  14   14   1.000000000000e+00
+   1   15   6.111087257716e-02
+   2   15   3.738979197853e-01
+   3   15   6.435570260814e-01
+  10   15   3.333333333333e-01
+  13   15   3.333333333333e-01
+  15   15   1.000000000000e+00
+   1   16   2.122870259016e-02
+   2   16   8.933604469184e-02
+   3   16   1.330876041146e-01
+   6   16   1.517284987815e-01
+   9   16   1.136230292020e-01
+  10   16   3.661059004591e-02
+  13   16   3.661059004591e-02
+  15   16   3.113618431571e-01
+  16   16   3.453140483256e-01
+   1   17   1.518845348249e-01
+   2   17   3.029195011034e-01
+   3   17   5.530798525107e-01
+   6   17   3.784156759728e-01
+   9   17   4.226976097304e-01
+  10   17   3.223839042442e-01
+  13   17   3.223839042442e-01
+  15   17   6.971907703330e-01
+  16   17   4.024473615409e-01
+  17   17   1.000000000000e+00
+   4   18   8.243364452249e-02
+   7   18   5.878583263933e-02
+  12   18   2.818737540320e-01
+  18   18   1.000000000000e+00
+   5   19   8.243364452249e-02
+   8   19   5.878583263933e-02
+  14   19   2.818737540320e-01
+  19   19   1.000000000000e+00
+   1   20  -4.507855171808e-02
+   2   20  -4.679063013530e-02
+   3   20  -3.523797537160e-02
+   6   20  -6.404549505991e-02
+   9   20   3.252087980480e-02
+  10   20  -2.773517482041e-02
+  13   20  -2.773517482041e-02
+  15   20  -3.315348957381e-02
+  20   20   1.000000000000e+00