Inserted overlap

Basis/Basis.ml

@@ -45,4 +45,21 @@ let to_string b =
     )
     ^ line
 
+let file : string option ref = ref None
+
+let set_file f =
+  file := Some f
+
+let general_basis = lazy(
+  match !file with
+  | None -> failwith "basis set file not defined"
+  | Some filename -> Gamess_reader.read ~filename
+)   
+
+
+let basis = lazy (
+  of_nuclei_and_general_basis
+    (Lazy.force Nuclei.nuclei) (Lazy.force general_basis)
+)
+
 

Basis/Basis.mli

@@ -1,3 +1,7 @@
 type t = Contracted_shell.t array
 val of_nuclei_and_general_basis : Nuclei.t -> General_basis.t list -> t
 val to_string : t -> string
+val set_file : string -> unit
+val general_basis :
+  (Element.t * General_basis.general_contracted_shell array) list lazy_t
+val basis : Contracted_shell.t array lazy_t

Basis/OneElectronRR.ml

@@ -129,27 +129,22 @@ let contracted_class_nuc ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t
   in
   ()
 
-(* TODO
   (* Compute all integrals in the shell for each pair of significant shell pairs *)
 
   for ab=0 to (Array.length shell_p - 1)
   do
-    let b = shell_p.(ab).Shell_pair.j in
-
-    for cd=0 to (Array.length shell_q - 1)
-    do
       let coef_prod =
-         shell_p.(ab).Shell_pair.coef *. shell_q.(cd).Shell_pair.coef
+         shell_p.(ab).Shell_pair.coef 
       in
       (** Screening on thr product of coefficients *)
       if (abs_float coef_prod) > 1.e-4*.cutoff then
         begin
 
             let expo_pq_inv =
-              shell_p.(ab).Shell_pair.expo_inv +. shell_q.(cd).Shell_pair.expo_inv
+              shell_p.(ab).Shell_pair.expo_inv 
             in
-            let center_pq = 
-              Coordinate.(shell_p.(ab).Shell_pair.center |- shell_q.(cd).Shell_pair.center)
+            let center_ab = 
+              Coordinate.shell_p.(ab).Shell_pair.center_ab 
             in
             let norm_pq_sq =
               Coordinate.dot center_pq center_pq

Basis/Overlap.ml

@@ -0,0 +1,111 @@
+open Util
+
+(** Computes all the overlap integrals of the contracted shell pair *)
+let contracted_class shell_a shell_b : float Zmap.t =
+
+  let shell_p =
+    Shell_pair.create_array shell_a shell_b
+  in
+
+  (* Pre-computation of integral class indices *)
+  let class_indices = 
+    Angular_momentum.zkey_array (Angular_momentum.Doublet
+         Contracted_shell.(totAngMom shell_a, totAngMom shell_b))
+  in
+
+  let contracted_class = 
+    Array.make (Array.length class_indices) 0.
+  in
+
+  (* Compute all integrals in the shell for each pair of significant shell pairs *)
+
+  for ab=0 to (Array.length shell_p - 1)
+  do
+      let coef_prod =
+         shell_p.(ab).Shell_pair.coef 
+      in
+      (** Screening on thr product of coefficients *)
+      if (abs_float coef_prod) > 1.e-4*.cutoff then
+        begin
+            let center_ab = 
+              shell_p.(ab).Shell_pair.center_ab 
+            in
+            let center_a = 
+              shell_p.(ab).Shell_pair.center_a
+            in
+            let expo_inv = 
+              shell_p.(ab).Shell_pair.expo_inv
+            in
+
+            Array.iteri (fun i key ->
+                      let (angMomA,angMomB) =
+                        let a = Zkey.to_int_array Zkey.Kind_6 key in
+                        ( [| a.(0) ;  a.(1) ;  a.(2) |], 
+                          [| a.(3) ;  a.(4) ;  a.(5) |] )
+                      in
+                      let norm = 
+                        shell_p.(ab).Shell_pair.norm_fun angMomA angMomB 
+                      in
+                      let f k = 
+                        Overlap_primitives.hvrr (angMomA.(k), angMomB.(k)) 
+                        expo_inv
+                        (Coordinate.coord center_ab k,
+                         Coordinate.coord center_a  k) 
+                      in
+                      let integral = chop norm (fun () -> (f 0)*.(f 1)*.(f 2)) in 
+                      contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral
+                ) class_indices
+        end
+  done;
+  let result = 
+    Zmap.create (Array.length contracted_class)
+  in
+  Array.iteri (fun i key -> Zmap.add result key contracted_class.(i)) class_indices;
+  result
+
+
+
+
+
+(** Write all overlap integrals to a file *)
+let to_file ~filename basis =
+  let to_int_tuple x = 
+    let open Zkey in
+    match to_int_tuple Kind_3 x with
+    | Three x -> x
+    | _ -> assert false
+  in
+
+  let oc = open_out filename in
+  for i=0 to (Array.length basis) - 1 do
+    print_int basis.(i).Contracted_shell.indice ; print_newline ();
+    for j=0 to i do
+      (* Compute all the integrals of the class *)
+      let cls =
+        contracted_class basis.(i) basis.(j) 
+      in
+
+      (* Write the data in the output file *)
+      Array.iteri (fun i_c powers_i ->
+          let i_c = basis.(i).Contracted_shell.indice + i_c + 1 in
+          let xi = to_int_tuple powers_i in
+          Array.iteri (fun j_c powers_j ->
+              let j_c = basis.(j).Contracted_shell.indice + j_c + 1 in
+              let xj = to_int_tuple powers_j in
+              let key = 
+                Zkey.of_int_tuple (Zkey.Six (xi,xj))
+              in
+              let value = 
+                try
+                  Zmap.find cls key 
+                with Not_found -> failwith "Bug in overlap integrals"
+              in
+              if (abs_float value > cutoff) then
+                Printf.fprintf oc "%4d %4d %20.12e\n"
+                  i_c j_c value
+            ) basis.(j).Contracted_shell.powers
+        ) basis.(i).Contracted_shell.powers;
+    done;
+  done;
+  close_out oc
+

Basis/Overlap_primitives.ml

@@ -0,0 +1,34 @@
+open Util
+
+
+(** Horizontal and Vertical Recurrence Relations (HVRR) *)
+let hvrr (angMom_a, angMom_b) (expo_inv_p) (center_ab, center_pa)
+  =
+
+  (** Vertical recurrence relations *)
+  let rec vrr angMom_a =
+
+    if angMom_a < 0 then 0.
+    else if angMom_a = 0 then 1.
+    else
+      chop center_pa (fun () -> vrr (angMom_a-1))
+        +. chop (0.5 *. (float_of_int (angMom_a-1)) *. expo_inv_p) 
+        (fun () -> vrr (angMom_a-2))
+
+
+  (** Horizontal recurrence relations *)
+  and hrr angMom_a angMom_b =
+
+    if angMom_b < 0 then 0.
+    else if angMom_b = 0 then vrr angMom_a 
+    else
+      hrr (angMom_a+1) (angMom_b-1) +. chop center_ab
+      (fun () -> hrr angMom_a (angMom_b-1) )
+
+  in
+  hrr angMom_a angMom_b 
+
+
+
+let overlap = hvrr
+

Basis/Shell_pair.ml

@@ -4,6 +4,7 @@ type t = {
   expo     : float;
   expo_inv : float;
   center_ab: Coordinate.t;
+  center_a : Coordinate.t;
   center   : Coordinate.t;
   norm_sq  : float;
   norm     : float;
@@ -38,6 +39,12 @@ let create_array ?cutoff p_a p_b =
 
       Array.init (Contracted_shell.size p_b) (fun j ->
           try
+            if Contracted_shell.center p_a = Contracted_shell.center p_b &&
+              ((Angular_momentum.to_int @@ Contracted_shell.totAngMom p_a) +
+               (Angular_momentum.to_int @@ Contracted_shell.totAngMom p_b)) land 1 = 1
+            then
+              raise Null_contribution;
+
             let f2 =
               Contracted_shell.norm_coef p_b j
             in
@@ -76,7 +83,10 @@ let create_array ?cutoff p_a p_b =
             in
             if (abs_float coef  < cutoff) then
               raise Null_contribution;
-            Some { i ; j ; norm_fun ; norm ; coef ; expo ; expo_inv ; center ; center_ab ; norm_sq }
+            let center_a =
+              Coordinate.(center |- Contracted_shell.center p_a)
+            in
+            Some { i ; j ; norm_fun ; norm ; coef ; expo ; expo_inv ; center ; center_a ; center_ab ; norm_sq }
           with
           | Null_contribution -> None
         )
@@ -88,3 +98,5 @@ let create_array ?cutoff p_a p_b =
   |> List.map (function Some x -> x | None -> assert false)
   |> Array.of_list
 
+open Util
+

Basis/TwoElectronRR.ml

@@ -1,10 +1,5 @@
 open Util
 
-(** In chop f g, evaluate g only if f is non zero, and return f *. (g ()) *)
-let chop f g =
-  if (abs_float f) < cutoff then 0.
-  else f *. (g ())
-
 
 (** Horizontal and Vertical Recurrence Relations (HVRR) *)
 let hvrr_two_e m (angMom_a, angMom_b, angMom_c, angMom_d)

Nuclei/Nuclei.ml

@@ -49,3 +49,14 @@ let to_string atoms =
 "
 
 
+let file : string option ref = ref None
+
+let set_file f =
+  file := Some f
+
+let nuclei = lazy(
+  match !file with
+  | None -> failwith "coordinate file not defined"
+  | Some filename -> of_xyz_file ~filename
+)   
+

Utils/Angular_momentum.ml

@@ -99,8 +99,8 @@ let zkey_array a =
     | Doublet (l1, l2) -> 
       List.map (fun a ->
           List.map (fun b ->
-              Zkey.of_int_tuple (Zkey.Six (a,b))) (keys_1d @@ to_int l1)
-        ) (keys_1d @@ to_int l2)
+              Zkey.of_int_tuple (Zkey.Six (a,b))) (keys_1d @@ to_int l2)
+        ) (keys_1d @@ to_int l1)
       |> List.concat
 
     | Triplet (l1, l2, l3) ->

Utils/Util.ml

@@ -69,3 +69,8 @@ let rec pow a = function
 ;;
 
 
+(** In chop f g, evaluate g only if f is non zero, and return f *. (g ()) *)
+let chop f g =
+  if (abs_float f) < cutoff then 0.
+  else f *. (g ())
+

Utils/Zkey.ml

@@ -10,18 +10,18 @@ type kind_array =
 | Kind_2
 | Kind_1
 
+let (<|) x a = 
+  Z.logor (Z.shift_left x 64) a
+
+let (<<) x a = 
+  Int64.logor (Int64.shift_left x 10) (Int64.of_int a)
+
+let (<+) x a = 
+  Int64.logor (Int64.shift_left x 16) (Int64.of_int a)
+
 
 (** Build a Zkey from an array or 1, 2, 3, 4, 6, 9, or 12 integers *)
 let of_int_array ~kind a = 
-  let (<|) x a = 
-    Z.logor (Z.shift_left x 64) a
-  in
-  let (<<) x a = 
-    Int64.logor (Int64.shift_left x 10) (Int64.of_int a)
-  in
-  let (<+) x a = 
-    Int64.logor (Int64.shift_left x 16) (Int64.of_int a)
-  in
   match kind with
   | Kind_3 -> (Int64.of_int a.(0)) << a.(1) << a.(2) |> Z.of_int64
   | Kind_6 -> (Int64.of_int a.(0)) << a.(1) << a.(2) << a.(3) << a.(4) << a.(5)
@@ -56,15 +56,6 @@ type kind =
   | Nine   of ((int*int*int)*(int*int*int)*(int*int*int))
   | Twelve of ((int*int*int)*(int*int*int)*(int*int*int)*(int*int*int))
 
-let (<|) x a = 
-  Z.logor (Z.shift_left x 64) a
-
-let (<<) x a = 
-  Int64.logor (Int64.shift_left x 10) (Int64.of_int a)
-
-let (<+) x a = 
-  Int64.logor (Int64.shift_left x 16) (Int64.of_int a)
-
 let of_int_tuple a = 
   match a with
   | One   (a)       -> Z.of_int a
@@ -143,16 +134,8 @@ let to_int_array ~kind a =
   | Kind_1  -> [| Z.to_int a |]
 
 
-let to_string ~kind a =
-  "< " ^ ( Z.to_string a ) ^ " | " ^ (
-  to_int_array kind a 
-  |> Array.map string_of_int
-  |> Array.to_list
-  |> String.concat ", "
-  ) ^ " >"
 
-
-(** Transform the Zkey into an int array *)
+(** Transform the Zkey into an int tuple *)
 let to_int_tuple ~kind a = 
   match kind with
   | Kind_3  -> Three  ( Z.to_int @@ Z.extract a  20 10 ,
@@ -202,6 +185,14 @@ let to_int_tuple ~kind a =
 
 include Z
 
+let to_string ~kind a =
+  "< " ^ ( Z.to_string a ) ^ " | " ^ (
+  to_int_array kind a 
+  |> Array.map string_of_int
+  |> Array.to_list
+  |> String.concat ", "
+  ) ^ " >"
+
 (*
 let debug () = 
   let k2  = of_int_array Kind_2  [| 1 ; 2 |]

run_integrals.ml

@@ -1,60 +1,36 @@
-let basis_file : string option ref = ref None
-let coord_file : string option ref = ref None
 let out_file   : string option ref = ref None
 
 
 let speclist = [
-  ( "-b" , Arg.String (fun x -> basis_file := Some x) ,
+  ( "-b" , Arg.String Basis.set_file ,
     "File containing the atomic basis set") ;
-  ( "-c" , Arg.String (fun x -> coord_file := Some x) ,
+  ( "-c" , Arg.String Nuclei.set_file ,
     "File containing the nuclear coordinates") ;
   ( "-o" , Arg.String (fun x -> out_file := Some x) ,
     "Output file") ;
 ]
 
-let run ~coord ~basis ~out =
-  let coord_file = 
-    match coord with 
-    | None -> raise (Invalid_argument "Coordinate file should be specified with -c")
-    | Some x -> x
-  and basis_file = 
-    match basis with 
-    | None -> raise (Invalid_argument "Basis set file should be specified with -b")
-    | Some x -> x
-  and out_file = 
+let run ~out =
+  let out_file = 
     match out with 
     | None -> raise (Invalid_argument "Output file should be specified with -o")
     | Some x -> x
   in
 
-  let nuclei = 
-    Nuclei.of_xyz_file ~filename:coord_file
-  in
+  let nuclei = Lazy.force Nuclei.nuclei in
   print_endline @@ Nuclei.to_string nuclei;
 
-  let basis =
-    let general_basis =
-      Gamess_reader.read ~filename:basis_file
-    in
-    Basis.of_nuclei_and_general_basis nuclei general_basis
-  in
+  let basis = Lazy.force Basis.basis in
   print_endline @@ Basis.to_string basis;
 
-  ERI.to_file ~filename:out_file basis
-
+  ERI.to_file ~filename:(out_file^".eri") basis
+(*
+  Overlap.to_file ~filename:(out_file^".overlap") basis
+  *)
 
 
 let () =
   let usage_msg = "Available options:" in
   Arg.parse speclist (fun _ -> ()) usage_msg;
-    run ~coord:!coord_file ~basis:!basis_file ~out:!out_file
-  (*
-  try
-  with
-  | Invalid_argument e ->
-    begin
-      print_string "Error: " ; print_endline e; print_newline ();
-      Arg.usage speclist usage_msg
-    end
-    *)
+  run ~out:!out_file