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

Fixed Overlap

This commit is contained in:
Anthony Scemama 2018-02-06 16:33:28 +01:00
parent 666a4f4f69
commit cccf8d40b0
4 changed files with 14 additions and 28 deletions

View File

@ -32,20 +32,6 @@ let contracted_class_shell_pair shell_p geometry: float Zmap.t =
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 geometry =
let to_int_tuple x =
let open Zkey in

View File

@ -31,7 +31,7 @@ let contracted_class shell_a shell_b : float Zmap.t =
let center_ab =
shell_p.(ab).Shell_pair.center_ab
in
let center_a =
let center_pa =
shell_p.(ab).Shell_pair.center_a
in
let expo_inv =
@ -51,10 +51,9 @@ let contracted_class shell_a shell_b : float Zmap.t =
Overlap_primitives.hvrr (angMomA.(k), angMomB.(k))
expo_inv
(Coordinate.coord center_ab k,
Coordinate.coord center_a k)
Coordinate.coord center_pa 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

View File

@ -11,8 +11,8 @@ let hvrr (angMom_a, angMom_b) (expo_inv_p) (center_ab, center_pa)
if angMom_a < 0 then 0.
else if angMom_a = 0 then 1.
else
chop center_pa (fun () -> vrr (angMom_a-1))
+. chop (0.5 *. (float_of_int (angMom_a-1)) *. expo_inv_p)
(chop center_pa (fun () -> vrr (angMom_a-1)))
+. chop (0.5 *. (float_of_int (angMom_a-1)) *. expo_inv_p)
(fun () -> vrr (angMom_a-2))

View File

@ -2,15 +2,16 @@ open Util
open Constants
type t = {
expo : float;
expo_inv : float;
center_ab: Coordinate.t;
center_a : Coordinate.t;
center : Coordinate.t;
norm_sq : float;
norm_coef: float;
coef : float;
norm_coef_scale : float array;
expo : float; (* alpha + beta *)
expo_inv : float; (* 1/(alpha + beta) *)
center : Coordinate.t; (* P = (alpha * A + beta B)/(alpha+beta) *)
center_a : Coordinate.t; (* P - A *)
center_ab: Coordinate.t; (* A - B *)
norm_sq : float; (* |A-B|^2 *)
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;