diff --git a/CI/CI.ml b/CI/CI.ml index c192bec..2f19a68 100644 --- a/CI/CI.ml +++ b/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 diff --git a/CI/F12CI.ml b/CI/F12CI.ml index cc8497d..f727fd6 100644 --- a/CI/F12CI.ml +++ b/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 = diff --git a/Utils/Davidson.ml b/Utils/Davidson.ml index d72d482..0af4e15 100644 --- a/Utils/Davidson.ml +++ b/Utils/Davidson.ml @@ -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 diff --git a/Utils/FourIdxStorage.ml b/Utils/FourIdxStorage.ml index c77494b..f3ede9a 100644 --- a/Utils/FourIdxStorage.ml +++ b/Utils/FourIdxStorage.ml @@ -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 +