Fixed 4-idx

MOBasis/MOBasis.ml

@@ -1,5 +1,6 @@
 open Lacaml.D
 open Util
+open Constants
 
 (** One-electron orthogonal basis set, corresponding to Molecular Orbitals. *)
 
@@ -33,8 +34,7 @@ type t =
 
 
 let mo_matrix_of_ao_matrix ~mo_coef ao_matrix = 
-  gemm ~transa:`T mo_coef   @@
-  gemm ao_matrix mo_coef  
+  xt_o_x ~x:mo_coef ~o:ao_matrix
 
 
 let ao_matrix_of_mo_matrix ~mo_coef ~ao_overlap mo_matrix =
@@ -56,14 +56,10 @@ let four_index_transform ~mo_coef eri_ao =
   let range_mo = list_range ~start:1 mo_num in
   let range_ao = list_range ~start:1 ao_num in
 
-  let u = 
-      Mat.create mo_num_2 mo_num
-  and o = 
-      Mat.create ao_num ao_num_2 
-  and p = 
-      Mat.create ao_num_2 mo_num
-  and q = 
-      Mat.create ao_mo_num mo_num
+  let u = Mat.create mo_num_2 mo_num
+  and o = Mat.create ao_num ao_num_2 
+  and p = Mat.create ao_num_2 mo_num
+  and q = Mat.create ao_mo_num mo_num
   in
   Printf.eprintf "Transforming %d integrals : %!" mo_num;
   List.iter (fun delta ->
@@ -71,7 +67,8 @@ let four_index_transform ~mo_coef eri_ao =
       Mat.fill u 0.;
 
       List.iter (fun l ->
-
+          if abs_float mo_coef.{l,delta} > epsilon then
+          begin
             let jk = ref 0 in
             List.iter (fun k ->
                 List.iter (fun j ->
@@ -101,9 +98,9 @@ let four_index_transform ~mo_coef eri_ao =
             in
 
             ignore @@
-            gemm ~transa:`T ~beta:1. ~c:u q' mo_coef 
+              gemm ~transa:`T ~beta:1. ~alpha:mo_coef.{l,delta} ~c:u q' mo_coef ;
               (* u_alphabeta_gamma = \sum_k q_k_alphabeta c_k_gamma *)
-
+          end
       ) range_ao;
       let u =
         Bigarray.reshape 
@@ -172,7 +169,35 @@ let of_rhf ~frozen_core hf =
       |> Vec.of_array
   in
   let mo_coef = hf.HF.eigenvectors in
+  let result = 
     make ~ao_basis ~mo_type ~mo_class ~mo_occupation ~mo_coef ()
+  in
+
+  let () =
+    let e = ref 0. in
+    let t = KinInt.matrix (Lazy.force result.kin_ints) in
+    let v = NucInt.matrix (Lazy.force result.eN_ints) in
+    let g = Lazy.force result.ee_ints in
+    for i = 1 to 5 do
+      e := !e +. 2. *. (t.{i,i} +. v.{i,i})
+    done;
+    Printf.printf "Energy Mono = %20.15f\n" !e;
+    let e2 = ref 0. in
+    for i = 1 to 5 do
+      for j = i+1 to 5 do
+        e2 := !e2 +. 2. *. (ERI.get_phys g i j i j -.
+                            ERI.get_phys g i j j i)
+      done;
+    done;
+    for i = 1 to 5 do
+      for j = 1 to 5 do
+        e2 := !e2 +. ERI.get_phys g i j i j 
+      done;
+    done;
+    Printf.printf "Energy bi = %20.15f\n" !e2;
+    Printf.printf "Energy = %20.15f\n" (simulation.Si.nuclear_repulsion +. !e +. !e2)
+  in
+  result
 
 
 let of_hartree_fock ~frozen_core = function

Utils/Util.ml

@@ -172,7 +172,7 @@ let list_some l =
 let stream_range ?(start=0) n = 
   Stream.from (fun i ->
     let result = i+start in
-    if result < n then
+    if result <= n then
       Some result
     else None
   )

Utils/Util.mli

@@ -64,11 +64,11 @@ val list_some : 'a option list -> 'a list
     the [Some]. *)
 
 val list_range : ?start:int -> int -> int list
-(** Returns a list [start ; start+1 ; ... ; start+(n-1)]. Default is [start=0]. *)
+(** Returns a list [start ; start+1 ; ... ; start+n]. Default is [start=0]. *)
 
 (** {2 Useful streams} *)
 val stream_range : ?start:int -> int -> int Stream.t
-(** Returns a stream <start ; start+1 ; ... ; start+(n-1)>. Default is [start=0]. *)
+(** Returns a stream <start ; start+1 ; ... ; start+n>. Default is [start=0]. *)
 
 
 (** {2 Linear algebra } *)