2021-01-04 23:15:27 +01:00
|
|
|
#+begin_src elisp tangle: no :results none :exports none
|
|
|
|
(setq pwd (file-name-directory buffer-file-name))
|
|
|
|
(setq name (file-name-nondirectory (substring buffer-file-name 0 -4)))
|
|
|
|
(setq lib (concat pwd "lib/"))
|
|
|
|
(setq testdir (concat pwd "test/"))
|
|
|
|
(setq mli (concat lib name ".mli"))
|
|
|
|
(setq ml (concat lib name ".ml"))
|
|
|
|
(setq test-ml (concat testdir name ".ml"))
|
|
|
|
(org-babel-tangle)
|
|
|
|
#+end_src
|
|
|
|
|
|
|
|
* Gaussian basis
|
|
|
|
:PROPERTIES:
|
|
|
|
:header-args: :noweb yes :comments both
|
|
|
|
:END:
|
|
|
|
|
|
|
|
Data structure for Gaussian Atomic Orbitals:
|
|
|
|
|
|
|
|
$$
|
2021-01-04 23:52:36 +01:00
|
|
|
\chi_i(\mathbf{r}) = P_i(\mathbf{r}) \sum_k c_k \exp\left( -\alpha_k (\mathbf{r-R_A})^2 \right)
|
2021-01-04 23:15:27 +01:00
|
|
|
$$
|
|
|
|
|
2021-01-04 23:52:36 +01:00
|
|
|
where the polynomial $P_Ai and the Gaussian part are both centered on
|
2021-01-04 23:15:27 +01:00
|
|
|
nucleus $A$.
|
|
|
|
|
|
|
|
** Type
|
|
|
|
|
|
|
|
#+NAME: types
|
|
|
|
#+begin_src ocaml :tangle (eval mli)
|
|
|
|
open Common
|
|
|
|
open Particles
|
|
|
|
open Linear_algebra
|
|
|
|
open Gaussian_integrals
|
|
|
|
open Operators
|
|
|
|
|
|
|
|
type t
|
|
|
|
#+end_src
|
|
|
|
|
|
|
|
#+begin_src ocaml :tangle (eval ml) :exports none
|
|
|
|
open Linear_algebra
|
|
|
|
open Gaussian
|
|
|
|
open Gaussian_integrals
|
|
|
|
open Operators
|
|
|
|
|
|
|
|
module Basis = Gaussian.Basis
|
|
|
|
|
|
|
|
type 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 ;
|
|
|
|
}
|
|
|
|
#+end_src
|
|
|
|
|
|
|
|
** Access
|
|
|
|
|
|
|
|
#+begin_src ocaml :tangle (eval mli)
|
|
|
|
val basis : t -> Gaussian.Basis.t
|
|
|
|
val cartesian : t -> bool
|
|
|
|
val ee_ints : t -> Eri.t
|
|
|
|
val ee_lr_ints : t -> Eri_long_range.t
|
|
|
|
val eN_ints : t -> Electron_nucleus.t
|
|
|
|
val f12_ints : t -> F12.t
|
|
|
|
val f12_over_r12_ints : t -> Screened_eri.t
|
|
|
|
val kin_ints : t -> Kinetic.t
|
|
|
|
val multipole : t -> Multipole.t
|
|
|
|
val ortho : t -> Orthonormalization.t
|
|
|
|
val overlap : t -> Overlap.t
|
|
|
|
val size : t -> int
|
|
|
|
#+end_src
|
|
|
|
|
|
|
|
| ~basis~ | One-electron basis set |
|
|
|
|
| ~cartesian~ | If true, use cartesian Gaussians (6d, 10f, ...) |
|
|
|
|
| ~ee_ints~ | Electron-electron potential integrals |
|
|
|
|
| ~ee_lr_ints~ | Electron-electron long-range potential integrals |
|
|
|
|
| ~eN_ints~ | Electron-nucleus potential integrals |
|
|
|
|
| ~f12_ints~ | Electron-electron potential integrals |
|
|
|
|
| ~f12_over_r12_ints~ | Electron-electron potential integrals |
|
|
|
|
| ~kin_ints~ | Kinetic energy integrals |
|
|
|
|
| ~multipole~ | Multipole matrices |
|
|
|
|
| ~ortho~ | Orthonormalization matrix of the overlap |
|
|
|
|
| ~overlap~ | Overlap matrix |
|
|
|
|
| ~size~ | Number of atomic orbitals |
|
|
|
|
|
|
|
|
#+begin_src ocaml :tangle (eval ml) :exports none
|
|
|
|
let basis t = t.basis
|
|
|
|
let cartesian t = t.cartesian
|
|
|
|
let ee_ints t = Lazy.force t.ee_ints
|
|
|
|
let ee_lr_ints t = Lazy.force t.ee_lr_ints
|
|
|
|
let eN_ints t = Lazy.force t.eN_ints
|
|
|
|
let f12_ints t = Lazy.force t.f12_ints
|
|
|
|
let f12_over_r12_ints t = Lazy.force t.f12_over_r12_ints
|
|
|
|
let kin_ints t = Lazy.force t.kin_ints
|
|
|
|
let multipole t = Lazy.force t.multipole
|
|
|
|
let ortho t = Lazy.force t.ortho
|
|
|
|
let overlap t = Lazy.force t.overlap
|
|
|
|
let size t = Matrix.dim1 (Lazy.force t.overlap)
|
|
|
|
#+end_src
|
|
|
|
|
|
|
|
** Computation
|
|
|
|
|
|
|
|
#+begin_src ocaml :tangle (eval mli)
|
|
|
|
val values : t -> Coordinate.t -> Gaussian.Basis.t Vector.t
|
|
|
|
#+end_src
|
|
|
|
|
|
|
|
| ~values~ | Returns the values of all the AOs evaluated at a given point |
|
|
|
|
|
|
|
|
#+begin_src ocaml :tangle (eval ml) :exports none
|
|
|
|
module Cs = Contracted_shell
|
|
|
|
|
|
|
|
let values t point =
|
|
|
|
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
|
|
|
|
resultx.{i} <- value)
|
|
|
|
) (Basis.contracted_shells t.basis);
|
|
|
|
result
|
|
|
|
#+end_src
|
|
|
|
|
|
|
|
** Creation
|
|
|
|
|
|
|
|
#+begin_src ocaml :tangle (eval mli)
|
|
|
|
val make : basis:Gaussian.Basis.t ->
|
|
|
|
?operators:Operator.t list ->
|
|
|
|
?cartesian:bool ->
|
|
|
|
Nuclei.t -> t
|
|
|
|
#+end_src
|
|
|
|
|
|
|
|
Creates the data structure for atomic orbitals from a Gaussian basis and the
|
|
|
|
molecular geometry ~Nuclei.t~.
|
|
|
|
|
|
|
|
Defaults:
|
|
|
|
- ~operators~ : ~[]~
|
|
|
|
- ~cartesian~ : ~false~
|
|
|
|
|
|
|
|
#+begin_example
|
|
|
|
let b = Ao.Basis_gaussian.make ~basis nuclei ;;
|
|
|
|
val b : Ao.Basis_gaussian.t = Gaussian Basis, spherical, 15 AOs
|
|
|
|
#+end_example
|
|
|
|
|
|
|
|
#+begin_src ocaml :tangle (eval ml) :exports none
|
|
|
|
let make ~basis ?(operators=[]) ?(cartesian=false) nuclei =
|
|
|
|
|
|
|
|
let overlap = lazy (
|
|
|
|
Overlap.of_basis basis
|
|
|
|
) 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 rec get_rs = function
|
|
|
|
| (Operator.Range_sep _ as r) :: _ -> r
|
|
|
|
| [] -> failwith "Missing range-separation operator"
|
|
|
|
| _ :: rest -> get_rs rest
|
|
|
|
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 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 }
|
|
|
|
#+end_src
|
|
|
|
|
|
|
|
|
|
|
|
** Printers
|
|
|
|
|
|
|
|
#+begin_src ocaml :tangle (eval mli)
|
|
|
|
val pp : Format.formatter -> t -> unit
|
|
|
|
#+end_src
|
|
|
|
|
|
|
|
#+begin_src ocaml :tangle (eval ml) :exports none
|
|
|
|
let pp ppf t =
|
|
|
|
let cart = if t.cartesian then "cartesian" else "spherical" in
|
|
|
|
let nao = size t in
|
|
|
|
Format.fprintf ppf "@[@[Gaussian Basis@], @[%s@], @[%d AOs@]@]"
|
|
|
|
cart nao
|
|
|
|
#+end_src
|
|
|
|
|