open Lacaml.D module Ds = DeterminantSpace 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 hf_ij_non_zero hf12_integrals deg_a deg_b ki kj = let integrals = [ let one_e _ _ _ = 0. in let hf12_s, hf12_o = hf12_integrals in let two_e i j k l s s' = if s' = Spin.other s then hf12_o.{i,j,k,l} else hf12_s.{i,j,k,l} in (one_e, two_e) ] in CIMatrixElement.non_zero integrals deg_a deg_b ki kj |> List.hd 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 m_HF = let f = match Ds.determinants det_space with | Ds.Arbitrary _ -> CI.create_matrix_arbitrary | Ds.Spin _ -> CI.create_matrix_spin_computed in f (fun deg_a deg_b ki kj -> hf_ij_non_zero hf12_integrals deg_a deg_b ki kj ) det_space in Matrix.mm (Lazy.force m_HF) (Matrix.dense_of_mat f12_amplitudes) 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 = 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 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 = (* Column dressing *) 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 (* 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 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); 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 }