LCPQ/QCaml

Fork 0

mirror of https://gitlab.com/scemama/QCaml.git synced 2024-06-19 19:52:06 +02:00

Anthony Scemama 537c5540c4 Notebook 0 2 22 and 1 2 22 OK

2019-12-04 18:13:54 +01:00

29 KiB

Raw Blame History

None <html> <head> </head>

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  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_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.

exc is the degree of excitation between $|I\rangle$ and $|J\rangle$
cabs is the number of electrons in the CABS
l is the degree of excitation between $|I\rangle$ and $|\alpha \rangle$
r is 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 ->
                  Determinant.single_excitation Spin.Alfa hole particle ki
              ) mos_cabs
          ) mos_a
      ;
        List.map (fun hole -> 
          List.map (fun particle ->
                  Determinant.single_excitation Spin.Beta hole particle ki
              ) mos_cabs
          ) mos_b 
      ]
      |> List.concat
      |> List.concat
      |> 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' > hole && particle' < particle then
                  Some (Determinant.double_excitation Spin.Alfa hole particle Spin.Alfa hole' particle' ki)
                else None
                      ) mos_cabs
                  ) mos_a 
              ) mos_cabs
          ) 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_cabs
                  ) mos_b 
              ) mos_cabs
          ) 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_cabs
                  ) mos_b 
              ) mos_cabs
          ) 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
    Printf.printf "Checking...\n%!";
    let result =
        Array.mapi (fun i ki -> 
            Array.mapi (fun j kj -> 
                if i > j then 0. else
                if Determinant.degree ki kj <> exc then 0. else
                let alphas = generate_alphas ki kj exc cabs lexc rexc in
                let det_value =
                  compute_HaaF ki alphas kj
                in
                integral_value ki kj -. det_value 
            ) det_list
        ) det_list
        |> Array.to_list
        |> Array.concat
  in
  let error = 
      Array.fold_left (fun accu x -> accu +. x *. x) 0. result
      /. (float_of_int @@ Array.length result)
  in
  if error < epsilon_float then
      Printf.printf "OK: %e\n%!" error
  else
      Printf.printf "Failed: %e\n%!" error

0 2 22¶

In [ ]:

png_image "0_2_22.png"

$i,j$ : orbital indices of MOs occupied in $|I\rangle$.

$$ \sum_{j} \sum_{i<j} \sum_{b}\sum_{a<b} \langle i j || a b \rangle \left[ ab || ij \right] $$

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 10 integral_value 0 2 2 2

1 2 22¶

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

   (* Alpha-Beta *)
   let s' = Spin.other s in
   List.fold_left (fun accu j -> 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 *)
     let s' = s in 
     List.fold_left (fun accu j -> 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 
     

let _ = check 10 integral_value 1 2 2 2

In [ ]:

</html>

29 KiB Raw Blame History

Test of F12 matrix elements¶

Initialization¶

Modules to load¶

Printers¶

Run¶

Simulation¶

Hartree-Fock¶

FCI-F12¶

Common functions¶

Integral-based functions¶

H integrals¶

F12 integrals¶

Determinant-based functions¶

Integrals¶

Determinant space¶

Matrices $\langle I | H | \alpha \rangle$ and $\langle I | F | \alpha \rangle$¶

0 2 22¶

1 2 22¶

29 KiB

Raw Blame History