LCPQ/QCaml
10
1
Fork 0
mirror of https://gitlab.com/scemama/QCaml.git synced 2024-11-19 04:22:21 +01:00
Code Releases Activity

Cleaning in ERI

Browse Source
This commit is contained in:
Anthony Scemama 2018-02-21 17:27:31 +01:00
parent 5fd94119f6
commit 0207cc7497
2 changed files with 135 additions and 152 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

170
Basis/TwoElectronRR.ml
View File

@ -32,7 +32,6 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
let maxm = angMom_a.tot + angMom_b.tot + angMom_c.tot + angMom_d.tot in let maxm = angMom_a.tot + angMom_b.tot + angMom_c.tot + angMom_d.tot in
let maxsze = maxm+1 in let maxsze = maxm+1 in
let empty = Array.make (maxm+1) 0. in
let get_xyz angMom = let get_xyz angMom =
match angMom with match angMom with
@ -57,36 +56,33 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
let am = Powers.decr xyz angMom_a in let am = Powers.decr xyz angMom_a in
let amxyz = Powers.get xyz am in let amxyz = Powers.get xyz am in
if amxyz >= 0 then let f1 = expo_inv_p *. (Coordinate.get xyz center_pq)
let f1 = expo_inv_p *. (Coordinate.get xyz center_pq) and f2 = expo_b *. expo_inv_p *. (Coordinate.get xyz center_ab)
and f2 = expo_b *. expo_inv_p *. (Coordinate.get xyz center_ab) in
in let result = Array.create_float maxsze in
let result = Array.create_float maxsze in if amxyz = 0 then
if amxyz < 1 then begin
begin let v1 =
let v1 = vrr0 am
vrr0 am in
in for m=0 to maxm-1 do
for m=0 to maxm-1 do result.(m) <- f1 *. v1.(m+1) -. f2 *. v1.(m)
result.(m) <- f1 *. v1.(m+1) -. f2 *. v1.(m) done;
done; result.(maxm) <- -. f2 *. v1.(maxm)
result.(maxm) <- -. f2 *. v1.(maxm) end
end
else
begin
let amm = Powers.decr xyz am in
let v3 = vrr0 amm in
let v1 = vrr0 am in
let f3 = (float_of_int amxyz) *. expo_inv_p *. 0.5 in
for m=0 to maxm-1 do
result.(m) <- f1 *. v1.(m+1) -. f2 *. v1.(m)
+. f3 *. (v3.(m) +. expo_inv_p *. v3.(m+1))
done;
result.(maxm) <- f3 *. v3.(maxm)
end;
result
else else
empty begin
let amm = Powers.decr xyz am in
let v3 = vrr0 amm in
let v1 = vrr0 am in
let f3 = (float_of_int amxyz) *. expo_inv_p *. 0.5 in
for m=0 to maxm-1 do
result.(m) <- f1 *. v1.(m+1) -. f2 *. v1.(m)
+. f3 *. (v3.(m) +. expo_inv_p *. v3.(m+1))
done;
result.(maxm) <- f3 *. v3.(maxm)
end;
result
in Zmap.add map_1d key result; in Zmap.add map_1d key result;
result result
@ -102,68 +98,84 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
try Zmap.find map_2d key with try Zmap.find map_2d key with
| Not_found -> | Not_found ->
let result = let result =
(* Invariant : angMom_c.tot > 0 so cm.tot >= 0 *)
let xyz = get_xyz angMom_c in let xyz = get_xyz angMom_c in
let cm = Powers.decr xyz angMom_c in let cm = Powers.decr xyz angMom_c in
let cmxyz = Powers.get xyz cm in let cmxyz = Powers.get xyz cm in
let axyz = Powers.get xyz angMom_a in let axyz = Powers.get xyz angMom_a in
if cmxyz >= 0 then let f1 =
let f1 = -. expo_d *. expo_inv_q *. (Coordinate.get xyz center_cd)
-. expo_d *. expo_inv_q *. (Coordinate.get xyz center_cd) and f2 =
and f2 = expo_inv_q *. (Coordinate.get xyz center_pq)
expo_inv_q *. (Coordinate.get xyz center_pq) in
in let result = Array.make maxsze 0. in
let result = Array.make maxsze 0. in if ( (abs_float f1 > cutoff) || (abs_float f2 > cutoff) ) then
if ( (abs_float f1 > cutoff) || (abs_float f2 > cutoff) ) then begin
let v1 =
vrr angMom_a cm
in
for m=0 to maxm-1 do
result.(m) <- f1 *. v1.(m) -. f2 *. v1.(m+1) ;
done;
result.(maxm) <- f1 *. v1.(maxm) ;
end;
if cmxyz > 0 then
begin
let f3 =
(float_of_int cmxyz) *. expo_inv_q *. 0.5
in
if (abs_float f3 > cutoff) ||
(abs_float (f3 *. expo_inv_q) > cutoff) then
begin begin
let v1 = let v3 =
vrr angMom_a cm let cmm = Powers.decr xyz cm in
vrr angMom_a cmm
in in
for m=0 to maxm-1 do for m=0 to maxm-1 do
result.(m) <- f1 *. v1.(m) -. f2 *. v1.(m+1) ; result.(m) <- result.(m) +.
f3 *. (v3.(m) +. expo_inv_q *. v3.(m+1))
done; done;
result.(maxm) <- f1 *. v1.(maxm) ; result.(maxm) <- result.(maxm) +. f3 *. v3.(maxm)
end; end
if cmxyz > 0 then end;
begin if axyz > 0 then
let f3 = begin
(float_of_int cmxyz) *. expo_inv_q *. 0.5 let am = Powers.decr xyz angMom_a in
in let f5 =
if (abs_float f3 > cutoff) || (float_of_int axyz) *. expo_inv_p *. expo_inv_q *. 0.5
(abs_float (f3 *. expo_inv_q) > cutoff) then in
begin if (abs_float f5 > cutoff) then
let v3 = let v5 =
let cmm = Powers.decr xyz cm in vrr am cm
vrr angMom_a cmm in
in for m=0 to maxm-1 do
for m=0 to maxm-1 do result.(m) <- result.(m) -. f5 *. v5.(m+1)
result.(m) <- result.(m) +. done
f3 *. (v3.(m) +. expo_inv_q *. v3.(m+1)) end;
done; result
result.(maxm) <- result.(maxm) +. f3 *. v3.(maxm)
end
end;
if axyz > 0 && cmxyz >= 0 then
begin
let am = Powers.decr xyz angMom_a in
let f5 =
(float_of_int axyz) *. expo_inv_p *. expo_inv_q *. 0.5
in
if (abs_float f5 > cutoff) then
let v5 =
vrr am cm
in
for m=0 to maxm-1 do
result.(m) <- result.(m) -. f5 *. v5.(m+1)
done
end;
result
else
empty
in Zmap.add map_2d key result; in Zmap.add map_2d key result;
result result
(*
and trr angMom_a angMom_c =
match (angMom_a.tot, angMom_c.tot) with
| (i,0) -> if (i>0) then vrr0 angMom_a
else zero_m_array
| (_,_) ->
let key = Zkey.of_powers (Zkey.Six (angMom_a, angMom_c)) in
try Zmap.find map_2d key with
| Not_found ->
let result =
let xyz = get_xyz angMom_c in
let cm = Powers.decr xyz angMom_c in
let cmxyz = Powers.get xyz cm in
let axyz = Powers.get xyz angMom_a in
*)

117
Basis/TwoElectronRRVectorized.ml
View File

@ -35,7 +35,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 = let rec vrr0_v m angMom_a =
match angMom_a.tot with match angMom_a.tot with
| 0 -> Some zero_m_array.(m) | 0 -> zero_m_array.(m)
| _ -> | _ ->
let key = Zkey.of_powers (Zkey.Three angMom_a) let key = Zkey.of_powers (Zkey.Three angMom_a)
in in
@ -45,83 +45,54 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
let result = let result =
let xyz = get_xyz angMom_a in let xyz = get_xyz angMom_a in
let am = Powers.decr xyz angMom_a in let am = Powers.decr xyz angMom_a in
let amxyz = Powers.get xyz am in let cab = Coordinate.get xyz center_ab in
if amxyz >= 0 then let p0 = vrr0_v (m+1) am in
begin
let cab = Coordinate.get xyz center_ab in
let v1_top, p1_top =
if abs_float cab < cutoff then
None,
vrr0_v (m+1) am
else
vrr0_v m am, vrr0_v (m+1) am
in
let v1_top2, p1_top2 =
if amxyz < 1 then None, None else
let amm = Powers.decr xyz am in
vrr0_v m amm, vrr0_v (m+1) amm
in
let result = Array.make_matrix np nq 0. in let result = Array.make_matrix np nq 0. in
let p0 = begin
match p1_top with if abs_float cab >= cutoff then
| Some p1_top -> p1_top let v0 = vrr0_v m am in
| _ -> assert false
in
begin
match v1_top with
| None -> ()
| Some v0 ->
Array.iteri (fun l result_l ->
let f0 = -. expo_b.(l) *. expo_inv_p.(l) *. cab
and v0_l = v0.(l)
in
Array.iteri (fun k v0_lk ->
result_l.(k) <- v0_lk *. f0) v0_l ) result
end;
let amxyz = Powers.get xyz am in
if amxyz < 1 then
Array.iteri (fun l result_l -> Array.iteri (fun l result_l ->
let expo_inv_p_l = expo_inv_p.(l) let f0 = -. expo_b.(l) *. expo_inv_p.(l) *. cab
and v0_l = v0.(l)
in
Array.iteri (fun k v0_lk ->
result_l.(k) <- v0_lk *. f0) v0_l ) result
end;
let amxyz = Powers.get xyz am in
if amxyz < 1 then
Array.iteri (fun l result_l ->
let expo_inv_p_l = expo_inv_p.(l)
and center_pq_xyz_l = (center_pq xyz).(l)
and result_l = result.(l)
and p0_l = p0.(l)
in
Array.iteri (fun k p0_lk ->
result_l.(k) <- result_l.(k)
+. expo_inv_p_l *. center_pq_xyz_l.(k) *. p0_lk
) p0_l ) result
else
begin
let amm = Powers.decr xyz am in
let v1 = vrr0_v m amm in
let v2 = vrr0_v (m+1) amm in
let amxyz = float_of_int amxyz in
Array.iteri (fun l result_l ->
let f = amxyz *. expo_inv_p.(l) *. 0.5
and expo_inv_p_l = expo_inv_p.(l)
and center_pq_xyz_l = (center_pq xyz).(l) and center_pq_xyz_l = (center_pq xyz).(l)
and v1_l = v1.(l)
and v2_l = v2.(l)
and result_l = result.(l) and result_l = result.(l)
and p0_l = p0.(l)
in in
Array.iteri (fun k p0_lk -> Array.iteri (fun k p0_lk ->
result_l.(k) <- result_l.(k) result_l.(k) <- result_l.(k) +.
+. expo_inv_p_l *. center_pq_xyz_l.(k) *. p0_lk expo_inv_p_l *. center_pq_xyz_l.(k) *. p0_lk +.
) p0_l ) result f *. (v1_l.(k) +. v2_l.(k) *. expo_inv_p_l)
else ) p0.(l)
begin ) result
let v1 = end;
match v1_top2 with result
| Some v1_top2 -> v1_top2
| None -> assert false
in
let v2 =
match p1_top2 with
| Some p1_top2 -> p1_top2
| None -> assert false
in
Array.iteri (fun l result_l ->
let f = float_of_int amxyz *. expo_inv_p.(l) *. 0.5
and expo_inv_p_l = expo_inv_p.(l)
and center_pq_xyz_l = (center_pq xyz).(l)
and v1_l = v1.(l)
and v2_l = v2.(l)
and result_l = result.(l)
in
Array.iteri (fun k p0_lk ->
result_l.(k) <- result_l.(k) +.
expo_inv_p_l *. center_pq_xyz_l.(k) *. p0_lk +.
f *. (v1_l.(k) +. v2_l.(k) *. expo_inv_p_l)
) p0.(l)
) result
end;
Some result
end
else
None
in in
Zmap.add map_1d.(m) key result; Zmap.add map_1d.(m) key result;
result result
@ -129,7 +100,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
and vrr_v m angMom_a angMom_c = and vrr_v m angMom_a angMom_c =
match (angMom_a.tot, angMom_c.tot) with match (angMom_a.tot, angMom_c.tot) with
| (i,0) -> vrr0_v m angMom_a | (i,0) -> Some (vrr0_v m angMom_a)
| (_,_) -> | (_,_) ->
let key = Zkey.of_powers (Zkey.Six (angMom_a, angMom_c)) let key = Zkey.of_powers (Zkey.Six (angMom_a, angMom_c))