Fixed f12 integrals

Basis/F12.ml

@@ -1,4 +1,342 @@
-(** Two-electron integrals over Slater geminals via a fit using Gaussian geminals.
+(** Two electron integral functor for operators that are separable among %{ $(x,y,z)$ %}.
+    It is parameterized by the [zero_m] function.
 *)
 
-include TwoElectronIntegralsSeparable
+open Constants
+let cutoff = integrals_cutoff
+
+module Bs  = Basis
+module Cs  = ContractedShell
+module Csp = ContractedShellPair
+module Cspc = ContractedShellPairCouple
+module Fis = FourIdxStorage
+
+
+include FourIdxStorage
+
+(** Exponent of the geminal *)
+let expo_s = 1.0
+
+(** Coefficients and exponents of the Gaussian fit of the Slater Geminal*) 
+let coef_g =
+    [| 0.3144 ; 0.3037 ; 0.1681 ; 0.09811 ; 0.06024 ; 0.03726 |]
+        
+let expo_sg_inv =
+  Array.map (fun x -> 1. /. (x *. expo_s *. expo_s))
+    [| 0.2209 ; 1.004 ;  3.622 ; 12.16 ;   45.87 ;  254.4 |]
+
+
+
+(*
+
+Fit of 1/r:
+
+let coef_g = [|
+  841.88478132     ;  70.590185207     ;  18.3616020768    ;  7.2608642093    ;
+  3.57483416444    ;  2.01376031082    ;  1.24216542801    ;  0.81754348620   ;
+  0.564546514023   ;  0.404228610699   ;  0.297458536575   ;  0.223321219537  ;
+  0.169933732064   ;  0.130190978230   ;  0.099652303426   ;  0.075428246546  ;
+  0.0555635614051  ;  0.0386791283055  ;  0.0237550435652  ;  0.010006278387  ;
+  |]
+
+let expo_sg_inv =
+  Array.map (fun x -> 1. /. (x *. expo_s *. expo_s))
+   [| 84135.654509    ;  2971.58727634    ;  474.716025959    ;  130.676724560    ;
+      47.3938388887   ;  20.2078651631    ;  9.5411021938     ;  4.8109546955     ;
+      2.52795733067   ;  1.35894103210    ;  0.73586710268    ;  0.39557629706    ;
+      0.20785895177   ;  0.104809693858   ;  0.049485682527   ;  0.021099788990   ;
+      0.007652472186  ;  0.0021065225215  ;  0.0003365204879  ;  0.0000118855674  |]
+*)
+
+
+let class_of_contracted_shell_pair_couple shell_pair_couple = 
+    F12RR.contracted_class_shell_pair_couple 
+      expo_sg_inv coef_g shell_pair_couple
+
+
+
+
+module Zero_m = struct
+  let name = "F12"
+end 
+
+
+
+let filter_contracted_shell_pairs ?(cutoff=integrals_cutoff) shell_pairs =
+  List.map (fun pair ->
+      match Cspc.make ~cutoff pair pair with
+      | Some cspc ->
+        let cls = class_of_contracted_shell_pair_couple cspc in                               
+        (pair, Zmap.fold (fun key value accu -> max (abs_float value) accu) cls 0. )
+      (* TODO \sum_k |coef_k * integral_k| *)
+      | None -> (pair, -1.)
+    ) shell_pairs
+  |> List.filter (fun (_, schwartz_p_max) -> schwartz_p_max >= cutoff)
+  |> List.map fst
+
+
+(*  TODO                                                                                      
+  let filter_contracted_shell_pair_couples
+  ?(cutoff=integrals_cutoff) shell_pair_couples = 
+  List.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
+
+*)
+
+
+let store_class basis ?(cutoff=integrals_cutoff) data contracted_shell_pair_couple cls = 
+  let to_powers x = 
+    let open Zkey in
+    match to_powers x with
+    | Three x -> x
+    | _ -> assert false
+  in
+
+  let shell_p = Cspc.shell_pair_p contracted_shell_pair_couple
+  and shell_q = Cspc.shell_pair_q contracted_shell_pair_couple
+  in
+
+  let s =
+    Overlap.of_basis basis
+    |> Overlap.matrix
+  in
+  let lambda_inv = 1. /. expo_s in
+  Array.iteri (fun i_c powers_i ->
+      let i_c = Cs.index (Csp.shell_a shell_p) + i_c + 1 in
+      let xi = to_powers powers_i in
+      Array.iteri (fun j_c powers_j ->
+          let j_c = Cs.index (Csp.shell_b shell_p) + j_c + 1 in
+          let xj = to_powers powers_j in
+          Array.iteri (fun k_c powers_k ->
+              let k_c = Cs.index (Csp.shell_a shell_q) + k_c + 1 in 
+              let xk = to_powers powers_k in
+              Array.iteri (fun l_c powers_l ->
+                  let l_c = Cs.index (Csp.shell_b shell_q) + l_c + 1 in
+                  let xl = to_powers powers_l in
+                  let key = Zkey.of_powers_twelve  xi xj xk xl in
+                  let value = Zmap.find cls key in
+                  lambda_inv *. s.{i_c,j_c} *. s.{k_c,l_c} -. value
+                  |> set_chem data i_c j_c k_c l_c
+                ) (Cs.zkey_array (Csp.shell_b shell_q))
+            ) (Cs.zkey_array (Csp.shell_a shell_q))
+        ) (Cs.zkey_array (Csp.shell_b shell_p))
+    ) (Cs.zkey_array (Csp.shell_a shell_p))
+
+
+
+
+
+
+let of_basis_serial basis = 
+
+  let n = Bs.size basis
+  and shell = Bs.contracted_shells basis
+  in
+
+  let eri_array = 
+    Fis.create ~size:n `Dense
+(*
+  Fis.create ~size:n `Sparse
+*)
+  in
+
+  let t0 = Unix.gettimeofday () in
+
+  let shell_pairs =
+    Csp.of_contracted_shell_array shell
+    |> filter_contracted_shell_pairs ~cutoff
+  in
+
+  Printf.printf "%d significant shell pairs computed in %f seconds\n"
+    (List.length shell_pairs) (Unix.gettimeofday () -. t0);
+
+
+  let t0 = Unix.gettimeofday () in
+  let ishell = ref 0 in
+
+  List.iter (fun shell_p ->
+      let () = 
+        if (Cs.index (Csp.shell_a shell_p) > !ishell) then
+          (ishell := Cs.index (Csp.shell_a shell_p) ; print_int !ishell ; print_newline ())
+      in
+
+      let sp =
+        Csp.shell_pairs shell_p
+      in
+
+      try
+        List.iter (fun shell_q ->
+            let () = 
+              if  Cs.index (Csp.shell_a shell_q) >
+                  Cs.index (Csp.shell_a shell_p) then
+                raise Exit
+            in
+            let sq = Csp.shell_pairs shell_q in
+            let cspc = 
+              if Array.length sp < Array.length sq then
+                Cspc.make ~cutoff shell_p shell_q
+              else
+                Cspc.make ~cutoff shell_q shell_p
+            in
+
+            match cspc with
+            | Some cspc ->
+              let cls =
+                class_of_contracted_shell_pair_couple cspc
+              in
+              store_class basis ~cutoff eri_array cspc cls
+            | None -> ()
+          ) shell_pairs
+      with Exit -> ()
+    ) shell_pairs ;
+  Printf.printf "Computed ERIs in %f seconds\n%!" (Unix.gettimeofday () -. t0);
+  eri_array
+
+
+
+
+
+
+
+(* Parallel functions *)
+
+
+
+let of_basis_parallel basis = 
+
+  let n = Bs.size basis
+  and shell = Bs.contracted_shells basis
+  in
+
+  let store_class_parallel
+      ?(cutoff=integrals_cutoff) contracted_shell_pair_couple cls = 
+    let to_powers x = 
+      let open Zkey in
+      match to_powers x with
+      | Three x -> x
+      | _ -> assert false
+    in
+
+    let shell_p = Cspc.shell_pair_p contracted_shell_pair_couple
+    and shell_q = Cspc.shell_pair_q contracted_shell_pair_couple
+    in
+
+    let result = ref [] in
+    Array.iteri (fun i_c powers_i ->
+        let i_c = Cs.index (Csp.shell_a shell_p) + i_c + 1 in
+        let xi = to_powers powers_i in
+        Array.iteri (fun j_c powers_j ->
+            let j_c = Cs.index (Csp.shell_b shell_p) + j_c + 1 in
+            let xj = to_powers powers_j in
+            Array.iteri (fun k_c powers_k ->
+                let k_c = Cs.index (Csp.shell_a shell_q) + k_c + 1 in 
+                let xk = to_powers powers_k in
+                Array.iteri (fun l_c powers_l ->
+                    let l_c = Cs.index (Csp.shell_b shell_q) + l_c + 1 in
+                    let xl = to_powers powers_l in
+                    let key = Zkey.of_powers_twelve  xi xj xk xl in
+                    let value = Zmap.find cls key in
+                    result := (i_c, j_c, k_c, l_c, value) :: !result
+                  ) (Cs.zkey_array (Csp.shell_b shell_q))
+              ) (Cs.zkey_array (Csp.shell_a shell_q))
+          ) (Cs.zkey_array (Csp.shell_b shell_p))
+      ) (Cs.zkey_array (Csp.shell_a shell_p));
+    !result
+  in
+
+
+
+  let t0 = Unix.gettimeofday () in
+
+  let shell_pairs =
+    Csp.of_contracted_shell_array shell
+    |> filter_contracted_shell_pairs ~cutoff
+  in
+
+  if Parallel.master then
+    Printf.printf "%d significant shell pairs computed in %f seconds\n"
+      (List.length shell_pairs) (Unix.gettimeofday () -. t0);
+
+
+  let t0 = Unix.gettimeofday () in
+  let ishell = ref max_int in
+
+  let input_stream = Stream.of_list (List.rev shell_pairs) in
+
+  let f shell_p =
+    let () = 
+      if Parallel.rank < 2 && Cs.index (Csp.shell_a shell_p) < !ishell then
+        (ishell := Cs.index (Csp.shell_a shell_p) ; print_int !ishell ; print_newline ())
+    in
+
+    let sp =
+      Csp.shell_pairs shell_p
+    in
+
+    let result = ref [] in
+    try
+      List.iter (fun shell_q ->
+          let () = 
+            if  Cs.index (Csp.shell_a shell_q) >
+                Cs.index (Csp.shell_a shell_p) then
+              raise Exit
+          in
+          let sq = Csp.shell_pairs shell_q in
+          let cspc = 
+            if Array.length sp < Array.length sq then
+              Cspc.make ~cutoff shell_p shell_q
+            else
+              Cspc.make ~cutoff shell_q shell_p
+          in
+
+          match cspc with
+          | Some cspc ->
+            let cls =
+              class_of_contracted_shell_pair_couple cspc
+            in
+            result := (store_class_parallel ~cutoff cspc cls) :: !result;
+          | None -> ()
+        ) shell_pairs;
+      raise Exit
+    with Exit -> List.concat !result |> Array.of_list
+  in
+
+  let eri_array = 
+    if Parallel.master then
+      Fis.create ~size:n `Dense
+    else
+      Fis.create ~size:0 `Dense
+  in
+  let s =
+    Overlap.of_basis basis
+    |> Overlap.matrix
+  in
+  let lambda_inv = 1. /. expo_s in
+  Farm.run ~ordered:false ~f input_stream
+  |> Stream.iter (fun l ->
+      Array.iter (fun (i_c,j_c,k_c,l_c,value) ->
+        lambda_inv *. s.{i_c,j_c} *. s.{k_c,l_c} -. value
+        |> set_chem eri_array i_c j_c k_c l_c) l);
+
+  if Parallel.master then
+    Printf.printf
+      "Computed %s Integrals in parallel in %f seconds\n%!" Zero_m.name (Unix.gettimeofday () -. t0);
+  Parallel.broadcast (lazy eri_array)
+
+
+
+let of_basis = 
+  match Parallel.size with
+  | 1 -> of_basis_serial
+  | _ -> of_basis_parallel
+
+
+
+
+

CI/CI.ml

@@ -410,14 +410,11 @@ let make ?(n_states=1) det_space =
 
 
 let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
-    list_holes1  list_particles1 ?(unique=true)
+    list_holes  list_particles ?(unique=true) is_internal
-    list_holes2  list_particles2 is_internal
     i_o1_alfa  alfa_o2_i  w_alfa  psi0 =
 
-  let list_holes1     = Array.of_list list_holes1 in
+  let list_holes     = Array.of_list list_holes in
-  let list_particles1 = Array.of_list list_particles1 in
+  let list_particles = Array.of_list list_particles in
-  let list_holes2     = Array.of_list list_holes2 in
-  let list_particles2 = Array.of_list list_particles2 in
 
   let psi0 =
     let stream = 
@@ -473,11 +470,7 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
           accu +. coef.{1} *. (i_o1_alfa det alfa)) 0. psi_filtered
     in
 
-    let alfa_h_psi =
+    let alfa_h_psi alfa =
-      if symmetric then 
-        psi_h_alfa 
-      else
-        fun alfa ->
       List.fold_left (fun accu (det, coef) ->
 (* Single state here *)
               accu +. coef.{1} *. (alfa_o2_i alfa det)) 0. psi_filtered
@@ -529,12 +522,12 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
                                         accu
                                       else
                                         accu +. w_alfa alfa *. psi_h_alfa_alfa_h_psi alfa 
-                                ) 0. list_holes1
+                                ) 0. list_holes
-                          ) 0. list_particles1
+                          ) 0. list_particles
                       in
                       accu +. single +. double
-                ) 0. list_holes2
+                ) 0. list_holes
-            ) 0. list_particles2
+            ) 0. list_particles
         ) 0. [ Spin.Alfa ; Spin.Beta ]
     in
 
@@ -563,13 +556,13 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
                                   accu
                                 else
                                   accu +. w_alfa alfa *. psi_h_alfa_alfa_h_psi alfa 
-                          ) 0. list_holes1
+                          ) 0. list_holes
-                      ) 0. list_particles1
+                      ) 0. list_particles
                   in
 
                   accu +. double_ab 
-            ) 0. list_holes2
+            ) 0. list_holes
-        ) 0. list_particles2
+        ) 0. list_particles
     in
     same_spin +. opposite_spin
   in
@@ -580,6 +573,68 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
 
 
 
+let second_order_sum2 { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
+    list_holes  list_particles i_o1_alfa e0 psi0 =
+
+  let psi0 =
+    let stream = 
+      Ds.determinant_stream det_space
+    in
+    Array.init (Ds.size det_space) (fun i ->
+        (Stream.next stream), (Mat.copy_row psi0 (i+1)) )
+  in
+
+  let determinants =
+
+    Ds.determinants_array det_space
+    |> Array.to_list
+    |> List.map (fun det_i ->
+        [ Spin.Alfa ; Spin.Beta ] 
+        |>  List.map (fun spin ->
+            List.map (fun particle ->
+              List.map (fun hole ->
+                  [ [ Determinant.single_excitation spin hole particle det_i ] ;
+                    List.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.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
+
+  List.fold_left (fun accu alfa ->
+    let alfa_o2 = i_o1_alfa alfa in
+    let a_h_psi = 
+      Array.fold_left (fun accu (det,ci) -> ci.{1} *. (alfa_o2 det)) 0. psi0
+    in
+    accu +. (a_h_psi *. a_h_psi) /. (e0 -. (alfa_o2 alfa))
+    ) 0. determinants
+
+
+
 
 let is_internal det_space =
     let m l =
@@ -608,7 +663,7 @@ let is_internal det_space =
         Z.logand neg_active_mask beta = occ_mask
 
 
-let pt2_en ci =
+let _pt2_en ci =
 
   let mo_basis  = Ds.mo_basis ci.det_space in
   let psi0, e0 = Parallel.broadcast ci.eigensystem in
@@ -668,11 +723,25 @@ let pt2_en ci =
     [ MOClass.active_mos mo_class ; MOClass.virtual_mos mo_class ]
   in
 
-  second_order_sum  ci  list_holes  list_particles  list_holes  list_particles
+  second_order_sum  ci  list_holes  list_particles  
     (is_internal ci.det_space) i_o1_alfa  i_o1_alfa  w_alfa  psi0
   |> List.fold_left (+.) 0.
 
 
+let pt2_en ci =
+
+  let mo_basis  = Ds.mo_basis ci.det_space in
+  let psi0, e0 = Parallel.broadcast ci.eigensystem in
+
+  let i_o1_alfa = h_ij mo_basis in
+
+  let mo_class = mo_class ci in
+  let list_holes = List.concat
+    [ MOClass.inactive_mos mo_class ; MOClass.active_mos mo_class ]
+  and list_particles = List.concat 
+    [ MOClass.active_mos mo_class ; MOClass.virtual_mos mo_class ]
+  in
+  second_order_sum2 ci list_holes  list_particles i_o1_alfa e0.{1} psi0 
 
 
 
@@ -701,7 +770,7 @@ let pt2_mp ci =
   in
 
   let psi0, _ = Parallel.broadcast ci.eigensystem in
-  second_order_sum  ci  list_holes  list_particles  list_holes  list_particles
+  second_order_sum  ci  list_holes  list_particles  
     (is_internal ci.det_space) i_o1_alfa  i_o1_alfa  w_alfa  psi0
   |> List.fold_left (+.) 0.
 
@@ -722,7 +791,7 @@ let variance ci =
     [ MOClass.active_mos mo_class ; MOClass.virtual_mos mo_class ]
   in
 
-  second_order_sum  ci  list_holes  list_particles  list_holes  list_particles
+  second_order_sum  ci  list_holes  list_particles  
     (is_internal ci.det_space) i_o1_alfa  i_o1_alfa  w_alfa  psi0
   |> List.fold_left (+.) 0.
 

CI/F12CI.ml

@@ -1,6 +1,3 @@
-let debug s =
-  Printf.printf "%s\n%!" s;
-
 open Lacaml.D
 
 type t = 
@@ -43,13 +40,22 @@ let h_ij mo_basis ki kj =
   |> List.hd
 
 
-let f_ij mo_basis ki kj =
+let f_ij gamma mo_basis ki kj =
   let integrals =
     List.map (fun f -> f mo_basis)
       [ f12_integrals ]
   in
+  let integral = 
     CIMatrixElement.make integrals ki kj 
     |> List.hd
+  in
+  gamma *. integral 
+  (*
+  match Determinant.degrees ki kj with
+  | (2,0) 
+  | (0,2) -> 0.5 *. gamma *. integral
+  | _     -> gamma *. integral 
+  *)
 
 
 let is_internal det_space =
@@ -75,9 +81,7 @@ let is_internal det_space =
         Z.logand aux_mask beta = Z.zero 
 
 
-let dressing_vector aux_basis f12_amplitudes ci =
+let dressing_vector gamma aux_basis f12_amplitudes ci =
-
-debug "Computing dressing vector";
 
 (*
   let i_o1_alfa = h_ij aux_basis in
@@ -93,9 +97,6 @@ debug "Computing dressing vector";
   and list_particles1 = List.concat
     [ MOClass.active_mos mo_class ; MOClass.virtual_mos mo_class ; MOClass.auxiliary_mos mo_class ]
   in
-(*
-Util.debug_matrix "f12 amplitudes" f12_amplitudes;
-*)
   (* Single state here *)
   let result = 
     CI.second_order_sum  ci  list_holes  list_particles1  list_holes  list_particles2
@@ -130,7 +131,7 @@ Util.debug_matrix "f12 amplitudes" f12_amplitudes;
   let m_F_aux = 
     Array.map (fun ki ->
       Array.map (fun kj ->
-        f_ij aux_basis ki kj
+        f_ij gamma aux_basis ki kj
       ) out_dets
     ) in_dets
     |> Mat.of_array
@@ -148,7 +149,7 @@ Util.debug_matrix "f12 amplitudes" f12_amplitudes;
 
 let make ~simulation ?(threshold=1.e-12) ?(frozen_core=true) ~mo_basis ~aux_basis_filename () =
 
-  let gamma = 1.0 in
+  let gamma = 0.5 in
 
   let mo_num = MOBasis.size mo_basis in
 
@@ -188,24 +189,19 @@ let make ~simulation ?(threshold=1.e-12) ?(frozen_core=true) ~mo_basis ~aux_basi
     (* While in a sequential region, initiate the parallel
        4-idx transformation to avoid nested parallel jobs
     *)
-debug "Four-idx transform of f12 intergals";
     ignore @@ MOBasis.f12_ints aux_basis;
 
     let f = fun ki kj ->
       if ki <> kj then
-        gamma *. (f_ij aux_basis ki kj)
+        (f_ij gamma aux_basis ki kj)
       else
-        1. +. gamma *. (f_ij aux_basis ki kj)
+        1. +. (f_ij gamma aux_basis ki kj)
     in
     let m_F = 
       CI.create_matrix_spin f det_space
       |> Lazy.force
     in
     fun ci_coef -> 
-(*
-Util.debug_matrix "F" (Matrix.to_mat m_F);
-debug "Solving linear system";
-*)
     Matrix.ax_eq_b m_F (Matrix.dense_of_mat ci_coef)
     |> Matrix.to_mat
   in
@@ -222,14 +218,10 @@ debug "Solving linear system";
     let m_H =
       Lazy.force ci.CI.m_H
     in
-(*
-Util.debug_matrix "H" (Matrix.to_mat m_H);
-*)
 
     let rec iteration ?(state=1) psi = 
-debug "Iteration";
       let delta = 
-         dressing_vector aux_basis (f12_amplitudes psi) ci
+         dressing_vector gamma aux_basis (f12_amplitudes psi) ci
       in
 
       let f = 1.0 /. psi.{1,1} in
@@ -273,7 +265,9 @@ TODO SINGLE STATE HERE
       in
 
       let eigenvectors, eigenvalues = 
+        Parallel.broadcast (lazy (
           Davidson.make ~threshold:1.e-6 ~guess:psi ~n_states diagonal matrix_prod
+          ))
       in
 
       let conv =
@@ -281,7 +275,8 @@ TODO SINGLE STATE HERE
           (Mat.to_col_vecs psi).(0)
           (Mat.to_col_vecs eigenvectors).(0) )
       in
-      Printf.printf "Convergence : %e  %f\n" conv (eigenvalues.{1} +. e_shift);
+      if Parallel.master then
+        Printf.printf "F12 Convergence : %e  %f\n" conv (eigenvalues.{1} +. e_shift);
 
       if conv > threshold then
         iteration eigenvectors