diff --git a/CI/F12CI.ml b/CI/F12CI.ml index f46d311..d3c8430 100644 --- a/CI/F12CI.ml +++ b/CI/F12CI.ml @@ -24,7 +24,11 @@ let eigensystem t = Lazy.force t.eigensystem let hf_ij_non_zero hf12_integrals deg_a deg_b ki kj = let integrals = [ let one_e _ _ _ = 0. in - let hf12_s, hf12_o, hf12m_s, hf12m_o = hf12_integrals in + let { HF12. + simulation ; aux_basis ; + hf12 ; hf12_anti ; + hf12_single ; hf12_single_anti } = hf12_integrals + in let kia = De.alfa ki and kib = De.beta ki and kja = De.alfa kj and kjb = De.beta kj in @@ -37,14 +41,14 @@ let hf_ij_non_zero hf12_integrals deg_a deg_b ki kj = in let two_e i j k l s s' = if s' = Spin.other s then - hf12_o.{i,j,k,l} -. 1. *. ( - (List.fold_left (fun accu m -> accu +. hf12m_o.{m,i,j,k,l}) 0. mo_a) +. - (List.fold_left (fun accu m -> accu +. hf12m_o.{m,j,i,l,k}) 0. mo_b) + hf12.{i,j,k,l} -. 1. *. ( + (List.fold_left (fun accu m -> accu +. hf12_single.{m,i,j,k,l}) 0. mo_a) +. + (List.fold_left (fun accu m -> accu +. hf12_single.{m,j,i,l,k}) 0. mo_b) ) else - hf12_s.{i,j,k,l} -. 1. *. ( - (List.fold_left (fun accu m -> accu +. hf12m_s.{m,i,j,k,l}) 0. mo_a) +. - (List.fold_left (fun accu m -> accu +. hf12m_s.{m,j,i,l,k}) 0. mo_b) + hf12_anti.{i,j,k,l} -. 1. *. ( + (List.fold_left (fun accu m -> accu +. hf12_single_anti.{m,i,j,k,l}) 0. mo_a) +. + (List.fold_left (fun accu m -> accu +. hf12_single_anti.{m,j,i,l,k}) 0. mo_b) ) in (one_e, two_e) @@ -128,8 +132,6 @@ let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filen let eigenvectors, eigenvalues = -(* Column dressing -*) let delta = lacpy delta in Mat.scal f delta; for k=1 to state-1 do @@ -166,65 +168,10 @@ let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filen in -(* Diagonal dressing *) -(* - let diagonal = - Vec.init (Matrix.dim1 m_H) (fun i -> - Matrix.get m_H i i +. - if (abs_float psi.{i,state} > 1.e-8) then - delta.{i,state} /. psi.{i,state} - else 0. - ) - in - - let matrix_prod c = - let w = - Matrix.mm ~transa:`T m_H c - |> Matrix.to_mat - in - for i=1 to (Mat.dim1 w) do - w.{i,state} <- w.{i,state} +. delta.{i,state} - done; - Matrix.dense_of_mat w - in -*) - - Parallel.broadcast (lazy ( Davidson.make ~threshold:1.e-9 ~guess:psi ~n_states:state diagonal matrix_prod )) -(* - let m_H = Matrix.to_mat m_H |> lacpy in -*) - -(* DIAGONAL TEST - for i=1 to Mat.dim1 m_H do - if (abs_float psi.{i,state} > 1.e-8) then - m_H.{i,i} <- m_H.{i,i} +. delta.{i,state} /. psi.{i,state}; - done; -*) - - -(* COLUMN TEST - for i=1 to Mat.dim1 m_H do - let d = delta.{i,state} /. psi.{column_idx,state} in - m_H.{i,column_idx} <- m_H.{i,column_idx} +. d; - if (i <> column_idx) then - begin - m_H.{column_idx,i} <- m_H.{column_idx,i} +. d; - m_H.{column_idx,column_idx} <- m_H.{column_idx,column_idx} -. - d *. psi.{i,state} - end - done; -*) - - - -(* - let m_v = syev m_H in - m_H, m_v -*) in @@ -240,6 +187,14 @@ let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filen Printf.printf "F12 Convergence : %e %f\n" conv (eigenvalues.{state} +. Simulation.nuclear_repulsion simulation); +(* + let cabs_singles = + let f = + Fock.make_rhf ~density ~ao_basis:large_ao_basis + in + in +*) + if conv > threshold then iteration ~state eigenvectors else diff --git a/MOBasis/HF12.ml b/MOBasis/HF12.ml index d514ec9..df8d77a 100644 --- a/MOBasis/HF12.ml +++ b/MOBasis/HF12.ml @@ -4,10 +4,25 @@ open Lacaml.D module Fis = FourIdxStorage -type t = (float, Bigarray.float64_elt, Bigarray.fortran_layout) Bigarray.Genarray.t - * (float, Bigarray.float64_elt, Bigarray.fortran_layout) Bigarray.Genarray.t - * (float, Bigarray.float64_elt, Bigarray.fortran_layout) Bigarray.Genarray.t - * (float, Bigarray.float64_elt, Bigarray.fortran_layout) Bigarray.Genarray.t + +type q = (float, Bigarray.float64_elt, Bigarray.fortran_layout) Bigarray.Genarray.t + +type t = { + simulation : Simulation.t ; + aux_basis : MOBasis.t ; + hf12 : q ; (* hf12.{i,j,k,l} = \sum_{ab} + + \sum_{am} + + \sum_{bm} + *) + hf12_anti: q ; (* hf12_anti.{i,j,k,l} = + \sum_{ab} ( - ) ( - ) + + \sum_{am} ( - ) ( - ) + + \sum_{bm} ( - ) ( - ) + *) + hf12_single : q ; (* hf12.{m,i,j,k,l} = \sum_{a} *) + hf12_single_anti: q ; (* hf12_anti.{m,i,j,k,l} = + \sum_{ab} ( - ) ( - ) *) +} let make ~simulation ~mo_basis ~aux_basis_filename () = @@ -16,22 +31,23 @@ let make ~simulation ~mo_basis ~aux_basis_filename () = let mo_num = MOBasis.size mo_basis in (* Add auxiliary basis set *) - let aux_basis = - let s = - 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 simulation = + let charge = Charge.to_int @@ Simulation.charge simulation + and multiplicity = Electrons.multiplicity @@ Simulation.electrons simulation + and nuclei = Simulation.nuclei simulation in - MOBasis.of_mo_basis s mo_basis + 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 + in + + let aux_basis = + MOBasis.of_mo_basis simulation mo_basis in let aux_num = MOBasis.size aux_basis in @@ -47,11 +63,17 @@ let make ~simulation ~mo_basis ~aux_basis_filename () = (* Compute the integrals *) if Parallel.master then Printf.eprintf "Computing HF12 integrals\n%!"; - let result_s, result_o, resultm_s, resultm_o = - Bigarray.Genarray.create Float64 Bigarray.fortran_layout [| mo_num ; mo_num ; mo_num ; mo_num |] , - Bigarray.Genarray.create Float64 Bigarray.fortran_layout [| mo_num ; mo_num ; mo_num ; mo_num |] , - Bigarray.Genarray.create Float64 Bigarray.fortran_layout [| mo_num ; mo_num ; mo_num ; mo_num ; mo_num |] , - Bigarray.Genarray.create Float64 Bigarray.fortran_layout [| mo_num ; mo_num ; mo_num ; mo_num ; mo_num |] + let hf12, hf12_anti, hf12_single, hf12_single_anti = + let create4 n = + Bigarray.Genarray.create Float64 Bigarray.fortran_layout [| n ; n ; n ; n |] + in + let create5 n = + Bigarray.Genarray.create Float64 Bigarray.fortran_layout [| n ; n ; n ; n ; n |] + in + create4 mo_num , + create4 mo_num , + create5 mo_num , + create5 mo_num in @@ -200,11 +222,11 @@ let make ~simulation ~mo_basis ~aux_basis_filename () = |> Stream.iter (fun (hf_s, hf_o, hfm_s, hfm_o, j, l) -> for k=1 to mo_num do for i=1 to mo_num do - result_s.{i,j,k,l} <- hf_s.{i,k} ; - result_o.{i,j,k,l} <- hf_o.{i,k} ; + hf12.{i,j,k,l} <- hf_o.{i,k} ; + hf12_anti.{i,j,k,l} <- hf_s.{i,k} ; for m=1 to mo_num do - resultm_s.{m,i,j,k,l} <- hfm_s.{m,i,k} ; - resultm_o.{m,i,j,k,l} <- hfm_o.{m,i,k} + hf12_single.{m,i,j,k,l} <- hfm_o.{m,i,k} ; + hf12_single_anti.{m,i,j,k,l} <- hfm_s.{m,i,k} done done done ); @@ -224,7 +246,13 @@ let make ~simulation ~mo_basis ~aux_basis_filename () = Printf.printf "%!"; *) - Parallel.broadcast (lazy (result_s, result_o, resultm_s, resultm_o) ) + let result = + { simulation ; aux_basis ; + hf12 ; hf12_anti ; + hf12_single ; hf12_single_anti + } + in + Parallel.broadcast (lazy result)