diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml index 346615f..07aad6f 100644 --- a/Basis/TwoElectronRRVectorized.ml +++ b/Basis/TwoElectronRRVectorized.ml @@ -558,47 +558,47 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q match Contracted_shell.(totAngMom shell_a, totAngMom shell_b, totAngMom shell_c, totAngMom shell_d) with | Angular_momentum.(S,S,S,S) -> - contracted_class.(0) <- - begin + contracted_class.(0) <- + begin try - let expo_inv_p = - Vec.init np (fun ab -> sp.(ab-1).ShellPair.expo_inv) - and expo_inv_q = - Vec.init nq (fun cd -> sq.(cd-1).ShellPair.expo_inv) - in + let expo_inv_p = + Vec.init np (fun ab -> sp.(ab-1).ShellPair.expo_inv) + and expo_inv_q = + Vec.init nq (fun cd -> sq.(cd-1).ShellPair.expo_inv) + in - let coef = - let result = Mat.make0 nq np in - Lacaml.D.ger - (Vec.of_array @@ filter_q shell_q.ContractedShellPair.coef) - (Vec.of_array @@ filter_p shell_p.ContractedShellPair.coef) - result; - result - in - let zm_array = Mat.init_cols np nq (fun i j -> - try - if (abs_float coef.{j,i} ) < 1.e-3*.cutoff then - raise NullQuartet; + let coef = + let result = Mat.make0 nq np in + Lacaml.D.ger + (Vec.of_array @@ filter_q shell_q.ContractedShellPair.coef) + (Vec.of_array @@ filter_p shell_p.ContractedShellPair.coef) + result; + result + in + let zm_array = Mat.init_cols np nq (fun i j -> + try + if (abs_float coef.{j,i} ) < 1.e-3*.cutoff then + raise NullQuartet; - let expo_pq_inv = - expo_inv_p.{i} +. expo_inv_q.{j} - in - let center_pq = - Coordinate.(sp.(i-1).ShellPair.center |- sq.(j-1).ShellPair.center) - in - let norm_pq_sq = - Coordinate.dot center_pq center_pq - in + let expo_pq_inv = + expo_inv_p.{i} +. expo_inv_q.{j} + in + let center_pq = + Coordinate.(sp.(i-1).ShellPair.center |- sq.(j-1).ShellPair.center) + in + let norm_pq_sq = + Coordinate.dot center_pq center_pq + in - let zero_m_array = - zero_m ~maxm:0 ~expo_pq_inv ~norm_pq_sq - in - zero_m_array.(0) - with NullQuartet -> 0. - ) in - Mat.gemm_trace zm_array coef + let zero_m_array = + zero_m ~maxm:0 ~expo_pq_inv ~norm_pq_sq + in + zero_m_array.(0) + with NullQuartet -> 0. + ) in + Mat.gemm_trace zm_array coef with (Invalid_argument _) -> 0. - end + end | _ -> let coef =