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

Array of m

Browse Source
This commit is contained in:
Anthony Scemama 2018-01-31 14:37:51 +01:00
parent 803fa3a0d9
commit 148052025c
4 changed files with 144 additions and 80 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

44
Basis/Contracted_shell.ml
View File

@ -39,34 +39,42 @@ let to_string s =
exponent and the angular momentum. Two conventions can be chosen : a single
normalisation factor for all functions of the class, or a coefficient which
depends on the powers of x,y and z.
Returns, for each contracted function, an array of functions taking as
argument the [|x;y;z|] powers.
*)
let compute_norm_coef s =
let compute_norm_coef expo totAngMom =
let atot =
Angular_momentum.to_int s.totAngMom
Angular_momentum.to_int totAngMom
in
Array.mapi (fun i alpha ->
let alpha_2 =
alpha +. alpha
in
let c =
(alpha_2 *. pi_inv)**(1.5) *. (pow (alpha_2 +. alpha_2) atot)
in
let result a =
let factor int_array =
let dfa = Array.map (fun j ->
( float_of_int (1 lsl j) *. fact j) /. fact (j+j)
) a
in sqrt (c *. dfa.(0) *.dfa.(1) *. dfa.(2))
) int_array
in
result
) s.expo
sqrt (dfa.(0) *.dfa.(1) *. dfa.(2))
in
let expo =
if atot mod 2 = 0 then
Array.map (fun alpha ->
let alpha_2 = alpha +. alpha in
(alpha_2 *. pi_inv)**(0.75) *. (pow (alpha_2 +. alpha_2) (atot/2))
) expo
else
Array.map (fun alpha ->
let alpha_2 = alpha +. alpha in
(alpha_2 *. pi_inv)**(0.75) *. sqrt (pow (alpha_2 +. alpha_2) atot)
) expo
in
Array.map (fun x -> let f a = x *. (factor a) in f) expo
let create ~indice ~expo ~coef ~center ~totAngMom =
assert (Array.length expo = Array.length coef);
assert (Array.length expo > 0);
let tmp =
{ indice ; expo ; coef ; center ; totAngMom ; size=Array.length expo ; norm_coef = [||];
powers = Angular_momentum.zkey_array (Angular_momentum.Singlet totAngMom) }
let norm_coef =
compute_norm_coef expo totAngMom
in
{ tmp with norm_coef = compute_norm_coef tmp }
{ indice ; expo ; coef ; center ; totAngMom ; size=Array.length expo ; norm_coef ;
powers = Angular_momentum.zkey_array (Angular_momentum.Singlet totAngMom) }

1
Basis/ERI.ml
View File

@ -32,6 +32,7 @@ let contracted_class shell_a shell_b shell_c shell_d : float Zmap.t =
(** Compute all the integrals of a contracted class *)
let contracted_class_shell_pairs_vec ?schwartz_p ?schwartz_q shell_p shell_q : float Zmap.t =
TwoElectronRRVectorized.contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
let contracted_class_shell_pairs ?schwartz_p ?schwartz_q shell_p shell_q : float Zmap.t =
TwoElectronRR.contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q

171
Basis/TwoElectronRR.ml
View File

@ -5,7 +5,7 @@ let cutoff2 = cutoff *. cutoff
exception NullQuartet
(** Horizontal and Vertical Recurrence Relations (HVRR) *)
let hvrr_two_e m (angMom_a, angMom_b, angMom_c, angMom_d)
let hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
(totAngMom_a, totAngMom_b, totAngMom_c, totAngMom_d)
(maxm, zero_m_array)
(expo_b, expo_d)
@ -19,21 +19,20 @@ let hvrr_two_e m (angMom_a, angMom_b, angMom_c, angMom_d)
and totAngMom_c = Angular_momentum.to_int totAngMom_c
and totAngMom_d = Angular_momentum.to_int totAngMom_d
in
let maxm = totAngMom_a + totAngMom_b + totAngMom_c + totAngMom_d in
let empty = Array.make (maxm+1) 0.
in
(** Vertical recurrence relations *)
let rec vrr0 m angMom_a = function
| 1 -> let i = if angMom_a.(0) = 1 then 0 else if angMom_a.(1) = 1 then 1 else 2
in expo_inv_p *.( (Coordinate.coord center_pq i) *. zero_m_array.(m+1)
-. expo_b *. (Coordinate.coord center_ab i) *. zero_m_array.(m) )
| 0 -> zero_m_array.(m)
let rec vrr0 angMom_a = function
| 0 -> zero_m_array
| totAngMom_a ->
let key = Zkey.of_int_tuple (Zkey.Three
(angMom_a.(0)+1, angMom_a.(1)+1, angMom_a.(2)+1) )
in
let (found, result) =
try (true, Zmap.find map.(m) key) with
try (true, Zmap.find map key) with
| Not_found -> (false,
let am = [| angMom_a.(0) ; angMom_a.(1) ; angMom_a.(2) |]
and amm = [| angMom_a.(0) ; angMom_a.(1) ; angMom_a.(2) |]
@ -45,27 +44,37 @@ let hvrr_two_e m (angMom_a, angMom_b, angMom_c, angMom_d)
in
am.(xyz) <- am.(xyz) - 1;
amm.(xyz) <- amm.(xyz) - 2;
if am.(xyz) < 0 then 0. else
chop (-. expo_b *. expo_inv_p *. (Coordinate.coord center_ab xyz))
(fun () -> vrr0 m am (totAngMom_a-1) )
+. chop (expo_inv_p *. (Coordinate.coord center_pq xyz))
(fun () -> vrr0 (m+1) am (totAngMom_a-1) )
+. (if amm.(xyz) < 0 then 0. else
chop ((float_of_int am.(xyz)) *. expo_inv_p *. 0.5)
(fun () -> vrr0 m amm (totAngMom_a-2)
+. chop expo_inv_p (fun () ->
vrr0 (m+1) amm (totAngMom_a-2) ) ) )
if am.(xyz) < 0 then empty else
let v1 =
vrr0 am (totAngMom_a-1)
in
let f1 = expo_inv_p *. (Coordinate.coord center_pq xyz)
and f2 = expo_b *. expo_inv_p *. (Coordinate.coord center_ab xyz)
in
if amm.(xyz) < 0 then
Array.init (maxm+1) (fun m ->
if m = maxm then 0. else (f1 *. v1.(m+1) ) -. f2 *. v1.(m) )
else
let f3 = (float_of_int am.(xyz)) *. expo_inv_p *. 0.5 in
let v3 =
vrr0 amm (totAngMom_a-2)
in
Array.init (maxm+1) (fun m ->
(if m = maxm then 0. else (f1 *. v1.(m+1) ) -. f2 *. v1.(m) )
+. f3 *. (v3.(m) +. if m = maxm then 0. else
expo_inv_p *. v3.(m+1))
)
)
in
if not found then
Zmap.add map.(m) key result;
Zmap.add map key result;
result
and vrr m angMom_a angMom_c totAngMom_a totAngMom_c =
and vrr angMom_a angMom_c totAngMom_a totAngMom_c =
match (totAngMom_a, totAngMom_c) with
| (0,0) -> zero_m_array.(m)
| (_,0) -> vrr0 m angMom_a totAngMom_a
| (0,0) -> zero_m_array
| (_,0) -> vrr0 angMom_a totAngMom_a
| (_,_) ->
let key = Zkey.of_int_tuple (Zkey.Six
@ -75,7 +84,7 @@ let hvrr_two_e m (angMom_a, angMom_b, angMom_c, angMom_d)
in
let (found, result) =
try (true, Zmap.find map.(m) key) with
try (true, Zmap.find map key) with
| Not_found -> (false,
let am = [| angMom_a.(0) ; angMom_a.(1) ; angMom_a.(2) |]
and cm = [| angMom_c.(0) ; angMom_c.(1) ; angMom_c.(2) |]
@ -89,39 +98,75 @@ let hvrr_two_e m (angMom_a, angMom_b, angMom_c, angMom_d)
am.(xyz) <- am.(xyz) - 1;
cm.(xyz) <- cm.(xyz) - 1;
cmm.(xyz) <- cmm.(xyz) - 2;
if cm.(xyz) < 0 then 0. else
chop (-. expo_d *. expo_inv_q *. (Coordinate.coord center_cd xyz) )
(fun () -> vrr m angMom_a cm totAngMom_a (totAngMom_c-1) )
-. chop (expo_inv_q *. (Coordinate.coord center_pq xyz))
(fun () -> vrr (m+1) angMom_a cm totAngMom_a (totAngMom_c-1) )
+. (if cmm.(xyz) < 0 then 0. else
chop ((float_of_int cm.(xyz)) *. expo_inv_q *. 0.5 )
(fun () -> vrr m angMom_a cmm totAngMom_a (totAngMom_c-2)
+. chop expo_inv_q
(fun () -> vrr (m+1) angMom_a cmm totAngMom_a (totAngMom_c-2) ) ) )
-. (if am.(xyz) lor cm.(xyz) < 0 then 0. else
chop ((float_of_int angMom_a.(xyz)) *. expo_inv_p *. expo_inv_q *. 0.5 )
(fun () -> vrr (m+1) am cm (totAngMom_a-1) (totAngMom_c-1) ) ))
if cm.(xyz) < 0 then empty else
let f1 =
-. expo_d *. expo_inv_q *. (Coordinate.coord center_cd xyz)
in
let f2 =
expo_inv_q *. (Coordinate.coord center_pq xyz)
in
let result =
if ( (abs_float f1 < cutoff) && (abs_float f2 < cutoff) ) then empty else
let v1 =
vrr angMom_a cm totAngMom_a (totAngMom_c-1)
in
Array.init (maxm+1) (fun m ->
f1 *. v1.(m) -. (if m = maxm then 0. else f2 *. v1.(m+1)) )
in
let result =
if cmm.(xyz) < 0 then result else
let f3 =
(float_of_int cm.(xyz)) *. expo_inv_q *. 0.5
in
if (abs_float f3 < cutoff) && (abs_float (f3 *. abs_float expo_inv_q) < cutoff) then result else
let v3 =
vrr angMom_a cmm totAngMom_a (totAngMom_c-2)
in
Array.init (maxm+1) (fun m -> result.(m) +.
f3 *. (v3.(m) +. (if m=maxm then 0. else expo_inv_q *. v3.(m+1)) ))
in
let result =
if am.(xyz) lor cm.(xyz) < 0 then result else
let f5 =
(float_of_int angMom_a.(xyz)) *. expo_inv_p *. expo_inv_q *. 0.5
in
if (abs_float f5 < cutoff) then result else
let v5 =
vrr am cm (totAngMom_a-1) (totAngMom_c-1)
in
Array.init (maxm+1) (fun m ->
result.(m) -. (if m = maxm then 0. else f5 *. v5.(m+1)))
in
result
)
in
if not found then
Zmap.add map.(m) key result;
Zmap.add map key result;
result
(** Horizontal recurrence relations *)
and hrr0 m angMom_a angMom_b angMom_c
and hrr0 angMom_a angMom_b angMom_c
totAngMom_a totAngMom_b totAngMom_c =
match totAngMom_b with
| 0 -> vrr m angMom_a angMom_c totAngMom_a totAngMom_c
| 0 -> (vrr angMom_a angMom_c totAngMom_a totAngMom_c).(0)
| 1 -> let xyz = if angMom_b.(0) = 1 then 0 else if angMom_b.(1) = 1 then 1 else 2 in
let ap = [| angMom_a.(0) ; angMom_a.(1) ; angMom_a.(2) |] in
ap.(xyz) <- ap.(xyz) + 1;
vrr m ap angMom_c (totAngMom_a+1) totAngMom_c
+. chop (Coordinate.coord center_ab xyz) (fun () ->
vrr m angMom_a angMom_c totAngMom_a totAngMom_c)
let v1 =
vrr ap angMom_c (totAngMom_a+1) totAngMom_c
in
let f2 =
(Coordinate.coord center_ab xyz)
in
if (abs_float f2 < cutoff) then v1.(0) else
let v2 =
vrr angMom_a angMom_c totAngMom_a totAngMom_c
in
v1.(0) +. f2 *. v2.(0)
| _ ->
let ap = [| angMom_a.(0) ; angMom_a.(1) ; angMom_a.(2) |]
and bm = [| angMom_b.(0) ; angMom_b.(1) ; angMom_b.(2) |]
@ -134,18 +179,24 @@ let hvrr_two_e m (angMom_a, angMom_b, angMom_c, angMom_d)
ap.(xyz) <- ap.(xyz) + 1;
bm.(xyz) <- bm.(xyz) - 1;
if (bm.(xyz) < 0) then 0. else
hrr0 m ap bm angMom_c (totAngMom_a+1) (totAngMom_b-1) totAngMom_c
+. chop (Coordinate.coord center_ab xyz) (fun () ->
hrr0 m angMom_a bm angMom_c totAngMom_a (totAngMom_b-1)
totAngMom_c )
let h1 =
hrr0 ap bm angMom_c (totAngMom_a+1) (totAngMom_b-1) totAngMom_c
in
let f2 =
(Coordinate.coord center_ab xyz)
in
if (abs_float f2 < cutoff) then h1 else
let h2 =
hrr0 angMom_a bm angMom_c totAngMom_a (totAngMom_b-1) totAngMom_c
in
h1 +. f2 *. h2
and hrr m angMom_a angMom_b angMom_c angMom_d
and hrr angMom_a angMom_b angMom_c angMom_d
totAngMom_a totAngMom_b totAngMom_c totAngMom_d =
match (totAngMom_b, totAngMom_d) with
| (0,0) ->
vrr m angMom_a angMom_c totAngMom_a totAngMom_c
| (_,0) -> hrr0 m angMom_a angMom_b angMom_c totAngMom_a totAngMom_b totAngMom_c
| (0,0) -> (vrr angMom_a angMom_c totAngMom_a totAngMom_c).(0)
| (_,0) -> hrr0 angMom_a angMom_b angMom_c totAngMom_a totAngMom_b totAngMom_c
| (_,_) ->
let cp = [| angMom_c.(0) ; angMom_c.(1) ; angMom_c.(2) |]
and dm = [| angMom_d.(0) ; angMom_d.(1) ; angMom_d.(2) |]
@ -157,13 +208,17 @@ let hvrr_two_e m (angMom_a, angMom_b, angMom_c, angMom_d)
in
cp.(xyz) <- cp.(xyz) + 1;
dm.(xyz) <- dm.(xyz) - 1;
hrr m angMom_a angMom_b cp dm totAngMom_a totAngMom_b
(totAngMom_c+1) (totAngMom_d-1)
+. chop (Coordinate.coord center_cd xyz) (fun () ->
hrr m angMom_a angMom_b angMom_c dm totAngMom_a totAngMom_b
totAngMom_c (totAngMom_d-1) )
let h1 =
hrr angMom_a angMom_b cp dm totAngMom_a totAngMom_b (totAngMom_c+1) (totAngMom_d-1)
in
hrr m angMom_a angMom_b angMom_c angMom_d totAngMom_a totAngMom_b
let f2 = Coordinate.coord center_cd xyz in
if (abs_float f2 < cutoff) then h1 else
let h2 =
hrr angMom_a angMom_b angMom_c dm totAngMom_a totAngMom_b totAngMom_c (totAngMom_d-1)
in
h1 +. f2 *. h2
in
hrr angMom_a angMom_b angMom_c angMom_d totAngMom_a totAngMom_b
totAngMom_c totAngMom_d
@ -231,7 +286,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
contracted_class.(0) <- contracted_class.(0) +. coef_prod *. integral
| _ ->
let d = shell_q.(cd).Shell_pair.j in
let map = Array.init maxm (fun _ -> Zmap.create (Array.length class_indices)) in
let map = Zmap.create (Array.length class_indices) in
(* Compute the integral class from the primitive shell quartet *)
Array.iteri (fun i key ->
let a = Zkey.to_int_array Zkey.Kind_12 key in
@ -277,7 +332,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
shell_p.(ab).Shell_pair.norm_fun angMomA angMomB *. shell_q.(cd).Shell_pair.norm_fun angMomC angMomD
in
let integral = chop norm (fun () ->
hvrr_two_e 0 (angMomA, angMomB, angMomC, angMomD)
hvrr_two_e (angMomA, angMomB, angMomC, angMomD)
(Contracted_shell.totAngMom shell_a, Contracted_shell.totAngMom shell_b,
Contracted_shell.totAngMom shell_c, Contracted_shell.totAngMom shell_d)
(maxm, zero_m_array)

2
Makefile
View File

@ -42,4 +42,4 @@ clean:
rm -rf _build $(ALL_EXE) $(ALL_TESTS) *.native *.byte
debug: run_integrals.native
time ./run_integrals -c h2o.xyz -b ~/quantum_package/data/basis/cc-pvtz -o /dev/shm/out ; sleep 2 ; diff /dev/shm/out.eri REF | head -30
time ./run_integrals -c h2o.xyz -b ~/quantum_package/data/basis/cc-pvtz -o /dev/shm/out ; sleep 2 ; diff /dev/shm/out.eri REF | head -50