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

Fixed multipoles

This commit is contained in:
Anthony Scemama 2020-04-19 23:24:52 +02:00
parent 28e6241346
commit d7be3ac86e
1 changed files with 21 additions and 37 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

58
Basis/Multipole.ml
View File

@ -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);