From d4242e406e303bfaaadaca35ef9a83b28805a216 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 29 Jan 2018 15:29:38 +0100 Subject: [PATCH] Works --- Basis/TwoElectronRRVectorized.ml | 152 ++++++++++++++++--------------- 1 file changed, 78 insertions(+), 74 deletions(-) diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml index 0a7595a..84b66d9 100644 --- a/Basis/TwoElectronRRVectorized.ml +++ b/Basis/TwoElectronRRVectorized.ml @@ -238,85 +238,87 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q | Angular_momentum.(S,S,S,S) -> contracted_class.(0) <- Array.fold_left (fun accu shell_ab -> accu +. - Array.fold_left (fun accu shell_cd -> - let coef_prod = - shell_ab.Shell_pair.coef *. shell_cd.Shell_pair.coef - in - (** Screening on the product of coefficients *) - try - if (abs_float coef_prod) < 1.e-4*.cutoff then - raise NullQuartet; + Array.fold_left (fun accu shell_cd -> + let coef_prod = + shell_ab.Shell_pair.coef *. shell_cd.Shell_pair.coef + in + (** Screening on the product of coefficients *) + try + if (abs_float coef_prod) < 1.e-3*.cutoff then + raise NullQuartet; - let expo_pq_inv = - shell_ab.Shell_pair.expo_inv +. shell_cd.Shell_pair.expo_inv - in - let center_pq = - Coordinate.(shell_ab.Shell_pair.center |- shell_cd.Shell_pair.center) - in - let norm_pq_sq = - Coordinate.dot center_pq center_pq - in + let expo_pq_inv = + shell_ab.Shell_pair.expo_inv +. shell_cd.Shell_pair.expo_inv + in + let center_pq = + Coordinate.(shell_ab.Shell_pair.center |- shell_cd.Shell_pair.center) + in + let norm_pq_sq = + Coordinate.dot center_pq center_pq + in - let zero_m_array = - zero_m ~maxm ~expo_pq_inv ~norm_pq_sq - in + let zero_m_array = + zero_m ~maxm ~expo_pq_inv ~norm_pq_sq + in - let coef_prod = - shell_ab.Shell_pair.coef *. shell_cd.Shell_pair.coef - in - let integral = - zero_m_array.(0) - in - accu +. coef_prod *. integral - with NullQuartet -> accu - ) 0. shell_q - ) 0. shell_p + let coef_prod = + shell_ab.Shell_pair.coef *. shell_cd.Shell_pair.coef + in + let integral = + zero_m_array.(0) + in + accu +. coef_prod *. integral + with NullQuartet -> accu + ) 0. shell_q + ) 0. shell_p | _ -> -Array.iter (fun shell_ab -> - let b = shell_ab.Shell_pair.j in - let common = - Array.map (fun shell_cd -> - let coef_prod = - shell_ab.Shell_pair.coef *. shell_cd.Shell_pair.coef - in - let coef_prod = - if (abs_float coef_prod) < 1.e-4*.cutoff then - 0. else coef_prod - in - let expo_pq_inv = - shell_ab.Shell_pair.expo_inv +. shell_cd.Shell_pair.expo_inv - in - let center_pq = - Coordinate.(shell_ab.Shell_pair.center |- shell_cd.Shell_pair.center) - in - let norm_pq_sq = - Coordinate.dot center_pq center_pq - in + Array.iter (fun shell_ab -> + let b = shell_ab.Shell_pair.j in + let common = + Array.mapi (fun idx shell_cd -> + let coef_prod = + shell_ab.Shell_pair.coef *. shell_cd.Shell_pair.coef + in + let expo_pq_inv = + shell_ab.Shell_pair.expo_inv +. shell_cd.Shell_pair.expo_inv + in + let center_pq = + Coordinate.(shell_ab.Shell_pair.center |- shell_cd.Shell_pair.center) + in + let norm_pq_sq = + Coordinate.dot center_pq center_pq + in - let zero_m_array = - zero_m ~maxm ~expo_pq_inv ~norm_pq_sq + let zero_m_array = + zero_m ~maxm ~expo_pq_inv ~norm_pq_sq + in + + let d = shell_cd.Shell_pair.j in + + (zero_m_array, shell_cd.Shell_pair.expo_inv, + Contracted_shell.expo shell_d d, shell_cd.Shell_pair.center_ab, + center_pq,coef_prod,idx) + ) shell_q + |> Array.to_list + |> List.filter (fun (zero_m_array, expo_inv, d, center_cd, + center_pq,coef_prod,idx) -> abs_float coef_prod >= 1.e-4 *. cutoff) + |> Array.of_list in - - let d = shell_cd.Shell_pair.j in - - (zero_m_array, shell_cd.Shell_pair.expo_inv, - Contracted_shell.expo shell_d d, shell_cd.Shell_pair.center_ab, - center_pq,coef_prod) - ) shell_q - in let zero_m_array = Array.map (fun (zero_m_array, expo_inv, d, center_cd, - center_pq,coef_prod) -> zero_m_array) common + center_pq,coef_prod,idx) -> zero_m_array) common and expo_inv = Array.map (fun (zero_m_array, expo_inv, d, center_cd, - center_pq,coef_prod) -> expo_inv ) common + center_pq,coef_prod,idx) -> expo_inv ) common and d = Array.map (fun (zero_m_array, expo_inv, d, center_cd, - center_pq,coef_prod) -> d) common + center_pq,coef_prod,idx) -> d) common and center_cd = Array.map (fun (zero_m_array, expo_inv, d, center_cd, - center_pq,coef_prod) -> center_cd) common + center_pq,coef_prod,idx) -> center_cd) common and center_pq = Array.map (fun (zero_m_array, expo_inv, d, center_cd, - center_pq,coef_prod) -> center_pq) common + center_pq,coef_prod,idx) -> center_pq) common and coef_prod = Array.map (fun (zero_m_array, expo_inv, d, center_cd, - center_pq,coef_prod) -> coef_prod) common + center_pq,coef_prod,idx) -> coef_prod) common + and idx = Array.map (fun (zero_m_array, expo_inv, d, center_cd, + center_pq,coef_prod,idx) -> idx) common in let map = Array.init maxm (fun _ -> Zmap.create (Array.length class_indices)) in (* Compute the integral class from the primitive shell quartet *) @@ -330,24 +332,26 @@ Array.iter (fun shell_ab -> in let norm = Array.map (fun shell_cd -> - shell_ab.Shell_pair.norm_fun angMomA angMomB *. shell_cd.Shell_pair.norm_fun angMomC angMomD + shell_ab.Shell_pair.norm_fun angMomA angMomB *. shell_cd.Shell_pair.norm_fun angMomC angMomD ) shell_q in - let integral = + try + let integral = hvrr_two_e 0 (angMomA, angMomB, angMomC, angMomD) (Contracted_shell.totAngMom shell_a, Contracted_shell.totAngMom shell_b, - Contracted_shell.totAngMom shell_c, Contracted_shell.totAngMom shell_d) + Contracted_shell.totAngMom shell_c, Contracted_shell.totAngMom shell_d) (maxm, zero_m_array) (Contracted_shell.expo shell_b b, d) (shell_ab.Shell_pair.expo_inv, expo_inv) (shell_ab.Shell_pair.center_ab, center_cd, center_pq) coef_prod map - |> Array.map2 (fun x y -> x *. y ) norm - in - let x = Array.fold_left (+.) 0. integral in - contracted_class.(i) <- contracted_class.(i) +. x + |> Array.mapi (fun i x -> x *. norm.(idx.(i)) ) + in + let x = Array.fold_left (+.) 0. integral in + contracted_class.(i) <- contracted_class.(i) +. x + with NullQuartet -> () ) class_indices - ) shell_p + ) shell_p end;