LCPQ/QCaml
10
1
Fork 0
mirror of https://gitlab.com/scemama/QCaml.git synced 2024-11-06 22:23:42 +01:00
Code Releases Activity

Introduced multipole integrals

Browse Source
This commit is contained in:
Anthony Scemama 2020-04-16 19:49:23 +02:00
parent a1c29fd3f1
commit 11151c997e
11 changed files with 383 additions and 85 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
Basis/AOBasis.ml
View File

@ -5,6 +5,7 @@ type t =
{
basis : Basis.t ;
overlap : Overlap.t lazy_t;
multipole : Multipole.t lazy_t;
ortho : Orthonormalization.t lazy_t;
eN_ints : NucInt.t lazy_t;
kin_ints : KinInt.t lazy_t;
@ -17,6 +18,7 @@ type t =
let basis t = t.basis
let overlap t = Lazy.force t.overlap
let multipole t = Lazy.force t.multipole
let ortho t = Lazy.force t.ortho
let eN_ints t = Lazy.force t.eN_ints
let kin_ints t = Lazy.force t.kin_ints
@ -62,7 +64,11 @@ let make ~cartesian ~basis ?f12 nuclei =
ScreenedERI.of_basis basis
) in
{ basis ; overlap ; ortho ; eN_ints ; kin_ints ; ee_ints ;
let multipole = lazy (
Multipole.of_basis basis
) in
{ basis ; overlap ; multipole ; ortho ; eN_ints ; kin_ints ; ee_ints ;
ee_lr_ints ; f12_ints ; f12_over_r12_ints ; cartesian ;
}

3
Basis/AOBasis.mli
View File

@ -10,6 +10,9 @@ val basis : t -> Basis.t
val overlap : t -> Overlap.t
(** Overlap matrix *)
val multipole : t -> Multipole.t
(** Multipole matrices *)
val ortho : t -> Orthonormalization.t
(** Orthonormalization matrix of the overlap *)

6
Basis/Basis.ml
View File

@ -93,6 +93,12 @@ let of_nuclei_and_basis_filename ~nuclei filename =
in
of_nuclei_and_general_basis nuclei general_basis
let of_nuclei_and_basis_string ~nuclei str =
let general_basis =
GeneralBasis.of_string str
in
of_nuclei_and_general_basis nuclei general_basis
let of_nuclei_and_basis_filenames ~nuclei filenames =
let general_basis =

5
Basis/Basis.mli
View File

@ -18,6 +18,11 @@ val of_nuclei_and_basis_filename : nuclei:Nuclei.t -> string -> t
from a file.
*)
val of_nuclei_and_basis_string : nuclei:Nuclei.t -> string -> t
(** Same as {!of_nuclei_and_general_basis}, but taking the {!GeneralBasis.t}
from a string.
*)
val of_nuclei_and_basis_filenames : nuclei:Nuclei.t -> string list -> t
(** Same as {!of_nuclei_and_general_basis}, but taking the {!GeneralBasis.t}
from multiple files.

40
Basis/GeneralBasis.ml
View File

@ -17,23 +17,24 @@ exception No_shell
exception Malformed_shell of string
let read_shell ic =
let read_shell line_stream =
try
let shell, n =
try
let line = input_line ic in
let line = Stream.next line_stream in
if String.trim line = "$END" then raise End_of_file;
Scanf.sscanf line " %c %d " (fun shell n -> shell, n)
with
| End_of_file -> raise No_shell
| Stream.Failure -> raise No_shell
| Scanf.Scan_failure m -> raise (Malformed_shell m)
in
let rec loop = function
| 0 -> []
| i -> let contraction =
let line = (input_line ic) in
let line = Stream.next line_stream in
try Scanf.sscanf line " %_d %f %f "
(fun exponent coefficient -> { exponent ; coefficient })
with _ -> raise (Malformed_shell (Printf.sprintf
@ -50,36 +51,46 @@ let read_shell ic =
let read_element ic =
let read_element line_stream =
try
let line = Stream.next line_stream in
let element =
Scanf.sscanf (input_line ic) " %s " Element.of_string
Scanf.sscanf line " %s " Element.of_string
in
let rec loop () =
match read_shell ic with
match read_shell line_stream with
| Some shell -> shell :: loop ()
| None -> []
in
Some (element, Array.of_list (loop ()) )
with
| End_of_file -> None
| Stream.Failure -> None
let read filename =
let ic = open_in filename in
let read_stream line_stream =
let rec loop accu =
try
match read_element ic with
match read_element line_stream with
| Some e -> loop (e :: accu)
| None -> accu
with
Element.ElementError _ -> loop accu
in
loop []
loop []
let read filename =
let ic = open_in filename in
let line_stream =
Stream.from (fun _ ->
try Some (input_line ic)
with End_of_file -> None )
in
let result = read_stream line_stream in
close_in ic;
result
let combine basis_list =
@ -124,3 +135,8 @@ let to_string (name, contracted_shell_array) =
Printf.sprintf "%s\n%s" name (string_of_contracted_shell_array contracted_shell_array)
let of_string input_string =
String.split_on_char '\n' input_string
|> Stream.of_list
|> read_stream

16
Basis/GeneralBasis.mli
View File

@ -46,7 +46,7 @@ exception Malformed_shell of string
val read : string -> t
(** Reads a basis set file and return an association list where
(** Reads a basis-set file and return an association list where
the key is an {!Element.t} and the value is the parsed basis set.
*)
@ -65,8 +65,10 @@ val read_many : string list -> t
val read_element : in_channel -> element_basis option
(** Reads an element from the input channel [ic]. For example,
val read_element : string Stream.t -> element_basis option
(** Reads an element from the input [string Stream]. The [string Stream] is a
stream of lines, like a text file split on the end-of-line character.
For example,
{[
HYDROGEN
S 3
@ -95,8 +97,10 @@ Some
*)
val read_shell : in_channel -> general_contracted_shell option
(** Reads a shell from the input channel [ic]. For example,
val read_shell : string Stream.t -> general_contracted_shell option
(** Reads a shell from the input [string Stream]. The [string Stream] is a
stream of lines, like a text file split on the end-of-line character.
For example,
{[
S 3
1 13.0100000 0.0196850
@ -119,4 +123,6 @@ Some
val to_string : string * (general_contracted_shell array) -> string
(** Pretty-prints the basis set of an {Element.t}. *)
val of_string : string -> t
(** Reads a GAMESS-formatted string. *)

199
Basis/Multipole.ml Normal file
View File

@ -0,0 +1,199 @@
open Util
open Constants
open Lacaml.D
type t = Mat.t array
(*
[| "x"; "y"; "z"; "x2"; "y2"; "z2" |]
*)
module Am = AngularMomentum
module Bs = Basis
module Co = Coordinate
module Cs = ContractedShell
module Csp = ContractedShellPair
module Po = Powers
module Psp = PrimitiveShellPair
let multiply a b =
let n = Mat.dim1 a in
let c = Mat.create n n in
Mat.cpab c a b;
c
let x0 t = t.(0)
let y0 t = t.(1)
let z0 t = t.(2)
let x1 t = t.(3)
let y1 t = t.(4)
let z1 t = t.(5)
let x2 t = t.(6)
let y2 t = t.(7)
let z2 t = t.(8)
let matrix_x t = multiply (x1 t) @@ multiply (y0 t) (z0 t)
let matrix_y t = multiply (x0 t) @@ multiply (y1 t) (z0 t)
let matrix_z t = multiply (x0 t) @@ multiply (y0 t) (z1 t)
let matrix_x2 t = multiply (x2 t) @@ multiply (y0 t) (z0 t)
let matrix_y2 t = multiply (x0 t) @@ multiply (y2 t) (z0 t)
let matrix_z2 t = multiply (x0 t) @@ multiply (y0 t) (z2 t)
let matrix_xy t = multiply (x1 t) @@ multiply (y1 t) (z0 t)
let matrix_yz t = multiply (x0 t) @@ multiply (y1 t) (z1 t)
let matrix_zx t = multiply (x1 t) @@ multiply (y0 t) (z1 t)
let cutoff = integrals_cutoff
let to_powers x =
let open Zkey in
match to_powers x with
| Six x -> x
| _ -> assert false
(** Computes all the integrals of the contracted shell pair *)
let contracted_class shell_a shell_b : float Zmap.t array =
match Csp.make shell_a shell_b with
| None -> Array.init 9 (fun _ -> Zmap.create 0)
| Some shell_p ->
begin
(* Pre-computation of integral class indices *)
let class_indices = Csp.zkey_array shell_p in
let contracted_class =
Array.init 9 (fun _ -> Array.make (Array.length class_indices) 0.)
in
let a_minus_b =
Csp.a_minus_b shell_p
in
let norm_coef_scales =
Csp.norm_scales shell_p
in
(* Compute all integrals in the shell for each pair of significant shell pairs *)
let xyz_of_int k =
match k with
| 0 -> Co.X
| 1 -> Co.Y
| _ -> Co.Z
in
List.iter (fun (coef_prod, psp) ->
(** Screening on the product of coefficients *)
if (abs_float coef_prod) > 1.e-6*.cutoff then
begin
let expo_inv = Psp.exponent_inv psp
and center_pa = Psp.center_minus_a psp
and xa = Co.get Co.X @@ Cs.center shell_a
and ya = Co.get Co.Y @@ Cs.center shell_a
and za = Co.get Co.Z @@ Cs.center shell_a
in
Array.iteri (fun i key ->
let (angMomA, angMomB) = to_powers key in
(* 1D Overlap <i|j> *)
let f k =
let xyz = xyz_of_int k in
Overlap_primitives.hvrr (Po.get xyz angMomA, Po.get xyz angMomB)
expo_inv
(Co.get xyz a_minus_b, Co.get xyz center_pa)
in
(* 1D <i|x-Xa|j> *)
let g k =
let xyz = xyz_of_int k in
Overlap_primitives.hvrr (Po.get xyz angMomA + 1, Po.get xyz angMomB)
expo_inv
(Co.get xyz a_minus_b, Co.get xyz center_pa)
in
(* 1D <i|(x-Xa)^2|j> *)
let h k =
let xyz = xyz_of_int k in
Overlap_primitives.hvrr (Po.get xyz angMomA + 2, Po.get xyz angMomB)
expo_inv
(Co.get xyz a_minus_b, Co.get xyz center_pa)
in
let norm = norm_coef_scales.(i) in
let f0, f1, f2, g0, g1, g2, h0, h1, h2 =
f 0, f 1, f 2, g 0, g 1, g 2, h 0, h 1, h 2
in
let x = g0 +. f0 *. xa in
let y = g1 +. f1 *. ya in
let z = g2 +. f2 *. za in
let x2 = h0 +. xa *. (2. *. x -. xa *. f0) in
let y2 = h1 +. ya *. (2. *. y -. ya *. f1) in
let z2 = h2 +. za *. (2. *. z -. za *. f2) in
let c = contracted_class in
let d = coef_prod *. norm in
c.(0).(i) <- c.(0).(i) +. d *. f0 ;
c.(1).(i) <- c.(1).(i) +. d *. f1 ;
c.(2).(i) <- c.(2).(i) +. d *. f2 ;
c.(3).(i) <- c.(3).(i) +. d *. x ;
c.(4).(i) <- c.(4).(i) +. d *. y ;
c.(5).(i) <- c.(5).(i) +. d *. z ;
c.(6).(i) <- c.(6).(i) +. d *. x2 ;
c.(7).(i) <- c.(7).(i) +. d *. y2 ;
c.(8).(i) <- c.(8).(i) +. d *. z2 ;
) class_indices
end
) (Csp.coefs_and_shell_pairs shell_p);
let result =
Array.map (fun c -> Zmap.create (Array.length c) ) contracted_class
in
for j=0 to Array.length result do
let rj = result.(j) in
let cj = contracted_class.(j) in
Array.iteri (fun i key -> Zmap.add rj key cj.(i)) class_indices
done;
result
end
(** Create multipole matrices *)
let of_basis basis =
let to_powers x =
let open Zkey in
match to_powers x with
| Three x -> x
| _ -> assert false
in
let n = Bs.size basis
and shell = Bs.contracted_shells basis
in
let result = Array.init 9 (fun _ -> 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
for k=0 to 8 do
Array.iteri (fun j_c powers_j ->
let j_c = Cs.index shell.(j) + j_c + 1 in
let xj = to_powers powers_j in
Array.iteri (fun i_c powers_i ->
let i_c = Cs.index shell.(i) + i_c + 1 in
let xi = to_powers powers_i in
let key =
Zkey.of_powers_six xi xj
in
let value =
try Zmap.find cls.(k) key
with Not_found -> 0.
in
result.(k).{i_c,j_c} <- value;
result.(k).{j_c,i_c} <- value;
) (Am.zkey_array (Singlet (Cs.ang_mom shell.(i))))
) (Am.zkey_array (Singlet (Cs.ang_mom shell.(j))))
done;
done;
done;
for k=0 to 8 do
Mat.detri result.(k);
done;
result

34
Basis/Multipole.mli Normal file
View File

@ -0,0 +1,34 @@
(** Multipole atomic integrals:
{% $$ \langle \chi_i | x | \chi_j \rangle $$ %}
{% $$ \langle \chi_i | y | \chi_j \rangle $$ %}
{% $$ \langle \chi_i | z | \chi_j \rangle $$ %}
{% $$ \langle \chi_i | x^2 | \chi_j \rangle $$ %}
{% $$ \langle \chi_i | y^2 | \chi_j \rangle $$ %}
{% $$ \langle \chi_i | z^2 | \chi_j \rangle $$ %}
*)
open Lacaml.D
type t
val matrix_x : t -> Mat.t
(** {% $$ \langle \chi_i | x | \chi_j \rangle $$ %} *)
val matrix_y : t -> Mat.t
(** {% $$ \langle \chi_i | y | \chi_j \rangle $$ %} *)
val matrix_z : t -> Mat.t
(** {% $$ \langle \chi_i | z | \chi_j \rangle $$ %} *)
val matrix_x2 : t -> Mat.t
(** {% $$ \langle \chi_i | x^2 | \chi_j \rangle $$ %} *)
val matrix_y2 : t -> Mat.t
(** {% $$ \langle \chi_i | y^2 | \chi_j \rangle $$ %} *)
val matrix_z2 : t -> Mat.t
(** {% $$ \langle \chi_i | z^2 | \chi_j \rangle $$ %} *)
val of_basis : Basis.t -> t

124
Basis/Overlap.ml
View File

@ -18,77 +18,77 @@ module Psp = PrimitiveShellPair
let cutoff = integrals_cutoff
let to_powers x =
let to_powers x =
let open Zkey in
match to_powers x with
| Six x -> x
| _ -> assert false
(** Computes all the overlap integrals of the contracted shell pair *)
let contracted_class shell_a shell_b : float Zmap.t =
let contracted_class shell_a shell_b : float Zmap.t =
match Csp.make shell_a shell_b with
| Some shell_p ->
| None -> Zmap.create 0
| Some shell_p ->
begin
(* Pre-computation of integral class indices *)
let class_indices = Csp.zkey_array shell_p in
(* Pre-computation of integral class indices *)
let class_indices = Csp.zkey_array shell_p in
let contracted_class =
Array.make (Array.length class_indices) 0.
in
let contracted_class =
Array.make (Array.length class_indices) 0.
in
let a_minus_b =
Csp.a_minus_b shell_p
in
let norm_coef_scales =
Csp.norm_scales shell_p
in
let a_minus_b =
Csp.a_minus_b shell_p
in
let norm_coef_scales =
Csp.norm_scales shell_p
in
(* Compute all integrals in the shell for each pair of significant shell pairs *)
(* Compute all integrals in the shell for each pair of significant shell pairs *)
let xyz_of_int k =
match k with
| 0 -> Co.X
| 1 -> Co.Y
| _ -> Co.Z
in
let xyz_of_int k =
match k with
| 0 -> Co.X
| 1 -> Co.Y
| _ -> Co.Z
in
List.iter (fun (coef_prod, psp) ->
(** Screening on thr product of coefficients *)
if (abs_float coef_prod) > 1.e-3*.cutoff then
begin
let expo_inv = Psp.exponent_inv psp
and center_pa = Psp.center_minus_a psp
in
Array.iteri (fun i key ->
let (angMomA,angMomB) = to_powers key in
let f k =
let xyz = xyz_of_int k in
Overlap_primitives.hvrr (Po.get xyz angMomA, Po.get xyz angMomB)
expo_inv
(Co.get xyz a_minus_b,
Co.get xyz center_pa)
List.iter (fun (coef_prod, psp) ->
(** Screening on the product of coefficients *)
if (abs_float coef_prod) > 1.e-6*.cutoff then
begin
let expo_inv = Psp.exponent_inv psp
and center_pa = Psp.center_minus_a psp
in
let norm = norm_coef_scales.(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
) (Csp.coefs_and_shell_pairs shell_p);
let result =
Zmap.create (Array.length contracted_class)
in
Array.iteri (fun i key -> Zmap.add result key contracted_class.(i)) class_indices;
result
end
| None -> Zmap.create 0
Array.iteri (fun i key ->
let (angMomA,angMomB) = to_powers key in
let f k =
let xyz = xyz_of_int k in
Overlap_primitives.hvrr (Po.get xyz angMomA, Po.get xyz angMomB)
expo_inv
(Co.get xyz a_minus_b,
Co.get xyz center_pa)
in
let norm = norm_coef_scales.(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
) (Csp.coefs_and_shell_pairs shell_p);
let result =
Zmap.create (Array.length contracted_class)
in
Array.iteri (fun i key -> Zmap.add result key contracted_class.(i)) class_indices;
result
end
(** Create overlap matrix *)
let of_basis basis =
let to_powers x =
let to_powers x =
let open Zkey in
match to_powers x with
| Three x -> x
@ -104,7 +104,7 @@ let of_basis basis =
for i=0 to j do
(* Compute all the integrals of the class *)
let cls =
contracted_class shell.(i) shell.(j)
contracted_class shell.(i) shell.(j)
in
Array.iteri (fun j_c powers_j ->
@ -113,11 +113,11 @@ let of_basis basis =
Array.iteri (fun i_c powers_i ->
let i_c = Cs.index shell.(i) + i_c + 1 in
let xi = to_powers powers_i in
let key =
let key =
Zkey.of_powers_six xi xj
in
let value =
try Zmap.find cls key
let value =
try Zmap.find cls key
with Not_found -> 0.
in
result.{i_c,j_c} <- value;
@ -132,7 +132,7 @@ let of_basis basis =
(** Create mixed overlap matrix *)
let of_basis_pair first_basis second_basis =
let to_powers x =
let to_powers x =
let open Zkey in
match to_powers x with
| Three x -> x
@ -140,7 +140,7 @@ let of_basis_pair first_basis second_basis =
in
let n = Bs.size first_basis
and m = Bs.size second_basis
and m = Bs.size second_basis
and first = Bs.contracted_shells first_basis
and second = Bs.contracted_shells second_basis
in
@ -150,7 +150,7 @@ let of_basis_pair first_basis second_basis =
for i=0 to (Array.length first) - 1 do
(* Compute all the integrals of the class *)
let cls =
contracted_class first.(i) second.(j)
contracted_class first.(i) second.(j)
in
Array.iteri (fun j_c powers_j ->
@ -159,11 +159,11 @@ let of_basis_pair first_basis second_basis =
Array.iteri (fun i_c powers_i ->
let i_c = Cs.index first.(i) + i_c + 1 in
let xi = to_powers powers_i in
let key =
let key =
Zkey.of_powers_six xi xj
in
let value =
try Zmap.find cls key
let value =
try Zmap.find cls key
with Not_found -> 0.
in
result.{i_c,j_c} <- value;
@ -180,7 +180,7 @@ let to_file ~filename overlap =
let oc = open_out filename in
let n =
Mat.dim1 overlap
Mat.dim1 overlap
in
for j=1 to n do

30
Nuclei/Nuclei.ml
View File

@ -3,11 +3,7 @@ open Xyz_ast
type t = (Element.t * Coordinate.t) array
let of_xyz_file filename =
let lexbuf =
let ic = open_in filename in
Lexing.from_channel ic
in
let of_xyz_lexbuf lexbuf =
let data =
Xyz_parser.input Nuclei_lexer.read_all lexbuf
in
@ -24,6 +20,30 @@ let of_xyz_file filename =
|> Array.of_list
let of_xyz_string buffer =
Zmatrix.of_string buffer
|> Zmatrix.to_xyz
|> Array.map (fun (e,x,y,z) ->
(e, Coordinate.(angstrom_to_bohr @@ make_angstrom { x ; y ; z} ))
)
let of_xyz_string input_string =
Lexing.from_string input_string
|> of_xyz_lexbuf
let of_xyz_file filename =
let ic = open_in filename in
let lexbuf =
Lexing.from_channel ic
in
let result =
of_xyz_lexbuf lexbuf
in
close_in ic;
result
let of_zmt_string buffer =
Zmatrix.of_string buffer

3
Nuclei/Nuclei.mli
View File

@ -5,6 +5,9 @@ of tuples ({!Element.t}, {!Coordinate.t}).
type t = (Element.t * Coordinate.t) array
val of_xyz_string : string -> t
(** Create from a string, in [xyz] format. *)
val of_xyz_file : string -> t
(** Create from a file, in [xyz] format. *)