From 1a9f1ac54b602e678cff67d3a71e996f7390efdf Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 4 Apr 2019 00:28:29 +0200 Subject: [PATCH] Accelerated FCI --- CI/CI.ml | 105 +++++++++++++++++++++++++++++-------------------------- 1 file changed, 55 insertions(+), 50 deletions(-) diff --git a/CI/CI.ml b/CI/CI.ml index e3791cd..1cf7487 100644 --- a/CI/CI.ml +++ b/CI/CI.ml @@ -56,6 +56,15 @@ let h_ij mo_basis ki kj = |> List.hd +let h_ij_non_zero mo_basis deg_a deg_b ki kj = + let integrals = + List.map (fun f -> f mo_basis) + [ h_integrals ] + in + CIMatrixElement.non_zero integrals deg_a deg_b ki kj + |> List.hd + + let create_matrix_arbitrary f det_space = lazy ( let ndet = Ds.size det_space in @@ -106,8 +115,8 @@ let create_matrix_arbitrary f det_space = in (** Update function when ki and kj are connected *) - let update i j ki kj = - let x = f ki kj in + let update i j deg_a deg_b ki kj = + let x = f deg_a deg_b ki kj in if abs_float x > Constants.epsilon then result.(i - shift) <- (j, x) :: result.(i - shift) ; in @@ -131,7 +140,7 @@ let create_matrix_arbitrary f det_space = let ki = Determinant.of_spindeterminants i_alfa i_beta in let kj = Determinant.of_spindeterminants j_alfa i_beta in update (index_start.(i) + i') - (index_start.(j) + j' + 1) ki kj); + (index_start.(j) + j' + 1) 2 0 ki kj); raise Exit) ) j_dets with Exit -> () @@ -153,7 +162,7 @@ let create_matrix_arbitrary f det_space = Determinant.of_spindeterminants j_alfa j_beta in (update (index_start.(i) + i') (index_start.(j) + j' + 1) - ki kj; + 1 (Determinant.degree_beta ki kj) ki kj; aux r_singles (j'+1);) | 1 -> if (j' < Array.length j_dets) then aux singles (j'+1) | _ -> assert false @@ -177,7 +186,7 @@ let create_matrix_arbitrary f det_space = Determinant.of_spindeterminants j_alfa j_beta in (update (index_start.(i) + i') (index_start.(j) + j' + 1) - ki kj; + 0 (Determinant.degree_beta ki kj) ki kj; aux r_doubles (j'+1);) | 1 -> if (j' < Array.length j_dets) then aux doubles (j'+1) | _ -> assert false @@ -268,8 +277,8 @@ let create_matrix_spin f det_space = in (** Update function when ki and kj are connected *) - let update i j ki kj = - let x = f ki kj in + let update i j deg_a deg_b ki kj = + let x = f deg_a deg_b ki kj in if abs_float x > Constants.epsilon then result.(i) <- (j, x) :: result.(i) ; in @@ -285,7 +294,7 @@ let create_matrix_spin f det_space = List.iteri (fun ib i_beta -> let ki = Determinant.of_spindeterminants i_alfa i_beta in let kj = Determinant.of_spindeterminants j_alfa i_beta in - update !i' (ib + !j) ki kj; + update !i' (ib + !j) 2 0 ki kj; incr i'; ) b; | 1 -> @@ -295,7 +304,7 @@ let create_matrix_spin f det_space = let singles, _ = degree_bb.(ib) in List.iter (fun (j', j_beta) -> let kj = Determinant.of_spindeterminants j_alfa j_beta in - update !i' (j' + !j) ki kj + update !i' (j' + !j) 1 (Determinant.degree_beta ki kj) ki kj ) singles; incr i'; ) b; @@ -306,7 +315,7 @@ let create_matrix_spin f det_space = let _singles, doubles = degree_bb.(ib) in List.iter (fun (j', j_beta) -> let kj = Determinant.of_spindeterminants j_alfa j_beta in - update !i' (j' + !j) ki kj + update !i' (j' + !j) 0 (Determinant.degree_beta ki kj) ki kj ) doubles; incr i'; ) b; @@ -337,15 +346,6 @@ let create_matrix_spin f det_space = ) -type mdata = - { i_a : int ref; - j_a : int ref; - j0 : int ref; - j_alfa_prev : int ref; - bi : int; - h1 : (int -> int -> int -> int -> float) ref - h123_prev : (int -> float) ref - } (* Create a matrix using the fact that the determinant space is made of the outer product of spindeterminants. *) @@ -361,31 +361,39 @@ let create_matrix_spin_computed f det_space = let n_alfa = Array.length a in let h i_alfa j_alfa = - match Spindeterminant.degree a.(i_alfa) a.(j_alfa) with + let deg_a = Spindeterminant.degree a.(i_alfa) a.(j_alfa) in + match deg_a with | 2 -> let ai, aj = a.(i_alfa), a.(j_alfa) in (fun i_beta j_beta -> if i_beta <> j_beta then 0. else - let ki = Determinant.of_spindeterminants ai b.(i_beta) in - let kj = Determinant.of_spindeterminants aj b.(j_beta) in - f ki kj + let deg_b = 0 in + let ki = Determinant.of_spindeterminants ai b.(i_beta) in + let kj = Determinant.of_spindeterminants aj b.(j_beta) in + f deg_a deg_b ki kj ) | 1 -> let ai, aj = a.(i_alfa), a.(j_alfa) in (fun i_beta j_beta -> - match Spindeterminant.degree b.(i_beta) b.(j_beta) with + let deg_b = Spindeterminant.degree b.(i_beta) b.(j_beta) in + match deg_b with | 0 | 1 -> let ki = Determinant.of_spindeterminants ai b.(i_beta) in let kj = Determinant.of_spindeterminants aj b.(j_beta) in - f ki kj + f deg_a deg_b ki kj | _ -> 0. ) | 0 -> let ai, aj = a.(i_alfa), a.(j_alfa) in (fun i_beta j_beta -> - let ki = Determinant.of_spindeterminants ai b.(i_beta) in - let kj = Determinant.of_spindeterminants aj b.(j_beta) in - f ki kj + let deg_b = Spindeterminant.degree b.(i_beta) b.(j_beta) in + match deg_b with + | 0 | 1 | 2 -> + let deg_b = Spindeterminant.degree b.(i_beta) b.(j_beta) in + let ki = Determinant.of_spindeterminants ai b.(i_beta) in + let kj = Determinant.of_spindeterminants aj b.(j_beta) in + f deg_a deg_b ki kj + | _ -> 0. ) | _ -> (fun _ _ -> 0.) in @@ -398,26 +406,24 @@ let create_matrix_spin_computed f det_space = if i <> !i_prev then begin i_prev := i; - let j_a = ref (-n_alfa) in - let j0 = ref (!j_a * n_beta) in let i_a = (i-1)/n_beta in let i_alfa = i_a + 1 in let h1 = h (i_alfa-1) in - let bi = - let i_beta = i - i_a*n_beta in - i_beta-1 - in - let j_alfa_prev = ref (-10) in + let i_beta = i - i_a*n_beta in + let bi = (i_beta-1) in let h123_prev = ref (fun _ -> 0.) in - let f j = - if j > !j0 + n_beta || j < !j0 then - begin - let x = (j-1)/n_beta in - j_a := x; - j0 := x * n_beta - end; + let j_a = ref (-n_alfa) in + let j_alfa_prev = ref (-10) in + let j0 = ref (!j_a * n_beta) in + result := fun j -> + if j > !j0 + n_beta + || j < !j0 + then begin + j_a := (j-1)/n_beta; + j0 := !j_a * n_beta + end; let j_alfa = !j_a + 1 in let h123 = if j_alfa <> !j_alfa_prev then @@ -429,8 +435,6 @@ let create_matrix_spin_computed f det_space = in let j_beta = j - !j0 in h123 (j_beta-1) - in - result := f ; end; !result in @@ -449,7 +453,7 @@ let make ?(n_states=1) ?(algo=`Direct) det_space = Ds.determinant_stream det_space |> Stream.next in - h_ij mo_basis d0 d0 + h_ij_non_zero mo_basis 0 0 d0 d0 in let m_H = @@ -468,11 +472,11 @@ let make ?(n_states=1) ?(algo=`Direct) det_space = else create_matrix_spin in - f (fun ki kj -> - if ki <> kj then - h_ij mo_basis ki kj + f (fun deg_a deg_b ki kj -> + if deg_a + deg_b > 0 then + h_ij_non_zero mo_basis deg_a deg_b ki kj else - h_ij mo_basis ki ki -. e_shift + h_ij_non_zero mo_basis 0 0 ki ki -. e_shift ) det_space in @@ -482,7 +486,8 @@ let make ?(n_states=1) ?(algo=`Direct) det_space = | Ds.Arbitrary _ -> create_matrix_arbitrary | Ds.Spin _ -> create_matrix_spin in - f (fun ki kj -> CIMatrixElement.make_s2 ki kj) det_space + f (fun _deg_a _deg_b ki kj -> CIMatrixElement.make_s2 ki kj) det_space + (*TODO*) in let eigensystem = lazy (