LCPQ/QCaml
10
1
Fork 0
mirror of https://gitlab.com/scemama/QCaml.git synced 2024-10-05 07:46:07 +02:00
Code Releases Activity

CAS+EN OK

Browse Source
This commit is contained in:
Anthony Scemama 2019-03-18 23:38:01 +01:00
parent f07dc33c73
commit e7780804bf
7 changed files with 294 additions and 176 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

420
CI/CI.ml
View File

@ -94,45 +94,45 @@ let create_matrix_arbitrary f det_space =
in in
let task (i,i_dets) = let task (i,i_dets) =
let shift = index_start.(i) in let shift = index_start.(i) in
let result = let result =
Array.init (index_start.(i+1) - shift) Array.init (index_start.(i+1) - shift)
(fun _ -> []) (fun _ -> [])
in in
(** Update function when ki and kj are connected *) (** Update function when ki and kj are connected *)
let update i j ki kj = let update i j ki kj =
let x = f ki kj in let x = f ki kj in
if abs_float x > Constants.epsilon then if abs_float x > Constants.epsilon then
result.(i - shift) <- (j, x) :: result.(i - shift) ; result.(i - shift) <- (j, x) :: result.(i - shift) ;
in in
let i_alfa = det_alfa.(i) in let i_alfa = det_alfa.(i) in
let deg_a = Spindeterminant.degree i_alfa in let deg_a = Spindeterminant.degree i_alfa in
Array.iteri (fun j j_dets -> Array.iteri (fun j j_dets ->
let j_alfa = det_alfa.(j) in let j_alfa = det_alfa.(j) in
let degree_a = deg_a j_alfa in let degree_a = deg_a j_alfa in
begin begin
match degree_a with match degree_a with
| 2 -> | 2 ->
Array.iteri (fun i' i_b -> Array.iteri (fun i' i_b ->
try try
Array.iteri (fun j' j_b -> Array.iteri (fun j' j_b ->
if j_b >= i_b then if j_b >= i_b then
( if j_b = i_b then ( if j_b = i_b then
( let i_beta = det_beta.(i_b) in ( let i_beta = det_beta.(i_b) in
let ki = Determinant.of_spindeterminants i_alfa i_beta in let ki = Determinant.of_spindeterminants i_alfa i_beta in
let kj = Determinant.of_spindeterminants j_alfa i_beta in let kj = Determinant.of_spindeterminants j_alfa i_beta in
update (index_start.(i) + i') update (index_start.(i) + i')
(index_start.(j) + j' + 1) ki kj); (index_start.(j) + j' + 1) ki kj);
raise Exit) raise Exit)
) j_dets ) j_dets
with Exit -> () with Exit -> ()
) i_dets ) i_dets
| 1 -> | 1 ->
Array.iteri (fun i' i_b -> Array.iteri (fun i' i_b ->
let i_beta = det_beta.(i_b) in let i_beta = det_beta.(i_b) in
let ki = Determinant.of_spindeterminants i_alfa i_beta in let ki = Determinant.of_spindeterminants i_alfa i_beta in
@ -145,18 +145,18 @@ let create_matrix_arbitrary f det_space =
match compare js j_dets.(j') with match compare js j_dets.(j') with
| -1 -> aux r_singles j' | -1 -> aux r_singles j'
| 0 -> | 0 ->
let kj = let kj =
Determinant.of_spindeterminants j_alfa j_beta Determinant.of_spindeterminants j_alfa j_beta
in (update in (update
(index_start.(i) + i') (index_start.(j) + j' + 1) (index_start.(i) + i') (index_start.(j) + j' + 1)
ki kj; ki kj;
aux r_singles (j'+1);) aux r_singles (j'+1);)
| 1 -> if (j' < Array.length j_dets) then aux singles (j'+1) | 1 -> if (j' < Array.length j_dets) then aux singles (j'+1)
| _ -> assert false | _ -> assert false
end end
in aux singles 0 in aux singles 0
) i_dets ) i_dets
| 0 -> | 0 ->
Array.iteri (fun i' i_b -> Array.iteri (fun i' i_b ->
let i_beta = det_beta.(i_b) in let i_beta = det_beta.(i_b) in
let ki = Determinant.of_spindeterminants i_alfa i_beta in let ki = Determinant.of_spindeterminants i_alfa i_beta in
@ -169,27 +169,27 @@ let create_matrix_arbitrary f det_space =
match compare js j_dets.(j') with match compare js j_dets.(j') with
| -1 -> aux r_doubles j' | -1 -> aux r_doubles j'
| 0 -> | 0 ->
let kj = let kj =
Determinant.of_spindeterminants j_alfa j_beta Determinant.of_spindeterminants j_alfa j_beta
in (update in (update
(index_start.(i) + i') (index_start.(j) + j' + 1) (index_start.(i) + i') (index_start.(j) + j' + 1)
ki kj; ki kj;
aux r_doubles (j'+1);) aux r_doubles (j'+1);)
| 1 -> if (j' < Array.length j_dets) then aux doubles (j'+1) | 1 -> if (j' < Array.length j_dets) then aux doubles (j'+1)
| _ -> assert false | _ -> assert false
end end
in aux doubles 0 in aux doubles 0
) i_dets ) i_dets
| _ -> (); | _ -> ();
end end
) det; ) det;
let r = let r =
Array.map (fun l -> Array.map (fun l ->
List.rev l List.rev l
|> Vector.sparse_of_assoc_list ndet |> Vector.sparse_of_assoc_list ndet
) result ) result
in (i,r) in (i,r)
in in
let result = let result =
@ -259,63 +259,63 @@ let create_matrix_spin f det_space =
in in
let task (i,i_alfa) = let task (i,i_alfa) =
let result = let result =
Array.init n_beta (fun _ -> []) Array.init n_beta (fun _ -> [])
in in
(** Update function when ki and kj are connected *) (** Update function when ki and kj are connected *)
let update i j ki kj = let update i j ki kj =
let x = f ki kj in let x = f ki kj in
if abs_float x > Constants.epsilon then if abs_float x > Constants.epsilon then
result.(i) <- (j, x) :: result.(i) ; result.(i) <- (j, x) :: result.(i) ;
in in
let j = ref 1 in let j = ref 1 in
let deg_a = Spindeterminant.degree i_alfa in let deg_a = Spindeterminant.degree i_alfa in
List.iter (fun j_alfa -> List.iter (fun j_alfa ->
let degree_a = deg_a j_alfa in let degree_a = deg_a j_alfa in
begin begin
match degree_a with match degree_a with
| 2 -> | 2 ->
let i' = ref 0 in let i' = ref 0 in
List.iteri (fun ib i_beta -> List.iteri (fun ib i_beta ->
let ki = Determinant.of_spindeterminants i_alfa i_beta in let ki = Determinant.of_spindeterminants i_alfa i_beta in
let kj = Determinant.of_spindeterminants j_alfa i_beta in let kj = Determinant.of_spindeterminants j_alfa i_beta in
update !i' (ib + !j) ki kj; update !i' (ib + !j) ki kj;
incr i'; incr i';
) b; ) b;
| 1 -> | 1 ->
let i' = ref 0 in let i' = ref 0 in
List.iteri (fun ib i_beta -> List.iteri (fun ib i_beta ->
let ki = Determinant.of_spindeterminants i_alfa i_beta in let ki = Determinant.of_spindeterminants i_alfa i_beta in
let singles, _ = degree_bb.(ib) in let singles, _ = degree_bb.(ib) in
List.iter (fun (j', j_beta) -> List.iter (fun (j', j_beta) ->
let kj = Determinant.of_spindeterminants j_alfa j_beta in let kj = Determinant.of_spindeterminants j_alfa j_beta in
update !i' (j' + !j) ki kj update !i' (j' + !j) ki kj
) singles; ) singles;
incr i'; incr i';
) b; ) b;
| 0 -> | 0 ->
let i' = ref 0 in let i' = ref 0 in
List.iteri (fun ib i_beta -> List.iteri (fun ib i_beta ->
let ki = Determinant.of_spindeterminants i_alfa i_beta in let ki = Determinant.of_spindeterminants i_alfa i_beta in
let _singles, doubles = degree_bb.(ib) in let _singles, doubles = degree_bb.(ib) in
List.iter (fun (j', j_beta) -> List.iter (fun (j', j_beta) ->
let kj = Determinant.of_spindeterminants j_alfa j_beta in let kj = Determinant.of_spindeterminants j_alfa j_beta in
update !i' (j' + !j) ki kj update !i' (j' + !j) ki kj
) doubles; ) doubles;
incr i'; incr i';
) b; ) b;
| _ -> (); | _ -> ();
end; end;
j := !j + n_beta j := !j + n_beta
) a; ) a;
let r = let r =
Array.map (fun l -> Array.map (fun l ->
List.rev l List.rev l
|> Vector.sparse_of_assoc_list ndet |> Vector.sparse_of_assoc_list ndet
) result ) result
in (i,r) in (i,r)
in in
let result = let result =
@ -348,7 +348,7 @@ let make ?(n_states=1) det_space =
(* While in a sequential region, initiate the parallel (* While in a sequential region, initiate the parallel
4-idx transformation to avoid nested parallel jobs 4-idx transformation to avoid nested parallel jobs
*) *)
ignore @@ MOBasis.two_e_ints mo_basis; ignore @@ MOBasis.two_e_ints mo_basis;
let f = let f =
@ -385,7 +385,9 @@ let make ?(n_states=1) det_space =
let pt2 { det_space ; m_H ; m_S2 ; eigensystem ; n_states } =
let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
i_o1_alfa alfa_o2_i w_alfa =
let mo_basis = Ds.mo_basis det_space in let mo_basis = Ds.mo_basis det_space in
@ -394,23 +396,22 @@ let pt2 { det_space ; m_H ; m_S2 ; eigensystem ; n_states } =
let cls = MOClass.mo_class_array mo_class in let cls = MOClass.mo_class_array mo_class in
Util.list_range 1 (MOBasis.size mo_basis) Util.list_range 1 (MOBasis.size mo_basis)
|> List.filter (fun i -> match cls.(i) with |> List.filter (fun i -> match cls.(i) with
| MOClass.Deleted _ | MOClass.Deleted _
| MOClass.Core _ -> false | MOClass.Core _ -> false
| _ -> true | _ -> true
) )
in in
(* Only the gournd state is computed here *)
let psi0, e0 = Lazy.force eigensystem in
let psi0 = let psi0 =
let psi0, _ = Lazy.force eigensystem in
let stream = let stream =
Ds.determinant_stream det_space Ds.determinant_stream det_space
in in
Array.init (Ds.size det_space) (fun i -> Array.init (Ds.size det_space) (fun i ->
Stream.next stream, psi0.{i+1,1}) Stream.next stream, psi0.{i+1,1})
in in
let e0 = e0.{1} in
(* (*
let is_internal = let is_internal =
@ -445,11 +446,11 @@ let pt2 { det_space ; m_H ; m_S2 ; eigensystem ; n_states } =
let psi_filtered_idx = let psi_filtered_idx =
let rec aux accu = function let rec aux accu = function
| j when j <= i -> List.rev accu | j when j < i -> List.rev accu
| j -> if Determinant.degree (fst psi0.(i)) (fst psi0.(j)) < 4 then | j -> if Determinant.degree (fst psi0.(i)) (fst psi0.(j)) < 4 then
aux (j::accu) (j-1) aux (j::accu) (j-1)
else else
aux accu (j-1) aux accu (j-1)
in aux [] (Array.length psi0 - 1) in aux [] (Array.length psi0 - 1)
in in
@ -457,18 +458,35 @@ let pt2 { det_space ; m_H ; m_S2 ; eigensystem ; n_states } =
List.map (fun i -> psi0.(i)) psi_filtered_idx List.map (fun i -> psi0.(i)) psi_filtered_idx
in in
let psi_h alfa = let symmetric = i_o1_alfa == alfa_o2_i in
let hij = h_ij mo_basis alfa in
let psi_h_alfa alfa =
List.fold_left (fun accu (det, coef) -> List.fold_left (fun accu (det, coef) ->
accu +. coef *. (hij det)) 0. psi_filtered accu +. coef *. (i_o1_alfa det alfa)) 0. psi_filtered
in
let alfa_h_psi =
if symmetric then
psi_h_alfa
else
fun alfa ->
List.fold_left (fun accu (det, coef) ->
accu +. coef *. (alfa_o2_i alfa det)) 0. psi_filtered
in
let psi_h_alfa_alfa_h_psi alfa =
if symmetric then
let x = psi_h_alfa alfa in x *. x
else
(psi_h_alfa alfa) *. (alfa_h_psi alfa)
in in
let is_internal alfa = let is_internal alfa =
let rec aux = function let rec aux = function
| -1 -> false | -1 -> false
| j -> Determinant.degree (fst psi0.(j)) alfa = 0 || aux (j-1) | j -> Determinant.degree (fst psi0.(j)) alfa = 0 || aux (j-1)
in in
aux (Array.length psi0 - 1) aux (Array.length psi0 - 1)
in in
let already_generated alfa = let already_generated alfa =
@ -476,60 +494,136 @@ let pt2 { det_space ; m_H ; m_S2 ; eigensystem ; n_states } =
true true
else else
let rec aux = function let rec aux = function
| -1 -> false | -1 -> false
| j -> Determinant.degree (fst psi0.(j)) alfa <= 2 || aux (j-1) | j -> Determinant.degree (fst psi0.(j)) alfa <= 2 || aux (j-1)
in in
aux (i-1) aux (i-1)
in in
let det_i = fst psi0.(i) in let det_i = fst psi0.(i) in
let w_alfa = w_alfa det_i in
List.fold_left (fun accu particle -> let same_spin =
accu +. List.fold_left (fun accu spin ->
List.fold_left (fun accu hole -> accu +.
if hole = particle then accu else List.fold_left (fun accu particle ->
accu +. accu +.
List.fold_left (fun accu spin -> List.fold_left (fun accu hole ->
let alfa = if hole = particle then accu else
Determinant.single_excitation spin hole particle det_i let alfa =
in Determinant.single_excitation spin hole particle det_i
if Determinant.is_none alfa then accu else
let single =
if already_generated alfa then 0. else
let h_aa = h_ij mo_basis alfa alfa in
let psi_h_alfa = psi_h alfa in
psi_h_alfa *. psi_h_alfa /. (e0 -. h_aa)
in in
if Determinant.is_none alfa then accu else
let double = let single =
List.fold_left (fun accu particle' -> if already_generated alfa then 0. else
accu +. w_alfa alfa *. psi_h_alfa_alfa_h_psi alfa
List.fold_left (fun accu hole' -> in
let double =
List.fold_left (fun accu particle' ->
if particle' > particle then
accu
else
accu +. accu +.
List.fold_left (fun accu spin' -> List.fold_left (fun accu hole' ->
let alfa = if hole' = particle' || hole' < hole then
Determinant.double_excitation
spin hole particle
spin' hole' particle' det_i
in
if Determinant.is_none alfa ||
already_generated alfa then
accu accu
else else
let h_aa = h_ij mo_basis alfa alfa in let alfa =
let psi_h_alfa = psi_h alfa in Determinant.double_excitation
accu +. psi_h_alfa *. psi_h_alfa /. (e0 -. h_aa) spin hole particle
) 0. [ Spin.Alfa ; Spin.Beta ] spin hole' particle' det_i
) 0. mo_indices in
if Determinant.is_none alfa ||
already_generated alfa then
accu
else
accu +. w_alfa alfa *. psi_h_alfa_alfa_h_psi alfa
) 0. mo_indices
) 0. mo_indices ) 0. mo_indices
in in
accu +. single +. double accu +. single +. double
) 0. [ Spin.Alfa ; Spin.Beta ] ) 0. mo_indices
) 0. mo_indices ) 0. mo_indices
) 0. mo_indices ) 0. [ Spin.Alfa ; Spin.Beta ]
in in
let opposite_spin =
List.fold_left (fun accu particle ->
accu +.
List.fold_left (fun accu hole ->
if hole = particle then accu else
let alfa =
Determinant.single_excitation Spin.Alfa hole particle det_i
in
if Determinant.is_none alfa then accu else
let double =
List.fold_left (fun accu particle' ->
accu +.
List.fold_left (fun accu hole' ->
if hole' = particle' then accu else
let alfa =
Determinant.double_excitation
Spin.Alfa hole particle
Spin.Beta hole' particle' det_i
in
if Determinant.is_none alfa ||
already_generated alfa then
accu
else
accu +. w_alfa alfa *. psi_h_alfa_alfa_h_psi alfa
) 0. mo_indices
) 0. mo_indices
in
accu +. double
) 0. mo_indices
) 0. mo_indices
in
same_spin +. opposite_spin
in
Array.mapi (fun i (_,c_i) -> c_i *. det_contribution i) psi0 Array.mapi (fun i (_,c_i) -> c_i *. det_contribution i) psi0
|> Array.fold_left (+.) 0. |> Array.fold_left (+.) 0.
let pt2_en ci =
let mo_basis = Ds.mo_basis ci.det_space in
let _psi0, e0 = Lazy.force ci.eigensystem in
let i_o1_alfa = h_ij mo_basis in
let w_alfa _ =
let e0 = e0.{1} in
let h_aa alfa = h_ij mo_basis alfa alfa in
fun alfa ->
1. /. (e0 -. h_aa alfa)
in
second_order_sum ci i_o1_alfa i_o1_alfa w_alfa
let pt2_mp ci =
let mo_basis = Ds.mo_basis ci.det_space in
let i_o1_alfa = h_ij mo_basis in
let eps = MOBasis.mo_energies mo_basis in
let w_alfa det_i alfa=
match Excitation.of_det det_i alfa with
| Excitation.Single (_, { hole ; particle ; spin })->
1./.(eps.{hole} -. eps.{particle})
| Excitation.Double (_, { hole=h ; particle=p ; spin=s },
{ hole=h'; particle=p'; spin=s'})->
1./.(eps.{h} +. eps.{h'} -. eps.{p} -. eps.{p'})
| _ -> assert false
in
second_order_sum ci i_o1_alfa i_o1_alfa w_alfa

10
CI/Excitation.ml
View File

@ -99,6 +99,16 @@ let double_of_det t t' =
| _ -> assert false | _ -> assert false
let of_det t t' =
match Determinant.degree t t' with
| 0 -> if Determinant.phase t = Determinant.phase t' then
Identity Phase.Pos
else
Identity Phase.Neg
| 1 -> single_of_det t t'
| 2 -> double_of_det t t'
| _ -> multiple_of_det t t'
let pp_s_exc ppf t = let pp_s_exc ppf t =
Format.fprintf ppf "@[T^{%s}_{%d->%d}@]" Format.fprintf ppf "@[T^{%s}_{%d->%d}@]"
(match t.spin with (match t.spin with

14
MOBasis/MOBasis.ml
View File

@ -40,19 +40,19 @@ let two_e_ints t = Lazy.force t.ee_ints
let one_e_ints t = Lazy.force t.one_e_ints let one_e_ints t = Lazy.force t.one_e_ints
let mo_energies t = let mo_energies t =
let m_C = mo_coef t in
let f = let f =
let m_C = mo_coef t in
let m_N = Mat.of_diag @@ mo_occupation t in let m_N = Mat.of_diag @@ mo_occupation t in
let m_P = let m_P = x_o_xt m_N m_C in
gemm m_C @@ (gemm m_N m_C ~transb:`T)
in
match t.mo_type with match t.mo_type with
| RHF -> Fock.make_rhf ~density:m_P (ao_basis t) | RHF -> Fock.make_rhf ~density:m_P (ao_basis t)
| ROHF -> Fock.make_uhf ~density_same:m_P ~density_other:m_P (ao_basis t) | ROHF -> (Mat.scal 0.5 m_P;
Fock.make_uhf ~density_same:m_P ~density_other:m_P (ao_basis t))
| _ -> failwith "Not implemented" | _ -> failwith "Not implemented"
in in
let m_F = Fock.fock f in let m_F0 = Fock.fock f in
Vec.init (size t) (fun i -> m_F.{i,i}) xt_o_x m_F0 m_C
|> Mat.copy_diag
let mo_matrix_of_ao_matrix ~mo_coef ao_matrix = let mo_matrix_of_ao_matrix ~mo_coef ao_matrix =

4
MOBasis/MOBasis.mli
View File

@ -49,6 +49,10 @@ val two_e_ints : t -> ERI.t
val size : t -> int val size : t -> int
(** Number of molecular orbitals in the basis *) (** Number of molecular orbitals in the basis *)
val mo_energies : t -> Vec.t
(** Fock MO energies *)
(** {1 Creators} *) (** {1 Creators} *)
val make : simulation:Simulation.t -> val make : simulation:Simulation.t ->

2
Perturbation/MP2.ml
View File

@ -5,7 +5,7 @@ let make ?(frozen_core=true) hf =
MOBasis.of_hartree_fock hf MOBasis.of_hartree_fock hf
in in
let epsilon = let epsilon =
HartreeFock.eigenvalues hf MOBasis.mo_energies mo_basis
in in
let mo_class = let mo_class =
MOClass.cas_sd mo_basis 0 0 MOClass.cas_sd mo_basis 0 0

14
SCF/HartreeFock.ml
View File

@ -95,8 +95,18 @@ let density t =
| _ -> failwith "Not implemented" | _ -> failwith "Not implemented"
let occupation t = let occupation t =
density t let n_alfa, n_beta =
|> Mat.copy_diag ~n:(Mat.dim2 @@ eigenvectors t) El.n_alfa @@ Simulation.electrons @@ simulation t,
El.n_beta @@ Simulation.electrons @@ simulation t
in
match kind t with
| RHF -> Vec.init (Mat.dim2 @@ eigenvectors t) (fun i ->
if i <= nocc t then 2.0 else 0.0)
| ROHF -> Vec.init (Mat.dim2 @@ eigenvectors t) (fun i ->
if i <= n_beta then 2.0 else
if i <= n_alfa then 1.0 else
0.0)
| _ -> failwith "Not implemented"
let energy t = let energy t =

2
run_cas.ml
View File

@ -72,7 +72,7 @@ let () =
let ci = CI.make space in let ci = CI.make space in
Format.fprintf ppf "CAS-CI energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s); Format.fprintf ppf "CAS-CI energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s);
let pt2 = CI.pt2 ci in 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 s2 = Util.xt_o_x ~o:(CI.s2_matrix ci) ~x:(CI.eigenvectors ci) in let s2 = Util.xt_o_x ~o:(CI.s2_matrix ci) ~x:(CI.eigenvectors ci) in