open Lacaml.D module Ds = DeterminantSpace module De = Determinant module Sp = Spindeterminant type t = { mo_basis : MOBasis.t ; det_space : DeterminantSpace.t ; ci : CI.t ; hf12_integrals : HF12.t ; eigensystem : (Mat.t * Vec.t) lazy_t; } let ci t = t.ci let mo_basis t = t.mo_basis let det_space t = t.det_space let mo_class t = Ds.mo_class @@ det_space t let eigensystem t = Lazy.force t.eigensystem let dressing_vector ~frozen_core hf12_integrals f12_amplitudes ci = if Parallel.master then Printf.printf "Building matrix\n%!"; let det_space = ci.CI.det_space in let { HF12. simulation ; aux_basis ; f_0 ; f_1 ; f_2 ; f_3 } = hf12_integrals in let m_HF = let f = match Ds.determinants det_space with | Ds.Arbitrary _ -> CI.create_matrix_arbitrary | Ds.Spin _ -> CI.create_matrix_spin_computed ~nmax:3 in f (fun deg_a deg_b ki kj -> match deg_a + deg_b with | 0 -> f_0 ki | 1 -> f_1 ki kj | 2 -> f_2 ki kj | 3 -> f_3 ki kj | _ -> assert false ) det_space in Matrix.mm (Lazy.force m_HF) (Matrix.dense_of_mat f12_amplitudes) let sum l f = List.fold_left (fun accu i -> accu +. f i) 0. l let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filename ?(state=1) () = let det_space = DeterminantSpace.fci_of_mo_basis mo_basis ~frozen_core in let ci = CI.make ~n_states:state det_space in let hf12_integrals = HF12.make ~simulation ~mo_basis ~aux_basis_filename () in let ci_coef, ci_energy = let x = Lazy.force ci.eigensystem in Parallel.broadcast (lazy x) in let eigensystem = lazy ( let m_H = Lazy.force ci.CI.m_H in let rec iteration ~state psi = (* Format.printf "%a@." DeterminantSpace.pp_det_space @@ CI.det_space ci; Format.printf "%a@." Matrix.pp_matrix @@ Matrix.dense_of_mat psi; *) let column_idx = iamax (Mat.to_col_vecs psi).(state-1) in let delta = (* delta_i = {% $\sum_j c_j H_{ij}$ %} *) dressing_vector ~frozen_core hf12_integrals psi ci |> Matrix.to_mat in (* Format.printf "%a@." Matrix.pp_matrix @@ Matrix.dense_of_mat delta; *) Printf.printf "Cmax : %e\n" psi.{column_idx,state}; Printf.printf "Norm : %e\n" (sqrt (gemm ~transa:`T delta delta).{state,state}); let f = 1.0 /. psi.{column_idx,state} in let delta_00 = (* Delta_00 = {% $\sum_{j \ne x} delta_j c_j / c_x$ %} *) f *. ( (gemm ~transa:`T delta psi).{state,state} -. delta.{column_idx,state} *. psi.{column_idx,state} ) in Printf.printf "Delta_00 : %e %e\n" delta.{column_idx,state} delta_00; delta.{column_idx,state} <- delta.{column_idx,state} -. delta_00; let eigenvectors, eigenvalues = let delta = lacpy delta in Mat.scal f delta; for k=1 to state-1 do for i=1 to Mat.dim1 delta do delta.{i,k} <- delta.{i,state} done; done; let diagonal = Vec.init (Matrix.dim1 m_H) (fun i -> if i = column_idx then Matrix.get m_H i i +. delta.{column_idx,state} 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 c = Matrix.to_mat c in for k=1 to state do for i=1 to (Mat.dim1 w) do w.{i,k} <- w.{i,k} +. delta.{i,k} *. c.{column_idx, k} ; w.{column_idx,k} <- w.{column_idx,k} +. delta.{i,k} *. c.{i,k}; done; w.{column_idx,k} <- w.{column_idx,k} -. delta.{column_idx,k} *. c.{column_idx,k}; done; Matrix.dense_of_mat w in Parallel.broadcast (lazy ( Davidson.make ~threshold:1.e-10 ~guess:psi ~n_states:state diagonal matrix_prod )) in let eigenvectors = Conventions.rephase @@ Util.remove_epsilons eigenvectors in Vec.iter (fun energy -> Printf.printf "%g\t" energy) eigenvalues; print_newline (); let conv = 1.0 -. abs_float ( dot (Mat.to_col_vecs psi).(0) (Mat.to_col_vecs eigenvectors).(0) ) in if Parallel.master then 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 let eigenvalues = Vec.map (fun x -> x +. ci.CI.e_shift) eigenvalues in eigenvectors, eigenvalues in iteration ~state ci_coef ) in { mo_basis ; det_space ; ci ; hf12_integrals ; eigensystem }