mirror of
https://gitlab.com/scemama/QCaml.git
synced 2024-11-14 10:03:39 +01:00
F12 almost OK
This commit is contained in:
parent
ec19195a6f
commit
f88409bfbd
33
CI/CI.ml
33
CI/CI.ml
@ -403,25 +403,34 @@ 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 max_int in
|
||||
let j1 = ref min_int in
|
||||
let j0 = ref (2*ndet) in
|
||||
let j1 = ref (2*ndet) in
|
||||
let j_a = ref 0 in
|
||||
result := fun j ->
|
||||
(* 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
|
||||
if (!j0 < j) && (j <= !j1) then
|
||||
()
|
||||
else
|
||||
begin
|
||||
if j < !j1 + n_beta then
|
||||
if (!j1 < j) && (j <= !j1 + n_beta) then
|
||||
begin
|
||||
incr j_a;
|
||||
j0 := !j1;
|
||||
@ -433,9 +442,9 @@ let create_matrix_spin_computed f det_space =
|
||||
end;
|
||||
j1 := !j0 + n_beta;
|
||||
h123 := h1 !j_a i_b;
|
||||
let j_b = j - !j0 - 1 in
|
||||
!h123 j_b
|
||||
end
|
||||
end;
|
||||
let j_b = j - !j0 - 1 in
|
||||
!h123 j_b
|
||||
end;
|
||||
!result
|
||||
in
|
||||
|
||||
48
CI/F12CI.ml
48
CI/F12CI.ml
@ -25,12 +25,13 @@ let f12_integrals mo_basis =
|
||||
else
|
||||
begin
|
||||
let ijkl = F12.get_phys two_e_ints i j k l
|
||||
and ijlk = F12.get_phys two_e_ints i j l k
|
||||
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.375 *. ijkl -. 0.125 *. ijlk
|
||||
0.5 *. ijkl
|
||||
else
|
||||
0.25 *. (ijkl -. ijlk)
|
||||
end
|
||||
@ -48,22 +49,17 @@ 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
|
||||
CIMatrixElement.make integrals ki kj
|
||||
CIMatrixElementF12.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_internal det_space =
|
||||
let mo_class = DeterminantSpace.mo_class det_space in
|
||||
@ -172,17 +168,6 @@ 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 ".%!";
|
||||
@ -192,6 +177,14 @@ 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
|
||||
@ -271,16 +264,17 @@ Printf.printf "det space\n%!";
|
||||
let rec iteration ~state psi =
|
||||
let column_idx = iamax (Mat.to_col_vecs psi).(state-1) in
|
||||
let delta =
|
||||
(* delta_i = {% $\sum_j c_j H_{ij}$ %} *)
|
||||
dressing_vector ~frozen_core aux_basis psi ci
|
||||
|> Matrix.to_mat
|
||||
in
|
||||
|
||||
let f = 1.0 /. psi.{column_idx,state} in
|
||||
let delta_00 =
|
||||
Util.list_range 2 (Mat.dim1 psi)
|
||||
|> List.fold_left (fun accu i -> accu +.
|
||||
(Matrix.get delta i state) *. psi.{i,state} *. f) 0.
|
||||
(* 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} )
|
||||
in
|
||||
let delta = Matrix.to_mat delta in
|
||||
delta.{column_idx,state} <- delta.{column_idx,state} -. delta_00;
|
||||
|
||||
let diagonal =
|
||||
@ -297,7 +291,7 @@ Printf.printf "det space\n%!";
|
||||
Matrix.mm ~transa:`T m_H c
|
||||
|> Matrix.to_mat
|
||||
in
|
||||
let c11 = Matrix.get c state state in
|
||||
let c11 = Matrix.get c column_idx state in
|
||||
Util.list_range 1 (Mat.dim1 w)
|
||||
|> List.iter (fun i ->
|
||||
let dci =
|
||||
|
||||
@ -16,7 +16,7 @@ let make
|
||||
|
||||
let m = n_states in
|
||||
|
||||
(* Create guess vectors u, with randomly initialized unknown vectors. *)
|
||||
(* Create guess vectors u, with unknown vectors initialized to unity. *)
|
||||
let init_vectors m =
|
||||
let init_vector k =
|
||||
Vector.sparse_of_assoc_list n [ (k,1.0) ]
|
||||
@ -63,7 +63,7 @@ let make
|
||||
|> pick_new
|
||||
in
|
||||
|
||||
(* Apply the operator the m last vectors *)
|
||||
(* Apply the operator to the m last vectors *)
|
||||
let w_new =
|
||||
matrix_prod (
|
||||
u_new_ortho
|
||||
|
||||
@ -10,6 +10,7 @@ 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 =
|
||||
@ -86,6 +87,7 @@ 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.
|
||||
@ -174,6 +176,7 @@ 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 =
|
||||
@ -301,6 +304,7 @@ 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
|
||||
|
||||
|
||||
|
||||
@ -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 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,3 +512,181 @@ let four_index_transform 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
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user