mirror of
https://gitlab.com/scemama/QCaml.git
synced 2024-07-25 12:17:25 +02:00
Removed Mat and Vec
This commit is contained in:
parent
032f1a0913
commit
9e9fe01f8d
@ -4,7 +4,7 @@ open Util
|
|||||||
open Constants
|
open Constants
|
||||||
open Bigarray
|
open Bigarray
|
||||||
|
|
||||||
type t = (float, float64_elt, fortran_layout) Bigarray.Genarray.t
|
type t = (float, float32_elt, fortran_layout) Bigarray.Genarray.t
|
||||||
|
|
||||||
(** (00|00)^m : Fundamental electron repulsion integral
|
(** (00|00)^m : Fundamental electron repulsion integral
|
||||||
$ \int \int \phi_p(r1) 1/r_{12} \phi_q(r2) dr_1 dr_2 $
|
$ \int \int \phi_p(r1) 1/r_{12} \phi_q(r2) dr_1 dr_2 $
|
||||||
@ -114,7 +114,7 @@ let of_basis basis =
|
|||||||
|
|
||||||
(* 4D data initialization *)
|
(* 4D data initialization *)
|
||||||
let eri_array =
|
let eri_array =
|
||||||
Genarray.create Float64 fortran_layout [| n ; n ; n ; n|]
|
Genarray.create Float32 fortran_layout [| n ; n ; n ; n|]
|
||||||
in
|
in
|
||||||
Genarray.fill eri_array 0.;
|
Genarray.fill eri_array 0.;
|
||||||
|
|
||||||
|
@ -8,6 +8,11 @@ let cutoff2 = cutoff *. cutoff
|
|||||||
exception NullQuartet
|
exception NullQuartet
|
||||||
exception Found
|
exception Found
|
||||||
|
|
||||||
|
let at_least_one_valid arr =
|
||||||
|
try
|
||||||
|
Array.iter (fun x -> if (abs_float x > cutoff) then raise Found) arr ; false
|
||||||
|
with Found -> true
|
||||||
|
|
||||||
|
|
||||||
(** Horizontal and Vertical Recurrence Relations (HVRR) *)
|
(** Horizontal and Vertical Recurrence Relations (HVRR) *)
|
||||||
let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
|
let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
|
||||||
@ -16,11 +21,9 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
|
|||||||
(expo_b, expo_d)
|
(expo_b, expo_d)
|
||||||
(expo_inv_p, expo_inv_q)
|
(expo_inv_p, expo_inv_q)
|
||||||
(center_ab, center_cd, center_pq)
|
(center_ab, center_cd, center_pq)
|
||||||
coef_prod map_1d map_2d
|
coef_prod map_1d map_2d np nq
|
||||||
=
|
=
|
||||||
|
|
||||||
let nq = Mat.dim1 coef_prod in
|
|
||||||
let np = Mat.dim2 coef_prod in
|
|
||||||
|
|
||||||
let empty =
|
let empty =
|
||||||
Array.make nq 0.
|
Array.make nq 0.
|
||||||
@ -35,23 +38,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
|
|||||||
|
|
||||||
(** Vertical recurrence relations *)
|
(** Vertical recurrence relations *)
|
||||||
let rec vrr0_v m angMom_a = function
|
let rec vrr0_v m angMom_a = function
|
||||||
(*
|
(*
|
||||||
| 1 ->
|
|
||||||
let xyz =
|
|
||||||
match angMom_a with
|
|
||||||
| (1,_,_) -> 0
|
|
||||||
| (_,1,_) -> 1
|
|
||||||
| _ -> 2
|
|
||||||
in
|
|
||||||
let a = Mat.mul center_pq.(xyz) zero_m_array.(m+1) in
|
|
||||||
let b = copy expo_b in
|
|
||||||
scal (Coordinate.coord center_ab xyz) b;
|
|
||||||
let c = Mat.map (fun x -> x) zero_m_array.(m) in
|
|
||||||
Mat.scal_cols c b;
|
|
||||||
let d = Mat.sub a c in
|
|
||||||
Mat.scal_cols d expo_inv_p;
|
|
||||||
Some (Mat.mul coef_prod d)
|
|
||||||
*)
|
|
||||||
| 1 ->
|
| 1 ->
|
||||||
let xyz =
|
let xyz =
|
||||||
match angMom_a with
|
match angMom_a with
|
||||||
@ -60,14 +47,18 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
|
|||||||
| _ -> 2
|
| _ -> 2
|
||||||
in
|
in
|
||||||
Some (
|
Some (
|
||||||
Array.init np (fun ab -> let l=ab+1 in
|
Array.init np (fun l ->
|
||||||
Array.init nq (fun k ->
|
Array.init nq (fun k ->
|
||||||
let f = expo_b.{l} *. (Coordinate.coord center_ab xyz) in
|
let f = expo_b.(l) *. (Coordinate.coord center_ab xyz) in
|
||||||
coef_prod.{k+1,l} *. expo_inv_p.{l} *.
|
coef_prod.(l).(k) *. expo_inv_p.(l) *.
|
||||||
(center_pq.(xyz).{k+1,l} *. zero_m_array.(m+1).{k+1,l}
|
(center_pq.(xyz).(l).(k) *. zero_m_array.(m+1).(l).(k)
|
||||||
-. f *. zero_m_array.(m).{k+1,l} ) ))
|
-. f *. zero_m_array.(m).(l).(k) ) ))
|
||||||
)
|
)
|
||||||
| 0 -> Some (Mat.mul zero_m_array.(m) coef_prod |> Mat.transpose_copy |> Mat.to_array)
|
*)
|
||||||
|
| 0 -> Some (
|
||||||
|
Array.init np (fun l ->
|
||||||
|
Array.init nq (fun k ->
|
||||||
|
zero_m_array.(m).(l).(k) *. coef_prod.(l).(k) ) ) )
|
||||||
| totAngMom_a ->
|
| totAngMom_a ->
|
||||||
let key = Zkey.of_int_tuple (Zkey.Three angMom_a)
|
let key = Zkey.of_int_tuple (Zkey.Three angMom_a)
|
||||||
in
|
in
|
||||||
@ -99,42 +90,42 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
|
|||||||
in
|
in
|
||||||
|
|
||||||
Some (
|
Some (
|
||||||
Array.init np (fun ab -> let l = ab+1 in
|
Array.init np (fun l ->
|
||||||
let v1 =
|
let v1 =
|
||||||
let f =
|
let f =
|
||||||
-. expo_b.{l} *. expo_inv_p.{l} *. (Coordinate.coord center_ab xyz)
|
-. expo_b.(l) *. expo_inv_p.(l) *. (Coordinate.coord center_ab xyz)
|
||||||
in
|
in
|
||||||
match v1_top with
|
match v1_top with
|
||||||
| Some v1_top ->
|
| Some v1_top ->
|
||||||
v1_top.(l-1)
|
v1_top.(l)
|
||||||
|> Array.map (fun x -> f *. x)
|
|> Array.map (fun x -> f *. x)
|
||||||
| None -> empty
|
| None -> empty
|
||||||
in
|
in
|
||||||
let p1 =
|
let p1 =
|
||||||
match p1_top with
|
match p1_top with
|
||||||
| Some p1_top ->
|
| Some p1_top ->
|
||||||
p1_top.(l-1)
|
p1_top.(l)
|
||||||
| _ -> assert false
|
| _ -> assert false
|
||||||
in
|
in
|
||||||
let p1 =
|
let p1 =
|
||||||
Array.init nq (fun k ->
|
Array.init nq (fun k ->
|
||||||
v1.(k) +. expo_inv_p.{l} *. center_pq.(xyz).{k+1,l} *. p1.(k))
|
v1.(k) +. expo_inv_p.(l) *. center_pq.(xyz).(l).(k) *. p1.(k))
|
||||||
in
|
in
|
||||||
if amxyz < 1 then p1 else
|
if amxyz < 1 then p1 else
|
||||||
let f = (float_of_int amxyz) *. expo_inv_p.{l} *. 0.5
|
let f = (float_of_int amxyz) *. expo_inv_p.(l) *. 0.5
|
||||||
in
|
in
|
||||||
let v1 =
|
let v1 =
|
||||||
match v1_top2 with
|
match v1_top2 with
|
||||||
| Some v1_top2 -> v1_top2.(l-1)
|
| Some v1_top2 -> v1_top2.(l)
|
||||||
| None -> assert false
|
| None -> assert false
|
||||||
in
|
in
|
||||||
let v2 =
|
let v2 =
|
||||||
match p1_top2 with
|
match p1_top2 with
|
||||||
| Some p1_top2 -> p1_top2.(l-1)
|
| Some p1_top2 -> p1_top2.(l)
|
||||||
| None -> assert false
|
| None -> assert false
|
||||||
in
|
in
|
||||||
Array.init nq (fun k ->
|
Array.init nq (fun k ->
|
||||||
p1.(k) +. f *. (v1.(k) +. v2.(k) *. expo_inv_p.{l} ) )
|
p1.(k) +. f *. (v1.(k) +. v2.(k) *. expo_inv_p.(l) ) )
|
||||||
)
|
)
|
||||||
)
|
)
|
||||||
in Zmap.add map_1d.(m) key result;
|
in Zmap.add map_1d.(m) key result;
|
||||||
@ -145,14 +136,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
|
|||||||
|
|
||||||
match (totAngMom_a, totAngMom_c) with
|
match (totAngMom_a, totAngMom_c) with
|
||||||
| (i,0) ->
|
| (i,0) ->
|
||||||
if (i>0) then
|
vrr0_v m angMom_a totAngMom_a
|
||||||
begin
|
|
||||||
match vrr0_v m angMom_a totAngMom_a with
|
|
||||||
| Some x -> Some (Mat.of_array x |> Mat.transpose_copy)
|
|
||||||
| None -> None
|
|
||||||
end
|
|
||||||
else
|
|
||||||
Some (Mat.mul zero_m_array.(m) coef_prod )
|
|
||||||
| (_,_) ->
|
| (_,_) ->
|
||||||
|
|
||||||
let key = Zkey.of_int_tuple (Zkey.Six (angMom_a, angMom_c))
|
let key = Zkey.of_int_tuple (Zkey.Six (angMom_a, angMom_c))
|
||||||
@ -183,34 +167,37 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
|
|||||||
in
|
in
|
||||||
let f1 =
|
let f1 =
|
||||||
let f = (Coordinate.coord center_cd xyz) in
|
let f = (Coordinate.coord center_cd xyz) in
|
||||||
Vec.init nq (fun k ->
|
Array.init nq (fun k -> expo_d.(k) *. expo_inv_q.(k) *. f)
|
||||||
expo_d.{k} *. expo_inv_q.{k} *. f)
|
|
||||||
in
|
in
|
||||||
let v1 =
|
let v1 =
|
||||||
if (abs_float @@ amax f1 > cutoff) then
|
if (at_least_one_valid f1) then
|
||||||
vrr_v m angMom_a cm totAngMom_a (totAngMom_c-1)
|
vrr_v m angMom_a cm totAngMom_a (totAngMom_c-1)
|
||||||
else None
|
else None
|
||||||
in
|
in
|
||||||
let f2 =
|
let f2 =
|
||||||
Mat.init_cols nq np (fun k l ->
|
Array.init np (fun l ->
|
||||||
expo_inv_q.{k} *. center_pq.(xyz).{k,l} )
|
Array.init nq (fun k ->
|
||||||
|
expo_inv_q.(k) *. center_pq.(xyz).(l).(k) ) )
|
||||||
in
|
in
|
||||||
let v2 =
|
let v2 =
|
||||||
if (Mat.as_vec f2 |> amax |> abs_float) < cutoff then
|
if (at_least_one_valid (Array.to_list f2 |> Array.concat)) then
|
||||||
None
|
|
||||||
else
|
|
||||||
vrr_v (m+1) angMom_a cm totAngMom_a (totAngMom_c-1)
|
vrr_v (m+1) angMom_a cm totAngMom_a (totAngMom_c-1)
|
||||||
|
else
|
||||||
|
None
|
||||||
in
|
in
|
||||||
|
|
||||||
let p1 =
|
let p1 =
|
||||||
match v1, v2 with
|
match v1, v2 with
|
||||||
| Some v1, Some v2 ->
|
| Some v1, Some v2 ->
|
||||||
Some (Mat.init_cols nq np (fun k l ->
|
Some (Array.init np (fun l ->
|
||||||
-. v1.{k,l} *. f1.{k} -. v2.{k,l} *. f2.{k,l}) )
|
Array.init nq (fun k ->
|
||||||
|
-. v1.(l).(k) *. f1.(k) -. v2.(l).(k) *. f2.(l).(k)) ) )
|
||||||
| None, Some v2 ->
|
| None, Some v2 ->
|
||||||
Some (Mat.init_cols nq np (fun k l -> -. v2.{k,l} *. f2.{k,l}) )
|
Some (Array.init np (fun l ->
|
||||||
|
Array.init nq (fun k -> -. v2.(l).(k) *. f2.(l).(k)) ) )
|
||||||
| Some v1, None ->
|
| Some v1, None ->
|
||||||
Some (Mat.init_cols nq np (fun k l -> -. v1.{k,l} *. f1.{k} ) )
|
Some (Array.init np (fun l ->
|
||||||
|
Array.init nq (fun k -> -. v1.(l).(k) *. f1.(k) ) ) )
|
||||||
| None, None -> None
|
| None, None -> None
|
||||||
in
|
in
|
||||||
|
|
||||||
@ -220,35 +207,43 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
|
|||||||
(float_of_int (cxyz-1)) *. 0.5
|
(float_of_int (cxyz-1)) *. 0.5
|
||||||
in
|
in
|
||||||
let f1 =
|
let f1 =
|
||||||
Vec.map (fun e -> fcm *. e) expo_inv_q
|
Array.init nq (fun k -> fcm *. expo_inv_q.(k))
|
||||||
in
|
in
|
||||||
let f2 =
|
let f2 =
|
||||||
Vec.mul f1 expo_inv_q
|
Array.mapi (fun k x -> x *. expo_inv_q.(k)) f1
|
||||||
in
|
in
|
||||||
let v1 =
|
let v1 =
|
||||||
if (abs_float @@ amax f1 > cutoff) then
|
if (at_least_one_valid f1) then
|
||||||
vrr_v m angMom_a cmm totAngMom_a (totAngMom_c-2)
|
vrr_v m angMom_a cmm totAngMom_a (totAngMom_c-2)
|
||||||
else None
|
else None
|
||||||
in
|
in
|
||||||
let v2 =
|
let v2 =
|
||||||
if (abs_float @@ amax f2 > cutoff) then
|
if (at_least_one_valid f2) then
|
||||||
vrr_v (m+1) angMom_a cmm totAngMom_a (totAngMom_c-2)
|
vrr_v (m+1) angMom_a cmm totAngMom_a (totAngMom_c-2)
|
||||||
else None
|
else None
|
||||||
in
|
in
|
||||||
match p1, v1, v2 with
|
match p1, v1, v2 with
|
||||||
| Some p1, Some v1, Some v2 ->
|
| Some p1, Some v1, Some v2 ->
|
||||||
Some (Mat.init_cols nq np (fun k l -> p1.{k,l} +. f1.{k} *. v1.{k,l} +. f2.{k} *. v2.{k,l}) )
|
Some (Array.init np (fun l ->
|
||||||
|
Array.init nq (fun k ->
|
||||||
|
p1.(l).(k) +. f1.(k) *. v1.(l).(k) +. f2.(k) *. v2.(l).(k)) ) )
|
||||||
| Some p1, Some v1, None ->
|
| Some p1, Some v1, None ->
|
||||||
Some (Mat.init_cols nq np (fun k l -> p1.{k,l} +. f1.{k} *. v1.{k,l} ))
|
Some (Array.init np (fun l ->
|
||||||
|
Array.init nq (fun k -> p1.(l).(k) +. f1.(k) *. v1.(l).(k) ) ) )
|
||||||
| Some p1, None, Some v2 ->
|
| Some p1, None, Some v2 ->
|
||||||
Some (Mat.init_cols nq np (fun k l -> p1.{k,l} +. f2.{k} *. v2.{k,l}) )
|
Some (Array.init np (fun l ->
|
||||||
|
Array.init nq (fun k -> p1.(l).(k) +. f2.(k) *. v2.(l).(k)) ) )
|
||||||
| None , Some v1, Some v2 ->
|
| None , Some v1, Some v2 ->
|
||||||
Some (Mat.init_cols nq np (fun k l -> f1.{k} *. v1.{k,l} +. f2.{k} *. v2.{k,l}) )
|
Some (Array.init np (fun l ->
|
||||||
|
Array.init nq (fun k ->
|
||||||
|
f1.(k) *. v1.(l).(k) +. f2.(k) *. v2.(l).(k)) ) )
|
||||||
| Some p1, None, None -> Some p1
|
| Some p1, None, None -> Some p1
|
||||||
| None , Some v1, None ->
|
| None , Some v1, None ->
|
||||||
Some (Mat.init_cols nq np (fun k l -> f1.{k} *. v1.{k,l}))
|
Some (Array.init np (fun l ->
|
||||||
|
Array.init nq (fun k -> f1.(k) *. v1.(l).(k) ) ) )
|
||||||
| None, None, Some v2 ->
|
| None, None, Some v2 ->
|
||||||
Some (Mat.init_cols nq np (fun k l -> f2.{k} *. v2.{k,l}) )
|
Some (Array.init np (fun l ->
|
||||||
|
Array.init nq (fun k -> f2.(k) *. v2.(l).(k)) ) )
|
||||||
| None, None, None -> None
|
| None, None, None -> None
|
||||||
in
|
in
|
||||||
if (axyz < 1) || (cxyz < 1) then p2 else
|
if (axyz < 1) || (cxyz < 1) then p2 else
|
||||||
@ -258,28 +253,20 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
|
|||||||
begin
|
begin
|
||||||
match (p2, v) with
|
match (p2, v) with
|
||||||
| Some p2, Some v -> Some (
|
| Some p2, Some v -> Some (
|
||||||
Array.init np (fun ab -> let l = ab+1 in
|
Array.init np (fun l ->
|
||||||
let fa =
|
let fa =
|
||||||
(float_of_int axyz) *. expo_inv_p.{l} *. 0.5
|
(float_of_int axyz) *. expo_inv_p.(l) *. 0.5
|
||||||
in
|
in
|
||||||
let f1 =
|
Array.init nq (fun k -> p2.(l).(k) -. expo_inv_q.(k) *. fa *. v.(l).(k))
|
||||||
Vec.map (fun e -> fa *. e ) expo_inv_q
|
) )
|
||||||
in
|
|
||||||
Vec.init nq (fun k -> p2.{k,l} -. f1.{k} *. v.{k,l})
|
|
||||||
)
|
|
||||||
|> Mat.of_col_vecs )
|
|
||||||
| Some p2, None -> Some p2
|
| Some p2, None -> Some p2
|
||||||
| None, Some v -> Some (
|
| None, Some v -> Some (
|
||||||
Array.init np (fun ab -> let l = ab+1 in
|
Array.init np (fun l ->
|
||||||
let fa =
|
let fa =
|
||||||
(float_of_int axyz) *. expo_inv_p.{l} *. 0.5
|
(float_of_int axyz) *. expo_inv_p.(l) *. 0.5
|
||||||
in
|
in
|
||||||
let f1 =
|
Array.init nq (fun k -> -. fa *. expo_inv_q.(k) *. v.(l).(k))
|
||||||
Vec.map (fun e -> fa *. e ) expo_inv_q
|
) )
|
||||||
in
|
|
||||||
Vec.init nq (fun k -> -. f1.{k} *. v.{k,l})
|
|
||||||
)
|
|
||||||
|> Mat.of_col_vecs )
|
|
||||||
| None, None -> None
|
| None, None -> None
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
@ -297,16 +284,25 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
|
|||||||
| 0 ->
|
| 0 ->
|
||||||
begin
|
begin
|
||||||
match (totAngMom_a, totAngMom_c) with
|
match (totAngMom_a, totAngMom_c) with
|
||||||
| (0,0) -> Mat.gemm_trace zero_m_array.(0) coef_prod
|
| (0,0) ->
|
||||||
|
begin
|
||||||
|
let result = ref 0. in
|
||||||
|
for l=0 to np-1 do
|
||||||
|
for k=0 to nq-1 do
|
||||||
|
result := !result +. zero_m_array.(0).(l).(k) *. coef_prod.(l).(k)
|
||||||
|
done
|
||||||
|
done;
|
||||||
|
!result
|
||||||
|
end
|
||||||
| (_,0) ->
|
| (_,0) ->
|
||||||
begin
|
begin
|
||||||
match vrr0_v 0 angMom_a totAngMom_a with
|
match vrr0_v 0 angMom_a totAngMom_a with
|
||||||
| Some matrix -> Mat.sum (Mat.of_array matrix)
|
| Some matrix -> Array.fold_left (fun accu c -> accu +. Array.fold_left (+.) 0. c) 0. matrix
|
||||||
| None -> 0.
|
| None -> 0.
|
||||||
end
|
end
|
||||||
| (_,_) ->
|
| (_,_) ->
|
||||||
match vrr_v 0 angMom_a angMom_c totAngMom_a totAngMom_c with
|
match vrr_v 0 angMom_a angMom_c totAngMom_a totAngMom_c with
|
||||||
| Some matrix -> Mat.sum matrix
|
| Some matrix -> Array.fold_left (fun accu c -> accu +. Array.fold_left (+.) 0. c) 0. matrix
|
||||||
| None -> 0.
|
| None -> 0.
|
||||||
end
|
end
|
||||||
| 1 ->
|
| 1 ->
|
||||||
@ -320,13 +316,13 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
|
|||||||
let f = Coordinate.coord center_ab xyz in
|
let f = Coordinate.coord center_ab xyz in
|
||||||
let v1 =
|
let v1 =
|
||||||
match vrr_v 0 ap angMom_c (totAngMom_a+1) totAngMom_c with
|
match vrr_v 0 ap angMom_c (totAngMom_a+1) totAngMom_c with
|
||||||
| Some matrix -> Mat.sum matrix
|
| Some matrix -> Array.fold_left (fun accu c -> accu +. Array.fold_left (+.) 0. c) 0. matrix
|
||||||
| None -> 0.
|
| None -> 0.
|
||||||
in
|
in
|
||||||
if (abs_float f < cutoff) then v1 else
|
if (abs_float f < cutoff) then v1 else
|
||||||
let v2 =
|
let v2 =
|
||||||
match vrr_v 0 angMom_a angMom_c totAngMom_a totAngMom_c with
|
match vrr_v 0 angMom_a angMom_c totAngMom_a totAngMom_c with
|
||||||
| Some matrix -> Mat.sum matrix
|
| Some matrix -> Array.fold_left (fun accu c -> accu +. Array.fold_left (+.) 0. c) 0. matrix
|
||||||
| None -> 0.
|
| None -> 0.
|
||||||
in
|
in
|
||||||
v1 +. v2 *. f
|
v1 +. v2 *. f
|
||||||
@ -364,7 +360,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
|
|||||||
| (_,0) -> if (totAngMom_b = 0) then
|
| (_,0) -> if (totAngMom_b = 0) then
|
||||||
begin
|
begin
|
||||||
match vrr_v 0 angMom_a angMom_c totAngMom_a totAngMom_c with
|
match vrr_v 0 angMom_a angMom_c totAngMom_a totAngMom_c with
|
||||||
| Some matrix -> Mat.sum matrix
|
| Some matrix -> Array.fold_left (fun accu c -> accu +. Array.fold_left (+.) 0. c) 0. matrix
|
||||||
| None -> 0.
|
| None -> 0.
|
||||||
end
|
end
|
||||||
else
|
else
|
||||||
@ -428,30 +424,30 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
|
|||||||
|
|
||||||
(* Compute all integrals in the shell for each pair of significant shell pairs *)
|
(* Compute all integrals in the shell for each pair of significant shell pairs *)
|
||||||
|
|
||||||
let expo_inv_p =
|
|
||||||
Array.map (fun shell_ab -> shell_ab.ShellPair.expo_inv) sp
|
|
||||||
|> Vec.of_array
|
|
||||||
and expo_inv_q =
|
|
||||||
Array.map (fun shell_cd -> shell_cd.ShellPair.expo_inv) sq
|
|
||||||
|> Vec.of_array
|
|
||||||
in
|
|
||||||
let np, nq =
|
|
||||||
Vec.dim expo_inv_p, Vec.dim expo_inv_q
|
|
||||||
in
|
|
||||||
let coef =
|
|
||||||
let result = Mat.make0 nq np in
|
|
||||||
Lacaml.D.ger
|
|
||||||
(Vec.of_array shell_q.ContractedShellPair.coef)
|
|
||||||
(Vec.of_array shell_p.ContractedShellPair.coef)
|
|
||||||
result;
|
|
||||||
result
|
|
||||||
in
|
|
||||||
|
|
||||||
begin
|
begin
|
||||||
match Contracted_shell.(totAngMom shell_a, totAngMom shell_b,
|
match Contracted_shell.(totAngMom shell_a, totAngMom shell_b,
|
||||||
totAngMom shell_c, totAngMom shell_d) with
|
totAngMom shell_c, totAngMom shell_d) with
|
||||||
| Angular_momentum.(S,S,S,S) ->
|
| Angular_momentum.(S,S,S,S) ->
|
||||||
contracted_class.(0) <-
|
contracted_class.(0) <-
|
||||||
|
let expo_inv_p =
|
||||||
|
Array.map (fun shell_ab -> shell_ab.ShellPair.expo_inv) sp
|
||||||
|
|> Vec.of_array
|
||||||
|
and expo_inv_q =
|
||||||
|
Array.map (fun shell_cd -> shell_cd.ShellPair.expo_inv) sq
|
||||||
|
|> Vec.of_array
|
||||||
|
in
|
||||||
|
let np, nq =
|
||||||
|
Vec.dim expo_inv_p, Vec.dim expo_inv_q
|
||||||
|
in
|
||||||
|
|
||||||
|
let coef =
|
||||||
|
let result = Mat.make0 nq np in
|
||||||
|
Lacaml.D.ger
|
||||||
|
(Vec.of_array shell_q.ContractedShellPair.coef)
|
||||||
|
(Vec.of_array shell_p.ContractedShellPair.coef)
|
||||||
|
result;
|
||||||
|
result
|
||||||
|
in
|
||||||
let zm_array = Mat.init_cols np nq (fun i j ->
|
let zm_array = Mat.init_cols np nq (fun i j ->
|
||||||
(** Screening on the product of coefficients *)
|
(** Screening on the product of coefficients *)
|
||||||
try
|
try
|
||||||
@ -476,47 +472,63 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
|
|||||||
) in
|
) in
|
||||||
Mat.gemm_trace zm_array coef
|
Mat.gemm_trace zm_array coef
|
||||||
| _ ->
|
| _ ->
|
||||||
|
let expo_inv_p =
|
||||||
|
Array.map (fun shell_ab -> shell_ab.ShellPair.expo_inv) sp
|
||||||
|
and expo_inv_q =
|
||||||
|
Array.map (fun shell_cd -> shell_cd.ShellPair.expo_inv) sq
|
||||||
|
in
|
||||||
|
let np, nq =
|
||||||
|
Array.length expo_inv_p, Array.length expo_inv_q
|
||||||
|
in
|
||||||
let expo_b =
|
let expo_b =
|
||||||
Array.map (fun shell_ab -> Contracted_shell.expo shell_b shell_ab.ShellPair.j) sp
|
Array.map (fun shell_ab -> Contracted_shell.expo shell_b shell_ab.ShellPair.j) sp
|
||||||
|> Vec.of_array
|
|
||||||
and expo_d =
|
and expo_d =
|
||||||
Array.map (fun shell_cd -> Contracted_shell.expo shell_d shell_cd.ShellPair.j) sq
|
Array.map (fun shell_cd -> Contracted_shell.expo shell_d shell_cd.ShellPair.j) sq
|
||||||
|> Vec.of_array
|
|
||||||
in
|
in
|
||||||
let norm_coef_scale_p = shell_p.ContractedShellPair.norm_coef_scale in
|
let norm_coef_scale_p = shell_p.ContractedShellPair.norm_coef_scale in
|
||||||
|
|
||||||
|
let coef =
|
||||||
|
Array.init np (fun l ->
|
||||||
|
Array.init nq (fun k ->
|
||||||
|
shell_q.ContractedShellPair.coef.(k) *.
|
||||||
|
shell_p.ContractedShellPair.coef.(l)
|
||||||
|
) )
|
||||||
|
in
|
||||||
|
|
||||||
let center_pq =
|
let center_pq =
|
||||||
Array.init 3 (fun xyz ->
|
Array.init 3 (fun xyz ->
|
||||||
Mat.init_cols nq np (fun cd ab ->
|
Array.init np (fun ab ->
|
||||||
let shell_ab = sp.(ab-1)
|
let shell_ab = sp.(ab) in
|
||||||
and shell_cd = sq.(cd-1)
|
Array.init nq (fun cd ->
|
||||||
in
|
let shell_cd = sq.(cd)
|
||||||
let cpq =
|
in
|
||||||
Coordinate.(shell_ab.ShellPair.center |- shell_cd.ShellPair.center)
|
let cpq =
|
||||||
in
|
Coordinate.(shell_ab.ShellPair.center |- shell_cd.ShellPair.center)
|
||||||
match xyz with
|
in
|
||||||
| 0 -> Coordinate.x cpq;
|
match xyz with
|
||||||
| 1 -> Coordinate.y cpq;
|
| 0 -> Coordinate.x cpq;
|
||||||
| 2 -> Coordinate.z cpq;
|
| 1 -> Coordinate.y cpq;
|
||||||
| _ -> assert false
|
| 2 -> Coordinate.z cpq;
|
||||||
)
|
| _ -> assert false
|
||||||
|
)
|
||||||
)
|
)
|
||||||
|
)
|
||||||
in
|
in
|
||||||
let zero_m_array =
|
let zero_m_array =
|
||||||
let result =
|
let result =
|
||||||
Array.init (maxm+1) (fun _ -> Mat.make0 nq np)
|
Array.init (maxm+1) (fun _ ->
|
||||||
|
Array.init np (fun _ -> Array.make nq 0. ) )
|
||||||
in
|
in
|
||||||
Array.iteri (fun ab shell_ab ->
|
Array.iteri (fun ab shell_ab ->
|
||||||
let zero_m_array_tmp =
|
let zero_m_array_tmp =
|
||||||
Array.mapi (fun cd shell_cd ->
|
Array.mapi (fun cd shell_cd ->
|
||||||
let expo_pq_inv =
|
let expo_pq_inv =
|
||||||
expo_inv_p.{ab+1} +. expo_inv_q.{cd+1}
|
expo_inv_p.(ab) +. expo_inv_q.(cd)
|
||||||
in
|
in
|
||||||
let norm_pq_sq =
|
let norm_pq_sq =
|
||||||
center_pq.(0).{cd+1,ab+1} *. center_pq.(0).{cd+1,ab+1} +.
|
center_pq.(0).(ab).(cd) *. center_pq.(0).(ab).(cd) +.
|
||||||
center_pq.(1).{cd+1,ab+1} *. center_pq.(1).{cd+1,ab+1} +.
|
center_pq.(1).(ab).(cd) *. center_pq.(1).(ab).(cd) +.
|
||||||
center_pq.(2).{cd+1,ab+1} *. center_pq.(2).{cd+1,ab+1}
|
center_pq.(2).(ab).(cd) *. center_pq.(2).(ab).(cd)
|
||||||
in
|
in
|
||||||
|
|
||||||
zero_m ~maxm ~expo_pq_inv ~norm_pq_sq
|
zero_m ~maxm ~expo_pq_inv ~norm_pq_sq
|
||||||
@ -530,8 +542,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
|
|||||||
in
|
in
|
||||||
(* Transpose result *)
|
(* Transpose result *)
|
||||||
for m=0 to maxm do
|
for m=0 to maxm do
|
||||||
for cd=1 to nq do
|
for cd=0 to nq-1 do
|
||||||
result.(m).{cd,ab+1} <- zero_m_array_tmp.(cd-1).(m)
|
result.(m).(ab).(cd) <- zero_m_array_tmp.(cd).(m)
|
||||||
done
|
done
|
||||||
done
|
done
|
||||||
) sp;
|
) sp;
|
||||||
@ -568,7 +580,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
|
|||||||
(expo_inv_p, expo_inv_q)
|
(expo_inv_p, expo_inv_q)
|
||||||
(shell_p.ContractedShellPair.center_ab,
|
(shell_p.ContractedShellPair.center_ab,
|
||||||
shell_q.ContractedShellPair.center_ab, center_pq)
|
shell_q.ContractedShellPair.center_ab, center_pq)
|
||||||
coef map_1d map_2d
|
coef map_1d map_2d np nq
|
||||||
in
|
in
|
||||||
contracted_class.(i) <- contracted_class.(i) +. integral *. norm.(i)
|
contracted_class.(i) <- contracted_class.(i) +. integral *. norm.(i)
|
||||||
) class_indices
|
) class_indices
|
||||||
|
Loading…
Reference in New Issue
Block a user