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

Introduced reference PT2 energy

Browse Source
This commit is contained in:
Anthony Scemama 2019-03-22 18:16:41 +01:00
parent e09a9011e4
commit dc50c99b43
4 changed files with 210 additions and 59 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

124
CI/CI.ml
View File

@ -411,13 +411,13 @@ let make ?(n_states=1) det_space =
let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
list_holes1 list_particles1 ?(unique=true)
list_holes2 list_particles2
list_holes2 list_particles2 is_internal
i_o1_alfa alfa_o2_i w_alfa psi0 =
let list_holes1 = Array.of_list list_holes1
and list_holes2 = Array.of_list list_holes2
and list_particles2 = Array.of_list list_particles1
and list_particles1 = Array.of_list list_particles2
and list_particles1 = Array.of_list list_particles1
and list_particles2 = Array.of_list list_particles2
in
let psi0 =
@ -429,32 +429,6 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
(Stream.next stream), (Mat.copy_row psi0 (i+1)) )
in
let is_internal =
let m l =
List.fold_left (fun accu i ->
let j = i-1 in Z.(logor accu (shift_left one j))
) Z.zero l
in
let mo_class = DeterminantSpace.mo_class det_space in
let active_mask = m (MOClass.active_mos mo_class) in
let occ_mask = m (MOClass.core_mos mo_class) in
let inactive_mask = m (MOClass.inactive_mos mo_class) in
let occ_mask = Z.logor occ_mask inactive_mask in
let neg_active_mask = Z.lognot active_mask in
fun a ->
let alfa =
Determinant.alfa a
|> Spindeterminant.bitstring
in
if Z.logand neg_active_mask alfa <> occ_mask then
false
else
let beta =
Determinant.beta a
|> Spindeterminant.bitstring
in
Z.logand neg_active_mask beta = occ_mask
in
let symmetric = i_o1_alfa == alfa_o2_i in
@ -606,6 +580,35 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
|> Util.stream_to_list
let is_internal det_space =
let m l =
List.fold_left (fun accu i ->
let j = i-1 in Z.(logor accu (shift_left one j))
) Z.zero l
in
let mo_class = DeterminantSpace.mo_class det_space in
let active_mask = m (MOClass.active_mos mo_class) in
let occ_mask = m (MOClass.core_mos mo_class) in
let inactive_mask = m (MOClass.inactive_mos mo_class) in
let occ_mask = Z.logor occ_mask inactive_mask in
let neg_active_mask = Z.lognot active_mask in
fun a ->
let alfa =
Determinant.alfa a
|> Spindeterminant.bitstring
in
if Z.logand neg_active_mask alfa <> occ_mask then
false
else
let beta =
Determinant.beta a
|> Spindeterminant.bitstring
in
Z.logand neg_active_mask beta = occ_mask
let pt2_en ci =
let mo_basis = Ds.mo_basis ci.det_space in
@ -652,7 +655,7 @@ let pt2_en ci =
+. (fock_diag.(p) -. two_e h p h p s s)
|> ignore;
h_ij mo_basis alfa alfa
| _ -> assert false
| _ -> e0.{1} -. 1.0
in
let e0 = e0.{1} in
fun alfa ->
@ -667,12 +670,13 @@ let pt2_en ci =
in
second_order_sum ci list_holes list_particles list_holes list_particles
i_o1_alfa i_o1_alfa w_alfa psi0
(is_internal ci.det_space) i_o1_alfa i_o1_alfa w_alfa psi0
|> List.fold_left (+.) 0.
let pt2_mp ci =
let mo_basis = Ds.mo_basis ci.det_space in
@ -699,13 +703,14 @@ let pt2_mp ci =
let psi0, _ = Parallel.broadcast ci.eigensystem in
second_order_sum ci list_holes list_particles list_holes list_particles
i_o1_alfa i_o1_alfa w_alfa psi0
(is_internal ci.det_space) i_o1_alfa i_o1_alfa w_alfa psi0
|> List.fold_left (+.) 0.
let variance ci =
let mo_basis = Ds.mo_basis ci.det_space in
let psi0, _ = Parallel.broadcast ci.eigensystem in
let i_o1_alfa = h_ij mo_basis in
@ -718,9 +723,60 @@ let variance ci =
[ MOClass.active_mos mo_class ; MOClass.virtual_mos mo_class ]
in
let psi0, _ = Parallel.broadcast ci.eigensystem in
second_order_sum ci list_holes list_particles list_holes list_particles
i_o1_alfa i_o1_alfa w_alfa psi0
(is_internal ci.det_space) i_o1_alfa i_o1_alfa w_alfa psi0
|> List.fold_left (+.) 0.
let pt2_en_reference ci =
let mo_basis = Ds.mo_basis ci.det_space in
let psi0, e0 = Parallel.broadcast ci.eigensystem in
let aux_basis = mo_basis in
let ds =
DeterminantSpace.fci_of_mo_basis ~frozen_core:false aux_basis
in
let out_dets =
ds
|> DeterminantSpace.determinants_array
|> Array.to_list
|> List.filter (fun i -> not (is_internal ci.det_space i))
|> Array.of_list
in
let in_dets =
DeterminantSpace.determinants_array ci.det_space
in
let m_H_aux =
let h_aa =
Array.map (fun ki -> h_ij aux_basis ki ki) out_dets
in
Array.map (fun ki ->
Array.mapi (fun j kj ->
(h_ij aux_basis ki kj) /. (e0.{1} -. h_aa.(j))
) out_dets
) in_dets
|> Mat.of_array
in
let m_F_aux =
Array.map (fun ki ->
Array.map (fun kj ->
h_ij aux_basis ki kj
) out_dets
) in_dets
|> Mat.of_array
in
let m_HF =
gemm m_H_aux m_F_aux ~transb:`T
in
(gemm ~transa:`T psi0 @@ gemm m_HF psi0).{1,1}

130
CI/F12CI.ml
View File

@ -5,6 +5,7 @@ open Lacaml.D
type t =
{
gamma : float;
mo_basis : MOBasis.t ;
aux_basis : MOBasis.t ;
det_space : DeterminantSpace.t ;
@ -20,14 +21,14 @@ let eigensystem t = Lazy.force t.eigensystem
let f12_integrals mo_basis =
let two_e_ints = MOBasis.two_e_ints mo_basis in
let two_e_ints = MOBasis.f12_ints mo_basis in
( (fun i j _ -> 0.),
(fun i j k l s s' ->
if s' = Spin.other s then
ERI.get_phys two_e_ints i j k l
F12.get_phys two_e_ints i j k l
else
(ERI.get_phys two_e_ints i j k l) -.
(ERI.get_phys two_e_ints i j l k)
(F12.get_phys two_e_ints i j k l) -.
(F12.get_phys two_e_ints i j l k)
) )
@ -51,13 +52,36 @@ let f_ij mo_basis ki kj =
|> List.hd
let dressing_vector f12_amplitudes ci =
let is_internal det_space =
let m l =
List.fold_left (fun accu i ->
let j = i-1 in Z.(logor accu (shift_left one j))
) Z.zero l
in
let mo_class = DeterminantSpace.mo_class det_space in
let aux_mask = m (MOClass.auxiliary_mos mo_class) in
fun a ->
let alfa =
Determinant.alfa a
|> Spindeterminant.bitstring
in
if Z.logand aux_mask alfa <> Z.zero then
false
else
let beta =
Determinant.beta a
|> Spindeterminant.bitstring
in
Z.logand aux_mask beta = Z.zero
let dressing_vector aux_basis f12_amplitudes ci =
debug "Computing dressing vector";
let mo_basis = DeterminantSpace.mo_basis ci.CI.det_space in
let i_o1_alfa = h_ij mo_basis in
let alfa_o2_i = f_ij mo_basis in
(*
let i_o1_alfa = h_ij aux_basis in
let alfa_o2_i = f_ij aux_basis in
let w_alfa _ _ = 1. in
@ -65,19 +89,59 @@ debug "Computing dressing vector";
let list_holes = List.concat
[ MOClass.inactive_mos mo_class ; MOClass.active_mos mo_class ]
and list_particles1 = MOClass.auxiliary_mos mo_class
and list_particles2 = List.concat
and list_particles2 = MOClass.auxiliary_mos mo_class
and list_particles1 = List.concat
[ MOClass.active_mos mo_class ; MOClass.virtual_mos mo_class ; MOClass.auxiliary_mos mo_class ]
in
(*
Util.debug_matrix "f12 amplitudes" f12_amplitudes;
*)
(* Single state here *)
let result =
CI.second_order_sum ci list_holes list_particles1 list_holes list_particles2
i_o1_alfa alfa_o2_i w_alfa f12_amplitudes ~unique:false
(is_internal ci.det_space) i_o1_alfa alfa_o2_i w_alfa f12_amplitudes ~unique:false
|> Vec.of_list
in
Matrix.sparse_of_vector_array [| Vector.sparse_of_vec result |]
*)
let out_dets =
DeterminantSpace.fci_of_mo_basis ~frozen_core:false aux_basis
|> DeterminantSpace.determinants_array
|> Array.to_list
|> List.filter (fun i -> not (is_internal ci.CI.det_space i))
|> Array.of_list
in
let in_dets =
DeterminantSpace.determinants_array ci.CI.det_space
in
let m_H_aux =
Array.map (fun ki ->
Array.map (fun kj ->
h_ij aux_basis ki kj
) out_dets
) in_dets
|> Mat.of_array
in
let m_F_aux =
Array.map (fun ki ->
Array.map (fun kj ->
f_ij aux_basis ki kj
) out_dets
) in_dets
|> Mat.of_array
in
Printf.printf "%d %d\n" (Mat.dim1 m_H_aux) (Mat.dim2 m_H_aux);
let m_HF =
gemm m_H_aux m_F_aux ~transb:`T
in
gemm m_HF f12_amplitudes
|> Matrix.sparse_of_mat
@ -85,6 +149,8 @@ Util.debug_matrix "f12 amplitudes" f12_amplitudes;
let make ~simulation ?(threshold=1.e-12) ?(frozen_core=true) ~mo_basis ~aux_basis_filename () =
let gamma = 1. in
let mo_num = MOBasis.size mo_basis in
(* Add auxiliary basis set *)
@ -118,6 +184,7 @@ let make ~simulation ?(threshold=1.e-12) ?(frozen_core=true) ~mo_basis ~aux_basi
Parallel.broadcast (lazy x)
in
let f12_amplitudes =
(* While in a sequential region, initiate the parallel
4-idx transformation to avoid nested parallel jobs
@ -127,9 +194,9 @@ debug "Four-idx transform of f12 intergals";
let f = fun ki kj ->
if ki <> kj then
f_ij aux_basis ki kj
gamma *. (f_ij aux_basis ki kj)
else
f_ij aux_basis ki kj +. 1.
1. +. gamma *. (f_ij aux_basis ki kj)
in
debug "Computing F matrix";
let m_F =
@ -137,7 +204,10 @@ debug "Computing F matrix";
|> Lazy.force
in
fun ci_coef ->
(*
Util.debug_matrix "F" (Matrix.to_mat m_F);
debug "Solving linear system";
*)
Matrix.ax_eq_b m_F (Matrix.dense_of_mat ci_coef)
|> Matrix.to_mat
in
@ -154,32 +224,43 @@ debug "Solving linear system";
let m_H =
Lazy.force ci.CI.m_H
in
(*
Util.debug_matrix "H" (Matrix.to_mat m_H);
*)
let rec iteration ?(state=1) psi =
(*
debug "Iteration";
Util.debug_matrix "T" (f12_amplitudes psi);
*)
let delta =
dressing_vector (f12_amplitudes psi) ci
dressing_vector aux_basis (f12_amplitudes psi) ci
in
Format.printf "Dressing vector : %a\n@." Matrix.pp_matrix delta;
(*------
(*TODO : enlever le double comptage de la symmetrisation*)
*)
let m_H_dressed = Matrix.to_mat m_H in
(Matrix.dim1 delta) (Matrix.dim2 delta) (Mat.dim1 psi) (Mat.dim2 psi) ;
let f = 1.0 /. psi.{1,1} in
Util.list_range 1 (Mat.dim1 psi)
|> List.iter (fun i -> m_H_dressed.{i,1} <- m_H_dressed.{i,1} +. (Matrix.get delta i 1) /. (psi.{1,1}));
|> List.iter (fun i ->
let delta = (Matrix.get delta i 1) *. f in
m_H_dressed.{i,1} <- m_H_dressed.{i,1} +. delta;
if i <> 1 then
begin
m_H_dressed.{1,i} <- m_H_dressed.{1,i} +. delta;
m_H_dressed.{1,1} <- m_H_dressed.{1,1} -. delta *. psi.{i,1} *. f;
end
);
let eigenvectors, eigenvalues =
Util.diagonalize_symm m_H_dressed
in
let m =
gemm ~transa:`T psi eigenvectors
|> Mat.abs
in
let conv =
Mat.sum m -. (Vec.sum (Mat.copy_diag m))
1.0 -. abs_float ( dot
(Mat.to_col_vecs psi).(0)
(Mat.to_col_vecs eigenvectors).(0) )
in
Printf.printf "Convergence : %f %f\n" conv eigenvalues.{1};
Printf.printf "Convergence : %f %f\n" conv (eigenvalues.{1} +. e_shift);
if conv > threshold then
iteration eigenvectors
else
@ -193,6 +274,7 @@ Format.printf "Dressing vector : %a\n@." Matrix.pp_matrix delta;
------- *)
(*
(*TODO : enlever le double comptage de la symmetrisation*)
let n_states = ci.CI.n_states in
let diagonal =
Vec.init (Matrix.dim1 m_H) (fun i -> Matrix.get m_H i i +. (if i=1 then Matrix.get delta 1 1 else 0.) )
@ -226,6 +308,6 @@ Format.printf "Dressing vector : %a\n@." Matrix.pp_matrix delta;
)
in
{ mo_basis ; aux_basis ; det_space ; ci ; eigensystem }
{ mo_basis ; aux_basis ; det_space ; ci ; eigensystem ; gamma }

5
Utils/FourIdxStorage.ml
View File

@ -303,6 +303,11 @@ let four_index_transform coef source =
incr jk;
get_chem_all_i source ~j ~k ~l
|> Array.iteri (fun i x -> o.{i+1,!jk} <- x)
(*
lacpy ~bc:!jk ~b:o
(Mat.of_col_vecs [| Vec.of_array (get_chem_all_i source ~j ~k ~l) |] )
|> ignore
*)
) range_ao
) range_ao;
(* o_i_jk *)

10
run_cas.ml
View File

@ -72,14 +72,22 @@ let () =
let ci = CI.make space in
Format.fprintf ppf "CAS-CI energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s);
(*
let pt2 = CI.pt2_mp ci in
Format.fprintf ppf "CAS-MP2 energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s +. pt2);
*)
let pt2 = CI.pt2_en ci in
Format.fprintf ppf "CAS-EN2 energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s +. pt2);
Format.fprintf ppf "CAS-EN2 energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s +. pt2);
let pt2 = CI.pt2_en_reference ci in
Format.fprintf ppf "CAS-EN2 energy (reference) : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s +. pt2);
(*
let variance = CI.variance ci in
Format.fprintf ppf "CAS variance : %20.16f@." variance
*)
(*
let s2 = Util.xt_o_x ~o:(CI.s2_matrix ci) ~x:(CI.eigenvectors ci) in
Util.list_range 1 (DeterminantSpace.size space)