diff --git a/Basis/Contracted_shell.ml b/Basis/Contracted_shell.ml index 9f54e46..4fcb9f4 100644 --- a/Basis/Contracted_shell.ml +++ b/Basis/Contracted_shell.ml @@ -45,14 +45,17 @@ let compute_norm_coef s = Angular_momentum.to_int s.totAngMom in Array.mapi (fun i alpha -> + let alpha_2 = + alpha +. alpha + in let c = - ((alpha +. alpha) *. pi_inv)**(1.5) *. (pow (4. *. alpha) atot) + (alpha_2 *. pi_inv)**(1.5) *. (pow (alpha_2 +. alpha_2) atot) in let result a = let dfa = Array.map (fun j -> - fact (j+j) /. ( float_of_int (1 lsl j) *. fact j) + ( float_of_int (1 lsl j) *. fact j) /. fact (j+j) ) a - in sqrt (c /. (dfa.(0) *.dfa.(1) *. dfa.(2))) + in sqrt (c *. dfa.(0) *.dfa.(1) *. dfa.(2)) in result ) s.expo diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml index 766bcfa..08ff18b 100644 --- a/Basis/TwoElectronRR.ml +++ b/Basis/TwoElectronRR.ml @@ -22,11 +22,11 @@ let hvrr_two_e m (angMom_a, angMom_b, angMom_c, angMom_d) (** Vertical recurrence relations *) let rec vrr0 m angMom_a = function - | 0 -> zero_m_array.(m) | 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) | totAngMom_a -> let key = Zkey.of_int_tuple (Zkey.Three (angMom_a.(0)+1, angMom_a.(1)+1, angMom_a.(2)+1) ) @@ -196,74 +196,40 @@ 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 *) - begin - match Contracted_shell.(totAngMom shell_a, totAngMom shell_b, - totAngMom shell_c, totAngMom shell_d) with - | Angular_momentum.(S,S,S,S) -> - begin - for ab=0 to (Array.length shell_p - 1) do - for cd=0 to (Array.length shell_q - 1) do - let coef_prod = - shell_p.(ab).Shell_pair.coef *. shell_q.(cd).Shell_pair.coef - in - (** Screening on the product of coefficients *) - try - if (abs_float coef_prod) < 1.e-4*.cutoff then - raise NullQuartet; + for ab=0 to (Array.length shell_p - 1) do + let cab = shell_p.(ab).Shell_pair.coef in + let b = shell_p.(ab).Shell_pair.j in + for cd=0 to (Array.length shell_q - 1) do + let coef_prod = + cab *. shell_q.(cd).Shell_pair.coef + in + (** Screening on the product of coefficients *) + try + if (abs_float coef_prod) < 1.e-4*.cutoff then + raise NullQuartet; - let expo_pq_inv = - shell_p.(ab).Shell_pair.expo_inv +. shell_q.(cd).Shell_pair.expo_inv - in - let center_pq = - Coordinate.(shell_p.(ab).Shell_pair.center |- shell_q.(cd).Shell_pair.center) - in - let norm_pq_sq = - Coordinate.dot center_pq center_pq - in + let expo_pq_inv = + shell_p.(ab).Shell_pair.expo_inv +. shell_q.(cd).Shell_pair.expo_inv + in + let center_pq = + Coordinate.(shell_p.(ab).Shell_pair.center |- shell_q.(cd).Shell_pair.center) + in + let norm_pq_sq = + Coordinate.dot center_pq center_pq + in - let zero_m_array = - zero_m ~maxm ~expo_pq_inv ~norm_pq_sq - in - - let coef_prod = - shell_p.(ab).Shell_pair.coef *. shell_q.(cd).Shell_pair.coef - in + let zero_m_array = + zero_m ~maxm ~expo_pq_inv ~norm_pq_sq + in + begin + match Contracted_shell.(totAngMom shell_a, totAngMom shell_b, + totAngMom shell_c, totAngMom shell_d) with + | Angular_momentum.(S,S,S,S) -> let integral = zero_m_array.(0) in contracted_class.(0) <- contracted_class.(0) +. coef_prod *. integral - with NullQuartet -> () - done - done; - end - - | _ -> - begin - for ab=0 to (Array.length shell_p - 1) do - let b = shell_p.(ab).Shell_pair.j in - for cd=0 to (Array.length shell_q - 1) do - try - let coef_prod = - shell_p.(ab).Shell_pair.coef *. shell_q.(cd).Shell_pair.coef - in - (** Screening on the product of coefficients *) - if (abs_float coef_prod) < 1.e-4*.cutoff then - raise NullQuartet; - - let expo_pq_inv = - shell_p.(ab).Shell_pair.expo_inv +. shell_q.(cd).Shell_pair.expo_inv - in - let center_pq = - Coordinate.(shell_p.(ab).Shell_pair.center |- shell_q.(cd).Shell_pair.center) - in - let norm_pq_sq = - Coordinate.dot center_pq center_pq - in - - let zero_m_array = - zero_m ~maxm ~expo_pq_inv ~norm_pq_sq - in - + | _ -> let d = shell_q.(cd).Shell_pair.j in let map = Array.init maxm (fun _ -> Zmap.create (Array.length class_indices)) in (* Compute the integral class from the primitive shell quartet *) @@ -323,11 +289,10 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral with NullQuartet -> () ) class_indices - with NullQuartet -> () - done - done; - end - end; + end + with NullQuartet -> () + done + done; let result = Zmap.create (Array.length contracted_class) diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml index d148bc0..21d1331 100644 --- a/Basis/TwoElectronRRVectorized.ml +++ b/Basis/TwoElectronRRVectorized.ml @@ -27,7 +27,6 @@ let hvrr_two_e_vector m (angMom_a, angMom_b, angMom_c, angMom_d) (** Vertical recurrence relations *) let rec vrr0_v m angMom_a = function - | 0 -> Array.mapi (fun k c -> c *. zero_m_array.(k).(m)) coef_prod | 1 -> let i = if angMom_a.(0) = 1 then 0 else if angMom_a.(1) = 1 then 1 else 2 in let f = expo_b *. (Coordinate.coord center_ab i) in @@ -35,6 +34,7 @@ let hvrr_two_e_vector m (angMom_a, angMom_b, angMom_c, angMom_d) ( (Coordinate.coord center_pq.(k) i) *. zero_m_array.(k).(m+1) -. f *. zero_m_array.(k).(m) ) ) coef_prod + | 0 -> Array.mapi (fun k c -> c *. zero_m_array.(k).(m)) coef_prod | totAngMom_a -> let key = Zkey.of_int_tuple (Zkey.Three (angMom_a.(0)+1, angMom_a.(1)+1, angMom_a.(2)+1) ) diff --git a/run_integrals.ml b/run_integrals.ml index 77be74a..96c1950 100644 --- a/run_integrals.ml +++ b/run_integrals.ml @@ -11,6 +11,10 @@ let speclist = [ ] let run ~out = + (* + let gc = Gc.get () in + Gc.set { gc with minor_heap_size=(262144 / 16) }; + *) let out_file = match out with | None -> raise (Invalid_argument "Output file should be specified with -o")