mirror of
https://gitlab.com/scemama/QCaml.git
synced 2024-12-22 20:33:36 +01:00
Working on F12
This commit is contained in:
parent
3874dd1d95
commit
8e7107cfba
55
CI/F12CI.ml
55
CI/F12CI.ml
@ -6,6 +6,7 @@ type t =
|
|||||||
aux_basis : MOBasis.t ;
|
aux_basis : MOBasis.t ;
|
||||||
det_space : DeterminantSpace.t ;
|
det_space : DeterminantSpace.t ;
|
||||||
ci : CI.t ;
|
ci : CI.t ;
|
||||||
|
f12_amplitudes : Mat.t;
|
||||||
}
|
}
|
||||||
|
|
||||||
let ci t = t.ci
|
let ci t = t.ci
|
||||||
@ -13,6 +14,7 @@ let mo_basis t = t.mo_basis
|
|||||||
let det_space t = t.det_space
|
let det_space t = t.det_space
|
||||||
let mo_class t = DeterminantSpace.mo_class @@ det_space t
|
let mo_class t = DeterminantSpace.mo_class @@ det_space t
|
||||||
|
|
||||||
|
|
||||||
let f12_integrals mo_basis =
|
let f12_integrals mo_basis =
|
||||||
let two_e_ints = MOBasis.two_e_ints mo_basis in
|
let two_e_ints = MOBasis.two_e_ints mo_basis in
|
||||||
( (fun i j _ -> 0.),
|
( (fun i j _ -> 0.),
|
||||||
@ -26,6 +28,7 @@ let f12_integrals mo_basis =
|
|||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
let h_ij mo_basis ki kj =
|
let h_ij mo_basis ki kj =
|
||||||
let integrals =
|
let integrals =
|
||||||
List.map (fun f -> f mo_basis)
|
List.map (fun f -> f mo_basis)
|
||||||
@ -54,15 +57,29 @@ let dressing_vector ci =
|
|||||||
let w_alfa _ _ = 1. in
|
let w_alfa _ _ = 1. in
|
||||||
|
|
||||||
let mo_class = CI.mo_class ci in
|
let mo_class = CI.mo_class ci in
|
||||||
let list_holes = List.concat
|
|
||||||
[ MOClass.inactive_mos mo_class ; MOClass.active_mos mo_class ]
|
let internal_term =
|
||||||
and list_particles = MOClass.auxiliary_mos mo_class
|
let list_holes = List.concat
|
||||||
|
[ MOClass.inactive_mos mo_class ; MOClass.active_mos mo_class ]
|
||||||
|
and list_particles = List.concat
|
||||||
|
[ MOClass.active_mos mo_class ; MOClass.virtual_mos mo_class ]
|
||||||
|
in
|
||||||
|
CI.second_order_sum ci list_holes list_particles
|
||||||
|
i_o1_alfa alfa_o2_i w_alfa
|
||||||
|
|> Vec.of_list
|
||||||
in
|
in
|
||||||
|
|
||||||
CI.second_order_sum ci list_holes list_particles
|
let external_term =
|
||||||
i_o1_alfa alfa_o2_i w_alfa
|
let list_holes = List.concat
|
||||||
|> Vec.of_list
|
[ MOClass.inactive_mos mo_class ; MOClass.active_mos mo_class ]
|
||||||
|
and list_particles = MOClass.auxiliary_mos mo_class
|
||||||
|
in
|
||||||
|
CI.second_order_sum ci list_holes list_particles
|
||||||
|
i_o1_alfa alfa_o2_i w_alfa
|
||||||
|
|> Vec.of_list
|
||||||
|
in
|
||||||
|
|
||||||
|
Vec.sub external_term internal_term
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
@ -97,6 +114,30 @@ let make ~simulation ?(frozen_core=true) ~mo_basis ~aux_basis_filename () =
|
|||||||
in
|
in
|
||||||
|
|
||||||
let ci = CI.make det_space in
|
let ci = CI.make det_space in
|
||||||
{ mo_basis ; aux_basis ; det_space ; ci }
|
|
||||||
|
let ci_coef, ci_energy = Parallel.broadcast ci.eigensystem in
|
||||||
|
|
||||||
|
let f12_amplitudes =
|
||||||
|
(* While in a sequential region, initiate the parallel
|
||||||
|
4-idx transformation to avoid nested parallel jobs
|
||||||
|
*)
|
||||||
|
ignore @@ MOBasis.f12_ints mo_basis;
|
||||||
|
|
||||||
|
let f = fun ki kj ->
|
||||||
|
if ki <> kj then
|
||||||
|
f_ij mo_basis ki kj
|
||||||
|
else
|
||||||
|
f_ij mo_basis ki kj +. 1.
|
||||||
|
in
|
||||||
|
let m_F =
|
||||||
|
CI.create_matrix_spin f det_space
|
||||||
|
|> Lazy.force
|
||||||
|
in
|
||||||
|
Matrix.ax_eq_b (Matrix.dense_of_sparse m_F) (Matrix.dense_of_mat ci_coef)
|
||||||
|
|> Matrix.to_mat
|
||||||
|
in
|
||||||
|
|
||||||
|
|
||||||
|
{ mo_basis ; aux_basis ; det_space ; ci ; f12_amplitudes }
|
||||||
|
|
||||||
|
|
||||||
|
@ -307,7 +307,22 @@ let mv ?(sparse=false) ?(trans=`N) ?(threshold=epsilon) a b =
|
|||||||
else
|
else
|
||||||
Vector.dense_of_vec dense_result
|
Vector.dense_of_vec dense_result
|
||||||
|
|
||||||
|
let iterative_ax_eq_b ~trans a b =
|
||||||
|
failwith "Not implemented"
|
||||||
|
|
||||||
|
let rec ax_eq_b ?(trans=`N) a b =
|
||||||
|
match a, b with
|
||||||
|
| (Dense a), (Dense b) ->
|
||||||
|
let x = lacpy a in
|
||||||
|
(getrs ~trans x b; Dense x)
|
||||||
|
| (Dense _), (Sparse _) ->
|
||||||
|
let b = dense_of_sparse b in
|
||||||
|
ax_eq_b ~trans a b
|
||||||
|
| _ -> iterative_ax_eq_b ~trans a b
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
(* ------------ Printers ------------ *)
|
||||||
|
|
||||||
let rec pp_matrix ppf = function
|
let rec pp_matrix ppf = function
|
||||||
| Dense m -> Util.pp_matrix ppf m
|
| Dense m -> Util.pp_matrix ppf m
|
||||||
|
@ -60,6 +60,9 @@ val mm : ?transa:trans3 -> ?transb:trans3 -> ?threshold:float -> t -> t -> t
|
|||||||
val mv : ?sparse:bool -> ?trans:trans3 -> ?threshold:float -> t -> Vector.t -> Vector.t
|
val mv : ?sparse:bool -> ?trans:trans3 -> ?threshold:float -> t -> Vector.t -> Vector.t
|
||||||
(** Matrix Vector product *)
|
(** Matrix Vector product *)
|
||||||
|
|
||||||
|
val ax_eq_b : ?trans:trans3 -> t -> t -> t
|
||||||
|
(** Solves A.X = B or A'.X = B *)
|
||||||
|
|
||||||
|
|
||||||
(** {1 Printers } *)
|
(** {1 Printers } *)
|
||||||
|
|
||||||
|
Loading…
Reference in New Issue
Block a user