diff --git a/Basis/Contracted_shell.ml b/Basis/Contracted_shell.ml index 4fcb9f4..798fabd 100644 --- a/Basis/Contracted_shell.ml +++ b/Basis/Contracted_shell.ml @@ -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 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)) - in - result - ) s.expo + let factor int_array = + let dfa = Array.map (fun j -> + ( float_of_int (1 lsl j) *. fact j) /. fact (j+j) + ) int_array + in + 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) } diff --git a/Basis/ERI.ml b/Basis/ERI.ml index e1e4e69..c69d9a4 100644 --- a/Basis/ERI.ml +++ b/Basis/ERI.ml @@ -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 diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml index b301c84..773512d 100644 --- a/Basis/TwoElectronRR.ml +++ b/Basis/TwoElectronRR.ml @@ -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 + 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 m angMom_a angMom_b angMom_c angMom_d totAngMom_a totAngMom_b + 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) diff --git a/Makefile b/Makefile index b9b4652..71911b6 100644 --- a/Makefile +++ b/Makefile @@ -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