Preparing tests for Erf integrals

Basis/AOBasis.ml

@@ -9,6 +9,7 @@ type t =
   eN_ints            :  NucInt.t              lazy_t;
   kin_ints           :  KinInt.t              lazy_t;
   ee_ints            :  ERI.t                 lazy_t;
+  ee_lr_ints         :  ERI_lr.t              lazy_t;
   f12_ints           :  F12.t                 lazy_t;
   f12_over_r12_ints  :  ScreenedERI.t         lazy_t;
   cartesian          :  bool                        ;
@@ -20,6 +21,7 @@ let ortho     t = Lazy.force t.ortho
 let eN_ints      t = Lazy.force t.eN_ints
 let kin_ints     t = Lazy.force t.kin_ints
 let ee_ints      t = Lazy.force t.ee_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
@@ -48,6 +50,10 @@ 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
   ) in
@@ -56,13 +62,13 @@ let make ~cartesian ~basis ?f12 nuclei =
     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
@@ -86,7 +92,7 @@ let test_case t =
 
   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 +102,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 +112,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,12 +122,24 @@ 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
     in
     check_eri "ee_ints" ee_ints reference
+    ;
+  in
+
+  let test_ee_lr_ints () =
+    let reference =
+      ERI.of_file ("test_files/"^name^"_eri_lr.ref") ~sparsity:`Dense
+          ~size:(Basis.size t.basis)
+    in
+    let ee_ints =
+      Lazy.force t.ee_ints
+    in
+    check_eri "ee_lr_ints" ee_ints reference
   in
 
   [
@@ -129,6 +147,7 @@ let test_case t =
     "eN_ints",     `Quick,  test_eN_ints;
     "kin_ints",    `Quick,  test_kin_ints;
     "ee_ints",     `Quick,  test_ee_ints;
+    "ee_lr_ints",  `Quick,  test_ee_lr_ints;
   ]
 
 

Basis/AOBasis.mli

@@ -19,6 +19,9 @@ val eN_ints : t -> NucInt.t
 val ee_ints : t -> ERI.t
 (** Electron-electron potential integrals *)
 
+val ee_lr_ints : t -> ERI_lr.t
+(** Electron-electron long-range potential integrals *)
+
 val f12_ints : t -> F12.t
 (** Electron-electron potential integrals *)
 
@@ -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/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/Zero_m_parameters.ml

@@ -9,13 +9,15 @@ type 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.;

INSTALL.md

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

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"
   },

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,18 @@ 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/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
-  fun n k ->
+    for n=0 to Array.length m - 1 do
-    if n < 64 then 
+      m.(n).(0) <- 1;
-      begin
+      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
+  let rec f n k =
+    assert (k >= 0);
     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

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;
@@ -39,7 +43,10 @@ let run () =
   NucInt.to_file ~filename:("Nuc.dat") eN_ints;
   KinInt.to_file ~filename:("Kin.dat") kin_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,6 +1,41 @@
   
 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-pvtz"
+    and nuclei_file  = "test_files/lih.xyz"
+    in
+
+    let simulation_closed_shell =
+      Simulation.of_filenames ~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 ();
       "Bitstring", Bitstring.test_case ();
@@ -10,23 +45,10 @@ let () =
       "Sparse vectors", Vector.test_case ();
       "Sparse matrices", Matrix.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;
     ]
 
 
+