From 2735688efee5aef632edbdc9ba68a61dc14946af Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 22 Mar 2019 23:42:35 +0100 Subject: [PATCH 1/4] Davidson F12CI OK --- CI/F12CI.ml | 76 +++++++++++++++++++---------------------------------- 1 file changed, 27 insertions(+), 49 deletions(-) diff --git a/CI/F12CI.ml b/CI/F12CI.ml index b338be3..051ccf2 100644 --- a/CI/F12CI.ml +++ b/CI/F12CI.ml @@ -241,47 +241,7 @@ Util.debug_matrix "psi" psi; Format.printf "Amplitude (1,1) : %f@." (f12_amplitudes psi).{1,1}; Format.printf "Dressing vector(1,1) : %f@." (Matrix.get delta 1 1); -(*------ -TODO SINGLE STATE HERE -*) - let m_H_dressed = Matrix.to_mat m_H in let f = 1.0 /. psi.{1,1} in - Util.list_range 1 (Mat.dim1 psi) - |> List.iter (fun i -> - let delta = (Matrix.get delta i 1) *. f in - m_H_dressed.{i,1} <- m_H_dressed.{i,1} +. delta; - if i <> 1 then - begin - m_H_dressed.{1,i} <- m_H_dressed.{1,i} +. delta; - m_H_dressed.{1,1} <- m_H_dressed.{1,1} -. delta *. psi.{i,1} *. f; - end - (* DIAGONAL DRESSING - m_H_dressed.{i,i} <- m_H_dressed.{i,i} +. (Matrix.get delta i 1) /. (psi.{i,1} +. 1.e-5); - *) - ); - let eigenvectors, eigenvalues = - Util.diagonalize_symm m_H_dressed - in - let conv = - 1.0 -. abs_float ( dot - (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 conv > threshold then - iteration eigenvectors - else - let eigenvalues = - Vec.map (fun x -> x +. e_shift) eigenvalues - in - eigenvectors, eigenvalues - in - iteration ci_coef -(* -------- *) -(* - let n_states = ci.CI.n_states in - let f = 1.0 /. (psi.{1,1} ) in let delta_00 = Util.list_range 2 (Mat.dim1 psi) |> List.fold_left (fun accu i -> accu +. @@ -289,19 +249,38 @@ TODO SINGLE STATE HERE in let delta = Matrix.to_mat delta in delta.{1,1} <- delta.{1,1} -. delta_00; + +(*------ +TODO SINGLE STATE HERE +*) + let n_states = ci.CI.n_states in let diagonal = - Vec.init (Matrix.dim1 m_H) (fun i -> Matrix.get m_H i i +. (if i=1 then delta.{1,1} *. psi.{1,1} else 0.) ) + Vec.init (Matrix.dim1 m_H) (fun i -> + if i = 1 then + Matrix.get m_H i i +. delta.{1,1} *. f + else + Matrix.get m_H i i + ) in - let delta = Matrix.dense_of_mat delta in let matrix_prod c = - Matrix.add - (Matrix.mm ~transa:`T m_H c) - delta + let w = + Matrix.mm ~transa:`T m_H c + |> Matrix.to_mat + in + let c11 = Matrix.get c 1 1 in + Util.list_range 1 (Mat.dim1 w) + |> List.iter (fun i -> + let dci = + delta.{i,1} *. f ; + in + w.{i,1} <- w.{i,1} +. dci *. c11; + if (i <> 1) then + w.{1,1} <- w.{1,1} +. dci *. (Matrix.get c i 1); + ); + Matrix.dense_of_mat w in let eigenvectors, eigenvalues = - Parallel.broadcast (lazy ( - Davidson.make ~threshold:1.e-6 ~guess:psi ~n_states diagonal matrix_prod - )) + Davidson.make ~threshold:1.e-6 ~guess:psi ~n_states diagonal matrix_prod in let conv = 1.0 -. abs_float ( dot @@ -318,7 +297,6 @@ TODO SINGLE STATE HERE eigenvectors, eigenvalues in iteration ci_coef -*) ) in From 494b0a2e3cad3e704e691c1ad3f077d283a9a328 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sat, 23 Mar 2019 01:06:38 +0100 Subject: [PATCH 2/4] Searching bug in CAS+EN2 --- CI/CI.ml | 55 ++++++++++++++++++++++++++--------------------------- CI/F12CI.ml | 13 ++++--------- 2 files changed, 31 insertions(+), 37 deletions(-) diff --git a/CI/CI.ml b/CI/CI.ml index ea31578..f8fa0a7 100644 --- a/CI/CI.ml +++ b/CI/CI.ml @@ -414,14 +414,12 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states } list_holes2 list_particles2 is_internal i_o1_alfa alfa_o2_i w_alfa psi0 = - let list_holes1 = Array.of_list list_holes1 - and list_holes2 = Array.of_list list_holes2 - and list_particles1 = Array.of_list list_particles1 - and list_particles2 = Array.of_list list_particles2 - in + 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 psi0 = - let stream = Ds.determinant_stream det_space in @@ -448,7 +446,8 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states } || aux (j-1) in aux (i-1) - ) else + ) + else is_internal in @@ -513,18 +512,17 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states } let double = Array.fold_left (fun accu particle' -> - if particle' > particle || particle' = hole then + if particle' >= particle || particle' = hole then accu else accu +. Array.fold_left (fun accu hole' -> - if hole' = particle' || hole' = particle || hole' < hole then + if hole' = particle' || hole' = particle || hole' <= hole then accu else let alfa = - Determinant.double_excitation - spin hole particle - spin hole' particle' det_i + Determinant.single_excitation + spin hole' particle' alfa in if Determinant.is_none alfa || already_generated alfa then @@ -550,25 +548,26 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states } in if Determinant.is_none alfa then accu else - let double = + let double_ab = Array.fold_left (fun accu particle' -> - accu +. - Array.fold_left (fun accu hole' -> - if hole' = particle' then accu else - let alfa = - Determinant.double_excitation - Spin.Alfa hole particle - Spin.Beta hole' particle' det_i - in - if Determinant.is_none alfa || - already_generated alfa then - accu - else - accu +. w_alfa alfa *. psi_h_alfa_alfa_h_psi alfa + accu +. + Array.fold_left (fun accu hole' -> + if hole' = particle' then accu else + let alfa = + Determinant.double_excitation + Spin.Alfa hole particle + Spin.Beta hole' particle' det_i + in + if Determinant.is_none alfa || + already_generated alfa then + accu + else + accu +. w_alfa alfa *. psi_h_alfa_alfa_h_psi alfa ) 0. list_holes1 ) 0. list_particles1 in - accu +. double + + accu +. double_ab ) 0. list_holes2 ) 0. list_particles2 in @@ -757,7 +756,7 @@ let pt2_en_reference ci = in Array.map (fun ki -> Array.mapi (fun j kj -> - (h_ij aux_basis ki kj) /. (e0.{1} -. h_aa.(j)) + (h_ij aux_basis ki kj) /. (e0.{1} -. h_aa.(j)) ) out_dets ) in_dets |> Mat.of_array diff --git a/CI/F12CI.ml b/CI/F12CI.ml index 051ccf2..f81f17e 100644 --- a/CI/F12CI.ml +++ b/CI/F12CI.ml @@ -197,7 +197,6 @@ debug "Four-idx transform of f12 intergals"; else 1. +. gamma *. (f_ij aux_basis ki kj) in -debug "Computing F matrix"; let m_F = CI.create_matrix_spin f det_space |> Lazy.force @@ -228,18 +227,10 @@ Util.debug_matrix "H" (Matrix.to_mat m_H); *) let rec iteration ?(state=1) psi = -(* debug "Iteration"; -Util.debug_matrix "T" (f12_amplitudes psi); -*) let delta = dressing_vector aux_basis (f12_amplitudes psi) ci in - (* -Util.debug_matrix "psi" psi; -*) -Format.printf "Amplitude (1,1) : %f@." (f12_amplitudes psi).{1,1}; -Format.printf "Dressing vector(1,1) : %f@." (Matrix.get delta 1 1); let f = 1.0 /. psi.{1,1} in let delta_00 = @@ -262,6 +253,7 @@ TODO SINGLE STATE HERE Matrix.get m_H i i ) in + let matrix_prod c = let w = Matrix.mm ~transa:`T m_H c @@ -279,15 +271,18 @@ TODO SINGLE STATE HERE ); Matrix.dense_of_mat w in + let eigenvectors, eigenvalues = Davidson.make ~threshold:1.e-6 ~guess:psi ~n_states diagonal matrix_prod in + let conv = 1.0 -. abs_float ( dot (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 conv > threshold then iteration eigenvectors else From 0a6b8e30a473555eddddee5f6947c31023a75384 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sat, 23 Mar 2019 14:54:59 +0100 Subject: [PATCH 3/4] 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 From 48939ed7cf45d0c55f22cc08f6fc90e15dc384a4 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sat, 23 Mar 2019 15:54:46 +0100 Subject: [PATCH 4/4] Fixed open shells CAS --- CI/CI.ml | 2 +- CI/DeterminantSpace.ml | 16 ++++++++-------- CI/DeterminantSpace.mli | 8 ++++---- CI/F12CI.ml | 2 +- CI/SpindeterminantSpace.ml | 6 +++--- CI/SpindeterminantSpace.mli | 4 ++-- MOBasis/MOClass.ml | 10 ++++++---- MOBasis/MOClass.mli | 2 +- Perturbation/MP2.ml | 4 ++-- run_cas.ml | 2 +- run_mp2.ml | 2 +- 11 files changed, 30 insertions(+), 28 deletions(-) diff --git a/CI/CI.ml b/CI/CI.ml index bfe8742..4a5cb1d 100644 --- a/CI/CI.ml +++ b/CI/CI.ml @@ -805,7 +805,7 @@ let pt2_en_reference ci = let aux_basis = mo_basis in let ds = - DeterminantSpace.fci_of_mo_basis ~frozen_core:true aux_basis + DeterminantSpace.fci_of_mo_basis ~frozen_core:false aux_basis in let out_dets = ds diff --git a/CI/DeterminantSpace.ml b/CI/DeterminantSpace.ml index d9cb913..48b1d02 100644 --- a/CI/DeterminantSpace.ml +++ b/CI/DeterminantSpace.ml @@ -292,21 +292,21 @@ let arbitrary_of_mo_basis mo_basis f = -let cas_of_mo_basis mo_basis n m = +let cas_of_mo_basis mo_basis ~frozen_core n m = let f n_alfa = - Ss.cas_of_mo_basis mo_basis n_alfa n m + Ss.cas_of_mo_basis mo_basis ~frozen_core n_alfa n m in spin_of_mo_basis mo_basis f -let fci_of_mo_basis ?(frozen_core=true) mo_basis = +let fci_of_mo_basis mo_basis ~frozen_core = let f n_alfa = - Ss.fci_of_mo_basis ~frozen_core mo_basis n_alfa + Ss.fci_of_mo_basis mo_basis ~frozen_core n_alfa in spin_of_mo_basis mo_basis f -let fci_f12_of_mo_basis ?(frozen_core=true) mo_basis mo_num = +let fci_f12_of_mo_basis mo_basis ~frozen_core mo_num = let s = MOBasis.simulation mo_basis in let e = Simulation.electrons s in let n_alfa = Electrons.n_alfa e @@ -321,7 +321,7 @@ let fci_f12_of_mo_basis ?(frozen_core=true) mo_basis mo_num = (mo_num - n_core) in let f n_alfa = - Ss.cas_of_mo_basis mo_basis n_alfa n m + Ss.cas_of_mo_basis mo_basis ~frozen_core n_alfa n m in let r = spin_of_mo_basis mo_basis f @@ -335,9 +335,9 @@ let fci_f12_of_mo_basis ?(frozen_core=true) mo_basis mo_num = |> MOClass.of_list } -let cas_f12_of_mo_basis mo_basis n m mo_num = +let cas_f12_of_mo_basis mo_basis ~frozen_core n m mo_num = let f n_alfa = - Ss.cas_of_mo_basis mo_basis n_alfa n m + Ss.cas_of_mo_basis mo_basis ~frozen_core n_alfa n m in let r = spin_of_mo_basis mo_basis f diff --git a/CI/DeterminantSpace.mli b/CI/DeterminantSpace.mli index e8aa4ee..1c189a3 100644 --- a/CI/DeterminantSpace.mli +++ b/CI/DeterminantSpace.mli @@ -50,20 +50,20 @@ val fock_diag : t -> Determinant.t -> float array * float array *) -val fci_of_mo_basis : ?frozen_core:bool -> MOBasis.t -> t +val fci_of_mo_basis : MOBasis.t -> frozen_core:bool -> t (** Creates a space of all possible ways to put [n_alfa] electrons in the {% $\alpha$ %} [Active] MOs and [n_beta] electrons in the {% $\beta$ %} [Active] MOs. All other MOs are untouched. *) -val cas_of_mo_basis : MOBasis.t -> int -> int -> t +val cas_of_mo_basis : MOBasis.t -> frozen_core:bool -> int -> int -> t (** Creates a CAS(n,m) space of determinants. *) -val fci_f12_of_mo_basis : ?frozen_core:bool -> MOBasis.t -> int -> t +val fci_f12_of_mo_basis : MOBasis.t -> frozen_core:bool -> int -> t (** Creates the active space to perform a FCI-F12 with an auxiliary basis set. *) -val cas_f12_of_mo_basis : MOBasis.t -> int -> int -> int -> t +val cas_f12_of_mo_basis : MOBasis.t -> frozen_core:bool -> int -> int -> int -> t (** [cas_of_mo_basis mo_basis m n mo_num] Creates a CAS(n,m) space of determinants with an auxiliary basis set defined as the MOs from [mo_num+1] to [MOBasis.size mo_basis]. diff --git a/CI/F12CI.ml b/CI/F12CI.ml index 3d69de1..e9e493f 100644 --- a/CI/F12CI.ml +++ b/CI/F12CI.ml @@ -147,7 +147,7 @@ let dressing_vector gamma aux_basis f12_amplitudes ci = -let make ~simulation ?(threshold=1.e-12) ?(frozen_core=true) ~mo_basis ~aux_basis_filename () = +let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filename () = let gamma = 0.5 in diff --git a/CI/SpindeterminantSpace.ml b/CI/SpindeterminantSpace.ml index a88fc65..913c94a 100644 --- a/CI/SpindeterminantSpace.ml +++ b/CI/SpindeterminantSpace.ml @@ -13,7 +13,7 @@ let mo_basis t = t.mo_basis let mo_class t = t.mo_class let size t = Array.length t.spin_determinants -let fci_of_mo_basis ?(frozen_core=true) mo_basis elec_num = +let fci_of_mo_basis ~frozen_core mo_basis elec_num = let mo_num = MOBasis.size mo_basis in let mo_class = MOClass.fci ~frozen_core mo_basis in let m l = @@ -35,9 +35,9 @@ let fci_of_mo_basis ?(frozen_core=true) mo_basis elec_num = { elec_num ; mo_basis ; mo_class ; spin_determinants } -let cas_of_mo_basis mo_basis elec_num n m = +let cas_of_mo_basis mo_basis ~frozen_core elec_num n m = let mo_num = MOBasis.size mo_basis in - let mo_class = MOClass.cas_sd mo_basis n m in + let mo_class = MOClass.cas_sd ~frozen_core mo_basis n m in let m l = List.fold_left (fun accu i -> let j = i-1 in Z.(logor accu (shift_left one j)) ) Z.zero l diff --git a/CI/SpindeterminantSpace.mli b/CI/SpindeterminantSpace.mli index e1c2653..555abcd 100644 --- a/CI/SpindeterminantSpace.mli +++ b/CI/SpindeterminantSpace.mli @@ -24,12 +24,12 @@ val mo_basis : t -> MOBasis.t (** {1 Creation} *) -val fci_of_mo_basis : ?frozen_core:bool -> MOBasis.t -> int -> t +val fci_of_mo_basis : frozen_core:bool -> MOBasis.t -> int -> t (** Create a space of all possible ways to put [n_elec-ncore] electrons in the [Active] MOs. All other MOs are untouched. *) -val cas_of_mo_basis : MOBasis.t -> int -> int -> int -> t +val cas_of_mo_basis : MOBasis.t -> frozen_core:bool -> int -> int -> int -> t (** [cas_of_mo_basis mo_basis n_elec n m] creates a CAS(n,m) space of [Active] MOs. The unoccupied MOs are [Virtual], and the occupied MOs are [Core] and [Inactive]. diff --git a/MOBasis/MOClass.ml b/MOBasis/MOClass.ml index 6959aab..242295d 100644 --- a/MOBasis/MOClass.ml +++ b/MOBasis/MOClass.ml @@ -104,17 +104,19 @@ let fci ?(frozen_core=true) mo_basis = ) -let cas_sd mo_basis n m = +let cas_sd mo_basis ~frozen_core n m = let mo_num = MOBasis.size mo_basis in let n_alfa = MOBasis.simulation mo_basis |> Simulation.electrons |> Electrons.n_alfa in let n_beta = MOBasis.simulation mo_basis |> Simulation.electrons |> Electrons.n_beta in let n_unpaired = n_alfa - n_beta in - let n_alfa_in_cas = (n - n_unpaired)/2 in + let n_alfa_in_cas = (n - n_unpaired)/2 + n_unpaired in let last_inactive = n_alfa - n_alfa_in_cas in let last_active = last_inactive + m in let ncore = - (Nuclei.small_core @@ Simulation.nuclei @@ MOBasis.simulation mo_basis) / 2 - |> min last_inactive + if frozen_core then + (Nuclei.small_core @@ Simulation.nuclei @@ MOBasis.simulation mo_basis) / 2 + |> min last_inactive + else 0 in of_list ( List.concat [ diff --git a/MOBasis/MOClass.mli b/MOBasis/MOClass.mli index adf5392..4e63c64 100644 --- a/MOBasis/MOClass.mli +++ b/MOBasis/MOClass.mli @@ -20,7 +20,7 @@ val fci : ?frozen_core:bool -> MOBasis.t -> t [n] lowest MOs are [Core] if [frozen_core = true]. *) -val cas_sd: MOBasis.t -> int -> int -> t +val cas_sd: MOBasis.t -> frozen_core:bool -> int -> int -> t (** [cas_sd mo_basis n m ] creates the MO classes for CAS(n,m) + SD calculations. lowest MOs are [Core], then all the next MOs are [Inactive], then [Active], then [Virtual]. diff --git a/Perturbation/MP2.ml b/Perturbation/MP2.ml index 98f481d..037cc94 100644 --- a/Perturbation/MP2.ml +++ b/Perturbation/MP2.ml @@ -1,6 +1,6 @@ type t = float -let make ?(frozen_core=true) hf = +let make ~frozen_core hf = let mo_basis = MOBasis.of_hartree_fock hf in @@ -8,7 +8,7 @@ let make ?(frozen_core=true) hf = MOBasis.mo_energies mo_basis in let mo_class = - MOClass.cas_sd mo_basis 0 0 + MOClass.cas_sd mo_basis ~frozen_core 0 0 |> MOClass.to_list in let eri = diff --git a/run_cas.ml b/run_cas.ml index a58a330..89a91e2 100644 --- a/run_cas.ml +++ b/run_cas.ml @@ -67,7 +67,7 @@ let () = in let space = - DeterminantSpace.cas_of_mo_basis mos n m + DeterminantSpace.cas_of_mo_basis mos ~frozen_core:true n m in let ci = CI.make space in Format.fprintf ppf "CAS-CI energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s); diff --git a/run_mp2.ml b/run_mp2.ml index 6be4bf8..458caa2 100644 --- a/run_mp2.ml +++ b/run_mp2.ml @@ -73,7 +73,7 @@ let () = let e_hf = HartreeFock.energy hf in - let mp2 = MP2.make hf in + let mp2 = MP2.make ~frozen_core:true hf in Format.fprintf ppf "@[MP2 = %15.10f@]@." mp2; Format.fprintf ppf "@[E+MP2 = %15.10f@]@." (mp2 +. e_hf)