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

Added AO basis tests

This commit is contained in:
Anthony Scemama 2020-10-03 00:54:17 +02:00
parent 7de85e6a44
commit 2881d69bb6
8 changed files with 3280 additions and 3607 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

215
ao_basis/lib/ao_basis_gaussian.ml
View File

@ -1,21 +1,23 @@
module GB = Qcaml_gaussian_basis
module GI = Qcaml_gaussian_integrals
open Qcaml_linear_algebra
open Qcaml_gaussian_basis
open Qcaml_gaussian_integrals
open Qcaml_operators
type t =
{
basis : GB.Basis.t ;
overlap : GI.Overlap.t lazy_t;
multipole : GI.Multipole.t lazy_t;
ortho : GI.Orthonormalization.t lazy_t;
eN_ints : GI.Electron_nucleus.t lazy_t;
kin_ints : GI.Kinetic.t lazy_t;
ee_ints : GI.Eri.t lazy_t;
ee_lr_ints : GI.Eri_long_range.t lazy_t;
f12_ints : GI.F12.t lazy_t;
f12_over_r12_ints : GI.Screened_eri.t lazy_t;
basis : Basis.t ;
overlap : Overlap.t lazy_t;
multipole : Multipole.t lazy_t;
ortho : Orthonormalization.t lazy_t;
eN_ints : Electron_nucleus.t lazy_t;
kin_ints : Kinetic.t lazy_t;
ee_ints : Eri.t lazy_t;
ee_lr_ints : Eri_long_range.t lazy_t;
f12_ints : F12.t lazy_t;
f12_over_r12_ints : Screened_eri.t lazy_t;
cartesian : bool ;
}
(*
let basis t = t.basis
let overlap t = Lazy.force t.overlap
let multipole t = Lazy.force t.multipole
@ -28,162 +30,71 @@ let f12_ints t = Lazy.force t.f12_ints
let f12_over_r12_ints t = Lazy.force t.f12_over_r12_ints
let cartesian t = t.cartesian
module Cs = ContractedShell
module Cs = Contracted_shell
let values t point =
let result = Vec.create (Basis.size t.basis) in
let result = Vector.create (Basis.size t.basis) in
let resultx = Vector.to_bigarray_inplace result in
Array.iter (fun shell ->
Cs.values shell point
|> Array.iteri
(fun i_c value ->
let i = Cs.index shell + i_c + 1 in
result.{i} <- value)
resultx.{i} <- value)
) (Basis.contracted_shells t.basis);
result
let make ~cartesian ~basis ?f12 nuclei =
let make ~basis ?(operators=[]) ?(cartesian=false) nuclei =
let overlap =
lazy (
Overlap.of_basis basis
) in
let overlap = lazy (
Overlap.of_basis basis
) in
let ortho =
lazy (
Orthonormalization.make ~cartesian ~basis (Lazy.force overlap)
) in
let ortho = lazy (
Lazy.force overlap
|> Orthonormalization.make ~cartesian ~basis
) in
let eN_ints = lazy (
Electron_nucleus.of_basis_nuclei ~basis nuclei
) in
let kin_ints = lazy (
Kinetic.of_basis basis
) in
let ee_ints = lazy (
Eri.of_basis basis
) in
let rec get_f12 = function
| (Operator.F12 _ as f) :: _ -> f
| [] -> failwith "Missing F12 operator"
| _ :: rest -> get_f12 rest
in
let eN_ints =
lazy (
Electron_nucleus.of_basis_nuclei ~basis nuclei
) in
let rec get_rs = function
| (Operator.Range_sep _ as r) :: _ -> r
| [] -> failwith "Missing range-separation operator"
| _ :: rest -> get_rs rest
in
let kin_ints =
lazy (
Kinetic.of_basis basis
) in
let ee_lr_ints = lazy (
Eri_long_range.of_basis basis~operator:(get_rs operators)
) in
let f12_ints = lazy (
F12.of_basis basis ~operator:(get_f12 operators)
) in
let f12_over_r12_ints = lazy (
Screened_eri.of_basis basis ~operator:(get_f12 operators)
) in
let ee_ints =
lazy (
ERI.of_basis basis
) in
let ee_lr_ints =
lazy (
ERI_lr.of_basis basis
) in
let f12_ints =
lazy (
F12.of_basis basis
) in
let f12_over_r12_ints =
lazy (
ScreenedERI.of_basis basis
) in
let multipole =
lazy (
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 ;
}
let test_case name t =
let check_matrix title a r =
let a = Mat.to_array a in
Mat.to_array r
|> Array.iteri (fun i x ->
let message =
Printf.sprintf "%s line %d" title i
in
Alcotest.(check (array (float 1.e-10))) message a.(i) x
)
in
let check_eri a r =
let f { ERI.i_r1 ; j_r2 ; k_r1 ; l_r2 ; value } =
(i_r1, (j_r2, (k_r1, (l_r2, value))))
in
let a = ERI.to_list a |> List.rev_map f |> List.rev in
let r = ERI.to_list r |> List.rev_map f |> List.rev in
Alcotest.(check (list (pair int (pair int (pair int (pair int (float 1.e-10))))))) "ERI" a r
in
let check_eri_lr a r =
let f { ERI_lr.i_r1 ; j_r2 ; k_r1 ; l_r2 ; value } =
(i_r1, (j_r2, (k_r1, (l_r2, value))))
in
let a = ERI_lr.to_list a |> List.rev_map f |> List.rev in
let r = ERI_lr.to_list r |> List.rev_map f |> List.rev in
Alcotest.(check (list (pair int (pair int (pair int (pair int (float 1.e-10))))))) "ERI_lr" a r
in
let test_overlap () =
let reference =
sym_matrix_of_file ("test_files/"^name^"_overlap.ref")
in
let overlap =
Lazy.force t.overlap |> Overlap.matrix
in
check_matrix "Overlap" overlap reference
in
let test_eN_ints () =
let reference =
sym_matrix_of_file ("test_files/"^name^"_nuc.ref")
in
let eN_ints =
Lazy.force t.eN_ints |> NucInt.matrix
in
check_matrix "eN_ints" eN_ints reference
in
let test_kin_ints () =
let reference =
sym_matrix_of_file ("test_files/"^name^"_kin.ref")
in
let kin_ints =
Lazy.force t.kin_ints |> KinInt.matrix
in
check_matrix "kin_ints" kin_ints reference
in
let test_ee_ints () =
let reference =
ERI.of_file ("test_files/"^name^"_eri.ref") ~sparsity:`Dense ~size:(Basis.size t.basis)
in
let ee_ints =
Lazy.force t.ee_ints
in
check_eri ee_ints reference
;
in
let test_ee_lr_ints () =
let reference =
ERI_lr.of_file ("test_files/"^name^"_eri_lr.ref") ~sparsity:`Dense
~size:(Basis.size t.basis)
in
let ee_lr_ints =
Lazy.force t.ee_lr_ints
in
check_eri_lr ee_lr_ints reference
in
[
"Overlap", `Quick, test_overlap;
"eN_ints", `Quick, test_eN_ints;
"kin_ints", `Quick, test_kin_ints;
"ee_ints", `Quick, test_ee_ints;
"ee_lr_ints", `Quick, test_ee_lr_ints;
]
*)
ee_lr_ints ; f12_ints ; f12_over_r12_ints ; cartesian }

53
ao_basis/lib/ao_basis_gaussian.mli
View File

@ -1,24 +1,27 @@
(** Data structure for Atomic Orbitals. *)
module GB = Qcaml_gaussian_basis
module GI = Qcaml_gaussian_integrals
open Qcaml_common
open Qcaml_particles
open Qcaml_linear_algebra
open Qcaml_gaussian_basis
open Qcaml_gaussian_integrals
open Qcaml_operators
type t =
{
basis : GB.Basis.t ;
overlap : GI.Overlap.t lazy_t;
multipole : GI.Multipole.t lazy_t;
ortho : GI.Orthonormalization.t lazy_t;
eN_ints : GI.Electron_nucleus.t lazy_t;
kin_ints : GI.Kinetic.t lazy_t;
ee_ints : GI.Eri.t lazy_t;
ee_lr_ints : GI.Eri_long_range.t lazy_t;
f12_ints : GI.F12.t lazy_t;
f12_over_r12_ints : GI.Screened_eri.t lazy_t;
basis : Basis.t ;
overlap : Overlap.t lazy_t;
multipole : Multipole.t lazy_t;
ortho : Orthonormalization.t lazy_t;
eN_ints : Electron_nucleus.t lazy_t;
kin_ints : Kinetic.t lazy_t;
ee_ints : Eri.t lazy_t;
ee_lr_ints : Eri_long_range.t lazy_t;
f12_ints : F12.t lazy_t;
f12_over_r12_ints : Screened_eri.t lazy_t;
cartesian : bool ;
}
(*
(** {1 Accessors} *)
val basis : t -> Basis.t
@ -33,37 +36,35 @@ val multipole : t -> Multipole.t
val ortho : t -> Orthonormalization.t
(** Orthonormalization matrix of the overlap *)
val eN_ints : t -> NucInt.t
val eN_ints : t -> Electron_nucleus.t
(** Electron-nucleus potential integrals *)
val ee_ints : t -> ERI.t
val ee_ints : t -> Eri.t
(** Electron-electron potential integrals *)
val ee_lr_ints : t -> ERI_lr.t
val ee_lr_ints : t -> Eri_long_range.t
(** Electron-electron long-range potential integrals *)
val f12_ints : t -> F12.t
(** Electron-electron potential integrals *)
val kin_ints : t -> KinInt.t
val f12_over_r12_ints : t -> Screened_eri.t
(** Electron-electron potential integrals *)
val kin_ints : t -> Kinetic.t
(** Kinetic energy integrals *)
val cartesian : t -> bool
(** If true, use cartesian Gaussians (6d, 10f, ...) *)
val values : t -> Coordinate.t -> Vec.t
val values : t -> Coordinate.t -> Basis.t Vector.t
(** Values of the AOs evaluated at a given point *)
(** {1 Creators} *)
val make : cartesian:bool -> basis:Basis.t -> ?f12:F12factor.t -> Nuclei.t -> t
val make : basis:Basis.t -> ?operators:Operator.t list -> ?cartesian:bool
-> Nuclei.t -> t
(** Creates the data structure for atomic orbitals from a {Basis.t} and the
molecular geometry {Nuclei.t} *)
(** {2 Tests} *)
val test_case : string -> t -> unit Alcotest.test_case list
*)

1
ao_basis/lib/dune
View File

@ -5,6 +5,7 @@
(public_name qcaml.ao_basis)
(libraries
qcaml.common
qcaml.particles
qcaml.gaussian_basis
qcaml.gaussian_integrals
qcaml.operators

11
test/caffeine.xyz
View File

@ -1,8 +1,8 @@
24
Caffeine test file
H -3.3804130 -1.1272367 0.5733036
N 0.9668296 -1.0737425 -0.8198227
C 0.0567293 0.8527195 0.3923156
Caffeine
1 -3.3804130 -1.1272367 0.5733036
7 0.9668296 -1.0737425 -0.8198227
6 0.0567293 0.8527195 0.3923156
N -1.3751742 -1.0212243 -0.0570552
C -1.2615018 0.2590713 0.5234135
C -0.3068337 -1.6836331 -0.7169344
@ -24,6 +24,3 @@ N 2.2861252 0.9968439 -0.2440298
H -0.1687028 4.0436553 0.9301094
H 0.3535322 3.2979060 2.5177747
H -1.2074498 2.7537592 1.7203047

6602
test/cc-pvdz
View File

File diff suppressed because it is too large Load Diff

2
test/dune
View File

@ -6,6 +6,8 @@
test_linear_algebra
test_particles
test_gaussian_basis
test_gaussian_integrals
test_ao_basis
))
(alias

1
test/run_tests.ml
View File

@ -9,6 +9,7 @@ let test_suites: unit Alcotest.test list = [
"Particles.Nuclei", Test_particles.Nuclei.tests;
"Particles.Electrons", Test_particles.Electrons.tests;
"Gaussian_basis.General_basis", Test_gaussian_basis.General_basis.tests;
"Ao_basis.Ao_basis_gaussian", Test_ao_basis.Ao_basis_gaussian.tests;
]
let () = Alcotest.run "QCaml" test_suites

2
test/water.xyz
View File

@ -3,5 +3,3 @@ Water
O 0. 0. 0.
H -0.756950272703377558 0. -0.585882234512562827
H 0.756950272703377558 0. -0.585882234512562827