From e7780804bf3e84e21b718d70bf68eb929513db0f Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 18 Mar 2019 23:38:01 +0100 Subject: [PATCH] CAS+EN OK --- CI/CI.ml | 424 +++++++++++++++++++++++++++----------------- CI/Excitation.ml | 10 ++ MOBasis/MOBasis.ml | 14 +- MOBasis/MOBasis.mli | 4 + Perturbation/MP2.ml | 2 +- SCF/HartreeFock.ml | 14 +- run_cas.ml | 2 +- 7 files changed, 294 insertions(+), 176 deletions(-) diff --git a/CI/CI.ml b/CI/CI.ml index bc6ebe6..b6531d0 100644 --- a/CI/CI.ml +++ b/CI/CI.ml @@ -94,45 +94,45 @@ let create_matrix_arbitrary f det_space = in let task (i,i_dets) = - let shift = index_start.(i) in + let shift = index_start.(i) in - let result = - Array.init (index_start.(i+1) - shift) - (fun _ -> []) - 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 + (** 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 + 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 + 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 -> + 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) + 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 -> () + with Exit -> () ) i_dets - | 1 -> + | 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 @@ -145,18 +145,18 @@ let create_matrix_arbitrary f det_space = match compare js j_dets.(j') with | -1 -> aux r_singles j' | 0 -> - let kj = - Determinant.of_spindeterminants j_alfa j_beta - in (update + 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);) + 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 -> + | 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 @@ -169,27 +169,27 @@ let create_matrix_arbitrary f det_space = match compare js j_dets.(j') with | -1 -> aux r_doubles j' | 0 -> - let kj = - Determinant.of_spindeterminants j_alfa j_beta - in (update + 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);) + 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 + | _ -> (); + end + ) det; + let r = + Array.map (fun l -> + List.rev l + |> Vector.sparse_of_assoc_list ndet + ) result + in (i,r) + in let result = @@ -259,63 +259,63 @@ 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 + (** 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 -> - let degree_a = deg_a j_alfa in - begin - match degree_a with - | 2 -> - let i' = ref 0 in - List.iteri (fun ib i_beta -> - let ki = Determinant.of_spindeterminants i_alfa i_beta in - let kj = Determinant.of_spindeterminants j_alfa i_beta in - update !i' (ib + !j) ki kj; - incr i'; - ) b; - | 1 -> - let i' = ref 0 in - List.iteri (fun ib i_beta -> - let ki = Determinant.of_spindeterminants i_alfa i_beta in - let singles, _ = degree_bb.(ib) in - List.iter (fun (j', j_beta) -> - let kj = Determinant.of_spindeterminants j_alfa j_beta in - update !i' (j' + !j) ki kj - ) singles; - incr i'; - ) b; - | 0 -> - let i' = ref 0 in - List.iteri (fun ib i_beta -> - let ki = Determinant.of_spindeterminants i_alfa i_beta in - let _singles, doubles = degree_bb.(ib) in - List.iter (fun (j', j_beta) -> - let kj = Determinant.of_spindeterminants j_alfa j_beta in - update !i' (j' + !j) ki kj - ) doubles; - incr i'; - ) b; - | _ -> (); - end; - j := !j + n_beta - ) a; - let r = - Array.map (fun l -> - List.rev l - |> Vector.sparse_of_assoc_list ndet - ) result - in (i,r) + let j = ref 1 in + let deg_a = Spindeterminant.degree i_alfa in + List.iter (fun j_alfa -> + let degree_a = deg_a j_alfa in + begin + match degree_a with + | 2 -> + let i' = ref 0 in + List.iteri (fun ib i_beta -> + let ki = Determinant.of_spindeterminants i_alfa i_beta in + let kj = Determinant.of_spindeterminants j_alfa i_beta in + update !i' (ib + !j) ki kj; + incr i'; + ) b; + | 1 -> + let i' = ref 0 in + List.iteri (fun ib i_beta -> + let ki = Determinant.of_spindeterminants i_alfa i_beta in + let singles, _ = degree_bb.(ib) in + List.iter (fun (j', j_beta) -> + let kj = Determinant.of_spindeterminants j_alfa j_beta in + update !i' (j' + !j) ki kj + ) singles; + incr i'; + ) b; + | 0 -> + let i' = ref 0 in + List.iteri (fun ib i_beta -> + let ki = Determinant.of_spindeterminants i_alfa i_beta in + let _singles, doubles = degree_bb.(ib) in + List.iter (fun (j', j_beta) -> + let kj = Determinant.of_spindeterminants j_alfa j_beta in + update !i' (j' + !j) ki kj + ) doubles; + incr i'; + ) b; + | _ -> (); + end; + j := !j + n_beta + ) a; + let r = + Array.map (fun l -> + List.rev l + |> Vector.sparse_of_assoc_list ndet + ) result + in (i,r) in let result = @@ -348,7 +348,7 @@ let make ?(n_states=1) det_space = (* While in a sequential region, initiate the parallel 4-idx transformation to avoid nested parallel jobs - *) + *) ignore @@ MOBasis.two_e_ints mo_basis; let f = @@ -385,32 +385,33 @@ 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_class = DeterminantSpace.mo_class det_space in let mo_indices = let cls = MOClass.mo_class_array mo_class in Util.list_range 1 (MOBasis.size mo_basis) |> List.filter (fun i -> match cls.(i) with - | MOClass.Deleted _ - | MOClass.Core _ -> false - | _ -> true + | MOClass.Deleted _ + | MOClass.Core _ -> false + | _ -> true ) in -(* Only the gournd state is computed here *) - let psi0, e0 = Lazy.force eigensystem in - let psi0 = + let psi0, _ = Lazy.force eigensystem in + let stream = Ds.determinant_stream det_space in Array.init (Ds.size det_space) (fun i -> - Stream.next stream, psi0.{i+1,1}) + Stream.next stream, psi0.{i+1,1}) in - let e0 = e0.{1} in + (* let is_internal = @@ -440,16 +441,16 @@ let pt2 { det_space ; m_H ; m_S2 ; eigensystem ; n_states } = in *) - + let det_contribution i = let psi_filtered_idx = let rec aux accu = function - | j when j <= i -> List.rev accu - | j -> if Determinant.degree (fst psi0.(i)) (fst psi0.(j)) < 4 then - aux (j::accu) (j-1) - else - aux accu (j-1) + | j when j < i -> List.rev accu + | j -> if Determinant.degree (fst psi0.(i)) (fst psi0.(j)) < 4 then + aux (j::accu) (j-1) + else + aux accu (j-1) in aux [] (Array.length psi0 - 1) 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 in - let psi_h alfa = - let hij = h_ij mo_basis alfa in + let symmetric = i_o1_alfa == alfa_o2_i in + + let psi_h_alfa alfa = 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 let is_internal alfa = - let rec aux = function + let rec aux = function | -1 -> false | j -> Determinant.degree (fst psi0.(j)) alfa = 0 || aux (j-1) - in - aux (Array.length psi0 - 1) + in + aux (Array.length psi0 - 1) in let already_generated alfa = @@ -476,60 +494,136 @@ let pt2 { det_space ; m_H ; m_S2 ; eigensystem ; n_states } = true else let rec aux = function - | -1 -> false - | j -> Determinant.degree (fst psi0.(j)) alfa <= 2 || aux (j-1) + | -1 -> false + | j -> Determinant.degree (fst psi0.(j)) alfa <= 2 || aux (j-1) in aux (i-1) in let det_i = fst psi0.(i) in + let w_alfa = w_alfa det_i in - List.fold_left (fun accu particle -> - accu +. - List.fold_left (fun accu hole -> - if hole = particle then accu else + let same_spin = + List.fold_left (fun accu spin -> + accu +. + List.fold_left (fun accu particle -> accu +. - List.fold_left (fun accu spin -> - let alfa = - Determinant.single_excitation spin hole particle det_i - in - 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) + List.fold_left (fun accu hole -> + if hole = particle then accu else + let alfa = + Determinant.single_excitation spin 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' -> + let single = + if already_generated alfa then 0. else + w_alfa alfa *. psi_h_alfa_alfa_h_psi alfa + in + + let double = + List.fold_left (fun accu particle' -> + if particle' > particle then + accu + else accu +. - List.fold_left (fun accu spin' -> - let alfa = - Determinant.double_excitation - spin hole particle - spin' hole' particle' det_i - in - if Determinant.is_none alfa || - already_generated alfa then + List.fold_left (fun accu hole' -> + if hole' = particle' || hole' < hole then accu else - let h_aa = h_ij mo_basis alfa alfa in - let psi_h_alfa = psi_h alfa in - accu +. psi_h_alfa *. psi_h_alfa /. (e0 -. h_aa) - ) 0. [ Spin.Alfa ; Spin.Beta ] - ) 0. mo_indices + let alfa = + Determinant.double_excitation + spin hole particle + spin 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 +. single +. double - ) 0. [ Spin.Alfa ; Spin.Beta ] - ) 0. mo_indices - ) 0. mo_indices + in + accu +. single +. double + ) 0. mo_indices + ) 0. mo_indices + ) 0. [ Spin.Alfa ; Spin.Beta ] 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.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 + + diff --git a/CI/Excitation.ml b/CI/Excitation.ml index 25fe73d..c1fdcc6 100644 --- a/CI/Excitation.ml +++ b/CI/Excitation.ml @@ -99,6 +99,16 @@ let double_of_det t t' = | _ -> 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 = Format.fprintf ppf "@[T^{%s}_{%d->%d}@]" (match t.spin with diff --git a/MOBasis/MOBasis.ml b/MOBasis/MOBasis.ml index d32d457..4e2e8e7 100644 --- a/MOBasis/MOBasis.ml +++ b/MOBasis/MOBasis.ml @@ -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 mo_energies t = + let m_C = mo_coef t in let f = - let m_C = mo_coef t in let m_N = Mat.of_diag @@ mo_occupation t in - let m_P = - gemm m_C @@ (gemm m_N m_C ~transb:`T) - in + let m_P = x_o_xt m_N m_C in match t.mo_type with | 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" in - let m_F = Fock.fock f in - Vec.init (size t) (fun i -> m_F.{i,i}) + let m_F0 = Fock.fock f in + xt_o_x m_F0 m_C + |> Mat.copy_diag let mo_matrix_of_ao_matrix ~mo_coef ao_matrix = diff --git a/MOBasis/MOBasis.mli b/MOBasis/MOBasis.mli index 124e531..25ef20c 100644 --- a/MOBasis/MOBasis.mli +++ b/MOBasis/MOBasis.mli @@ -49,6 +49,10 @@ val two_e_ints : t -> ERI.t val size : t -> int (** Number of molecular orbitals in the basis *) +val mo_energies : t -> Vec.t +(** Fock MO energies *) + + (** {1 Creators} *) val make : simulation:Simulation.t -> diff --git a/Perturbation/MP2.ml b/Perturbation/MP2.ml index eb3ef7e..98f481d 100644 --- a/Perturbation/MP2.ml +++ b/Perturbation/MP2.ml @@ -5,7 +5,7 @@ let make ?(frozen_core=true) hf = MOBasis.of_hartree_fock hf in let epsilon = - HartreeFock.eigenvalues hf + MOBasis.mo_energies mo_basis in let mo_class = MOClass.cas_sd mo_basis 0 0 diff --git a/SCF/HartreeFock.ml b/SCF/HartreeFock.ml index 1c8cfe7..e2039a6 100644 --- a/SCF/HartreeFock.ml +++ b/SCF/HartreeFock.ml @@ -95,8 +95,18 @@ let density t = | _ -> failwith "Not implemented" let occupation t = - density t - |> Mat.copy_diag ~n:(Mat.dim2 @@ eigenvectors t) + let n_alfa, n_beta = + 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 = diff --git a/run_cas.ml b/run_cas.ml index 05b7f90..6c77ab9 100644 --- a/run_cas.ml +++ b/run_cas.ml @@ -72,7 +72,7 @@ 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 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) (* let s2 = Util.xt_o_x ~o:(CI.s2_matrix ci) ~x:(CI.eigenvectors ci) in