From 9e5933f0d0d927a6724ce583247910d82be2b44f Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 13 Jun 2018 17:49:58 +0200 Subject: [PATCH] Introduced AOBasis --- .merlin | 2 ++ Basis/AOBasis.ml | 43 +++++++++++++++++++++++++++++++++++++++++++ Basis/AOBasis.mli | 19 +++++++++++++++++++ Basis/KinInt.mli | 17 +++++++++++++++++ Basis/NucInt.ml | 2 +- Basis/NucInt.mli | 15 +++++++++++++++ SCF/Fock.ml | 11 ++++++----- SCF/Guess.ml | 47 ++++++++++++++++++++++++----------------------- SCF/Guess.mli | 2 +- SCF/RHF.ml | 25 ++++++++++++++++--------- Simulation.ml | 36 ++++++++++++++++-------------------- run_integrals.ml | 9 +++++---- 12 files changed, 165 insertions(+), 63 deletions(-) create mode 100644 Basis/AOBasis.ml create mode 100644 Basis/AOBasis.mli create mode 100644 Basis/KinInt.mli create mode 100644 Basis/NucInt.mli diff --git a/.merlin b/.merlin index f2649f0..83ecc10 100644 --- a/.merlin +++ b/.merlin @@ -1,7 +1,9 @@ PKG str unix bigarray lacaml S . S Nuclei +S CI S Utils S Basis S SCF +S MOBasis B _build/** diff --git a/Basis/AOBasis.ml b/Basis/AOBasis.ml new file mode 100644 index 0000000..e018ab0 --- /dev/null +++ b/Basis/AOBasis.ml @@ -0,0 +1,43 @@ +open Lacaml.D + +type t = +{ + basis : Basis.t ; + overlap : Overlap.t lazy_t; + ortho : Orthonormalization.t lazy_t; + eN_ints : NucInt.t lazy_t; + ee_ints : ERI.t lazy_t; + kin_ints : KinInt.t lazy_t; + cartesian : bool; +} + +let make ~cartesian ~basis nuclei = + + let overlap = lazy ( + Overlap.of_basis basis + ) in + + let ortho = lazy ( + Orthonormalization.make ~cartesian ~basis (Lazy.force overlap) + ) in + + let eN_ints = lazy ( + NucInt.of_basis_nuclei ~basis nuclei + ) in + + let kin_ints = lazy ( + KinInt.of_basis basis + ) in + + let ee_ints = lazy ( + ERI.of_basis basis + ) in + + { basis ; overlap ; ortho ; eN_ints ; kin_ints ; ee_ints ; + cartesian ; + } + + + + + diff --git a/Basis/AOBasis.mli b/Basis/AOBasis.mli new file mode 100644 index 0000000..1369dfb --- /dev/null +++ b/Basis/AOBasis.mli @@ -0,0 +1,19 @@ +(** Data structure for Atomic Orbitals. *) + +type t = +{ + basis : Basis.t ; (* One-electron basis set *) + overlap : Overlap.t lazy_t; (* Overlap matrix *) + ortho : Orthonormalization.t lazy_t; (* Orthonormalization matrix of the overlap *) + eN_ints : NucInt.t lazy_t; (* Electron-nucleus potential integrals *) + ee_ints : ERI.t lazy_t; (* Electron-electron potential integrals *) + kin_ints : KinInt.t lazy_t; (* Kinetic energy integrals *) + cartesian : bool ; (* If true, use cartesian Gaussians (6d, 10f, ...) *) +} + + +val make : cartesian:bool -> basis:Basis.t -> Nuclei.t -> t +(** Creates the data structure for atomic orbitals from a {Basis.t} and the + molecular geometry {Nuclei.t} *) + + diff --git a/Basis/KinInt.mli b/Basis/KinInt.mli new file mode 100644 index 0000000..dbccc0b --- /dev/null +++ b/Basis/KinInt.mli @@ -0,0 +1,17 @@ +(** Electronic kinetic energy integrals, expressed as a matrix in a {Basis.t}. + +{% +$$ +T_{ij} = \left \langle \chi_i \left| -\frac{1}{2} \Delta \right| \chi_j \right \rangle +$$ +%} +*) + +type t = Lacaml.D.Mat.t + +val of_basis : Basis.t -> t +(** Build from a {Basis.t}. *) + +val to_file : filename:string -> t -> unit +(** Write the integrals in a file. *) + diff --git a/Basis/NucInt.ml b/Basis/NucInt.ml index 9b555e7..dc92009 100644 --- a/Basis/NucInt.ml +++ b/Basis/NucInt.ml @@ -41,7 +41,7 @@ let cutoff2 = cutoff *. cutoff exception NullIntegral -let of_basis_nuclei basis nuclei = +let of_basis_nuclei ~basis nuclei = let to_powers x = let open Zkey in match to_powers x with diff --git a/Basis/NucInt.mli b/Basis/NucInt.mli new file mode 100644 index 0000000..28cab01 --- /dev/null +++ b/Basis/NucInt.mli @@ -0,0 +1,15 @@ +(** Electron-Nucleus attractive potential integrals, expressed as a matrix in a {Basis.t}. + +{% $$ +V_{ij} = \left \langle \chi_i \left| \sum_A \frac{-Z_A}{ | \mathbf{r} - \mathbf{R}_A |} \right| \chi_j \right \rangle +$$ %} + +*) + +type t = Lacaml.D.Mat.t + +val of_basis_nuclei : basis:Basis.t -> Nuclei.t -> t +(** Build from a {Basis.t} and the nuclei (coordinates and charges). *) + +val to_file : filename:string -> t -> unit +(** Write the integrals in a file. *) diff --git a/SCF/Fock.ml b/SCF/Fock.ml index 20a5d8e..b4c4d9b 100644 --- a/SCF/Fock.ml +++ b/SCF/Fock.ml @@ -3,7 +3,6 @@ open Simulation open Constants - type t = { fock : Mat.t ; @@ -12,11 +11,13 @@ type t = exchange : Mat.t ; } -let make ~density simulation = +module Ao = AOBasis + +let make ~density ao_basis = let m_P = density - and m_T = Lazy.force simulation.kin_ints - and m_V = Lazy.force simulation.eN_ints - and m_G = Lazy.force simulation.ee_ints + and m_T = Lazy.force ao_basis.Ao.kin_ints + and m_V = Lazy.force ao_basis.Ao.eN_ints + and m_G = Lazy.force ao_basis.Ao.ee_ints in let nBas = Mat.dim1 m_T in diff --git a/SCF/Guess.ml b/SCF/Guess.ml index 653f193..3b5c4a4 100644 --- a/SCF/Guess.ml +++ b/SCF/Guess.ml @@ -7,43 +7,44 @@ type guess = type t = guess -module Si = Simulation +module Ao = AOBasis module El = Electrons -let hcore_guess simulation = - let eN_ints = Lazy.force simulation.Si.eN_ints - and kin_ints = Lazy.force simulation.Si.kin_ints +let hcore_guess ao_basis = + let eN_ints = Lazy.force ao_basis.Ao.eN_ints + and kin_ints = Lazy.force ao_basis.Ao.kin_ints in Mat.add eN_ints kin_ints -let huckel_guess simulation = +let huckel_guess ao_basis = let c = 0.5 *. 1.75 in - let ao_num = Basis.size simulation.Si.basis in - let eN_ints = Lazy.force simulation.Si.eN_ints - and kin_ints = Lazy.force simulation.Si.kin_ints - and overlap = Lazy.force simulation.Si.overlap - and m_X = Lazy.force simulation.Si.overlap_ortho + let ao_num = Basis.size ao_basis.Ao.basis in + let eN_ints = Lazy.force ao_basis.Ao.eN_ints + and kin_ints = Lazy.force ao_basis.Ao.kin_ints + and overlap = Lazy.force ao_basis.Ao.overlap + and m_X = Lazy.force ao_basis.Ao.ortho in let diag = Array.init (ao_num+1) (fun i -> if i=0 then 0. else eN_ints.{i,i} +. kin_ints.{i,i}) in - let nocc = - simulation.Si.electrons.El.n_alpha - in - let m_F = (Fock.make ~density:(gemm ~alpha:2. ~transb:`T ~k:nocc m_X m_X) simulation).Fock.fock in - for j=1 to ao_num do - for i=1 to ao_num do - if (i <> j) then - m_F.{i,j} <- c *. overlap.{i,j} *. (diag.(i) +. diag.(j)) (*TODO Pseudo *) + function + | 0 -> invalid_arg "Huckel guess needs a non-zero number of occupied MOs." + | nocc -> + let m_F = (Fock.make ~density:(gemm ~alpha:2. ~transb:`T ~k:nocc m_X m_X) ao_basis).Fock.fock in + for j=1 to ao_num do + for i=1 to ao_num do + if (i <> j) then + m_F.{i,j} <- c *. overlap.{i,j} *. (diag.(i) +. diag.(j)) (*TODO Pseudo *) + done; done; - done; - m_F + m_F -let make ~guess simulation = + +let make ?(nocc=0) ~guess ao_basis = match guess with - | `Hcore -> Hcore (hcore_guess simulation) - | `Huckel -> Huckel (huckel_guess simulation) + | `Hcore -> Hcore (hcore_guess ao_basis) + | `Huckel -> Huckel (huckel_guess ao_basis nocc) diff --git a/SCF/Guess.mli b/SCF/Guess.mli index 2fa54c2..be7fdd1 100644 --- a/SCF/Guess.mli +++ b/SCF/Guess.mli @@ -7,7 +7,7 @@ type guess = type t = guess -val make : guess:[ `Hcore | `Huckel ] -> Simulation.t -> t +val make : ?nocc:int -> guess:[ `Hcore | `Huckel ] -> AOBasis.t -> t diff --git a/SCF/RHF.ml b/SCF/RHF.ml index bb492de..5200afb 100644 --- a/SCF/RHF.ml +++ b/SCF/RHF.ml @@ -1,37 +1,44 @@ open Util open Lacaml.D -open Simulation + +module Si = Simulation +module El = Electrons +module Ao = AOBasis let make ?guess:(guess=`Huckel) ?max_scf:(max_scf=64) ?level_shift:(level_shift=0.1) ?threshold_SCF:(threshold_SCF=1.e-5) simulation = (* Number of occupied MOs *) let nocc = - simulation.electrons.Electrons.n_alpha + simulation.Si.electrons.El.n_alpha in let nuclear_repulsion = - simulation.nuclear_repulsion + simulation.Si.nuclear_repulsion + in + + let ao_basis = + simulation.Si.ao_basis in (* Initial guess *) let guess = - Guess.make ~guess simulation + Guess.make ~nocc ~guess ao_basis in (* Orthogonalization matrix *) let m_X = - Lazy.force simulation.overlap_ortho + Lazy.force ao_basis.Ao.ortho in (* Overlap matrix *) let m_S = - Lazy.force simulation.overlap + Lazy.force ao_basis.Ao.overlap in - let m_T = Lazy.force simulation.kin_ints - and m_V = Lazy.force simulation.eN_ints + let m_T = Lazy.force ao_basis.Ao.kin_ints + and m_V = Lazy.force ao_basis.Ao.eN_ints in (* Level shift in MO basis *) @@ -53,7 +60,7 @@ let make ?guess:(guess=`Huckel) ?max_scf:(max_scf=64) ?level_shift:(level_shift= (* Fock matrix in AO basis *) let m_F, m_Hc, m_J, m_K = let x = - Fock.make ~density:m_P simulation + Fock.make ~density:m_P ao_basis in x.Fock.fock, x.Fock.core, x.Fock.coulomb, x.Fock.exchange in diff --git a/Simulation.ml b/Simulation.ml index f404774..3b7bc43 100644 --- a/Simulation.ml +++ b/Simulation.ml @@ -1,15 +1,10 @@ type t = { - charge : Charge.t; - electrons : Electrons.t ; - basis : Basis.t; - nuclei : Nuclei.t; - overlap : Overlap.t lazy_t; - overlap_ortho : Orthonormalization.t lazy_t; - eN_ints : NucInt.t lazy_t; - kin_ints : KinInt.t lazy_t; - ee_ints : ERI.t lazy_t; + charge : Charge.t; + electrons : Electrons.t ; + nuclei : Nuclei.t; + basis : Basis.t; + ao_basis : AOBasis.t; nuclear_repulsion : float; - cartesian : bool; } let make ?cartesian:(cartesian=false) @@ -26,21 +21,22 @@ let make ?cartesian:(cartesian=false) let electrons = Electrons.make ~multiplicity ~charge nuclei in + let charge = Charge.(Nuclei.charge nuclei + Electrons.charge electrons) in - let overlap = - lazy (Overlap.of_basis basis) + + let ao_basis = + AOBasis.make ~basis ~cartesian nuclei in + + let nuclear_repulsion = + Nuclei.repulsion nuclei + in + { - charge ; - basis ; nuclei ; electrons ; overlap ; - overlap_ortho = lazy (Orthonormalization.make ~cartesian ~basis (Lazy.force overlap)); - eN_ints = lazy (NucInt.of_basis_nuclei basis nuclei); - kin_ints = lazy (KinInt.of_basis basis); - ee_ints = lazy (ERI.of_basis basis); - nuclear_repulsion = Nuclei.repulsion nuclei; - cartesian ; + charge ; basis ; nuclei ; electrons ; ao_basis ; + nuclear_repulsion ; } diff --git a/run_integrals.ml b/run_integrals.ml index 85472f8..4d16b54 100644 --- a/run_integrals.ml +++ b/run_integrals.ml @@ -36,10 +36,11 @@ let run ~out = print_float s.Simulation.nuclear_repulsion; print_newline (); print_endline @@ Basis.to_string s.Simulation.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 + let ao_basis = s.Simulation.ao_basis in + let overlap = Lazy.force ao_basis.AOBasis.overlap in + let eN_ints = Lazy.force ao_basis.AOBasis.eN_ints in + let kin_ints = Lazy.force ao_basis.AOBasis.kin_ints in + let ee_ints = Lazy.force ao_basis.AOBasis.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;