From cccf8d40b0a48c02d92724515c2c39e01510f79a Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Tue, 6 Feb 2018 16:33:28 +0100
Subject: [PATCH] Fixed Overlap

---
 Basis/NucInt.ml             | 14 --------------
 Basis/Overlap.ml            |  5 ++---
 Basis/Overlap_primitives.ml |  4 ++--
 Basis/Shell_pair.ml         | 19 ++++++++++---------
 4 files changed, 14 insertions(+), 28 deletions(-)

diff --git a/Basis/NucInt.ml b/Basis/NucInt.ml
index 110a419..12c1a59 100644
--- a/Basis/NucInt.ml
+++ b/Basis/NucInt.ml
@@ -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
diff --git a/Basis/Overlap.ml b/Basis/Overlap.ml
index 276aa3c..e3d7da7 100644
--- a/Basis/Overlap.ml
+++ b/Basis/Overlap.ml
@@ -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
diff --git a/Basis/Overlap_primitives.ml b/Basis/Overlap_primitives.ml
index 1fe7a24..ea6dbfc 100644
--- a/Basis/Overlap_primitives.ml
+++ b/Basis/Overlap_primitives.ml
@@ -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))
 
 
diff --git a/Basis/Shell_pair.ml b/Basis/Shell_pair.ml
index 2f2ae58..a3ae733 100644
--- a/Basis/Shell_pair.ml
+++ b/Basis/Shell_pair.ml
@@ -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;