10
1
mirror of https://gitlab.com/scemama/QCaml.git synced 2024-11-07 22:53:41 +01:00

Fixed bug in F12 integrals

This commit is contained in:
Anthony Scemama 2019-03-25 10:54:23 +01:00
parent e21aea2b16
commit 075455fc9c
3 changed files with 26 additions and 12 deletions

View File

@ -14,8 +14,14 @@ module Fis = FourIdxStorage
include FourIdxStorage include FourIdxStorage
(*
type gaussian_geminal =
{
expo_s : float ;
coef_g :
*)
(** Exponent of the geminal *) (** Exponent of the geminal *)
let expo_s = 2.5 let expo_s = 1.
(** Coefficients and exponents of the Gaussian fit of the Slater Geminal*) (** Coefficients and exponents of the Gaussian fit of the Slater Geminal*)
let coef_g = let coef_g =
@ -121,7 +127,7 @@ let store_class basis ?(cutoff=integrals_cutoff) data contracted_shell_pair_coup
let xl = to_powers powers_l in let xl = to_powers powers_l in
let key = Zkey.of_powers_twelve xi xj xk xl in let key = Zkey.of_powers_twelve xi xj xk xl in
let value = Zmap.find cls key in let value = Zmap.find cls key in
lambda_inv *. s.{i_c,j_c} *. s.{k_c,l_c} -. value lambda_inv *. (s.{i_c,j_c} *. s.{k_c,l_c} -. value)
|> set_chem data i_c j_c k_c l_c |> set_chem data i_c j_c k_c l_c
) (Cs.zkey_array (Csp.shell_b shell_q)) ) (Cs.zkey_array (Csp.shell_b shell_q))
) (Cs.zkey_array (Csp.shell_a shell_q)) ) (Cs.zkey_array (Csp.shell_a shell_q))
@ -321,7 +327,7 @@ let of_basis_parallel basis =
Farm.run ~ordered:false ~f input_stream Farm.run ~ordered:false ~f input_stream
|> Stream.iter (fun l -> |> Stream.iter (fun l ->
Array.iter (fun (i_c,j_c,k_c,l_c,value) -> 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 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); |> set_chem eri_array i_c j_c k_c l_c) l);
if Parallel.master then if Parallel.master then

View File

@ -128,6 +128,7 @@ let dressing_vector gamma aux_basis f12_amplitudes ci =
in in
Printf.printf "Building matrix\n%!"; Printf.printf "Building matrix\n%!";
let m_H_aux, m_F_aux = let m_H_aux, m_F_aux =
let rec col_vecs_list accu_H accu_F = function let rec col_vecs_list accu_H accu_F = function
| [] -> | [] ->
@ -137,35 +138,39 @@ let dressing_vector gamma aux_basis f12_amplitudes ci =
let h, f = let h, f =
List.map (fun kj -> List.map (fun kj ->
match hf_ij aux_basis ki kj with match hf_ij aux_basis ki kj with
| [ a ; b ] -> a, b | [ a ; b ] -> a, gamma *. b
| _ -> assert false ) out_dets | _ -> assert false ) in_dets
|> List.split |> List.split
in in
col_vecs_list (h::accu_H) (f::accu_F) rest col_vecs_list (h::accu_H) (f::accu_F) rest
in in
let h, f = let h, f =
col_vecs_list [] [] in_dets col_vecs_list [] [] out_dets
in in
Mat.of_col_vecs_list h, Mat.of_col_vecs_list h,
Mat.of_col_vecs_list f Mat.of_col_vecs_list f
in in
(* (*
let m_H_aux = let m_H_aux =
Array.map (fun ki -> List.map (fun ki ->
Array.map (fun kj -> List.map (fun kj ->
h_ij aux_basis ki kj h_ij aux_basis ki kj
) out_dets ) out_dets
|> Array.of_list
) in_dets ) in_dets
|> Array.of_list
|> Mat.of_array |> Mat.of_array
in in
let m_F_aux = let m_F_aux =
Array.map (fun ki -> List.map (fun ki ->
Array.map (fun kj -> List.map (fun kj ->
f_ij gamma aux_basis ki kj f_ij gamma aux_basis ki kj
) out_dets ) out_dets
|> Array.of_list
) in_dets ) in_dets
|> Array.of_list
|> Mat.of_array |> Mat.of_array
in in
*) *)
@ -175,6 +180,9 @@ let dressing_vector gamma aux_basis f12_amplitudes ci =
gemm m_H_aux m_F_aux ~transb:`T gemm m_H_aux m_F_aux ~transb:`T
in in
Printf.printf "Done\n%!"; Printf.printf "Done\n%!";
Printf.printf "%d %d %d %d %d %d\n" (Mat.dim1 m_H_aux) (Mat.dim2 m_H_aux)
(Mat.dim1 m_F_aux) (Mat.dim2 m_F_aux)
(Mat.dim1 m_HF) (Mat.dim2 m_HF);
gemm m_HF f12_amplitudes gemm m_HF f12_amplitudes
|> Matrix.sparse_of_mat |> Matrix.sparse_of_mat

View File

@ -69,7 +69,7 @@ let () =
let space = let space =
DeterminantSpace.cas_of_mo_basis mos ~frozen_core:true n m DeterminantSpace.cas_of_mo_basis mos ~frozen_core:true n m
in in
let ci = CI.make ~frozen_core:true space in let ci = CI.make space in
Format.fprintf ppf "CAS-CI energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s); Format.fprintf ppf "CAS-CI energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s);
(* (*