From ac45b0b0cb511abd87780529ad21b38eaba59eb5 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 28 Feb 2019 19:37:54 +0100 Subject: [PATCH] Cleaning --- CI/CI.ml | 212 ++++++++++++++++++++++++++++++------------------------- 1 file changed, 115 insertions(+), 97 deletions(-) diff --git a/CI/CI.ml b/CI/CI.ml index ef35da9..f670f53 100644 --- a/CI/CI.ml +++ b/CI/CI.ml @@ -51,47 +51,61 @@ let h_ij mo_basis ki kj = |> List.hd +let create_matrix_arbitrary f det_space = + lazy ( + let det = + match Ds.determinants det_space with + | Ds.Arbitrary a -> a + | _ -> assert false + in -let make ?(n_states=1) det_space = - let ndet = Ds.size det_space in - let mo_basis = Ds.mo_basis det_space in + let ndet = Ds.size det_space in - let m_H = - let m_H_arbitrary det = lazy ( - let v = Vec.make0 ndet in - Array.init ndet (fun i -> let ki = det.(i) in + let v = Vec.make0 ndet in + + Array.init ndet + (fun i -> let ki = det.(i) in Printf.eprintf "%8d / %8d\r%!" i ndet; let j = ref 1 in Ds.determinant_stream det_space - |> Stream.iter (fun kj -> - v.{!j} <- h_ij mo_basis ki kj ; incr j); + |> Stream.iter (fun kj -> v.{!j} <- f ki kj ; incr j); Vector.sparse_of_vec v) - |> Matrix.sparse_of_vector_array) + |> Matrix.sparse_of_vector_array + ) + + +(* Create a matrix using the fact that the determinant space is made of + the outer product of spindeterminants. *) +let create_matrix_spin f det_space = + lazy ( + let ndet = Ds.size det_space in + let a, b = + match Ds.determinants det_space with + | Ds.Spin (a,b) -> (a,b) + | _ -> assert false + in + let n_alfa = Array.length a in + let n_beta = Array.length b in + let result = Array.init ndet (fun _ -> []) in + + (** Update function when ki and kj are connected *) + let update i j ki kj = + let x = f ki kj in + if x <> 0. then + result.(i) <- (j, x) :: result.(i) ; in - let m_H_spin a b = lazy ( - let n_alfa = Array.length a in - let n_beta = Array.length b in - let result = Array.init ndet (fun _ -> []) in - - (** Update function when ki and kj are connected *) - let update i j ki kj = - let x = h_ij mo_basis ki kj in - if x <> 0. then - result.(i) <- (j, x) :: result.(i) ; - in - - (** Array of (list of singles, list of doubles) in the beta spin *) - let degree_bb = - Array.map (fun det_i -> + (** Array of (list of singles, list of doubles) in the beta spin *) + let degree_bb = + Array.map (fun det_i -> let deg = Spindeterminant.degree det_i in let doubles = Array.mapi (fun i det_j -> - let d = deg det_j in - if d < 3 then - Some (i,d,det_j) - else - None + let d = deg det_j in + if d < 3 then + Some (i,d,det_j) + else + None ) b |> Array.to_list |> Util.list_some @@ -105,79 +119,83 @@ let make ?(n_states=1) det_space = in (singles, doubles) ) b - in - let a = Array.to_list a - and b = Array.to_list b - in - let i = ref 0 in - List.iteri (fun ia i_alfa -> - Printf.eprintf "%8d / %8d\r%!" ia n_alfa; - let j = ref 1 in - let deg_a = Spindeterminant.degree i_alfa in - List.iter (fun j_alfa -> - let degree_a = deg_a j_alfa in - begin - match degree_a with - | 2 -> - let i' = ref !i in - 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; - incr i'; - ) b; - | 1 -> - let i' = ref !i in - List.iteri (fun ib i_beta -> - let ki = Determinant.of_spindeterminants i_alfa i_beta in - 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 - ) singles; - incr i'; - ) b; - | 0 -> - let i' = ref !i in - List.iteri (fun ib i_beta -> - let ki = Determinant.of_spindeterminants i_alfa i_beta in - 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 - ) doubles; - incr i'; - ) b; - | _ -> (); - end; - j := !j + n_beta - ) a; - i := !i + n_beta - ) a; + in + let a = Array.to_list a + and b = Array.to_list b + in + let i = ref 0 in + List.iteri (fun ia i_alfa -> + Printf.eprintf "%8d / %8d\r%!" ia n_alfa; + let j = ref 1 in + let deg_a = Spindeterminant.degree i_alfa in + List.iter (fun j_alfa -> + let degree_a = deg_a j_alfa in + begin + match degree_a with + | 2 -> + let i' = ref !i in + 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; + incr i'; + ) b; + | 1 -> + let i' = ref !i in + List.iteri (fun ib i_beta -> + let ki = Determinant.of_spindeterminants i_alfa i_beta in + 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 + ) singles; + incr i'; + ) b; + | 0 -> + let i' = ref !i in + List.iteri (fun ib i_beta -> + let ki = Determinant.of_spindeterminants i_alfa i_beta in + 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 + ) doubles; + incr i'; + ) b; + | _ -> (); + end; + j := !j + n_beta + ) a; + i := !i + n_beta + ) a; - Array.map (fun l -> + Array.map (fun l -> List.sort compare l |> Vector.sparse_of_assoc_list ndet ) result - |> Matrix.sparse_of_vector_array - ) - in + |> Matrix.sparse_of_vector_array + ) - match Ds.determinants det_space with - | Ds.Arbitrary a -> m_H_arbitrary a - | Ds.Spin (a,b) -> m_H_spin a b + + +let make ?(n_states=1) det_space = + + let m_H = + let mo_basis = Ds.mo_basis det_space in + let f = + match Ds.determinants det_space with + | Ds.Arbitrary _ -> create_matrix_arbitrary + | Ds.Spin _ -> create_matrix_spin + in + f (fun ki kj -> h_ij mo_basis ki kj) det_space in - let m_S2 = lazy ( - match Ds.determinants det_space with - | Ds.Spin (a,b) -> failwith "Not implemented" - | Ds.Arbitrary det -> - let v = Vec.make0 ndet in - Array.init ndet (fun i -> let ki = det.(i) in - Array.iteri (fun j kj -> - v.{j+1} <- CIMatrixElement.make_s2 ki kj) det; - Vector.sparse_of_vec v) - |> Matrix.sparse_of_vector_array - ) + let m_S2 = + let f = + match Ds.determinants det_space with + | Ds.Arbitrary _ -> create_matrix_arbitrary + | Ds.Spin _ -> create_matrix_spin + in + f (fun ki kj -> CIMatrixElement.make_s2 ki kj) det_space in let eigensystem = lazy (