10
1
mirror of https://gitlab.com/scemama/QCaml.git synced 2024-06-13 16:55:18 +02:00
QCaml/Basis/Overlap.ml
2018-02-24 23:57:38 +01:00

154 lines
3.7 KiB
OCaml

open Util
open Constants
open Lacaml.D
type t = Mat.t
module Am = AngularMomentum
module Bs = Basis
module Co = Coordinate
module Cs = ContractedShell
module Csp = ContractedShellPair
module Sp = ShellPair
let cutoff = integrals_cutoff
(** Computes all the overlap integrals of the contracted shell pair *)
let contracted_class shell_a shell_b : float Zmap.t =
let shell_p =
Csp.create shell_a shell_b
in
(* Pre-computation of integral class indices *)
let class_indices =
Am.zkey_array (Am.Doublet (shell_a.totAngMom, shell_b.totAngMom))
in
let contracted_class =
Array.make (Array.length class_indices) 0.
in
let sp =
shell_p.Csp.shell_pairs
in
let center_ab =
shell_p.Csp.center_ab
in
let norm_coef_scale =
shell_p.Csp.norm_coef_scale
in
(* Compute all integrals in the shell for each pair of significant shell pairs *)
for ab=0 to (Array.length sp - 1)
do
let coef_prod =
shell_p.Csp.coef.(ab)
in
(** Screening on thr product of coefficients *)
if (abs_float coef_prod) > 1.e-3*.cutoff then
begin
let expo_inv =
shell_p.Csp.expo_inv.(ab)
in
let center_pa =
sp.(ab).Sp.center_a
in
Array.iteri (fun i key ->
let (angMomA,angMomB) =
let a = Zkey.to_int_array Zkey.Kind_6 key in
( [| a.(0) ; a.(1) ; a.(2) |],
[| a.(3) ; a.(4) ; a.(5) |] )
in
let f k =
let xyz = match k with
| 0 -> Co.X
| 1 -> Co.Y
| _ -> Co.Z
in
Overlap_primitives.hvrr (angMomA.(k), angMomB.(k))
expo_inv
(Co.get xyz center_ab,
Co.get xyz center_pa)
in
let norm = norm_coef_scale.(i) in
let integral = chop norm (fun () -> (f 0)*.(f 1)*.(f 2)) in
contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral
) class_indices
end
done;
let result =
Zmap.create (Array.length contracted_class)
in
Array.iteri (fun i key -> Zmap.add result key contracted_class.(i)) class_indices;
result
(** Create overlap matrix *)
let of_basis basis =
let to_powers x =
let open Zkey in
match to_powers Kind_3 x with
| Three x -> x
| _ -> assert false
in
let n = basis.Bs.size
and shell = basis.Bs.contracted_shells
in
let result = Mat.create n n in
for j=0 to (Array.length shell) - 1 do
for i=0 to j do
(* Compute all the integrals of the class *)
let cls =
contracted_class shell.(i) shell.(j)
in
Array.iteri (fun j_c powers_j ->
let j_c = shell.(j).index + j_c + 1 in
let xj = to_powers powers_j in
Array.iteri (fun i_c powers_i ->
let i_c = shell.(i).index + i_c + 1 in
let xi = to_powers powers_i in
let key =
Zkey.of_powers (Zkey.Six (xi,xj))
in
let value =
try Zmap.find cls key
with Not_found -> failwith "Bug in overlap integrals"
in
result.{i_c,j_c} <- value;
result.{j_c,i_c} <- value;
) shell.(i).powers
) shell.(j).powers
done;
done;
Mat.detri result;
result
(** Write all overlap integrals to a file *)
let to_file ~filename overlap =
let oc = open_out filename in
let n =
Mat.dim1 overlap
in
for j=1 to n do
for i=1 to j do
if (abs_float overlap.{i,j} > cutoff) then
Printf.fprintf oc "%4d %4d %20.12e\n" i j overlap.{i,j}
done;
done;
close_out oc