mirror of
https://gitlab.com/scemama/QCaml.git
synced 2024-08-17 02:01:44 +02:00
337 lines
8.3 KiB
OCaml
337 lines
8.3 KiB
OCaml
open Lacaml.D
|
|
|
|
type t =
|
|
{
|
|
mo_basis : MOBasis.t ;
|
|
aux_basis : MOBasis.t ;
|
|
det_space : DeterminantSpace.t ;
|
|
ci : CI.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 = DeterminantSpace.mo_class @@ det_space t
|
|
let eigensystem t = Lazy.force t.eigensystem
|
|
|
|
|
|
let f12_integrals mo_basis =
|
|
let two_e_ints = MOBasis.f12_ints mo_basis in
|
|
( (fun _ _ _ -> 0.),
|
|
(fun i j k l s s' ->
|
|
if (i=k && j<>l) || (j=l && i<>k) then
|
|
0.
|
|
else
|
|
begin
|
|
if s' = Spin.other s then
|
|
0.5 *. (F12.get_phys two_e_ints i j k l)
|
|
else
|
|
0.25 *. ((F12.get_phys two_e_ints i j k l) -.
|
|
(F12.get_phys two_e_ints i j l k))
|
|
end
|
|
) )
|
|
|
|
|
|
|
|
|
|
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 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
|
|
|
|
|
|
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 is_internal 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 a ->
|
|
let alfa =
|
|
Determinant.alfa a
|
|
|> Spindeterminant.bitstring
|
|
in
|
|
let beta =
|
|
Determinant.beta a
|
|
|> 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
|
|
| 1 | 2 -> false
|
|
| _ -> true
|
|
|
|
|
|
let dressing_vector ~frozen_core aux_basis f12_amplitudes ci =
|
|
|
|
if Parallel.master then
|
|
Printf.printf "Building matrix\n%!";
|
|
|
|
(* Determinants of the FCI space as a list *)
|
|
let in_dets =
|
|
DeterminantSpace.determinant_stream ci.CI.det_space
|
|
|> Util.stream_to_list
|
|
in
|
|
|
|
(* Stream that generates only singly and doubly excited determinants
|
|
wrt FCI space *)
|
|
let out_dets_stream =
|
|
(* Stream that generates all determinants of FCI space *)
|
|
let s =
|
|
DeterminantSpace.fci_of_mo_basis ~frozen_core aux_basis
|
|
|> DeterminantSpace.determinant_stream
|
|
|
|
in
|
|
(* Select only singly and doubly excited determinants
|
|
wrt FCI space *)
|
|
Stream.from (fun _ ->
|
|
try
|
|
let rec result () =
|
|
let ki = Stream.next s in
|
|
if is_internal ci.CI.det_space ki then
|
|
result ()
|
|
else
|
|
Some ki
|
|
in
|
|
result ()
|
|
with Stream.Failure -> None
|
|
)
|
|
in
|
|
|
|
let make_h_and_f alpha_list =
|
|
|
|
let rec col_vecs_list accu_H accu_F = function
|
|
| [] ->
|
|
List.rev accu_H,
|
|
List.rev accu_F
|
|
| ki :: rest ->
|
|
let h, f =
|
|
List.map (fun kj ->
|
|
match hf_ij aux_basis ki kj with
|
|
| [ a ; b ] -> a, b
|
|
| _ -> assert false ) in_dets
|
|
|> List.split
|
|
in
|
|
let h =
|
|
Vec.of_list h
|
|
and f =
|
|
Vec.of_list f
|
|
in
|
|
col_vecs_list (h::accu_H) (f::accu_F) rest
|
|
in
|
|
let h, f =
|
|
col_vecs_list [] [] alpha_list
|
|
in
|
|
Mat.of_col_vecs_list h,
|
|
Mat.of_col_vecs_list f
|
|
in
|
|
|
|
|
|
let m_HF =
|
|
let batch_size = 1 + 1_000_000 / (Mat.dim1 f12_amplitudes) in
|
|
let input_stream =
|
|
Stream.from (fun i ->
|
|
let rec make_batch accu = function
|
|
| 0 -> accu
|
|
| n -> try
|
|
let alpha = Stream.next out_dets_stream in
|
|
let accu = alpha :: accu in
|
|
make_batch accu (n-1)
|
|
with Stream.Failure -> accu
|
|
in
|
|
let result = make_batch [] batch_size in
|
|
if result = [] then None else Some result
|
|
)
|
|
in
|
|
let result =
|
|
let m_H_aux, m_F_aux = make_h_and_f [(Stream.next out_dets_stream)] in
|
|
let m_HF =
|
|
gemm m_H_aux m_F_aux ~transb:`T
|
|
in
|
|
gemm m_HF f12_amplitudes
|
|
in
|
|
|
|
let iteration input =
|
|
Printf.printf ".%!";
|
|
let m_H_aux, m_F_aux = make_h_and_f input in
|
|
let m_HF =
|
|
gemm m_H_aux m_F_aux ~transb:`T
|
|
in
|
|
gemm m_HF f12_amplitudes
|
|
in
|
|
|
|
input_stream
|
|
|> Farm.run ~ordered:false ~f:iteration
|
|
|> Stream.iter (fun hf ->
|
|
ignore @@ Mat.add result hf ~c:result );
|
|
Printf.printf "\n";
|
|
Parallel.broadcast (lazy result)
|
|
in
|
|
|
|
Printf.printf "Done\n%!";
|
|
Matrix.dense_of_mat m_HF
|
|
|
|
|
|
|
|
|
|
|
|
let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filename () =
|
|
|
|
let f12 = Util.of_some @@ Simulation.f12 simulation in
|
|
let mo_num = MOBasis.size mo_basis in
|
|
|
|
Printf.printf "Add aux basis\n%!";
|
|
(* Add auxiliary basis set *)
|
|
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
|
|
in
|
|
|
|
let aux_basis =
|
|
MOBasis.of_mo_basis s mo_basis
|
|
in
|
|
let () =
|
|
Printf.printf "F12 ints\n%!";
|
|
ignore @@ MOBasis.f12_ints aux_basis
|
|
in
|
|
let () =
|
|
Printf.printf "2e ints\n%!";
|
|
ignore @@ MOBasis.two_e_ints aux_basis
|
|
in
|
|
|
|
Printf.printf "det space\n%!";
|
|
let det_space =
|
|
DeterminantSpace.fci_f12_of_mo_basis aux_basis ~frozen_core mo_num
|
|
in
|
|
|
|
let ci = CI.make det_space in
|
|
|
|
let ci_coef, ci_energy =
|
|
let x = Lazy.force ci.eigensystem in
|
|
Parallel.broadcast (lazy x)
|
|
in
|
|
|
|
let e_shift =
|
|
let det =
|
|
DeterminantSpace.determinant_stream det_space
|
|
|> Stream.next
|
|
in
|
|
h_ij aux_basis det det
|
|
in
|
|
|
|
let eigensystem = lazy (
|
|
let m_H =
|
|
Lazy.force ci.CI.m_H
|
|
in
|
|
|
|
let rec iteration ?(state=1) psi =
|
|
let delta =
|
|
dressing_vector ~frozen_core aux_basis psi ci
|
|
in
|
|
|
|
let f = 1.0 /. psi.{1,1} in
|
|
let delta_00 =
|
|
Util.list_range 2 (Mat.dim1 psi)
|
|
|> List.fold_left (fun accu i -> accu +.
|
|
(Matrix.get delta i 1) *. psi.{i,1} *. f) 0.
|
|
in
|
|
let delta = Matrix.to_mat delta in
|
|
delta.{1,1} <- delta.{1,1} -. delta_00;
|
|
|
|
(*------
|
|
TODO SINGLE STATE HERE
|
|
*)
|
|
let n_states = ci.CI.n_states in
|
|
let diagonal =
|
|
Vec.init (Matrix.dim1 m_H) (fun i ->
|
|
if i = 1 then
|
|
Matrix.get m_H i i +. delta.{1,1} *. f
|
|
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 c11 = Matrix.get c 1 1 in
|
|
Util.list_range 1 (Mat.dim1 w)
|
|
|> List.iter (fun i ->
|
|
let dci =
|
|
delta.{i,1} *. f ;
|
|
in
|
|
w.{i,1} <- w.{i,1} +. dci *. c11;
|
|
if (i <> 1) then
|
|
w.{1,1} <- w.{1,1} +. dci *. (Matrix.get c i 1);
|
|
);
|
|
Matrix.dense_of_mat w
|
|
in
|
|
|
|
let eigenvectors, eigenvalues =
|
|
Parallel.broadcast (lazy (
|
|
Davidson.make ~threshold:1.e-6 ~guess:psi ~n_states diagonal matrix_prod
|
|
))
|
|
in
|
|
|
|
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.{1} +. e_shift
|
|
+. Simulation.nuclear_repulsion simulation);
|
|
|
|
if conv > threshold then
|
|
iteration eigenvectors
|
|
else
|
|
let eigenvalues =
|
|
Vec.map (fun x -> x +. e_shift) eigenvalues
|
|
in
|
|
eigenvectors, eigenvalues
|
|
in
|
|
iteration ci_coef
|
|
|
|
)
|
|
in
|
|
{ mo_basis ; aux_basis ; det_space ; ci ; eigensystem }
|
|
|
|
|