diff --git a/CI/F12CI.ml b/CI/F12CI.ml index 2a19acf..afa9506 100644 --- a/CI/F12CI.ml +++ b/CI/F12CI.ml @@ -70,7 +70,7 @@ let single_matrices hf12_integrals density = done; *) -let hf_ij_non_zero hf12_integrals deg_a deg_b ki kj = +let hf_ij_non_zero mo_basis hf12_integrals deg_a deg_b ki kj = let integrals = [ let { HF12. @@ -85,26 +85,20 @@ let hf_ij_non_zero hf12_integrals deg_a deg_b ki kj = let mo_a = Bitstring.logand (Sp.bitstring kia) (Sp.bitstring kja) |> Bitstring.to_list + |> Array.of_list and mo_b = Bitstring.logand (Sp.bitstring kib) (Sp.bitstring kjb) |> Bitstring.to_list + |> Array.of_list + in + + let aux_mos = + Util.list_range (MOBasis.size mo_basis) (MOBasis.size aux_basis) + |> Array.of_list in let one_e _ _ _ = 0. in - let two_e i j k l s s' = - if s = s' then - hf12_anti.{i,j,k,l} -. ( - (List.fold_left (fun accu m -> accu +. hf12_single_anti.{m,i,j,k,l}) 0. mo_a) +. - (List.fold_left (fun accu m -> accu +. hf12_single_anti.{m,j,i,l,k}) 0. mo_b) - ) - else - hf12.{i,j,k,l} -. ( - (List.fold_left (fun accu m -> accu +. hf12_single.{m,i,j,k,l}) 0. mo_a) +. - (List.fold_left (fun accu m -> accu +. hf12_single.{m,j,i,l,k}) 0. mo_b) - ) - in - let h = MOBasis.ee_ints aux_basis in let two_e_h i j k l s s' = if s' <> s then @@ -115,24 +109,48 @@ let hf_ij_non_zero hf12_integrals deg_a deg_b ki kj = let f = MOBasis.f12_ints aux_basis in let two_e_f i j k l s s' = + let fijkl = F12.get_phys f i j k l + and fijlk = F12.get_phys f i j l k + in if s' <> s then - F12.get_phys f i j k l + 0.325 *. fijkl +. 0.125 *. fijlk else - (F12.get_phys f i j k l) -. (F12.get_phys f i j l k) + 0.25 *. (fijkl -. fijlk) in + (* let mo_of_s = function | Spin.Alfa -> mo_a | Spin.Beta -> mo_b in + *) + + let two_e i j k l s s' = + ( + if s = s' then + hf12_anti.{i,j,k,l} -. ( + (Array.fold_left (fun accu m -> accu +. hf12_single_anti.{m,i,j,k,l}) 0. mo_a) +. + (Array.fold_left (fun accu m -> accu +. hf12_single_anti.{m,j,i,l,k}) 0. mo_b) ) + else + hf12.{i,j,k,l} -. ( + (Array.fold_left (fun accu m -> accu +. hf12_single.{m,i,j,k,l}) 0. mo_a) +. + (Array.fold_left (fun accu m -> accu +. hf12_single.{m,j,i,l,k}) 0. mo_b) ) + ) + (* + +. Array.fold_left ( fun accu a -> accu +. + (Array.fold_left ( fun accu m -> accu +. two_e_h m i m a s s) 0. (mo_of_s s) +. + Array.fold_left ( fun accu m -> accu +. two_e_h m i m a s s) 0. (mo_of_s @@ Spin.other s) ) + *. (two_e_f a j k l s s') ) 0. aux_mos + *) + in let three_e i j k l m n s s' s'' = - List.fold_left (fun accu a -> accu +. two_e_h i j l a s s' *. two_e_f a k m n s' s'') 0. (mo_of_s s' ) - -. List.fold_left (fun accu a -> accu +. two_e_h j i m a s' s *. two_e_f a k l n s s'') 0. (mo_of_s s ) - +. List.fold_left (fun accu a -> accu +. two_e_h j k m a s' s'' *. two_e_f a i n l s'' s ) 0. (mo_of_s s'') - -. List.fold_left (fun accu a -> accu +. two_e_h k j n a s'' s' *. two_e_f a i m l s' s ) 0. (mo_of_s s' ) - +. List.fold_left (fun accu a -> accu +. two_e_h k i n a s'' s *. two_e_f a j l m s s' ) 0. (mo_of_s s ) - -. List.fold_left (fun accu a -> accu +. two_e_h i k l a s s'' *. two_e_f a j n m s'' s' ) 0. (mo_of_s s'') + Array.fold_left (fun accu a -> accu +. two_e_h i j l a s s' *. two_e_f a k m n s' s'') 0. aux_mos + -. Array.fold_left (fun accu a -> accu +. two_e_h j i m a s' s *. two_e_f a k l n s s'') 0. aux_mos + +. Array.fold_left (fun accu a -> accu +. two_e_h j k m a s' s'' *. two_e_f a i n l s'' s ) 0. aux_mos + -. Array.fold_left (fun accu a -> accu +. two_e_h k j n a s'' s' *. two_e_f a i m l s' s ) 0. aux_mos + +. Array.fold_left (fun accu a -> accu +. two_e_h k i n a s'' s *. two_e_f a j l m s s' ) 0. aux_mos + -. Array.fold_left (fun accu a -> accu +. two_e_h i k l a s s'' *. two_e_f a j n m s'' s' ) 0. aux_mos in (one_e, two_e, Some three_e) ] @@ -148,6 +166,9 @@ let dressing_vector ~frozen_core hf12_integrals f12_amplitudes ci = let det_space = ci.CI.det_space in + let mo_basis = + Ds.mo_basis det_space + in let m_HF = @@ -157,7 +178,7 @@ let dressing_vector ~frozen_core hf12_integrals f12_amplitudes ci = | Ds.Spin _ -> CI.create_matrix_spin_computed ~nmax:3 in f (fun deg_a deg_b ki kj -> - hf_ij_non_zero hf12_integrals deg_a deg_b ki kj + hf_ij_non_zero mo_basis hf12_integrals deg_a deg_b ki kj ) det_space in