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

Provided example

This commit is contained in:
Anthony Scemama 2020-10-19 18:33:02 +02:00
parent a907572217
commit 718d39924f
18 changed files with 318 additions and 139 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
.gitignore vendored
View File

@ -3,3 +3,4 @@ _build
.merlin
*.install
qcaml.opam
#*#

39
ao/lib/basis.ml
View File

@ -1,11 +1,8 @@
open Linear_algebra
type basis =
| Unknown
| Gaussian of Basis_gaussian.t
open Common
type t =
{ ao_basis : basis ;
{ ao_basis : Basis_poly.t ;
cartesian : bool
}
@ -19,25 +16,24 @@ let of_nuclei_and_basis_filename ?(kind=`Gaussian) ?operators ?(cartesian=false)
Gaussian.Basis.of_nuclei_and_basis_filename ~nuclei filename
in
let ao_basis =
Gaussian (Basis_gaussian.make ~basis ?operators ~cartesian nuclei )
Basis_poly.Gaussian (Basis_gaussian.make ~basis ?operators ~cartesian nuclei )
in
{ ao_basis ; cartesian }
| _ -> failwith "of_nuclei_and_basis_filename needs to be called with `Gaussian"
let not_implemented () =
failwith "Not implemented"
let not_implemented = Util.not_implemented
let ao_basis t = t.ao_basis
let size t =
match t.ao_basis with
| Gaussian b -> Basis_gaussian.size b
| Basis_poly.Gaussian b -> Basis_gaussian.size b
| _ -> not_implemented ()
let overlap t =
begin
match t.ao_basis with
| Gaussian b -> Basis_gaussian.overlap b
| Basis_poly.Gaussian b -> Basis_gaussian.overlap b
| _ -> not_implemented ()
end
|> Matrix.relabel
@ -45,15 +41,18 @@ let overlap t =
let multipole t =
begin
match t.ao_basis with
| Gaussian b -> Basis_gaussian.multipole b
| Basis_poly.Gaussian b ->
let m = Basis_gaussian.multipole b in
fun s ->
Gaussian_integrals.Multipole.matrix m s
|> Matrix.relabel
| _ -> not_implemented ()
end
|> Array.map Matrix.relabel
let ortho t =
begin
match t.ao_basis with
| Gaussian b -> Basis_gaussian.ortho b
| Basis_poly.Gaussian b -> Basis_gaussian.ortho b
| _ -> not_implemented ()
end
|> Matrix.relabel
@ -61,7 +60,7 @@ let ortho t =
let eN_ints t =
begin
match t.ao_basis with
| Gaussian b -> Basis_gaussian.eN_ints b
| Basis_poly.Gaussian b -> Basis_gaussian.eN_ints b
| _ -> not_implemented ()
end
|> Matrix.relabel
@ -69,7 +68,7 @@ let eN_ints t =
let kin_ints t =
begin
match t.ao_basis with
| Gaussian b -> Basis_gaussian.kin_ints b
| Basis_poly.Gaussian b -> Basis_gaussian.kin_ints b
| _ -> not_implemented ()
end
|> Matrix.relabel
@ -77,7 +76,7 @@ let kin_ints t =
let ee_ints t =
begin
match t.ao_basis with
| Gaussian b -> Basis_gaussian.ee_ints b
| Basis_poly.Gaussian b -> Basis_gaussian.ee_ints b
| _ -> not_implemented ()
end
|> Four_idx_storage.relabel
@ -85,7 +84,7 @@ let ee_ints t =
let ee_lr_ints t =
begin
match t.ao_basis with
| Gaussian b -> Basis_gaussian.ee_lr_ints b
| Basis_poly.Gaussian b -> Basis_gaussian.ee_lr_ints b
| _ -> not_implemented ()
end
|> Four_idx_storage.relabel
@ -93,7 +92,7 @@ let ee_lr_ints t =
let f12_ints t =
begin
match t.ao_basis with
| Gaussian b -> Basis_gaussian.f12_ints b
| Basis_poly.Gaussian b -> Basis_gaussian.f12_ints b
| _ -> not_implemented ()
end
|> Four_idx_storage.relabel
@ -101,7 +100,7 @@ let f12_ints t =
let f12_over_r12_ints t =
begin
match t.ao_basis with
| Gaussian b -> Basis_gaussian.f12_over_r12_ints b
| Basis_poly.Gaussian b -> Basis_gaussian.f12_over_r12_ints b
| _ -> not_implemented ()
end
|> Four_idx_storage.relabel
@ -112,7 +111,7 @@ let cartesian t = t.cartesian
let values t point =
begin
match t.ao_basis with
| Gaussian b -> Basis_gaussian.values b point
| Basis_poly.Gaussian b -> Basis_gaussian.values b point
| _ -> not_implemented ()
end
|> Vector.relabel

8
ao/lib/basis.mli
View File

@ -5,10 +5,6 @@ open Particles
open Operators
open Linear_algebra
type basis =
| Unknown
| Gaussian of Basis_gaussian.t
type t
type ao = Ao_dim.t
@ -17,13 +13,13 @@ type ao = Ao_dim.t
val size : t -> int
(** Number of atomic orbitals in the AO basis set *)
val ao_basis : t -> basis
val ao_basis : t -> Basis_poly.t
(** One-electron basis set *)
val overlap : t -> (ao,ao) Matrix.t
(** Overlap matrix *)
val multipole : t -> (ao,ao) Matrix.t array
val multipole : t -> string -> (ao,ao) Matrix.t
(** Multipole matrices *)
val ortho : t -> (ao,'a) Matrix.t

3
ao/lib/basis_poly.ml Normal file
View File

@ -0,0 +1,3 @@
type t =
| Unknown
| Gaussian of Basis_gaussian.t

6
ao/lib/basis_poly.mli Normal file
View File

@ -0,0 +1,6 @@
(** Polymorphic type for the basis (Gaussian, Slater, etc) *)
type t =
| Unknown
| Gaussian of Basis_gaussian.t

2
common/lib/util.ml
View File

@ -43,6 +43,8 @@ let () =
Sys.set_signal Sys.sigint (Sys.Signal_handle f)
;;
let not_implemented () =
failwith "Not implemented"
let memo_float_of_int =
Array.init 64 float_of_int

3
common/lib/util.mli
View File

@ -41,6 +41,9 @@ val leadz : int64 -> int
(** {2 General functions} *)
val not_implemented : unit -> 'a
(** Fails with error message if some functionality is not implemented. *)
val fact : int -> float
(** Factorial function.
@raise Invalid_argument for negative arguments or arguments >100.

6
examples/dune Normal file
View File

@ -0,0 +1,6 @@
(executable
(name ex_integrals)
(libraries
qcaml
))

93
examples/ex_integrals.ml Normal file
View File

@ -0,0 +1,93 @@
(* [[file:~/QCaml/examples/ex_integrals.org::*Header][Header:1]] *)
module Command_line = Qcaml.Common.Command_line
module Util = Qcaml.Common.Util
open Qcaml.Linear_algebra
let () =
(* Header:1 ends here *)
(* [[file:~/QCaml/examples/ex_integrals.org::*Definition][Definition:1]] *)
let open Command_line in
begin
set_header_doc (Sys.argv.(0) ^ " - QuAcK command");
set_description_doc "Computes the one- and two-electron integrals on the Gaussian atomic basis set.";
set_specs
[ { short='b' ; long="basis" ; opt=Mandatory;
arg=With_arg "<string>";
doc="Name of the file containing the basis set"; } ;
{ short='x' ; long="xyz" ; opt=Mandatory;
arg=With_arg "<string>";
doc="Name of the file containing the nuclear coordinates in xyz format"; } ;
{ short='u' ; long="range-separation" ; opt=Optional;
arg=With_arg "<float>";
doc="Range-separation parameter."; } ;
]
end;
(* Definition:1 ends here *)
(* [[file:~/QCaml/examples/ex_integrals.org::*Interpretation][Interpretation:1]] *)
let basis_file = Util.of_some @@ Command_line.get "basis" in
let nuclei_file = Util.of_some @@ Command_line.get "xyz" in
let range_separation =
match Command_line.get "range-separation" with
| None -> None
| Some mu -> Some (float_of_string mu)
in
let operators =
match range_separation with
| None -> []
| Some mu -> [ Qcaml.Operators.Operator.of_range_separation mu ]
in
(* Interpretation:1 ends here *)
(* [[file:~/QCaml/examples/ex_integrals.org::*Computation][Computation:1]] *)
let nuclei =
Qcaml.Particles.Nuclei.of_xyz_file nuclei_file
in
(* Computation:1 ends here *)
(* [[file:~/QCaml/examples/ex_integrals.org::*Computation][Computation:2]] *)
let ao_basis =
Qcaml.Ao.Basis.of_nuclei_and_basis_filename ~kind:`Gaussian
~operators ~cartesian:true ~nuclei basis_file
in
(* Computation:2 ends here *)
(* [[file:~/QCaml/examples/ex_integrals.org::*Computation][Computation:3]] *)
let overlap = Qcaml.Ao.Basis.overlap ao_basis in
let eN_ints = Qcaml.Ao.Basis.eN_ints ao_basis in
let kin_ints = Qcaml.Ao.Basis.kin_ints ao_basis in
let multipole = Qcaml.Ao.Basis.multipole ao_basis in
let x_mat = multipole "x" in
let y_mat = multipole "y" in
let z_mat = multipole "z" in
(* Computation:3 ends here *)
(* [[file:~/QCaml/examples/ex_integrals.org::*Computation][Computation:4]] *)
let ee_ints = Qcaml.Ao.Basis.ee_ints ao_basis in
let lr_ints =
match range_separation with
| Some _mu -> Some (Qcaml.Ao.Basis.ee_lr_ints ao_basis)
| None -> None
in
(* Computation:4 ends here *)
(* [[file:~/QCaml/examples/ex_integrals.org::*Output][Output:1]] *)
Matrix.to_file ~filename:"overlap.dat" ~sym:true overlap;
Matrix.to_file ~filename:"eN.dat" ~sym:true eN_ints;
Matrix.to_file ~filename:"kinetic.dat" ~sym:true kin_ints;
Matrix.to_file ~filename:"x.dat" ~sym:true x_mat;
Matrix.to_file ~filename:"y.dat" ~sym:true y_mat;
Matrix.to_file ~filename:"z.dat" ~sym:true z_mat;
(* Output:1 ends here *)
(* [[file:~/QCaml/examples/ex_integrals.org::*Output][Output:2]] *)
Four_idx_storage.to_file ~filename:"eri.dat" ee_ints;
match lr_ints with
| Some integrals -> Four_idx_storage.to_file ~filename:"eri_lr.dat" integrals;
| None -> ()
(* Output:2 ends here *)

130
examples/ex_integrals.org Normal file
View File

@ -0,0 +1,130 @@
#+TITLE: Integrals
#+PROPERTY
In this example, we write a program that reads the geometry of a
molecule in in =xyz= format and a Gaussian atomic basis set in GAMESS
format. The output is a set of files containing the one- and two-
electron integrals.
* Header
#+BEGIN_SRC ocaml :comments link :exports code :tangle ex_integrals.ml
module Command_line = Qcaml.Common.Command_line
module Util = Qcaml.Common.Util
open Qcaml.Linear_algebra
let () =
#+END_SRC
* Command-line arguments
We use the =Command_line= module to define the following possible
arguments:
- =-b --basis= : The name of the file containing the basis set
- =-x --xyz= : The name of the file containing the atomic coordinates
- =-u --range-separation= : The value of $\mu$, the range-separation
parameter in range-separated DFT. If this option is not present,
no output file will be generated for the range-separated integrals.
** Definition
#+BEGIN_SRC ocaml :comments link :exports code :tangle ex_integrals.ml
let open Command_line in
begin
set_header_doc (Sys.argv.(0) ^ " - QuAcK command");
set_description_doc "Computes the one- and two-electron integrals on the Gaussian atomic basis set.";
set_specs
[ { short='b' ; long="basis" ; opt=Mandatory;
arg=With_arg "<string>";
doc="Name of the file containing the basis set"; } ;
{ short='x' ; long="xyz" ; opt=Mandatory;
arg=With_arg "<string>";
doc="Name of the file containing the nuclear coordinates in xyz format"; } ;
{ short='u' ; long="range-separation" ; opt=Optional;
arg=With_arg "<float>";
doc="Range-separation parameter."; } ;
]
end;
#+END_SRC
** Interpretation
#+BEGIN_SRC ocaml :comments link :exports code :tangle ex_integrals.ml
let basis_file = Util.of_some @@ Command_line.get "basis" in
let nuclei_file = Util.of_some @@ Command_line.get "xyz" in
let range_separation =
match Command_line.get "range-separation" with
| None -> None
| Some mu -> Some (float_of_string mu)
in
let operators =
match range_separation with
| None -> []
| Some mu -> [ Qcaml.Operators.Operator.of_range_separation mu ]
in
#+END_SRC
* Computation
We first read the =xyz= file to create a molecule:
#+BEGIN_SRC ocaml :comments link :exports code :tangle ex_integrals.ml
let nuclei =
Qcaml.Particles.Nuclei.of_xyz_file nuclei_file
in
#+END_SRC
Then we create an Gaussian AO basis using the atomic coordinates,
and we optionally introduce the range-separation parameter via the
=operators=:
#+BEGIN_SRC ocaml :comments link :exports code :tangle ex_integrals.ml
let ao_basis =
Qcaml.Ao.Basis.of_nuclei_and_basis_filename ~kind:`Gaussian
~operators ~cartesian:true ~nuclei basis_file
in
#+END_SRC
We compute the required one-electron integrals:
#+BEGIN_SRC ocaml :comments link :exports code :tangle ex_integrals.ml
let overlap = Qcaml.Ao.Basis.overlap ao_basis in
let eN_ints = Qcaml.Ao.Basis.eN_ints ao_basis in
let kin_ints = Qcaml.Ao.Basis.kin_ints ao_basis in
let multipole = Qcaml.Ao.Basis.multipole ao_basis in
let x_mat = multipole "x" in
let y_mat = multipole "y" in
let z_mat = multipole "z" in
#+END_SRC
and the two-electron integrals (1/r and long range):
#+BEGIN_SRC ocaml :comments link :exports code :tangle ex_integrals.ml
let ee_ints = Qcaml.Ao.Basis.ee_ints ao_basis in
let lr_ints =
match range_separation with
| Some _mu -> Some (Qcaml.Ao.Basis.ee_lr_ints ao_basis)
| None -> None
in
#+END_SRC
* Output
We write the one-electron integrals:
#+BEGIN_SRC ocaml :comments link :exports code :tangle ex_integrals.ml
Matrix.to_file ~filename:"overlap.dat" ~sym:true overlap;
Matrix.to_file ~filename:"eN.dat" ~sym:true eN_ints;
Matrix.to_file ~filename:"kinetic.dat" ~sym:true kin_ints;
Matrix.to_file ~filename:"x.dat" ~sym:true x_mat;
Matrix.to_file ~filename:"y.dat" ~sym:true y_mat;
Matrix.to_file ~filename:"z.dat" ~sym:true z_mat;
#+END_SRC
and the the two-electron integrals:
#+BEGIN_SRC ocaml :comments link :exports code :tangle ex_integrals.ml
Four_idx_storage.to_file ~filename:"eri.dat" ee_ints;
match lr_ints with
| Some integrals -> Four_idx_storage.to_file ~filename:"eri_lr.dat" integrals;
| None -> ()
#+END_SRC

15
gaussian_integrals/lib/electron_nucleus.ml
View File

@ -91,21 +91,6 @@ let of_basis_nuclei ~basis nuclei =
Matrix.detri_inplace eni_array;
eni_array
let to_file ~filename eni_array =
let n = Matrix.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 j value;
done;
done;
close_out oc
let of_basis _basis =
invalid_arg "of_basis_nuclei should be called for NucInt"

17
gaussian_integrals/lib/kinetic.ml
View File

@ -166,20 +166,3 @@ let of_basis basis =
let of_basis_pair _first_basis _second_basis =
failwith "Not implemented"
(** Write all kinetic integrals to a file *)
let to_file ~filename kinetic =
let oc = open_out filename in
let n =
Matrix.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

3
gaussian_integrals/lib/matrix_on_basis.mli
View File

@ -10,6 +10,3 @@ val of_basis : Basis.t -> t
val of_basis_pair : Basis.t -> Basis.t -> t
(** Build from a pair of {!Basis.t}. *)
val to_file : filename:string -> t -> unit
(** Write the integrals in a file. *)

33
gaussian_integrals/lib/multipole.ml
View File

@ -16,22 +16,23 @@ module Csp = Contracted_shell_pair
module Po = Powers
module Psp = Primitive_shell_pair
let matrix_x t = t.(0)
let matrix_y t = t.(1)
let matrix_z t = t.(2)
let matrix_x2 t = t.(3)
let matrix_y2 t = t.(4)
let matrix_z2 t = t.(5)
let matrix_xy t = t.(6)
let matrix_xz t = t.(8)
let matrix_yz t = t.(7)
let matrix_x3 t = t.(9)
let matrix_y3 t = t.(10)
let matrix_z3 t = t.(11)
let matrix_x4 t = t.(12)
let matrix_y4 t = t.(13)
let matrix_z4 t = t.(14)
let matrix t = function
| "x" -> t.(0)
| "y" -> t.(1)
| "z" -> t.(2)
| "x2" -> t.(3)
| "y2" -> t.(4)
| "z2" -> t.(5)
| "xy" -> t.(6)
| "xz" -> t.(8)
| "yz" -> t.(7)
| "x3" -> t.(9)
| "y3" -> t.(10)
| "z3" -> t.(11)
| "x4" -> t.(12)
| "y4" -> t.(13)
| "z4" -> t.(14)
| _ -> Util.not_implemented ()

48
gaussian_integrals/lib/multipole.mli
View File

@ -14,50 +14,10 @@ open Gaussian
type t = (Basis.t, Basis.t) Matrix.t array
val matrix_x : t -> (Basis.t, Basis.t) Matrix.t
(** {% $$ \langle \chi_i | x | \chi_j \rangle $$ %} *)
val matrix_y : t -> (Basis.t, Basis.t) Matrix.t
(** {% $$ \langle \chi_i | y | \chi_j \rangle $$ %} *)
val matrix_z : t -> (Basis.t, Basis.t) Matrix.t
(** {% $$ \langle \chi_i | z | \chi_j \rangle $$ %} *)
val matrix_x2 : t -> (Basis.t, Basis.t) Matrix.t
(** {% $$ \langle \chi_i | x^2 | \chi_j \rangle $$ %} *)
val matrix_xy : t -> (Basis.t, Basis.t) Matrix.t
(** {% $$ \langle \chi_i | xy | \chi_j \rangle $$ %} *)
val matrix_yz : t -> (Basis.t, Basis.t) Matrix.t
(** {% $$ \langle \chi_i | yz | \chi_j \rangle $$ %} *)
val matrix_xz : t -> (Basis.t, Basis.t) Matrix.t
(** {% $$ \langle \chi_i | xz | \chi_j \rangle $$ %} *)
val matrix_y2 : t -> (Basis.t, Basis.t) Matrix.t
(** {% $$ \langle \chi_i | y^2 | \chi_j \rangle $$ %} *)
val matrix_z2 : t -> (Basis.t, Basis.t) Matrix.t
(** {% $$ \langle \chi_i | z^2 | \chi_j \rangle $$ %} *)
val matrix_x3 : t -> (Basis.t, Basis.t) Matrix.t
(** {% $$ \langle \chi_i | x^3 | \chi_j \rangle $$ %} *)
val matrix_y3 : t -> (Basis.t, Basis.t) Matrix.t
(** {% $$ \langle \chi_i | y^3 | \chi_j \rangle $$ %} *)
val matrix_z3 : t -> (Basis.t, Basis.t) Matrix.t
(** {% $$ \langle \chi_i | z^3 | \chi_j \rangle $$ %} *)
val matrix_x4 : t -> (Basis.t, Basis.t) Matrix.t
(** {% $$ \langle \chi_i | x^4 | \chi_j \rangle $$ %} *)
val matrix_y4 : t -> (Basis.t, Basis.t) Matrix.t
(** {% $$ \langle \chi_i | y^4 | \chi_j \rangle $$ %} *)
val matrix_z4 : t -> (Basis.t, Basis.t) Matrix.t
(** {% $$ \langle \chi_i | z^4 | \chi_j \rangle $$ %} *)
val matrix : t -> string -> (Basis.t, Basis.t) Matrix.t
(** Returns the requested matrix. Choices:
[x, y, z, x2, xy, yz, xz, y2, z2, z3, y3, z3, x4, y4, z4]
*)
val of_basis : Basis.t -> t

18
gaussian_integrals/lib/overlap.ml
View File

@ -174,21 +174,3 @@ let of_basis_pair first_basis second_basis =
(** Write all overlap integrals to a file *)
let to_file ~filename overlap =
let oc = open_out filename in
let n =
Matrix.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

29
linear_algebra/lib/matrix.ml
View File

@ -372,8 +372,37 @@ let qr a =
let q = result in
q, r
let (%:) t (i,j) = t.{i,j}
let to_file ~filename ?(sym=false) ?(cutoff=0.) t =
let oc = open_out filename in
let n = Mat.dim1 t in
let m = Mat.dim2 t in
if sym then
for j=1 to n do
for i=1 to j do
if (abs_float (t.{i,j}) > cutoff) then
Printf.fprintf oc "%4d %4d %20.12e\n" i j (t.{i,j})
done;
done
else
for j=1 to n do
for i=1 to m do
if (abs_float (t.{i,j}) > cutoff) then
Printf.fprintf oc "%4d %4d %20.12e\n" i j (t.{i,j})
done;
done;
close_out oc
let set t i j v = t.{i,j} <- v

3
linear_algebra/lib/matrix.mli
View File

@ -319,5 +319,8 @@ val (%:) : ('a,'b) t -> (int*int) -> float
val set : ('a,'b) t -> int -> int -> float -> unit
(** [set t i j v] sets the (i,j)-th element to v *)
val to_file : filename:string -> ?sym:bool -> ?cutoff:float -> ('a,'b) t -> unit
(** Write the matrix to a file. *)
val pp : Format.formatter -> ('a,'b) t -> unit