LCPQ/QCaml
10
1
Fork 0
mirror of https://gitlab.com/scemama/QCaml.git synced 2024-12-22 04:13:33 +01:00
Code Releases Activity

Introduces Simulation

Browse Source
This commit is contained in:
Anthony Scemama 2018-02-09 00:37:25 +01:00
parent 659a9a3fcf
commit 730f6a9e3f
10 changed files with 236 additions and 211 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

41
Basis/Basis.ml
View File

@ -1,4 +1,11 @@
type t = Contracted_shell.t array
type t =
{
size : int;
contracted_shells : Contracted_shell.t array;
}
let size b = b.size
let contracted_shells b = b.contracted_shells
(** Returns an array of the basis set per atom *)
let of_nuclei_and_general_basis n b =
@ -22,10 +29,19 @@ let of_nuclei_and_general_basis n b =
(Contracted_shell.index result.(i-1)) +
(Array.length (Contracted_shell.powers result.(i-1))))
) result ;
result
let size =
let n = ref 0 in
for i=0 to (Array.length result) - 1 do
n := !n + (Array.length (Contracted_shell.powers (result.(i))))
done; !n
in
{ contracted_shells = result ; size }
let to_string b =
let b = b.contracted_shells in
let line ="
-----------------------------------------------------------------------
" in
@ -46,21 +62,10 @@ let to_string b =
)
^ line
let file : string option ref = ref None
let set_file f =
file := Some f
let general_basis = lazy(
match !file with
| None -> failwith "basis set file not defined"
| Some filename -> Gamess_reader.read ~filename
)
let basis = lazy (
of_nuclei_and_general_basis
(Lazy.force Nuclei.nuclei) (Lazy.force general_basis)
)
let of_nuclei_and_basis_filename ~nuclei ~filename =
let general_basis =
Gamess_reader.read ~filename
in
of_nuclei_and_general_basis nuclei general_basis

29
Basis/Basis.mli
View File

@ -1,4 +1,10 @@
type t = Contracted_shell.t array
type t
(** Number of contracted Gaussians *)
val size : t -> int
(** Array of contracted shells *)
val contracted_shells : t -> Contracted_shell.t array
(** Returns an array of the basis set per atom *)
@ -9,22 +15,7 @@ val of_nuclei_and_general_basis : Nuclei.t -> General_basis.t list -> t
val to_string : t -> string
(** Mutates the state of the file variable to Some f. Required for command-line
interface.
*)
val set_file : string -> unit
(** Create a basis using the coordinates of Nuclei and a the filename of
the general basis set *)
val of_nuclei_and_basis_filename : nuclei:Nuclei.t -> filename:string -> t
(** Basis set file read and parsed. Requires the set_file function to have
been called before.
*)
val general_basis :
(Element.t * General_basis.general_contracted_shell array) list lazy_t
(** Global variable which sets the basis set of the current run.
Lazy evaluated from the nuclear coordinates (Nuclei.nuclei) and
the basis set file (general_basis). The set_file function has
to be called before the basis is used.
*)
val basis : Contracted_shell.t array lazy_t

106
Basis/ERI.ml
View File

@ -4,6 +4,8 @@ open Util
open Constants
open Bigarray
type t = (float, float64_elt, fortran_layout) Bigarray.Genarray.t
(** (00|00)^m : Fundamental electron repulsion integral
$ \int \int \phi_p(r1) 1/r_{12} \phi_q(r2) dr_1 dr_2 $
@ -59,8 +61,7 @@ let index i j k l =
*)
(** Write all integrals to a file with the <ij|kl> convention *)
let to_file ~filename basis =
let of_basis basis =
let to_int_tuple x =
let open Zkey in
match to_int_tuple Kind_3 x with
@ -68,19 +69,16 @@ let to_file ~filename basis =
| _ -> assert false
in
let oc = open_out filename in
let n =
Basis.size basis
and shell =
Basis.contracted_shells basis
in
Printf.printf "%d shells\n" (Array.length basis);
(* Pre-compute all shell pairs *)
let shell_pairs =
ContractedShellPair.shell_pairs basis
in
let indices_of_shell_pairs =
ContractedShellPair.indices shell_pairs
in
let unique_shell_pairs =
ContractedShellPair.unique shell_pairs
ContractedShellPair.shell_pairs shell
in
(* Pre-compute diagonal integrals for Schwartz *)
@ -96,8 +94,8 @@ let to_file ~filename basis =
in
let icount = ref 0 in
for i=0 to (Array.length basis) - 1 do
print_int (Contracted_shell.index basis.(i)) ; print_newline ();
for i=0 to (Array.length shell) - 1 do
print_int (Contracted_shell.index shell.(i)) ; print_newline ();
for j=0 to i do
let schwartz_p, schwartz_p_max = schwartz.(i).(j) in
if (schwartz_p_max >= cutoff) then
@ -108,35 +106,26 @@ let to_file ~filename basis =
(* 4D data initialization *)
let eri_array =
let n = ref 0 in
for i=0 to (Array.length basis) - 1 do
n := !n + (Array.length (Contracted_shell.powers (basis.(i))))
done;
let n = !n in
Genarray.create Float64 c_layout [| n ; n ; n ; n|]
Genarray.create Float64 fortran_layout [| n ; n ; n ; n|]
in
Genarray.fill eri_array 0.;
(* Compute ERIs *)
let t0 = Unix.gettimeofday () in
let t0 = Unix.gettimeofday () in
let inn = ref 0 and out = ref 0 in
(*
for i=0 to (Array.length basis) - 1 do
print_int (Contracted_shell.index basis.(i)) ; print_newline ();
for i=0 to (Array.length shell) - 1 do
print_int (Contracted_shell.index shell.(i)) ; print_newline ();
for j=0 to i do
*)
List.iter (fun shell_p ->
let i,j =
Hashtbl.find indices_of_shell_pairs (ContractedShellPair.hash shell_p)
in
assert (ContractedShellPair.equivalent (ContractedShellPair.create basis.(i) basis.(j)) shell_p);
let schwartz_p, schwartz_p_max = schwartz.(i).(j) in
try
if (schwartz_p_max < cutoff) then raise NullIntegral;
let
shell_p = shell_pairs.(i).(j)
in
for k=0 to i do
for l=0 to k do
let schwartz_q, schwartz_q_max = schwartz.(k).(l) in
@ -165,21 +154,18 @@ let to_file ~filename basis =
contracted_class_shell_pairs_vec ~schwartz_p ~schwartz_q shell_p shell_q
in
Hashtbl.find_all indices_of_shell_pairs (ContractedShellPair.hash shell_p)
|> List.iter (fun (i,j) ->
Printf.printf "%d %d\n" i j;
(* Write the data in the output file *)
Array.iteri (fun i_c powers_i ->
let i_c = (Contracted_shell.index basis.(i)) + i_c in
let i_c = (Contracted_shell.index shell.(i)) + i_c + 1 in
let xi = to_int_tuple powers_i in
Array.iteri (fun j_c powers_j ->
let j_c = (Contracted_shell.index basis.(j)) + j_c in
let j_c = (Contracted_shell.index shell.(j)) + j_c + 1 in
let xj = to_int_tuple powers_j in
Array.iteri (fun k_c powers_k ->
let k_c = (Contracted_shell.index basis.(k)) + k_c in
let k_c = (Contracted_shell.index shell.(k)) + k_c + 1 in
let xk = to_int_tuple powers_k in
Array.iteri (fun l_c powers_l ->
let l_c = (Contracted_shell.index basis.(l)) + l_c in
let l_c = (Contracted_shell.index shell.(l)) + l_c + 1 in
let xl = to_int_tuple powers_l in
let key =
if swap then
@ -203,36 +189,36 @@ let to_file ~filename basis =
)
else
out := !out + 1;
) (Contracted_shell.powers basis.(l))
) (Contracted_shell.powers basis.(k))
) (Contracted_shell.powers basis.(j))
) (Contracted_shell.powers basis.(i));
);
) (Contracted_shell.powers shell.(l))
) (Contracted_shell.powers shell.(k))
) (Contracted_shell.powers shell.(j))
) (Contracted_shell.powers shell.(i));
with NullIntegral -> ()
done;
done;
with NullIntegral -> ()
) unique_shell_pairs;
(*
done;
done;
*)
Printf.printf "Computed %d non-zero ERIs in %f seconds\n" !inn (Unix.gettimeofday () -. t0);
(* Print ERIs *)
for i_c=1 to (Genarray.nth_dim eri_array 0) do
for j_c=1 to (Genarray.nth_dim eri_array 2) do
for k_c=1 to (Genarray.nth_dim eri_array 1) do
for l_c=1 to (Genarray.nth_dim eri_array 3) do
let value = eri_array.{(i_c-1),(k_c-1),(j_c-1),(l_c-1)} in
if (abs_float value > cutoff) then
Printf.fprintf oc " %5d %5d %5d %5d%20.15f\n" i_c k_c j_c l_c value;
done;
done;
done;
done;
Printf.printf "In: %d Out:%d\n" !inn !out ;
Printf.printf "Computed ERIs in %f seconds\n" (Unix.gettimeofday () -. t0);
eri_array
(** Write all integrals to a file with the <ij|kl> convention *)
let to_file ~filename eri_array =
let oc = open_out filename in
(* Print ERIs *)
for l_c=1 to (Genarray.nth_dim eri_array 3) do
for k_c=1 to l_c do
for j_c=1 to l_c do
for i_c=1 to k_c do
let value = eri_array.{i_c,j_c,k_c,l_c} in
if (abs_float value > cutoff) then
Printf.fprintf oc " %5d %5d %5d %5d%20.15f\n" i_c j_c k_c l_c value;
done;
done;
done;
done;
close_out oc
(*

77
Basis/KinInt.ml
View File

@ -1,5 +1,9 @@
open Util
open Constants
open Lacaml.D
type t = Mat.t
(** Computes all the kinetic integrals of the contracted shell pair *)
let contracted_class shell_a shell_b : float Zmap.t =
@ -94,48 +98,65 @@ let contracted_class shell_a shell_b : float Zmap.t =
result
(** Write all kinetic integrals to a file *)
let to_file ~filename basis =
let to_int_tuple x =
(** Create kinetic energy matrix *)
let of_basis basis =
let to_int_tuple x =
let open Zkey in
match to_int_tuple Kind_3 x with
| Three x -> x
| _ -> assert false
in
let oc = open_out filename in
for i=0 to (Array.length basis) - 1 do
print_int (Contracted_shell.index basis.(i)) ; print_newline ();
for j=0 to i do
let n =
Basis.size basis
and shell =
Basis.contracted_shells basis
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 basis.(i) basis.(j)
contracted_class shell.(i) shell.(j)
in
(* Write the data in the output file *)
Array.iteri (fun i_c powers_i ->
let i_c = Contracted_shell.index basis.(i) + i_c + 1 in
let xi = to_int_tuple powers_i in
Array.iteri (fun j_c powers_j ->
let j_c = Contracted_shell.index basis.(j) + j_c + 1 in
let xj = to_int_tuple powers_j in
let key =
Array.iteri (fun j_c powers_j ->
let j_c = Contracted_shell.index shell.(j) + j_c + 1 in
let xj = to_int_tuple powers_j in
Array.iteri (fun i_c powers_i ->
let i_c = Contracted_shell.index shell.(i) + i_c + 1 in
let xi = to_int_tuple powers_i in
let key =
Zkey.of_int_tuple (Zkey.Six (xi,xj))
in
let value =
try
Zmap.find cls key
with Not_found -> failwith "Bug in Kinetic integrals"
let value =
try Zmap.find cls key
with Not_found -> failwith "Bug in kinetic integrals"
in
if (abs_float value > cutoff) then
Printf.fprintf oc "%4d %4d %20.12e\n"
i_c j_c value
) (Contracted_shell.powers basis.(j))
) (Contracted_shell.powers basis.(i));
result.{i_c,j_c} <- value
) (Contracted_shell.powers shell.(i));
) (Contracted_shell.powers shell.(j))
done;
done;
Mat.detri result;
result
(** Write all kinetic integrals to a file *)
let to_file ~filename kinetic =
let oc = open_out filename in
let n =
Mat.dim1 kinetic
in
for j=1 to n do
for i=1 to j do
if (abs_float kinetic.{i,j} > cutoff) then
Printf.fprintf oc "%4d %4d %20.12e\n" i j kinetic.{i,j}
done;
done;
close_out oc

69
Basis/NucInt.ml
View File

@ -2,7 +2,9 @@
open Util
open Constants
open Bigarray
open Lacaml.D
type t = Mat.t
(** (0|0)^m : Fundamental electron-nucleus repulsion integral
$ \int \phi_p(r1) 1/r_{C} dr_1 $
@ -32,7 +34,8 @@ let contracted_class_shell_pair shell_p geometry: float Zmap.t =
let cutoff2 = cutoff *. cutoff
exception NullIntegral
let to_file ~filename basis geometry =
let of_basis_nuclei basis nuclei =
let to_int_tuple x =
let open Zkey in
match to_int_tuple Kind_3 x with
@ -40,32 +43,22 @@ let to_file ~filename basis geometry =
| _ -> assert false
in
let oc = open_out filename in
let n =
Basis.size basis
and shell =
Basis.contracted_shells basis
in
let eni_array = Mat.create n n in
(* Pre-compute all shell pairs *)
let shell_pairs =
Array.mapi (fun i shell_a -> Array.map (fun shell_b ->
ContractedShellPair.create shell_a shell_b) (Array.sub basis 0 (i+1)) ) basis
ContractedShellPair.create shell_a shell_b) (Array.sub shell 0 (i+1)) ) shell
in
Printf.printf "%d shells\n" (Array.length basis);
let eni_array =
let n = ref 0 in
for i=0 to (Array.length basis) - 1 do
n := !n + (Array.length (Contracted_shell.powers basis.(i)))
done;
let n = !n in
Array2.create Float64 c_layout n n
in
Array2.fill eni_array 0.;
(* Compute Integrals *)
let t0 = Unix.gettimeofday () in
let inn = ref 0 and out = ref 0 in
for i=0 to (Array.length basis) - 1 do
print_int (Contracted_shell.index basis.(i)) ; print_newline ();
for i=0 to (Array.length shell) - 1 do
for j=0 to i do
let
shell_p = shell_pairs.(i).(j)
@ -73,15 +66,15 @@ let to_file ~filename basis geometry =
(* Compute all the integrals of the class *)
let cls =
contracted_class_shell_pair shell_p geometry
contracted_class_shell_pair shell_p nuclei
in
(* Write the data in the output file *)
Array.iteri (fun i_c powers_i ->
let i_c = (Contracted_shell.index basis.(i)) + i_c + 1 in
let i_c = (Contracted_shell.index shell.(i)) + i_c + 1 in
let xi = to_int_tuple powers_i in
Array.iteri (fun j_c powers_j ->
let j_c = (Contracted_shell.index basis.(j)) + j_c + 1 in
let j_c = (Contracted_shell.index shell.(j)) + j_c + 1 in
let xj = to_int_tuple powers_j in
let key =
Zkey.of_int_tuple (Zkey.Six (xi,xj))
@ -89,27 +82,25 @@ let to_file ~filename basis geometry =
let value =
Zmap.find cls key
in
if (abs_float value > cutoff) then
(inn := !inn + 1;
eni_array.{(i_c-1),(j_c-1)} <- value;
)
else
out := !out + 1;
) (Contracted_shell.powers basis.(j))
) (Contracted_shell.powers basis.(i));
eni_array.{j_c,i_c} <- value;
) (Contracted_shell.powers shell.(j))
) (Contracted_shell.powers shell.(i));
done;
done;
Mat.detri eni_array;
eni_array
Printf.printf "Computed %d non-zero ENIs in %f seconds\n" !inn (Unix.gettimeofday () -. t0);
(* Print ENIs *)
for i_c=1 to (Array2.dim1 eni_array) do
for j_c=1 to (Array2.dim2 eni_array) do
let value = eni_array.{(i_c-1),(j_c-1)} in
let to_file ~filename eni_array =
let n = Mat.dim1 eni_array in
let oc = open_out filename in
for j=1 to n do
for i=1 to j do
let value = eni_array.{i,j} in
if (value <> 0.) then
Printf.fprintf oc " %5d %5d %20.15f\n" i_c j_c value;
Printf.fprintf oc " %5d %5d %20.15f\n" i j value;
done;
done;
Printf.printf "In: %d Out:%d\n" !inn !out ;
close_out oc

70
Basis/Overlap.ml
View File

@ -1,5 +1,9 @@
open Util
open Constants
open Lacaml.D
type t = Mat.t
(** Computes all the overlap integrals of the contracted shell pair *)
let contracted_class shell_a shell_b : float Zmap.t =
@ -66,11 +70,8 @@ let contracted_class shell_a shell_b : float Zmap.t =
result
(** Write all overlap integrals to a file *)
let to_file ~filename basis =
(** Create overlap matrix *)
let of_basis basis =
let to_int_tuple x =
let open Zkey in
match to_int_tuple Kind_3 x with
@ -78,35 +79,56 @@ let to_file ~filename basis =
| _ -> assert false
in
let oc = open_out filename in
for i=0 to (Array.length basis) - 1 do
print_int (Contracted_shell.index basis.(i)) ; print_newline ();
for j=0 to i do
let n =
Basis.size basis
and shell =
Basis.contracted_shells basis
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 basis.(i) basis.(j)
contracted_class shell.(i) shell.(j)
in
(* Write the data in the output file *)
Array.iteri (fun i_c powers_i ->
let i_c = Contracted_shell.index basis.(i) + i_c + 1 in
let xi = to_int_tuple powers_i in
Array.iteri (fun j_c powers_j ->
let j_c = Contracted_shell.index basis.(j) + j_c + 1 in
let xj = to_int_tuple powers_j in
Array.iteri (fun j_c powers_j ->
let j_c = Contracted_shell.index shell.(j) + j_c + 1 in
let xj = to_int_tuple powers_j in
Array.iteri (fun i_c powers_i ->
let i_c = Contracted_shell.index shell.(i) + i_c + 1 in
let xi = to_int_tuple powers_i in
let key =
Zkey.of_int_tuple (Zkey.Six (xi,xj))
in
let value =
try
Zmap.find cls key
try Zmap.find cls key
with Not_found -> failwith "Bug in overlap integrals"
in
if (abs_float value > cutoff) then
Printf.fprintf oc "%4d %4d %20.12e\n"
i_c j_c value
) (Contracted_shell.powers basis.(j))
) (Contracted_shell.powers basis.(i));
result.{i_c,j_c} <- value
) (Contracted_shell.powers shell.(i));
) (Contracted_shell.powers shell.(j))
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

14
Nuclei/Nuclei.ml
View File

@ -23,6 +23,7 @@ let of_zmt_file ~filename =
|> Zmatrix.to_xyz
|> Array.map (fun (e,x,y,z) -> (e, Coordinate.of_3_floats x y z ))
let to_string atoms =
"
Nuclear Coordinates (Angstrom)
@ -49,14 +50,7 @@ let to_string atoms =
"
let file : string option ref = ref None
let set_file f =
file := Some f
let nuclei = lazy(
match !file with
| None -> failwith "coordinate file not defined"
| Some filename -> of_xyz_file ~filename
)
let of_filename ~filename =
of_xyz_file filename

2
README.rst
View File

@ -2,8 +2,8 @@ Requirements
------------
* gmp : GNU Multiple Precision arithmetic library
* zarith : Arbitrary-precision integers
* BLAS/LAPACK : Linear algebra
* Zarith : Arbitrary-precision integers
* LaCaml : LAPACK OCaml interface
* SklMl : Parallel skeletons for OCaml

2
_tags
View File

@ -1,3 +1,3 @@
true: package(str,unix,bigarray,zarith)
true: package(str,unix,bigarray,zarith,lacaml)
<*.byte> : linkdep(Utils/math_functions.o), custom
<*.native>: linkdep(Utils/math_functions.o)

37
run_integrals.ml
View File

@ -1,10 +1,12 @@
let out_file : string option ref = ref None
let out_file : string option ref = ref None
let basis_file : string option ref = ref None
let nuclei_file : string option ref = ref None
let speclist = [
( "-b" , Arg.String Basis.set_file ,
( "-b" , Arg.String (fun x -> basis_file := Some x),
"File containing the atomic basis set") ;
( "-c" , Arg.String Nuclei.set_file ,
( "-c" , Arg.String (fun x -> nuclei_file := Some x),
"File containing the nuclear coordinates") ;
( "-o" , Arg.String (fun x -> out_file := Some x) ,
"Output file") ;
@ -19,18 +21,31 @@ let run ~out =
match out with
| None -> raise (Invalid_argument "Output file should be specified with -o")
| Some x -> x
and basis_file =
match !basis_file with
| None -> raise (Invalid_argument "Basis set file should be specified with -b")
| Some x -> x
and nuclei_file =
match !nuclei_file with
| None -> raise (Invalid_argument "Basis set file should be specified with -c")
| Some x -> x
in
let nuclei = Lazy.force Nuclei.nuclei in
print_endline @@ Nuclei.to_string nuclei;
let s =
Simulation.of_filenames ~nuclei:nuclei_file ~basis:basis_file
in
let basis = Lazy.force Basis.basis in
print_endline @@ Basis.to_string basis;
print_endline @@ Nuclei.to_string s.Simulation.nuclei;
print_endline @@ Basis.to_string s.Simulation.basis;
Overlap.to_file ~filename:(out_file^".overlap") basis;
NucInt.to_file ~filename:(out_file^".nuc") basis nuclei;
KinInt.to_file ~filename:(out_file^".kin") basis;
ERI.to_file ~filename:(out_file^".eri") basis
let overlap = Lazy.force s.Simulation.overlap in
let eN_ints = Lazy.force s.Simulation.eN_ints in
let kin_ints = Lazy.force s.Simulation.kin_ints in
let ee_ints = Lazy.force s.Simulation.ee_ints in
Overlap.to_file ~filename:(out_file^".overlap") overlap;
NucInt.to_file ~filename:(out_file^".nuc") eN_ints;
KinInt.to_file ~filename:(out_file^".kin") kin_ints;
ERI.to_file ~filename:(out_file^".eri") ee_ints
let () =