LCPQ/QCaml
10
1
Fork 0
mirror of https://gitlab.com/scemama/QCaml.git synced 2024-08-29 07:51:59 +02:00
Code Releases Activity

Works

Browse Source
This commit is contained in:
Anthony Scemama 2018-01-29 15:29:38 +01:00
parent 2d759e736e
commit d4242e406e
1 changed files with 78 additions and 74 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

152
Basis/TwoElectronRRVectorized.ml
View File

@ -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) -> | Angular_momentum.(S,S,S,S) ->
contracted_class.(0) <- contracted_class.(0) <-
Array.fold_left (fun accu shell_ab -> accu +. Array.fold_left (fun accu shell_ab -> accu +.
Array.fold_left (fun accu shell_cd -> Array.fold_left (fun accu shell_cd ->
let coef_prod = let coef_prod =
shell_ab.Shell_pair.coef *. shell_cd.Shell_pair.coef shell_ab.Shell_pair.coef *. shell_cd.Shell_pair.coef
in in
(** Screening on the product of coefficients *) (** Screening on the product of coefficients *)
try try
if (abs_float coef_prod) < 1.e-4*.cutoff then if (abs_float coef_prod) < 1.e-3*.cutoff then
raise NullQuartet; raise NullQuartet;
let expo_pq_inv = let expo_pq_inv =
shell_ab.Shell_pair.expo_inv +. shell_cd.Shell_pair.expo_inv shell_ab.Shell_pair.expo_inv +. shell_cd.Shell_pair.expo_inv
in in
let center_pq = let center_pq =
Coordinate.(shell_ab.Shell_pair.center |- shell_cd.Shell_pair.center) Coordinate.(shell_ab.Shell_pair.center |- shell_cd.Shell_pair.center)
in in
let norm_pq_sq = let norm_pq_sq =
Coordinate.dot center_pq center_pq Coordinate.dot center_pq center_pq
in in
let zero_m_array = let zero_m_array =
zero_m ~maxm ~expo_pq_inv ~norm_pq_sq zero_m ~maxm ~expo_pq_inv ~norm_pq_sq
in in
let coef_prod = let coef_prod =
shell_ab.Shell_pair.coef *. shell_cd.Shell_pair.coef shell_ab.Shell_pair.coef *. shell_cd.Shell_pair.coef
in in
let integral = let integral =
zero_m_array.(0) zero_m_array.(0)
in in
accu +. coef_prod *. integral accu +. coef_prod *. integral
with NullQuartet -> accu with NullQuartet -> accu
) 0. shell_q ) 0. shell_q
) 0. shell_p ) 0. shell_p
| _ -> | _ ->
Array.iter (fun shell_ab -> Array.iter (fun shell_ab ->
let b = shell_ab.Shell_pair.j in let b = shell_ab.Shell_pair.j in
let common = let common =
Array.map (fun shell_cd -> Array.mapi (fun idx shell_cd ->
let coef_prod = let coef_prod =
shell_ab.Shell_pair.coef *. shell_cd.Shell_pair.coef shell_ab.Shell_pair.coef *. shell_cd.Shell_pair.coef
in in
let coef_prod = let expo_pq_inv =
if (abs_float coef_prod) < 1.e-4*.cutoff then shell_ab.Shell_pair.expo_inv +. shell_cd.Shell_pair.expo_inv
0. else coef_prod in
in let center_pq =
let expo_pq_inv = Coordinate.(shell_ab.Shell_pair.center |- shell_cd.Shell_pair.center)
shell_ab.Shell_pair.expo_inv +. shell_cd.Shell_pair.expo_inv in
in let norm_pq_sq =
let center_pq = Coordinate.dot center_pq center_pq
Coordinate.(shell_ab.Shell_pair.center |- shell_cd.Shell_pair.center) in
in
let norm_pq_sq =
Coordinate.dot center_pq center_pq
in
let zero_m_array = let zero_m_array =
zero_m ~maxm ~expo_pq_inv ~norm_pq_sq 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 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, 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, 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, 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, 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, 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, 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 in
let map = Array.init maxm (fun _ -> Zmap.create (Array.length class_indices)) in let map = Array.init maxm (fun _ -> Zmap.create (Array.length class_indices)) in
(* Compute the integral class from the primitive shell quartet *) (* Compute the integral class from the primitive shell quartet *)
@ -330,24 +332,26 @@ Array.iter (fun shell_ab ->
in in
let norm = let norm =
Array.map (fun shell_cd -> 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 ) shell_q
in in
let integral = try
let integral =
hvrr_two_e 0 (angMomA, angMomB, angMomC, angMomD) hvrr_two_e 0 (angMomA, angMomB, angMomC, angMomD)
(Contracted_shell.totAngMom shell_a, Contracted_shell.totAngMom shell_b, (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) (maxm, zero_m_array)
(Contracted_shell.expo shell_b b, d) (Contracted_shell.expo shell_b b, d)
(shell_ab.Shell_pair.expo_inv, expo_inv) (shell_ab.Shell_pair.expo_inv, expo_inv)
(shell_ab.Shell_pair.center_ab, center_cd, center_pq) (shell_ab.Shell_pair.center_ab, center_cd, center_pq)
coef_prod map coef_prod map
|> Array.map2 (fun x y -> x *. y ) norm |> Array.mapi (fun i x -> x *. norm.(idx.(i)) )
in in
let x = Array.fold_left (+.) 0. integral in let x = Array.fold_left (+.) 0. integral in
contracted_class.(i) <- contracted_class.(i) +. x contracted_class.(i) <- contracted_class.(i) +. x
with NullQuartet -> ()
) class_indices ) class_indices
) shell_p ) shell_p
end; end;