LCPQ/QCaml
10
1
Fork 0
mirror of https://gitlab.com/scemama/QCaml.git synced 2024-12-22 20:33:36 +01:00
Code Releases Activity

Merge branch 'master' of gitlab.com:scemama/QCaml

Browse Source
This commit is contained in:
Anthony Scemama 2019-08-21 23:16:45 +02:00
commit 16a0fbacaa
10 changed files with 581 additions and 138 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
Basis/F12.ml
View File

@ -373,9 +373,6 @@ let of_basis_parallel f12 basis =
Fis.create ~size:n `Dense Fis.create ~size:n `Dense
in in
(*
let lambda_inv = -. 1. /. f12.expo_s in
*)
Farm.run ~ordered:false ~f input_stream Farm.run ~ordered:false ~f input_stream
|> Stream.iter (fun l -> |> Stream.iter (fun l ->
Array.iter (fun (i_c,j_c,k_c,l_c,value) -> Array.iter (fun (i_c,j_c,k_c,l_c,value) ->

29
CI/CI.ml
View File

@ -403,25 +403,34 @@ let create_matrix_spin_computed f det_space =
let h123 = ref (fun _ -> 0.) in let h123 = ref (fun _ -> 0.) in
let g i = let g i =
(*
i : index of the i-th determinant. 1 <= i <= ndet
i_prev : index of the i-th determinant in the previous function call.
1 <= i_prev <= ndet
i_a : index of the i_a-th alpha determinant. 0 <= i_a < n_alfa
i_b : index of the i_b-th beta determinant. 0 <= i_b < n_beta
j0 : index - 1 of the first determinant with the same alfa component
0 <= j0 < n_beta*(n_alfa-1)
j1 : index - 1 of the next determinant with the 1st beta determinant
n_beta <= j1 <= ndet
j_a : index of the j_a-th alpha determinant. 0 <= j_a < n_alfa
j_b : index of the j_b-th beta determinant. 0 <= j_b < n_beta
*)
if i <> !i_prev then if i <> !i_prev then
begin begin
i_prev := i; i_prev := i;
let i_a = (i-1)/n_beta in let i_a = (i-1)/n_beta in
let h1 = h i_a in let h1 = h i_a in
let i_b = i - i_a*n_beta - 1 in let i_b = i - i_a*n_beta - 1 in
let j0 = ref max_int in let j0 = ref (2*ndet) in
let j1 = ref min_int in let j1 = ref (2*ndet) in
let j_a = ref 0 in let j_a = ref 0 in
result := fun j -> result := fun j ->
(* The following test checks if !j0 < j < !j1 *) if (!j0 < j) && (j <= !j1) then
if (!j1 - j) lor (j - !j0) > 0 then ()
begin
let j_b = j - !j0 - 1 in
!h123 j_b
end
else else
begin begin
if j < !j1 + n_beta then if (!j1 < j) && (j <= !j1 + n_beta) then
begin begin
incr j_a; incr j_a;
j0 := !j1; j0 := !j1;
@ -433,9 +442,9 @@ let create_matrix_spin_computed f det_space =
end; end;
j1 := !j0 + n_beta; j1 := !j0 + n_beta;
h123 := h1 !j_a i_b; h123 := h1 !j_a i_b;
end;
let j_b = j - !j0 - 1 in let j_b = j - !j0 - 1 in
!h123 j_b !h123 j_b
end
end; end;
!result !result
in in

147
CI/CIMatrixElementF12.ml Normal file
View File

@ -0,0 +1,147 @@
open Lacaml.D
module De = Determinant
module Ex = Excitation
module Sp = Spindeterminant
type t = float list
let non_zero integrals degree_a degree_b ki kj =
let kia = De.alfa ki and kib = De.beta ki
and kja = De.alfa kj and kjb = De.beta kj
in
let diag_element =
let mo_a = Sp.to_list kia
and mo_b = Sp.to_list kib
in
fun one_e two_e ->
let one =
(List.fold_left (fun accu i -> accu +. one_e i i Spin.Alfa) 0. mo_a)
+.
(List.fold_left (fun accu i -> accu +. one_e i i Spin.Beta) 0. mo_b)
in
let two =
let rec aux_same spin accu = function
| [] -> accu
| i :: rest ->
let new_accu =
List.fold_left (fun accu j -> accu +. two_e i j i j spin spin) accu rest
in
aux_same spin new_accu rest
in
let rec aux_opposite accu other = function
| [] -> accu
| i :: rest ->
let new_accu =
List.fold_left (fun accu j -> accu +. two_e i j i j Spin.Alfa Spin.Beta) accu other
in
aux_opposite new_accu other rest
in
(aux_same Spin.Alfa 0. mo_a) +. (aux_same Spin.Beta 0. mo_b) +.
(aux_opposite 0. mo_a mo_b)
in
one +. two
in
let result =
match degree_a, degree_b with
| 1, 1 -> (* alpha-beta double *)
begin
let ha, pa, phase_a = Ex.single_of_spindet kia kja in
let hb, pb, phase_b = Ex.single_of_spindet kib kjb in
let s1 =
match phase_a, phase_b with
| Phase.Pos, Phase.Pos
| Phase.Neg, Phase.Neg -> fun _ two_e -> two_e ha hb pa pb Spin.Alfa Spin.Beta
| Phase.Neg, Phase.Pos
| Phase.Pos, Phase.Neg -> fun _ two_e -> -. two_e ha hb pa pb Spin.Alfa Spin.Beta
in
let kk = Determinant.double_excitation Spin.Alfa ha pb Spin.Beta hb pa ki in
let kka = De.alfa kk and kkb = De.beta kk in
match Spindeterminant.(degree kia kka, degree kib kkb) with
| (1,1) ->
let s2 =
begin
let ha, pa, phase_a = Ex.single_of_spindet kia kka in
let hb, pb, phase_b = Ex.single_of_spindet kib kkb in
match phase_a, phase_b with
| Phase.Pos, Phase.Pos
| Phase.Neg, Phase.Neg -> fun _ two_e -> two_e ha hb pa pb Spin.Alfa Spin.Beta
| Phase.Neg, Phase.Pos
| Phase.Pos, Phase.Neg -> fun _ two_e -> -. two_e ha hb pa pb Spin.Alfa Spin.Beta
end
in fun x two_e -> 0.75 *. (s1 x two_e) +. 0.25 *. (s2 x two_e)
| _ -> fun x two_e -> 0.
end
| 2, 0 -> (* alpha double *)
begin
let h1, p1, h2, p2, phase = Ex.double_of_spindet kia kja in
match phase with
| Phase.Pos -> fun _ two_e -> two_e h1 h2 p1 p2 Spin.Alfa Spin.Alfa
| Phase.Neg -> fun _ two_e -> -. two_e h1 h2 p1 p2 Spin.Alfa Spin.Alfa
end
| 0, 2 -> (* beta double *)
begin
let h1, p1, h2, p2, phase = Ex.double_of_spindet kib kjb in
match phase with
| Phase.Pos -> fun _ two_e -> two_e h1 h2 p1 p2 Spin.Beta Spin.Beta
| Phase.Neg -> fun _ two_e -> -. two_e h1 h2 p1 p2 Spin.Beta Spin.Beta
end
| 1, 0 (* alpha single *)
| 0, 1 (* beta single *)
-> fun _ _ -> 0.
| 0, 0 -> (* diagonal element *)
diag_element
| _ -> assert false
in
List.map (fun (one_e, two_e) -> result one_e two_e) integrals
let make integrals ki kj =
let degree_a, degree_b =
De.degrees ki kj
in
if degree_a+degree_b > 2 then
List.map (fun _ -> 0.) integrals
else
non_zero integrals degree_a degree_b ki kj
let make_s2 ki kj =
let degree_a = De.degree_alfa ki kj in
let kia = De.alfa ki in
let kja = De.alfa kj in
if degree_a > 1 then 0.
else
let degree_b = De.degree_beta ki kj in
let kib = De.beta ki in
let kjb = De.beta kj in
match degree_a, degree_b with
| 1, 1 -> (* alpha-beta double *)
let ha, pa, phase_a = Ex.single_of_spindet kia kja in
let hb, pb, phase_b = Ex.single_of_spindet kib kjb in
if ha = pb && hb = pa then
begin
match phase_a, phase_b with
| Phase.Pos, Phase.Pos
| Phase.Neg, Phase.Neg -> -1.
| Phase.Neg, Phase.Pos
| Phase.Pos, Phase.Neg -> 1.
end
else 0.
| 0, 0 ->
let ba = Sp.bitstring kia and bb = Sp.bitstring kib in
let tmp = Bitstring.logxor ba bb in
let n_a = Bitstring.logand ba tmp |> Bitstring.popcount in
let n_b = Bitstring.logand bb tmp |> Bitstring.popcount in
let s_z = 0.5 *. float_of_int (n_a - n_b) in
float_of_int n_a +. s_z *. (s_z -. 1.)
| _ -> 0.

108
CI/F12CI.ml
View File

@ -24,11 +24,16 @@ let f12_integrals mo_basis =
0. 0.
else else
begin begin
let ijkl = F12.get_phys two_e_ints i j k l
in
let ijlk = F12.get_phys two_e_ints i j l k
in
if s' = Spin.other s then if s' = Spin.other s then
0.5 *. (F12.get_phys two_e_ints i j k l) (* Minus sign because we swap spin variables
instead of orbital variables *)
0.375 *. ijkl +. 0.125 *. ijlk
else else
0.25 *. ((F12.get_phys two_e_ints i j k l) -. 0.25 *. (ijkl -. ijlk)
(F12.get_phys two_e_ints i j l k))
end end
) ) ) )
@ -44,14 +49,6 @@ let h_ij mo_basis ki kj =
|> List.hd |> List.hd
let hf_ij mo_basis ki kj =
let integrals =
List.map (fun f -> f mo_basis)
[ CI.h_integrals ; f12_integrals ]
in
CIMatrixElement.make integrals ki kj
let f_ij mo_basis ki kj = let f_ij mo_basis ki kj =
let integrals = let integrals =
List.map (fun f -> f mo_basis) List.map (fun f -> f mo_basis)
@ -60,8 +57,11 @@ let f_ij mo_basis ki kj =
CIMatrixElement.make integrals ki kj CIMatrixElement.make integrals ki kj
|> List.hd |> List.hd
let hf_ij mo_basis ki kj =
[ h_ij mo_basis ki kj ; f_ij mo_basis ki kj ]
let is_internal det_space =
let is_a_double det_space =
let mo_class = DeterminantSpace.mo_class det_space in let mo_class = DeterminantSpace.mo_class det_space in
let mo_num = Array.length @@ MOClass.mo_class_array mo_class in let mo_num = Array.length @@ MOClass.mo_class_array mo_class in
let m l = let m l =
@ -83,8 +83,8 @@ let is_internal det_space =
and b = Bitstring.logand aux_mask beta and b = Bitstring.logand aux_mask beta
in in
match Bitstring.popcount a + Bitstring.popcount b with match Bitstring.popcount a + Bitstring.popcount b with
| 1 | 2 -> false | 2 -> true
| _ -> true | _ -> false
let dressing_vector ~frozen_core aux_basis f12_amplitudes ci = let dressing_vector ~frozen_core aux_basis f12_amplitudes ci =
@ -107,13 +107,12 @@ let dressing_vector ~frozen_core aux_basis f12_amplitudes ci =
|> DeterminantSpace.determinant_stream |> DeterminantSpace.determinant_stream
in in
(* Select only singly and doubly excited determinants (* Select only doubly excited determinants wrt FCI space *)
wrt FCI space *)
Stream.from (fun _ -> Stream.from (fun _ ->
try try
let rec result () = let rec result () =
let ki = Stream.next s in let ki = Stream.next s in
if is_internal ci.CI.det_space ki then if not (is_a_double ci.CI.det_space ki) then
result () result ()
else else
Some ki Some ki
@ -168,17 +167,6 @@ let dressing_vector ~frozen_core aux_basis f12_amplitudes ci =
if result = [] then None else Some result if result = [] then None else Some result
) )
in in
let result =
let x =
try [ Stream.next out_dets_stream ]
with Stream.Failure -> failwith "Auxiliary basis set does not produce any excited determinant"
in
let m_H_aux, m_F_aux = make_h_and_f x in
let m_HF =
gemm m_H_aux m_F_aux ~transb:`T
in
gemm m_HF f12_amplitudes
in
let iteration input = let iteration input =
Printf.printf ".%!"; Printf.printf ".%!";
@ -189,6 +177,14 @@ let dressing_vector ~frozen_core aux_basis f12_amplitudes ci =
gemm m_HF f12_amplitudes gemm m_HF f12_amplitudes
in in
let result =
let x =
try [ Stream.next out_dets_stream ]
with Stream.Failure -> failwith "Auxiliary basis set does not produce any excited determinant"
in
iteration x
in
input_stream input_stream
|> Farm.run ~ordered:false ~f:iteration |> Farm.run ~ordered:false ~f:iteration
|> Stream.iter (fun hf -> |> Stream.iter (fun hf ->
@ -204,7 +200,7 @@ let dressing_vector ~frozen_core aux_basis f12_amplitudes ci =
let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filename () = let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filename ?(state=1) () =
let f12 = Util.of_some @@ Simulation.f12 simulation in let f12 = Util.of_some @@ Simulation.f12 simulation in
let mo_num = MOBasis.size mo_basis in let mo_num = MOBasis.size mo_basis in
@ -230,20 +226,17 @@ Printf.printf "Add aux basis\n%!";
MOBasis.of_mo_basis s mo_basis MOBasis.of_mo_basis s mo_basis
in in
let () = let () =
Printf.printf "F12 ints\n%!";
ignore @@ MOBasis.f12_ints aux_basis ignore @@ MOBasis.f12_ints aux_basis
in in
let () = let () =
Printf.printf "2e ints\n%!";
ignore @@ MOBasis.two_e_ints aux_basis ignore @@ MOBasis.two_e_ints aux_basis
in in
Printf.printf "det space\n%!";
let det_space = let det_space =
DeterminantSpace.fci_f12_of_mo_basis aux_basis ~frozen_core mo_num DeterminantSpace.fci_f12_of_mo_basis aux_basis ~frozen_core mo_num
in in
let ci = CI.make det_space in let ci = CI.make ~n_states:state det_space in
let ci_coef, ci_energy = let ci_coef, ci_energy =
let x = Lazy.force ci.eigensystem in let x = Lazy.force ci.eigensystem in
@ -263,28 +256,31 @@ Printf.printf "det space\n%!";
Lazy.force ci.CI.m_H Lazy.force ci.CI.m_H
in in
let rec iteration ?(state=1) psi =
let rec iteration ~state psi =
let column_idx = iamax (Mat.to_col_vecs psi).(state-1) in
let delta = let delta =
(* delta_i = {% $\sum_j c_j H_{ij}$ %} *)
dressing_vector ~frozen_core aux_basis psi ci dressing_vector ~frozen_core aux_basis psi ci
|> Matrix.to_mat
in in
let f = 1.0 /. psi.{1,1} in Printf.printf "Cmax : %e\n" psi.{column_idx,state};
Printf.printf "Norm : %e\n" (sqrt (gemm ~transa:`T delta delta).{state,state});
let f = 1.0 /. psi.{column_idx,state} in
let delta_00 = let delta_00 =
Util.list_range 2 (Mat.dim1 psi) (* Delta_00 = {% $\sum_{j \ne x} delta_j c_j / c_x$ %} *)
|> List.fold_left (fun accu i -> accu +. f *. ( (gemm ~transa:`T delta psi).{state,state} -.
(Matrix.get delta i 1) *. psi.{i,1} *. f) 0. delta.{column_idx,state} *. psi.{column_idx,state} )
in in
let delta = Matrix.to_mat delta in Printf.printf "Delta_00 : %e %e\n" delta.{column_idx,state} delta_00;
delta.{1,1} <- delta.{1,1} -. delta_00; delta.{column_idx,state} <- delta.{column_idx,state} -. delta_00;
(*------
TODO SINGLE STATE HERE
*)
let n_states = ci.CI.n_states in
let diagonal = let diagonal =
Vec.init (Matrix.dim1 m_H) (fun i -> Vec.init (Matrix.dim1 m_H) (fun i ->
if i = 1 then if i = column_idx then
Matrix.get m_H i i +. delta.{1,1} *. f Matrix.get m_H i i +. delta.{column_idx,state} *. f
else else
Matrix.get m_H i i Matrix.get m_H i i
) )
@ -295,24 +291,26 @@ TODO SINGLE STATE HERE
Matrix.mm ~transa:`T m_H c Matrix.mm ~transa:`T m_H c
|> Matrix.to_mat |> Matrix.to_mat
in in
let c11 = Matrix.get c 1 1 in let c11 = Matrix.get c column_idx state in
Util.list_range 1 (Mat.dim1 w) Util.list_range 1 (Mat.dim1 w)
|> List.iter (fun i -> |> List.iter (fun i ->
let dci = let dci =
delta.{i,1} *. f ; delta.{i,state} *. f ;
in in
w.{i,1} <- w.{i,1} +. dci *. c11; w.{i,state} <- w.{i,state} +. dci *. c11;
if (i <> 1) then if (i <> column_idx) then
w.{1,1} <- w.{1,1} +. dci *. (Matrix.get c i 1); w.{column_idx,state} <- w.{column_idx,state} +. dci *. (Matrix.get c i state);
); );
Matrix.dense_of_mat w Matrix.dense_of_mat w
in in
let eigenvectors, eigenvalues = let eigenvectors, eigenvalues =
Parallel.broadcast (lazy ( Parallel.broadcast (lazy (
Davidson.make ~threshold:1.e-6 ~guess:psi ~n_states diagonal matrix_prod Davidson.make ~threshold:1.e-6 ~guess:psi ~n_states:state diagonal matrix_prod
)) ))
in in
Vec.iter (fun energy -> Printf.printf "%g\t" energy) eigenvalues;
print_newline ();
let conv = let conv =
1.0 -. abs_float ( dot 1.0 -. abs_float ( dot
@ -320,18 +318,18 @@ TODO SINGLE STATE HERE
(Mat.to_col_vecs eigenvectors).(0) ) (Mat.to_col_vecs eigenvectors).(0) )
in in
if Parallel.master then if Parallel.master then
Printf.printf "F12 Convergence : %e %f\n" conv (eigenvalues.{1} +. e_shift Printf.printf "F12 Convergence : %e %f\n" conv (eigenvalues.{state} +. e_shift
+. Simulation.nuclear_repulsion simulation); +. Simulation.nuclear_repulsion simulation);
if conv > threshold then if conv > threshold then
iteration eigenvectors iteration ~state eigenvectors
else else
let eigenvalues = let eigenvalues =
Vec.map (fun x -> x +. e_shift) eigenvalues Vec.map (fun x -> x +. e_shift) eigenvalues
in in
eigenvectors, eigenvalues eigenvectors, eigenvalues
in in
iteration ci_coef iteration ~state ci_coef
) )
in in

112
Utils/Cholesky.ml
View File

@ -52,6 +52,10 @@ let full_ldl m_A =
let pivoted_ldl threshold m_A = let pivoted_ldl threshold m_A =
(** {% $P A P^\dagger = L D L^\dagger$. %}
Input : Matrix $A$
Output : Matrices $L, D, P$.
*)
let n = Mat.dim1 m_A in let n = Mat.dim1 m_A in
assert (Mat.dim2 m_A = n); assert (Mat.dim2 m_A = n);
@ -63,14 +67,14 @@ let pivoted_ldl threshold m_A =
let v_D_inv = Vec.make0 n in let v_D_inv = Vec.make0 n in
let compute_d j = let compute_d i =
let l_jk = let l_ik =
Mat.col m_Lt j Mat.col m_Lt i
in in
let l_jk__d_k = let l_ik__d_k =
Vec.mul ~n:(j-1) l_jk v_D Vec.mul ~n:(i-1) l_ik v_D
in in
m_A.{pi.(j),pi.(j)} -. dot ~n:(j-1) l_jk l_jk__d_k m_A.{pi.(i),pi.(i)} -. dot ~n:(i-1) l_ik l_ik__d_k
in in
let compute_l i = let compute_l i =
@ -86,8 +90,9 @@ let pivoted_ldl threshold m_A =
v_D_inv.{j} *. ( m_A.{pi.(j),pi.(i)} -. dot ~n:(j-1) l_ik__d_k l_jk ) v_D_inv.{j} *. ( m_A.{pi.(j),pi.(i)} -. dot ~n:(j-1) l_ik__d_k l_jk )
in in
let pivot i =
let rec aux (pos,value) i = let maxloc v i =
let rec aux pos value i =
if i > n then if i > n then
pos pos
else if v_D.{i} > value then else if v_D.{i} > value then
@ -95,23 +100,36 @@ let pivoted_ldl threshold m_A =
else else
aux pos value (i+1) aux pos value (i+1)
in in
let j = aux i v_D.{i} (i+1) in aux i v.{i} (i+1)
let p_i, p_j = pi.(i), pi.(j) in
pi.(i) <- pj;
pi.(j) <- pi;
in in
let pivot i =
let j = maxloc v_D i in
let p_i, p_j = pi.(i), pi.(j) in
pi.(i) <- p_j;
pi.(j) <- p_i;
in
let () =
try
for i=1 to n do for i=1 to n do
pivot i; pivot i;
for j=1 to (i-1) do for j=1 to (i-1) do
m_Lt.{j,i} <- compute_l i j; m_Lt.{j,i} <- compute_l i j;
done; done;
let d_i = compute_d i in let d_i = compute_d i in
if abs_float d_i < threshold then
raise Exit;
v_D.{i} <- d_i; v_D.{i} <- d_i;
v_D_inv.{i} <- 1. /. d_i; v_D_inv.{i} <- 1. /. d_i;
done; done
m_Lt, v_D with Exit -> ()
in
m_Lt, v_D, pi
@ -121,6 +139,11 @@ let make_ldl ?(threshold=Constants.epsilon) m_A =
let test_case () = let test_case () =
let matrix_diff m_A m_B = let matrix_diff m_A m_B =
@ -132,7 +155,10 @@ let test_case () =
dot v v dot v v
in in
let test_full () = let m_A = Mat.random 1000 1000 in
let m_A = Mat.add m_A (Mat.transpose_copy m_A) in
let test_full_small () =
let m_A = let m_A =
Mat.of_array [| [| 4. ; 12. ; -16. |] ; Mat.of_array [| [| 4. ; 12. ; -16. |] ;
[| 12. ; 37. ; -43. |] ; [| 12. ; 37. ; -43. |] ;
@ -154,11 +180,61 @@ let test_case () =
in in
let m_Lt, v_D = full_ldl m_A in let m_Lt, v_D = full_ldl m_A in
Alcotest.(check (float 1.e-15)) "full L" 0. (matrix_diff m_Lt m_Lt_ref); Alcotest.(check (float 1.e-15)) "full L" 0. (matrix_diff m_Lt m_Lt_ref);
Alcotest.(check (float 1.e-15)) "full D" 0. (vector_diff v_D v_D_ref) Alcotest.(check (float 1.e-15)) "full D" 0. (vector_diff v_D v_D_ref);
let m_D = Mat.of_diag v_D in
let m_B = gemm ~transa:`T m_Lt @@ gemm m_D m_Lt in
Alcotest.(check (float 1.e-15)) "full L" 0. (matrix_diff m_A m_B);
()
in
let test_full () =
let m_Lt, v_D = full_ldl m_A in
let m_D = Mat.of_diag v_D in
let m_B = gemm ~transa:`T m_Lt @@ gemm m_D m_Lt in
Alcotest.(check (float 1.e-15)) "full D" 0. (matrix_diff m_A m_B);
()
in
let test_pivoted () =
let m_Lt, v_D, pi = pivoted_ldl 0. m_A in
let n = Mat.dim1 m_A in
let m_P = Mat.make0 n n in
for i=1 to n do
m_P.{i,pi.(i)} <- 1.
done;
let m_D = Mat.of_diag v_D in
let m_B =
gemm ~transa:`T m_P @@
gemm ~transa:`T m_Lt @@
gemm m_D @@
gemm m_Lt m_P
in
Alcotest.(check (float 1.e-14)) "pivoted D" 0. (matrix_diff m_A m_B);
()
in
let test_truncated () =
let m_Lt, v_D, pi = pivoted_ldl 0.001 m_A in
let n = Mat.dim1 m_Lt in
let m_P = Mat.make0 n n in
for i=1 to n do
m_P.{i,pi.(i)} <- 1.
done;
let m_D = Mat.of_diag v_D in
let m_B =
gemm ~transa:`T m_P @@
gemm ~transa:`T m_Lt @@
gemm m_D @@
gemm m_Lt m_P
in
Alcotest.(check (float 1.e-3)) "full D" 0. (matrix_diff m_A m_B);
()
in in
[ [
"Small", `Quick, test_full_small;
"Full", `Quick, test_full; "Full", `Quick, test_full;
"Pivoted", `Quick, test_pivoted ;
"Truncated", `Quick, test_truncated ;
] ]

70
Utils/Davidson.ml
View File

@ -14,15 +14,10 @@ let make
(* Size of the matrix to diagonalize *) (* Size of the matrix to diagonalize *)
let n = Vec.dim diagonal in let n = Vec.dim diagonal in
let m = n_states in
let m = (* Number of requested states *) (* Create guess vectors u, with unknown vectors initialized to unity. *)
match guess with let init_vectors m =
| Some vectors -> (Mat.dim2 vectors) * n_states
| None -> n_states
in
(* Create guess vectors u, with randomly initialized unknown vectors. *)
let init_vectors =
let init_vector k = let init_vector k =
Vector.sparse_of_assoc_list n [ (k,1.0) ] Vector.sparse_of_assoc_list n [ (k,1.0) ]
in in
@ -31,6 +26,20 @@ let make
|> Matrix.to_mat |> Matrix.to_mat
in in
let guess =
match guess with
| Some vectors -> vectors
| None -> init_vectors m
in
let guess =
if Mat.dim2 guess = n_states then guess
else
(Mat.to_col_vecs_list guess) @
(Mat.to_col_vecs_list (init_vectors (m-(Mat.dim2 guess))) )
|> Mat.of_col_vecs_list
in
let pick_new u = let pick_new u =
Mat.to_col_vecs_list u Mat.to_col_vecs_list u
|> Util.list_pack m |> Util.list_pack m
@ -38,19 +47,14 @@ let make
|> List.hd |> List.hd
in in
let u_new = let u_new = Mat.to_col_vecs_list guess in
match guess with
| Some vectors -> Mat.to_col_vecs_list vectors
| None -> Mat.to_col_vecs_list init_vectors
in
let rec iteration u u_new w iter = let rec iteration u u_new w iter macro =
(* u is a list of orthonormal vectors, on which the operator has (* u is a list of orthonormal vectors, on which the operator has
been applied : w = op.u been applied : w = op.u
u_new is a list of vector which will increase the size of the u_new is a list of vectors which will increase the size of the
space. space.
*) *)
(* Orthonormalize input vectors u_new *) (* Orthonormalize input vectors u_new *)
let u_new_ortho = let u_new_ortho =
List.concat [u ; u_new] List.concat [u ; u_new]
@ -59,7 +63,7 @@ let make
|> pick_new |> pick_new
in in
(* Apply the operator the m last vectors *) (* Apply the operator to the m last vectors *)
let w_new = let w_new =
matrix_prod ( matrix_prod (
u_new_ortho u_new_ortho
@ -101,17 +105,29 @@ let make
(* Compute the residual as proposed new vectors *) (* Compute the residual as proposed new vectors *)
let u_proposed = let u_proposed =
Mat.init_cols n m (fun i k -> Mat.init_cols n m (fun i k ->
(lambda.{k} *. m_new_U.{i,k} -. m_new_W.{i,k}) /. let delta = diagonal.{i} -. lambda.{k} in
(max (diagonal.{i} -. lambda.{k}) 0.01) ) let delta =
if abs_float delta > 0.01 then delta
else if delta < 0. then -.0.01
else 0.01
in
(lambda.{k} *. m_new_U.{i,k} -. m_new_W.{i,k}) /. delta )
|> Mat.to_col_vecs_list |> Mat.to_col_vecs_list
in in
let residual_norms = List.map nrm2 u_proposed in let residual_norms = List.map nrm2 u_proposed in
let residual_norm = List.fold_left (fun accu i -> max accu i) 0. residual_norms in let residual_norm =
List.fold_left (fun accu i -> accu +. i *. i) 0. residual_norms
|> sqrt
in
if Parallel.master then if Parallel.master then
Printf.printf "%3d %16.10f %16.8e%!\n" iter lambda.{1} residual_norm; begin
Printf.printf "%3d " iter;
Vec.iteri (fun i x -> if (i<=m) then Printf.printf "%16.10f " x) lambda;
Printf.printf "%16.8e%!\n" residual_norm
end;
(* Make new vectors sparse *) (* Make new vectors sparse *)
let u_proposed = let u_proposed =
@ -125,20 +141,20 @@ let make
in in
if residual_norm > threshold then if residual_norm > threshold && macro > 0 then
let u_next, w_next, iter = let u_next, w_next, iter, macro =
if iter = n_iter then if iter = n_iter then
m_new_U |> pick_new, m_new_U |> pick_new,
m_new_W |> pick_new, m_new_W |> pick_new,
0 0, (macro-1)
else else
u_next, w_next, iter u_next, w_next, (iter+1), macro
in in
iteration u_next u_proposed w_next (iter+1) iteration u_next u_proposed w_next iter macro
else else
(m_new_U |> pick_new |> Mat.of_col_vecs_list), lambda (m_new_U |> pick_new |> Mat.of_col_vecs_list), lambda
in in
iteration [] u_new [] 1 iteration [] u_new [] 1 5

194
Utils/FourIdxStorage.ml
View File

@ -10,6 +10,7 @@ type array2 = (float, Bigarray.float64_elt, Bigarray.fortran_layout) Bigarray.Ar
type storage_t = type storage_t =
| Dense of array2 | Dense of array2
| Svd of Vec.t * array2 * array2
| Sparse of (int, float) Hashtbl.t | Sparse of (int, float) Hashtbl.t
type t = type t =
@ -86,6 +87,7 @@ let unsafe_get_four_index ~r1 ~r2 t =
else else
match t.four_index with match t.four_index with
| Dense a -> unsafe_get a (dense_index i k t.size) (sym_index j l) | Dense a -> unsafe_get a (dense_index i k t.size) (sym_index j l)
| Svd _ -> assert false
| Sparse a -> let key = key_of_indices ~r1 ~r2 in | Sparse a -> let key = key_of_indices ~r1 ~r2 in
try Hashtbl.find a key try Hashtbl.find a key
with Not_found -> 0. with Not_found -> 0.
@ -174,6 +176,7 @@ let unsafe_set_four_index ~r1 ~r2 ~value t =
end end
| Sparse a -> let key = key_of_indices ~r1 ~r2 in | Sparse a -> let key = key_of_indices ~r1 ~r2 in
Hashtbl.replace a key value Hashtbl.replace a key value
| Svd _ -> assert false
let set_four_index ~r1 ~r2 ~value t = let set_four_index ~r1 ~r2 ~value t =
@ -301,6 +304,7 @@ let get_chem_all_ij d ~k ~l =
Bigarray.reshape_2 result d.size d.size Bigarray.reshape_2 result d.size d.size
| Sparse a -> | Sparse a ->
Mat.init_cols d.size d.size (fun i j -> get_chem d i j k l) Mat.init_cols d.size d.size (fun i j -> get_chem d i j k l)
| Svd _ -> assert false
@ -419,16 +423,12 @@ t
let four_index_transform coef source = let four_index_transform_dense_sparse ds coef source =
let mo_num = Mat.dim2 coef in
let destination = create ~size:mo_num ds in
let ao_num = Mat.dim1 coef in let ao_num = Mat.dim1 coef in
let mo_num = Mat.dim2 coef in
let destination =
match source.four_index with
| Dense _ -> create ~size:mo_num `Dense
| Sparse _ -> create ~size:mo_num `Sparse
in
let mo_num_2 = mo_num * mo_num in let mo_num_2 = mo_num * mo_num in
let ao_num_2 = ao_num * ao_num in let ao_num_2 = ao_num * ao_num in
let ao_mo_num = ao_num * mo_num in let ao_mo_num = ao_num * mo_num in
@ -512,3 +512,181 @@ let four_index_transform coef source =
if Parallel.master then Printf.eprintf "\n%!"; if Parallel.master then Printf.eprintf "\n%!";
broadcast destination broadcast destination
let svd_of_dense t =
let four_index =
let f =
match t.four_index with
| Dense a -> a
| _ -> invalid_arg "svd_of_dense expects a Dense storage"
in
let size = t.size in
let b = Mat.create (size*size) (size*size) in
for k = 1 to size do
for l = 1 to size do
let ac = sym_index k l in
lacpy f ~n:1 ~ac ~b ~bc:(dense_index k l size)
|> ignore;
lacpy f ~n:1 ~ac ~b ~bc:(dense_index l k size)
|> ignore
done
done;
let d, u, v = gesdd b in
let rank =
let w = copy d in
scal (1. /. d.{1}) w;
let rec aux i =
if i = Vec.dim w then i else
if w.{i} < 1.e-15 then i else
aux (i+1)
in aux 1
in
let d = copy d ~n:rank in
let u = lacpy ~n:rank u in
let v = lacpy ~m:rank v in
Svd (d,u,v)
in
{ t with four_index }
let dense_of_svd t =
let four_index =
let d, u, v =
match t.four_index with
| Svd (d, u, v) -> d, u, v
| _ -> invalid_arg "dense_of_svd expects a Svd storage"
in
let f =
gemm u @@ gemm (Mat.of_diag d) v
in
let size = t.size in
let b = Mat.create (size*size) ((size*(size+1))/2) in
for k = 1 to size do
for l = 1 to k do
let bc = sym_index k l in
lacpy f ~ac:(dense_index l k size) ~n:1 ~bc ~b
|> ignore
done
done;
Dense b
in
{ t with four_index }
let four_index_transform_svd coef source =
let t0 = Unix.gettimeofday () in
let svd = svd_of_dense source in
Printf.printf "Computed SVD in %f seconds\n%!" (Unix.gettimeofday () -. t0);
let u0, d0, v0 =
match svd.four_index with
| Svd (d, u, v) -> u, d, (Mat.transpose_copy v)
| _ -> assert false
in
let t0 = Unix.gettimeofday () in
let ao_num = Mat.dim1 coef in
let mo_num = Mat.dim2 coef in
let svd_num = Vec.dim d0 in
let destination = create ~size:mo_num `Dense in
let mo_num_2 = mo_num * mo_num in
if Parallel.master then Printf.eprintf "4-idx transformation \n%!";
let u_vecs = Mat.to_col_vecs u0
and v_vecs = Mat.to_col_vecs v0
in
Printf.printf "%d %d %d %d\n" (Mat.dim1 u0) (Mat.dim2 u0) (Mat.dim1 v0) (Mat.dim2 v0);
let task a =
let o =
Bigarray.reshape_2 (Bigarray.genarray_of_array1 u_vecs.(a-1)) ao_num ao_num
in
let u =
let p = gemm ~transa:`T coef @@ gemm o coef in
Bigarray.reshape_1 (Bigarray.genarray_of_array2 p) mo_num_2
in
let o =
Bigarray.reshape_2 (Bigarray.genarray_of_array1 v_vecs.(a-1)) ao_num ao_num
in
let v =
let p = gemm ~transa:`T coef @@ gemm o coef in
Bigarray.reshape_1 (Bigarray.genarray_of_array2 p) mo_num_2
in
(u, v)
in
(*
let four_index = Svd (d,u'',v'') in
let t = make_from_four_index
{ source with
size = mo_num;
two_index;
two_index_anti;
three_index;
three_index_anti;
four_index
}
*)
let n = ref 0 in
let u_list = ref []
and v_list = ref []
in
Util.list_range 1 svd_num
|> Stream.of_list
|> Farm.run ~f:task ~ordered:true ~comm:Parallel.Node.comm
|> Stream.iter (fun (u,v) ->
if Parallel.master then (incr n ; Printf.eprintf "\r%d / %d%!" !n mo_num);
u_list := u :: !u_list;
v_list := v :: !v_list);
let u = Mat.of_col_vecs_list (List.rev !u_list)
and v = Mat.of_col_vecs_list (List.rev !v_list)
in
Printf.printf "Computed 4-idx in %f seconds\n%!" (Unix.gettimeofday () -. t0);
let result =
let p = gemm u @@ gemm ~transb:`T (Mat.of_diag d0) v in
Bigarray.reshape (Bigarray.genarray_of_array2 p) [| mo_num ; mo_num ; mo_num ; mo_num |]
in
for a=1 to mo_num do
for b=a to mo_num do
for c=1 to mo_num do
for d=c to mo_num do
let x = result.{a,b,c,d} in
if abs_float x > epsilon then
set_chem destination a b c d x;
done
done
done
done;
Printf.printf "Computed 4-idx in %f seconds\n%!" (Unix.gettimeofday () -. t0);
if Parallel.master then Printf.eprintf "\n%!";
broadcast destination
let four_index_transform coef source =
match source.four_index with
| Sparse _ -> four_index_transform_dense_sparse `Sparse coef source
| Svd _ -> assert false (* four_index_transform_svd coef source *)
| Dense _ -> four_index_transform_dense_sparse `Dense coef source

3
Utils/Vector.ml
View File

@ -171,7 +171,7 @@ let axpy ?(threshold=epsilon) ?(alpha=1.) x y =
in aux new_accu r1 r2 in aux new_accu r1 r2
| _ -> assert false | _ -> assert false
end end
| ({index=i ; value=x}::r1), [] -> aux ({index=i ; value=x}::accu) r1 [] | ({index=i ; value=x}::r1), [] -> aux ({index=i ; value=alpha *. x}::accu) r1 []
| [] , ({index=j ; value=y}::r2) -> aux ({index=j ; value=y}::accu) [] r2 | [] , ({index=j ; value=y}::r2) -> aux ({index=j ; value=y}::accu) [] r2
| [] , [] -> {n ; v=List.rev accu} | [] , [] -> {n ; v=List.rev accu}
in in
@ -297,6 +297,7 @@ let test_case () =
Alcotest.(check bool) "dense dense axpy" true (axpy ~alpha:3. v1 v2 = v6); Alcotest.(check bool) "dense dense axpy" true (axpy ~alpha:3. v1 v2 = v6);
Alcotest.(check bool) "dense sparse axpy" true (sub ~threshold:1.e-12 (axpy ~alpha:3. v1 v2_s) v6_s = zero_s); Alcotest.(check bool) "dense sparse axpy" true (sub ~threshold:1.e-12 (axpy ~alpha:3. v1 v2_s) v6_s = zero_s);
Alcotest.(check bool) "sparse dense axpy" true (sub ~threshold:1.e-12 (axpy ~alpha:3. v1_s v2) v6_s = zero_s); Alcotest.(check bool) "sparse dense axpy" true (sub ~threshold:1.e-12 (axpy ~alpha:3. v1_s v2) v6_s = zero_s);
Alcotest.(check bool) "sparse sparse axpy" true (sub ~threshold:1.e-12 (axpy ~alpha:3. v1_s v2_s) v6_s = zero_s); Alcotest.(check bool) "sparse sparse axpy" true (sub ~threshold:1.e-12 (axpy ~alpha:3. v1_s v2_s) v6_s = zero_s);
in in
let test_dot () = let test_dot () =

28
run_fci_f12.ml
View File

@ -1,3 +1,5 @@
open Lacaml.D
let () = let () =
let open Command_line in let open Command_line in
begin begin
@ -27,6 +29,10 @@ let () =
{ short='e' ; long="expo" ; opt=Optional; { short='e' ; long="expo" ; opt=Optional;
arg=With_arg "<float>"; arg=With_arg "<float>";
doc="Exponent of the Gaussian geminal"; } ; doc="Exponent of the Gaussian geminal"; } ;
{ short='s' ; long="state" ; opt=Optional;
arg=With_arg "<int>";
doc="Requested state. Default is 1."; } ;
] ]
end; end;
@ -49,6 +55,12 @@ let () =
| None -> 1.0 | None -> 1.0
in in
let state =
match Command_line.get "state" with
| Some x -> int_of_string x
| None -> 1
in
let multiplicity = let multiplicity =
match Command_line.get "multiplicity" with match Command_line.get "multiplicity" with
| Some x -> int_of_string x | Some x -> int_of_string x
@ -76,11 +88,21 @@ let () =
in in
let fcif12 = let fcif12 =
F12CI.make ~simulation ~frozen_core:false ~mo_basis ~aux_basis_filename () F12CI.make ~simulation ~frozen_core:false ~mo_basis ~aux_basis_filename ~state ()
in in
let ci = F12CI.ci fcif12 in let ci = F12CI.ci fcif12 in
Format.fprintf ppf "FCI energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion simulation); Format.fprintf ppf "FCI energy : ";
Vec.iteri (fun i x -> if i <= state then
Format.fprintf ppf "%20.16f@; " (x +. Simulation.nuclear_repulsion simulation) )
(CI.eigenvalues ci);
Format.fprintf ppf "@.";
let _, e_cif12 = F12CI.eigensystem fcif12 in let _, e_cif12 = F12CI.eigensystem fcif12 in
Format.fprintf ppf "FCI-F12 energy : %20.16f@." (e_cif12.{1} +. Simulation.nuclear_repulsion simulation); Format.fprintf ppf "FCI-F12 energy : ";
Vec.iteri (fun i x -> if i <= state then
Format.fprintf ppf "%20.16f@; " (x +. Simulation.nuclear_repulsion simulation) )
e_cif12;
Format.fprintf ppf "@."

5
run_integrals.ml
View File

@ -27,8 +27,9 @@ let run ~out =
| Some x -> x | Some x -> x
in in
let f12 = F12.gaussian_geminal 1.0 in
let s = let s =
Simulation.of_filenames ~nuclei:nuclei_file basis_file Simulation.of_filenames ~nuclei:nuclei_file ~f12 basis_file
in in
print_endline @@ Nuclei.to_string @@ Simulation.nuclei s; print_endline @@ Nuclei.to_string @@ Simulation.nuclei s;
@ -45,10 +46,8 @@ let run ~out =
NucInt.to_file ~filename:(out_file^".nuc") eN_ints; NucInt.to_file ~filename:(out_file^".nuc") eN_ints;
KinInt.to_file ~filename:(out_file^".kin") kin_ints; KinInt.to_file ~filename:(out_file^".kin") kin_ints;
ERI.to_file ~filename:(out_file^".eri") ee_ints; ERI.to_file ~filename:(out_file^".eri") ee_ints;
(*
let f12_ints = AOBasis.f12_ints ao_basis in let f12_ints = AOBasis.f12_ints ao_basis in
F12.to_file ~filename:(out_file^".f12") f12_ints; F12.to_file ~filename:(out_file^".f12") f12_ints;
*)
() ()