#+TITLE: Integrals #+PROPERTY In this example, we write a program that reads the geometry of a molecule in =xyz= format and a Gaussian atomic basis set in GAMESS format. The output is written in a trexio file. * 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 #+END_SRC #+RESULTS: : module Command_line = Qcaml.Common.Command_line : module Util = Qcaml.Common.Util #+BEGIN_SRC ocaml :comments link :exports code :tangle ex_integrals.ml 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)); 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 ""; doc="Name of the file containing the basis set"; } ; { short='x' ; long="xyz" ; opt=Mandatory; arg=With_arg ""; doc="Name of the file containing the nuclear coordinates in xyz format"; } ; { short='u' ; long="range-separation" ; opt=Optional; arg=With_arg ""; 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 * Interactive test :noexport: #+begin_src ocaml :results drawer :session :cache no :exports none #require "qcaml.top" ;; #require "trexio" ;; open Qcaml ;; #+end_src #+RESULTS: :results: :end: #+begin_src ocaml :results drawer :session :cache no :exports none let nuclei_file = "/dev/shm/f2.xyz" let nuclei = Qcaml.Particles.Nuclei.of_xyz_file nuclei_file #+end_src #+RESULTS: :results: val nuclei_file : string = "/dev/shm/f2.xyz" val nuclei : Qcaml.Particles.Nuclei.t = Nuclear Coordinates (Angstrom) ------------------------------ ----------------------------------------------------------------------- Center Atomic Element Coordinates (Angstroms) Number X Y Z ----------------------------------------------------------------------- 1 9 F 0.000000 0.000000 0.000000 2 9 F 0.000000 0.000000 1.411900 ----------------------------------------------------------------------- :end: #+begin_src ocaml :results drawer :session :cache no :exports none let trexio_file = Trexio.open_file "test.trexio" 'w' Trexio.HDF5 #+end_src #+RESULTS: :results: val trexio_file : Trexio.trexio_file = :end: #+begin_src ocaml :results drawer :session :cache no :exports none Qcaml.Particles.Nuclei.to_trexio trexio_file nuclei;; Trexio.close_file trexio_file;; #+end_src #+RESULTS: :results: - : unit = () :end: #+begin_src ocaml :results drawer :session :cache no :exports none let basis_file = "/home/scemama/qp2/data/basis/cc-pvdz";; let ao_basis = Qcaml.Ao.Basis.of_nuclei_and_basis_filename ~nuclei basis_file #+end_src #+RESULTS: :results: val ao_basis : Qcaml.Ao.Basis.t = Gaussian Basis, spherical, 30 AOs :end: #+begin_src ocaml :results drawer :session :cache no :exports none Trexio.close_file trexio_file;; #+end_src