mirror of https://gitlab.com/scemama/QCaml.git
61 KiB
61 KiB
None
<html>
<head>
</head>
</html>
Test of F12 matrix elements¶
Initialization¶
In [ ]:
#cd "/home/scemama/QCaml";;
#use "topfind";;
#require "jupyter.notebook";;
let png_image name =
Jupyter_notebook.display_file ~base64:true "image/png" ("Notebooks/images/"^name)
;;
#require "lacaml.top";;
#require "alcotest";;
#require "str";;
#require "bigarray";;
#require "zarith";;
#require "getopt";;
#directory "_build";;
#directory "_build/Basis";;
#directory "_build/CI";;
#directory "_build/MOBasis";;
#directory "_build/Nuclei";;
#directory "_build/Parallel";;
#directory "_build/Perturbation";;
#directory "_build/SCF";;
#directory "_build/Utils";;
Modules to load¶
In [ ]:
#load "Constants.cmo";;
#load_rec "Util.cma";;
#load_rec "Matrix.cmo";;
#load_rec "Simulation.cmo";;
#load_rec "Spindeterminant.cmo";;
#load_rec "Determinant.cmo";;
#load_rec "HartreeFock.cmo";;
#load_rec "MOBasis.cmo";;
#load_rec "F12CI.cmo";;
open Lacaml.D;;
Printers¶
In [ ]:
#install_printer AngularMomentum.pp_string ;;
#install_printer Basis.pp ;;
#install_printer Charge.pp ;;
#install_printer Coordinate.pp ;;
#install_printer Vector.pp;;
#install_printer Matrix.pp;;
#install_printer Util.pp_float_2darray;;
#install_printer Util.pp_float_array;;
#install_printer Util.pp_matrix;;
#install_printer HartreeFock.pp ;;
#install_printer Fock.pp ;;
#install_printer MOClass.pp ;;
let pp_mo ppf t = MOBasis.pp ~start:1 ~finish:0 ppf t ;;
#install_printer pp_mo;;
#install_printer DeterminantSpace.pp;;
#install_printer SpindeterminantSpace.pp;;
#install_printer Bitstring.pp;;
Run¶
Simulation¶
In [ ]:
let basis_filename = "/home/scemama/qp2/data/basis/6-31g"
let aux_basis_filename = "/home/scemama/qp2/data/basis/cc-pvdz"
let nuclei = Nuclei.of_zmt_string "be"
let frozen_core = false
let multiplicity = 1
let state = 1
let basis = Basis.of_nuclei_and_basis_filenames ~nuclei [basis_filename]
let aux_basis = Basis.of_nuclei_and_basis_filenames ~nuclei (basis_filename :: aux_basis_filename :: [])
let f12 = F12factor.gaussian_geminal 1.0
let charge = 0
let simulation =
Simulation.make
~f12 ~charge ~multiplicity ~nuclei
~cartesian:true
basis
Hartree-Fock¶
In [ ]:
let hf = HartreeFock.make ~guess:`Hcore ~max_scf:1 simulation ;;
let mo_basis = MOBasis.of_hartree_fock hf
FCI-F12¶
Notations:
- $\langle ij || kl \rangle = \int \phi_i(r_1) \phi_j(r_2) \frac{1}{r_{12}} \phi_k(r1) \phi_l(r2) $
- $\left[ ij || kl \right] = \int \phi_i(r_1) \phi_j(r_2) f_{12} \phi_k(r1) \phi_l(r2) $
Common functions¶
In [ ]:
let f12 = Util.of_some @@ Simulation.f12 simulation
let mo_num = MOBasis.size mo_basis
let pp_spindet = Spindeterminant.pp mo_num
let pp_det = Determinant.pp mo_num
;;
#install_printer pp_spindet ;;
#install_printer pp_det ;;
In [ ]:
let simulation_aux =
let charge = Charge.to_int @@ Simulation.charge simulation
and multiplicity = Electrons.multiplicity @@ Simulation.electrons simulation
and nuclei = Simulation.nuclei simulation
in
let general_basis =
Basis.general_basis @@ Simulation.basis simulation
in
GeneralBasis.combine [
general_basis ; GeneralBasis.read aux_basis_filename
]
|> Basis.of_nuclei_and_general_basis nuclei
|> Simulation.make ~f12 ~charge ~multiplicity ~nuclei
let aux_basis =
MOBasis.of_mo_basis simulation_aux mo_basis
let aux_num =
MOBasis.size aux_basis
In [ ]:
let () = ignore @@ MOBasis.f12_ints aux_basis
let () = ignore @@ MOBasis.two_e_ints aux_basis
Integral-based functions¶
\begin{equation}
\langle I | \hat{H} | J \rangle = \begin{cases}
\sum_i h_{ii} + \frac{1}{2} \sum_{ij} \langle ij || ij \rangle \text{ when } |J\rangle = |I\rangle \\
h_{ik} + \sum_{j} \langle ij || kj \rangle \text{ when } |J\rangle = \hat{T}_i^k |I\rangle \\
\langle ij || kl \rangle \text{ when } |J\rangle = \hat{T}_{ij}^{kl} |I\rangle \\
\end{cases}
\end{equation}\begin{equation}
\langle I | \hat{F} | J \rangle = \begin{cases}
\sum_i f_{ii} + \frac{1}{2} \sum_{ij} \langle ij || ij \rangle \text{ when } |J\rangle = |I\rangle \\
f_{ik} + \sum_{j} \langle ij || kj \rangle \text{ when } |J\rangle = \hat{T}_i^k |I\rangle \\
\langle ij || kl \rangle \text{ when } |J\rangle = \hat{T}_{ij}^{kl} |I\rangle \\
\end{cases}
\end{equation}
In [ ]:
let cancel_singles = false
let mos_cabs =
Util.list_range (mo_num+1) aux_num
let mos_abs =
Util.list_range 1 aux_num
let mos_in =
Util.list_range 1 mo_num
let mos_a k =
Determinant.alfa k
|> Spindeterminant.to_list
let mos_b k =
Determinant.beta k
|> Spindeterminant.to_list
H integrals¶
In [ ]:
let h_one =
let h =
MOBasis.one_e_ints aux_basis
in fun i j _ -> h.{i,j}
let h_two =
let two_e_ints = MOBasis.two_e_ints aux_basis in
let h2 i j k l (s:Spin.t) (s':Spin.t) =
if s' <> s then
ERI.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)
in
h2
F12 integrals¶
In [ ]:
let f_two =
let two_e_ints = MOBasis.f12_ints aux_basis in
let f2 i j k l (s:Spin.t) (s':Spin.t) =
if s' <> s then
0.5 *. F12.get_phys two_e_ints i j k l
else
0.5 *. (
(F12.get_phys two_e_ints i j k l) -.
(F12.get_phys two_e_ints i j l k) )
in
let f3 i j k l (s:Spin.t) (s':Spin.t) =
if (i=k && j<>l) || (j=l && i<>k) then
0.
else
f2 i j k l s s'
in
if cancel_singles then f3 else f2
let f_one = fun _ _ _ -> 0.
(*
let f_two = h_two
let f_one = h_one
*)
Determinant-based functions¶
Integrals¶
In [ ]:
let f12_integrals mo_basis =
( f_one, f_two, None )
let h_ij mo_basis ki kj =
let integrals =
List.map (fun f -> f mo_basis)
[ CI.h_integrals ]
in
CIMatrixElement.make integrals ki kj
|> List.hd
let f_ij mo_basis ki kj =
let integrals =
List.map (fun f -> f mo_basis)
[ f12_integrals ]
in
CIMatrixElement.make integrals ki kj
|> List.hd
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
Determinant space¶
In [ ]:
let is_a_double det_space =
let mo_class = DeterminantSpace.mo_class det_space in
let mo_num = Array.length @@ MOClass.mo_class_array mo_class in
let m l =
List.fold_left (fun accu i ->
let j = i-1 in Bitstring.logor accu (Bitstring.shift_left_one mo_num j)
) (Bitstring.zero mo_num) l
in
let aux_mask = m (MOClass.auxiliary_mos mo_class) in
fun k ->
let alfa =
Determinant.alfa k
|> Spindeterminant.bitstring
in
let beta =
Determinant.beta k
|> Spindeterminant.bitstring
in
let a = Bitstring.logand aux_mask alfa
and b = Bitstring.logand aux_mask beta
in
match Bitstring.popcount a + Bitstring.popcount b with
| 2 | 1 -> true
| 0 | _ -> false
In [ ]:
let in_space =
DeterminantSpace.fci_f12_of_mo_basis aux_basis ~frozen_core mo_num
let aux_space =
DeterminantSpace.fci_of_mo_basis aux_basis ~frozen_core
let det_space_in () =
DeterminantSpace.determinant_stream in_space
let det_space_out () =
let s =
DeterminantSpace.determinant_stream aux_space
in
Stream.from (fun _ ->
try
let is_a_double = is_a_double in_space in
let rec result () =
let ki = Stream.next s in
if is_a_double ki then
Some (ki,ki)
else
result ()
in
result ()
with Stream.Failure -> None
)
In [ ]:
let ci = CI.make ~n_states:state in_space
let ci_coef, ci_energy = Lazy.force ci.eigensystem
let _ = print_newline ()
Permutation operator $p_{12}$ that generates a new determinant with electrons 1 and 2 swapped.
let p12 det_space = let mo_class = DeterminantSpace.mo_class det_space in let mo_num = Array.length @@ MOClass.mo_class_array mo_class in let m l = List.fold_left (fun accu i -> let j = i-1 in Bitstring.logor accu (Bitstring.shift_left_one mo_num j) ) (Bitstring.zero mo_num) l in let aux_mask = m (MOClass.auxiliary_mos mo_class) in let not_aux_mask = Bitstring.(shift_left_one mo_num (mo_num-1) |> minus_one |> logxor aux_mask) in fun k -> let alfa = Determinant.alfa k |> Spindeterminant.bitstring in let beta = Determinant.beta k |> Spindeterminant.bitstring in let a = Bitstring.logand aux_mask alfa and b = Bitstring.logand aux_mask beta in match Bitstring.popcount a, Bitstring.popcount b with | 2, 0 | 0, 2 -> Some (Determinant.negate_phase k) | 1, 1 -> Some (Determinant.of_spindeterminants (Spindeterminant.of_bitstring @@ Bitstring.(logor b (logand not_aux_mask alfa)) ) (Spindeterminant.of_bitstring @@ Bitstring.(logor a (logand not_aux_mask beta)) ) ) | 1, 0 | 0, 1 -> Some k | _ -> None
Matrices $\langle I | H | \alpha \rangle$ and $\langle I | F | \alpha \rangle$¶
In [ ]:
let out_list =
Util.stream_to_list (det_space_out ())
let in_list =
Util.stream_to_list (det_space_in ())
let det_a = Array.of_list out_list
|> Array.map (fun (i,_) -> i)
let det_I = Array.of_list in_list
let m_H_aux, m_F_aux = let rec col_vecs_list accu_H accu_F = function | [] -> List.rev accu_H, List.rev accu_F | (ki, ki') :: rest when ki = ki' -> begin let h, f = List.map (fun kj -> match hf_ij aux_basis kj ki with | [ a ; b ] -> a, b | _ -> assert false ) in_list |> List.split in let h = Vec.of_list h in let f = Vec.of_list f in scal 0.5 f; (* Specific to He singlet *) col_vecs_list (h::accu_H) (f::accu_F) rest end | (ki, ki') :: rest -> begin assert false; end in let h, f = col_vecs_list [] [] out_list in Mat.of_col_vecs_list h, Mat.of_col_vecs_list f
Function to generate all intermediate external determinants $|\alpha \rangle$ between $|I\rangle$ and $|J\rangle$, with a positive phase factor.
excis the degree of excitation between $|I\rangle$ and $|J\rangle$cabsis the number of electrons in the CABSlis the degree of excitation between $|I\rangle$ and $|\alpha \rangle$ris the degree of excitation between $|\alpha \rangle$ and $|J\rangle$
In [ ]:
let generate_alphas ki kj exc cabs l r =
(* Check input excitation degree *)
let d = Determinant.degree ki kj in
if d <> exc then
Printf.sprintf "Invalid excitation degree. Expected %d and found %d." exc d
|> failwith;
(* Generate single excitations *)
let all_singles ki =
let mos_a, mos_b = Determinant.to_lists ki in
[ List.map (fun hole ->
List.map (fun particle ->
if hole = particle then None else
Some (Determinant.single_excitation Spin.Alfa hole particle ki)
) mos_abs
) mos_a
;
List.map (fun hole ->
List.map (fun particle ->
if hole = particle then None else
Some (Determinant.single_excitation Spin.Beta hole particle ki)
) mos_abs
) mos_b
]
|> List.concat
|> List.concat
|> Util.list_some
|> List.filter (fun x -> Determinant.is_none x = false)
|> List.map (fun x -> Determinant.set_phase Phase.Pos x)
in
(* Generate double excitations *)
let all_doubles ki =
let mos_a, mos_b = Determinant.to_lists ki in
[ List.map (fun hole -> (* Alpha-Alpha *)
List.map (fun particle ->
List.map (fun hole' ->
List.map (fun particle' ->
if hole = particle || hole' = particle' then None else
if hole' > hole && particle' < particle then
Some (Determinant.double_excitation Spin.Alfa hole particle Spin.Alfa hole' particle' ki)
else None
) mos_abs
) mos_a
) mos_abs
) mos_a
;
List.map (fun hole -> (* Beta-Beta *)
List.map (fun particle ->
List.map (fun hole' ->
List.map (fun particle' ->
if hole' > hole && particle' < particle then
Some (Determinant.double_excitation Spin.Beta hole particle Spin.Beta hole' particle' ki)
else None
) mos_abs
) mos_b
) mos_abs
) mos_b
;
List.map (fun hole -> (* Alpha-Beta *)
List.map (fun particle ->
List.map (fun hole' ->
List.map (fun particle' ->
Some (Determinant.double_excitation Spin.Alfa hole particle Spin.Beta hole' particle' ki)
) mos_abs
) mos_b
) mos_abs
) mos_a
]
|> List.concat
|> List.concat
|> List.concat
|> List.concat
|> Util.list_some
|> List.filter (fun x -> Determinant.is_none x = false)
|> List.map (fun x -> Determinant.set_phase Phase.Pos x)
in
(* Generate left and right excitations *)
let al =
match l with
| 1 -> all_singles ki
| 2 -> all_doubles ki
| _ -> assert false
in
let ar =
match r with
| 1 -> all_singles kj
| 2 -> all_doubles kj
| _ -> assert false
in
let mo_class = DeterminantSpace.mo_class in_space in
let m l =
List.fold_left (fun accu i ->
let j = i-1 in Bitstring.logor accu (Bitstring.shift_left_one mo_num j)
) (Bitstring.zero mo_num) l
in
let aux_mask = m (MOClass.auxiliary_mos mo_class) in
let good_cabs k =
let alfa =
Determinant.alfa k
|> Spindeterminant.bitstring
in
let beta =
Determinant.beta k
|> Spindeterminant.bitstring
in
let a = Bitstring.logand aux_mask alfa
and b = Bitstring.logand aux_mask beta
in
Bitstring.popcount a + Bitstring.popcount b = cabs
in
let good_lr k =
Determinant.degree ki k = l &&
Determinant.degree k kj = r
in
(* Merge lists in a set of unique determinants *)
List.concat [ al; ar ]
|> List.sort_uniq compare
(* Filter out all determinants with incorrect numbers of electrons in the CABS *)
|> List.filter good_cabs
(* Filter out all determinants with incorrect excitation numbers *)
|> List.filter good_lr
let compute_HaaF ki alphas kj =
List.fold_left (fun accu alpha -> accu
+. h_ij aux_basis ki alpha
*. f_ij aux_basis alpha kj
) 0. alphas
In [ ]:
let check n integral_value exc cabs lexc rexc =
let det_list =
match n with
| 0 -> det_I
| n -> Array.sub det_I 0 n
in
let result =
Printf.printf "Checking ... \n%!";
let percent = ref 0 in
Array.mapi (fun i ki ->
let p = (10 * (i+1))/(Array.length det_list) in
if p <> !percent then
( percent := p ; Printf.printf " - %3d %%\n%!" (p*10) );
Array.mapi (fun j kj ->
if i > j then (0,0,0.,0.) else
if Determinant.degree ki kj <> exc then (0,0,0.,0.) else
let alphas = generate_alphas ki kj exc cabs lexc rexc in
let det_value =
compute_HaaF ki alphas kj
in
let int_value =
integral_value ki kj
in
(i,j,det_value,int_value)
) det_list
) det_list
|> Array.to_list
|> Array.concat
in
let i,j,d,v =
let rec aux k imax jmax emax dmax vmax =
if k = -1 then
imax, jmax, dmax, vmax
else
let i, j, d, v = result.(k) in
let e = abs_float (d -. v) in
if e > emax then
aux (k-1) i j e d v
else
aux (k-1) imax jmax emax dmax vmax
in
aux (n-1) 0 0 0. 0. 0.
in
let error = abs_float (d -. v) in
if error < epsilon_float then
Printf.printf "OK: %e\n%!" error
else
Printf.printf "Failed: (%d, %d) | %e %e | %e\n%!" i j d v error
1. 0 1 11¶
In [ ]:
png_image "0_1_11.png"
$i,j$ : orbital indices of MOs occupied in $|I\rangle$.
$\sum_i \sum_a \hat{T}_{i}^{a} \hat{T}_{a}^{i}$
\begin{align} \sum_{a} \sum_{i} \left( h_{ia} + \sum_{j} \langle i j || a j \rangle + \sum_{\bar{j}} \langle i \bar{j} | a \bar{j} \rangle \right) \left( f_{ai} + \sum_{j} \left[ a j || i j \right] + \sum_{\bar{j}} \left[ a \bar{j} | i \bar{j} \right] \right) + \\ \sum_{\bar{a}} \sum_{\bar{i}} \left( h_{ia} + \sum_{j} \langle \bar{i} j | \bar{a} j \rangle + \sum_{\bar{j}} \langle \bar{i} \bar{j} || \bar{a} \bar{j} \rangle \right) \left( f_{ai} + \sum_{j} \left[ \bar{a} j | \bar{i} j \right] + \sum_{\bar{j}} \left[ \bar{a} \bar{j} || \bar{i} \bar{j} \right] \right) \end{align}In [ ]:
let integral_value ki kj =
(* Alpha-Beta *)
let s, s' = Spin.(Alfa, Beta) in
let mos_a, mos_b = mos_a ki, mos_b ki in
List.fold_left (fun accu a -> accu +.
List.fold_left (fun accu i -> accu +.
(h_one i a s +.
List.fold_left (fun accu j -> accu +. h_two i j a j s s ) 0. mos_a +.
List.fold_left (fun accu j -> accu +. h_two i j a j s s') 0. mos_b ) *.
(f_one a i s +.
List.fold_left (fun accu j -> accu +. f_two a j i j s s ) 0. mos_a +.
List.fold_left (fun accu j -> accu +. f_two a j i j s s') 0. mos_b )
) 0. mos_a
) 0. mos_cabs
+.
List.fold_left (fun accu a -> accu +.
List.fold_left (fun accu i -> accu +.
(h_one i a s +.
List.fold_left (fun accu j -> accu +. h_two i j a j s' s') 0. mos_b +.
List.fold_left (fun accu j -> accu +. h_two i j a j s' s ) 0. mos_a ) *.
(f_one a i s +.
List.fold_left (fun accu j -> accu +. f_two a j i j s' s ) 0. mos_a +.
List.fold_left (fun accu j -> accu +. f_two a j i j s' s') 0. mos_b )
) 0. mos_b
) 0. mos_cabs
let _ =
check 0 integral_value 0 1 1 1
2. 0 1 22¶
In [ ]:
png_image "0_1_22.png"
$i,j$ : orbital indices of MOs occupied in $|I\rangle$.
$k$ : orbital indices of MOs unoccupied in $|I\rangle$.
$\sum_i \sum_a \sum_j \sum_k \hat{T}_{ij}^{ak} \hat{T}_{ak}^{ij}$
\begin{align} \sum_{a} \sum_{k} \sum_{i} \sum_{j<i} \langle i j || a k \rangle \left[ a k || i j \right] + \sum_{a} \sum_{\bar{k}} \sum_{i} \sum_{\bar{j}} \langle i \bar{j} | a \bar{k} \rangle \left[ a \bar{k} | i \bar{j} \right] + \\ \sum_{\bar{a}} \sum_{k} \sum_{\bar{i}} \sum_{j} \langle \bar{i} j | \bar{a} k \rangle \left[ \bar{a} k | \bar{i} j \right] + \sum_{\bar{a}} \sum_{\bar{k}} \sum_{\bar{i}} \sum_{\bar{j}<\bar{i}} \langle \bar{i} \bar{j} || \bar{a} \bar{k} \rangle \left[ \bar{a} \bar{k} || \bar{i} \bar{j} \right] \end{align}In [ ]:
let integral_value ki kj =
(* mos unoccupied in both I and J *)
let mos_virt_a, mos_virt_b =
Array.init mo_num (fun i -> Some (i+1)) ,
Array.init mo_num (fun i -> Some (i+1))
in
List.iter (fun i -> mos_virt_a.(i-1) <- None) (mos_a ki);
List.iter (fun i -> mos_virt_a.(i-1) <- None) (mos_a kj);
List.iter (fun i -> mos_virt_b.(i-1) <- None) (mos_b ki);
List.iter (fun i -> mos_virt_b.(i-1) <- None) (mos_b kj);
let mos_virt_a, mos_virt_b =
Array.to_list mos_virt_a |> Util.list_some,
Array.to_list mos_virt_b |> Util.list_some
in
(* Alpha-Beta *)
let s, s' = Spin.(Alfa, Beta) in
List.fold_left (fun accu j -> accu +.
List.fold_left (fun accu i -> accu +.
List.fold_left (fun accu k -> accu +.
List.fold_left (fun accu a ->
accu +. h_two i j a k s s' *. f_two a k i j s s'
) 0. mos_cabs
) 0. mos_virt_b
) 0. (mos_a ki)
) 0. (mos_b ki)
+.
List.fold_left (fun accu j -> accu +.
List.fold_left (fun accu i -> accu +.
List.fold_left (fun accu k -> accu +.
List.fold_left (fun accu a ->
accu +. h_two i j a k s s' *. f_two a k i j s s'
) 0. mos_cabs
) 0. mos_virt_a
) 0. (mos_b ki)
) 0. (mos_a ki)
+.
(* Alpha-Alpha / Beta-Beta *)
List.fold_left (fun accu (mos_virt,mos,s) ->
let s' = s in accu +.
List.fold_left (fun accu j -> accu +.
List.fold_left (fun accu i -> accu +.
if i < j then accu else
List.fold_left (fun accu k -> accu +.
List.fold_left (fun accu a ->
accu +. h_two i j a k s s' *. f_two a k i j s s'
) 0. mos_cabs
) 0. mos_virt
) 0. mos
) 0. mos
) 0. [ (mos_virt_a,mos_a ki,Spin.Alfa) ; (mos_virt_b,mos_b ki,Spin.Beta)]
let _ =
check 0 integral_value 0 1 2 2
3. 0 2 22¶
In [ ]:
png_image "0_2_22.png"
$i,j$ : orbital indices of MOs occupied in $|I\rangle$.
$\sum_i \sum_a \sum_j \sum_b \hat{T}_{ij}^{ab} \hat{T}_{ab}^{ij}$
\begin{align} \sum_{a} \sum_{b>a} \sum_{i} \sum_{j<i} \langle i j || a b \rangle \left[ a b || i j \right] + \sum_{a} \sum_{\bar{b}} \sum_{i} \sum_{\bar{j}} \langle i \bar{j} | a \bar{b} \rangle \left[ a \bar{b} | i \bar{j} \right] + \\ \sum_{\bar{a}} \sum_{b} \sum_{\bar{i}} \sum_{j} \langle \bar{i} j | \bar{a} b \rangle \left[ \bar{a} b | \bar{i} j \right] + \sum_{\bar{a}} \sum_{\bar{b}>\bar{a}} \sum_{\bar{i}} \sum_{\bar{j}<\bar{i}} \langle \bar{i} \bar{j} || \bar{a} \bar{b} \rangle \left[ \bar{a} \bar{b} || \bar{i} \bar{j} \right] \end{align}In [ ]:
let integral_value ki kj =
(* Alpha-Beta *)
let s, s' = Spin.(Alfa, Beta) in
List.fold_left (fun accu j -> accu +.
List.fold_left (fun accu i -> accu +.
List.fold_left (fun accu b -> accu +.
List.fold_left (fun accu a ->
accu +. h_two i j a b s s' *. f_two a b i j s s'
) 0. mos_cabs
) 0. mos_cabs
) 0. (mos_a ki)
) 0. (mos_b ki)
+.
(* Alpha-Alpha / Beta-Beta *)
List.fold_left (fun accu (mos,s) ->
let s' = s in accu +.
List.fold_left (fun accu j -> accu +.
List.fold_left (fun accu i -> accu +.
if i < j then accu else
List.fold_left (fun accu b -> accu +.
List.fold_left (fun accu a ->
if b > a then accu else
accu +. h_two i j a b s s' *. f_two a b i j s s'
) 0. mos_cabs
) 0. mos_cabs
) 0. mos
) 0. mos
) 0. [ (mos_a ki,Spin.Alfa) ; (mos_b ki,Spin.Beta)]
let _ =
check 0 integral_value 0 2 2 2
4. 1 1 11¶
In [ ]:
png_image "1_1_11.png"
$|J\rangle = \hat{T}_{i}^{k}|I\rangle$
$j$ : orbital indices of MOs occupied both in $|I\rangle$ and $|J\rangle$.
$\sum_a \hat{T}_{i}^{a} \hat{T}_{a}^{k}$
\begin{align} \sum_{a} \left( h_{ia} +\sum_{j} \langle i j || a j \rangle + \sum_{\bar{j}} \langle i \bar{j} | a \bar{j} \rangle \right) \left( f_{ak} +\sum_{j} \left[ a j || k j \right] + \sum_{\bar{j}} \left[ a \bar{j} | k \bar{j} \right]\right) \\ \end{align}In [ ]:
let integral_value ki kj =
let i, k, s, phase =
match Excitation.of_det ki kj with
| Excitation.(Single (phase, { hole ; particle ; spin })) ->
hole, particle, spin, phase
| _ -> assert false
in
let mos_i, mos_i' =
match s with
| Spin.Alfa -> mos_a ki, mos_b ki
| Spin.Beta -> mos_b ki, mos_a ki
in
let mos_j, mos_j' =
match s with
| Spin.Alfa -> mos_a kj, mos_b kj
| Spin.Beta -> mos_b kj, mos_a kj
in
let result =
let s' = Spin.other s in
List.fold_left (fun accu a -> accu +.
(h_one i a s +.
(List.fold_left (fun accu j -> accu +. h_two i j a j s s ) 0. mos_i ) +.
(List.fold_left (fun accu j -> accu +. h_two i j a j s s') 0. mos_i') ) *.
(f_one a k s +.
(List.fold_left (fun accu j -> accu +. f_two a j k j s s ) 0. mos_j ) +.
(List.fold_left (fun accu j -> accu +. f_two a j k j s s') 0. mos_j') )
) 0. mos_cabs
in
match phase with
| Phase.Pos -> result
| Phase.Neg -> -. result
let _ = check 0 integral_value 1 1 1 1
5. 1 1 12¶
In [ ]:
png_image "1_1_12.png"
$|J\rangle = \hat{T}_{i}^{k}|I\rangle$
$\sum_{j} \sum_a \hat{T}_{j}^{a} \hat{T}_{ai}^{jk}$
$j$ : orbital indices of MOs occupied both in $|I\rangle$ and $|J\rangle$.
\begin{align} \sum_{a} \sum_{j} \left( h_{ja} +\sum_{m} \langle j m || a m \rangle + \sum_{\bar{m}} \langle j \bar{m} | a \bar{m} \rangle \right) \left[ a i || j k \right] + \\ \sum_{\bar{a}} \sum_{\bar{j}} \left( h_{\bar{j}\bar{a}} +\sum_{m} \langle \bar{j} m | \bar{a} m \rangle + \sum_{\bar{m}} \langle \bar{j} \bar{m} || \bar{a} \bar{m} \rangle \right) \left[ \bar{a} i | \bar{j} k \right] \end{align}In [ ]:
let integral_value ki kj =
let i, k, s, phase =
match Excitation.of_det ki kj with
| Excitation.(Single (phase, { hole ; particle ; spin })) ->
hole, particle, spin, phase
| _ -> assert false
in
let mos_j, mos_j' =
let alfa =
let i = Spindeterminant.bitstring @@ Determinant.alfa ki in
let j = Spindeterminant.bitstring @@ Determinant.alfa kj in
Bitstring.to_list (Bitstring.logand i j)
in
let beta =
let i = Spindeterminant.bitstring @@ Determinant.beta ki in
let j = Spindeterminant.bitstring @@ Determinant.beta kj in
Bitstring.to_list (Bitstring.logand i j)
in
match s with
| Spin.Alfa -> alfa, beta
| Spin.Beta -> beta, alfa
in
let mos_i, mos_i' =
match s with
| Spin.Alfa -> mos_a ki, mos_b ki
| Spin.Beta -> mos_b ki, mos_a ki
in
let result =
let s' = Spin.other s in
List.fold_left (fun accu a -> accu +.
List.fold_left (fun accu j -> accu +.
(h_one j a s +.
(List.fold_left (fun accu m -> accu +. h_two j m a m s s ) 0. mos_i ) +.
(List.fold_left (fun accu m -> accu +. h_two j m a m s s') 0. mos_i') ) *.
(f_two a i j k s s )
) 0. mos_j
+.
List.fold_left (fun accu j -> accu +.
(h_one j a s +.
(List.fold_left (fun accu m -> accu +. h_two j m a m s' s ) 0. mos_i ) +.
(List.fold_left (fun accu m -> accu +. h_two j m a m s' s') 0. mos_i') ) *.
(f_two a i j k s' s )
) 0. mos_j'
) 0. mos_cabs
in
match phase with
| Phase.Pos -> result
| Phase.Neg -> -. result
let _ = check 0 integral_value 1 1 1 2
6. 1 1 21¶
In [ ]:
png_image "1_1_21.png"
$|J\rangle = \hat{T}_{i}^{k}|I\rangle$
$\sum_{j} \sum_a \hat{T}_{ij}^{ka} \hat{T}_{a}^{j}$
$j$ : orbital indices of MOs occupied both in $|I\rangle$ and $|J\rangle$.
\begin{align} \sum_{a} \sum_{j} \langle i j || k a \rangle \left( f_{aj} +\sum_{m} \left[ m a || m j \right] + \sum_{\bar{m}} \left[ \bar{m} a | \bar{m} j \right] \right) + \\ \sum_{\bar{a}} \sum_{\bar{j}} \langle i \bar{j} | k \bar{a} \rangle \left( f_{\bar{a}\bar{j}} +\sum_{m} \left[ m \bar{a} | m \bar{j} \right] + \sum_{\bar{m}} \left[ \bar{m} \bar{a} || \bar{m} \bar{j} \right] \right) \end{align}In [ ]:
let integral_value ki kj =
let i, k, s, phase =
match Excitation.of_det ki kj with
| Excitation.(Single (phase, { hole ; particle ; spin })) ->
hole, particle, spin, phase
| _ -> assert false
in
let mos_j, mos_j' =
let alfa =
let i = Spindeterminant.bitstring @@ Determinant.alfa ki in
let j = Spindeterminant.bitstring @@ Determinant.alfa kj in
Bitstring.to_list (Bitstring.logand i j)
in
let beta =
let i = Spindeterminant.bitstring @@ Determinant.beta ki in
let j = Spindeterminant.bitstring @@ Determinant.beta kj in
Bitstring.to_list (Bitstring.logand i j)
in
match s with
| Spin.Alfa -> alfa, beta
| Spin.Beta -> beta, alfa
in
let mos_i, mos_i' =
match s with
| Spin.Alfa -> mos_a kj, mos_b kj
| Spin.Beta -> mos_b kj, mos_a kj
in
let result =
let s' = Spin.other s in
List.fold_left (fun accu a -> accu +.
List.fold_left (fun accu j -> accu +.
(h_two i j k a s s ) *.
(f_one a j s +.
(List.fold_left (fun accu m -> accu +. f_two m a m j s s) 0. mos_i ) +.
(List.fold_left (fun accu m -> accu +. f_two m a m j s' s) 0. mos_i') )
) 0. mos_j
+.
List.fold_left (fun accu j -> accu +.
(h_two i j k a s s') *.
(f_one a j s' +.
(List.fold_left (fun accu m -> accu +. f_two m a m j s s') 0. mos_i ) +.
(List.fold_left (fun accu m -> accu +. f_two m a m j s' s') 0. mos_i') )
) 0. mos_j'
) 0. mos_cabs
in
match phase with
| Phase.Pos -> result
| Phase.Neg -> -. result
let _ = check 0 integral_value 1 1 2 1
7. 1 1 22¶
In [ ]:
png_image "1_1_22.png"
$|J\rangle = \hat{T}_{i}^{k}|I\rangle$
$\sum_{j} \sum_{a} \sum_{m} \hat{T}_{ij}^{ma} \hat{T}_{ma}^{kj} + \hat{T}_{lj}^{ka} \hat{T}_{ia}^{lj} + \hat{T}_{ij}^{am} \hat{T}_{am}^{kj} $
$j,l$ : orbital indices of MOs occupied both in $|I\rangle$ and $|J\rangle$.
$m$ : orbital indices of virtual MOs, unoccupied both in $|I\rangle$ and $|J\rangle$.
\begin{align} \sum_{a} \sum_{j} \sum_{m} \langle i j || m a \rangle \left[ m a || k j \right] + \sum_{\bar{a}} \sum_{\bar{j}} \sum_{m} \langle i \bar{j} | m \bar{a} \rangle \left[ m \bar{a} | k \bar{j} \right] + \\ \sum_{a} \sum_{j} \sum_{l} \langle l j || k a \rangle \left[ i a || l j \right] + \sum_{\bar{a}} \sum_{\bar{j}} \sum_{l} \langle l \bar{j} | k \bar{a} \rangle \left[ i \bar{a} | l \bar{j} \right] + \\ \sum_{a} \sum_{\bar{j}} \sum_{\bar{m}} \langle i \bar{j} | a \bar{m} \rangle \left[ a \bar{m} | k \bar{j} \right] \\ \end{align}In [ ]:
let integral_value ki kj =
let i, k, s, phase =
match Excitation.of_det ki kj with
| Excitation.(Single (phase, { hole ; particle ; spin })) ->
hole, particle, spin, phase
| _ -> assert false
in
(* mos unoccupied in both I and J *)
let mos_virt_a, mos_virt_b =
Array.init mo_num (fun i -> Some (i+1)) ,
Array.init mo_num (fun i -> Some (i+1))
in
List.iter (fun i -> mos_virt_a.(i-1) <- None) (mos_a ki);
List.iter (fun i -> mos_virt_a.(i-1) <- None) (mos_a kj);
List.iter (fun i -> mos_virt_b.(i-1) <- None) (mos_b ki);
List.iter (fun i -> mos_virt_b.(i-1) <- None) (mos_b kj);
let mos_virt_a, mos_virt_b =
Array.to_list mos_virt_a |> Util.list_some,
Array.to_list mos_virt_b |> Util.list_some
in
let mos_virt =
match s with
| Spin.Alfa -> mos_virt_a
| Spin.Beta -> mos_virt_b
in
let mos_j, mos_j' =
let alfa =
let i = Spindeterminant.bitstring @@ Determinant.alfa ki in
let j = Spindeterminant.bitstring @@ Determinant.alfa kj in
Bitstring.to_list (Bitstring.logand i j)
in
let beta =
let i = Spindeterminant.bitstring @@ Determinant.beta ki in
let j = Spindeterminant.bitstring @@ Determinant.beta kj in
Bitstring.to_list (Bitstring.logand i j)
in
match s with
| Spin.Alfa -> alfa, beta
| Spin.Beta -> beta, alfa
in
let result =
let s' = Spin.other s in
List.fold_left (fun accu a -> accu +.
List.fold_left (fun accu m -> accu +.
List.fold_left (fun accu j -> accu +.
h_two i j m a s s *. f_two m a k j s s ) 0. mos_j
+.
List.fold_left (fun accu j -> accu +.
h_two i j m a s s' *. f_two m a k j s s' ) 0. mos_j'
) 0. mos_virt
-.
List.fold_left (fun accu l -> accu +.
List.fold_left (fun accu j -> accu +.
h_two l j k a s s *. f_two i a l j s s ) 0. mos_j
+.
List.fold_left (fun accu j -> accu +.
h_two l j k a s s' *. f_two i a l j s s' ) 0. mos_j'
) 0. mos_j
+.
List.fold_left (fun accu m -> accu +.
List.fold_left (fun accu j -> accu +.
h_two i j a m s s' *. f_two a m k j s s' ) 0. mos_j'
) 0. mos_virt
) 0. mos_cabs
in
match phase with
| Phase.Pos -> result
| Phase.Neg -> -. result
let _ = check 20 integral_value 1 1 2 2
In [ ]:
let ki = det_I.(0)
let kj = det_I.(11)
let alpha =
generate_alphas ki kj 1 1 2 2
|> List.filter (fun alpha -> abs_float (compute_HaaF ki [alpha] kj) < 2.0e-05)
|> List.fold_left (fun (a,mx) alpha ->
let x = abs_float (compute_HaaF ki [alpha] kj) in
if x > mx then (alpha, x) else (a,mx) ) (ki,-1.)
|> fst
let _ = Determinant.to_lists ki
let _ = Determinant.to_lists alpha
let _ = Determinant.to_lists kj
let _ =
Format.printf "Excitation: %a\n@." Excitation.pp_exc (Excitation.of_det ki alpha);
Format.printf "Excitation: %a\n@." Excitation.pp_exc (Excitation.of_det alpha kj);
Format.printf "Excitation: %a\n@." Excitation.pp_exc (Excitation.of_det ki kj)
let integral_value ki kj =
let i, k, s, phase =
match Excitation.of_det ki kj with
| Excitation.(Single (phase, { hole ; particle ; spin })) ->
hole, particle, spin, phase
| _ -> assert false
in
(* mos unoccupied in both I and J *)
let mos_virt_a, mos_virt_b =
Array.init mo_num (fun i -> Some (i+1)) ,
Array.init mo_num (fun i -> Some (i+1))
in
List.iter (fun i -> mos_virt_a.(i-1) <- None) (mos_a ki);
List.iter (fun i -> mos_virt_a.(i-1) <- None) (mos_a kj);
List.iter (fun i -> mos_virt_b.(i-1) <- None) (mos_b ki);
List.iter (fun i -> mos_virt_b.(i-1) <- None) (mos_b kj);
let mos_virt_a, mos_virt_b =
Array.to_list mos_virt_a |> Util.list_some,
Array.to_list mos_virt_b |> Util.list_some
in
let mos_virt =
match s with
| Spin.Alfa -> mos_virt_a
| Spin.Beta -> mos_virt_b
in
let mos_j, mos_j' =
let alfa =
let i = Spindeterminant.bitstring @@ Determinant.alfa ki in
let j = Spindeterminant.bitstring @@ Determinant.alfa kj in
Bitstring.to_list (Bitstring.logand i j)
in
let beta =
let i = Spindeterminant.bitstring @@ Determinant.beta ki in
let j = Spindeterminant.bitstring @@ Determinant.beta kj in
Bitstring.to_list (Bitstring.logand i j)
in
match s with
| Spin.Alfa -> alfa, beta
| Spin.Beta -> beta, alfa
in
let result =
let s' = Spin.other s in
let a = 13 in
(*
List.fold_left (fun accu m -> accu +.
List.fold_left (fun accu j -> accu +.
h_two i j m a s s *. f_two m a k j s s ) 0. mos_j
+.
List.fold_left (fun accu j -> accu +.
h_two i j m a s s' *. f_two m a k j s s' ) 0. mos_j'
) 0. mos_virt
+.
*)
let a = 16
and j = 1
and m = 7
in
(*
h_two i j a m s s' *. f_two a m k j s s'
*)
h_two i j m a s s' *. f_two m a k j s s'
(*
let a = 18
and j = 1
and l = 2
in
-. h_two l j k a s s' *. f_two i a l j s s'
*)
in
match phase with
| Phase.Pos -> result
| Phase.Neg -> -. result
let _ = compute_HaaF ki [alpha] kj
let _ = integral_value ki kj
8. 1 2 22a¶
In [ ]:
png_image "1_2_22.png"
$|J\rangle = \hat{T}_{i}^{k}|I\rangle$
$j$ : orbital indices of MOs occupied both in $|I\rangle$ and $|J\rangle$.
$$ \sum_{j} \sum_{b}\sum_{a<b} \langle i j || a b \rangle \left[ ab || kj \right] $$In [ ]:
let integral_value ki kj =
let h, p, s, phase =
match Excitation.of_det ki kj with
| Excitation.(Single (phase, { hole ; particle ; spin })) ->
hole, particle, spin, phase
| _ -> assert false
in
let mos, mos' =
match s with
| Spin.Alfa -> mos_a ki, mos_b ki
| Spin.Beta -> mos_b ki, mos_a ki
in
let result =
(* Alpha-Beta *)
List.fold_left (fun accu j ->
let s' = Spin.other s in
accu +. List.fold_left (fun accu b ->
accu +. List.fold_left (fun accu a ->
accu +. h_two h j a b s s' *. f_two a b p j s s'
) 0. mos_cabs
) 0. mos_cabs
) 0. mos'
+.
(* Alpha-Alpha / Beta-Beta *)
List.fold_left (fun accu j ->
let s' = s in
accu +. List.fold_left (fun accu b ->
accu +. List.fold_left (fun accu a ->
if b > a then accu else
accu +. h_two h j a b s s' *. f_two a b p j s s'
) 0. mos_cabs
) 0. mos_cabs
) 0. mos
in
match phase with
| Phase.Pos -> result
| Phase.Neg -> -. result
let _ = check 50 integral_value 1 2 2 2
9. 2 1 12¶
In [ ]:
png_image "2_1_12.png"
10. 2 1 21¶
In [ ]:
png_image "2_1_21.png"
11. 2 1 22¶
In [ ]:
png_image "2_1_22.png"
12. 2 2 22¶
In [ ]:
png_image "2_2_22.png"
$|J\rangle = \hat{T}_{ij}^{kl}|I\rangle$
$$ \sum_{b}\sum_{a<b} \langle i j || a b \rangle \left[ ab || kj \right] $$In [ ]:
let integral_value ki kj =
let h, h', p, p', s, s', phase =
match Excitation.of_det ki kj with
| Excitation.(Double (phase,
{ hole=h ; particle=p ; spin=s },
{ hole=h'; particle=p'; spin=s'}) ) -> h, h', p, p', s, s', phase
| _ -> assert false
in
let result =
if s <> s' then (* Alpha-Beta *)
List.fold_left (fun accu b ->
accu +. List.fold_left (fun accu a ->
accu +. h_two h h' a b s s' *. f_two a b p p' s s'
) 0. mos_cabs
) 0. mos_cabs
else (* Alpha-Alpha / Beta-Beta *)
List.fold_left (fun accu b ->
accu +. List.fold_left (fun accu a ->
if b > a then accu else
accu +. h_two h h' a b s s' *. f_two a b p p' s s'
) 0. mos_cabs
) 0. mos_cabs
in
match phase with
| Phase.Pos -> result
| Phase.Neg -> -. result
let _ = check 100 integral_value 2 2 2 2
13. 3 1 2 2¶
In [ ]:
png_image "3_1_22.png"
$|J\rangle = \hat{T}_{ijm}^{kln}|I\rangle$
$$ \sum_{a} \langle i j || k a \rangle \left[ m a || n l \right] $$In [ ]:
let integral_value ki kj =
let h1, h2, h3, p1, p2, p3, s1, s2, s3, phase =
match Excitation.of_det ki kj with
| Excitation.(Triple (phase,
{ hole=h1 ; particle=p1 ; spin=s1 },
{ hole=h2 ; particle=p2 ; spin=s2 },
{ hole=h3 ; particle=p3 ; spin=s3 }) ) -> h1, h2, h3, p1, p2, p3, s1, s2, s3, phase
| _ -> assert false
in
(*
123
213
231
321
312
132
*)
let result =
List.fold_left (fun accu a -> accu +.
h_two h1 h2 p1 a s1 s2 *. f_two h3 a p3 p2 s3 s2 -.
h_two h2 h1 p2 a s2 s1 *. f_two h3 a p3 p1 s3 s1 +.
h_two h2 h3 p2 a s2 s3 *. f_two h1 a p1 p3 s1 s3 -.
h_two h3 h2 p3 a s3 s2 *. f_two h1 a p1 p2 s1 s2 +.
h_two h3 h1 p3 a s3 s1 *. f_two h2 a p2 p1 s2 s1 -.
h_two h1 h3 p1 a s1 s3 *. f_two h2 a p2 p3 s2 s3
) 0. mos_cabs
in
match phase with
| Phase.Pos -> result
| Phase.Neg -> -. result
let _ = check 100 integral_value 3 1 2 2
In [ ]: