diff --git a/ao/lib/ao_dim.ml b/ao/lib/ao_dim.ml new file mode 100644 index 0000000..63c57c4 --- /dev/null +++ b/ao/lib/ao_dim.ml @@ -0,0 +1 @@ +type t diff --git a/ao/lib/ao_dim.mli b/ao/lib/ao_dim.mli new file mode 100644 index 0000000..63c57c4 --- /dev/null +++ b/ao/lib/ao_dim.mli @@ -0,0 +1 @@ +type t diff --git a/ao/lib/basis.ml b/ao/lib/basis.ml index 81b46f4..d541f50 100644 --- a/ao/lib/basis.ml +++ b/ao/lib/basis.ml @@ -9,6 +9,8 @@ type t = cartesian : bool } +type ao = Ao_dim.t + let of_nuclei_and_basis_filename ?(kind=`Gaussian) ?operators ?(cartesian=false) ~nuclei filename = match kind with diff --git a/ao/lib/basis.mli b/ao/lib/basis.mli index 4fdd054..38021b0 100644 --- a/ao/lib/basis.mli +++ b/ao/lib/basis.mli @@ -10,6 +10,7 @@ type basis = | Gaussian of Basis_gaussian.t type t +type ao = Ao_dim.t (** {1 Accessors} *) @@ -19,37 +20,37 @@ val size : t -> int val ao_basis : t -> basis (** One-electron basis set *) -val overlap : t -> ('a,'a) Matrix.t +val overlap : t -> (ao,ao) Matrix.t (** Overlap matrix *) -val multipole : t -> ('a,'a) Matrix.t array +val multipole : t -> (ao,ao) Matrix.t array (** Multipole matrices *) -val ortho : t -> ('a,'a) Matrix.t +val ortho : t -> (ao,'a) Matrix.t (** Orthonormalization matrix of the overlap *) -val eN_ints : t -> ('a,'a) Matrix.t +val eN_ints : t -> (ao,ao) Matrix.t (** Electron-nucleus potential integrals *) -val kin_ints : t -> ('a,'a) Matrix.t +val kin_ints : t -> (ao,ao) Matrix.t (** Kinetic energy integrals *) -val ee_ints : t -> 'a Four_idx_storage.t +val ee_ints : t -> ao Four_idx_storage.t (** Electron-electron potential integrals *) -val ee_lr_ints : t -> 'a Four_idx_storage.t +val ee_lr_ints : t -> ao Four_idx_storage.t (** Electron-electron long-range potential integrals *) -val f12_ints : t -> 'a Four_idx_storage.t +val f12_ints : t -> ao Four_idx_storage.t (** Electron-electron potential integrals *) -val f12_over_r12_ints : t -> 'a Four_idx_storage.t +val f12_over_r12_ints : t -> ao Four_idx_storage.t (** Electron-electron potential integrals *) val cartesian : t -> bool (** If true, use cartesian Gaussians (6d, 10f, ...) *) -val values : t -> Coordinate.t -> 'a Vector.t +val values : t -> Coordinate.t -> ao Vector.t (** Values of the AOs evaluated at a given point *) diff --git a/ao/test/ao_basis.ml b/ao/test/ao_basis.ml index 2a8d359..8215d68 100644 --- a/ao/test/ao_basis.ml +++ b/ao/test/ao_basis.ml @@ -33,21 +33,21 @@ let make_tests name t = let test_overlap () = let reference = - Matrix.sym_matrix_of_file (wd ^ Filename.dir_sep ^ name ^ "_overlap.ref") + Matrix.sym_of_file (wd ^ Filename.dir_sep ^ name ^ "_overlap.ref") in check_matrix "Overlap" (overlap t) reference in let test_eN_ints () = let reference = - Matrix.sym_matrix_of_file (wd ^ Filename.dir_sep ^ name ^ "_nuc.ref") + Matrix.sym_of_file (wd ^ Filename.dir_sep ^ name ^ "_nuc.ref") in check_matrix "eN_ints" (eN_ints t) reference in let test_kin_ints () = let reference = - Matrix.sym_matrix_of_file (wd ^ Filename.dir_sep ^ name ^ "_kin.ref") + Matrix.sym_of_file (wd ^ Filename.dir_sep ^ name ^ "_kin.ref") in check_matrix "kin_ints" (kin_ints t) reference in diff --git a/ao/test/ao_basis_gaussian.ml b/ao/test/ao_basis_gaussian.ml index b68d948..5130cda 100644 --- a/ao/test/ao_basis_gaussian.ml +++ b/ao/test/ao_basis_gaussian.ml @@ -35,21 +35,21 @@ let make_tests name t = let test_overlap () = let reference = - Matrix.sym_matrix_of_file (wd ^ Filename.dir_sep ^ name ^ "_overlap.ref") + Matrix.sym_of_file (wd ^ Filename.dir_sep ^ name ^ "_overlap.ref") in check_matrix "Overlap" (overlap t) reference in let test_eN_ints () = let reference = - Matrix.sym_matrix_of_file (wd ^ Filename.dir_sep ^ name ^ "_nuc.ref") + Matrix.sym_of_file (wd ^ Filename.dir_sep ^ name ^ "_nuc.ref") in check_matrix "eN_ints" (eN_ints t) reference in let test_kin_ints () = let reference = - Matrix.sym_matrix_of_file (wd ^ Filename.dir_sep ^ name ^ "_kin.ref") + Matrix.sym_of_file (wd ^ Filename.dir_sep ^ name ^ "_kin.ref") in check_matrix "kin_ints" (kin_ints t) reference in @@ -69,7 +69,6 @@ let make_tests name t = in check_eri (ee_lr_ints t) reference in - [ "Overlap", `Quick, test_overlap; "eN_ints", `Quick, test_eN_ints; diff --git a/linear_algebra/lib/four_idx_storage.ml b/linear_algebra/lib/four_idx_storage.ml index 5c97ea8..d73b74b 100644 --- a/linear_algebra/lib/four_idx_storage.ml +++ b/linear_algebra/lib/four_idx_storage.ml @@ -426,13 +426,13 @@ let four_index_transform_dense_sparse ds coef source = Matrix.gemm_inplace ~transa:`T ~c:p o coef; (* p_j_kalpha *) - let p' = Matrix.reshape p ao_num ao_mo_num in + let p' = Matrix.reshape_inplace ao_num ao_mo_num p in (* q_kalpha_beta = \sum_j p_j_kalpha c_j_beta *) Matrix.gemm_inplace ~transa:`T ~c:q p' coef; (* q_k_alphabeta = \sum_j p_j_kalpha c_j_beta *) - let q' = Matrix.reshape q ao_num mo_num_2 in + let q' = Matrix.reshape_inplace ao_num mo_num_2 q in (* u_alphabeta_gamma = \sum_k q_k_alphabeta c_k_gamma *) Matrix.gemm_inplace ~transa:`T ~beta:1. ~alpha:coefx.{l,delta} ~c:u q' coef ; diff --git a/linear_algebra/lib/linear_algebra.ml b/linear_algebra/lib/linear_algebra.ml new file mode 100644 index 0000000..c0a3d3b --- /dev/null +++ b/linear_algebra/lib/linear_algebra.ml @@ -0,0 +1,10 @@ +let (%.) = Vector.(%.) +let (%:) = Matrix.(%:) + +module Conventions = Conventions +module Diis = Diis +module Four_idx_storage = Four_idx_storage +module Matrix = Matrix +module Orthonormalization = Orthonormalization +module Spherical_to_cartesian = Spherical_to_cartesian +module Vector = Vector diff --git a/linear_algebra/lib/matrix.ml b/linear_algebra/lib/matrix.ml index 0e8af4a..8a4927e 100644 --- a/linear_algebra/lib/matrix.ml +++ b/linear_algebra/lib/matrix.ml @@ -48,7 +48,7 @@ external of_bigarray_inplace : (float, Stdlib.Bigarray.float64_elt, Stdlib.Bigar let to_bigarray t = lacpy t let of_bigarray t = lacpy t -let reshape a m n = +let reshape_inplace m n a = assert ( (dim1 a) * (dim2 a) = m*n); let b = to_bigarray a @@ -56,6 +56,14 @@ let reshape a m n = in Bigarray.reshape_2 b m n +let reshape_vec_inplace m n v = + assert ( Vector.dim v = m*n); + let b = + Vector.to_bigarray_inplace v + |> Bigarray.genarray_of_array1 + in + Bigarray.reshape_2 b m n + let col_inplace t j = Mat.col t j |> Vector.of_bigarray_inplace @@ -200,6 +208,19 @@ let of_col_vecs_list a = |> List.rev |> Mat.of_col_vecs_list +let of_array a = + Bigarray.Array2.of_array Bigarray.Float64 Bigarray.fortran_layout a + +let to_array a = + let result = Array.make_matrix (Mat.dim1 a) (Mat.dim2 a) 0. in + for i=1 to Mat.dim1 a do + let result_i = result.(i-1) in + for j=1 to Mat.dim2 a do + result_i.(j-1) <- a.{i,j} + done + done; + result + let normalize_mat_inplace t = let norm = Mat.as_vec t |> nrm2 in Mat.scal norm t @@ -227,6 +248,9 @@ let x_o_xt ~o ~x = gemm o x ~transb:`T |> gemm x +let amax t = + Mat.as_vec t |> amax + let pp ppf m = let rows = Mat.dim1 m and cols = Mat.dim2 m @@ -263,7 +287,7 @@ let outer_product ?(alpha=1.0) u v = outer_product_inplace m ~alpha u v; m -let matrix_of_file filename = +let of_file filename = let ic = Scanf.Scanning.open_in filename in let rec read_line accu = let result = @@ -290,9 +314,9 @@ let matrix_of_file filename = result -let sym_matrix_of_file filename = +let sym_of_file filename = let result = - matrix_of_file filename + of_file filename in for j=1 to Mat.dim1 result do for i=1 to j do diff --git a/linear_algebra/lib/matrix.mli b/linear_algebra/lib/matrix.mli index fb74db1..9b3d92e 100644 --- a/linear_algebra/lib/matrix.mli +++ b/linear_algebra/lib/matrix.mli @@ -17,9 +17,12 @@ val make0 : int -> int -> ('a,'b) t val create : int -> int -> ('a,'b) t (** Creates an uninitialized matrix. *) -val reshape : ('a,'b) t -> int -> int -> ('c,'d) t +val reshape_inplace : int -> int -> ('a,'b) t -> ('c,'d) t (** Changes the dimensions of the matrix *) +val reshape_vec_inplace : int -> int -> ('a*'b) Vector.t -> ('a,'b) t +(** Reshapres a vector into a matrix *) + val init_cols : int -> int -> (int -> int -> float) -> ('a,'b) t (** Creates an uninitialized matrix. *) @@ -59,6 +62,9 @@ val mul : ('a,'b) t -> ('a,'b) t -> ('a,'b) t val div : ('a,'b) t -> ('a,'b) t -> ('a,'b) t (** Divides two matrices element-wise *) +val amax : ('a,'b) t -> float +(** Maximum of the absolute values of the elements of the matrix. *) + val add_inplace : c:('a,'b) t -> ('a,'b) t -> ('a,'b) t -> unit (** [add_inplace c a b] : performs [c = a+b] in-place. *) @@ -99,6 +105,12 @@ val of_col_vecs : 'a Vector.t array -> ('a,'b) t val of_col_vecs_list : 'a Vector.t list -> ('a,'b) t (** Converts a list of vectors into a matrix *) +val to_array : ('a,'b) t -> float array array +(** Converts the matrix into an array of arrays *) + +val of_array : float array array -> ('a,'b) t +(** Converts an array of arrays into a matrix *) + val copy: ?m:int -> ?n:int -> ?br:int -> ?bc:int -> ?ar:int -> ?ac:int -> ('a,'b) t -> ('a,'b) t (** Copies all or part of a two-dimensional matrix A to a new matrix B *) @@ -228,7 +240,7 @@ val gemm_nn: ?m:int -> ?n:int -> ?k:int -> ?beta:float -> (** Performs gemm with [~transa=`N] and [~transb=`N]. *) val gemm_nt: ?m:int -> ?n:int -> ?k:int -> ?beta:float -> - ?c:('a,'c) t -> ?alpha:float -> ('a,'b) t -> ('c,'b) t -> ('c,'a) t + ?c:('a,'c) t -> ?alpha:float -> ('a,'b) t -> ('c,'b) t -> ('a,'c) t (** Performs gemm with [~transa=`N] and [~transb=`T]. *) val gemm_tn: ?m:int -> ?n:int -> ?k:int -> ?beta:float -> @@ -236,7 +248,7 @@ val gemm_tn: ?m:int -> ?n:int -> ?k:int -> ?beta:float -> (** Performs gemm with [~transa=`T] and [~transb=`N]. *) val gemm_tt: ?m:int -> ?n:int -> ?k:int -> ?beta:float -> - ?c:('a,'c) t -> ?alpha:float -> ('b,'a) t -> ('c,'b) t -> ('c,'b) t + ?c:('a,'c) t -> ?alpha:float -> ('b,'a) t -> ('c,'b) t -> ('a,'c) t (** Performs gemm with [~transa=`T] and [~transb=`T]. *) @@ -285,13 +297,13 @@ val x_o_xt : o:('b,'b) t -> x:('a,'b) t -> ('a,'a) t val debug_matrix: string -> ('a,'b) t -> unit (** Prints a matrix in stdout for debug *) -val matrix_of_file : string -> ('a,'b) t +val of_file : string -> ('a,'b) t (** Reads a matrix from a file with format "%d %d %f" corresponding to [i, j, A.{i,j}]. *) val relabel : ('a,'b) t -> ('c,'d) t -val sym_matrix_of_file : string -> ('a,'b) t +val sym_of_file : string -> ('a,'b) t (** Reads a symmetric matrix from a file with format "%d %d %f" corresponding to [i, j, A.{i,j}]. *) diff --git a/linear_algebra/lib/vector.ml b/linear_algebra/lib/vector.ml index 47b68a5..938ee2b 100644 --- a/linear_algebra/lib/vector.ml +++ b/linear_algebra/lib/vector.ml @@ -32,6 +32,7 @@ let sub t1 t2 = Vec.sub t1 t2 let mul t1 t2 = Vec.mul t1 t2 let div t1 t2 = Vec.div t1 t2 let dot t1 t2 = dot t1 t2 +let amax t = amax t let create n = Vec.create n let make0 n = Vec.make0 n diff --git a/linear_algebra/lib/vector.mli b/linear_algebra/lib/vector.mli index 3b6ec62..1db391e 100644 --- a/linear_algebra/lib/vector.mli +++ b/linear_algebra/lib/vector.mli @@ -67,6 +67,9 @@ val asin : 'a t -> 'a t val atan : 'a t -> 'a t (** [atan t = map (f x -> atan x) t *) +val amax : 'a t -> float +(** Maximum of the absolute values of the elements of the vector. *) + val normalize : 'a t -> 'a t (** Returns the vector normalized *) diff --git a/linear_algebra/test/vector.ml b/linear_algebra/test/vector.ml index 99fee9e..3a17163 100644 --- a/linear_algebra/test/vector.ml +++ b/linear_algebra/test/vector.ml @@ -41,8 +41,8 @@ let test_all () = check (float 1.e-14) "norm" (sqrt (dot u2 u2)) (Vector.norm v2); check (float 1.e-14) "sum" (Vec.sum u1) (Vector.sum v1); check (float 1.e-14) "sum" (Vec.sum u2) (Vector.sum v2); - check (float 1.e-14) "at" (u1.{n/2}) (Vector.at v1 (n/2)); - check (float 1.e-14) "at" (u2.{n/2}) (Vector.at v2 (n/2)); + check (float 1.e-14) "at" (u1.{n/2}) (v1%.(n/2)); + check (float 1.e-14) "at" (u2.{n/2}) (v2%.(n/2)); check (bool) "of_list" true (v1 = Vector.of_list @@ Array.to_list a1); check (bool) "of_list" true (v2 = Vector.of_list @@ Array.to_list a2); check (bool) "to_list" true (Vector.to_list v1 = Array.to_list a1); diff --git a/mo/lib/basis.ml b/mo/lib/basis.ml index d9b21aa..7e8c619 100644 --- a/mo/lib/basis.ml +++ b/mo/lib/basis.ml @@ -1,12 +1,13 @@ open Linear_algebra -open Common.Util -open Common.Constants (** One-electron orthogonal basis set, corresponding to Molecular Orbitals. *) -module HF = HartreeFock +module HF = Hartree_fock module Si = Simulation +type ao = Ao.Ao_dim.t +type mo = Mo_dim.t + type mo_type = | RHF | ROHF | UHF | CASSCF | Projected | Natural of string @@ -16,12 +17,12 @@ type t = { simulation : Simulation.t; (* Simulation which produced the MOs *) mo_type : mo_type; (* Kind of MOs (RHF, CASSCF, Localized...) *) - mo_occupation : t Vector.t; (* Occupation numbers *) - mo_coef : (Ao.Basis.t,t) Matrix.t; (* Matrix of the MO coefficients in the AO basis *) - eN_ints : (t,t) Matrix.t lazy_t; (* Electron-nucleus potential integrals *) - ee_ints : t Four_idx_storage.t lazy_t; (* Electron-electron potential integrals *) - kin_ints : (t,t) Matrix.t lazy_t; (* Kinetic energy integrals *) - one_e_ints : (t,t) Matrix.t lazy_t; (* One-electron integrals *) + mo_occupation : mo Vector.t; (* Occupation numbers *) + mo_coef : (ao,mo) Matrix.t; (* Matrix of the MO coefficients in the AO basis *) + eN_ints : (mo,mo) Matrix.t lazy_t; (* Electron-nucleus potential integrals *) + ee_ints : mo Four_idx_storage.t lazy_t; (* Electron-electron potential integrals *) + kin_ints : (mo,mo) Matrix.t lazy_t; (* Kinetic energy integrals *) + one_e_ints : (mo,mo) Matrix.t lazy_t; (* One-electron integrals *) (* TODO f12_ints : F12.t lazy_t; (* F12 integrals *) *) @@ -50,16 +51,16 @@ let mo_energies t = let m_C = mo_coef t in let f = let m_N = Matrix.of_diag @@ mo_occupation t in - let m_P = Matrix.x_o_xt m_N m_C in + let m_P = Matrix.x_o_xt ~o:m_N ~x:m_C in match t.mo_type with | RHF -> Fock.make_rhf ~density:m_P (ao_basis t) | Projected - | ROHF -> (Matrix.scal 0.5 m_P; + | ROHF -> (Matrix.scale_inplace 0.5 m_P; Fock.make_uhf ~density_same:m_P ~density_other:m_P (ao_basis t)) | _ -> failwith "Not implemented" in let m_F0 = Fock.fock f in - Matrix.xt_o_x m_F0 m_C + Matrix.xt_o_x ~o:m_F0 ~x:m_C |> Matrix.diag @@ -86,7 +87,7 @@ let make ~simulation ~mo_type ~mo_occupation ~mo_coef () = ) and ee_ints = lazy ( Ao.Basis.ee_ints ao_basis - |> Eri.four_index_transform mo_coef + |> Four_idx_storage.four_index_transform mo_coef ) (* and f12_ints = lazy ( @@ -106,7 +107,7 @@ let make ~simulation ~mo_type ~mo_occupation ~mo_coef () = let values t point = let c = mo_coef t in let a = Ao.Basis.values (Simulation.ao_basis t.simulation) point in - Matrix.gemv ~trans:`T c a + Matrix.gemv_t c a let of_hartree_fock hf = let mo_coef = HF.eigenvectors hf in @@ -114,24 +115,24 @@ let of_hartree_fock hf = let mo_occupation = HF.occupation hf in let mo_type = match HF.kind hf with - | HartreeFock.RHF -> RHF - | HartreeFock.ROHF -> ROHF - | HartreeFock.UHF -> UHF + | HF.RHF -> RHF + | HF.ROHF -> ROHF + | HF.UHF -> UHF in make ~simulation ~mo_type ~mo_occupation ~mo_coef () - +(* let of_mo_basis simulation other = let mo_coef = let basis = Simulation.ao_basis simulation in let basis_other = ao_basis other in let m_S = - Overlap.(matrix @@ of_basis_pair - (AOBasis.basis basis) - (AOBasis.basis basis_other) ) + Ao.Overlap.(of_basis_pair + (Ao.Basis.ao_basis basis) + (Ao.Basis.ao_basis basis_other) ) in - let m_X = AOBasis.ortho basis in + let m_X = Ao.Basis.ortho basis in (* Project other vectors in the current basis *) let m_C = gemm m_S @@ mo_coef other @@ -156,15 +157,15 @@ let of_mo_basis simulation other = else 0.) in make ~simulation ~mo_type:Projected ~mo_occupation ~mo_coef () +*) let pp ?(start=1) ?(finish=0) ppf t = - let open Lacaml.Io in - let rows = Mat.dim1 t.mo_coef - and cols = Mat.dim2 t.mo_coef + let rows = Matrix.dim1 t.mo_coef + and cols = Matrix.dim2 t.mo_coef in let finish = match finish with @@ -182,7 +183,7 @@ let pp ?(start=1) ?(finish=0) ppf t = Array.iteri (fun i x -> if (i+1 >= first) && (i+1 <= first+4 ) then Format.fprintf ppf "%12f@ " x) - (Vec.to_array @@ mo_energies t); + (Vector.to_array @@ mo_energies t); Format.fprintf ppf "@]@;"; Format.fprintf ppf "@[%a@]" @@ -193,7 +194,8 @@ let pp ?(start=1) ?(finish=0) ppf t = (Array.init (min 5 (cols-first+1)) (fun i -> Printf.sprintf "-- %d --" (i + first) )) ~print_right:false ~print_foot:false - () ) (lacpy ~ac:first ~n:(min 5 (cols-first+1)) (t.mo_coef)) ; + () ) (Matrix.to_bigarray_inplace t.mo_coef + |> Lacaml.D.lacpy ~ac:first ~n:(min 5 (cols-first+1)) ) ; Format.fprintf ppf "@]@;@;@]"; (aux [@tailcall]) (first+5) end diff --git a/mo/lib/basis.mli b/mo/lib/basis.mli index b243765..01243ca 100644 --- a/mo/lib/basis.mli +++ b/mo/lib/basis.mli @@ -13,6 +13,8 @@ type mo_type = | Localized of string type t +type mo = Mo_dim.t +type ao = Ao.Ao_dim.t (** {1 Accessors} *) @@ -25,25 +27,25 @@ val mo_type : t -> mo_type val ao_basis : t -> Ao.Basis.t (** Matrix of the MO coefficients in the AO basis *) -val mo_occupation : t -> t Vector.t +val mo_occupation : t -> mo Vector.t (** Occupation numbers *) -val mo_coef : t -> (Ao.Basis.t, t) Matrix.t +val mo_coef : t -> (ao, mo) Matrix.t (** Molecular orbitcal coefficients *) -val eN_ints : t -> (t,t) Matrix.t +val eN_ints : t -> (mo,mo) Matrix.t (** Electron-nucleus potential integrals *) -val ee_ints : t -> t Four_idx_storage.t +val ee_ints : t -> mo Four_idx_storage.t (** Electron-electron repulsion integrals *) -val kin_ints : t -> (t,t) Matrix.t +val kin_ints : t -> (mo,mo) Matrix.t (** Kinetic energy integrals *) -val one_e_ints : t -> (t,t) Matrix.t +val one_e_ints : t -> (mo,mo) Matrix.t (** One-electron integrals {% $\hat{T} + V$ %} *) -val two_e_ints : t -> t Four_idx_storage.t +val two_e_ints : t -> mo Four_idx_storage.t (** Electron-electron repulsion integrals *) (* TODO @@ -54,39 +56,38 @@ val f12_ints : t -> F12.t val size : t -> int (** Number of molecular orbitals in the basis *) -val mo_energies : t -> t Vector.t +val mo_energies : t -> mo Vector.t (** Fock MO energies *) -val values : t -> Coordinate.t -> t Vector.t +val values : t -> Coordinate.t -> mo Vector.t (** Values of the MOs evaluated at a given coordinate. *) (** {1 Creators} *) val make : simulation:Simulation.t -> mo_type:mo_type -> - mo_occupation:t Vector.t -> - mo_coef:(Ao.Basis.t,t) Matrix.t -> + mo_occupation:mo Vector.t -> + mo_coef:(ao,mo) Matrix.t -> unit -> t (** Function to build a data structure representing the molecular orbitals. *) -val of_hartree_fock : HartreeFock.t -> t +val of_hartree_fock : Hartree_fock.t -> t (** Build MOs from a Restricted Hartree-Fock calculation. *) +(* val of_mo_basis : Simulation.t -> t -> t (** Project the MOs of the other basis on the current one. *) +*) val mo_matrix_of_ao_matrix : - mo_coef:(Ao.Basis.t,t) Matrix.t -> - (Ao.Basis.t,Ao.Basis.t) Matrix.t -> - (t,t) Matrix.t + mo_coef:(ao,mo) Matrix.t -> + (ao,ao) Matrix.t -> (mo,mo) Matrix.t (** Build a matrix in MO basis from a matrix in AO basis. *) val ao_matrix_of_mo_matrix : - mo_coef:(Ao.Basis.t,t) Matrix.t -> - ao_overlap:(Ao.Basis.t,Ao.Basis.t) Matrix.t -> - (t,t) Matrix.t -> - (Ao.Basis.t,Ao.Basis.t) Matrix.t + mo_coef:(ao,mo) Matrix.t -> ao_overlap:(ao,ao) Matrix.t -> + (mo,mo) Matrix.t -> (ao,ao) Matrix.t (** Build a matrix in AO basis from a matrix in MO basis. *) (** {1 Printers} *) diff --git a/mo/lib/class.ml b/mo/lib/class.ml index 265fe4b..51591ac 100644 --- a/mo/lib/class.ml +++ b/mo/lib/class.ml @@ -1,4 +1,5 @@ open Particles +open Common type mo_class = | Core of int (* Always doubly occupied *) diff --git a/mo/lib/fock.ml b/mo/lib/fock.ml new file mode 100644 index 0000000..27192e3 --- /dev/null +++ b/mo/lib/fock.ml @@ -0,0 +1,240 @@ +open Linear_algebra +open Common + +type ao = Ao.Ao_dim.t + +type t = + { + fock : (ao,ao) Matrix.t ; + core : (ao,ao) Matrix.t ; + coulomb : (ao,ao) Matrix.t ; + exchange : (ao,ao) Matrix.t ; + } + + +let fock t = t.fock +let core t = t.core +let coulomb t = t.coulomb +let exchange t = t.exchange + + +let make_rhf ~density ?(threshold=Constants.epsilon) ao_basis = + let m_P = density + and m_T = Ao.Basis.kin_ints ao_basis + and m_V = Ao.Basis.eN_ints ao_basis + and m_G = Ao.Basis.ee_ints ao_basis + in + let nBas = Matrix.dim1 m_T + in + + let m_Hc = Matrix.add m_T m_V + and m_J = Array.make_matrix nBas nBas 0. + and m_K = Array.make_matrix nBas nBas 0. + in + + for sigma = 1 to nBas do + let m_Ksigma = m_K.(sigma-1) in + for nu = 1 to nBas do + let m_Jnu = m_J.(nu-1) in + for lambda = 1 to nBas do + let pJ = m_P%:(lambda,sigma) + and pK = 0.5 *. (m_P%:(lambda,nu)) + in + match (abs_float pJ > threshold , abs_float pK > threshold, nu < sigma) with + | (false, false, _) -> () + | (true , true , true) -> + begin + for mu = 1 to nu do + let integral = + Four_idx_storage.get_phys m_G mu lambda nu sigma + in + if (integral <> 0.) then begin + m_Jnu.(mu-1) <- m_Jnu.(mu-1) +. pJ *. integral; + m_Ksigma.(mu-1) <- m_Ksigma.(mu-1) +. pK *. integral + end + done; + for mu = nu+1 to sigma do + m_Ksigma.(mu-1) <- m_Ksigma.(mu-1) +. + pK *. Four_idx_storage.get_phys m_G mu lambda nu sigma + done + end + | (true , true , false) -> + begin + for mu = 1 to sigma do + let integral = + Four_idx_storage.get_phys m_G mu lambda nu sigma + in + if (integral <> 0.) then begin + m_Jnu.(mu-1) <- m_Jnu.(mu-1) +. pJ *. integral; + m_Ksigma.(mu-1) <- m_Ksigma.(mu-1) +. pK *. integral + end + done; + for mu = sigma+1 to nu do + m_Jnu.(mu-1) <- + m_Jnu.(mu-1) +. pJ *. Four_idx_storage.get_phys m_G mu lambda nu sigma + done + end + | (false, true , _) -> + for mu = 1 to sigma do + m_Ksigma.(mu-1) <- + m_Ksigma.(mu-1) +. pK *. Four_idx_storage.get_phys m_G mu lambda nu sigma + done + | (true , false, _) -> + for mu = 1 to nu do + m_Jnu.(mu-1) <- + m_Jnu.(mu-1) +. pJ *. Four_idx_storage.get_phys m_G mu lambda nu sigma + done + done + done; + for mu = 1 to sigma-1 do + m_K.(mu-1).(sigma-1) <- m_Ksigma.(mu-1); + done + done; + for nu = 1 to nBas do + let m_Jnu = m_J.(nu-1) in + for mu = 1 to nu-1 do + m_J.(mu-1).(nu-1) <- m_Jnu.(mu-1) + done + done; + + let m_J = Matrix.of_array m_J + and m_K = Matrix.of_array m_K + in + { fock = Matrix.add m_Hc (Matrix.sub m_J m_K) ; + core = m_Hc ; coulomb = m_J ; exchange = m_K } + + + +let make_uhf ~density_same ~density_other ?(threshold=Constants.epsilon) ao_basis = + let m_P_a = density_same + and m_P_b = density_other + and m_T = Ao.Basis.kin_ints ao_basis + and m_V = Ao.Basis.eN_ints ao_basis + and m_G = Ao.Basis.ee_ints ao_basis + in + let nBas = Matrix.dim1 m_T + in + + let m_Hc = Matrix.add m_T m_V + and m_J = Array.make_matrix nBas nBas 0. + and m_K = Array.make_matrix nBas nBas 0. + in + + for sigma = 1 to nBas do + let m_Ksigma = m_K.(sigma-1) in + for nu = 1 to nBas do + let m_Jnu = m_J.(nu-1) in + for lambda = 1 to nBas do + let pJ = m_P_a%:(lambda,sigma) +. m_P_b%:(lambda,sigma) + and pK = m_P_a%:(lambda,nu) + in + match (abs_float pJ > threshold , abs_float pK > threshold, nu < sigma) with + | (false, false, _) -> () + | (true , true , true) -> + begin + for mu = 1 to nu do + let integral = + Four_idx_storage.get_phys m_G mu lambda nu sigma + in + if (integral <> 0.) then begin + m_Jnu.(mu-1) <- m_Jnu.(mu-1) +. pJ *. integral; + m_Ksigma.(mu-1) <- m_Ksigma.(mu-1) +. pK *. integral + end + done; + for mu = nu+1 to sigma do + m_Ksigma.(mu-1) <- m_Ksigma.(mu-1) +. + pK *.Four_idx_storage.get_phys m_G mu lambda nu sigma + done + end + | (true , true , false) -> + begin + for mu = 1 to sigma do + let integral = + Four_idx_storage.get_phys m_G mu lambda nu sigma + in + if (integral <> 0.) then begin + m_Jnu.(mu-1) <- m_Jnu.(mu-1) +. pJ *. integral; + m_Ksigma.(mu-1) <- m_Ksigma.(mu-1) +. pK *. integral + end + done; + for mu = sigma+1 to nu do + m_Jnu.(mu-1) <- + m_Jnu.(mu-1) +. pJ *. Four_idx_storage.get_phys m_G mu lambda nu sigma + done + end + | (false, true , _) -> + for mu = 1 to sigma do + m_Ksigma.(mu-1) <- + m_Ksigma.(mu-1) +. pK *. Four_idx_storage.get_phys m_G mu lambda nu sigma + done + | (true , false, _) -> + for mu = 1 to nu do + m_Jnu.(mu-1) <- + m_Jnu.(mu-1) +. pJ *. Four_idx_storage.get_phys m_G mu lambda nu sigma + done + done + done; + for mu = 1 to sigma-1 do + m_K.(mu-1).(sigma-1) <- m_Ksigma.(mu-1); + done + done; + for nu = 1 to nBas do + let m_Jnu = m_J.(nu-1) in + for mu = 1 to nu-1 do + m_J.(mu-1).(nu-1) <- m_Jnu.(mu-1) + done + done; + + let m_J = Matrix.of_array m_J + and m_K = Matrix.of_array m_K + in + { fock = Matrix.add m_Hc (Matrix.sub m_J m_K) ; + core = m_Hc ; coulomb = m_J ; exchange = m_K } + + + + + + + +let op ~f f1 f2 = + assert (f1.core = f2.core); + let m_Hc = f1.core + and m_J = f f1.coulomb f2.coulomb + and m_K = f f1.exchange f2.exchange + in + { + fock = Matrix.add m_Hc (Matrix.sub m_J m_K); + core = m_Hc; + coulomb = m_J; + exchange = m_K; + } + + +let add = op ~f:(fun a b -> Matrix.add a b) + +let sub = op ~f:(fun a b -> Matrix.sub a b) + +let scale alpha f1 = + let m_Hc = f1.core + and m_J = Matrix.copy f1.coulomb + and m_K = Matrix.copy f1.exchange + in + Matrix.scale_inplace alpha m_J; + Matrix.scale_inplace alpha m_K; + { + fock = Matrix.add m_Hc (Matrix.sub m_J m_K); + core = m_Hc; + coulomb = m_J; + exchange = m_K; + } + + + +let pp ppf a = + Format.fprintf ppf "@[<2>"; + Format.fprintf ppf "@[ Fock matrix:@[<2>@[%a@]@.]@]" Matrix.pp a.fock; + Format.fprintf ppf "@[ Core Hamiltonian:@[<2>@[%a@]@.]@]" Matrix.pp a.core; + Format.fprintf ppf "@[ Coulomb matrix:@[<2>@[%a@]@.]@]" Matrix.pp a.coulomb; + Format.fprintf ppf "@[ Exchange matrix:@[<2>@[%a@]@.]@]" Matrix.pp a.exchange; + Format.fprintf ppf "@]" diff --git a/mo/lib/fock.mli b/mo/lib/fock.mli new file mode 100644 index 0000000..bd6d408 --- /dev/null +++ b/mo/lib/fock.mli @@ -0,0 +1,54 @@ +(** Type for the Fock operator in AO basis. *) + +open Linear_algebra + +type t +type ao = Ao.Ao_dim.t + +(** {1 Accessors} *) + +val fock : t -> (ao,ao) Matrix.t +(** Fock matrix in AO basis *) + +val core : t -> (ao,ao) Matrix.t +(** Core Hamiltonian : {% $\langle i | \hat{h} | j \rangle$ %} *) + +val coulomb : t -> (ao,ao) Matrix.t +(** Coulomb matrix : {% $\langle i | J | j \rangle$ %} *) + +val exchange : t -> (ao,ao) Matrix.t +(** Exchange matrix : {% $\langle i | K | j \rangle$ %} *) + + +(** {1 Creators} *) + +val make_rhf : density:(ao,ao) Matrix.t -> ?threshold:float -> Ao.Basis.t -> t +(** Create a Fock operator in the RHF formalism. Expected density is + {% $2 \mathbf{C\, C}^\dagger$ %}. [threshold] is a threshold on the + integrals. *) + +val make_uhf : + density_same:(ao,ao) Matrix.t -> + density_other:(ao,ao) Matrix.t -> + ?threshold:float -> + Ao.Basis.t -> t +(** Create a Fock operator in the UHF formalism. Expected density is + {% $\mathbf{C\, C}^\dagger$ %}. When building the {% $\alpha$ %} Fock + operator, [density_same] is the {% $\alpha$ %} density and [density_other] + is the {% $\beta$ %} density. [threshold] is a threshold on the integrals. *) + + +(** {1 Operations} *) + +val add : t -> t -> t +(** Add two Fock operators sharing the same core Hamiltonian. *) + +val sub : t -> t -> t +(** Subtract two Fock operators sharing the same core Hamiltonian. *) + +val scale : float -> t -> t +(** Scale the Fock matrix by a constant *) + +(** {1 Printers} *) + +val pp : Format.formatter -> t -> unit diff --git a/mo/lib/guess.ml b/mo/lib/guess.ml new file mode 100644 index 0000000..15733d3 --- /dev/null +++ b/mo/lib/guess.ml @@ -0,0 +1,56 @@ +open Linear_algebra + +type ao = Ao.Ao_dim.t +type mo = Mo_dim.t + +type guess = +| Hcore of (ao,ao) Matrix.t +| Huckel of (ao,ao) Matrix.t +| Matrix of (ao,mo) Matrix.t + +type t = guess + +module El = Particles.Electrons + +let hcore_guess ao_basis = + let eN_ints = Ao.Basis.eN_ints ao_basis + and kin_ints = Ao.Basis.kin_ints ao_basis + in + Matrix.add eN_ints kin_ints + + +let huckel_guess ao_basis = + let c = 0.5 *. 1.75 in + let eN_ints = Ao.Basis.eN_ints ao_basis + and kin_ints = Ao.Basis.kin_ints ao_basis + in + let m_F = + Matrix.add eN_ints kin_ints + in + let ao_num = Ao.Basis.size ao_basis + and overlap = Ao.Basis.overlap ao_basis + in + let diag = Vector.init ao_num (fun i -> m_F%:(i,i) ) in + + function + | 0 -> invalid_arg "Huckel guess needs a non-zero number of occupied MOs." + | _nocc -> + Matrix.init_cols ao_num ao_num (fun i j -> + if (i<>j) then + if (diag%.(i) +. diag%.(j)) < 0. then + c *. (overlap%:(i,j)) *. (diag%.(i) +. diag%.(j)) +. m_F%:(i,j) (*TODO Pseudo *) + else + m_F%:(i,j) (*TODO Pseudo *) + else + diag%.(i) + ) + + +let make ?(nocc=0) ~guess ao_basis = + match guess with + | `Hcore -> Hcore (hcore_guess ao_basis) + | `Huckel -> Huckel (huckel_guess ao_basis nocc) + | `Matrix m -> Matrix m + + + diff --git a/mo/lib/guess.mli b/mo/lib/guess.mli new file mode 100644 index 0000000..967a611 --- /dev/null +++ b/mo/lib/guess.mli @@ -0,0 +1,20 @@ +open Linear_algebra + +(** Guess for Hartree-Fock calculations. *) + +type ao = Ao.Ao_dim.t +type mo = Mo_dim.t + +type guess = +| Hcore of (ao,ao) Matrix.t (* Core Hamiltonian Matrix *) +| Huckel of (ao,ao) Matrix.t (* Huckel Hamiltonian Matrix *) +| Matrix of (ao,mo) Matrix.t (* Guess Eigenvectors *) + +type t = guess + + +val make : + ?nocc:int -> + guess:[ `Hcore | `Huckel | `Matrix of (ao,mo) Matrix.t ] -> + Ao.Basis.t -> t + diff --git a/mo/lib/hartree_fock.ml b/mo/lib/hartree_fock.ml new file mode 100644 index 0000000..8e307a2 --- /dev/null +++ b/mo/lib/hartree_fock.ml @@ -0,0 +1,721 @@ +open Linear_algebra +open Particles +open Common + +type ao = Ao.Ao_dim.t +type mo = Mo_dim.t + +type hartree_fock_data = + { + iteration : int ; + coefficients : (ao, mo) Matrix.t option ; + eigenvalues : mo Vector.t option ; + error : float option ; + diis : (mo*mo) Diis.t option ; + energy : float option ; + density : (ao,ao) Matrix.t option ; + density_a : (ao,ao) Matrix.t option ; + density_b : (ao,ao) Matrix.t option ; + fock : Fock.t option ; + fock_a : Fock.t option ; + fock_b : Fock.t option ; + } + +and hartree_fock_kind = +| RHF (** Restricted Hartree-Fock *) +| ROHF (** Restricted Open-shell Hartree-Fock *) +| UHF (** Unrestricted Hartree-Fock *) + +and t = + { + kind : hartree_fock_kind; + simulation : Simulation.t; + guess : Guess.t; + data : hartree_fock_data option lazy_t array; + nocc : int ; + } + + +let empty = + { + iteration = 0 ; + coefficients = None ; + eigenvalues = None ; + error = None ; + diis = None ; + energy = None ; + density = None ; + density_a = None ; + density_b = None ; + fock = None ; + fock_a = None ; + fock_b = None ; + } + + + + + + +module Si = Simulation +module El = Electrons + + +let kind t = t.kind +let simulation t = t.simulation +let guess t = t.guess +let nocc t = t.nocc + + + +let n_iterations t = + Array.fold_left (fun accu x -> + match Lazy.force x with + | Some _ -> accu + 1 + | None -> accu + ) 0 t.data + + +let last_iteration t = + Util.of_some @@ Lazy.force (t.data.(n_iterations t - 1)) + +let eigenvectors t = + let data = last_iteration t in + Util.of_some data.coefficients + +let eigenvalues t = + let data = last_iteration t in + Util.of_some data.eigenvalues + +let density t = + let data = last_iteration t in + match kind t with + | RHF -> Util.of_some data.density + | ROHF -> Matrix.add (Util.of_some data.density_a) (Util.of_some data.density_b) + | _ -> failwith "Not implemented" + +let occupation t = + let n_alfa, n_beta = + El.n_alfa @@ Simulation.electrons @@ simulation t, + El.n_beta @@ Simulation.electrons @@ simulation t + in + match kind t with + | RHF -> Vector.init (Matrix.dim2 @@ eigenvectors t) (fun i -> + if i <= nocc t then 2.0 else 0.0) + | ROHF -> Vector.init (Matrix.dim2 @@ eigenvectors t) (fun i -> + if i <= n_beta then 2.0 else + if i <= n_alfa then 1.0 else + 0.0) + | _ -> failwith "Not implemented" + + +let energy t = + let data = last_iteration t in + Util.of_some data.energy + + +let nuclear_repulsion t = + Si.nuclear_repulsion (simulation t) + + +let ao_basis t = + Si.ao_basis (simulation t) + + +let kin_energy t = + let m_T = + ao_basis t + |> Ao.Basis.kin_ints + in + let m_P = density t in + Matrix.gemm_trace m_P m_T + + +let eN_energy t = + let m_V = + ao_basis t + |> Ao.Basis.eN_ints + in + let m_P = density t in + Matrix.gemm_trace m_P m_V + + +let coulomb_energy t = + let data = + last_iteration t + in + match kind t with + | RHF -> let m_P = Util.of_some data.density in + let fock = Util.of_some data.fock in + let m_J = Fock.coulomb fock in + 0.5 *. Matrix.gemm_trace m_P m_J + + | ROHF -> let m_P_a = Util.of_some data.density_a in + let m_P_b = Util.of_some data.density_b in + let fock_a = Util.of_some data.fock_a in + let fock_b = Util.of_some data.fock_b in + let m_J_a = Fock.coulomb fock_a in + let m_J_b = Fock.coulomb fock_b in + 0.5 *. ( (Matrix.gemm_trace m_P_a m_J_a) +. (Matrix.gemm_trace m_P_b m_J_b) ) + + | _ -> failwith "Not implemented" + + +let exchange_energy t = + let data = + last_iteration t + in + match kind t with + | RHF -> let m_P = Util.of_some data.density in + let fock = Util.of_some data.fock in + let m_K = Fock.exchange fock in + -. 0.5 *. Matrix.gemm_trace m_P m_K + + | ROHF -> let m_P_a = Util.of_some data.density_a in + let m_P_b = Util.of_some data.density_b in + let fock_a = Util.of_some data.fock_a in + let fock_b = Util.of_some data.fock_b in + let m_K_a = Fock.exchange fock_a in + let m_K_b = Fock.exchange fock_b in + -. 0.5 *. ( (Matrix.gemm_trace m_P_a m_K_a) +. (Matrix.gemm_trace m_P_b m_K_b) ) + + | _ -> failwith "Not implemented" + + + + + + +let make + ?kind + ?guess:(guess=`Huckel) + ?max_scf:(max_scf=64) + ?level_shift:(level_shift=0.2) + ?threshold_SCF:(threshold_SCF=1.e-8) + simulation = + + + (* Number of occupied MOs *) + let n_alfa, n_beta = + El.n_alfa @@ Si.electrons simulation, + El.n_beta @@ Si.electrons simulation + in + + let nocc = n_alfa in + + let kind = + match kind with + | Some kind -> kind + | None -> if (n_alfa = n_beta) then RHF else ROHF + in + + let nuclear_repulsion = + Si.nuclear_repulsion simulation + in + + let ao_basis = + Si.ao_basis simulation + in + + + (* Orthogonalization matrix *) + let m_X = + Ao.Basis.ortho ao_basis + in + + (* Overlap matrix *) + let m_S = + Ao.Basis.overlap ao_basis + in + + (* Level shift in MO basis *) + let m_LSmo = + Array.init (Matrix.dim2 m_X) (fun i -> + if i > nocc then level_shift else 0.) + |> Vector.of_array + |> Matrix.of_diag + in + + (* Guess coefficients *) + let guess = + Guess.make ~nocc ~guess ao_basis + in + + let m_C : (ao,mo) Matrix.t = + let c_of_h m_H = + let m_Hmo = Matrix.xt_o_x ~o:m_H ~x:m_X in + let m_C', _ = Matrix.diagonalize_symm m_Hmo in + Matrix.gemm_nn m_X m_C' + in + match guess with + | Guess.Hcore m_H -> c_of_h m_H + | Guess.Huckel m_H -> c_of_h m_H + | Guess.Matrix m_C -> m_C + in + + (* A single SCF iteration *) + let scf_iteration_rhf data = + + let nSCF = data.iteration + 1 + and m_C = Util.of_some data.coefficients + and m_P_prev = data.density + and fock_prev = data.fock + and diis = + match data.diis with + | Some diis -> diis + | None -> Diis.make () + and threshold = + match data.error with + | Some error -> error + | None -> threshold_SCF *. 2. + in + + (* Density matrix over nocc occupied MOs *) + let m_P = + Matrix.gemm_nt ~alpha:2. ~k:nocc m_C m_C + in + + (* Fock matrix in AO basis *) + let fock = + match fock_prev, m_P_prev, threshold > 100. *. threshold_SCF with + | Some fock_prev, Some m_P_prev, true -> + let threshold = 1.e-8 in + Fock.make_rhf ~density:(Matrix.sub m_P m_P_prev) ~threshold ao_basis + |> Fock.add fock_prev + | _ -> Fock.make_rhf ~density:m_P ao_basis + in + + let m_F0, m_Hc, _, _ = + let x = fock in + Fock.(fock x, core x, coulomb x, exchange x) + in + + (* Add level shift in AO basis *) + let m_F = + let m_SC = + Matrix.gemm_nn m_S m_C + in + Matrix.gemm_nn m_SC (Matrix.gemm_nt m_LSmo m_SC) + |> Matrix.add m_F0 + in + + + (* Fock matrix in orthogonal basis *) + let m_F_ortho = + Matrix.xt_o_x ~o:m_F ~x:m_X + in + + let error_fock = + let fps = + Matrix.gemm_nn m_F (Matrix.gemm_nn m_P m_S) + and spf = + Matrix.gemm_nn m_S (Matrix.gemm_nn m_P m_F) + in + Matrix.xt_o_x ~o:(Matrix.sub fps spf) ~x:m_X + in + + let diis, m_F_diis = + let diis = + Diis.append + ~p:(Matrix.as_vec_inplace m_F_ortho) + ~e:(Matrix.as_vec_inplace error_fock) diis + in + + try + let m_F_diis = + Diis.next diis + |> Matrix.reshape_vec_inplace (Matrix.dim1 m_F_ortho) (Matrix.dim2 m_F_ortho) + in + diis, m_F_diis + + with Failure _ -> (* Failure in DIIS.next *) + Diis.make (), m_F_ortho + in + let diis = + Diis.append + ~p:(Matrix.as_vec_inplace m_F_ortho) + ~e:(Matrix.as_vec_inplace error_fock) diis + in + + + (* MOs in orthogonal MO basis *) + let m_C', _ = + Matrix.diagonalize_symm m_F_diis + in + + (* Re-compute eigenvalues to remove level-shift *) + let eigenvalues = + let m_F_ortho = + Matrix.xt_o_x ~o:m_F0 ~x:m_X + in + Matrix.xt_o_x ~o:m_F_ortho ~x:m_C' + |> Matrix.diag + in + + (* MOs in AO basis *) + let m_C = + Matrix.gemm_nn m_X m_C' + |> Conventions.rephase + in + + (* Hartree-Fock energy *) + let energy = + nuclear_repulsion +. 0.5 *. + Matrix.gemm_trace m_P (Matrix.add m_Hc m_F) + in + + (* Convergence criterion *) + let error = + error_fock + |> Matrix.amax + |> abs_float + in + + { empty with + iteration = nSCF ; + eigenvalues = Some eigenvalues ; + coefficients = Some m_C ; + error = Some error ; + diis = Some diis ; + energy = Some energy ; + density = Some m_P ; + fock = Some fock ; + } + + in + + let scf_iteration_rohf data = + + let nSCF = data.iteration + 1 + and m_C = Util.of_some data.coefficients + and m_P_a_prev = data.density_a + and m_P_b_prev = data.density_b + and fock_a_prev = data.fock_a + and fock_b_prev = data.fock_b + and diis = + match data.diis with + | Some diis -> diis + | None -> Diis.make () + and threshold = + match data.error with + | Some error -> error + | None -> threshold_SCF *. 2. + in + + (* Density matrix *) + let m_P_a = + Matrix.gemm_nt ~alpha:1. ~k:n_alfa m_C m_C + in + + let m_P_b = + Matrix.gemm_nt ~alpha:1. ~k:n_beta m_C m_C + in + + let m_P = + Matrix.add m_P_a m_P_b + in + + (* Fock matrix in AO basis *) + let fock_a = + match fock_a_prev, threshold > 100. *. threshold_SCF with + | Some fock_a_prev, true -> + let threshold = 1.e-8 in + Fock.make_uhf ~density_same:(Matrix.sub m_P_a @@ Util.of_some m_P_a_prev) ~density_other:(Matrix.sub m_P_b @@ Util.of_some m_P_b_prev) ~threshold ao_basis + |> Fock.add fock_a_prev + | _ -> Fock.make_uhf ~density_same:m_P_a ~density_other:m_P_b ao_basis + in + + let fock_b = + match fock_b_prev, threshold > 100. *. threshold_SCF with + | Some fock_b_prev, true -> + let threshold = 1.e-8 in + Fock.make_uhf ~density_same:(Matrix.sub m_P_b @@ Util.of_some m_P_b_prev) ~density_other:(Matrix.sub m_P_a @@ Util.of_some m_P_a_prev) ~threshold ao_basis + |> Fock.add fock_b_prev + | _ -> Fock.make_uhf ~density_same:m_P_b ~density_other:m_P_a ao_basis + in + + let m_F_a, m_Hc_a, _, _ = + let x = fock_a in + Fock.(fock x, core x, coulomb x, exchange x) + in + + let m_F_b, m_Hc_b, _, _ = + let x = fock_b in + Fock.(fock x, core x, coulomb x, exchange x) + in + + + let m_F_mo_a = + Matrix.xt_o_x ~o:m_F_a ~x:m_C + in + + let m_F_mo_b = + Matrix.xt_o_x ~o:m_F_b ~x:m_C + in + + let m_F_mo = + let space k = + if k <= n_beta then + `Core + else if k <= n_alfa then + `Active + else + `Virtual + in + Array.init (Matrix.dim2 m_F_mo_a) (fun i -> + let i = i+1 in + Array.init (Matrix.dim1 m_F_mo_a) (fun j -> + let j = j+1 in + match (space i), (space j) with + | `Core , `Core -> + 0.5 *. (m_F_mo_a%:(i,j) +. m_F_mo_b%:(i,j)) -. + (m_F_mo_b%:(i,j) -. m_F_mo_a%:(i,j)) + + | `Active , `Core + | `Core , `Active -> + m_F_mo_b%:(i,j) + + | `Core , `Virtual + | `Virtual , `Core + | `Active , `Active -> + 0.5 *. (m_F_mo_a%:(i,j) +. m_F_mo_b%:(i,j)) + + | `Virtual , `Active + | `Active , `Virtual -> + m_F_mo_a%:(i,j) + + | `Virtual , `Virtual -> + 0.5 *. (m_F_mo_a%:(i,j) +. m_F_mo_b%:(i,j)) +. + (m_F_mo_b%:(i,j) -. m_F_mo_a%:(i,j)) + ) ) + |> Matrix.of_array + in + + let m_SC = + Matrix.gemm_nn m_S m_C + in + + let m_F0 = + Matrix.x_o_xt ~x:m_SC ~o:m_F_mo + in + + + (* Add level shift in AO basis *) + let m_F = + Matrix.x_o_xt ~x:m_SC ~o:m_LSmo + |> Matrix.add m_F0 + in + + (* Fock matrix in orthogonal basis *) + let m_F_ortho = + Matrix.xt_o_x ~o:m_F ~x:m_X + in + + let error_fock = + let fps = + Matrix.gemm_nn m_F (Matrix.gemm_nn m_P m_S) + and spf = + Matrix.gemm_nn m_S (Matrix.gemm_nn m_P m_F) + in + Matrix.xt_o_x ~o:(Matrix.sub fps spf) ~x:m_X + in + + let diis, m_F_diis = + let diis = + Diis.append + ~p:(Matrix.as_vec_inplace m_F_ortho) + ~e:(Matrix.as_vec_inplace error_fock) diis + in + + try + let m_F_diis = + Diis.next diis + |> Matrix.reshape_vec_inplace (Matrix.dim1 m_F_ortho) (Matrix.dim2 m_F_ortho) + in + diis, m_F_diis + + with Failure _ -> (* Failure in DIIS.next *) + Diis.make (), m_F_ortho + in + + + (* MOs in orthogonal MO basis *) + let m_C', _ = + Matrix.diagonalize_symm m_F_diis + in + + (* Re-compute eigenvalues to remove level-shift *) + let eigenvalues = + let m_F_ortho = + Matrix.xt_o_x ~o:m_F0 ~x:m_X + in + Matrix.xt_o_x ~o:m_F_ortho ~x:m_C' + |> Matrix.diag + in + + (* MOs in AO basis *) + let m_C = + Matrix.gemm_nn m_X m_C' + |> Conventions.rephase + in + + (* Hartree-Fock energy *) + let energy = + nuclear_repulsion +. 0.5 *. ( Matrix.gemm_trace m_P_a (Matrix.add m_Hc_a m_F_a) +. + Matrix.gemm_trace m_P_b (Matrix.add m_Hc_b m_F_b) ) + in + + (* Convergence criterion *) + let error = + error_fock + |> Matrix.amax + |> abs_float + in + { empty with + iteration = nSCF ; + eigenvalues = Some eigenvalues ; + coefficients = Some m_C ; + error = Some error ; + diis = Some diis ; + energy = Some energy ; + density_a = Some m_P_a ; + density_b = Some m_P_b ; + fock_a = Some fock_a ; + fock_b = Some fock_b ; + } + + in + + + let scf_iteration data = + match kind with + | RHF -> scf_iteration_rhf data + | ROHF -> scf_iteration_rohf data + | _ -> failwith "Not implemented" + in + + let array_data = + + let storage = + Array.make max_scf None + in + + let rec iteration = function + | 0 -> Some (scf_iteration { empty with coefficients = Some m_C }) + | n -> begin + match storage.(n) with + | Some result -> Some result + | None -> + begin + let data = iteration (n-1) in + match data with + | None -> None + | Some data -> + begin + (* Check convergence *) + let converged, _error = + match data.error with + | None -> false, 0. + | Some error -> (data.iteration = max_scf || error < threshold_SCF), error + in + if converged then + None + else + begin + storage.(n-1) <- Some { empty with + iteration = data.iteration; + energy = data.energy ; + eigenvalues = data.eigenvalues ; + error = data.error ; + }; + storage.(n) <- Some (scf_iteration data); + storage.(n); + end + end + end + end + in + Array.init max_scf (fun i -> lazy (iteration i)) + in + { + kind; + simulation; + guess ; + data = array_data; + nocc; + } + + + + + + + + + + +let linewidth = 60 + +let pp_iterations ppf t = + let line = (String.make linewidth '-') in + Format.fprintf ppf "@[%4s%s@]@." "" line; + Format.fprintf ppf "@[%4s@[%5s@]@,@[%16s@]@,@[%16s@]@,@[%11s@]@]@." + "" "#" "HF energy " "Convergence" "HOMO-LUMO"; + Format.fprintf ppf "@[%4s%s@]@." "" line; + let nocc = nocc t in + Array.iter (fun data -> + let data = Lazy.force data in + match data with + | None -> () + | Some data -> + let e = Util.of_some data.eigenvalues in + let gap = e%.(nocc+1) -. e%.(nocc) in + begin + Format.fprintf ppf "@[%4s@[%5d@]@,@[%16.8f@]@,@[%16.4e@]@,@[%11.4f@]@]@." "" + (data.iteration) (Util.of_some data.energy) (Util.of_some data.error) gap; + end + ) t.data; + Format.fprintf ppf "@[%4s%s@]@." "" line + + +let pp_summary ppf t = + let print text value = + Format.fprintf ppf "@[@[%30s@]@,@[%16.10f@]@]@;" text value; + and line () = + Format.fprintf ppf "@[ %s @]@;" (String.make (linewidth-4) '-'); + in + let ha_to_ev = Constants.ha_to_ev in + let e = eigenvalues t in + + Format.fprintf ppf "@[%s@]@;" (String.make linewidth '=') + ; Format.fprintf ppf "@[" + ; print "One-electron energy" (kin_energy t +. eN_energy t) + ; print "Kinetic" (kin_energy t) + ; print "Potential" (eN_energy t) + ; line () + ; print "Two-electron energy" (coulomb_energy t +. exchange_energy t) + ; print "Coulomb" (coulomb_energy t) + ; print "Exchange" (exchange_energy t) + ; line () + ; print "HF HOMO" (ha_to_ev *. (e%.(nocc t))) + ; print "HF LUMO" (ha_to_ev *. (e%.(nocc t + 1))) + ; print "HF LUMO-LUMO" (ha_to_ev *. (e%.(nocc t + 1) -. e%.(nocc t))) + ; line () + ; print "Electronic energy" (energy t -. nuclear_repulsion t) + ; print "Nuclear repulsion" (nuclear_repulsion t) + ; print "Hartree-Fock energy" (energy t) + ; Format.fprintf ppf "@]" + ; Format.fprintf ppf "@[%s@]@;" (String.make linewidth '=') + + + + +let pp ppf t = + Format.fprintf ppf "@.@[%s@]@." (String.make 70 '='); + Format.fprintf ppf "@[%34s %-34s@]@." (match t.kind with + | UHF -> "Unrestricted" + | RHF -> "Restricted" + | ROHF -> "Restricted Open-shell") "Hartree-Fock"; + Format.fprintf ppf "@[%s@]@.@." (String.make 70 '='); + Format.fprintf ppf "@[%a@]@." pp_iterations t; + Format.fprintf ppf "@[@;%a@]@." pp_summary t diff --git a/mo/lib/hartree_fock.mli b/mo/lib/hartree_fock.mli new file mode 100644 index 0000000..0c2f6cb --- /dev/null +++ b/mo/lib/hartree_fock.mli @@ -0,0 +1,75 @@ +open Linear_algebra + +(** Data structure representing the output of a Hartree-Fock caculation *) + +type hartree_fock_data + +type hartree_fock_kind = +| RHF (** Restricted Hartree-Fock *) +| ROHF (** Restricted Open-shell Hartree-Fock *) +| UHF (** Unrestricted Hartree-Fock *) + + +type t +type mo = Mo_dim.t +type ao = Ao.Ao_dim.t + +val kind : t -> hartree_fock_kind +(** Kind of simulation : RHF, ROHF or UHF. *) + +val simulation : t -> Simulation.t +(** Simulation which was used for HF calculation *) + +val guess : t -> Guess.t +(** Initial guess *) + +val eigenvectors : t -> (ao, mo) Matrix.t +(** Final eigenvectors *) + +val eigenvalues : t -> mo Vector.t +(** Final eigenvalues *) + +val occupation : t -> mo Vector.t +(** Diagonal of the density matrix *) + +val energy : t -> float +(** Final energy *) + +val nuclear_repulsion : t -> float +(** Nucleus-Nucleus potential energy *) + +val kin_energy : t -> float +(** Kinetic energy *) + +val eN_energy : t -> float +(** Electron-nucleus potential energy *) + +val coulomb_energy : t -> float +(** Electron-Electron potential energy *) + +val exchange_energy : t -> float +(** Exchange energy *) + +val nocc : t -> int +(** Number of occupied MOs *) + +val empty: hartree_fock_data +(** Empty data *) + + +val make : + ?kind:hartree_fock_kind -> + ?guess:[ `Hcore | `Huckel | `Matrix of (ao,mo) Matrix.t ] -> + ?max_scf:int -> + ?level_shift:float -> ?threshold_SCF:float -> + Simulation.t -> t + + +(** {1 Printers} *) + +val pp : Format.formatter -> t -> unit + +val pp_iterations : Format.formatter -> t -> unit + +val pp_summary : Format.formatter -> t -> unit + diff --git a/mo/lib/mo_dim.ml b/mo/lib/mo_dim.ml new file mode 100644 index 0000000..63c57c4 --- /dev/null +++ b/mo/lib/mo_dim.ml @@ -0,0 +1 @@ +type t diff --git a/mo/lib/mo_dim.mli b/mo/lib/mo_dim.mli new file mode 100644 index 0000000..63c57c4 --- /dev/null +++ b/mo/lib/mo_dim.mli @@ -0,0 +1 @@ +type t diff --git a/mo/test/dune b/mo/test/dune new file mode 100644 index 0000000..16aec77 --- /dev/null +++ b/mo/test/dune @@ -0,0 +1,12 @@ +(library + (name test_mo_basis) + (libraries + alcotest + qcaml.mo + ) + (synopsis "Tests for the MO basis")) + + + + + diff --git a/mo/test/guess.ml b/mo/test/guess.ml new file mode 100644 index 0000000..313394d --- /dev/null +++ b/mo/test/guess.ml @@ -0,0 +1,50 @@ +open Alcotest +open Particles +open Operators +open Mo.Guess +open Linear_algebra + +let wd = Common.Qcaml.root ^ Filename.dir_sep ^ "test" + +let test_case ao_basis = + + let test_hcore () = + match make ~guess:`Hcore ao_basis with + | Hcore matrix -> + let a = Matrix.to_array matrix in + let reference = + Matrix.add + (Ao.Basis.eN_ints ao_basis) + (Ao.Basis.kin_ints ao_basis) + |> Matrix.to_array + in + Array.iteri (fun i x -> + let message = + Printf.sprintf "Guess line %d" (i) + in + check (array (float 1.e-15)) message a.(i) x) reference + | _ -> assert false + + in + [ + "HCore", `Quick, test_hcore; + ] + +let tests = + List.concat [ + let nuclei = + wd ^ Filename.dir_sep ^ "water.xyz" + |> Nuclei.of_xyz_file + in + let basis_filename = + wd ^ Filename.dir_sep ^ "cc-pvdz" + in + let rs = 0.5 in + let operators = [ Operator.of_range_separation rs ] in + let ao_basis = + Ao.Basis.of_nuclei_and_basis_filename ~kind:`Gaussian + ~operators ~cartesian:false ~nuclei basis_filename + in + test_case ao_basis ; + ] + diff --git a/test/dune b/test/dune index 2b0d37a..e0c8096 100644 --- a/test/dune +++ b/test/dune @@ -8,6 +8,7 @@ test_gaussian_basis test_gaussian_integrals test_ao_basis + test_mo_basis )) (alias diff --git a/test/run_tests.ml b/test/run_tests.ml index 4783a91..fecebd5 100644 --- a/test/run_tests.ml +++ b/test/run_tests.ml @@ -11,6 +11,7 @@ let test_suites: unit Alcotest.test list = [ "Gaussian_basis.General_basis", Test_gaussian_basis.General_basis.tests; "Ao_basis.Ao_basis_gaussian", Test_ao_basis.Ao_basis_gaussian.tests; "Ao_basis.Ao_basis", Test_ao_basis.Ao_basis.tests; + "Mo_basis.Guess", Test_mo_basis.Guess.tests; ] let () = Alcotest.run "QCaml" test_suites