From fe3fcd7a127a235e76737f908d84baa154d5b994 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 12 Feb 2018 00:56:32 +0100 Subject: [PATCH] Schwartz in VectorRR --- Basis/TwoElectronRR.ml | 94 +++++++++++++++----------------- Basis/TwoElectronRRVectorized.ml | 38 +++++++++---- Utils/Zkey.ml | 2 + 3 files changed, 73 insertions(+), 61 deletions(-) diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml index da5f37f..f4ed9da 100644 --- a/Basis/TwoElectronRR.ml +++ b/Basis/TwoElectronRR.ml @@ -350,7 +350,9 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q let monocentric = shell_p.ContractedShellPair.monocentric && - shell_q.ContractedShellPair.monocentric + shell_q.ContractedShellPair.monocentric && + Contracted_shell.center shell_p.ContractedShellPair.shell_a = + Contracted_shell.center shell_q.ContractedShellPair.shell_a in (* Compute all integrals in the shell for each pair of significant shell pairs *) @@ -369,10 +371,6 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q if (abs_float coef_prod) < 1.e-3*.cutoff then raise NullQuartet; - let expo_pq_inv = - shell_p.ContractedShellPair.expo_inv.(ab) +. - shell_q.ContractedShellPair.expo_inv.(cd) - in let center_pq = Coordinate.(sp.(ab).ShellPair.center |- sq.(cd).ShellPair.center) in @@ -380,6 +378,11 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q Coordinate.dot center_pq center_pq in + let expo_pq_inv = + shell_p.ContractedShellPair.expo_inv.(ab) +. + shell_q.ContractedShellPair.expo_inv.(cd) + in + let zero_m_array = zero_m ~maxm ~expo_pq_inv ~norm_pq_sq in @@ -403,14 +406,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q Array.map (fun v2 -> v1 *. v2) norm_coef_scale_q) |> Array.concat in - let () = - if debug then ( - if monocentric then - Printf.printf "Mono-centric\n" - else - Printf.printf "Multi-centric\n" - ) - in + (* Compute the integral class from the primitive shell quartet *) class_indices |> Array.iteri (fun i key -> @@ -420,49 +416,45 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q | _ -> assert false in try - if monocentric then - begin - let ax,ay,az = angMomA - and bx,by,bz = angMomB - and cx,cy,cz = angMomC - and dx,dy,dz = angMomD - in - if ( ((1 land ax+bx+cx+dx)=1) || - ((1 land ay+by+cy+dy)=1) || - ((1 land az+bz+cz+dz)=1) - ) then - raise NullQuartet - end; - (* Schwartz screening *) + if monocentric then + begin + let ax,ay,az = angMomA + and bx,by,bz = angMomB + and cx,cy,cz = angMomC + and dx,dy,dz = angMomD + in + if ( ((1 land ax+bx+cx+dx)=1) || + ((1 land ay+by+cy+dy)=1) || + ((1 land az+bz+cz+dz)=1) + ) then + raise NullQuartet + end; (* - let schwartz_p = - let key = - Zkey.of_int_array Zkey.Kind_12 - [| a.(0) ; a.(1) ; a.(2) ; - a.(3) ; a.(4) ; a.(5) ; - a.(0) ; a.(1) ; a.(2) ; - a.(3) ; a.(4) ; a.(5) |] + (* Schwartz screening *) + if (maxm > 2) then + ( + let schwartz_p = + let key = Zkey.of_int_tuple (Zkey.Twelve + (angMomA, angMomB, angMomA, angMomB) ) + in + match schwartz_p with + | None -> 1. + | Some schwartz_p -> Zmap.find schwartz_p key in - match schwartz_p with - | None -> 1. - | Some schwartz_p -> Zmap.find schwartz_p key - in - if schwartz_p < cutoff then raise NullQuartet; - let schwartz_q = - let key = - Zkey.of_int_array Zkey.Kind_12 - [| a.(6) ; a.(7) ; a.(8) ; - a.(9) ; a.(10) ; a.(11) ; - a.(6) ; a.(7) ; a.(8) ; - a.(9) ; a.(10) ; a.(11) |] + if schwartz_p < cutoff then raise NullQuartet; + let schwartz_q = + let key = Zkey.of_int_tuple (Zkey.Twelve + (angMomC, angMomD, angMomC, angMomD) ) + in + match schwartz_q with + | None -> 1. + | Some schwartz_q -> Zmap.find schwartz_q key in - match schwartz_q with - | None -> 1. - | Some schwartz_q -> Zmap.find schwartz_q key - in - if schwartz_p *. schwartz_q < cutoff2 then raise NullQuartet; + if schwartz_p *. schwartz_q < cutoff2 then raise NullQuartet; + ); *) + let norm = norm_coef_scale.(i) in let coef_prod = coef_prod *. norm in let integral = diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml index e555d14..6b046d4 100644 --- a/Basis/TwoElectronRRVectorized.ml +++ b/Basis/TwoElectronRRVectorized.ml @@ -495,7 +495,9 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q let monocentric = shell_p.ContractedShellPair.monocentric && - shell_q.ContractedShellPair.monocentric + shell_q.ContractedShellPair.monocentric && + Contracted_shell.center shell_p.ContractedShellPair.shell_a = + Contracted_shell.center shell_q.ContractedShellPair.shell_a in (** Screening on the product of coefficients *) @@ -566,6 +568,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q let expo_pq_inv = expo_inv_p.{i} +. expo_inv_q.{j} in + let center_pq = Coordinate.(sp.(i-1).ShellPair.center |- sq.(j-1).ShellPair.center) in @@ -691,15 +694,30 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q ) then raise NullQuartet end; - (* - let a = Zkey.to_int_array Zkey.Kind_12 key in - let (angMomA,angMomB,angMomC,angMomD) = - ( [| a.(0) ; a.(1) ; a.(2) |], - [| a.(3) ; a.(4) ; a.(5) |], - [| a.(6) ; a.(7) ; a.(8) |], - [| a.(9) ; a.(10) ; a.(11) |] ) - in - *) + + (* Schwartz screening *) + if (np+nq> 24) then + ( + let schwartz_p = + let key = Zkey.of_int_tuple (Zkey.Twelve + (angMomA, angMomB, angMomA, angMomB) ) + in + match schwartz_p with + | None -> 1. + | Some schwartz_p -> Zmap.find schwartz_p key + in + if schwartz_p < cutoff then raise NullQuartet; + let schwartz_q = + let key = Zkey.of_int_tuple (Zkey.Twelve + (angMomC, angMomD, angMomC, angMomD) ) + in + match schwartz_q with + | None -> 1. + | Some schwartz_q -> Zmap.find schwartz_q key + in + if schwartz_p *. schwartz_q < cutoff2 then raise NullQuartet; + ); + let integral = hvrr_two_e_vector (angMomA, angMomB, angMomC, angMomD) (Contracted_shell.totAngMom shell_a, Contracted_shell.totAngMom shell_b, diff --git a/Utils/Zkey.ml b/Utils/Zkey.ml index eb3cb62..49eeae0 100644 --- a/Utils/Zkey.ml +++ b/Utils/Zkey.ml @@ -187,6 +187,8 @@ let to_int_tuple ~kind a = | Kind_1 -> One ( Z.to_int a ) +let zero = Z.of_int 0 + include Z let to_string ~kind a =