LCPQ/QCaml
10
1
Fork 0
mirror of https://gitlab.com/scemama/QCaml.git synced 2024-12-22 04:13:33 +01:00
Code Releases Activity

Nuclear Integrals OK

Browse Source
This commit is contained in:
Anthony Scemama 2018-02-06 12:02:00 +01:00
parent db46395e4e
commit 2877a9d492
2 changed files with 17 additions and 10 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
Basis/NucInt.ml
View File

@ -121,7 +121,7 @@ let to_file ~filename basis geometry =
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;
Printf.fprintf oc " %5d %5d %20.15f\n" i_c j_c value;
done;
done;
Printf.printf "In: %d Out:%d\n" !inn !out ;

23
Basis/OneElectronRR.ml
View File

@ -13,7 +13,7 @@ let chop f g =
(** Horizontal and Vertical Recurrence Relations (HVRR) *)
let hvrr_one_e
(angMom_a, angMom_b) (totAngMom_a, totAngMom_b)
(maxm, zero_m_array) (expo_inv_p) (center_ab, center_pa, center_pc)
(maxm, zero_m_array) (expo_b) (expo_inv_p) (center_ab, center_pa, center_pc)
map
=
@ -53,7 +53,8 @@ let hvrr_one_e
vrr am (totAngMom_a-1)
in
Array.init maxsze (fun m ->
if m = maxm then 0. else (f1 *. v1.(m) ) -. f2 *. v1.(m+1) )
if m = maxm then (f1 *. v1.(m) ) else
(f1 *. v1.(m) ) -. f2 *. v1.(m+1) )
else
let v3 =
vrr amm (totAngMom_a-2)
@ -62,10 +63,10 @@ let hvrr_one_e
vrr am (totAngMom_a-1)
in
let f3 = (float_of_int amxyz) *. expo_inv_p *. 0.5 in
Array.init maxsze (fun m ->
Array.init maxsze (fun m -> f1 *. v1.(m) -.
(if m = maxm then 0. else
(f1 *. v1.(m+1) ) -. f2 *. v1.(m) )
+. f3 *. (v3.(m) +. if m = maxm then 0. else
f2 *. v1.(m+1) )
+. f3 *. (v3.(m) -. if m = maxm then 0. else
expo_inv_p *. v3.(m+1))
)
in Zmap.add map key result;
@ -148,6 +149,7 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t =
for ab=0 to (Array.length shell_p - 1)
do
let b = shell_p.(ab).Shell_pair.j in
try
begin
let coef_prod = shell_p.(ab).Shell_pair.coef in
@ -165,15 +167,18 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t =
let center_ab =
shell_p.(ab).Shell_pair.center_ab
in
let center_p =
shell_p.(ab).Shell_pair.center
in
let center_pa =
Coordinate.(center_ab |- shell_a.Contracted_shell.center)
Coordinate.(center_p |- shell_a.Contracted_shell.center)
in
for c=0 to Array.length geometry - 1 do
let element, nucl_coord = geometry.(c) in
let charge = Element.to_charge element |> Charge.to_float in
let center_pc =
Coordinate.(shell_p.(ab).Shell_pair.center |- nucl_coord )
Coordinate.(center_p |- nucl_coord )
in
let norm_pq_sq =
Coordinate.dot center_pc center_pc
@ -204,7 +209,9 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t =
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
(maxm, zero_m_array)
(Contracted_shell.expo shell_b b)
(shell_p.(ab).Shell_pair.expo_inv)
(center_ab, center_pa, center_pc)
map
in