10
1
mirror of https://gitlab.com/scemama/QCaml.git synced 2024-12-22 12:23:31 +01:00
QCaml/Basis/PrimitiveShellPair.ml

165 lines
3.3 KiB
OCaml
Raw Normal View History

2018-03-15 15:25:49 +01:00
open Util
open Constants
type t = {
2018-03-20 15:16:24 +01:00
exponent : float; (* alpha + beta *)
exponent_inv : float; (* 1/(alpha + beta) *)
2018-03-15 15:25:49 +01:00
center : Coordinate.t; (* P = (alpha * A + beta B)/(alpha+beta) *)
center_minus_a : Coordinate.t; (* P - A *)
a_minus_b : Coordinate.t; (* A - B *)
a_minus_b_sq : float; (* |A-B|^2 *)
2018-03-20 15:16:24 +01:00
norm_scales : float array lazy_t;
normalization : float; (* norm_coef_a * norm_coef_b * g, with
2018-03-15 15:25:49 +01:00
g = (pi/(alpha+beta))^(3/2) exp (-|A-B|^2 * alpha*beta/(alpha+beta)) *)
2018-03-21 15:01:39 +01:00
ang_mom : AngularMomentum.t;
2018-03-15 15:25:49 +01:00
shell_a : PrimitiveShell.t;
shell_b : PrimitiveShell.t;
}
exception Null_contribution
module Am = AngularMomentum
module Co = Coordinate
module Ps = PrimitiveShell
2018-03-15 19:11:59 +01:00
2018-03-15 15:25:49 +01:00
let hash a =
Hashtbl.hash a
2018-03-15 19:11:59 +01:00
2018-03-15 15:25:49 +01:00
let equivalent a b =
2018-03-20 15:16:24 +01:00
a.exponent = b.exponent &&
2018-03-21 15:01:39 +01:00
a.ang_mom = b.ang_mom &&
2018-03-20 15:16:24 +01:00
a.normalization = b.normalization &&
2018-03-15 19:11:59 +01:00
a.center = b.center &&
a.center_minus_a = b.center_minus_a &&
a.a_minus_b = b.a_minus_b
2018-03-15 15:25:49 +01:00
let cmp a b =
hash a - hash b
let create_make_of p_a p_b =
let a_minus_b =
Co.( Ps.center p_a |- Ps.center p_b )
in
let a_minus_b_sq =
Co.dot a_minus_b a_minus_b
in
2018-03-20 15:16:24 +01:00
let norm_scales = lazy (
2018-03-15 15:25:49 +01:00
Array.map (fun v1 ->
2018-03-20 15:16:24 +01:00
Array.map (fun v2 -> v1 *. v2) (Ps.norm_scales p_b)
) (Ps.norm_scales p_a)
2018-03-15 15:25:49 +01:00
|> Array.to_list
|> Array.concat
) in
2018-03-21 15:01:39 +01:00
let ang_mom =
Am.( Ps.ang_mom p_a + Ps.ang_mom p_b )
2018-03-15 15:25:49 +01:00
in
function p_a ->
2018-03-22 00:29:14 +01:00
let norm_coef_a =
Ps.normalization p_a
in
2018-03-15 15:25:49 +01:00
2018-03-22 00:29:14 +01:00
let alpha_a =
Co.( Ps.exponent p_a |. Ps.center p_a )
in
2018-03-15 15:25:49 +01:00
2018-03-22 00:29:14 +01:00
function p_b ->
2018-03-15 15:25:49 +01:00
2018-03-22 00:29:14 +01:00
let normalization =
norm_coef_a *. Ps.normalization p_b
in
2018-03-15 15:25:49 +01:00
2018-03-22 00:29:14 +01:00
let exponent =
Ps.exponent p_a +. Ps.exponent p_b
in
2018-03-15 15:25:49 +01:00
2018-03-22 00:29:14 +01:00
let exponent_inv = 1. /. exponent in
2018-03-15 15:25:49 +01:00
2018-03-22 00:29:14 +01:00
let normalization =
let argexpo =
Ps.exponent p_a *. Ps.exponent p_b *. a_minus_b_sq *. exponent_inv
in
normalization *. (pi *. exponent_inv)**1.5 *. exp (-. argexpo)
in
2018-03-15 15:25:49 +01:00
2018-03-22 00:29:14 +01:00
function cutoff ->
2018-03-15 15:25:49 +01:00
2018-03-22 00:29:14 +01:00
if abs_float normalization > cutoff then (
2018-03-15 15:25:49 +01:00
2018-03-22 00:29:14 +01:00
let beta_b =
Co.( Ps.exponent p_b |. Ps.center p_b )
in
2018-03-15 15:25:49 +01:00
2018-03-22 00:29:14 +01:00
let center =
Co.(exponent_inv |. (alpha_a |+ beta_b))
in
2018-03-15 15:25:49 +01:00
2018-03-22 00:29:14 +01:00
let center_minus_a =
Co.(center |- Ps.center p_a)
in
2018-03-15 15:25:49 +01:00
2018-03-22 00:29:14 +01:00
Some {
ang_mom ;
exponent ; exponent_inv ; center ; center_minus_a ; a_minus_b ;
a_minus_b_sq ; normalization ; norm_scales ; shell_a = p_a;
shell_b = p_b }
2018-03-15 15:25:49 +01:00
2018-03-22 00:29:14 +01:00
)
else None
2018-03-15 15:25:49 +01:00
let make p_a p_b =
let f =
create_make_of p_a p_b
in
2018-03-15 19:11:59 +01:00
match f p_a p_b 0. with
2018-03-15 15:25:49 +01:00
| Some result -> result
| None -> assert false
2018-03-20 15:16:24 +01:00
let norm_scales x =
Lazy.force x.norm_scales
2018-03-15 15:25:49 +01:00
2018-03-20 15:16:24 +01:00
let exponent_inv x = x.exponent_inv
2018-03-15 15:25:49 +01:00
let monocentric x =
Ps.center x.shell_a = Ps.center x.shell_b
2018-03-21 15:01:39 +01:00
let ang_mom x = x.ang_mom
2018-03-15 15:25:49 +01:00
let a_minus_b x = x.a_minus_b
let a_minus_b_sq x = x.a_minus_b_sq
let center_minus_a x = x.center_minus_a
2018-03-20 15:16:24 +01:00
let normalization x = x.normalization
2018-03-15 15:25:49 +01:00
2018-03-20 15:16:24 +01:00
let exponent x = x.exponent
2018-03-15 15:25:49 +01:00
let center x = x.center
2018-03-15 16:03:43 +01:00
let shell_a x = x.shell_a
let shell_b x = x.shell_b
2018-03-21 15:01:39 +01:00
let zkey_array x =
Am.zkey_array (Am.Doublet
Ps.(ang_mom x.shell_a, ang_mom x.shell_b)
)