From 0f2285b5571a5a37cb7b30d9b14576ff51b9beda Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 4 Mar 2019 22:56:03 +0100 Subject: [PATCH] Improved arbitrary space --- CI/CI.ml | 175 ++++++++++++++++-- ...terminant_space.ml => DeterminantSpace.ml} | 32 +++- ...rminant_space.mli => DeterminantSpace.mli} | 0 ...inant_space.ml => SpindeterminantSpace.ml} | 0 ...ant_space.mli => SpindeterminantSpace.mli} | 0 Makefile.include | 2 +- run_fci.ml | 4 +- 7 files changed, 185 insertions(+), 28 deletions(-) rename CI/{Determinant_space.ml => DeterminantSpace.ml} (86%) rename CI/{Determinant_space.mli => DeterminantSpace.mli} (100%) rename CI/{Spindeterminant_space.ml => SpindeterminantSpace.ml} (100%) rename CI/{Spindeterminant_space.mli => SpindeterminantSpace.mli} (100%) diff --git a/CI/CI.ml b/CI/CI.ml index 5380175..189c4d0 100644 --- a/CI/CI.ml +++ b/CI/CI.ml @@ -1,6 +1,6 @@ open Lacaml.D -module Ds = Determinant_space +module Ds = DeterminantSpace type t = { @@ -53,24 +53,165 @@ let h_ij mo_basis ki kj = let create_matrix_arbitrary f det_space = lazy ( - let det = + let ndet = Ds.size det_space in + let data = match Ds.determinants det_space with - | Ds.Arbitrary _ -> Ds.determinants_array det_space + | Ds.Arbitrary a -> a | _ -> assert false in + let det_alfa = data.Ds.det_alfa + and det_beta = data.Ds.det_beta + and det = data.Ds.det + and index_start = data.Ds.index_start + in - let ndet = Ds.size det_space in - let v = Vec.make0 ndet in + (** Array of (list of singles, list of doubles) in the beta spin *) + let degree_bb = + Array.map (fun det_i -> + let deg = Spindeterminant.degree det_i in + let doubles = + Array.mapi (fun i det_j -> + let d = deg det_j in + if d < 3 then + Some (i,d,det_j) + else + None + ) det_beta + |> Array.to_list + |> Util.list_some + in + let singles = + List.filter (fun (i,d,det_j) -> d < 2) doubles + |> List.map (fun (i,_,det_j) -> (i,det_j)) + in + let doubles = + List.map (fun (i,_,det_j) -> (i,det_j)) doubles + in + (singles, doubles) + ) det_beta + in - Array.init ndet - (fun i -> let ki = det.(i) in - Printf.eprintf "%8d / %8d\r%!" i ndet; - let j = ref 1 in - Ds.determinant_stream det_space - |> Stream.iter (fun kj -> v.{!j} <- f ki kj ; incr j); - Vector.sparse_of_vec v) - |> Matrix.sparse_of_vector_array + let task (i,i_dets) = + let shift = index_start.(i) in + + let result = + Array.init (index_start.(i+1) - shift) + (fun _ -> []) + in + + (** Update function when ki and kj are connected *) + let update i j ki kj = + let x = f ki kj in + if abs_float x > Constants.epsilon then + result.(i - shift) <- (j, x) :: result.(i - shift) ; + in + + let i_alfa = det_alfa.(i) in + let deg_a = Spindeterminant.degree i_alfa in + + Array.iteri (fun j j_dets -> + let j_alfa = det_alfa.(j) in + let degree_a = deg_a j_alfa in + + begin + match degree_a with + | 2 -> + Array.iteri (fun i' i_b -> + try + Array.iteri (fun j' j_b -> + if j_b >= i_b then + ( if j_b = i_b then + ( let i_beta = det_beta.(i_b) in + let ki = Determinant.of_spindeterminants i_alfa i_beta in + let kj = Determinant.of_spindeterminants j_alfa i_beta in + update (index_start.(i) + i') + (index_start.(j) + j' + 1) ki kj); + raise Exit) + ) j_dets + with Exit -> () + ) i_dets + | 1 -> + Array.iteri (fun i' i_b -> + let i_beta = det_beta.(i_b) in + let ki = Determinant.of_spindeterminants i_alfa i_beta in + let singles, _ = degree_bb.(i_b) in + let rec aux singles j' = + match singles with + | [] -> () + | (js, j_beta)::r_singles -> + begin + match compare js j_dets.(j') with + | -1 -> aux r_singles j' + | 0 -> + let kj = + Determinant.of_spindeterminants j_alfa j_beta + in (update + (index_start.(i) + i') (index_start.(j) + j' + 1) + ki kj; + aux r_singles (j'+1);) + | 1 -> if (j' < Array.length j_dets) then aux singles (j'+1) + | _ -> assert false + end + in aux singles 0 + ) i_dets + | 0 -> + Array.iteri (fun i' i_b -> + let i_beta = det_beta.(i_b) in + let ki = Determinant.of_spindeterminants i_alfa i_beta in + let _, doubles = degree_bb.(i_b) in + let rec aux doubles j' = + match doubles with + | [] -> () + | (js, j_beta)::r_doubles -> + begin + match compare js j_dets.(j') with + | -1 -> aux r_doubles j' + | 0 -> + let kj = + Determinant.of_spindeterminants j_alfa j_beta + in (update + (index_start.(i) + i') (index_start.(j) + j' + 1) + ki kj; + aux r_doubles (j'+1);) + | 1 -> if (j' < Array.length j_dets) then aux doubles (j'+1) + | _ -> assert false + end + in aux doubles 0 + ) i_dets + | _ -> (); + end + ) det; + let r = + Array.map (fun l -> + List.rev l + |> Vector.sparse_of_assoc_list ndet + ) result + in (i,r) + in + + + let result = + if Parallel.master then + Array.init ndet (fun _ -> Vector.sparse_of_assoc_list ndet []) + else + Array.init ndet (fun _ -> Vector.sparse_of_assoc_list ndet []) + in + + let n_det_alfa = + Array.length det_alfa + in + Array.mapi (fun i i_dets -> i, i_dets) det + |> Array.to_list + |> Stream.of_list + |> Farm.run ~ordered:false ~f:task + |> Stream.iter (fun (k, r) -> + let shift = index_start.(k) in + let det_k = det.(k) in + Array.iteri (fun j r_j -> result.(shift+det_k.(j)) <- r_j) r; + Printf.eprintf "%8d / %8d\r%!" (k+1) n_det_alfa; + ) ; + Matrix.sparse_of_vector_array result ) @@ -117,13 +258,17 @@ let create_matrix_spin f det_space = in let task (i,i_alfa) = - let result = Array.init n_beta (fun _ -> []) in + let result = + Array.init n_beta (fun _ -> []) + in + (** Update function when ki and kj are connected *) let update i j ki kj = let x = f ki kj in if abs_float x > Constants.epsilon then result.(i) <- (j, x) :: result.(i) ; in + let j = ref 1 in let deg_a = Spindeterminant.degree i_alfa in List.iter (fun j_alfa -> @@ -198,7 +343,7 @@ let make ?(n_states=1) det_space = let mo_basis = Ds.mo_basis det_space in (* While in a sequential region, initiate the parallel - 4-idx transformation + 4-idx transformation to avoid nested parallel jobs *) ignore @@ MOBasis.two_e_ints mo_basis; diff --git a/CI/Determinant_space.ml b/CI/DeterminantSpace.ml similarity index 86% rename from CI/Determinant_space.ml rename to CI/DeterminantSpace.ml index aa55036..7345e46 100644 --- a/CI/Determinant_space.ml +++ b/CI/DeterminantSpace.ml @@ -32,7 +32,7 @@ type t = } -module Ss = Spindeterminant_space +module Ss = SpindeterminantSpace let n_alfa t = t.n_alfa let n_beta t = t.n_beta @@ -44,8 +44,8 @@ let size t = match t.determinants with | Spin (a,b) -> (Array.length a) * (Array.length b) | Arbitrary a -> - let ndet_a = Array.length a.det_alfa in - a.index_start.(ndet_a - 1) + Array.length a.det.(ndet_a - 1) + let ndet_a = Array.length a.det_alfa in + a.index_start.(ndet_a) let determinant_stream t = @@ -150,15 +150,27 @@ let fci_of_mo_basis ?(frozen_core=true) mo_basis = let det_alfa = Ss.spin_determinants det_a and det_beta = Ss.spin_determinants det_b in - Spin (det_alfa, det_beta) - (* let n_det_beta = Array.length det_beta in - Arbitrary { - det_alfa ; det_beta ; - det = Array.make (Array.length det_alfa) (Array.init (Array.length det_beta) (fun i -> i)); - index_start = Array.mapi (fun i _ -> i*n_det_beta) det_alfa; - } + let n_det_alfa = Array.length det_alfa in + + (* + let ndet = n_det_alfa * n_det_beta in + Format.printf "Number of determinants : %d %d %d\n%!" + n_det_alfa n_det_beta ndet; + Spin (det_alfa, det_beta) *) + + let det = Array.make n_det_alfa + (Array.init n_det_beta (fun i -> i)) + in + let index_start = Array.init (n_det_alfa+1) (fun i -> i*n_det_beta) in + let ndet = (index_start.(n_det_alfa)) in + + Format.printf "Number of determinants : %d %d %d\n%!" + n_det_alfa n_det_beta ndet; + Arbitrary { + det_alfa ; det_beta ; det ; index_start + } in { n_alfa ; n_beta ; mo_class ; mo_basis ; determinants } diff --git a/CI/Determinant_space.mli b/CI/DeterminantSpace.mli similarity index 100% rename from CI/Determinant_space.mli rename to CI/DeterminantSpace.mli diff --git a/CI/Spindeterminant_space.ml b/CI/SpindeterminantSpace.ml similarity index 100% rename from CI/Spindeterminant_space.ml rename to CI/SpindeterminantSpace.ml diff --git a/CI/Spindeterminant_space.mli b/CI/SpindeterminantSpace.mli similarity index 100% rename from CI/Spindeterminant_space.mli rename to CI/SpindeterminantSpace.mli diff --git a/Makefile.include b/Makefile.include index 4adbe57..7063e7c 100644 --- a/Makefile.include +++ b/Makefile.include @@ -49,7 +49,7 @@ doc: QCaml.odocl $(OCAMLBUILD) -ocamlc ocamlcp $*.byte -use-ocamlfind $(PKGS) clean: - $(OCAMLBUILD) -clean # rm -rf _build $(ALL_EXE) *.native *.byte + rm QCaml.odocl && $(OCAMLBUILD) -clean debug: run_integrals.native ./debug.sh diff --git a/run_fci.ml b/run_fci.ml index 3d6a248..b95b7e4 100644 --- a/run_fci.ml +++ b/run_fci.ml @@ -56,13 +56,13 @@ let () = in let space = - Determinant_space.fci_of_mo_basis ~frozen_core:false mos + DeterminantSpace.fci_of_mo_basis ~frozen_core:false mos in let ci = CI.make space in Format.fprintf ppf "FCI energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s); (* let s2 = Util.xt_o_x ~o:(CI.s2_matrix ci) ~x:(CI.eigenvectors ci) in - Util.list_range 1 (Determinant_space.size space) + Util.list_range 1 (DeterminantSpace.size space) |> List.iter (fun i -> Format.printf "@[%f@]@;" s2.{i,i}); *)