mirror of
https://gitlab.com/scemama/QCaml.git
synced 2024-09-01 17:23:41 +02:00
205 lines
5.0 KiB
OCaml
205 lines
5.0 KiB
OCaml
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 dressing\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
|
|
|
|
let m_HF =
|
|
Lazy.force m_HF
|
|
in
|
|
|
|
let result =
|
|
Matrix.parallel_mm m_HF (Matrix.dense_of_mat f12_amplitudes)
|
|
in
|
|
|
|
if Parallel.master then
|
|
Printf.printf "dressing done\n%!";
|
|
|
|
Parallel.broadcast (lazy result)
|
|
|
|
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 ~frozen_core ~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 c m_H
|
|
|> Matrix.transpose
|
|
|> 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);
|
|
|
|
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 }
|
|
|
|
|
|
|
|
|