Cleaned code

CI/F12CI.ml

@@ -24,7 +24,11 @@ let eigensystem t =  Lazy.force t.eigensystem
 let hf_ij_non_zero hf12_integrals deg_a deg_b ki kj =
   let integrals = [
     let one_e _ _ _ = 0. in
-    let hf12_s, hf12_o, hf12m_s, hf12m_o = hf12_integrals in
+    let { HF12.
+          simulation ; aux_basis ;
+          hf12 ; hf12_anti ;
+          hf12_single ; hf12_single_anti } = hf12_integrals
+    in
     let kia = De.alfa ki and kib = De.beta ki
     and kja = De.alfa kj and kjb = De.beta kj
     in
@@ -37,14 +41,14 @@ let hf_ij_non_zero hf12_integrals deg_a deg_b ki kj =
     in
     let two_e i j k l s s' =
       if s' = Spin.other s then
-	hf12_o.{i,j,k,l}  -. 1. *. (
+	hf12.{i,j,k,l}  -. 1. *. (
-	(List.fold_left (fun accu m -> accu +. hf12m_o.{m,i,j,k,l}) 0. mo_a) +.
+	(List.fold_left (fun accu m -> accu +. hf12_single.{m,i,j,k,l}) 0. mo_a) +.
-	(List.fold_left (fun accu m -> accu +. hf12m_o.{m,j,i,l,k}) 0. mo_b) 
+	(List.fold_left (fun accu m -> accu +. hf12_single.{m,j,i,l,k}) 0. mo_b) 
         )
       else
-	hf12_s.{i,j,k,l} -. 1. *. (
+	hf12_anti.{i,j,k,l} -. 1. *. (
-	(List.fold_left (fun accu m -> accu +. hf12m_s.{m,i,j,k,l}) 0. mo_a) +.
+	(List.fold_left (fun accu m -> accu +. hf12_single_anti.{m,i,j,k,l}) 0. mo_a) +.
-	(List.fold_left (fun accu m -> accu +. hf12m_s.{m,j,i,l,k}) 0. mo_b) 
+	(List.fold_left (fun accu m -> accu +. hf12_single_anti.{m,j,i,l,k}) 0. mo_b) 
         )
     in
     (one_e, two_e)
@@ -128,8 +132,6 @@ let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filen
 
       let eigenvectors, eigenvalues = 
 
-(* Column dressing 
-*)
         let delta = lacpy delta in
         Mat.scal f delta;
         for k=1 to state-1 do
@@ -166,65 +168,10 @@ let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filen
         in
 
 
-(* Diagonal dressing *)
-(*
-        let diagonal =
-          Vec.init (Matrix.dim1 m_H) (fun i ->
-            Matrix.get m_H i i +.
-              if (abs_float psi.{i,state} > 1.e-8) then
-                delta.{i,state} /. psi.{i,state}
-              else  0.
-            )
-        in
-
-        let matrix_prod c = 
-          let w = 
-            Matrix.mm ~transa:`T m_H c
-            |> Matrix.to_mat
-          in
-          for i=1 to (Mat.dim1 w) do
-            w.{i,state} <- w.{i,state} +. delta.{i,state}
-          done;
-          Matrix.dense_of_mat w
-        in
-*)
-
-
         Parallel.broadcast (lazy (
           Davidson.make ~threshold:1.e-9 ~guess:psi ~n_states:state diagonal matrix_prod
           ))
 
-(*
-        let m_H = Matrix.to_mat m_H |> lacpy in
-*)
-
-(* DIAGONAL TEST 
-        for i=1 to Mat.dim1 m_H do
-          if (abs_float psi.{i,state} > 1.e-8) then
-            m_H.{i,i} <- m_H.{i,i} +. delta.{i,state} /. psi.{i,state};
-        done;
-*)
-
-
-(* COLUMN TEST
-        for i=1 to Mat.dim1 m_H do
-          let d = delta.{i,state} /. psi.{column_idx,state} in
-          m_H.{i,column_idx} <- m_H.{i,column_idx} +. d;
-          if (i <> column_idx) then
-            begin
-              m_H.{column_idx,i} <- m_H.{column_idx,i} +. d;
-              m_H.{column_idx,column_idx} <- m_H.{column_idx,column_idx} -.
-                d *. psi.{i,state}
-            end
-        done;
-*)
-
-
-
-(*
-        let m_v = syev m_H in
-        m_H, m_v
-*)
       in
 
 
@@ -240,6 +187,14 @@ let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filen
         Printf.printf "F12 Convergence : %e  %f\n" conv (eigenvalues.{state} 
           +. Simulation.nuclear_repulsion simulation);
 
+(*
+      let cabs_singles = 
+        let f =
+          Fock.make_rhf ~density ~ao_basis:large_ao_basis
+        in
+      in
+*)
+
       if conv > threshold then
         iteration ~state eigenvectors
       else

MOBasis/HF12.ml

@@ -4,10 +4,25 @@ open Lacaml.D
 
 module Fis = FourIdxStorage
 
-type t = (float, Bigarray.float64_elt, Bigarray.fortran_layout) Bigarray.Genarray.t 
+
-       * (float, Bigarray.float64_elt, Bigarray.fortran_layout) Bigarray.Genarray.t 
+type q = (float, Bigarray.float64_elt, Bigarray.fortran_layout) Bigarray.Genarray.t 
-       * (float, Bigarray.float64_elt, Bigarray.fortran_layout) Bigarray.Genarray.t 
+
-       * (float, Bigarray.float64_elt, Bigarray.fortran_layout) Bigarray.Genarray.t 
+type t = {
+  simulation : Simulation.t ;
+  aux_basis  : MOBasis.t ;
+  hf12     : q ; (* hf12.{i,j,k,l} = \sum_{ab} <ij|H|ab> <ab|F12|kl>
+                                   + \sum_{am} <ij|H|am> <am|F12|kl>
+                                   + \sum_{bm} <ij|H|mb> <mb|F12|kl>
+                 *)
+  hf12_anti: q ; (* hf12_anti.{i,j,k,l} = 
+                      \sum_{ab} (<ij|H|ab> - <ij|H|ba>) (<ab|F12|kl> - <ab|F12|lk>)
+                    + \sum_{am} (<ij|H|am> - <ij|H|ma>) (<am|F12|kl> - <am|F12|lk>)
+                    + \sum_{bm} (<ij|H|mb> - <ij|H|bm>) (<mb|F12|kl> - <mb|F12|lk>)
+                    *)
+  hf12_single     : q ; (* hf12.{m,i,j,k,l} = \sum_{a} <ij|H|am> <am|F12|kl> *)
+  hf12_single_anti: q ; (* hf12_anti.{m,i,j,k,l} = 
+                           \sum_{ab} (<ij|H|am> - <ij|H|ma>) (<am|F12|kl> - <am|F12|lk>) *)
+}
 
 
 let make ~simulation ~mo_basis ~aux_basis_filename () =
@@ -16,8 +31,7 @@ let make ~simulation ~mo_basis ~aux_basis_filename () =
   let mo_num = MOBasis.size mo_basis in
 
   (* Add auxiliary basis set *)                                                               
-  let aux_basis = 
+  let simulation =
-    let s =
     let charge        = Charge.to_int @@ Simulation.charge simulation 
     and multiplicity  = Electrons.multiplicity @@ Simulation.electrons simulation
     and nuclei        = Simulation.nuclei simulation
@@ -31,7 +45,9 @@ let make ~simulation ~mo_basis ~aux_basis_filename () =
     |> Basis.of_nuclei_and_general_basis nuclei
     |> Simulation.make ~f12 ~charge ~multiplicity ~nuclei 
   in
-    MOBasis.of_mo_basis s mo_basis
+
+  let aux_basis = 
+    MOBasis.of_mo_basis simulation mo_basis
   in
 
   let aux_num = MOBasis.size aux_basis in
@@ -47,11 +63,17 @@ let make ~simulation ~mo_basis ~aux_basis_filename () =
   (* Compute the <ij|QHF|kl> integrals *)
   if Parallel.master then Printf.eprintf "Computing HF12 integrals\n%!";
 
-  let result_s, result_o, resultm_s, resultm_o =
+  let hf12, hf12_anti, hf12_single, hf12_single_anti =
-    Bigarray.Genarray.create Float64 Bigarray.fortran_layout [| mo_num ; mo_num ; mo_num ; mo_num |] , 
+    let create4 n =
-    Bigarray.Genarray.create Float64 Bigarray.fortran_layout [| mo_num ; mo_num ; mo_num ; mo_num |] , 
+      Bigarray.Genarray.create Float64 Bigarray.fortran_layout [| n ; n ; n ; n |] 
-    Bigarray.Genarray.create Float64 Bigarray.fortran_layout [| mo_num ; mo_num ; mo_num ; mo_num ; mo_num |] , 
+    in
-    Bigarray.Genarray.create Float64 Bigarray.fortran_layout [| mo_num ; mo_num ; mo_num ; mo_num ; mo_num |] 
+    let create5 n =
+      Bigarray.Genarray.create Float64 Bigarray.fortran_layout [| n ; n ; n ; n ; n |] 
+    in
+    create4 mo_num ,
+    create4 mo_num ,
+    create5 mo_num ,
+    create5 mo_num 
   in
 
 
@@ -200,11 +222,11 @@ let make ~simulation ~mo_basis ~aux_basis_filename () =
   |> Stream.iter (fun (hf_s, hf_o, hfm_s, hfm_o, j, l) ->
     for k=1 to mo_num do
       for i=1 to mo_num do
-        result_s.{i,j,k,l} <- hf_s.{i,k} ;
+        hf12.{i,j,k,l}      <- hf_o.{i,k} ;
-        result_o.{i,j,k,l} <- hf_o.{i,k} ;
+        hf12_anti.{i,j,k,l} <- hf_s.{i,k} ;
         for m=1 to mo_num do
-          resultm_s.{m,i,j,k,l} <- hfm_s.{m,i,k} ;
+          hf12_single.{m,i,j,k,l}      <- hfm_o.{m,i,k} ;
-          resultm_o.{m,i,j,k,l} <- hfm_o.{m,i,k}
+          hf12_single_anti.{m,i,j,k,l} <- hfm_s.{m,i,k}
         done
       done
     done );
@@ -224,7 +246,13 @@ let make ~simulation ~mo_basis ~aux_basis_filename () =
   Printf.printf "%!";
   *)
 
-  Parallel.broadcast (lazy (result_s, result_o, resultm_s, resultm_o) )
+  let result = 
+    { simulation ; aux_basis ;
+      hf12 ; hf12_anti ;
+      hf12_single ; hf12_single_anti
+    } 
+  in
+  Parallel.broadcast (lazy result)