LCPQ/QCaml
10
1
Fork 0
mirror of https://gitlab.com/scemama/QCaml.git synced 2024-06-26 15:12:05 +02:00
Code Releases Activity

Working on contraction

Browse Source
This commit is contained in:
Anthony Scemama 2018-02-09 19:41:22 +01:00
parent 5420d41503
commit 365735d111
8 changed files with 366 additions and 296 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

132
Basis/ContractedShellPair.ml
View File

@ -6,7 +6,17 @@ exception Null_contribution
(** Array of shell pairs obtained from combining two contracted shells. The
two integers correspond to the indices of the combined shells.
*)
type t = ShellPair.t array
type t =
{
shell_a : Contracted_shell.t;
shell_b : Contracted_shell.t;
shell_pairs : ShellPair.t array;
coef : float array;
expo_inv : float array;
center_ab : Coordinate.t; (* A-B *)
norm_sq : float; (* |A-B|^2 *)
norm_coef_scale : float array; (* norm_coef.(i) / norm_coef.(0) *)
}
(** Creates an contracted shell pair : an array of pairs of primitive shells.
@ -38,62 +48,72 @@ let create ?cutoff p_a p_b =
|> Array.to_list
|> Array.concat
in
Array.init (Contracted_shell.size p_a) (fun i ->
let p_a_expo_center = Coordinate.(
Contracted_shell.expo p_a i |. Contracted_shell.center p_a )
in
let norm_coef_a =
Contracted_shell.norm_coef p_a i
in
let shell_pairs =
Array.init (Contracted_shell.size p_a) (fun i ->
let p_a_expo_center = Coordinate.(
Contracted_shell.expo p_a i |. Contracted_shell.center p_a )
in
let norm_coef_a =
Contracted_shell.norm_coef p_a i
in
Array.init (Contracted_shell.size p_b) (fun j ->
try
let norm_coef_b =
Contracted_shell.norm_coef p_b j
in
let norm_coef =
norm_coef_a *. norm_coef_b
in
if (norm_coef < cutoff) then
raise Null_contribution;
let p_b_expo_center = Coordinate.(
Contracted_shell.expo p_b j |. Contracted_shell.center p_b )
in
let expo = Contracted_shell.(expo p_a i +. expo p_b j) in
let expo_inv = 1. /. expo in
let center = Coordinate.(
expo_inv |. (p_a_expo_center |+ p_b_expo_center ) )
in
let argexpo =
Contracted_shell.(expo p_a i *. expo p_b j) *. norm_sq *. expo_inv
in
if (argexpo > log_cutoff) then
raise Null_contribution;
let g =
(pi *. expo_inv)**(1.5) *. exp(-. argexpo)
in
let coef =
norm_coef *. Contracted_shell.(coef p_a i *. coef p_b j) *. g
in
if (abs_float coef < cutoff) then
raise Null_contribution;
let center_a =
Coordinate.(center |- Contracted_shell.center p_a)
in
let monocentric =
Contracted_shell.center p_a = Contracted_shell.center p_b
in
Some ShellPair.{ i ; j ; shell_a=p_a ; shell_b=p_b ; norm_coef ; coef ; expo ; expo_inv ; center ; center_a ; center_ab ; norm_sq ; norm_coef_scale ; monocentric }
with
| Null_contribution -> None
)
)
|> Array.to_list
|> Array.concat
|> Array.to_list
|> List.filter (function Some _ -> true | None -> false)
|> List.map (function Some x -> x | None -> assert false)
|> Array.of_list
Array.init (Contracted_shell.size p_b) (fun j ->
try
let norm_coef_b =
Contracted_shell.norm_coef p_b j
in
let norm_coef =
norm_coef_a *. norm_coef_b
in
if (norm_coef < cutoff) then
raise Null_contribution;
let p_b_expo_center = Coordinate.(
Contracted_shell.expo p_b j |. Contracted_shell.center p_b )
in
let expo = Contracted_shell.(expo p_a i +. expo p_b j) in
let expo_inv = 1. /. expo in
let center = Coordinate.(
expo_inv |. (p_a_expo_center |+ p_b_expo_center ) )
in
let argexpo =
Contracted_shell.(expo p_a i *. expo p_b j) *. norm_sq *. expo_inv
in
if (argexpo > log_cutoff) then
raise Null_contribution;
let g =
(pi *. expo_inv)**(1.5) *. exp(-. argexpo)
in
let coef =
norm_coef *. Contracted_shell.(coef p_a i *. coef p_b j) *. g
in
if (abs_float coef < cutoff) then
raise Null_contribution;
let center_a =
Coordinate.(center |- Contracted_shell.center p_a)
in
let monocentric =
Contracted_shell.center p_a = Contracted_shell.center p_b
in
Some ShellPair.{ i ; j ; shell_a=p_a ; shell_b=p_b ; norm_coef ; coef ; expo ; expo_inv ; center ; center_a ; center_ab ; norm_sq ; monocentric }
with
| Null_contribution -> None
)
)
|> Array.to_list
|> Array.concat
|> Array.to_list
|> List.filter (function Some _ -> true | None -> false)
|> List.map (function Some x -> x | None -> assert false)
|> Array.of_list
in
let coef = Array.map (fun x -> (fun y -> y.ShellPair.coef) x) shell_pairs
and expo_inv = Array.map (fun x -> (fun y -> y.ShellPair.expo_inv) x) shell_pairs
in
{
shell_a = p_a ; shell_b = p_b ; coef ; expo_inv ;
shell_pairs ; center_ab=shell_pairs.(0).center_ab;
norm_coef_scale ; norm_sq=shell_pairs.(0).norm_sq
}
(** Returns an integer characteristic of a contracted shell pair *)

13
Basis/ERI.ml
View File

@ -134,6 +134,10 @@ let of_basis basis =
shell_p = shell_pairs.(i).(j)
in
let sp =
shell_p.ContractedShellPair.shell_pairs
in
for k=0 to i do
for l=0 to k do
let schwartz_q, schwartz_q_max = schwartz.(k).(l) in
@ -143,20 +147,23 @@ let of_basis basis =
let
shell_q = shell_pairs.(k).(l)
in
let sq =
shell_q.ContractedShellPair.shell_pairs
in
let swap =
Array.length shell_q < Array.length shell_p
Array.length sp > Array.length sq
in
(* Compute all the integrals of the class *)
let cls =
if swap then
if Array.length shell_p < 2 then
if Array.length sp < 2 then
contracted_class_shell_pairs ~schwartz_p:schwartz_q ~schwartz_q:schwartz_p shell_q shell_p
else
contracted_class_shell_pairs_vec ~schwartz_p:schwartz_q ~schwartz_q:schwartz_p shell_q shell_p
else
if Array.length shell_q < 2 then
if Array.length sq < 2 then
contracted_class_shell_pairs ~schwartz_p ~schwartz_q shell_p shell_q
else
contracted_class_shell_pairs_vec ~schwartz_p ~schwartz_q shell_p shell_q

28
Basis/KinInt.ml
View File

@ -24,33 +24,35 @@ let contracted_class shell_a shell_b : float Zmap.t =
(* Compute all integrals in the shell for each pair of significant shell pairs *)
for ab=0 to (Array.length shell_p - 1)
let sp = shell_p.ContractedShellPair.shell_pairs in
let center_ab =
shell_p.ContractedShellPair.center_ab
in
let norm_coef_scale =
shell_p.ContractedShellPair.norm_coef_scale
in
for ab=0 to (Array.length sp - 1)
do
let coef_prod =
shell_p.(ab).ShellPair.coef
shell_p.ContractedShellPair.coef.(ab)
in
(** Screening on thr product of coefficients *)
if (abs_float coef_prod) > 1.e-4*.cutoff then
begin
let center_ab =
shell_p.(ab).ShellPair.center_ab
in
let center_pa =
shell_p.(ab).ShellPair.center_a
sp.(ab).ShellPair.center_a
in
let expo_inv =
shell_p.(ab).ShellPair.expo_inv
in
let norm_coef_scale =
shell_p.(ab).ShellPair.norm_coef_scale
shell_p.ContractedShellPair.expo_inv.(ab)
in
let i, j =
shell_p.(ab).ShellPair.i, shell_p.(ab).ShellPair.j
sp.(ab).ShellPair.i, sp.(ab).ShellPair.j
in
let expo_a =
Contracted_shell.expo shell_p.(ab).ShellPair.shell_a i
Contracted_shell.expo sp.(ab).ShellPair.shell_a i
and expo_b =
Contracted_shell.expo shell_p.(ab).ShellPair.shell_b j
Contracted_shell.expo sp.(ab).ShellPair.shell_b j
in
Array.iteri (fun i key ->

21
Basis/OneElectronRR.ml
View File

@ -123,8 +123,8 @@ let hvrr_one_e
(** Computes all the one-electron integrals of the contracted shell pair *)
let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t =
let shell_a = shell_p.(0).ShellPair.shell_a
and shell_b = shell_p.(0).ShellPair.shell_b
let shell_a = shell_p.ContractedShellPair.shell_a
and shell_b = shell_p.ContractedShellPair.shell_b
in
let maxm =
let open Angular_momentum in
@ -144,13 +144,14 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t =
(* Compute all integrals in the shell for each pair of significant shell pairs *)
for ab=0 to (Array.length shell_p - 1)
let norm_coef_scale_p = shell_p.ContractedShellPair.norm_coef_scale
in
for ab=0 to (Array.length shell_p.ContractedShellPair.shell_pairs - 1)
do
let b = shell_p.(ab).ShellPair.j in
let b = shell_p.ContractedShellPair.shell_pairs.(ab).ShellPair.j in
try
begin
let coef_prod = shell_p.(ab).ShellPair.coef in
let norm_coef_scale_p = shell_p.(ab).ShellPair.norm_coef_scale in
let coef_prod = shell_p.ContractedShellPair.coef.(ab) in
(** Screening on the product of coefficients *)
if (abs_float coef_prod) < 1.e-4*.cutoff then
@ -158,14 +159,14 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t =
let expo_pq_inv =
shell_p.(ab).ShellPair.expo_inv
shell_p.ContractedShellPair.expo_inv.(ab)
in
let center_ab =
shell_p.(ab).ShellPair.center_ab
shell_p.ContractedShellPair.center_ab
in
let center_p =
shell_p.(ab).ShellPair.center
shell_p.ContractedShellPair.shell_pairs.(ab).ShellPair.center
in
let center_pa =
Coordinate.(center_p |- Contracted_shell.center shell_a)
@ -208,7 +209,7 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t =
(Contracted_shell.totAngMom shell_a, Contracted_shell.totAngMom shell_b)
(maxm, zero_m_array)
(Contracted_shell.expo shell_b b)
(shell_p.(ab).ShellPair.expo_inv)
(shell_p.ContractedShellPair.expo_inv.(ab))
(center_ab, center_pa, center_pc)
map
in

26
Basis/Overlap.ml
View File

@ -22,27 +22,31 @@ let contracted_class shell_a shell_b : float Zmap.t =
Array.make (Array.length class_indices) 0.
in
let sp =
shell_p.ContractedShellPair.shell_pairs
in
let center_ab =
shell_p.ContractedShellPair.center_ab
in
let norm_coef_scale =
shell_p.ContractedShellPair.norm_coef_scale
in
(* Compute all integrals in the shell for each pair of significant shell pairs *)
for ab=0 to (Array.length shell_p - 1)
for ab=0 to (Array.length sp - 1)
do
let coef_prod =
shell_p.(ab).ShellPair.coef
shell_p.ContractedShellPair.coef.(ab)
in
(** Screening on thr product of coefficients *)
if (abs_float coef_prod) > 1.e-4*.cutoff then
begin
let center_ab =
shell_p.(ab).ShellPair.center_ab
let expo_inv =
shell_p.ContractedShellPair.expo_inv.(ab)
in
let center_pa =
shell_p.(ab).ShellPair.center_a
in
let expo_inv =
shell_p.(ab).ShellPair.expo_inv
in
let norm_coef_scale =
shell_p.(ab).ShellPair.norm_coef_scale
sp.(ab).ShellPair.center_a
in
Array.iteri (fun i key ->

1
Basis/ShellPair.ml
View File

@ -11,7 +11,6 @@ type t = {
norm_coef: float; (* norm_coef_a * norm_coef_b *)
coef : float; (* norm_coef * coef_a * coef_b * g, with
g = (pi/(alpha+beta))^(3/2) exp (-|A-B|^2 * alpha*beta/(alpha+beta)) *)
norm_coef_scale : float array; (* norm_coef.(i) / norm_coef.(0) *)
i : int;
j : int;
shell_a : Contracted_shell.t;

36
Basis/TwoElectronRR.ml
View File

@ -310,10 +310,12 @@ let rec hvrr_two_e (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 shell_a = shell_p.(0).ShellPair.shell_a
and shell_b = shell_p.(0).ShellPair.shell_b
and shell_c = shell_q.(0).ShellPair.shell_a
and shell_d = shell_q.(0).ShellPair.shell_b
let shell_a = shell_p.ContractedShellPair.shell_a
and shell_b = shell_p.ContractedShellPair.shell_b
and shell_c = shell_q.ContractedShellPair.shell_a
and shell_d = shell_q.ContractedShellPair.shell_b
and sp = shell_p.ContractedShellPair.shell_pairs
and sq = shell_q.ContractedShellPair.shell_pairs
in
let maxm =
let open Angular_momentum in
@ -335,13 +337,14 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
(* Compute all integrals in the shell for each pair of significant shell pairs *)
for ab=0 to (Array.length shell_p - 1) do
let cab = shell_p.(ab).ShellPair.coef in
let b = shell_p.(ab).ShellPair.j in
let norm_coef_scale_p = shell_p.(ab).ShellPair.norm_coef_scale in
for cd=0 to (Array.length shell_q - 1) do
let norm_coef_scale_p = shell_p.ContractedShellPair.norm_coef_scale in
let norm_coef_scale_q = shell_q.ContractedShellPair.norm_coef_scale in
for ab=0 to (Array.length sp - 1) do
let cab = shell_p.ContractedShellPair.coef.(ab) in
let b = sp.(ab).ShellPair.j in
for cd=0 to (Array.length shell_q.ContractedShellPair.shell_pairs - 1) do
let coef_prod =
cab *. shell_q.(cd).ShellPair.coef
cab *. shell_q.ContractedShellPair.coef.(cd)
in
(** Screening on the product of coefficients *)
try
@ -349,10 +352,11 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
raise NullQuartet;
let expo_pq_inv =
shell_p.(ab).ShellPair.expo_inv +. shell_q.(cd).ShellPair.expo_inv
shell_p.ContractedShellPair.expo_inv.(ab) +.
shell_q.ContractedShellPair.expo_inv.(cd)
in
let center_pq =
Coordinate.(shell_p.(ab).ShellPair.center |- shell_q.(cd).ShellPair.center)
Coordinate.(sp.(ab).ShellPair.center |- sq.(cd).ShellPair.center)
in
let norm_pq_sq =
Coordinate.dot center_pq center_pq
@ -370,12 +374,11 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
in
contracted_class.(0) <- contracted_class.(0) +. coef_prod *. integral
| _ ->
let d = shell_q.(cd).ShellPair.j in
let d = shell_q.ContractedShellPair.shell_pairs.(cd).ShellPair.j in
let map_1d = Zmap.create (4*maxm) in
let map_2d = Zmap.create (Array.length class_indices) in
let map_1d' = Zmap.create (4*maxm) in
let map_2d' = Zmap.create (Array.length class_indices) in
let norm_coef_scale_q = shell_q.(cd).ShellPair.norm_coef_scale in
let norm_coef_scale =
Array.map (fun v1 ->
Array.map (fun v2 -> v1 *. v2) norm_coef_scale_q
@ -446,8 +449,9 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
Contracted_shell.totAngMom shell_c, Contracted_shell.totAngMom shell_d)
(maxm, zero_m_array)
(Contracted_shell.expo shell_b b, Contracted_shell.expo shell_d d)
(shell_p.(ab).ShellPair.expo_inv, shell_q.(cd).ShellPair.expo_inv)
(shell_p.(ab).ShellPair.center_ab, shell_q.(cd).ShellPair.center_ab, center_pq)
(shell_p.ContractedShellPair.expo_inv.(ab),
shell_q.ContractedShellPair.expo_inv.(cd) )
(sp.(ab).ShellPair.center_ab, sq.(cd).ShellPair.center_ab, center_pq)
map_1d map_2d map_1d' map_2d'
in
contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral

405
Basis/TwoElectronRRVectorized.ml
View File

@ -1,4 +1,6 @@
open Util
open Lacaml.D
open Bigarray
let cutoff = Constants.cutoff
let cutoff2 = cutoff *. cutoff
@ -8,9 +10,13 @@ exception Found
let at_least_one_valid arr =
try
Array.fold_left (fun _ x -> if (abs_float x > cutoff) then raise Found else false ) false arr
Vec.iter (fun x -> if (abs_float x > cutoff) then raise Found) arr ; false
with Found -> true
(*TODO : REMOVE *)
let sum integral =
Array.fold_left (+.) 0. integral
(** Horizontal and Vertical Recurrence Relations (HVRR) *)
let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
(totAngMom_a, totAngMom_b, totAngMom_c, totAngMom_d)
@ -18,12 +24,14 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
(expo_b, expo_d)
(expo_inv_p, expo_inv_q)
(center_ab, center_cd, center_pq)
coef_prod map_1d map_2d
coef_prod map_1d map_2d
=
let ncoef = (Array.length coef_prod) in
let nq = Mat.dim1 coef_prod in
let np = Mat.dim2 coef_prod in
let empty =
Array.make ncoef 0.
Array.make nq 0.
in
let totAngMom_a = Angular_momentum.to_int totAngMom_a
@ -33,7 +41,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
in
(** Vertical recurrence relations *)
let rec vrr0_v m angMom_a = function
let rec vrr0_v l m angMom_a = function
| 1 ->
let xyz =
match angMom_a with
@ -41,16 +49,16 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
| (_,1,_) -> 1
| _ -> 2
in
let f = expo_b *. (Coordinate.coord center_ab xyz) in
Array.init ncoef (fun k -> coef_prod.(k) *. expo_inv_p *.
( (Coordinate.coord center_pq.(k) xyz) *. zero_m_array.(m+1).(k)
-. f *. zero_m_array.(m).(k) ) )
| 0 -> Array.map2 ( *. ) zero_m_array.(m) coef_prod
let f = expo_b.{l} *. (Coordinate.coord center_ab xyz) in
Array.init nq (fun k -> coef_prod.{k+1,l} *. expo_inv_p.{l} *.
(center_pq.{xyz+1,k+1,l} *. zero_m_array.(m+1).{k+1,l}
-. f *. zero_m_array.(m).{k+1,l} ) )
| 0 -> Array.init nq (fun k -> zero_m_array.(m).{k+1,l} *. coef_prod.{k+1,l})
| totAngMom_a ->
let key = Zkey.of_int_tuple (Zkey.Three angMom_a)
in
try Zmap.find map_1d.(m) key with
try Zmap.find map_1d.(m).(l-1) key with
| Not_found ->
let result =
let am, amm, amxyz, xyz =
@ -62,43 +70,43 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
if amxyz < 0 then empty else
let v1 =
let f =
-. expo_b *. expo_inv_p *. (Coordinate.coord center_ab xyz)
-. expo_b.{l} *. expo_inv_p.{l} *. (Coordinate.coord center_ab xyz)
in
if (abs_float f < cutoff) then empty else
Array.map (fun v1k -> f *. v1k) (vrr0_v m am (totAngMom_a-1) )
Array.map (fun v1k -> f *. v1k) (vrr0_v l m am (totAngMom_a-1) )
in
let p1 =
Array.mapi (fun k v2k -> v1.(k) +. expo_inv_p *. (Coordinate.coord center_pq.(k) xyz) *. v2k) (vrr0_v (m+1) am (totAngMom_a-1))
Array.mapi (fun k v2k -> v1.(k) +. expo_inv_p.{l} *. center_pq.{xyz+1,k+1,l} *. v2k) (vrr0_v l (m+1) am (totAngMom_a-1))
in
if amxyz < 1 then p1 else
let f = (float_of_int amxyz) *. expo_inv_p *. 0.5
let f = (float_of_int amxyz) *. expo_inv_p.{l} *. 0.5
in
if (abs_float f < cutoff) then empty else
let v1 = vrr0_v m amm (totAngMom_a-2)
let v1 = vrr0_v l m amm (totAngMom_a-2)
in
let v2 =
if (abs_float (f *. expo_inv_p)) < cutoff then empty else
vrr0_v (m+1) amm (totAngMom_a-2)
if (abs_float (f *. expo_inv_p.{l})) < cutoff then empty else
vrr0_v l (m+1) amm (totAngMom_a-2)
in
Array.init ncoef (fun k -> p1.(k) +.
f *. (v1.(k) +. v2.(k) *. expo_inv_p ) )
in Zmap.add map_1d.(m) key result;
Array.init nq (fun k -> p1.(k) +.
f *. (v1.(k) +. v2.(k) *. expo_inv_p.{l} ) )
in Zmap.add map_1d.(m).(l-1) key result;
result
and vrr_v m angMom_a angMom_c totAngMom_a totAngMom_c =
and vrr_v l m angMom_a angMom_c totAngMom_a totAngMom_c =
match (totAngMom_a, totAngMom_c) with
| (i,0) -> if (i>0) then
vrr0_v m angMom_a totAngMom_a
vrr0_v l m angMom_a totAngMom_a
else
Array.map2 ( *. ) zero_m_array.(m) coef_prod
Array.init nq (fun k -> zero_m_array.(m).{k+1,l} *. coef_prod.{k+1,l})
| (_,_) ->
let key = Zkey.of_int_tuple (Zkey.Six (angMom_a, angMom_c))
in
try Zmap.find map_2d.(m) key with
try Zmap.find map_2d.(m).(l-1) key with
| Not_found ->
let result =
begin
@ -126,24 +134,26 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
*)
let f1 =
let f = (Coordinate.coord center_cd xyz) in
Array.init ncoef (fun k ->
expo_d.(k) *. expo_inv_q.(k) *. f)
Array.init nq (fun k ->
expo_d.{k+1} *. expo_inv_q.{k+1} *. f)
|> Vec.of_array
in
let f2 =
Array.init ncoef (fun k ->
expo_inv_q.(k) *. (Coordinate.coord center_pq.(k) xyz) )
Array.init nq (fun k ->
expo_inv_q.{k+1} *. center_pq.{xyz+1,k+1,l} )
|> Vec.of_array
in
let v1 =
if (at_least_one_valid f1) then
vrr_v m angMom_a cm totAngMom_a (totAngMom_c-1)
vrr_v l m angMom_a cm totAngMom_a (totAngMom_c-1)
else empty
and v2 =
if (at_least_one_valid f2) then
vrr_v (m+1) angMom_a cm totAngMom_a (totAngMom_c-1)
vrr_v l (m+1) angMom_a cm totAngMom_a (totAngMom_c-1)
else empty
in
let p1 =
Array.init ncoef (fun k -> -. v1.(k) *. f1.(k) -. v2.(k) *. f2.(k))
Array.init nq (fun k -> -. v1.(k) *. f1.{k+1} -. v2.(k) *. f2.{k+1})
in
let p2 =
if cxyz < 2 then p1 else
@ -151,54 +161,56 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
(float_of_int (cxyz-1)) *. 0.5
in
let f1 =
Array.map (fun e -> fcm *. e) expo_inv_q
Vec.map (fun e -> fcm *. e) expo_inv_q
in
let f2 =
Array.map2 ( *. ) f1 expo_inv_q
Vec.mul f1 expo_inv_q
in
let v1 =
if (at_least_one_valid f1) then
vrr_v m angMom_a cmm totAngMom_a (totAngMom_c-2)
vrr_v l m angMom_a cmm totAngMom_a (totAngMom_c-2)
else empty
in
let v2 =
if (at_least_one_valid f2) then
vrr_v (m+1) angMom_a cmm totAngMom_a (totAngMom_c-2)
vrr_v l (m+1) angMom_a cmm totAngMom_a (totAngMom_c-2)
else empty
in
Array.init ncoef (fun k -> p1.(k) +. f1.(k) *. v1.(k) +. f2.(k) *. v2.(k))
Array.init nq (fun k -> p1.(k) +. f1.{k+1} *. v1.(k) +. f2.{k+1} *. v2.(k))
in
if (axyz < 1) || (cxyz < 1) then p2 else
let fa =
(float_of_int axyz) *. expo_inv_p *. 0.5
(float_of_int axyz) *. expo_inv_p.{l} *. 0.5
in
let f1 =
Array.map (fun e -> fa *. e ) expo_inv_q
Vec.map (fun e -> fa *. e ) expo_inv_q
in
if (at_least_one_valid f1) then
let v =
vrr_v (m+1) am cm (totAngMom_a-1) (totAngMom_c-1)
vrr_v l (m+1) am cm (totAngMom_a-1) (totAngMom_c-1)
in
Array.init ncoef (fun k -> p2.(k) -. f1.(k) *. v.(k))
Array.init nq (fun k -> p2.(k) -. f1.{k+1} *. v.(k))
else p2
end
in Zmap.add map_2d.(m) key result;
in Zmap.add map_2d.(m).(l-1) key result;
result
(** Horizontal recurrence relations *)
and hrr0_v angMom_a angMom_b angMom_c
and hrr0_v l angMom_a angMom_b angMom_c
totAngMom_a totAngMom_b totAngMom_c =
match totAngMom_b with
| 0 ->
begin
match (totAngMom_a, totAngMom_c) with
| (0,0) -> Array.map2 ( *. ) zero_m_array.(0) coef_prod
| (_,0) -> vrr0_v 0 angMom_a totAngMom_a
| (_,_) -> vrr_v 0 angMom_a angMom_c totAngMom_a totAngMom_c
| (0,0) ->
Array.init nq (fun k -> zero_m_array.(0).{k+1,l} *. coef_prod.{k+1,l})
|> sum
| (_,0) -> vrr0_v l 0 angMom_a totAngMom_a |> sum
| (_,_) -> vrr_v l 0 angMom_a angMom_c totAngMom_a totAngMom_c |> sum
end
| 1 ->
let (aax, aay, aaz) = angMom_a in
@ -210,13 +222,13 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
in
let f = Coordinate.coord center_ab xyz in
let v1 =
vrr_v 0 ap angMom_c (totAngMom_a+1) totAngMom_c
vrr_v l 0 ap angMom_c (totAngMom_a+1) totAngMom_c
in
if (abs_float f < cutoff) then v1 else
if (abs_float f < cutoff) then sum v1 else
let v2 =
vrr_v 0 angMom_a angMom_c totAngMom_a totAngMom_c
vrr_v l 0 angMom_a angMom_c totAngMom_a totAngMom_c
in
Array.map2 (fun v1 v2 -> v1 +. v2 *. f) v1 v2
Array.map2 (fun v1 v2 -> v1 +. v2 *. f) v1 v2 |> sum
| _ ->
let (aax, aay, aaz) = angMom_a
and (abx, aby, abz) = angMom_b in
@ -226,7 +238,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
| (0,_,_) -> aby, 1
| _ -> abx, 0
in
if (bxyz < 1) then empty else
if (bxyz < 1) then 0. else
let ap, bm =
match xyz with
| 0 -> (aax+1,aay,aaz),(abx-1,aby,abz)
@ -235,22 +247,24 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
in
let h1 =
hrr0_v ap bm angMom_c (totAngMom_a+1) (totAngMom_b-1) totAngMom_c
hrr0_v l ap bm angMom_c (totAngMom_a+1) (totAngMom_b-1) totAngMom_c
in
let f = (Coordinate.coord center_ab xyz) in
if (abs_float f < cutoff) then h1 else
let h2 =
hrr0_v angMom_a bm angMom_c totAngMom_a (totAngMom_b-1) totAngMom_c
in Array.map2 (fun h1 h2 -> h1 +. h2 *. f) h1 h2
hrr0_v l angMom_a bm angMom_c totAngMom_a (totAngMom_b-1) totAngMom_c
in
h1 +. h2 *. f
and hrr_v angMom_a angMom_b angMom_c angMom_d
and hrr_v l angMom_a angMom_b angMom_c angMom_d
totAngMom_a totAngMom_b totAngMom_c totAngMom_d =
match (totAngMom_b, totAngMom_d) with
| (_,0) -> if (totAngMom_b = 0) then
vrr_v 0 angMom_a angMom_c totAngMom_a totAngMom_c
vrr_v l 0 angMom_a angMom_c totAngMom_a totAngMom_c
|> sum
else
hrr0_v angMom_a angMom_b angMom_c totAngMom_a totAngMom_b totAngMom_c
hrr0_v l angMom_a angMom_b angMom_c totAngMom_a totAngMom_b totAngMom_c
| (_,_) ->
let (acx, acy, acz) = angMom_c
and (adx, ady, adz) = angMom_d in
@ -261,19 +275,21 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
| _ -> (acx, acy, acz+1), (adx, ady, adz-1), 2
in
let h1 =
hrr_v angMom_a angMom_b cp dm totAngMom_a totAngMom_b (totAngMom_c+1) (totAngMom_d-1)
hrr_v l angMom_a angMom_b cp dm totAngMom_a totAngMom_b (totAngMom_c+1) (totAngMom_d-1)
and h2 =
hrr_v angMom_a angMom_b angMom_c dm totAngMom_a totAngMom_b totAngMom_c (totAngMom_d-1)
hrr_v l angMom_a angMom_b angMom_c dm totAngMom_a totAngMom_b totAngMom_c (totAngMom_d-1)
in
let f = (Coordinate.coord center_cd xyz) in
Array.init ncoef (fun k -> h1.(k) +. h2.(k) *. f)
h1 +. f *. h2
in
hrr_v
(angMom_a.(0),angMom_a.(1),angMom_a.(2))
(angMom_b.(0),angMom_b.(1),angMom_b.(2))
(angMom_c.(0),angMom_c.(1),angMom_c.(2))
(angMom_d.(0),angMom_d.(1),angMom_d.(2))
totAngMom_a totAngMom_b totAngMom_c totAngMom_d
Array.init np (fun ab ->
hrr_v (ab+1)
(angMom_a.(0),angMom_a.(1),angMom_a.(2))
(angMom_b.(0),angMom_b.(1),angMom_b.(2))
(angMom_c.(0),angMom_c.(1),angMom_c.(2))
(angMom_d.(0),angMom_d.(1),angMom_d.(2))
totAngMom_a totAngMom_b totAngMom_c totAngMom_d
) |> sum
@ -283,10 +299,12 @@ 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 shell_a = shell_p.(0).ShellPair.shell_a
and shell_b = shell_p.(0).ShellPair.shell_b
and shell_c = shell_q.(0).ShellPair.shell_a
and shell_d = shell_q.(0).ShellPair.shell_b
let shell_a = shell_p.ContractedShellPair.shell_a
and shell_b = shell_p.ContractedShellPair.shell_b
and shell_c = shell_q.ContractedShellPair.shell_a
and shell_d = shell_q.ContractedShellPair.shell_b
and sp = shell_p.ContractedShellPair.shell_pairs
and sq = shell_q.ContractedShellPair.shell_pairs
in
let maxm =
let open Angular_momentum in
@ -308,134 +326,149 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
(* Compute all integrals in the shell for each pair of significant shell pairs *)
let expo_inv_p =
Array.map (fun shell_ab -> shell_ab.ShellPair.expo_inv) sp
|> Vec.of_array
and expo_inv_q =
Array.map (fun shell_cd -> shell_cd.ShellPair.expo_inv) sq
|> Vec.of_array
in
let np, nq =
Vec.dim expo_inv_p, Vec.dim expo_inv_q
in
let coef =
let result = Mat.make0 nq np in
Lacaml.D.ger
(Vec.of_array shell_q.ContractedShellPair.coef)
(Vec.of_array shell_p.ContractedShellPair.coef)
result;
result
in
begin
match Contracted_shell.(totAngMom shell_a, totAngMom shell_b,
totAngMom shell_c, totAngMom shell_d) with
| 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.ShellPair.coef *. shell_cd.ShellPair.coef
in
(** Screening on the product of coefficients *)
try
if (abs_float coef_prod) < 1.e-3*.cutoff then
raise NullQuartet;
contracted_class.(0) <-
let zm_array = Mat.init_rows np nq (fun i j ->
(** Screening on the product of coefficients *)
try
if (abs_float coef.{j,i} ) < 1.e-3*.cutoff then
raise NullQuartet;
let expo_pq_inv =
shell_ab.ShellPair.expo_inv +. shell_cd.ShellPair.expo_inv
in
let center_pq =
Coordinate.(shell_ab.ShellPair.center |- shell_cd.ShellPair.center)
in
let norm_pq_sq =
Coordinate.dot center_pq center_pq
in
let zero_m_array =
zero_m ~maxm:0 ~expo_pq_inv ~norm_pq_sq
in
accu +. coef_prod *. zero_m_array.(0)
with NullQuartet -> accu
) 0. shell_q
) 0. shell_p
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
let norm_pq_sq =
Coordinate.dot center_pq center_pq
in
let zero_m_array =
zero_m ~maxm:0 ~expo_pq_inv ~norm_pq_sq
in
zero_m_array.(0)
with NullQuartet -> 0.
) in
Mat.gemm_trace zm_array coef
| _ ->
Array.iter (fun shell_ab ->
let norm_coef_scale_p = shell_ab.ShellPair.norm_coef_scale in
let b = shell_ab.ShellPair.j in
let common =
Array.map (fun shell_cd ->
let coef_prod =
shell_ab.ShellPair.coef *. shell_cd.ShellPair.coef
in
let expo_pq_inv =
shell_ab.ShellPair.expo_inv +. shell_cd.ShellPair.expo_inv
in
let center_pq =
Coordinate.(shell_ab.ShellPair.center |- shell_cd.ShellPair.center)
in
let norm_pq_sq =
Coordinate.dot center_pq center_pq
in
let expo_b =
Array.map (fun shell_ab -> Contracted_shell.expo shell_b shell_ab.ShellPair.j) sp
|> Vec.of_array
and expo_d =
Array.map (fun shell_cd -> Contracted_shell.expo shell_d shell_cd.ShellPair.j) sq
|> Vec.of_array
in
let norm_coef_scale_p = shell_p.ContractedShellPair.norm_coef_scale in
let center_pq =
let result =
Array3.create Float64 fortran_layout 3 nq np
in
Array.iteri (fun ab shell_ab ->
Array.iteri (fun cd shell_cd ->
let cpq =
Coordinate.(shell_ab.ShellPair.center |- shell_cd.ShellPair.center)
in
result.{1,cd+1,ab+1} <- Coordinate.x cpq;
result.{2,cd+1,ab+1} <- Coordinate.y cpq;
result.{3,cd+1,ab+1} <- Coordinate.z cpq;
) sq
) sp;
result
in
let zero_m_array =
let result =
Array.init (maxm+1) (fun _ -> Mat.make0 nq np)
in
Array.iteri (fun ab shell_ab ->
let zero_m_array_tmp =
Array.mapi (fun cd shell_cd ->
let expo_pq_inv =
expo_inv_p.{ab+1} +. expo_inv_q.{cd+1}
in
let norm_pq_sq =
center_pq.{1,cd+1,ab+1} *. center_pq.{1,cd+1,ab+1} +.
center_pq.{2,cd+1,ab+1} *. center_pq.{2,cd+1,ab+1} +.
center_pq.{3,cd+1,ab+1} *. center_pq.{3,cd+1,ab+1}
in
let zero_m_array =
zero_m ~maxm ~expo_pq_inv ~norm_pq_sq
in
let d = shell_cd.ShellPair.j in
(zero_m_array, shell_cd.ShellPair.expo_inv,
Contracted_shell.expo shell_d d,
center_pq,coef_prod)
) shell_q
|> Array.to_list
|> List.filter (fun (zero_m_array, expo_inv, d,
center_pq,coef_prod) -> abs_float coef_prod >= 1.e-4 *. cutoff)
|> Array.of_list
in
let zero_m_array = Array.map (fun (zero_m_array, expo_inv, d,
center_pq,coef_prod) -> zero_m_array) common
and expo_inv = Array.map (fun (zero_m_array, expo_inv, d,
center_pq,coef_prod) -> expo_inv ) common
and d = Array.map (fun (zero_m_array, expo_inv, d,
center_pq,coef_prod) -> d) common
and center_pq = Array.map (fun (zero_m_array, expo_inv, d,
center_pq,coef_prod) -> center_pq) common
and coef_prod = Array.map (fun (zero_m_array, expo_inv, d,
center_pq,coef_prod) -> coef_prod) common
in
(* Transpose zero_m_array
*)
let zero_m_array =
let result = Array.init (maxm+1) (fun _ ->
Array.make (Array.length coef_prod) 0.)
) sq
(*
|> Array.to_list
|> List.filter (fun (zero_m_array, d,
center_pq,coef_prod) -> abs_float coef_prod >= 1.e-4 *. cutoff)
|> Array.of_list
*)
in
for m=0 to maxm do
for k=0 to (Array.length coef_prod-1) do
result.(m).(k) <- zero_m_array.(k).(m)
done;
done;
result
in
(* Transpose result *)
for m=0 to maxm do
for cd=1 to nq do
result.(m).{cd,ab+1} <- zero_m_array_tmp.(cd-1).(m)
done
done
) sp;
result
in
let norm =
let norm_coef_scale_q = shell_q.ContractedShellPair.norm_coef_scale in
Array.map (fun v1 ->
Array.map (fun v2 -> v1 *. v2) norm_coef_scale_q
) norm_coef_scale_p
|> Array.to_list
|> Array.concat
in
let map_1d = Array.init maxm (fun _ -> Array.init np (fun _ -> Zmap.create (4*maxm)))
and map_2d = Array.init maxm (fun _ -> Array.init np (fun _ -> Zmap.create (Array.length class_indices)))
in
(* Compute the integral class from the primitive shell quartet *)
let map_1d = Array.init maxm (fun _ -> Zmap.create (4*maxm)) in
let map_2d = Array.init maxm (fun _ -> Zmap.create (Array.length class_indices)) in
let norm =
let norm_coef_scale_q = shell_q.(0).ShellPair.norm_coef_scale in
Array.map (fun v1 ->
Array.map (fun v2 -> v1 *. v2) norm_coef_scale_q
) norm_coef_scale_p
|> Array.to_list
|> Array.concat
in
Array.iteri (fun i key ->
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
let integral =
hvrr_two_e_vector (angMomA, angMomB, angMomC, angMomD)
(Contracted_shell.totAngMom shell_a, Contracted_shell.totAngMom shell_b,
Contracted_shell.totAngMom shell_c, Contracted_shell.totAngMom shell_d)
(maxm, zero_m_array)
(Contracted_shell.expo shell_b b, d)
(shell_ab.ShellPair.expo_inv, expo_inv)
(shell_p.(0).ShellPair.center_ab,
shell_q.(0).ShellPair.center_ab, center_pq)
coef_prod map_1d map_2d
in
let x = Array.fold_left (+.) 0. integral in
contracted_class.(i) <- contracted_class.(i) +. x *. norm.(i)
) class_indices
) shell_p
Array.iteri (fun i key ->
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
let integral =
hvrr_two_e_vector (angMomA, angMomB, angMomC, angMomD)
(Contracted_shell.totAngMom shell_a, Contracted_shell.totAngMom shell_b,
Contracted_shell.totAngMom shell_c, Contracted_shell.totAngMom shell_d)
(maxm, zero_m_array)
(expo_b, expo_d)
(expo_inv_p, expo_inv_q)
(shell_p.ContractedShellPair.center_ab,
shell_q.ContractedShellPair.center_ab, center_pq)
coef map_1d map_2d
in
contracted_class.(i) <- contracted_class.(i) +. integral *. norm.(i)
) class_indices
end;