LCPQ/QCaml
10
1
Fork 0
mirror of https://gitlab.com/scemama/QCaml.git synced 2024-06-01 19:05:18 +02:00
Code Releases Activity

Accelerated FCI

Browse Source
This commit is contained in:
Anthony Scemama 2019-04-04 00:28:29 +02:00
parent 34033829e4
commit 1a9f1ac54b
1 changed files with 55 additions and 50 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

105
CI/CI.ml
View File

@ -56,6 +56,15 @@ let h_ij mo_basis ki kj =
|> List.hd |> 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 = let create_matrix_arbitrary f det_space =
lazy ( lazy (
let ndet = Ds.size det_space in let ndet = Ds.size det_space in
@ -106,8 +115,8 @@ let create_matrix_arbitrary f det_space =
in in
(** Update function when ki and kj are connected *) (** Update function when ki and kj are connected *)
let update i j ki kj = let update i j deg_a deg_b ki kj =
let x = f ki kj in let x = f deg_a deg_b ki kj in
if abs_float x > Constants.epsilon then if abs_float x > Constants.epsilon then
result.(i - shift) <- (j, x) :: result.(i - shift) ; result.(i - shift) <- (j, x) :: result.(i - shift) ;
in in
@ -131,7 +140,7 @@ let create_matrix_arbitrary f det_space =
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 (index_start.(i) + i') update (index_start.(i) + i')
(index_start.(j) + j' + 1) ki kj); (index_start.(j) + j' + 1) 2 0 ki kj);
raise Exit) raise Exit)
) j_dets ) j_dets
with Exit -> () with Exit -> ()
@ -153,7 +162,7 @@ let create_matrix_arbitrary f det_space =
Determinant.of_spindeterminants j_alfa j_beta Determinant.of_spindeterminants j_alfa j_beta
in (update in (update
(index_start.(i) + i') (index_start.(j) + j' + 1) (index_start.(i) + i') (index_start.(j) + j' + 1)
ki kj; 1 (Determinant.degree_beta ki kj) ki kj;
aux r_singles (j'+1);) aux r_singles (j'+1);)
| 1 -> if (j' < Array.length j_dets) then aux singles (j'+1) | 1 -> if (j' < Array.length j_dets) then aux singles (j'+1)
| _ -> assert false | _ -> assert false
@ -177,7 +186,7 @@ let create_matrix_arbitrary f det_space =
Determinant.of_spindeterminants j_alfa j_beta Determinant.of_spindeterminants j_alfa j_beta
in (update in (update
(index_start.(i) + i') (index_start.(j) + j' + 1) (index_start.(i) + i') (index_start.(j) + j' + 1)
ki kj; 0 (Determinant.degree_beta ki kj) ki kj;
aux r_doubles (j'+1);) aux r_doubles (j'+1);)
| 1 -> if (j' < Array.length j_dets) then aux doubles (j'+1) | 1 -> if (j' < Array.length j_dets) then aux doubles (j'+1)
| _ -> assert false | _ -> assert false
@ -268,8 +277,8 @@ let create_matrix_spin f det_space =
in in
(** Update function when ki and kj are connected *) (** Update function when ki and kj are connected *)
let update i j ki kj = let update i j deg_a deg_b ki kj =
let x = f ki kj in let x = f deg_a deg_b ki kj in
if abs_float x > Constants.epsilon then if abs_float x > Constants.epsilon then
result.(i) <- (j, x) :: result.(i) ; result.(i) <- (j, x) :: result.(i) ;
in in
@ -285,7 +294,7 @@ let create_matrix_spin f det_space =
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) 2 0 ki kj;
incr i'; incr i';
) b; ) b;
| 1 -> | 1 ->
@ -295,7 +304,7 @@ let create_matrix_spin f det_space =
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) 1 (Determinant.degree_beta ki kj) ki kj
) singles; ) singles;
incr i'; incr i';
) b; ) b;
@ -306,7 +315,7 @@ let create_matrix_spin f det_space =
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) 0 (Determinant.degree_beta ki kj) ki kj
) doubles; ) doubles;
incr i'; incr i';
) b; ) 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 (* Create a matrix using the fact that the determinant space is made of
the outer product of spindeterminants. *) 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 n_alfa = Array.length a in
let h i_alfa j_alfa = 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 -> | 2 ->
let ai, aj = a.(i_alfa), a.(j_alfa) in let ai, aj = a.(i_alfa), a.(j_alfa) in
(fun i_beta j_beta -> (fun i_beta j_beta ->
if i_beta <> j_beta then 0. else if i_beta <> j_beta then 0. else
let ki = Determinant.of_spindeterminants ai b.(i_beta) in let deg_b = 0 in
let kj = Determinant.of_spindeterminants aj b.(j_beta) in let ki = Determinant.of_spindeterminants ai b.(i_beta) in
f ki kj let kj = Determinant.of_spindeterminants aj b.(j_beta) in
f deg_a deg_b ki kj
) )
| 1 -> | 1 ->
let ai, aj = a.(i_alfa), a.(j_alfa) in let ai, aj = a.(i_alfa), a.(j_alfa) in
(fun i_beta j_beta -> (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 -> | 0 | 1 ->
let ki = Determinant.of_spindeterminants ai b.(i_beta) in let ki = Determinant.of_spindeterminants ai b.(i_beta) in
let kj = Determinant.of_spindeterminants aj b.(j_beta) in let kj = Determinant.of_spindeterminants aj b.(j_beta) in
f ki kj f deg_a deg_b ki kj
| _ -> 0. | _ -> 0.
) )
| 0 -> | 0 ->
let ai, aj = a.(i_alfa), a.(j_alfa) in let ai, aj = a.(i_alfa), a.(j_alfa) in
(fun i_beta j_beta -> (fun i_beta j_beta ->
let ki = Determinant.of_spindeterminants ai b.(i_beta) in let deg_b = Spindeterminant.degree b.(i_beta) b.(j_beta) in
let kj = Determinant.of_spindeterminants aj b.(j_beta) in match deg_b with
f ki kj | 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.) | _ -> (fun _ _ -> 0.)
in in
@ -398,26 +406,24 @@ let create_matrix_spin_computed f det_space =
if i <> !i_prev then if i <> !i_prev then
begin begin
i_prev := i; 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_a = (i-1)/n_beta in
let i_alfa = i_a + 1 in let i_alfa = i_a + 1 in
let h1 = let h1 =
h (i_alfa-1) h (i_alfa-1)
in in
let bi = let i_beta = i - i_a*n_beta in
let i_beta = i - i_a*n_beta in let bi = (i_beta-1) in
i_beta-1
in
let j_alfa_prev = ref (-10) in
let h123_prev = ref (fun _ -> 0.) in let h123_prev = ref (fun _ -> 0.) in
let f j = let j_a = ref (-n_alfa) in
if j > !j0 + n_beta || j < !j0 then let j_alfa_prev = ref (-10) in
begin let j0 = ref (!j_a * n_beta) in
let x = (j-1)/n_beta in result := fun j ->
j_a := x; if j > !j0 + n_beta
j0 := x * n_beta || j < !j0
end; then begin
j_a := (j-1)/n_beta;
j0 := !j_a * n_beta
end;
let j_alfa = !j_a + 1 in let j_alfa = !j_a + 1 in
let h123 = let h123 =
if j_alfa <> !j_alfa_prev then if j_alfa <> !j_alfa_prev then
@ -429,8 +435,6 @@ let create_matrix_spin_computed f det_space =
in in
let j_beta = j - !j0 in let j_beta = j - !j0 in
h123 (j_beta-1) h123 (j_beta-1)
in
result := f ;
end; end;
!result !result
in in
@ -449,7 +453,7 @@ let make ?(n_states=1) ?(algo=`Direct) det_space =
Ds.determinant_stream det_space Ds.determinant_stream det_space
|> Stream.next |> Stream.next
in in
h_ij mo_basis d0 d0 h_ij_non_zero mo_basis 0 0 d0 d0
in in
let m_H = let m_H =
@ -468,11 +472,11 @@ let make ?(n_states=1) ?(algo=`Direct) det_space =
else else
create_matrix_spin create_matrix_spin
in in
f (fun ki kj -> f (fun deg_a deg_b ki kj ->
if ki <> kj then if deg_a + deg_b > 0 then
h_ij mo_basis ki kj h_ij_non_zero mo_basis deg_a deg_b ki kj
else else
h_ij mo_basis ki ki -. e_shift h_ij_non_zero mo_basis 0 0 ki ki -. e_shift
) det_space ) det_space
in in
@ -482,7 +486,8 @@ let make ?(n_states=1) ?(algo=`Direct) det_space =
| Ds.Arbitrary _ -> create_matrix_arbitrary | Ds.Arbitrary _ -> create_matrix_arbitrary
| Ds.Spin _ -> create_matrix_spin | Ds.Spin _ -> create_matrix_spin
in 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 in
let eigensystem = lazy ( let eigensystem = lazy (