10
1
mirror of https://gitlab.com/scemama/QCaml.git synced 2024-06-02 03:15:19 +02:00
This commit is contained in:
Anthony Scemama 2019-02-28 19:37:54 +01:00
parent bb677d4f3a
commit ac45b0b0cb

212
CI/CI.ml
View File

@ -51,47 +51,61 @@ let h_ij mo_basis ki kj =
|> List.hd |> 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 ndet = Ds.size det_space in
let mo_basis = Ds.mo_basis det_space in
let m_H = let v = Vec.make0 ndet in
let m_H_arbitrary det = lazy (
let v = Vec.make0 ndet in Array.init ndet
Array.init ndet (fun i -> let ki = det.(i) in (fun i -> let ki = det.(i) in
Printf.eprintf "%8d / %8d\r%!" i ndet; Printf.eprintf "%8d / %8d\r%!" i ndet;
let j = ref 1 in let j = ref 1 in
Ds.determinant_stream det_space Ds.determinant_stream det_space
|> Stream.iter (fun kj -> |> Stream.iter (fun kj -> v.{!j} <- f ki kj ; incr j);
v.{!j} <- h_ij mo_basis ki kj ; incr j);
Vector.sparse_of_vec v) 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 in
let m_H_spin a b = lazy ( (** Array of (list of singles, list of doubles) in the beta spin *)
let n_alfa = Array.length a in let degree_bb =
let n_beta = Array.length b in Array.map (fun det_i ->
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 ->
let deg = Spindeterminant.degree det_i in let deg = Spindeterminant.degree det_i in
let doubles = let doubles =
Array.mapi (fun i det_j -> Array.mapi (fun i det_j ->
let d = deg det_j in let d = deg det_j in
if d < 3 then if d < 3 then
Some (i,d,det_j) Some (i,d,det_j)
else else
None None
) b ) b
|> Array.to_list |> Array.to_list
|> Util.list_some |> Util.list_some
@ -105,79 +119,83 @@ let make ?(n_states=1) det_space =
in in
(singles, doubles) (singles, doubles)
) b ) b
in in
let a = Array.to_list a let a = Array.to_list a
and b = Array.to_list b and b = Array.to_list b
in in
let i = ref 0 in let i = ref 0 in
List.iteri (fun ia i_alfa -> List.iteri (fun ia i_alfa ->
Printf.eprintf "%8d / %8d\r%!" ia n_alfa; Printf.eprintf "%8d / %8d\r%!" ia n_alfa;
let j = ref 1 in let j = ref 1 in
let deg_a = Spindeterminant.degree i_alfa in let deg_a = Spindeterminant.degree i_alfa in
List.iter (fun j_alfa -> List.iter (fun j_alfa ->
let degree_a = deg_a j_alfa in let degree_a = deg_a j_alfa in
begin begin
match degree_a with match degree_a with
| 2 -> | 2 ->
let i' = ref !i in let i' = ref !i in
List.iteri (fun ib i_beta -> List.iteri (fun ib i_beta ->
let ki = Determinant.of_spindeterminants i_alfa i_beta in let ki = Determinant.of_spindeterminants i_alfa i_beta in
let kj = Determinant.of_spindeterminants j_alfa i_beta in let kj = Determinant.of_spindeterminants j_alfa i_beta in
update !i' (ib + !j) ki kj; update !i' (ib + !j) ki kj;
incr i'; incr i';
) b; ) b;
| 1 -> | 1 ->
let i' = ref !i in let i' = ref !i in
List.iteri (fun ib i_beta -> List.iteri (fun ib i_beta ->
let ki = Determinant.of_spindeterminants i_alfa i_beta in let ki = Determinant.of_spindeterminants i_alfa i_beta in
let singles, _ = degree_bb.(ib) in let singles, _ = degree_bb.(ib) in
List.iter (fun (j', j_beta) -> List.iter (fun (j', j_beta) ->
let kj = Determinant.of_spindeterminants j_alfa j_beta in let kj = Determinant.of_spindeterminants j_alfa j_beta in
update !i' (j' + !j) ki kj update !i' (j' + !j) ki kj
) singles; ) singles;
incr i'; incr i';
) b; ) b;
| 0 -> | 0 ->
let i' = ref !i in let i' = ref !i in
List.iteri (fun ib i_beta -> List.iteri (fun ib i_beta ->
let ki = Determinant.of_spindeterminants i_alfa i_beta in let ki = Determinant.of_spindeterminants i_alfa i_beta in
let _singles, doubles = degree_bb.(ib) in let _singles, doubles = degree_bb.(ib) in
List.iter (fun (j', j_beta) -> List.iter (fun (j', j_beta) ->
let kj = Determinant.of_spindeterminants j_alfa j_beta in let kj = Determinant.of_spindeterminants j_alfa j_beta in
update !i' (j' + !j) ki kj update !i' (j' + !j) ki kj
) doubles; ) doubles;
incr i'; incr i';
) b; ) b;
| _ -> (); | _ -> ();
end; end;
j := !j + n_beta j := !j + n_beta
) a; ) a;
i := !i + n_beta i := !i + n_beta
) a; ) a;
Array.map (fun l -> Array.map (fun l ->
List.sort compare l List.sort compare l
|> Vector.sparse_of_assoc_list ndet ) result |> Vector.sparse_of_assoc_list ndet ) result
|> Matrix.sparse_of_vector_array |> Matrix.sparse_of_vector_array
) )
in
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 in
let m_S2 = lazy ( let m_S2 =
match Ds.determinants det_space with let f =
| Ds.Spin (a,b) -> failwith "Not implemented" match Ds.determinants det_space with
| Ds.Arbitrary det -> | Ds.Arbitrary _ -> create_matrix_arbitrary
let v = Vec.make0 ndet in | Ds.Spin _ -> create_matrix_spin
Array.init ndet (fun i -> let ki = det.(i) in in
Array.iteri (fun j kj -> f (fun ki kj -> CIMatrixElement.make_s2 ki kj) det_space
v.{j+1} <- CIMatrixElement.make_s2 ki kj) det;
Vector.sparse_of_vec v)
|> Matrix.sparse_of_vector_array
)
in in
let eigensystem = lazy ( let eigensystem = lazy (