From 938935425f809fa641f0f18353f5ebc038f40291 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 26 Mar 2020 19:49:29 +0100 Subject: [PATCH] Erf integrals correct, but mu needs to be checked --- Basis/AOBasis.ml | 5 ++--- Basis/TwoElectronRR.ml | 14 ++++++-------- Basis/TwoElectronRRVectorized.ml | 17 +++++------------ Basis/Zero_m_parameters.ml | 4 ++-- Simulation.ml | 1 + run_tests.ml | 4 ++-- 6 files changed, 18 insertions(+), 27 deletions(-) diff --git a/Basis/AOBasis.ml b/Basis/AOBasis.ml index 7f128a4..aec1689 100644 --- a/Basis/AOBasis.ml +++ b/Basis/AOBasis.ml @@ -87,8 +87,7 @@ let test_case name t = in let a = ERI.to_list a |> List.rev_map f |> List.rev in let r = ERI.to_list r |> List.rev_map f |> List.rev in - Printf.eprintf "test \n%!"; - Alcotest.(check (list (pair int (pair int (pair int (pair int (float 1.e-12))))))) "ERI" a r + Alcotest.(check (list (pair int (pair int (pair int (pair int (float 1.e-06))))))) "ERI" a r in let check_eri_lr title a r = @@ -97,7 +96,7 @@ let test_case name t = in let a = ERI_lr.to_list a |> List.rev_map f |> List.rev in let r = ERI_lr.to_list r |> List.rev_map f |> List.rev in - Alcotest.(check (list (pair int (pair int (pair int (pair int (float 1.e-12))))))) "ERI_lr" a r + Alcotest.(check (list (pair int (pair int (pair int (pair int (float 1.e-06))))))) "ERI_lr" a r in let test_overlap () = diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml index ebc73b6..62b5fd3 100644 --- a/Basis/TwoElectronRR.ml +++ b/Basis/TwoElectronRR.ml @@ -301,7 +301,7 @@ let rec hvrr_two_e -let contracted_class_shell_pair_couple ~zero_m shell_pair_couple : float Zmap.t = +let contracted_class_shell_pair_couple ~zero_m ?mu_erf shell_pair_couple : float Zmap.t = let maxm = Am.to_int (Cspc.ang_mom shell_pair_couple) in @@ -340,12 +340,11 @@ let contracted_class_shell_pair_couple ~zero_m shell_pair_couple : float Zmap.t let norm_pq_sq = Co.dot center_pq center_pq in let expo_p_inv = Psp.exponent_inv sp_ab in let expo_q_inv = Psp.exponent_inv sp_cd in - let normalization = Psp.normalization sp_ab *. Psp.normalization sp_cd in - let zero = Zp.zero zero_m in + let zero = Zp.zero ?mu_erf zero_m in let zero_m_array = zero_m { zero with maxm ; expo_p_inv ; expo_q_inv ; norm_pq_sq ; - center_pq ; center_pa ; center_qc ; normalization ; + center_pq ; center_pa ; center_qc ; } in @@ -412,7 +411,7 @@ let contracted_class_shell_pair_couple ~zero_m shell_pair_couple : float Zmap.t -let contracted_class_atomic_shell_pair_couple ~zero_m atomic_shell_pair_couple : float Zmap.t = +let contracted_class_atomic_shell_pair_couple ~zero_m ?mu_erf atomic_shell_pair_couple : float Zmap.t = let maxm = Am.to_int (Aspc.ang_mom atomic_shell_pair_couple) in @@ -454,13 +453,12 @@ let contracted_class_atomic_shell_pair_couple ~zero_m atomic_shell_pair_couple : let center_pa = Psp.center_minus_a sp_ab in let norm_pq_sq = Co.dot center_pq center_pq in let expo_q_inv = Psp.exponent_inv sp_cd in - let normalization = Psp.normalization sp_ab *. Psp.normalization sp_cd in - let zero = Zp.zero zero_m in + let zero = Zp.zero ?mu_erf zero_m in let zero_m_array = zero_m { zero with maxm ; expo_p_inv ; expo_q_inv ; norm_pq_sq ; - center_pq ; center_pa ; center_qc ; normalization ; + center_pq ; center_pa ; center_qc ; } in diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml index 293f45c..4442f51 100644 --- a/Basis/TwoElectronRRVectorized.ml +++ b/Basis/TwoElectronRRVectorized.ml @@ -24,7 +24,7 @@ let at_least_one_valid arr = Array.iter (fun x -> if (abs_float x > cutoff) then raise Found) arr ; false with Found -> true -type fou_idx_intermediate = +type four_idx_intermediate = { expo_b : float array; expo_d : float array; @@ -577,7 +577,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) -let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q : float Zmap.t = +let contracted_class_shell_pairs ~zero_m ?mu_erf ?schwartz_p ?schwartz_q shell_p shell_q : float Zmap.t = let sp = Csp.shell_pairs shell_p and sq = Csp.shell_pairs shell_q @@ -651,16 +651,13 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q let norm_pq_sq = Co.dot center_pq center_pq in - let normalization = Psp.normalization sp.(i-1) *. - Psp.normalization sq.(i-1) - in - let zero = Zp.zero zero_m in + let zero = Zp.zero ?mu_erf zero_m in let zero_m_array = zero_m {zero with expo_p_inv ; expo_q_inv ; norm_pq_sq ; - center_pq ; center_pa ; center_qc ; normalization ; + center_pq ; center_pa ; center_qc ; } in zero_m_array.(0) @@ -786,15 +783,11 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q let norm_pq_sq = x *. x +. y *. y +. z *. z in - let normalization = Psp.normalization shell_ab *. - Psp.normalization shell_cd - in - let zero = Zp.zero zero_m in + let zero = Zp.zero ?mu_erf zero_m in zero_m {zero with maxm ; expo_p_inv ; expo_q_inv ; norm_pq_sq ; center_pq = Coordinate.make Coordinate.{x ; y ; z} ; center_pa ; center_qc = center_qc_tmp.(cd) ; - normalization ; } ) sq in diff --git a/Basis/Zero_m_parameters.ml b/Basis/Zero_m_parameters.ml index ebff133..b5e1bdf 100644 --- a/Basis/Zero_m_parameters.ml +++ b/Basis/Zero_m_parameters.ml @@ -3,7 +3,6 @@ type t = expo_p_inv : float ; expo_q_inv : float ; norm_pq_sq : float ; - normalization : float ; maxm : int ; center_pq : Coordinate.t ; center_pa : Coordinate.t ; @@ -22,8 +21,9 @@ let zero ?f12_factor ?mu_erf zero_m_func = expo_p_inv = 0.; expo_q_inv = 0.; norm_pq_sq = 0.; - normalization = 0.; center_pq = Coordinate.zero ; center_pa = Coordinate.zero ; center_qc = Coordinate.zero ; } + + diff --git a/Simulation.ml b/Simulation.ml index 03f599d..6a010b4 100644 --- a/Simulation.ml +++ b/Simulation.ml @@ -60,6 +60,7 @@ let of_filenames ?(cartesian=false) ?(multiplicity=1) ?(charge=0) ?f12 ?mu_erf ~ let basis = Basis.of_nuclei_and_basis_filenames ~nuclei (basis_filename :: aux_basis_filenames) in + lazy (make ?mu_erf ?f12 ~cartesian ~charge ~multiplicity ~nuclei basis) |> Parallel.broadcast diff --git a/run_tests.ml b/run_tests.ml index f364b4d..bd64ced 100644 --- a/run_tests.ml +++ b/run_tests.ml @@ -22,12 +22,12 @@ let () = (* LiH / cc-pVTZ *) let test_integrals_lih = - let basis_file = "test_files/cc-pvtz" + let basis_file = "test_files/cc-pvdz" and nuclei_file = "test_files/lih.xyz" in let simulation_closed_shell = - Simulation.of_filenames ~charge:0 ~cartesian:true ~multiplicity:1 ~nuclei:nuclei_file basis_file + Simulation.of_filenames ~mu_erf:0.5 ~charge:0 ~cartesian:true ~multiplicity:1 ~nuclei:nuclei_file basis_file in let ao_basis =