From 0a6b8e30a473555eddddee5f6947c31023a75384 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sat, 23 Mar 2019 14:54:59 +0100 Subject: [PATCH] Fixed f12 integrals --- Basis/F12.ml | 342 ++++++++++++++++++++++++++++++++++++++++++++++++++- CI/CI.ml | 117 ++++++++++++++---- CI/F12CI.ml | 49 ++++---- 3 files changed, 455 insertions(+), 53 deletions(-) diff --git a/Basis/F12.ml b/Basis/F12.ml index 29d832f..2c3caf4 100644 --- a/Basis/F12.ml +++ b/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 + + + + + diff --git a/CI/CI.ml b/CI/CI.ml index f8fa0a7..bfe8742 100644 --- a/CI/CI.ml +++ b/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_holes2 list_particles2 is_internal + list_holes list_particles ?(unique=true) is_internal i_o1_alfa alfa_o2_i w_alfa psi0 = - let list_holes1 = Array.of_list list_holes1 in - let list_particles1 = Array.of_list list_particles1 in - let list_holes2 = Array.of_list list_holes2 in - let list_particles2 = Array.of_list list_particles2 in + let list_holes = Array.of_list list_holes in + let list_particles = Array.of_list list_particles in let psi0 = let stream = @@ -473,12 +470,8 @@ 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 = - if symmetric then - psi_h_alfa - else - fun alfa -> - List.fold_left (fun accu (det, coef) -> + let alfa_h_psi alfa = + List.fold_left (fun accu (det, coef) -> (* Single state here *) accu +. coef.{1} *. (alfa_o2_i alfa det)) 0. psi_filtered in @@ -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_particles1 + ) 0. list_holes + ) 0. list_particles in accu +. single +. double - ) 0. list_holes2 - ) 0. list_particles2 + ) 0. list_holes + ) 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_particles1 + ) 0. list_holes + ) 0. list_particles in accu +. double_ab - ) 0. list_holes2 - ) 0. list_particles2 + ) 0. list_holes + ) 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. diff --git a/CI/F12CI.ml b/CI/F12CI.ml index f81f17e..3d69de1 100644 --- a/CI/F12CI.ml +++ b/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 - CIMatrixElement.make integrals ki kj - |> List.hd + 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 = - -debug "Computing dressing vector"; +let dressing_vector gamma aux_basis f12_amplitudes ci = (* 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