LCPQ/QCaml
10
1
Fork 0
mirror of https://gitlab.com/scemama/QCaml.git synced 2024-11-08 23:23:41 +01:00
Code Releases Activity

Compare commits

base: LCPQ:e9a835d7b8c72d201691c6ca1b82e24f55a69891
Branches Tags
LCPQ:master
LCPQ:dev
LCPQ:smallarray
..
compare: LCPQ:c517a6d6a59a612bbc085f1dd6148e448bcb8159
Branches Tags
LCPQ:dev
LCPQ:smallarray
LCPQ:master

No commits in common. "e9a835d7b8c72d201691c6ca1b82e24f55a69891" and "c517a6d6a59a612bbc085f1dd6148e448bcb8159" have entirely different histories.
e9a835d7b8 ... c517a6d6a5

10 changed files with 143 additions and 579 deletions
Show Stats Expand all files Collapse all files

3
Basis/F12.ml
View File

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

33
CI/CI.ml
View File

@ -403,34 +403,25 @@ let create_matrix_spin_computed f det_space =
let h123 = ref (fun _ -> 0.) in
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
begin
i_prev := i;
let i_a = (i-1)/n_beta in
let h1 = h i_a in
let i_b = i - i_a*n_beta - 1 in
let j0 = ref (2*ndet) in
let j1 = ref (2*ndet) in
let j0 = ref max_int in
let j1 = ref min_int in
let j_a = ref 0 in
result := fun j ->
if (!j0 < j) && (j <= !j1) then
()
(* The following test checks if !j0 < j < !j1 *)
if (!j1 - j) lor (j - !j0) > 0 then
begin
let j_b = j - !j0 - 1 in
!h123 j_b
end
else
begin
if (!j1 < j) && (j <= !j1 + n_beta) then
if j < !j1 + n_beta then
begin
incr j_a;
j0 := !j1;
@ -442,9 +433,9 @@ let create_matrix_spin_computed f det_space =
end;
j1 := !j0 + n_beta;
h123 := h1 !j_a i_b;
end;
let j_b = j - !j0 - 1 in
!h123 j_b
let j_b = j - !j0 - 1 in
!h123 j_b
end
end;
!result
in

147
CI/CIMatrixElementF12.ml
View File

@ -1,147 +0,0 @@
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.

110
CI/F12CI.ml
View File

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

127
Utils/Cholesky.ml
View File

@ -47,10 +47,6 @@ let full_ldl 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
assert (Mat.dim2 m_A = n);
@ -62,14 +58,14 @@ let pivoted_ldl threshold m_A =
let v_D_inv = Vec.make0 n in
let compute_d i =
let l_ik =
Mat.col m_Lt i
let compute_d j =
let l_jk =
Mat.col m_Lt j
in
let l_ik__d_k =
Vec.mul ~n:(i-1) l_ik v_D
let l_jk__d_k =
Vec.mul ~n:(j-1) l_jk v_D
in
m_A.{pi.(i),pi.(i)} -. dot ~n:(i-1) l_ik l_ik__d_k
m_A.{pi.(j),pi.(j)} -. dot ~n:(j-1) l_jk l_jk__d_k
in
let compute_l i =
@ -85,51 +81,39 @@ 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 )
in
let maxloc v i =
let rec aux pos value i =
let pivot i =
let rec aux (pos,value) i =
if i > n then
pos
else if v_D.{i} > value then
aux i v_D.{i} (i+1)
else if m_D.{i} > value then
aux i m_D.{i} (i+1)
else
aux pos value (i+1)
in
aux i v.{i} (i+1)
in
let pivot i =
let j = maxloc v_D i in
let j = aux i m_D.{i} (i+1) in
let p_i, p_j = pi.(i), pi.(j) in
pi.(i) <- p_j;
pi.(j) <- p_i;
pi.(i) <- pj;
pi.(j) <- pi;
in
let () =
try
for i=1 to n do
pivot i;
for j=1 to (i-1) do
m_Lt.{j,i} <- compute_l i j;
done;
let d_i = compute_d i in
if abs_float d_i < threshold then
raise Exit;
v_D.{i} <- d_i;
v_D_inv.{i} <- 1. /. d_i;
done
with Exit -> ()
in
m_Lt, v_D, pi
for i=1 to n do
pivot i;
for j=1 to (i-1) do
m_Lt.{j,i} <- compute_l i j;
done;
let d_i = compute_d i in
v_D.{i} <- d_i;
v_D_inv.{i} <- 1. /. d_i;
done;
m_Lt, v_D
let make_ldl ?(threshold=Constants.epsilon) m_A =
pivoted_ldl m_A
let test_case () =
@ -143,10 +127,7 @@ let test_case () =
dot v v
in
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 test_full () =
let m_A =
Mat.of_array [| [| 4. ; 12. ; -16. |] ;
[| 12. ; 37. ; -43. |] ;
@ -168,61 +149,11 @@ let test_case () =
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 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);
()
Alcotest.(check (float 1.e-15)) "full D" 0. (vector_diff v_D v_D_ref)
in
[
"Small", `Quick, test_full_small;
"Full", `Quick, test_full;
"Pivoted", `Quick, test_pivoted ;
"Truncated", `Quick, test_truncated ;
]

70
Utils/Davidson.ml
View File

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

194
Utils/FourIdxStorage.ml
View File

@ -10,7 +10,6 @@ type array2 = (float, Bigarray.float64_elt, Bigarray.fortran_layout) Bigarray.Ar
type storage_t =
| Dense of array2
| Svd of Vec.t * array2 * array2
| Sparse of (int, float) Hashtbl.t
type t =
@ -87,7 +86,6 @@ let unsafe_get_four_index ~r1 ~r2 t =
else
match t.four_index with
| 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
try Hashtbl.find a key
with Not_found -> 0.
@ -176,7 +174,6 @@ let unsafe_set_four_index ~r1 ~r2 ~value t =
end
| Sparse a -> let key = key_of_indices ~r1 ~r2 in
Hashtbl.replace a key value
| Svd _ -> assert false
let set_four_index ~r1 ~r2 ~value t =
@ -304,7 +301,6 @@ let get_chem_all_ij d ~k ~l =
Bigarray.reshape_2 result d.size d.size
| Sparse a ->
Mat.init_cols d.size d.size (fun i j -> get_chem d i j k l)
| Svd _ -> assert false
@ -423,12 +419,16 @@ t
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 four_index_transform coef source =
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 ao_num_2 = ao_num * ao_num in
let ao_mo_num = ao_num * mo_num in
@ -512,181 +512,3 @@ let four_index_transform_dense_sparse ds coef source =
if Parallel.master then Printf.eprintf "\n%!";
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
| _ -> assert false
end
| ({index=i ; value=x}::r1), [] -> aux ({index=i ; value=alpha *. x}::accu) r1 []
| ({index=i ; value=x}::r1), [] -> aux ({index=i ; value=x}::accu) r1 []
| [] , ({index=j ; value=y}::r2) -> aux ({index=j ; value=y}::accu) [] r2
| [] , [] -> {n ; v=List.rev accu}
in
@ -297,7 +297,6 @@ let test_case () =
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) "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);
in
let test_dot () =

28
run_fci_f12.ml
View File

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

7
run_integrals.ml
View File

@ -27,9 +27,8 @@ let run ~out =
| Some x -> x
in
let f12 = F12.gaussian_geminal 1.0 in
let s =
Simulation.of_filenames ~nuclei:nuclei_file ~f12 basis_file
Simulation.of_filenames ~nuclei:nuclei_file basis_file
in
print_endline @@ Nuclei.to_string @@ Simulation.nuclei s;
@ -46,8 +45,10 @@ let run ~out =
NucInt.to_file ~filename:(out_file^".nuc") eN_ints;
KinInt.to_file ~filename:(out_file^".kin") kin_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;
*)
()