mirror of
https://gitlab.com/scemama/QCaml.git
synced 2024-06-19 19:52:06 +02:00
108 lines
3.5 KiB
OCaml
108 lines
3.5 KiB
OCaml
open Util
|
|
open Constants
|
|
|
|
type t = {
|
|
expo : float; (* alpha + beta *)
|
|
expo_inv : float; (* 1/(alpha + beta) *)
|
|
center : Coordinate.t; (* P = (alpha * A + beta B)/(alpha+beta) *)
|
|
center_a : Coordinate.t; (* P - A *)
|
|
center_ab: Coordinate.t; (* A - B *)
|
|
norm_sq : float; (* |A-B|^2 *)
|
|
norm_coef: float; (* norm_coef_a * norm_coef_b *)
|
|
coef : float; (* norm_coef * coef_a * coef_b * g, with
|
|
g = (pi/(alpha+beta))^(3/2) exp (-|A-B|^2 * alpha*beta/(alpha+beta)) *)
|
|
norm_coef_scale : float array; (* norm_coef.(i) / norm_coef.(0) *)
|
|
i : int;
|
|
j : int;
|
|
shell_a : Contracted_shell.t;
|
|
shell_b : Contracted_shell.t;
|
|
monocentric : bool
|
|
}
|
|
|
|
exception Null_contribution
|
|
|
|
let create_array ?cutoff p_a p_b =
|
|
let cutoff, log_cutoff =
|
|
match cutoff with
|
|
| None -> -1., max_float
|
|
| Some cutoff -> cutoff, -. (log cutoff)
|
|
in
|
|
|
|
let center_ab = Coordinate.(
|
|
Contracted_shell.center p_a |- Contracted_shell.center p_b )
|
|
in
|
|
let norm_sq =
|
|
Coordinate.dot center_ab center_ab
|
|
in
|
|
let norm_coef_scale_a =
|
|
Contracted_shell.norm_coef_scale p_a
|
|
and norm_coef_scale_b =
|
|
Contracted_shell.norm_coef_scale p_b
|
|
in
|
|
let norm_coef_scale =
|
|
Array.map (fun v1 ->
|
|
Array.map (fun v2 -> v1 *. v2) norm_coef_scale_b
|
|
) norm_coef_scale_a
|
|
|> Array.to_list
|
|
|> Array.concat
|
|
in
|
|
Array.init (Contracted_shell.size p_a) (fun i ->
|
|
let p_a_expo_center = Coordinate.(
|
|
Contracted_shell.expo p_a i |. Contracted_shell.center p_a )
|
|
in
|
|
let norm_coef_a =
|
|
Contracted_shell.norm_coef p_a i
|
|
in
|
|
|
|
Array.init (Contracted_shell.size p_b) (fun j ->
|
|
try
|
|
let norm_coef_b =
|
|
Contracted_shell.norm_coef p_b j
|
|
in
|
|
let norm_coef =
|
|
norm_coef_a *. norm_coef_b
|
|
in
|
|
if (norm_coef < cutoff) then
|
|
raise Null_contribution;
|
|
let p_b_expo_center = Coordinate.(
|
|
Contracted_shell.expo p_b j |. Contracted_shell.center p_b )
|
|
in
|
|
let expo = Contracted_shell.(expo p_a i +. expo p_b j) in
|
|
let expo_inv = 1. /. expo in
|
|
let center = Coordinate.(
|
|
expo_inv |. (p_a_expo_center |+ p_b_expo_center ) )
|
|
in
|
|
let argexpo =
|
|
Contracted_shell.(expo p_a i *. expo p_b j) *. norm_sq *. expo_inv
|
|
in
|
|
if (argexpo > log_cutoff) then
|
|
raise Null_contribution;
|
|
let g =
|
|
(pi *. expo_inv)**(1.5) *. exp(-. argexpo)
|
|
in
|
|
let coef =
|
|
norm_coef *. Contracted_shell.(coef p_a i *. coef p_b j) *. g
|
|
in
|
|
if (abs_float coef < cutoff) then
|
|
raise Null_contribution;
|
|
let center_a =
|
|
Coordinate.(center |- Contracted_shell.center p_a)
|
|
in
|
|
let monocentric =
|
|
Contracted_shell.center p_a = Contracted_shell.center p_b
|
|
in
|
|
Some { i ; j ; shell_a=p_a ; shell_b=p_b ; norm_coef ; coef ; expo ; expo_inv ; center ; center_a ; center_ab ; norm_sq ; norm_coef_scale ; monocentric }
|
|
with
|
|
| Null_contribution -> None
|
|
)
|
|
)
|
|
|> Array.to_list
|
|
|> Array.concat
|
|
|> Array.to_list
|
|
|> List.filter (function Some _ -> true | None -> false)
|
|
|> List.map (function Some x -> x | None -> assert false)
|
|
|> Array.of_list
|
|
|
|
open Util
|
|
|