LCPQ/QCaml
10
1
Fork 0
mirror of https://gitlab.com/scemama/QCaml.git synced 2024-12-22 04:13:33 +01:00
Code Releases Activity

Merged

Browse Source
This commit is contained in:
Anthony Scemama 2019-03-25 09:17:13 +01:00
commit e21aea2b16
12 changed files with 561 additions and 154 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

342
Basis/F12.ml
View File

@ -1,4 +1,342 @@
(** Two-electron integrals over Slater geminals via a fit using Gaussian geminals.
(** Two electron integral functor for operators that are separable among %{ $(x,y,z)$ %}.
It is parameterized by the [zero_m] function.
*)
include TwoElectronIntegralsSeparable
open Constants
let cutoff = integrals_cutoff
module Bs = Basis
module Cs = ContractedShell
module Csp = ContractedShellPair
module Cspc = ContractedShellPairCouple
module Fis = FourIdxStorage
include FourIdxStorage
(** Exponent of the geminal *)
let expo_s = 2.5
(** Coefficients and exponents of the Gaussian fit of the Slater Geminal*)
let coef_g =
[| 0.3144 ; 0.3037 ; 0.1681 ; 0.09811 ; 0.06024 ; 0.03726 |]
let expo_sg_inv =
Array.map (fun x -> 1. /. (x *. expo_s *. expo_s))
[| 0.2209 ; 1.004 ; 3.622 ; 12.16 ; 45.87 ; 254.4 |]
(*
Fit of 1/r:
let coef_g = [|
841.88478132 ; 70.590185207 ; 18.3616020768 ; 7.2608642093 ;
3.57483416444 ; 2.01376031082 ; 1.24216542801 ; 0.81754348620 ;
0.564546514023 ; 0.404228610699 ; 0.297458536575 ; 0.223321219537 ;
0.169933732064 ; 0.130190978230 ; 0.099652303426 ; 0.075428246546 ;
0.0555635614051 ; 0.0386791283055 ; 0.0237550435652 ; 0.010006278387 ;
|]
let expo_sg_inv =
Array.map (fun x -> 1. /. (x *. expo_s *. expo_s))
[| 84135.654509 ; 2971.58727634 ; 474.716025959 ; 130.676724560 ;
47.3938388887 ; 20.2078651631 ; 9.5411021938 ; 4.8109546955 ;
2.52795733067 ; 1.35894103210 ; 0.73586710268 ; 0.39557629706 ;
0.20785895177 ; 0.104809693858 ; 0.049485682527 ; 0.021099788990 ;
0.007652472186 ; 0.0021065225215 ; 0.0003365204879 ; 0.0000118855674 |]
*)
let class_of_contracted_shell_pair_couple shell_pair_couple =
F12RR.contracted_class_shell_pair_couple
expo_sg_inv coef_g shell_pair_couple
module Zero_m = struct
let name = "F12"
end
let filter_contracted_shell_pairs ?(cutoff=integrals_cutoff) shell_pairs =
List.map (fun pair ->
match Cspc.make ~cutoff pair pair with
| Some cspc ->
let cls = class_of_contracted_shell_pair_couple cspc in
(pair, Zmap.fold (fun key value accu -> max (abs_float value) accu) cls 0. )
(* TODO \sum_k |coef_k * integral_k| *)
| None -> (pair, -1.)
) shell_pairs
|> List.filter (fun (_, schwartz_p_max) -> schwartz_p_max >= cutoff)
|> List.map fst
(* TODO
let filter_contracted_shell_pair_couples
?(cutoff=integrals_cutoff) shell_pair_couples =
List.map (fun pair ->
let cls =
class_of_contracted_shell_pairs pair pair
in
(pair, Zmap.fold (fun key value accu -> max (abs_float value) accu) cls 0. )
) shell_pairs
|> List.filter (fun (_, schwartz_p_max) -> schwartz_p_max >= cutoff)
|> List.map fst
*)
let store_class basis ?(cutoff=integrals_cutoff) data contracted_shell_pair_couple cls =
let to_powers x =
let open Zkey in
match to_powers x with
| Three x -> x
| _ -> assert false
in
let shell_p = Cspc.shell_pair_p contracted_shell_pair_couple
and shell_q = Cspc.shell_pair_q contracted_shell_pair_couple
in
let s =
Overlap.of_basis basis
|> Overlap.matrix
in
let lambda_inv = 1. /. expo_s in
Array.iteri (fun i_c powers_i ->
let i_c = Cs.index (Csp.shell_a shell_p) + i_c + 1 in
let xi = to_powers powers_i in
Array.iteri (fun j_c powers_j ->
let j_c = Cs.index (Csp.shell_b shell_p) + j_c + 1 in
let xj = to_powers powers_j in
Array.iteri (fun k_c powers_k ->
let k_c = Cs.index (Csp.shell_a shell_q) + k_c + 1 in
let xk = to_powers powers_k in
Array.iteri (fun l_c powers_l ->
let l_c = Cs.index (Csp.shell_b shell_q) + l_c + 1 in
let xl = to_powers powers_l in
let key = Zkey.of_powers_twelve xi xj xk xl in
let value = Zmap.find cls key in
lambda_inv *. s.{i_c,j_c} *. s.{k_c,l_c} -. value
|> set_chem data i_c j_c k_c l_c
) (Cs.zkey_array (Csp.shell_b shell_q))
) (Cs.zkey_array (Csp.shell_a shell_q))
) (Cs.zkey_array (Csp.shell_b shell_p))
) (Cs.zkey_array (Csp.shell_a shell_p))
let of_basis_serial basis =
let n = Bs.size basis
and shell = Bs.contracted_shells basis
in
let eri_array =
Fis.create ~size:n `Dense
(*
Fis.create ~size:n `Sparse
*)
in
let t0 = Unix.gettimeofday () in
let shell_pairs =
Csp.of_contracted_shell_array shell
|> filter_contracted_shell_pairs ~cutoff
in
Printf.printf "%d significant shell pairs computed in %f seconds\n"
(List.length shell_pairs) (Unix.gettimeofday () -. t0);
let t0 = Unix.gettimeofday () in
let ishell = ref 0 in
List.iter (fun shell_p ->
let () =
if (Cs.index (Csp.shell_a shell_p) > !ishell) then
(ishell := Cs.index (Csp.shell_a shell_p) ; print_int !ishell ; print_newline ())
in
let sp =
Csp.shell_pairs shell_p
in
try
List.iter (fun shell_q ->
let () =
if Cs.index (Csp.shell_a shell_q) >
Cs.index (Csp.shell_a shell_p) then
raise Exit
in
let sq = Csp.shell_pairs shell_q in
let cspc =
if Array.length sp < Array.length sq then
Cspc.make ~cutoff shell_p shell_q
else
Cspc.make ~cutoff shell_q shell_p
in
match cspc with
| Some cspc ->
let cls =
class_of_contracted_shell_pair_couple cspc
in
store_class basis ~cutoff eri_array cspc cls
| None -> ()
) shell_pairs
with Exit -> ()
) shell_pairs ;
Printf.printf "Computed ERIs in %f seconds\n%!" (Unix.gettimeofday () -. t0);
eri_array
(* Parallel functions *)
let of_basis_parallel basis =
let n = Bs.size basis
and shell = Bs.contracted_shells basis
in
let store_class_parallel
?(cutoff=integrals_cutoff) contracted_shell_pair_couple cls =
let to_powers x =
let open Zkey in
match to_powers x with
| Three x -> x
| _ -> assert false
in
let shell_p = Cspc.shell_pair_p contracted_shell_pair_couple
and shell_q = Cspc.shell_pair_q contracted_shell_pair_couple
in
let result = ref [] in
Array.iteri (fun i_c powers_i ->
let i_c = Cs.index (Csp.shell_a shell_p) + i_c + 1 in
let xi = to_powers powers_i in
Array.iteri (fun j_c powers_j ->
let j_c = Cs.index (Csp.shell_b shell_p) + j_c + 1 in
let xj = to_powers powers_j in
Array.iteri (fun k_c powers_k ->
let k_c = Cs.index (Csp.shell_a shell_q) + k_c + 1 in
let xk = to_powers powers_k in
Array.iteri (fun l_c powers_l ->
let l_c = Cs.index (Csp.shell_b shell_q) + l_c + 1 in
let xl = to_powers powers_l in
let key = Zkey.of_powers_twelve xi xj xk xl in
let value = Zmap.find cls key in
result := (i_c, j_c, k_c, l_c, value) :: !result
) (Cs.zkey_array (Csp.shell_b shell_q))
) (Cs.zkey_array (Csp.shell_a shell_q))
) (Cs.zkey_array (Csp.shell_b shell_p))
) (Cs.zkey_array (Csp.shell_a shell_p));
!result
in
let t0 = Unix.gettimeofday () in
let shell_pairs =
Csp.of_contracted_shell_array shell
|> filter_contracted_shell_pairs ~cutoff
in
if Parallel.master then
Printf.printf "%d significant shell pairs computed in %f seconds\n"
(List.length shell_pairs) (Unix.gettimeofday () -. t0);
let t0 = Unix.gettimeofday () in
let ishell = ref max_int in
let input_stream = Stream.of_list (List.rev shell_pairs) in
let f shell_p =
let () =
if Parallel.rank < 2 && Cs.index (Csp.shell_a shell_p) < !ishell then
(ishell := Cs.index (Csp.shell_a shell_p) ; print_int !ishell ; print_newline ())
in
let sp =
Csp.shell_pairs shell_p
in
let result = ref [] in
try
List.iter (fun shell_q ->
let () =
if Cs.index (Csp.shell_a shell_q) >
Cs.index (Csp.shell_a shell_p) then
raise Exit
in
let sq = Csp.shell_pairs shell_q in
let cspc =
if Array.length sp < Array.length sq then
Cspc.make ~cutoff shell_p shell_q
else
Cspc.make ~cutoff shell_q shell_p
in
match cspc with
| Some cspc ->
let cls =
class_of_contracted_shell_pair_couple cspc
in
result := (store_class_parallel ~cutoff cspc cls) :: !result;
| None -> ()
) shell_pairs;
raise Exit
with Exit -> List.concat !result |> Array.of_list
in
let eri_array =
if Parallel.master then
Fis.create ~size:n `Dense
else
Fis.create ~size:0 `Dense
in
let s =
Overlap.of_basis basis
|> Overlap.matrix
in
let lambda_inv = 1. /. expo_s in
Farm.run ~ordered:false ~f input_stream
|> Stream.iter (fun l ->
Array.iter (fun (i_c,j_c,k_c,l_c,value) ->
lambda_inv *. s.{i_c,j_c} *. s.{k_c,l_c} -. value
|> set_chem eri_array i_c j_c k_c l_c) l);
if Parallel.master then
Printf.printf
"Computed %s Integrals in parallel in %f seconds\n%!" Zero_m.name (Unix.gettimeofday () -. t0);
Parallel.broadcast (lazy eri_array)
let of_basis =
match Parallel.size with
| 1 -> of_basis_serial
| _ -> of_basis_parallel

166
CI/CI.ml
View File

@ -410,18 +410,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 is_internal
list_holes list_particles ?(unique=true) 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_particles1 = Array.of_list list_particles1
and list_particles2 = Array.of_list list_particles2
in
let list_holes = Array.of_list list_holes in
let list_particles = Array.of_list list_particles in
let psi0 =
let stream =
Ds.determinant_stream det_space
in
@ -448,7 +443,8 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
|| aux (j-1)
in
aux (i-1)
) else
)
else
is_internal
in
@ -474,12 +470,8 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
accu +. coef.{1} *. (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) ->
let alfa_h_psi alfa =
List.fold_left (fun accu (det, coef) ->
(* Single state here *)
accu +. coef.{1} *. (alfa_o2_i alfa det)) 0. psi_filtered
in
@ -513,30 +505,29 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
let double =
Array.fold_left (fun accu particle' ->
if particle' > particle || particle' = hole then
if particle' >= particle || particle' = hole then
accu
else
accu +.
Array.fold_left (fun accu hole' ->
if hole' = particle' || hole' = particle || hole' < hole then
if hole' = particle' || hole' = particle || hole' <= hole then
accu
else
let alfa =
Determinant.double_excitation
spin hole particle
spin hole' particle' det_i
Determinant.single_excitation
spin hole' particle' alfa
in
if Determinant.is_none alfa ||
already_generated alfa then
accu
else
accu +. w_alfa alfa *. psi_h_alfa_alfa_h_psi alfa
) 0. list_holes1
) 0. list_particles1
) 0. list_holes
) 0. list_particles
in
accu +. single +. double
) 0. list_holes2
) 0. list_particles2
) 0. list_holes
) 0. list_particles
) 0. [ Spin.Alfa ; Spin.Beta ]
in
@ -550,27 +541,28 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
in
if Determinant.is_none alfa then accu else
let double =
let double_ab =
Array.fold_left (fun accu particle' ->
accu +.
Array.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. list_holes1
) 0. list_particles1
accu +.
Array.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. list_holes
) 0. list_particles
in
accu +. double
) 0. list_holes2
) 0. list_particles2
accu +. double_ab
) 0. list_holes
) 0. list_particles
in
same_spin +. opposite_spin
in
@ -581,6 +573,68 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
let second_order_sum2 { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
list_holes list_particles i_o1_alfa e0 psi0 =
let psi0 =
let stream =
Ds.determinant_stream det_space
in
Array.init (Ds.size det_space) (fun i ->
(Stream.next stream), (Mat.copy_row psi0 (i+1)) )
in
let determinants =
Ds.determinants_array det_space
|> Array.to_list
|> List.map (fun det_i ->
[ Spin.Alfa ; Spin.Beta ]
|> List.map (fun spin ->
List.map (fun particle ->
List.map (fun hole ->
[ [ Determinant.single_excitation spin hole particle det_i ] ;
List.map (fun particle' ->
List.map (fun hole' ->
Determinant.double_excitation
spin hole particle
spin hole' particle' det_i
) list_holes
) list_particles
|> List.concat
;
List.map (fun particle' ->
List.map (fun hole' ->
Determinant.double_excitation
spin hole particle
(Spin.other spin) hole' particle' det_i
) list_holes
) list_particles
|> List.concat
]
|> List.concat
) list_holes
) list_particles
|> List.concat
)
|> List.concat
)
|> List.concat
|> List.concat
|> List.filter (fun alfa -> not (Determinant.is_none alfa))
|> List.sort_uniq compare
in
List.fold_left (fun accu alfa ->
let alfa_o2 = i_o1_alfa alfa in
let a_h_psi =
Array.fold_left (fun accu (det,ci) -> ci.{1} *. (alfa_o2 det)) 0. psi0
in
accu +. (a_h_psi *. a_h_psi) /. (e0 -. (alfa_o2 alfa))
) 0. determinants
let is_internal det_space =
let m l =
@ -609,7 +663,7 @@ let is_internal det_space =
Z.logand neg_active_mask beta = occ_mask
let pt2_en ci =
let _pt2_en ci =
let mo_basis = Ds.mo_basis ci.det_space in
let psi0, e0 = Parallel.broadcast ci.eigensystem in
@ -669,11 +723,25 @@ let pt2_en ci =
[ MOClass.active_mos mo_class ; MOClass.virtual_mos mo_class ]
in
second_order_sum ci list_holes list_particles list_holes list_particles
second_order_sum ci list_holes list_particles
(is_internal ci.det_space) i_o1_alfa i_o1_alfa w_alfa psi0
|> List.fold_left (+.) 0.
let pt2_en ci =
let mo_basis = Ds.mo_basis ci.det_space in
let psi0, e0 = Parallel.broadcast ci.eigensystem in
let i_o1_alfa = h_ij mo_basis in
let mo_class = mo_class ci in
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
second_order_sum2 ci list_holes list_particles i_o1_alfa e0.{1} psi0
@ -702,7 +770,7 @@ let pt2_mp ci =
in
let psi0, _ = Parallel.broadcast ci.eigensystem in
second_order_sum ci list_holes list_particles list_holes list_particles
second_order_sum ci list_holes list_particles
(is_internal ci.det_space) i_o1_alfa i_o1_alfa w_alfa psi0
|> List.fold_left (+.) 0.
@ -723,7 +791,7 @@ let variance ci =
[ MOClass.active_mos mo_class ; MOClass.virtual_mos mo_class ]
in
second_order_sum ci list_holes list_particles list_holes list_particles
second_order_sum ci list_holes list_particles
(is_internal ci.det_space) i_o1_alfa i_o1_alfa w_alfa psi0
|> List.fold_left (+.) 0.
@ -737,7 +805,7 @@ let pt2_en_reference ci =
let aux_basis = mo_basis in
let ds =
DeterminantSpace.fci_of_mo_basis ~frozen_core:true aux_basis
DeterminantSpace.fci_of_mo_basis ~frozen_core:false aux_basis
in
let out_dets =
ds
@ -757,7 +825,7 @@ let pt2_en_reference ci =
in
Array.map (fun ki ->
Array.mapi (fun j kj ->
(h_ij aux_basis ki kj) /. (e0.{1} -. h_aa.(j))
(h_ij aux_basis ki kj) /. (e0.{1} -. h_aa.(j))
) out_dets
) in_dets
|> Mat.of_array

16
CI/DeterminantSpace.ml
View File

@ -292,21 +292,21 @@ let arbitrary_of_mo_basis mo_basis f =
let cas_of_mo_basis mo_basis n m =
let cas_of_mo_basis mo_basis ~frozen_core n m =
let f n_alfa =
Ss.cas_of_mo_basis mo_basis n_alfa n m
Ss.cas_of_mo_basis mo_basis ~frozen_core n_alfa n m
in
spin_of_mo_basis mo_basis f
let fci_of_mo_basis ?(frozen_core=true) mo_basis =
let fci_of_mo_basis mo_basis ~frozen_core =
let f n_alfa =
Ss.fci_of_mo_basis ~frozen_core mo_basis n_alfa
Ss.fci_of_mo_basis mo_basis ~frozen_core n_alfa
in
spin_of_mo_basis mo_basis f
let fci_f12_of_mo_basis ?(frozen_core=true) mo_basis mo_num =
let fci_f12_of_mo_basis mo_basis ~frozen_core mo_num =
let s = MOBasis.simulation mo_basis in
let e = Simulation.electrons s in
let n_alfa = Electrons.n_alfa e
@ -321,7 +321,7 @@ let fci_f12_of_mo_basis ?(frozen_core=true) mo_basis mo_num =
(mo_num - n_core)
in
let f n_alfa =
Ss.cas_of_mo_basis mo_basis n_alfa n m
Ss.cas_of_mo_basis mo_basis ~frozen_core n_alfa n m
in
let r =
spin_of_mo_basis mo_basis f
@ -335,9 +335,9 @@ let fci_f12_of_mo_basis ?(frozen_core=true) mo_basis mo_num =
|> MOClass.of_list }
let cas_f12_of_mo_basis mo_basis n m mo_num =
let cas_f12_of_mo_basis mo_basis ~frozen_core n m mo_num =
let f n_alfa =
Ss.cas_of_mo_basis mo_basis n_alfa n m
Ss.cas_of_mo_basis mo_basis ~frozen_core n_alfa n m
in
let r =
spin_of_mo_basis mo_basis f

8
CI/DeterminantSpace.mli
View File

@ -50,20 +50,20 @@ val fock_diag : t -> Determinant.t -> float array * float array
*)
val fci_of_mo_basis : ?frozen_core:bool -> MOBasis.t -> t
val fci_of_mo_basis : MOBasis.t -> frozen_core:bool -> t
(** Creates a space of all possible ways to put [n_alfa] electrons in the {% $\alpha$ %}
[Active] MOs and [n_beta] electrons in the {% $\beta$ %} [Active] MOs.
All other MOs are untouched.
*)
val cas_of_mo_basis : MOBasis.t -> int -> int -> t
val cas_of_mo_basis : MOBasis.t -> frozen_core:bool -> int -> int -> t
(** Creates a CAS(n,m) space of determinants. *)
val fci_f12_of_mo_basis : ?frozen_core:bool -> MOBasis.t -> int -> t
val fci_f12_of_mo_basis : MOBasis.t -> frozen_core:bool -> int -> t
(** Creates the active space to perform a FCI-F12 with an
auxiliary basis set. *)
val cas_f12_of_mo_basis : MOBasis.t -> int -> int -> int -> t
val cas_f12_of_mo_basis : MOBasis.t -> frozen_core:bool -> int -> int -> int -> t
(** [cas_of_mo_basis mo_basis m n mo_num] Creates a CAS(n,m) space
of determinants with an auxiliary basis set defined as the MOs from
[mo_num+1] to [MOBasis.size mo_basis].

151
CI/F12CI.ml
View File

@ -1,6 +1,3 @@
let debug s =
Printf.printf "%s\n%!" s;
open Lacaml.D
type t =
@ -22,7 +19,7 @@ let eigensystem t = Lazy.force t.eigensystem
let f12_integrals mo_basis =
let two_e_ints = MOBasis.f12_ints mo_basis in
( (fun i j _ -> 0.),
( (fun _ _ _ -> 0.),
(fun i j k l s s' ->
if s' = Spin.other s then
F12.get_phys two_e_ints i j k l
@ -43,13 +40,30 @@ let h_ij mo_basis ki kj =
|> List.hd
let f_ij mo_basis ki kj =
let hf_ij mo_basis ki kj =
let integrals =
List.map (fun f -> f mo_basis)
[ CI.h_integrals ; f12_integrals ]
in
CIMatrixElement.make integrals ki kj
let f_ij gamma mo_basis ki kj =
let integrals =
List.map (fun f -> f mo_basis)
[ f12_integrals ]
in
CIMatrixElement.make integrals ki kj
|> List.hd
let integral =
CIMatrixElement.make integrals ki kj
|> List.hd
in
gamma *. integral
(*
match Determinant.degrees ki kj with
| (2,0)
| (0,2) -> 0.5 *. gamma *. integral
| _ -> gamma *. integral
*)
let is_internal det_space =
@ -75,9 +89,7 @@ let is_internal det_space =
Z.logand aux_mask beta = Z.zero
let dressing_vector aux_basis f12_amplitudes ci =
debug "Computing dressing vector";
let dressing_vector gamma aux_basis f12_amplitudes ci =
(*
let i_o1_alfa = h_ij aux_basis in
@ -93,9 +105,6 @@ debug "Computing dressing vector";
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
@ -111,13 +120,37 @@ Util.debug_matrix "f12 amplitudes" f12_amplitudes;
|> 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
|> Array.to_list
in
Printf.printf "Building matrix\n%!";
let m_H_aux, m_F_aux =
let rec col_vecs_list accu_H accu_F = function
| [] ->
List.rev_map Vec.of_list accu_H,
List.rev_map Vec.of_list accu_F
| ki :: rest ->
let h, f =
List.map (fun kj ->
match hf_ij aux_basis ki kj with
| [ a ; b ] -> a, b
| _ -> assert false ) out_dets
|> List.split
in
col_vecs_list (h::accu_H) (f::accu_F) rest
in
let h, f =
col_vecs_list [] [] in_dets
in
Mat.of_col_vecs_list h,
Mat.of_col_vecs_list f
in
(*
let m_H_aux =
Array.map (fun ki ->
Array.map (fun kj ->
@ -130,15 +163,18 @@ Util.debug_matrix "f12 amplitudes" f12_amplitudes;
let m_F_aux =
Array.map (fun ki ->
Array.map (fun kj ->
f_ij aux_basis ki kj
f_ij gamma aux_basis ki kj
) out_dets
) in_dets
|> Mat.of_array
in
*)
Printf.printf "Matrix product\n%!";
let m_HF =
gemm m_H_aux m_F_aux ~transb:`T
in
Printf.printf "Done\n%!";
gemm m_HF f12_amplitudes
|> Matrix.sparse_of_mat
@ -146,9 +182,9 @@ Util.debug_matrix "f12 amplitudes" f12_amplitudes;
let make ~simulation ?(threshold=1.e-12) ?(frozen_core=true) ~mo_basis ~aux_basis_filename () =
let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filename () =
let gamma = 1.0 in
let gamma = 0.5 in
let mo_num = MOBasis.size mo_basis in
@ -188,25 +224,19 @@ let make ~simulation ?(threshold=1.e-12) ?(frozen_core=true) ~mo_basis ~aux_basi
(* While in a sequential region, initiate the parallel
4-idx transformation to avoid nested parallel jobs
*)
debug "Four-idx transform of f12 intergals";
ignore @@ MOBasis.f12_ints aux_basis;
let f = fun ki kj ->
if ki <> kj then
gamma *. (f_ij aux_basis ki kj)
(f_ij gamma aux_basis ki kj)
else
1. +. gamma *. (f_ij aux_basis ki kj)
1. +. (f_ij gamma aux_basis ki kj)
in
debug "Computing F matrix";
let m_F =
CI.create_matrix_spin f det_space
|> 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
@ -223,23 +253,11 @@ 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 aux_basis (f12_amplitudes psi) ci
dressing_vector gamma aux_basis (f12_amplitudes psi) ci
in
(*
Util.debug_matrix "psi" psi;
*)
Format.printf "Amplitude (1,1) : %f@." (f12_amplitudes psi).{1,1};
Format.printf "Dressing vector(1,1) : %f@." (Matrix.get delta 1 1);
let f = 1.0 /. psi.{1,1} in
let delta_00 =
@ -253,67 +271,48 @@ Format.printf "Dressing vector(1,1) : %f@." (Matrix.get delta 1 1);
(*------
TODO SINGLE STATE HERE
*)
(*
let m_H_dressed = Matrix.to_mat m_H in
let f = 1.0 /. psi.{1,1} in
Util.list_range 1 (Mat.dim1 psi)
|> List.iter (fun i ->
let delta = 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;
end
);
let eigenvectors, eigenvalues =
Util.diagonalize_symm m_H_dressed
in
let conv =
1.0 -. abs_float ( dot
(Mat.to_col_vecs psi).(0)
(Mat.to_col_vecs eigenvectors).(0) )
in
Printf.printf "Convergence : %e %f\n" conv (eigenvalues.{1} +. e_shift);
if conv > threshold then
iteration eigenvectors
else
let eigenvalues =
Vec.map (fun x -> x +. e_shift) eigenvalues
in
eigenvectors, eigenvalues
in
iteration ci_coef
*)
(*
------- *)
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 delta.{1,1} *. psi.{1,1} else 0.) )
Vec.init (Matrix.dim1 m_H) (fun i ->
if i = 1 then
Matrix.get m_H i i +. delta.{1,1} *. f
else
Matrix.get m_H i i
)
in
let matrix_prod c =
let w =
Matrix.mm ~transa:`T m_H c
|> Matrix.to_mat
in
let c11 = Matrix.get c 1 1 in
Util.list_range 1 (Mat.dim1 w)
|> List.iter (fun i ->
w.{i,1} <- w.{i,1} +. delta.{i,1} ;
let dci =
delta.{i,1} *. f ;
in
w.{i,1} <- w.{i,1} +. dci *. c11;
if (i <> 1) then
w.{1,1} <- w.{1,1} +. delta.{i,1} *. f *. (Matrix.get c i 1);
w.{1,1} <- w.{1,1} +. dci *. (Matrix.get c i 1);
);
Matrix.dense_of_mat w
in
let eigenvectors, eigenvalues =
Parallel.broadcast (lazy (
Davidson.make ~threshold:1.e-6 ~guess:psi ~n_states diagonal matrix_prod
Davidson.make ~threshold:1.e-6 ~guess:psi ~n_states diagonal matrix_prod
))
in
let conv =
1.0 -. abs_float ( dot
(Mat.to_col_vecs psi).(0)
(Mat.to_col_vecs eigenvectors).(0) )
in
Printf.printf "Convergence : %e %f\n" conv (eigenvalues.{1} +. e_shift);
if Parallel.master then
Printf.printf "F12 Convergence : %e %f\n" conv (eigenvalues.{1} +. e_shift);
if conv > threshold then
iteration eigenvectors
else

6
CI/SpindeterminantSpace.ml
View File

@ -13,7 +13,7 @@ let mo_basis t = t.mo_basis
let mo_class t = t.mo_class
let size t = Array.length t.spin_determinants
let fci_of_mo_basis ?(frozen_core=true) mo_basis elec_num =
let fci_of_mo_basis ~frozen_core mo_basis elec_num =
let mo_num = MOBasis.size mo_basis in
let mo_class = MOClass.fci ~frozen_core mo_basis in
let m l =
@ -35,9 +35,9 @@ let fci_of_mo_basis ?(frozen_core=true) mo_basis elec_num =
{ elec_num ; mo_basis ; mo_class ; spin_determinants }
let cas_of_mo_basis mo_basis elec_num n m =
let cas_of_mo_basis mo_basis ~frozen_core elec_num n m =
let mo_num = MOBasis.size mo_basis in
let mo_class = MOClass.cas_sd mo_basis n m in
let mo_class = MOClass.cas_sd ~frozen_core mo_basis n m in
let m l =
List.fold_left (fun accu i -> let j = i-1 in Z.(logor accu (shift_left one j))
) Z.zero l

4
CI/SpindeterminantSpace.mli
View File

@ -24,12 +24,12 @@ val mo_basis : t -> MOBasis.t
(** {1 Creation} *)
val fci_of_mo_basis : ?frozen_core:bool -> MOBasis.t -> int -> t
val fci_of_mo_basis : frozen_core:bool -> MOBasis.t -> int -> t
(** Create a space of all possible ways to put [n_elec-ncore] electrons in the
[Active] MOs. All other MOs are untouched.
*)
val cas_of_mo_basis : MOBasis.t -> int -> int -> int -> t
val cas_of_mo_basis : MOBasis.t -> frozen_core:bool -> int -> int -> int -> t
(** [cas_of_mo_basis mo_basis n_elec n m] creates a CAS(n,m) space of
[Active] MOs. The unoccupied MOs are [Virtual], and the occupied MOs
are [Core] and [Inactive].

10
MOBasis/MOClass.ml
View File

@ -104,17 +104,19 @@ let fci ?(frozen_core=true) mo_basis =
)
let cas_sd mo_basis n m =
let cas_sd mo_basis ~frozen_core n m =
let mo_num = MOBasis.size mo_basis in
let n_alfa = MOBasis.simulation mo_basis |> Simulation.electrons |> Electrons.n_alfa in
let n_beta = MOBasis.simulation mo_basis |> Simulation.electrons |> Electrons.n_beta in
let n_unpaired = n_alfa - n_beta in
let n_alfa_in_cas = (n - n_unpaired)/2 in
let n_alfa_in_cas = (n - n_unpaired)/2 + n_unpaired in
let last_inactive = n_alfa - n_alfa_in_cas in
let last_active = last_inactive + m in
let ncore =
(Nuclei.small_core @@ Simulation.nuclei @@ MOBasis.simulation mo_basis) / 2
|> min last_inactive
if frozen_core then
(Nuclei.small_core @@ Simulation.nuclei @@ MOBasis.simulation mo_basis) / 2
|> min last_inactive
else 0
in
of_list (
List.concat [

2
MOBasis/MOClass.mli
View File

@ -20,7 +20,7 @@ val fci : ?frozen_core:bool -> MOBasis.t -> t
[n] lowest MOs are [Core] if [frozen_core = true].
*)
val cas_sd: MOBasis.t -> int -> int -> t
val cas_sd: MOBasis.t -> frozen_core:bool -> int -> int -> t
(** [cas_sd mo_basis n m ] creates the MO classes for CAS(n,m) + SD
calculations. lowest MOs are [Core], then all the next MOs are [Inactive],
then [Active], then [Virtual].

4
Perturbation/MP2.ml
View File

@ -1,6 +1,6 @@
type t = float
let make ?(frozen_core=true) hf =
let make ~frozen_core hf =
let mo_basis =
MOBasis.of_hartree_fock hf
in
@ -8,7 +8,7 @@ let make ?(frozen_core=true) hf =
MOBasis.mo_energies mo_basis
in
let mo_class =
MOClass.cas_sd mo_basis 0 0
MOClass.cas_sd mo_basis ~frozen_core 0 0
|> MOClass.to_list
in
let eri =

4
run_cas.ml
View File

@ -67,9 +67,9 @@ let () =
in
let space =
DeterminantSpace.cas_of_mo_basis mos n m
DeterminantSpace.cas_of_mo_basis mos ~frozen_core:true n m
in
let ci = CI.make space in
let ci = CI.make ~frozen_core:true space in
Format.fprintf ppf "CAS-CI energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s);
(*

2
run_mp2.ml
View File

@ -73,7 +73,7 @@ let () =
let e_hf = HartreeFock.energy hf in
let mp2 = MP2.make hf in
let mp2 = MP2.make ~frozen_core:true hf in
Format.fprintf ppf "@[MP2 = %15.10f@]@." mp2;
Format.fprintf ppf "@[E+MP2 = %15.10f@]@." (mp2 +. e_hf)