LCPQ
/
QCaml
mirror of https://gitlab.com/scemama/QCaml.git
10
1
Fork 0
Code Releases Activity

Working on Nuclear integrals

This commit is contained in:
Anthony Scemama 2018-02-03 23:26:20 +01:00
parent c836e3a74c
commit 925cdf693f
6 changed files with 400 additions and 216 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

179
Basis/ERI.ml
View File

@ -25,6 +25,10 @@ let zero_m ~maxm ~expo_pq_inv ~norm_pq_sq =
(** Compute all the integrals of a contracted class *)
(*
let contracted_class shell_a shell_b shell_c shell_d : float Zmap.t =
TwoElectronRRVectorized.contracted_class ~zero_m shell_a shell_b shell_c shell_d
*)
let contracted_class shell_a shell_b shell_c shell_d : float Zmap.t =
TwoElectronRR.contracted_class ~zero_m shell_a shell_b shell_c shell_d
@ -37,6 +41,9 @@ let contracted_class_shell_pairs ?schwartz_p ?schwartz_q shell_p shell_q : float
let cutoff2 = cutoff *. cutoff
(*
type n_cls = { n : int ; cls : Z.t array }
*)
exception NullIntegral
(*
@ -66,7 +73,10 @@ let to_file ~filename basis =
(* Pre-compute all shell pairs *)
let shell_pairs =
Array.mapi (fun i shell_a -> Array.map (fun shell_b ->
(*
Shell_pair.create_array shell_a shell_b) (Array.sub basis 0 (i+1)) ) basis
*)
Shell_pair.create_array shell_a shell_b) (basis) ) basis
in
Printf.printf "%d shells\n" (Array.length basis);
@ -87,8 +97,8 @@ let to_file ~filename basis =
print_int basis.(i).Contracted_shell.indice ; print_newline ();
for j=0 to i do
let schwartz_p, schwartz_p_max = schwartz.(i).(j) in
if (schwartz_p_max >= cutoff) then
icount := !icount + 1;
if (schwartz_p_max >= cutoff) then
icount := !icount + 1;
done;
done;
Printf.printf "%d shell pairs computed in %f seconds\n" !icount (Unix.gettimeofday () -. t0);
@ -111,83 +121,111 @@ let to_file ~filename basis =
for i=0 to (Array.length basis) - 1 do
print_int basis.(i).Contracted_shell.indice ; print_newline ();
(*
for j=0 to i do
*)
for j=0 to (Array.length basis) - 1 do
(*
let schwartz_p, schwartz_p_max = schwartz.(i).(j) in
*)
try
(*
if (schwartz_p_max < cutoff) then raise NullIntegral;
*)
let
shell_p = shell_pairs.(i).(j)
in
(*
for k=0 to i do
for l=0 to k do
*)
for k=0 to (Array.length basis) - 1 do
for l=0 to (Array.length basis) - 1 do
(*
let schwartz_q, schwartz_q_max = schwartz.(k).(l) in
*)
try
if schwartz_p_max *. schwartz_q_max < cutoff2 then
raise NullIntegral;
let
shell_q = shell_pairs.(k).(l)
in
(*
if schwartz_p_max *. schwartz_q_max < cutoff2 then
raise NullIntegral;
*)
let
shell_q = shell_pairs.(k).(l)
in
let swap =
Array.length shell_q < Array.length shell_p
in
let swap =
Array.length shell_q < Array.length shell_p
in
(* Compute all the integrals of the class *)
let cls =
if swap then
if Array.length shell_p < 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
(* Compute all the integrals of the class *)
let cls =
if swap then
if Array.length shell_p < 2 then
(*
contracted_class_shell_pairs ~schwartz_p:schwartz_q ~schwartz_q:schwartz_p shell_q shell_p
*)
contracted_class_shell_pairs shell_q shell_p
else
(*
contracted_class_shell_pairs_vec ~schwartz_p:schwartz_q ~schwartz_q:schwartz_p shell_q shell_p
*)
contracted_class_shell_pairs_vec shell_q shell_p
else
if Array.length shell_q < 2 then
contracted_class_shell_pairs
~schwartz_p ~schwartz_q
shell_p shell_q
(*
contracted_class_shell_pairs ~schwartz_p ~schwartz_q shell_p shell_q
*)
contracted_class_shell_pairs shell_p shell_q
else
contracted_class_shell_pairs_vec
~schwartz_p ~schwartz_q
shell_p shell_q
in
(*
contracted_class_shell_pairs_vec ~schwartz_p ~schwartz_q shell_p shell_q
*)
contracted_class_shell_pairs_vec shell_p shell_q
in
(* Write the data in the output file *)
Array.iteri (fun i_c powers_i ->
let i_c = basis.(i).Contracted_shell.indice + i_c + 1 in
let xi = to_int_tuple powers_i in
Array.iteri (fun j_c powers_j ->
let j_c = basis.(j).Contracted_shell.indice + j_c + 1 in
let xj = to_int_tuple powers_j in
Array.iteri (fun k_c powers_k ->
let k_c = basis.(k).Contracted_shell.indice + k_c + 1 in
let xk = to_int_tuple powers_k in
Array.iteri (fun l_c powers_l ->
let l_c = basis.(l).Contracted_shell.indice + l_c + 1 in
let xl = to_int_tuple powers_l in
let key =
if swap then
Zkey.of_int_tuple (Zkey.Twelve (xk,xl,xi,xj))
else
Zkey.of_int_tuple (Zkey.Twelve (xi,xj,xk,xl))
in
let value =
Zmap.find cls key
in
if (abs_float value > cutoff) then
(inn := !inn + 1;
eri_array.{(i_c-1),(k_c-1),(j_c-1),(l_c-1)} <- value;
)
(* Write the data in the output file *)
Array.iteri (fun i_c powers_i ->
let i_c = basis.(i).Contracted_shell.indice + i_c + 1 in
let xi = to_int_tuple powers_i in
Array.iteri (fun j_c powers_j ->
let j_c = basis.(j).Contracted_shell.indice + j_c + 1 in
let xj = to_int_tuple powers_j in
Array.iteri (fun k_c powers_k ->
let k_c = basis.(k).Contracted_shell.indice + k_c + 1 in
let xk = to_int_tuple powers_k in
Array.iteri (fun l_c powers_l ->
let l_c = basis.(l).Contracted_shell.indice + l_c + 1 in
let xl = to_int_tuple powers_l in
let key =
if swap then
Zkey.of_int_tuple (Zkey.Twelve (xk,xl,xi,xj))
else
out := !out + 1;
) basis.(l).Contracted_shell.powers
) basis.(k).Contracted_shell.powers
) basis.(j).Contracted_shell.powers
) basis.(i).Contracted_shell.powers;
with NullIntegral -> ()
Zkey.of_int_tuple (Zkey.Twelve (xi,xj,xk,xl))
in
let value =
Zmap.find cls key
in
eri_array.{(i_c-1),(k_c-1),(j_c-1),(l_c-1)} <- value;
(*
eri_array.{(j_c-1),(k_c-1),(i_c-1),(l_c-1)} <- value;
eri_array.{(i_c-1),(l_c-1),(j_c-1),(k_c-1)} <- value;
eri_array.{(j_c-1),(l_c-1),(i_c-1),(k_c-1)} <- value;
eri_array.{(k_c-1),(i_c-1),(l_c-1),(j_c-1)} <- value;
eri_array.{(k_c-1),(j_c-1),(l_c-1),(i_c-1)} <- value;
eri_array.{(l_c-1),(i_c-1),(k_c-1),(j_c-1)} <- value;
eri_array.{(l_c-1),(j_c-1),(k_c-1),(i_c-1)} <- value;
*)
if (abs_float value > cutoff) then
(inn := !inn + 1;
)
else
out := !out + 1;
) basis.(l).Contracted_shell.powers
) basis.(k).Contracted_shell.powers
) basis.(j).Contracted_shell.powers
) basis.(i).Contracted_shell.powers;
with NullIntegral -> ()
done;
done;
with NullIntegral -> ()
@ -198,15 +236,15 @@ let to_file ~filename basis =
(* Print ERIs *)
for i_c=1 to (Genarray.nth_dim eri_array 0) do
for j_c=1 to (Genarray.nth_dim eri_array 2) do
for k_c=1 to (Genarray.nth_dim eri_array 1) do
for l_c=1 to (Genarray.nth_dim eri_array 3) do
let value = eri_array.{(i_c-1),(k_c-1),(j_c-1),(l_c-1)} in
if (value <> 0.) then
Printf.fprintf oc " %5d %5d %5d %5d%20.15f\n" i_c k_c j_c l_c value;
done;
done;
for j_c=1 to (Genarray.nth_dim eri_array 2) do
for k_c=1 to (Genarray.nth_dim eri_array 1) do
for l_c=1 to (Genarray.nth_dim eri_array 3) do
let value = eri_array.{(i_c-1),(k_c-1),(j_c-1),(l_c-1)} in
if (abs_float value > cutoff) then
Printf.fprintf oc " %5d %5d %5d %5d%20.15f\n" i_c k_c j_c l_c value;
done;
done;
done;
done;
Printf.printf "In: %d Out:%d\n" !inn !out ;
close_out oc
@ -398,8 +436,3 @@ let xto_file ~filename basis =
*)

129
Basis/NucInt.ml Normal file
View File

@ -0,0 +1,129 @@
(** Electron-nucleus repulsion integrals *)
open Util
open Constants
open Bigarray
(** (0|0)^m : Fundamental electron-nucleus repulsion integral
$ \int \phi_p(r1) 1/r_{C} dr_1 $
maxm : Maximum total angular momentum
expo_pq_inv : $1./p + 1./q$ where $p$ and $q$ are the exponents of
$\phi_p$ and $\phi_q$
norm_pq_sq : square of the distance between the centers of $\phi_p$
and $\phi_q$
*)
let zero_m ~maxm ~expo_pq_inv ~norm_pq_sq =
let exp_pq = 1. /. expo_pq_inv in
let t = norm_pq_sq *. exp_pq in
boys_function ~maxm t
|> Array.mapi (fun m fm ->
two_over_sq_pi *. (if m mod 2 = 0 then fm else -.fm) *.
(pow exp_pq m) *. (sqrt exp_pq)
)
(** Compute all the integrals of a contracted class *)
let contracted_class_shell_pair shell_p nucl_coord : float Zmap.t =
OneElectronRR.contracted_class_shell_pair ~zero_m shell_p nucl_coord
let cutoff2 = cutoff *. cutoff
exception NullIntegral
(*
(** Unique index for integral <ij|kl> *)
let index i j k l =
let f i k =
let (p,r) =
if i <= k then (i,k) else (k,i)
in p+ (r*r-r)/2
in
let p = f i k and q = f j l in
f p q
*)
(** Write all integrals to a file with the <ij|kl> convention *)
let to_file ~filename basis nucl_coord =
let to_int_tuple x =
let open Zkey in
match to_int_tuple Kind_3 x with
| Three x -> x
| _ -> assert false
in
let oc = open_out filename in
(* Pre-compute all shell pairs *)
let shell_pairs =
Array.mapi (fun i shell_a -> Array.map (fun shell_b ->
Shell_pair.create_array shell_a shell_b) (Array.sub basis 0 (i+1)) ) basis
in
Printf.printf "%d shells\n" (Array.length basis);
let eni_array =
let n = ref 0 in
for i=0 to (Array.length basis) - 1 do
n := !n + (Array.length (basis.(i).Contracted_shell.powers))
done;
let n = !n in
Array2.create Float64 c_layout n n
in
Array2.fill eni_array 0.;
(* Compute Integrals *)
let t0 = Unix.gettimeofday () in
let inn = ref 0 and out = ref 0 in
for i=0 to (Array.length basis) - 1 do
print_int basis.(i).Contracted_shell.indice ; print_newline ();
for j=0 to i do
let
shell_p = shell_pairs.(i).(j)
in
(* Compute all the integrals of the class *)
let cls =
contracted_class_shell_pair shell_p nucl_coord
in
(* Write the data in the output file *)
Array.iteri (fun i_c powers_i ->
let i_c = basis.(i).Contracted_shell.indice + i_c + 1 in
let xi = to_int_tuple powers_i in
Array.iteri (fun j_c powers_j ->
let j_c = basis.(j).Contracted_shell.indice + j_c + 1 in
let xj = to_int_tuple powers_j in
let key =
Zkey.of_int_tuple (Zkey.Six (xi,xj))
in
let value =
Zmap.find cls key
in
if (abs_float value > cutoff) then
(inn := !inn + 1;
eni_array.{(i_c-1),(j_c-1)} <- value;
)
else
out := !out + 1;
) basis.(j).Contracted_shell.powers
) basis.(i).Contracted_shell.powers;
done;
done;
Printf.printf "Computed %d non-zero ENIs in %f seconds\n" !inn (Unix.gettimeofday () -. t0);
(* Print ENIs *)
for i_c=1 to (Array2.dim1 eni_array) do
for j_c=1 to (Array2.dim2 eni_array) do
let value = eni_array.{(i_c-1),(j_c-1)} in
if (value <> 0.) then
Printf.fprintf oc " %5d %5d%20.15f\n" i_c j_c value;
done;
done;
Printf.printf "In: %d Out:%d\n" !inn !out ;
close_out oc

290
Basis/OneElectronRR.ml
View File

@ -1,4 +1,7 @@
open Util
open Constants
exception NullPair
(** In chop f g, evaluate g only if f is non zero, and return f *. (g ()) *)
let chop f g =
@ -9,7 +12,7 @@ let chop f g =
(** Horizontal and Vertical Recurrence Relations (HVRR) *)
let hvrr_one_e
m (angMom_a, angMom_b) (totAngMom_a, totAngMom_b)
(angMom_a, angMom_b) (totAngMom_a, totAngMom_b)
(maxm, zero_m_array) (expo_inv_p) (center_ab, center_pa, center_pc)
map
=
@ -17,88 +20,99 @@ let hvrr_one_e
let totAngMom_a = Angular_momentum.to_int totAngMom_a
and totAngMom_b = Angular_momentum.to_int totAngMom_b
in
let maxm = totAngMom_a+totAngMom_b in
let maxsze = maxm+1 in
let empty = Array.make maxsze 0. in
(** Vertical recurrence relations *)
let rec vrr m angMom_a totAngMom_a =
if angMom_a.(0) < 0 || angMom_a.(1) < 0 || angMom_a.(2) < 0 then 0.
let rec vrr angMom_a totAngMom_a =
let ax,ay,az = angMom_a in
if (ax < 0) || (ay < 0) || (az < 0) then
empty
else
match totAngMom_a with
| 0 -> zero_m_array.(m)
| 0 -> zero_m_array
| _ ->
let key = Zkey.of_int_tuple (Zkey.Three angMom_a) in
let key = [| angMom_a.(0)+1; angMom_a.(1)+1; angMom_a.(2)+1; |]
|> Zkey.(of_int_array ~kind:Kind_3)
in
let (found, result) =
try (true, Zmap.find map.(m) key) with
| Not_found -> (false,
let am = [| angMom_a.(0) ; angMom_a.(1) ; angMom_a.(2) |]
and amm = [| angMom_a.(0) ; angMom_a.(1) ; angMom_a.(2) |]
and xyz =
match angMom_a with
| [|0;0;_|] -> 2
| [|0;_;_|] -> 1
| _ -> 0
in
am.(xyz) <- am.(xyz) - 1;
amm.(xyz) <- amm.(xyz) - 2;
chop (Coordinate.coord center_pa xyz)
(fun () -> vrr m am (totAngMom_a-1))
+. chop (0.5 *. (float_of_int am.(xyz)) *. expo_inv_p)
(fun () -> vrr m amm (totAngMom_a-2))
-. chop ((Coordinate.coord center_pc xyz) *. expo_inv_p)
(fun () -> vrr (m+1) am (totAngMom_a-1))
-. chop (0.5 *. (float_of_int am.(xyz)) *. expo_inv_p *. expo_inv_p)
(fun () -> vrr (m+1) amm (totAngMom_a-2))
in
if not found then
Zmap.add map.(m) key result;
result
try Zmap.find map key with
| Not_found ->
let result =
let am, amm, amxyz, xyz =
match angMom_a with
| (x,0,0) -> (x-1,0,0),(x-2,0,0), x-1, 0
| (x,y,0) -> (x,y-1,0),(x,y-2,0), y-1, 1
| (x,y,z) -> (x,y,z-1),(x,y,z-2), z-1, 2
in
if amxyz < 0 then empty else
let f1 = Coordinate.coord center_pa xyz
and f2 = expo_inv_p *. (Coordinate.coord center_pc xyz)
in
if amxyz < 1 then
let v1 =
vrr am (totAngMom_a-1)
in
Array.init maxsze (fun m ->
if m = maxm then 0. else (f1 *. v1.(m) ) -. f2 *. v1.(m+1) )
else
let v3 =
vrr amm (totAngMom_a-2)
in
let v1 =
vrr am (totAngMom_a-1)
in
let f3 = (float_of_int amxyz) *. expo_inv_p *. 0.5 in
Array.init maxsze (fun m ->
(if m = maxm then 0. else
(f1 *. v1.(m+1) ) -. f2 *. v1.(m) )
+. f3 *. (v3.(m) +. if m = maxm then 0. else
expo_inv_p *. v3.(m+1))
)
in Zmap.add map key result;
result
(** Horizontal recurrence relations *)
and hrr angMom_a angMom_b totAngMom_a totAngMom_b =
if angMom_b.(0) < 0 || angMom_b.(1) < 0 || angMom_b.(2) < 0 then 0.
let bx,by,bz = angMom_b in
if (bx < 0) || (by < 0) || (bz < 0) then 0.
else
match totAngMom_b with
| 0 -> vrr 0 angMom_a
| 0 -> (vrr angMom_a totAngMom_a).(0)
| _ ->
let key = [| angMom_a.(0)+1; angMom_a.(1)+1; angMom_a.(2)+1;
angMom_b.(0)+1; angMom_b.(1)+1; angMom_b.(2)+1 |]
|> Zkey.(of_int_array ~kind:Kind_6)
let angMom_ax, angMom_ay, angMom_az = angMom_a
and angMom_bx, angMom_by, angMom_bz = angMom_b in
let bxyz, xyz =
match angMom_b with
| (_,0,0) -> angMom_bx, 0
| (_,_,0) -> angMom_by, 1
| (_,_,_) -> angMom_bz, 2
in
let (found, result) =
try (true, Zmap.find map.(m) key) with
| Not_found -> (false,
begin
let ap = [| angMom_a.(0) ; angMom_a.(1) ; angMom_a.(2) |]
and bm = [| angMom_b.(0) ; angMom_b.(1) ; angMom_b.(2) |]
and xyz =
match angMom_b with
| [|0;0;_|] -> 2
| [|0;_;_|] -> 1
| _ -> 0
in
ap.(xyz) <- ap.(xyz) + 1;
bm.(xyz) <- bm.(xyz) - 1;
hrr ap bm (totAngMom_a+1) (totAngMom_b-1)
+. chop (Coordinate.coord center_ab xyz) (fun () ->
hrr angMom_a bm totAngMom_a (totAngMom_b-1) )
end)
in
if not found then
Zmap.add map.(m) key result;
result
if (bxyz < 1) then 0. else
let ap, bm =
match xyz with
| 0 -> (angMom_ax+1,angMom_ay,angMom_az),(angMom_bx-1,angMom_by,angMom_bz)
| 1 -> (angMom_ax,angMom_ay+1,angMom_az),(angMom_bx,angMom_by-1,angMom_bz)
| _ -> (angMom_ax,angMom_ay,angMom_az+1),(angMom_bx,angMom_by,angMom_bz-1)
in
let h1 =
hrr ap bm (totAngMom_a+1) (totAngMom_b-1)
in
let f2 =
(Coordinate.coord center_ab xyz)
in
if (abs_float f2 < cutoff) then h1 else
let h2 =
hrr angMom_a bm totAngMom_a (totAngMom_b-1)
in
h1 +. f2 *. h2
in
hrr m angMom_a angMom_b angMom_c angMom_d totAngMom_a totAngMom_b
totAngMom_c totAngMom_d
hrr
(angMom_a.(0),angMom_a.(1),angMom_a.(2))
(angMom_b.(0),angMom_b.(1),angMom_b.(2))
totAngMom_a totAngMom_b
@ -109,10 +123,12 @@ let hvrr_one_e
(** Computes all the one-electron integrals of the contracted shell pair *)
let contracted_class_nuc ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t =
let contracted_class_shell_pair ~zero_m shell_p nucl_coord : float Zmap.t =
let shell_p = Shell_pair.create_array shell_a shell_b
and maxm =
let shell_a = shell_p.(0).Shell_pair.shell_a
and shell_b = shell_p.(0).Shell_pair.shell_b
in
let maxm =
let open Angular_momentum in
(to_int @@ Contracted_shell.totAngMom shell_a) + (to_int @@ Contracted_shell.totAngMom shell_b)
in
@ -127,80 +143,80 @@ let contracted_class_nuc ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t
let contracted_class =
Array.make (Array.length class_indices) 0.;
in
()
(* Compute all integrals in the shell for each pair of significant shell pairs *)
for ab=0 to (Array.length shell_p - 1)
do
let coef_prod =
shell_p.(ab).Shell_pair.coef
in
(** Screening on thr product of coefficients *)
if (abs_float coef_prod) > 1.e-4*.cutoff then
begin
try
begin
let coef_prod = shell_p.(ab).Shell_pair.coef in
let norm_coef_scale_p = shell_p.(ab).Shell_pair.norm_coef_scale in
let expo_pq_inv =
shell_p.(ab).Shell_pair.expo_inv
(** Screening on the product of coefficients *)
if (abs_float coef_prod) < 1.e-4*.cutoff then
raise NullPair;
let expo_pq_inv =
shell_p.(ab).Shell_pair.expo_inv
in
let center_ab =
shell_p.(ab).Shell_pair.center_ab
in
let center_pa =
Coordinate.(center_ab |- shell_a.Contracted_shell.center)
in
for c=0 to Array.length nucl_coord - 1 do
let center_pc =
Coordinate.(shell_p.(ab).Shell_pair.center |- nucl_coord.(c) )
in
let norm_pq_sq =
Coordinate.dot center_pc center_pc
in
let zero_m_array =
zero_m ~maxm ~expo_pq_inv ~norm_pq_sq
in
match Contracted_shell.(totAngMom shell_a, totAngMom shell_b) with
| Angular_momentum.(S,S) ->
let integral =
zero_m_array.(0)
in
let center_ab =
Coordinate.shell_p.(ab).Shell_pair.center_ab
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
match Contracted_shell.(totAngMom shell_a, totAngMom shell_b,
totAngMom shell_c, totAngMom shell_d) with
| Angular_momentum.(S,S,S,S) -> Array.iteri (fun i key ->
let coef_prod =
shell_p.(ab).Shell_pair.coef *. shell_q.(cd).Shell_pair.coef
in
let integral =
zero_m_array.(0)
in
contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral
) class_indices
| _ ->
let d = shell_q.(cd).Shell_pair.j in
let map = Array.init maxm (fun _ -> Zmap.create (Array.length class_indices)) in
(* Compute the integral class from the primitive shell quartet *)
Array.iteri (fun i key ->
let (angMomA,angMomB,angMomC,angMomD) =
let a = Zkey.to_int_array Zkey.Kind_12 key in
( [| 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 norm =
shell_p.(ab).Shell_pair.norm_fun angMomA angMomB *. shell_q.(cd).Shell_pair.norm_fun angMomC angMomD
in
let integral = chop norm (fun () ->
ghvrr 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)
(maxm, zero_m_array)
(Contracted_shell.expo shell_b b, Contracted_shell.expo shell_d d)
(shell_p.(ab).Shell_pair.expo_inv, shell_q.(cd).Shell_pair.expo_inv)
(shell_p.(ab).Shell_pair.center_ab, shell_q.(cd).Shell_pair.center_ab, center_pq)
map )
in
contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral
) class_indices
end
done
done;
let result =
Zmap.create (Array.length contracted_class)
in
Array.iteri (fun i key -> Zmap.add result key contracted_class.(i)) class_indices;
result
*)
contracted_class.(0) <- contracted_class.(0) +. coef_prod *. integral
| _ ->
let map = Zmap.create (2*maxm) in
let norm_coef_scale = norm_coef_scale_p in
(* Compute the integral class from the primitive shell quartet *)
class_indices
|> Array.iteri (fun i key ->
let (angMomA,angMomB) =
let a = Zkey.to_int_array Zkey.Kind_6 key in
( [| a.(0) ; a.(1) ; a.(2) |],
[| a.(3) ; a.(4) ; a.(5) |])
in
let norm = norm_coef_scale.(i) in
let coef_prod = coef_prod *. norm in
let integral =
hvrr_one_e (angMomA, angMomB)
(Contracted_shell.totAngMom shell_a, Contracted_shell.totAngMom shell_b)
(maxm, zero_m_array) shell_p.(ab).Shell_pair.expo_inv
(center_ab, center_pa, center_pc)
map
in
contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral
)
done
end
with NullPair -> ()
done;
let result =
Zmap.create (Array.length contracted_class)
in
Array.iteri (fun i key -> Zmap.add result key contracted_class.(i)) class_indices;
result

9
Basis/Overlap.ml
View File

@ -1,4 +1,5 @@
open Util
open Constants
(** Computes all the overlap integrals of the contracted shell pair *)
let contracted_class shell_a shell_b : float Zmap.t =
@ -36,6 +37,9 @@ let contracted_class shell_a shell_b : float Zmap.t =
let expo_inv =
shell_p.(ab).Shell_pair.expo_inv
in
let norm_coef_scale =
shell_p.(ab).Shell_pair.norm_coef_scale
in
Array.iteri (fun i key ->
let (angMomA,angMomB) =
@ -43,15 +47,14 @@ let contracted_class shell_a shell_b : float Zmap.t =
( [| a.(0) ; a.(1) ; a.(2) |],
[| a.(3) ; a.(4) ; a.(5) |] )
in
let norm =
shell_p.(ab).Shell_pair.norm_fun angMomA angMomB
in
let f k =
Overlap_primitives.hvrr (angMomA.(k), angMomB.(k))
expo_inv
(Coordinate.coord center_ab k,
Coordinate.coord center_a k)
in
let norm = norm_coef_scale.(i) in
let coef_prod = coef_prod *. norm in
let integral = chop norm (fun () -> (f 0)*.(f 1)*.(f 2)) in
contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral
) class_indices

1
Basis/TwoElectronRR.ml
View File

@ -50,7 +50,6 @@ let hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
match totAngMom_a with
| 0 -> zero_m_array
| _ ->
let maxsze = maxm+1 in
let key = Zkey.of_int_tuple (Zkey.Three angMom_a) in
try Zmap.find map_1d key with

8
run_integrals.ml
View File

@ -27,10 +27,14 @@ let run ~out =
let basis = Lazy.force Basis.basis in
print_endline @@ Basis.to_string basis;
ERI.to_file ~filename:(out_file^".eri") basis
(*
(*
Overlap.to_file ~filename:(out_file^".overlap") basis
ERI.to_file ~filename:(out_file^".eri") basis
*)
let nucl_coord =
Array.map (fun (_,c) -> c) nuclei
in
NucInt.to_file ~filename:(out_file^".nuc") basis nucl_coord
let () =