LCPQ
/
QCaml
mirror of https://gitlab.com/scemama/QCaml.git
10
1
Fork 0
Code Releases Activity

A lot of cleaning

This commit is contained in:
Anthony Scemama 2018-02-13 17:36:25 +01:00
parent a892f2f125
commit ade45cd49d
22 changed files with 450 additions and 210 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
Basis/Contracted_shell.ml
View File

@ -1,6 +1,7 @@
open Util
open Constants
open Contracted_shell_type
open Coordinate
type t = Contracted_shell_type.t
@ -16,16 +17,14 @@ let with_index = Contracted_shell_type.with_index
let powers a = a.powers
let to_string s =
let coord =
Coordinate.to_Bohr s.center
in
let coord = s.center in
let open Printf in
(match s.totAngMom with
| Angular_momentum.S -> sprintf "%3d " (s.index+1)
| _ -> sprintf "%3d-%-3d" (s.index+1) (s.index+(Array.length s.powers))
) ^
( sprintf "%1s %8.3f %8.3f %8.3f " (Angular_momentum.to_string s.totAngMom)
(Coordinate.x coord) (Coordinate.y coord) (Coordinate.z coord) ) ^
(get X coord) (get Y coord) (get Z coord) ) ^
(Array.map2 (fun e c -> sprintf "%16.8e %16.8e" e c) s.expo s.coef
|> Array.to_list |> String.concat (sprintf "\n%36s" " ") )

73
Basis/Contracted_shell_type.ml Normal file
View File

@ -0,0 +1,73 @@
open Util
open Constants
type shell_contracted = {
expo : float array;
coef : float array;
center : Coordinate.t;
totAngMom : Angular_momentum.t;
size : int;
norm_coef : float array;
norm_coef_scale : float array;
index : int;
powers : Zkey.t array;
}
type t = shell_contracted
(** Normalization coefficient of contracted function i, which depends on the
exponent and the angular momentum. Two conventions can be chosen : a single
normalisation factor for all functions of the class, or a coefficient which
depends on the powers of x,y and z.
Returns, for each contracted function, an array of functions taking as
argument the [|x;y;z|] powers.
*)
let compute_norm_coef expo totAngMom =
let atot =
Angular_momentum.to_int totAngMom
in
let factor int_array =
let dfa = Array.map (fun j ->
( float_of_int (1 lsl j) *. fact j) /. fact (j+j)
) int_array
in
sqrt (dfa.(0) *.dfa.(1) *. dfa.(2))
in
let expo =
if atot mod 2 = 0 then
Array.map (fun alpha ->
let alpha_2 = alpha +. alpha in
(alpha_2 *. pi_inv)**(0.75) *. (pow (alpha_2 +. alpha_2) (atot/2))
) expo
else
Array.map (fun alpha ->
let alpha_2 = alpha +. alpha in
(alpha_2 *. pi_inv)**(0.75) *. sqrt (pow (alpha_2 +. alpha_2) atot)
) expo
in
Array.map (fun x -> let f a = x *. (factor a) in f) expo
let make ~index ~expo ~coef ~center ~totAngMom =
assert (Array.length expo = Array.length coef);
assert (Array.length expo > 0);
let norm_coef_func =
compute_norm_coef expo totAngMom
in
let powers =
Angular_momentum.zkey_array (Angular_momentum.Singlet totAngMom)
in
let norm_coef =
Array.map (fun f -> f [| Angular_momentum.to_int totAngMom ; 0 ; 0 |]) norm_coef_func
in
let norm_coef_scale =
Array.map (fun a ->
(norm_coef_func.(0) (Zkey.to_int_array ~kind:Zkey.Kind_3 a)) /. norm_coef.(0)
) powers
in
{ index ; expo ; coef ; center ; totAngMom ; size=Array.length expo ; norm_coef ;
norm_coef_scale ; powers }
let with_index a i = { a with index = i }

24
Basis/Contracted_shell_type.mli Normal file
View File

@ -0,0 +1,24 @@
type shell_contracted = private {
expo : float array;
coef : float array;
center : Coordinate.t;
totAngMom : Angular_momentum.t;
size : int;
norm_coef : float array;
norm_coef_scale : float array;
index : int;
powers : Zkey.t array;
}
type t = shell_contracted
val make :
index:int ->
expo:float array ->
coef:float array ->
center:Coordinate.t -> totAngMom:Angular_momentum.t -> t
(** Returns a copy of the contracted shell with a modified index *)
val with_index : t -> int -> t

10
Basis/KinInt.ml
View File

@ -1,6 +1,7 @@
open Util
open Constants
open Lacaml.D
open Coordinate
type t = Mat.t
@ -62,10 +63,15 @@ let contracted_class shell_a shell_b : float Zmap.t =
[| a.(3) ; a.(4) ; a.(5) |] )
in
let ov a b k =
let xyz = match k with
| 0 -> X
| 1 -> Y
| _ -> Z
in
Overlap_primitives.hvrr (a, b)
expo_inv
(Coordinate.coord center_ab k,
Coordinate.coord center_pa k)
(Coordinate.get xyz center_ab,
Coordinate.get xyz center_pa)
in
let f k =
ov angMomA.(k) angMomB.(k) k

31
Basis/OneElectronRR.ml
View File

@ -1,5 +1,6 @@
open Util
open Constants
open Coordinate
exception NullPair
@ -27,7 +28,7 @@ let hvrr_one_e
(** Vertical recurrence relations *)
let rec vrr angMom_a totAngMom_a =
let ax,ay,az = angMom_a in
if (ax < 0) || (ay < 0) || (az < 0) then
if ax < 0 || ay < 0 || az < 0 then
empty
else
match totAngMom_a with
@ -40,13 +41,13 @@ let hvrr_one_e
let result =
let am, amm, amxyz, xyz =
match angMom_a with
| (x,0,0) -> (x-1,0,0),(x-2,0,0), x-1, 0
| (x,y,0) -> (x,y-1,0),(x,y-2,0), y-1, 1
| (x,y,z) -> (x,y,z-1),(x,y,z-2), z-1, 2
| (x,0,0) -> (x-1,0,0),(x-2,0,0), x-1, X
| (x,y,0) -> (x,y-1,0),(x,y-2,0), y-1, Y
| (x,y,z) -> (x,y,z-1),(x,y,z-2), z-1, Z
in
if amxyz < 0 then empty else
let f1 = Coordinate.coord center_pa xyz
and f2 = expo_inv_p *. (Coordinate.coord center_pc xyz)
let f1 = Coordinate.get xyz center_pa
and f2 = expo_inv_p *. (Coordinate.get xyz center_pc)
in
if amxyz < 1 then
let v1 =
@ -77,7 +78,7 @@ let hvrr_one_e
and hrr angMom_a angMom_b totAngMom_a totAngMom_b =
let bx,by,bz = angMom_b in
if (bx < 0) || (by < 0) || (bz < 0) then 0.
if bx < 0 || by < 0 || bz < 0 then 0.
else
match totAngMom_b with
| 0 -> (vrr angMom_a totAngMom_a).(0)
@ -86,24 +87,24 @@ let hvrr_one_e
and angMom_bx, angMom_by, angMom_bz = angMom_b in
let bxyz, xyz =
match angMom_b with
| (_,0,0) -> angMom_bx, 0
| (_,_,0) -> angMom_by, 1
| (_,_,_) -> angMom_bz, 2
| (_,0,0) -> angMom_bx, X
| (_,_,0) -> angMom_by, Y
| (_,_,_) -> angMom_bz, Z
in
if (bxyz < 1) then 0. else
let ap, bm =
match xyz with
| 0 -> (angMom_ax+1,angMom_ay,angMom_az),(angMom_bx-1,angMom_by,angMom_bz)
| 1 -> (angMom_ax,angMom_ay+1,angMom_az),(angMom_bx,angMom_by-1,angMom_bz)
| _ -> (angMom_ax,angMom_ay,angMom_az+1),(angMom_bx,angMom_by,angMom_bz-1)
| X -> (angMom_ax+1,angMom_ay,angMom_az),(angMom_bx-1,angMom_by,angMom_bz)
| Y -> (angMom_ax,angMom_ay+1,angMom_az),(angMom_bx,angMom_by-1,angMom_bz)
| Z -> (angMom_ax,angMom_ay,angMom_az+1),(angMom_bx,angMom_by,angMom_bz-1)
in
let h1 =
hrr ap bm (totAngMom_a+1) (totAngMom_b-1)
in
let f2 =
(Coordinate.coord center_ab xyz)
Coordinate.get xyz center_ab
in
if (abs_float f2 < cutoff) then h1 else
if abs_float f2 < cutoff then h1 else
let h2 =
hrr angMom_a bm totAngMom_a (totAngMom_b-1)
in

10
Basis/Overlap.ml
View File

@ -1,6 +1,7 @@
open Util
open Constants
open Lacaml.D
open Coordinate
type t = Mat.t
@ -56,10 +57,15 @@ let contracted_class shell_a shell_b : float Zmap.t =
[| a.(3) ; a.(4) ; a.(5) |] )
in
let f k =
let xyz = match k with
| 0 -> X
| 1 -> Y
| _ -> Z
in
Overlap_primitives.hvrr (angMomA.(k), angMomB.(k))
expo_inv
(Coordinate.coord center_ab k,
Coordinate.coord center_pa k)
(Coordinate.get xyz center_ab,
Coordinate.get xyz center_pa)
in
let norm = norm_coef_scale.(i) in
let integral = chop norm (fun () -> (f 0)*.(f 1)*.(f 2)) in

76
Basis/TwoElectronRR.ml
View File

@ -1,5 +1,6 @@
open Util
open Constants
open Coordinate
let debug=false
@ -45,9 +46,9 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
let (x,y,z) = angMom_d in Printf.printf "%d %d %d\n" x y z;
Printf.printf "%f %f %f %f\n%f %f %f\n%f %f %f\n%f %f %f\n" expo_b expo_d
expo_inv_p expo_inv_q
(Coordinate.coord center_ab 0) (Coordinate.coord center_ab 1) (Coordinate.coord center_ab 2)
(Coordinate.coord center_cd 0) (Coordinate.coord center_cd 1) (Coordinate.coord center_cd 2)
(Coordinate.coord center_pq 0) (Coordinate.coord center_pq 1) (Coordinate.coord center_pq 2)
(get X center_ab) (get Y center_ab) (get Z center_ab)
(get X center_cd) (get Y center_cd) (get Z center_cd)
(get X center_pq) (get Y center_pq) (get Z center_pq)
end;
(** Vertical recurrence relations *)
@ -68,13 +69,13 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
let result =
let am, amm, amxyz, xyz =
match angMom_a with
| (x,0,0) -> (x-1,0,0),(x-2,0,0), x-1, 0
| (x,y,0) -> (x,y-1,0),(x,y-2,0), y-1, 1
| (x,y,z) -> (x,y,z-1),(x,y,z-2), z-1, 2
| (x,0,0) -> (x-1,0,0),(x-2,0,0), x-1, X
| (x,y,0) -> (x,y-1,0),(x,y-2,0), y-1, Y
| (x,y,z) -> (x,y,z-1),(x,y,z-2), z-1, Z
in
if amxyz < 0 then empty else
let f1 = expo_inv_p *. (Coordinate.coord center_pq xyz)
and f2 = expo_b *. expo_inv_p *. (Coordinate.coord center_ab xyz)
let f1 = expo_inv_p *. (Coordinate.get xyz center_pq)
and f2 = expo_b *. expo_inv_p *. (Coordinate.get xyz center_ab)
in
let result = Array.create_float maxsze in
if amxyz < 1 then
@ -140,24 +141,24 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
(angMom_ax-1, angMom_ay, angMom_az),
(angMom_cx-1, angMom_cy, angMom_cz),
(angMom_cx-2, angMom_cy, angMom_cz),
angMom_ax,angMom_cx-1, 0
angMom_ax,angMom_cx-1, X
| (_,_,0) ->
(angMom_ax, angMom_ay-1, angMom_az),
(angMom_cx, angMom_cy-1, angMom_cz),
(angMom_cx, angMom_cy-2, angMom_cz),
angMom_ay,angMom_cy-1, 1
angMom_ay,angMom_cy-1, Y
| _ ->
(angMom_ax, angMom_ay, angMom_az-1),
(angMom_cx, angMom_cy, angMom_cz-1),
(angMom_cx, angMom_cy, angMom_cz-2),
angMom_az,angMom_cz-1, 2
angMom_az,angMom_cz-1, Z
in
if cmxyz < 0 then empty
else
let f1 =
-. expo_d *. expo_inv_q *. (Coordinate.coord center_cd xyz)
-. expo_d *. expo_inv_q *. (Coordinate.get xyz center_cd)
and f2 =
expo_inv_q *. (Coordinate.coord center_pq xyz)
expo_inv_q *. (Coordinate.get xyz center_pq)
in
let result = Array.make maxsze 0. in
if ( (abs_float f1 > cutoff) || (abs_float f2 > cutoff) ) then
@ -213,6 +214,7 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
and hrr0 angMom_a angMom_b angMom_c
totAngMom_a totAngMom_b totAngMom_c =
(*
if debug then
begin
let angMom_ax, angMom_ay, angMom_az = angMom_a
@ -224,6 +226,7 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
angMom_bx angMom_by angMom_bz
angMom_cx angMom_cy angMom_cz
end;
*)
match totAngMom_b with
| 0 -> (vrr angMom_a angMom_c totAngMom_a totAngMom_c).(0)
@ -231,15 +234,15 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
let angMom_ax, angMom_ay, angMom_az = angMom_a in
let ap, xyz =
match angMom_b with
| (1,_,_) -> (angMom_ax+1,angMom_ay,angMom_az), 0
| (_,1,_) -> (angMom_ax,angMom_ay+1,angMom_az), 1
| _ -> (angMom_ax,angMom_ay,angMom_az+1), 2
| (1,_,_) -> (angMom_ax+1,angMom_ay,angMom_az), X
| (_,1,_) -> (angMom_ax,angMom_ay+1,angMom_az), Y
| _ -> (angMom_ax,angMom_ay,angMom_az+1), Z
in
let v1 =
vrr ap angMom_c (totAngMom_a+1) totAngMom_c
in
let f2 =
(Coordinate.coord center_ab xyz)
(Coordinate.get xyz center_ab)
in
if (abs_float f2 < cutoff) then v1.(0) else
let v2 =
@ -251,22 +254,22 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
and angMom_bx, angMom_by, angMom_bz = angMom_b in
let bxyz, xyz =
match angMom_b with
| (_,0,0) -> angMom_bx, 0
| (_,_,0) -> angMom_by, 1
| (_,_,_) -> angMom_bz, 2
| (_,0,0) -> angMom_bx, X
| (_,_,0) -> angMom_by, Y
| (_,_,_) -> angMom_bz, Z
in
if (bxyz < 1) then 0. else
let ap, bm =
match xyz with
| 0 -> (angMom_ax+1,angMom_ay,angMom_az),(angMom_bx-1,angMom_by,angMom_bz)
| 1 -> (angMom_ax,angMom_ay+1,angMom_az),(angMom_bx,angMom_by-1,angMom_bz)
| _ -> (angMom_ax,angMom_ay,angMom_az+1),(angMom_bx,angMom_by,angMom_bz-1)
| X -> (angMom_ax+1,angMom_ay,angMom_az),(angMom_bx-1,angMom_by,angMom_bz)
| Y -> (angMom_ax,angMom_ay+1,angMom_az),(angMom_bx,angMom_by-1,angMom_bz)
| Z -> (angMom_ax,angMom_ay,angMom_az+1),(angMom_bx,angMom_by,angMom_bz-1)
in
let h1 =
hrr0 ap bm angMom_c (totAngMom_a+1) (totAngMom_b-1) totAngMom_c
in
let f2 =
(Coordinate.coord center_ab xyz)
(Coordinate.get xyz center_ab)
in
if (abs_float f2 < cutoff) then h1 else
let h2 =
@ -274,23 +277,10 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
in
h1 +. f2 *. h2
and hrr angMom_a angMom_b angMom_c angMom_d
totAngMom_a totAngMom_b totAngMom_c totAngMom_d =
if debug then
begin
let angMom_ax, angMom_ay, angMom_az = angMom_a in
let angMom_bx, angMom_by, angMom_bz = angMom_b in
let angMom_cx, angMom_cy, angMom_cz = angMom_c in
let angMom_dx, angMom_dy, angMom_dz = angMom_d in
Printf.printf "hrr : %d %d %d %d : %d %d %d %d %d %d %d %d %d %d %d %d\n"
totAngMom_a totAngMom_b totAngMom_c totAngMom_d
angMom_ax angMom_ay angMom_az
angMom_bx angMom_by angMom_bz
angMom_cx angMom_cy angMom_cz
angMom_dx angMom_dy angMom_dz
end;
match (totAngMom_b, totAngMom_d) with
| (_,0) -> if (totAngMom_b = 0) then
(vrr angMom_a angMom_c totAngMom_a totAngMom_c).(0)
@ -301,14 +291,14 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
and (angMom_dx, angMom_dy, angMom_dz) = angMom_d in
let cp, dm, xyz =
match angMom_d with
| (_,0,0) -> (angMom_cx+1, angMom_cy, angMom_cz), (angMom_dx-1, angMom_dy, angMom_dz), 0
| (_,_,0) -> (angMom_cx, angMom_cy+1, angMom_cz), (angMom_dx, angMom_dy-1, angMom_dz), 1
| _ -> (angMom_cx, angMom_cy, angMom_cz+1), (angMom_dx, angMom_dy, angMom_dz-1), 2
| (_,0,0) -> (angMom_cx+1, angMom_cy, angMom_cz), (angMom_dx-1, angMom_dy, angMom_dz), X
| (_,_,0) -> (angMom_cx, angMom_cy+1, angMom_cz), (angMom_dx, angMom_dy-1, angMom_dz), Y
| _ -> (angMom_cx, angMom_cy, angMom_cz+1), (angMom_dx, angMom_dy, angMom_dz-1), Z
in
let h1 =
hrr angMom_a angMom_b cp dm totAngMom_a totAngMom_b (totAngMom_c+1) (totAngMom_d-1)
in
let f2 = Coordinate.coord center_cd xyz in
let f2 = Coordinate.get xyz center_cd in
if (abs_float f2 < cutoff) then h1 else
let h2 =
hrr angMom_a angMom_b angMom_c dm totAngMom_a totAngMom_b totAngMom_c (totAngMom_d-1)
@ -372,7 +362,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
raise NullQuartet;
let center_pq =
Coordinate.(sp.(ab).ShellPair.center |- sq.(cd).ShellPair.center)
sp.(ab).ShellPair.center |- sq.(cd).ShellPair.center
in
let norm_pq_sq =
Coordinate.dot center_pq center_pq

82
Basis/TwoElectronRRVectorized.ml
View File

@ -1,6 +1,7 @@
open Util
open Lacaml.D
open Bigarray
open Coordinate
let cutoff = Constants.cutoff
let cutoff2 = cutoff *. cutoff
@ -42,15 +43,15 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
let result =
let am, amm, amxyz, xyz =
match angMom_a with
| (x,0,0) -> (x-1,0,0),(x-2,0,0), x-1, 0
| (x,y,0) -> (x,y-1,0),(x,y-2,0), y-1, 1
| (x,y,z) -> (x,y,z-1),(x,y,z-2), z-1, 2
| (x,0,0) -> (x-1,0,0),(x-2,0,0), x-1, X
| (x,y,0) -> (x,y-1,0),(x,y-2,0), y-1, Y
| (x,y,z) -> (x,y,z-1),(x,y,z-2), z-1, Z
in
if amxyz < 0 then
None
else
begin
let cab = Coordinate.coord center_ab xyz in
let cab = Coordinate.get xyz center_ab in
let v1_top, p1_top =
if abs_float cab < cutoff then
None,
@ -89,7 +90,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
for l=0 to np-1 do
for k=0 to nq-1 do
result.(l).(k) <- result.(l).(k)
+. expo_inv_p.(l) *. center_pq.(xyz).(l).(k) *. p0.(l).(k)
+. expo_inv_p.(l) *. (center_pq xyz).(l).(k) *. p0.(l).(k)
done
done
end
@ -126,7 +127,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
for k=0 to nq-1
do
result.(l).(k) <- result.(l).(k) +.
expo_inv_p.(l) *. center_pq.(xyz).(l).(k) *. p0.(l).(k) +.
expo_inv_p.(l) *. (center_pq xyz).(l).(k) *. p0.(l).(k) +.
f *. (v1.(l).(k) +. v2.(l).(k) *. expo_inv_p.(l))
done
done
@ -157,24 +158,24 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
(aax, aay, aaz-1),
(acx, acy, acz-1),
(acx, acy, acz-2),
aaz, acz, 2
aaz, acz, Z
else if (acy > 0) then
(aax, aay-1,aaz),
(acx, acy-1,acz),
(acx, acy-2,acz),
aay,acy, 1
aay,acy, Y
else
(aax-1,aay,aaz),
(acx-1,acy,acz),
(acx-2,acy,acz),
aax,acx, 0
aax,acx, X
in
(*
let result = Array.make_matrix np nq 0. in
*)
let do_compute = ref false in
let v1 =
let f = -. (Coordinate.coord center_cd xyz) in
let f = -. (Coordinate.get xyz center_cd) in
let f1 =
Array.init nq (fun k ->
@ -203,7 +204,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
let f2 =
Array.init np (fun l ->
Array.init nq (fun k ->
let x = expo_inv_q.(k) *. center_pq.(xyz).(l).(k) in
let x = expo_inv_q.(k) *. (center_pq xyz).(l).(k) in
if (!do_compute) then x
else (if abs_float x > cutoff then do_compute := true ; x)
) )
@ -382,11 +383,11 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
let (aax, aay, aaz) = angMom_a in
let ap, xyz =
match angMom_b with
| (_,_,1) -> (aax,aay,aaz+1), 2
| (_,1,_) -> (aax,aay+1,aaz), 1
| (_,_,_) -> (aax+1,aay,aaz), 0
| (_,_,1) -> (aax,aay,aaz+1), Z
| (_,1,_) -> (aax,aay+1,aaz), Y
| (_,_,_) -> (aax+1,aay,aaz), X
in
let f = Coordinate.coord center_ab xyz in
let f = Coordinate.get xyz center_ab in
let v1 =
match vrr_v 0 ap angMom_c (totAngMom_a+1) totAngMom_c with
| Some matrix -> Array.fold_left (fun accu c -> accu +. Array.fold_left (+.) 0. c) 0. matrix
@ -404,23 +405,23 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
and (abx, aby, abz) = angMom_b in
let bxyz, xyz =
match angMom_b with
| (0,0,_) -> abz, 2
| (0,_,_) -> aby, 1
| _ -> abx, 0
| (0,0,_) -> abz, Z
| (0,_,_) -> aby, Y
| _ -> abx, X
in
if (bxyz < 1) then 0. else
let ap, bm =
match xyz with
| 0 -> (aax+1,aay,aaz),(abx-1,aby,abz)
| 1 -> (aax,aay+1,aaz),(abx,aby-1,abz)
| _ -> (aax,aay,aaz+1),(abx,aby,abz-1)
| X -> (aax+1,aay,aaz),(abx-1,aby,abz)
| Y -> (aax,aay+1,aaz),(abx,aby-1,abz)
| Z -> (aax,aay,aaz+1),(abx,aby,abz-1)
in
let h1 =
hrr0_v ap bm angMom_c (totAngMom_a+1) (totAngMom_b-1) totAngMom_c
in
let f = (Coordinate.coord center_ab xyz) in
if (abs_float f < cutoff) then h1 else
let f = Coordinate.get xyz center_ab in
if abs_float f < cutoff then h1 else
let h2 =
hrr0_v angMom_a bm angMom_c totAngMom_a (totAngMom_b-1) totAngMom_c
in
@ -430,7 +431,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
totAngMom_a totAngMom_b totAngMom_c totAngMom_d =
match (totAngMom_b, totAngMom_d) with
| (_,0) -> if (totAngMom_b = 0) then
| (_,0) -> if totAngMom_b = 0 then
begin
match vrr_v 0 angMom_a angMom_c totAngMom_a totAngMom_c with
| Some matrix -> Array.fold_left (fun accu c -> accu +. Array.fold_left (+.) 0. c) 0. matrix
@ -443,15 +444,15 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
and (adx, ady, adz) = angMom_d in
let cp, dm, xyz =
match angMom_d with
| (_,0,0) -> (acx+1, acy, acz), (adx-1, ady, adz), 0
| (_,_,0) -> (acx, acy+1, acz), (adx, ady-1, adz), 1
| _ -> (acx, acy, acz+1), (adx, ady, adz-1), 2
| (_,0,0) -> (acx+1, acy, acz), (adx-1, ady, adz), X
| (_,_,0) -> (acx, acy+1, acz), (adx, ady-1, adz), Y
| _ -> (acx, acy, acz+1), (adx, ady, adz-1), Z
in
let h1 =
hrr_v angMom_a angMom_b cp dm totAngMom_a totAngMom_b (totAngMom_c+1) (totAngMom_d-1)
in
let f = (Coordinate.coord center_cd xyz) in
if (abs_float f < cutoff) then
let f = Coordinate.get xyz center_cd in
if abs_float f < cutoff then
h1
else
let h2 =
@ -570,7 +571,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
in
let center_pq =
Coordinate.(sp.(i-1).ShellPair.center |- sq.(j-1).ShellPair.center)
sp.(i-1).ShellPair.center |- sq.(j-1).ShellPair.center
in
let norm_pq_sq =
Coordinate.dot center_pq center_pq
@ -608,6 +609,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
let norm_coef_scale_p = shell_p.ContractedShellPair.norm_coef_scale in
let center_pq =
let result =
Array.init 3 (fun xyz ->
Array.init np (fun ab ->
let shell_ab = sp.(ab) in
@ -615,16 +617,20 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
let shell_cd = sq.(cd)
in
let cpq =
Coordinate.(shell_ab.ShellPair.center |- shell_cd.ShellPair.center)
shell_ab.ShellPair.center |- shell_cd.ShellPair.center
in
match xyz with
| 0 -> Coordinate.x cpq;
| 1 -> Coordinate.y cpq;
| 2 -> Coordinate.z cpq;
| 0 -> Coordinate.get X cpq;
| 1 -> Coordinate.get Y cpq;
| 2 -> Coordinate.get Z cpq;
| _ -> assert false
)
)
)
in function
| X -> result.(0)
| Y -> result.(1)
| Z -> result.(2)
in
let zero_m_array =
let result =
@ -642,9 +648,11 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
expo_inv_p.(ab) +. expo_inv_q.(cd)
in
let norm_pq_sq =
center_pq.(0).(ab).(cd) *. center_pq.(0).(ab).(cd) +.
center_pq.(1).(ab).(cd) *. center_pq.(1).(ab).(cd) +.
center_pq.(2).(ab).(cd) *. center_pq.(2).(ab).(cd)
let x = (center_pq X).(ab).(cd)
and y = (center_pq Y).(ab).(cd)
and z = (center_pq Z).(ab).(cd)
in
x *. x +. y *. y +. z *. z
in
zero_m ~maxm ~expo_pq_inv ~norm_pq_sq

24
Nuclei/Nuclei.ml
View File

@ -1,3 +1,5 @@
open Xyz_ast
type t = (Element.t * Coordinate.t) array
let of_xyz_file ~filename =
@ -5,7 +7,21 @@ let of_xyz_file ~filename =
let ic = open_in filename in
Lexing.from_channel ic
in
Xyz_parser.input Nuclei_lexer.read_all lexbuf
let data =
Xyz_parser.input Nuclei_lexer.read_all lexbuf
in
let len = List.length data.nuclei in
if len <> data.number_of_atoms then
Printf.sprintf "Error: expected %d atoms but %d read"
data.number_of_atoms len
|> failwith;
List.map (fun nucleus ->
nucleus.element, Coordinate.angstrom_to_bohr nucleus.coord
) data.nuclei
|> Array.of_list
let of_zmt_file ~filename =
@ -21,7 +37,7 @@ let of_zmt_file ~filename =
in aux []
|> Zmatrix.of_string
|> Zmatrix.to_xyz
|> Array.map (fun (e,x,y,z) -> (e, Coordinate.of_3_floats x y z ))
|> Array.map (fun (e,x,y,z) -> (e, (Angstrom.make {Point.x ; y ; z} )))
let to_string atoms =
@ -36,11 +52,11 @@ let to_string atoms =
" ^
(Array.mapi (fun i (e, coord) ->
let coord =
Coordinate.to_Angstrom coord
Coordinate.bohr_to_angstrom coord
in
Printf.sprintf " %5d %5d %5s %12.6f %12.6f %12.6f"
(i+1) (Element.to_int e) (Element.to_string e)
(Coordinate.x coord) (Coordinate.y coord) (Coordinate.z coord)
coord.Angstrom.x coord.Angstrom.y coord.Angstrom.z
) atoms
|> Array.to_list
|> String.concat "\n" ) ^

13
Nuclei/Xyz_ast.mli Normal file
View File

@ -0,0 +1,13 @@
type nucleus =
{
element: Element.t ;
coord: Angstrom.t;
}
type xyz_file =
{
number_of_atoms : int;
file_title : string;
nuclei : nucleus list;
}

59
Nuclei/Xyz_parser.mly
View File

@ -2,6 +2,24 @@
%{
exception InputError of string
let make_angstrom x y z =
Angstrom.make {
Point.
x ; y ; z
}
let output_of f x y z =
let a = make_angstrom x y z in
fun e ->
{
Xyz_ast.
element = f e;
coord = a ;
}
let output_of_string = output_of Element.of_string
let output_of_int = output_of Element.of_int
%}
%token EOL
@ -12,54 +30,59 @@ exception InputError of string
%token EOF
%start input
%type <(Element.t * Coordinate.t) array> input
%type <Xyz_ast.xyz_file> input
%% /* Grammar rules and actions follow */
input:
| integer title atoms_xyz {
let len = List.length $3 in
if len <> $1 then
let error_msg = Printf.sprintf "%d atoms entered, expected %d" len $1 in
raise (InputError error_msg)
else
Array.of_list $3
}
{
number_of_atoms = $1;
file_title = $2;
nuclei = $3;
}
}
;
integer:
| INTEGER EOL { $1 }
| INTEGER SPACE EOL { $1 }
| SPACE INTEGER EOL { $2 }
| SPACE INTEGER SPACE EOL { $2 }
;
title:
| title_list EOL { $1 }
;
text:
| WORD { $1 }
| SPACE { $1 }
| FLOAT { (string_of_float $1)}
| INTEGER { (string_of_int $1)}
;
title_list:
| { "" }
| title_list text { ($1 ^ $2) }
;
atoms_xyz:
| atoms_list EOL { List.rev $1 }
| atoms_list EOF { List.rev $1 }
;
atoms_list:
| { [] }
| atoms_list WORD SPACE FLOAT SPACE FLOAT SPACE FLOAT EOL { (Element.of_string $2, Coordinate.of_3_floats $4 $6 $8 `Angstrom ) :: $1 }
| atoms_list WORD SPACE FLOAT SPACE FLOAT SPACE FLOAT SPACE EOL { (Element.of_string $2, Coordinate.of_3_floats $4 $6 $8 `Angstrom ) :: $1 }
| atoms_list INTEGER SPACE FLOAT SPACE FLOAT SPACE FLOAT EOL { (Element.of_int $2, Coordinate.of_3_floats $4 $6 $8 `Angstrom ) :: $1 }
| atoms_list INTEGER SPACE FLOAT SPACE FLOAT SPACE FLOAT SPACE EOL { (Element.of_int $2, Coordinate.of_3_floats $4 $6 $8 `Angstrom ) :: $1 }
| atoms_list SPACE WORD SPACE FLOAT SPACE FLOAT SPACE FLOAT EOL { (Element.of_string $3, Coordinate.of_3_floats $5 $7 $9 `Angstrom ) :: $1 }
| atoms_list SPACE WORD SPACE FLOAT SPACE FLOAT SPACE FLOAT SPACE EOL { (Element.of_string $3, Coordinate.of_3_floats $5 $7 $9 `Angstrom ) :: $1 }
| atoms_list SPACE INTEGER SPACE FLOAT SPACE FLOAT SPACE FLOAT EOL { (Element.of_int $3, Coordinate.of_3_floats $5 $7 $9 `Angstrom ) :: $1 }
| atoms_list SPACE INTEGER SPACE FLOAT SPACE FLOAT SPACE FLOAT SPACE EOL { (Element.of_int $3, Coordinate.of_3_floats $5 $7 $9 `Angstrom ) :: $1 }
| atoms_list WORD SPACE FLOAT SPACE FLOAT SPACE FLOAT EOL { output_of_string $4 $6 $8 $2 :: $1 }
| atoms_list WORD SPACE FLOAT SPACE FLOAT SPACE FLOAT SPACE EOL { output_of_string $4 $6 $8 $2 :: $1 }
| atoms_list INTEGER SPACE FLOAT SPACE FLOAT SPACE FLOAT EOL { output_of_int $4 $6 $8 $2 :: $1 }
| atoms_list INTEGER SPACE FLOAT SPACE FLOAT SPACE FLOAT SPACE EOL { output_of_int $4 $6 $8 $2 :: $1 }
| atoms_list SPACE WORD SPACE FLOAT SPACE FLOAT SPACE FLOAT EOL { output_of_string $5 $7 $9 $3 :: $1 }
| atoms_list SPACE WORD SPACE FLOAT SPACE FLOAT SPACE FLOAT SPACE EOL { output_of_string $5 $7 $9 $3 :: $1 }
| atoms_list SPACE INTEGER SPACE FLOAT SPACE FLOAT SPACE FLOAT EOL { output_of_int $5 $7 $9 $3 :: $1 }
| atoms_list SPACE INTEGER SPACE FLOAT SPACE FLOAT SPACE FLOAT SPACE EOL { output_of_int $5 $7 $9 $3 :: $1 }
;

12
Utils/Angstrom.ml Normal file
View File

@ -0,0 +1,12 @@
type t = {
x : float ;
y : float ;
z : float ;
}
let make { Point.x ; y ; z } = { x ; y ; z }
(* = %identity *)

10
Utils/Angstrom.mli Normal file
View File

@ -0,0 +1,10 @@
type t = private {
x : float ;
y : float ;
z : float ;
}
val make : Point.t -> t

12
Utils/Bohr.ml Normal file
View File

@ -0,0 +1,12 @@
type t = {
x : float ;
y : float ;
z : float ;
}
let make { Point.x ; y ; z } = { x ; y ; z }
(* = %identity *)

10
Utils/Bohr.mli Normal file
View File

@ -0,0 +1,10 @@
type t = private {
x : float ;
y : float ;
z : float ;
}
val make : Point.t -> t

1
Utils/Constants.ml
View File

@ -8,3 +8,4 @@ let pi_inv = 1. /. pi
let two_over_sq_pi = 2. /. (sqrt pi)
let a0 = 0.529_177_210_67
let a0_inv = 1. /. a0

8
Utils/Constants.mli Normal file
View File

@ -0,0 +1,8 @@
val cutoff : float
val pi : float
val sq_pi : float
val sq_pi_over_two : float
val pi_inv : float
val two_over_sq_pi : float
val a0 : float
val a0_inv : float

121
Utils/Coordinate.ml
View File

@ -1,89 +1,76 @@
open Util
type bohr = Bohr.t
type angstrom = Angstrom.t
type t =
| Bohr of (float * float * float)
| Angstrom of (float * float * float)
type t = bohr
let zero = Bohr (0., 0., 0.)
type axis = X | Y | Z
let of_float_triplet (x,y,z) = function
| `Bohr -> Bohr (x,y,z)
| `Angstrom -> Angstrom (x,y,z)
let a_to_b a = Constants.a0_inv *. a
let b_to_a b = Constants.a0 *. b
let of_3_floats x y z =
of_float_triplet (x,y,z)
let bohr_to_angstrom { Bohr.x ; y ; z } =
Angstrom.make
{
Point.
x = b_to_a x ;
y = b_to_a y ;
z = b_to_a z ;
}
let to_string t =
let result (x,y,z) =
(string_of_float x)^" "^(string_of_float y)^" "^(string_of_float z)
in
match t with
| Bohr x -> (result x) ^ " Bohr"
| Angstrom x -> (result x) ^ " Angstrom"
let angstrom_to_bohr { Angstrom.x ; y ; z } =
Bohr.make
{
Point.
x = a_to_b x ;
y = a_to_b y ;
z = a_to_b z ;
}
let extract_float_tuple = function
| Bohr a
| Angstrom a -> a
let zero =
Bohr.make { Point.x = 0. ; y = 0. ; z = 0. }
(** Linear algebra *)
let (|.) s a =
match a with
| Bohr (x,y,z) -> Bohr ( s*.x, s*.y, s*.z )
| Angstrom (x,y,z) -> Angstrom ( s*.x, s*.y, s*.z )
let ( |. ) s { Bohr.x ; y ; z } =
Bohr.make {
Point.
x = s *. x ;
y = s *. y ;
z = s *. z ;
}
let to_Angstrom = function
| Angstrom a -> Angstrom a
| Bohr a -> Angstrom (Constants.a0 |. Bohr a |> extract_float_tuple)
let ( |+ ) { Bohr.x = x1 ; y = y1 ; z = z1 } { Bohr.x = x2 ; y = y2 ; z = z2 } =
Bohr.make {
Point.
x = x1 +. x2 ;
y = y1 +. y2 ;
z = z1 +. z2 ;
}
let to_Bohr = function
| Angstrom a -> Bohr (1./.Constants.a0 |. Angstrom a |> extract_float_tuple)
| Bohr a -> Bohr a
let (|-), (|+) =
let rec op f p q =
match (p, q) with
| (Angstrom a, Angstrom b) -> Angstrom (f a b)
| (Bohr a, Bohr b) -> Bohr (f a b)
| (Angstrom a, Bohr b) -> op f (to_Bohr p) q
| (Bohr a, Angstrom b) -> op f p (to_Bohr q)
in
(op (fun (x,y,z) (x',y',z') -> ( x-.x', y-.y', z-.z' )) ,
op (fun (x,y,z) (x',y',z') -> ( x+.x', y+.y', z+.z' ))
)
let ( |- ) { Bohr.x = x1 ; y = y1 ; z = z1 } { Bohr.x = x2 ; y = y2 ; z = z2 } =
Bohr.make {
Point.
x = x1 -. x2 ;
y = y1 -. y2 ;
z = z1 -. z2 ;
}
let neg a = -1. |. a
let rec dot p q =
match (p,q) with
| Bohr (x,y,z), Bohr (x',y',z') -> x*.x' +. y*.y' +. z*.z'
| _ -> dot (to_Bohr p) (to_Bohr q)
let dot { Bohr.x = x1 ; y = y1 ; z = z1 } { Bohr.x = x2 ; y = y2 ; z = z2 } =
x1 *. x2 +. y1 *. y2 +. z1 *. z2
let norm u =
sqrt @@ dot u u
sqrt ( dot u u )
let rec to_tuple a =
to_Bohr a |> extract_float_tuple
let x a =
let (result, _, _) = extract_float_tuple @@ to_Bohr a in
result
let y a =
let (_, result, _) = extract_float_tuple @@ to_Bohr a in
result
let z a =
let (_, _, result) = extract_float_tuple @@ to_Bohr a in
result
let coord a = function
| 0 -> x a
| 1 -> y a
| 2 -> z a
| _ -> raise (Invalid_argument "Coordinate")
let get axis { Bohr.x ; y ; z } =
match axis with
| X -> x
| Y -> y
| Z -> z

36
Utils/Coordinate.mli
View File

@ -1,19 +1,19 @@
type t
val to_Angstrom : t -> t
val to_Bohr : t -> t
val zero : t
val of_float_triplet : (float * float * float) -> [< `Angstrom | `Bohr ] -> t
val of_3_floats : float -> float -> float -> [< `Angstrom | `Bohr ] -> t
val ( |. ) : float -> t -> t
val ( |- ) : t -> t -> t
val ( |+ ) : t -> t -> t
val neg : t -> t
val dot : t -> t -> float
val norm : t -> float
val to_string : t -> string
val to_tuple : t -> (float * float * float)
val x : t -> float
val y : t -> float
val z : t -> float
val coord : t -> int -> float
type bohr = Bohr.t
type angstrom = Angstrom.t
type t = bohr
type axis = X | Y | Z
val bohr_to_angstrom : bohr -> angstrom
val angstrom_to_bohr : angstrom -> bohr
val zero : bohr
val ( |. ) : float -> bohr -> bohr
val ( |+ ) : bohr -> bohr -> bohr
val ( |- ) : bohr -> bohr -> bohr
val neg : bohr -> bohr
val dot : bohr -> bohr -> float
val norm : bohr -> float
val get : axis -> bohr -> float

26
Utils/Coordinate_type.ml Normal file
View File

@ -0,0 +1,26 @@
type point = Point.t
type bohr = Bohr.t
type angstrom = Angstrom.t
let a_to_b a = Constants.a0_inv *. a
let b_to_a b = Constants.a0 *. b
let bohr_to_angstrom { Bohr.x ; y ; z } =
Angstrom.make
{
Point.
x = b_to_a x ;
y = b_to_a y ;
z = b_to_a z ;
}
let angstrom_to_bohr { Angstrom.x ; y ; z } =
Bohr.make
{
Point.
x = a_to_b x ;
y = a_to_b y ;
z = a_to_b z ;
}

8
Utils/Coordinate_type.mli Normal file
View File

@ -0,0 +1,8 @@
type bohr = Bohr.t
type angstrom = Angstrom.t
val bohr_to_angstrom : bohr -> angstrom
val angstrom_to_bohr : angstrom -> bohr

7
Utils/Point.mli Normal file
View File

@ -0,0 +1,7 @@
type t = {
x : float ;
y : float ;
z : float ;
}