From d7be3ac86e21492de91cfb4bff43f7edd5b3165f Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sun, 19 Apr 2020 23:24:52 +0200 Subject: [PATCH] Fixed multipoles --- Basis/Multipole.ml | 58 +++++++++++++++++----------------------------- 1 file changed, 21 insertions(+), 37 deletions(-) diff --git a/Basis/Multipole.ml b/Basis/Multipole.ml index 3b4f38f..cbc9d44 100644 --- a/Basis/Multipole.ml +++ b/Basis/Multipole.ml @@ -15,31 +15,15 @@ module Csp = ContractedShellPair module Po = Powers module Psp = PrimitiveShellPair -let multiply a b = - let n = Mat.dim1 a in - let c = Mat.create n n in - Mat.cpab c a b; - c - -let x0 t = t.(0) -let y0 t = t.(1) -let z0 t = t.(2) -let x1 t = t.(3) -let y1 t = t.(4) -let z1 t = t.(5) -let x2 t = t.(6) -let y2 t = t.(7) -let z2 t = t.(8) - -let matrix_x t = multiply (x1 t) @@ multiply (y0 t) (z0 t) -let matrix_y t = multiply (x0 t) @@ multiply (y1 t) (z0 t) -let matrix_z t = multiply (x0 t) @@ multiply (y0 t) (z1 t) -let matrix_x2 t = multiply (x2 t) @@ multiply (y0 t) (z0 t) -let matrix_y2 t = multiply (x0 t) @@ multiply (y2 t) (z0 t) -let matrix_z2 t = multiply (x0 t) @@ multiply (y0 t) (z2 t) -let matrix_xy t = multiply (x1 t) @@ multiply (y1 t) (z0 t) -let matrix_yz t = multiply (x0 t) @@ multiply (y1 t) (z1 t) -let matrix_zx t = multiply (x1 t) @@ multiply (y0 t) (z1 t) +let matrix_x t = t.(0) +let matrix_y t = t.(1) +let matrix_z t = t.(2) +let matrix_x2 t = t.(3) +let matrix_y2 t = t.(4) +let matrix_z2 t = t.(5) +let matrix_xy t = t.(6) +let matrix_yz t = t.(7) +let matrix_zx t = t.(8) let cutoff = integrals_cutoff @@ -86,9 +70,9 @@ let contracted_class shell_a shell_b : float Zmap.t array = begin let expo_inv = Psp.exponent_inv psp and center_pa = Psp.center_minus_a psp - and xa = Co.get Co.X @@ Cs.center shell_a - and ya = Co.get Co.Y @@ Cs.center shell_a - and za = Co.get Co.Z @@ Cs.center shell_a + and xa = Co.(get X) @@ Cs.center shell_a + and ya = Co.(get Y) @@ Cs.center shell_a + and za = Co.(get Z) @@ Cs.center shell_a in Array.iteri (fun i key -> @@ -126,15 +110,15 @@ let contracted_class shell_a shell_b : float Zmap.t array = let z2 = h2 +. za *. (2. *. z -. za *. f2) in let c = contracted_class in let d = coef_prod *. norm in - c.(0).(i) <- c.(0).(i) +. d *. f0 ; - c.(1).(i) <- c.(1).(i) +. d *. f1 ; - c.(2).(i) <- c.(2).(i) +. d *. f2 ; - c.(3).(i) <- c.(3).(i) +. d *. x ; - c.(4).(i) <- c.(4).(i) +. d *. y ; - c.(5).(i) <- c.(5).(i) +. d *. z ; - c.(6).(i) <- c.(6).(i) +. d *. x2 ; - c.(7).(i) <- c.(7).(i) +. d *. y2 ; - c.(8).(i) <- c.(8).(i) +. d *. z2 ; + c.(0).(i) <- c.(0).(i) +. d *. x *. f1 *. f2; + c.(1).(i) <- c.(1).(i) +. d *. f0 *. y *. f2; + c.(2).(i) <- c.(2).(i) +. d *. f0 *. f1 *. z; + c.(3).(i) <- c.(3).(i) +. d *. x2 *. f1 *. f2; + c.(4).(i) <- c.(4).(i) +. d *. f0 *. y2 *. f2; + c.(5).(i) <- c.(5).(i) +. d *. f0 *. f1 *. z2; + c.(6).(i) <- c.(6).(i) +. d *. x *. y *. f2; + c.(7).(i) <- c.(7).(i) +. d *. f0 *. y *. z; + c.(8).(i) <- c.(8).(i) +. d *. x *. f1 *. z; ) class_indices end ) (Csp.coefs_and_shell_pairs shell_p);