diff --git a/gaussian_integrals/lib/electron_nucleus.ml b/gaussian_integrals/lib/electron_nucleus.ml index 022db9b..7ca56fd 100644 --- a/gaussian_integrals/lib/electron_nucleus.ml +++ b/gaussian_integrals/lib/electron_nucleus.ml @@ -51,7 +51,6 @@ let of_basis_nuclei ~basis nuclei = in let eni_array = Matrix.create n n in - let eni_array_x = Matrix.to_bigarray_inplace eni_array in (* Pre-compute all shell pairs *) let shell_pairs = @@ -83,8 +82,8 @@ let of_basis_nuclei ~basis nuclei = let value = Zmap.find cls key in - eni_array_x.{j_c,i_c} <- value; - eni_array_x.{i_c,j_c} <- value; + Matrix.set eni_array j_c i_c value; + Matrix.set eni_array i_c j_c value; ) (Am.zkey_array (Singlet (Cs.ang_mom shell.(j)))) ) (Am.zkey_array (Singlet (Cs.ang_mom shell.(i)))) done; @@ -96,11 +95,10 @@ let of_basis_nuclei ~basis nuclei = let to_file ~filename eni_array = let n = Matrix.dim1 eni_array in let oc = open_out filename in - let eni_array_x = Matrix.to_bigarray_inplace eni_array in for j=1 to n do for i=1 to j do - let value = eni_array_x.{i,j} in + let value = eni_array%:(i,j) in if (value <> 0.) then Printf.fprintf oc " %5d %5d %20.15f\n" i j value; done; diff --git a/gaussian_integrals/lib/kinetic.ml b/gaussian_integrals/lib/kinetic.ml index 5bbc75d..9bd9849 100644 --- a/gaussian_integrals/lib/kinetic.ml +++ b/gaussian_integrals/lib/kinetic.ml @@ -134,7 +134,6 @@ let of_basis basis = in let result = Matrix.create n n in - let result_x = Matrix.to_bigarray_inplace result in for j=0 to (Array.length shell) - 1 do for i=0 to j do (* Compute all the integrals of the class *) @@ -155,8 +154,8 @@ let of_basis basis = try Zmap.find cls key with Not_found -> 0. in - result_x.{i_c,j_c} <- value; - result_x.{j_c,i_c} <- value; + Matrix.set result i_c j_c value; + Matrix.set result j_c i_c value; ) (Am.zkey_array (Singlet (Cs.ang_mom shell.(i)))) ) (Am.zkey_array (Singlet (Cs.ang_mom shell.(j)))) done; @@ -175,11 +174,10 @@ let to_file ~filename kinetic = Matrix.dim1 kinetic in - let kinetic_x = Matrix.to_bigarray_inplace kinetic in for j=1 to n do for i=1 to j do - if (abs_float kinetic_x.{i,j} > cutoff) then - Printf.fprintf oc "%4d %4d %20.12e\n" i j kinetic_x.{i,j} + if (abs_float (kinetic%:(i,j)) > cutoff) then + Printf.fprintf oc "%4d %4d %20.12e\n" i j (kinetic%:(i,j)) done; done; close_out oc diff --git a/gaussian_integrals/lib/orthonormalization.ml b/gaussian_integrals/lib/orthonormalization.ml index 1d1c1c6..dcc5b9b 100644 --- a/gaussian_integrals/lib/orthonormalization.ml +++ b/gaussian_integrals/lib/orthonormalization.ml @@ -50,7 +50,6 @@ let make_lowdin ~thresh ~overlap = let u_vec, u_val = Matrix.diagonalize_symm overlap in - let u_vec_x = Matrix.to_bigarray_inplace u_vec in Vector.iter (fun x -> if x < thresh then invalid_arg (__FILE__^": make_lowdin") ) u_val; @@ -59,7 +58,7 @@ let make_lowdin ~thresh ~overlap = let u_vec' = Matrix.init_cols (Matrix.dim1 u_vec) (Matrix.dim2 u_vec) - (fun i j -> u_vec_x.{i,j} *. (Vector.at u_val j)) + (fun i j -> u_vec%:(i,j) *. (u_val%.(j)) ) in Matrix.gemm_nt u_vec' u_vec diff --git a/gaussian_integrals/lib/overlap.ml b/gaussian_integrals/lib/overlap.ml index b22bd89..747a2a7 100644 --- a/gaussian_integrals/lib/overlap.ml +++ b/gaussian_integrals/lib/overlap.ml @@ -99,7 +99,6 @@ let of_basis basis = in let result = Matrix.create n n in - let result_x = Matrix.to_bigarray_inplace result in for j=0 to (Array.length shell) - 1 do for i=0 to j do (* Compute all the integrals of the class *) @@ -120,8 +119,8 @@ let of_basis basis = try Zmap.find cls key with Not_found -> 0. in - result_x.{i_c,j_c} <- value; - result_x.{j_c,i_c} <- value; + Matrix.set result i_c j_c value; + Matrix.set result j_c i_c value; ) (Am.zkey_array (Singlet (Cs.ang_mom shell.(i)))) ) (Am.zkey_array (Singlet (Cs.ang_mom shell.(j)))) done; @@ -146,7 +145,6 @@ let of_basis_pair first_basis second_basis = in let result = Matrix.create n m in - let result_x = Matrix.to_bigarray_inplace result in for j=0 to (Array.length second) - 1 do for i=0 to (Array.length first) - 1 do (* Compute all the integrals of the class *) @@ -167,7 +165,7 @@ let of_basis_pair first_basis second_basis = try Zmap.find cls key with Not_found -> 0. in - result_x.{i_c,j_c} <- value; + Matrix.set result i_c j_c value; ) (Am.zkey_array (Singlet (Cs.ang_mom first.(i)))) ) (Am.zkey_array (Singlet (Cs.ang_mom second.(j)))) done; @@ -184,11 +182,10 @@ let to_file ~filename overlap = Matrix.dim1 overlap in - let overlap_x = Matrix.to_bigarray_inplace overlap in for j=1 to n do for i=1 to j do - if (abs_float overlap_x.{i,j} > cutoff) then - Printf.fprintf oc "%4d %4d %20.12e\n" i j overlap_x.{i,j} + if (abs_float (overlap%:(i,j)) > cutoff) then + Printf.fprintf oc "%4d %4d %20.12e\n" i j (overlap%:(i,j)) done; done; close_out oc diff --git a/gaussian_integrals/lib/two_electron_rr_vectorized.ml b/gaussian_integrals/lib/two_electron_rr_vectorized.ml index 789d9ea..35e0cf3 100644 --- a/gaussian_integrals/lib/two_electron_rr_vectorized.ml +++ b/gaussian_integrals/lib/two_electron_rr_vectorized.ml @@ -638,8 +638,8 @@ let contracted_class_shell_pairs ?operator ~zero_m ?schwartz_p ?schwartz_q shell raise NullQuartet; let expo_p_inv, expo_q_inv = - (Vector.at expo_p_inv i), - (Vector.at expo_q_inv j) + (expo_p_inv%.(i)), + (expo_q_inv%.(j)) in let center_pq = diff --git a/linear_algebra/lib/conventions.ml b/linear_algebra/lib/conventions.ml index 234394d..682f429 100644 --- a/linear_algebra/lib/conventions.ml +++ b/linear_algebra/lib/conventions.ml @@ -5,18 +5,20 @@ *) +open Vector + let in_phase vec = let s = Vector.sum vec in if s = 0. then let rec first_non_zero k = if k > Vector.dim vec then k-1 - else if Vector.at vec k = 0. then + else if vec%.(k) = 0. then first_non_zero (k+1) else k in let k = first_non_zero 1 in - Vector.at vec k >= 0. + vec%.(k) >= 0. else s > 0. diff --git a/linear_algebra/lib/matrix.ml b/linear_algebra/lib/matrix.ml index 641119c..0e8af4a 100644 --- a/linear_algebra/lib/matrix.ml +++ b/linear_algebra/lib/matrix.ml @@ -41,8 +41,6 @@ let add_const_inplace x a = let add_const x a = Mat.add_const x a -let at t i j = t.{i,j} - external to_bigarray_inplace : ('a,'b) t -> (float, Stdlib.Bigarray.float64_elt, Stdlib.Bigarray.fortran_layout) Stdlib.Bigarray.Array2.t = "%identity" external of_bigarray_inplace : (float, Stdlib.Bigarray.float64_elt, Stdlib.Bigarray.fortran_layout) Stdlib.Bigarray.Array2.t -> ('a,'b) t = "%identity" @@ -108,6 +106,43 @@ let scale_inplace x t = let scale x t = out_of_place (fun t -> scale_inplace x t) t +let of_diag v = + Vector.to_bigarray_inplace v + |> Mat.of_diag + +let diag t = + Mat.copy_diag t + |> Vector.of_bigarray_inplace + +let gemv_n_inplace ?m ?n ?(beta=0.) y ?(alpha=1.) ?(ar=1) ?(ac=1) t v = + let y = Vector.to_bigarray_inplace y in + let v = Vector.to_bigarray_inplace v in + ignore @@ gemv ?m ?n ~beta ~trans:`N ~y ~alpha ~ar ~ac t v + +let gemv_t_inplace ?m ?n ?(beta=0.) y ?(alpha=1.) ?(ar=1) ?(ac=1) t v = + let y = Vector.to_bigarray_inplace y in + let v = Vector.to_bigarray_inplace v in + ignore @@ gemv ?m ?n ~beta ~trans:`T ~y ~alpha ~ar ~ac t v + +let gemv_n ?m ?n ?(beta=0.) ?y ?(alpha=1.) ?(ar=1) ?(ac=1) t v = + let v = Vector.to_bigarray_inplace v in + let y = + match y with + | None -> None + | Some y -> Some (Vector.to_bigarray_inplace y) + in + gemv ?m ?n ~beta ?y ~trans:`N ~alpha ~ar ~ac t v + |> Vector.of_bigarray_inplace + +let gemv_t ?m ?n ?(beta=0.) ?y ?(alpha=1.) ?(ar=1) ?(ac=1) t v = + let v = Vector.to_bigarray_inplace v in + let y = + match y with + | None -> None + | Some y -> Some (Vector.to_bigarray_inplace y) + in + gemv ?m ?n ~beta ?y ~trans:`T ~alpha ~ar ~ac t v + |> Vector.of_bigarray_inplace let gemm_inplace ?m ?n ?k ?(beta=0.) ~c ?(transa=`N) ?(alpha=1.0) a ?(transb=`N) b = ignore @@ gemm ?m ?n ?k ~beta ~c ~transa ~alpha a ~transb b @@ -312,3 +347,9 @@ let qr a = orgqr ~tau result; let q = result in q, r + +let (%:) t (i,j) = t.{i,j} + +let set t i j v = t.{i,j} <- v + + diff --git a/linear_algebra/lib/matrix.mli b/linear_algebra/lib/matrix.mli index c24ec21..fb74db1 100644 --- a/linear_algebra/lib/matrix.mli +++ b/linear_algebra/lib/matrix.mli @@ -26,6 +26,12 @@ val init_cols : int -> int -> (int -> int -> float) -> ('a,'b) t val identity: int -> ('a,'b) t (** Creates an identity matrix. *) +val of_diag: 'a Vector.t -> ('a,'a) t +(** Creates a diagonal matrix. *) + +val diag: ('a,'a) t -> 'a Vector.t +(** Returns the diagonal of a matrix. *) + val fill_inplace: ('a,'b) t -> float -> unit (** Fills the matrix with the give value. *) @@ -65,9 +71,6 @@ val mul_inplace : c:('a,'b) t -> ('a,'b) t -> ('a,'b) t -> unit val div_inplace : c:('a,'b) t -> ('a,'b) t -> ('a,'b) t -> unit (** [div_inplace c a b] : performs [c = a/b] element-wise in-place. *) -val at : ('a,'b) t -> int -> int -> float -(** [at i j] returns the element at i,j. *) - (* val to_bigarray : ('a,'b) t -> (float, Stdlib.Bigarray.float64_elt, Stdlib.Bigarray.fortran_layout) Stdlib.Bigarray.Array2.t (** Converts the matrix into a Bigarray in Fortran layout *) @@ -156,6 +159,39 @@ val outer_product_inplace : ('a,'b) t -> ?alpha:float -> 'a Vector.t -> 'b Vecto (** Computes M = %{ $\alpha u.v^t$ %} *) +val gemv_n_inplace : ?m:int -> ?n:int -> ?beta:float -> 'a Vector.t -> + ?alpha:float -> ?ar:int -> ?ac:int -> ('a,'b) t -> 'b Vector.t -> + unit +(** Performs the Lapack GEMV operation. Default values: +[beta=0.] [alpha=1.0]. +[gemv ~beta y ~alpha m v]: %{ $Y = \beta Y + \alpha M V$ +The vector Y is updated in-place. +*) + +val gemv_t_inplace : ?m:int -> ?n:int -> ?beta:float -> 'b Vector.t -> + ?alpha:float -> ?ar:int -> ?ac:int -> ('a,'b) t -> 'a Vector.t -> + unit +(** Performs the Lapack GEMV operation. Default values: +[beta=0.] [alpha=1.0]. +[gemv ~beta y ~alpha m v]: %{ $Y = \beta Y + \alpha M^\dagger V$ +The vector Y is updated in-place. +*) + +val gemv_n : ?m:int -> ?n:int -> ?beta:float -> ?y:'a Vector.t -> + ?alpha:float -> ?ar:int -> ?ac:int -> ('a,'b) t -> 'b Vector.t -> + 'a Vector.t +(** Performs the Lapack GEMV operation. Default values: +[beta=0.] [alpha=1.0]. +[gemv ~beta y ~alpha m v]: %{ $Y = \beta Y + \alpha M^\dagger V$ *) + +val gemv_t : ?m:int -> ?n:int -> ?beta:float -> ?y:'b Vector.t -> + ?alpha:float -> ?ar:int -> ?ac:int -> ('a,'b) t -> 'a Vector.t -> + 'b Vector.t +(** Performs the Lapack GEMV operation. Default values: +[beta=0.] [alpha=1.0]. +[gemv ~beta y ~alpha m v]: %{ $Y = \beta Y + \alpha M^\dagger V$ +*) + val gemm_inplace : ?m:int -> ?n:int -> ?k:int -> ?beta:float -> c:('a,'b) t -> ?transa:[`N | `T] -> ?alpha:float -> ('c,'d) t -> ?transb:[`N | `T] -> ('e,'f) t -> unit @@ -265,5 +301,11 @@ val sysv_inplace : b:('a,'b) t -> ('a,'a) t -> unit val sysv : b:('a,'b) t -> ('a,'a) t -> ('a,'b) t (** Solves %{ $AX=B$ %} when A is symmetric *) +val (%:) : ('a,'b) t -> (int*int) -> float +(** [t%.(i,j)] returns the element at i,j. *) + +val set : ('a,'b) t -> int -> int -> float -> unit +(** [set t i j v] sets the (i,j)-th element to v *) + val pp : Format.formatter -> ('a,'b) t -> unit diff --git a/linear_algebra/lib/vector.ml b/linear_algebra/lib/vector.ml index fa36f8b..47b68a5 100644 --- a/linear_algebra/lib/vector.ml +++ b/linear_algebra/lib/vector.ml @@ -65,5 +65,6 @@ let normalize v = result -let at t i = t.{i} +let (%.) t i = t.{i} +let set t i v = t.{i} <- v diff --git a/linear_algebra/lib/vector.mli b/linear_algebra/lib/vector.mli index fe2865d..3b6ec62 100644 --- a/linear_algebra/lib/vector.mli +++ b/linear_algebra/lib/vector.mli @@ -76,9 +76,6 @@ val init : int -> (int -> float) -> 'a t val sum : 'a t -> float (** Returns the sum of the elements of the vector *) -val at : 'a t -> int -> float -(** Returns t.{i} *) - val copy : ?n:int -> ?ofsy:int -> ?incy:int -> ?y:vec -> ?ofsx:int -> ?incx:int -> 'a t -> 'a t (** Returns a copy of the vector X into Y. [ofs] controls the offset and [inc] the increment. *) @@ -134,3 +131,11 @@ val of_bigarray_inplace : (float, Stdlib.Bigarray.float64_elt, Stdlib.Bigarray.f val to_bigarray_inplace : 'a t -> (float, Stdlib.Bigarray.float64_elt, Stdlib.Bigarray.fortran_layout) Stdlib.Bigarray.Array1.t (** Converts the vector into a Fortran bigarray *) + + +val (%.) : 'a t -> int -> float +(** [t%.(i)] Returns the i-th element of the vector *) + +val set : 'a t -> int -> float -> unit +(** Modifies the value in-place at the i-th position *) + diff --git a/mo/lib/basis.ml b/mo/lib/basis.ml new file mode 100644 index 0000000..d9b21aa --- /dev/null +++ b/mo/lib/basis.ml @@ -0,0 +1,201 @@ +open Linear_algebra +open Common.Util +open Common.Constants + +(** One-electron orthogonal basis set, corresponding to Molecular Orbitals. *) + +module HF = HartreeFock +module Si = Simulation + +type mo_type = + | RHF | ROHF | UHF | CASSCF | Projected + | Natural of string + | Localized of string + +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 *) + (* TODO + f12_ints : F12.t lazy_t; (* F12 integrals *) + *) + } + + +let size t = + Matrix.dim2 t.mo_coef + +let simulation t = t.simulation +let mo_type t = t.mo_type +let ao_basis t = Si.ao_basis t.simulation +let mo_occupation t = t.mo_occupation +let mo_coef t = t.mo_coef +let eN_ints t = Lazy.force t.eN_ints +let ee_ints t = Lazy.force t.ee_ints +let kin_ints t = Lazy.force t.kin_ints +let two_e_ints t = Lazy.force t.ee_ints +(* TODO +let f12_ints t = Lazy.force t.f12_ints +*) +let one_e_ints t = Lazy.force t.one_e_ints + + +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 + match t.mo_type with + | RHF -> Fock.make_rhf ~density:m_P (ao_basis t) + | Projected + | ROHF -> (Matrix.scal 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.diag + + +let mo_matrix_of_ao_matrix ~mo_coef ao_matrix = + Matrix.xt_o_x ~x:mo_coef ~o:ao_matrix + + +let ao_matrix_of_mo_matrix ~mo_coef ~ao_overlap mo_matrix = + let sc = Matrix.gemm ao_overlap mo_coef in + Matrix.x_o_xt ~x:sc ~o:mo_matrix + + +let make ~simulation ~mo_type ~mo_occupation ~mo_coef () = + let ao_basis = + Si.ao_basis simulation + in + let eN_ints = lazy ( + Ao.Basis.eN_ints ao_basis + |> mo_matrix_of_ao_matrix ~mo_coef + ) + and kin_ints = lazy ( + Ao.Basis.kin_ints ao_basis + |> mo_matrix_of_ao_matrix ~mo_coef + ) + and ee_ints = lazy ( + Ao.Basis.ee_ints ao_basis + |> Eri.four_index_transform mo_coef + ) +(* + and f12_ints = lazy ( + Ao.Basis.f12_ints ao_basis + |> F12.four_index_transform mo_coef + ) +*) + in + let one_e_ints = lazy ( + Matrix.add (Lazy.force eN_ints) (Lazy.force kin_ints) ) + in + { simulation ; mo_type ; mo_occupation ; mo_coef ; + eN_ints ; ee_ints ; kin_ints ; one_e_ints ; + } + + +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 + +let of_hartree_fock hf = + let mo_coef = HF.eigenvectors hf in + let simulation = HF.simulation hf in + let mo_occupation = HF.occupation hf in + let mo_type = + match HF.kind hf with + | HartreeFock.RHF -> RHF + | HartreeFock.ROHF -> ROHF + | HartreeFock.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) ) + in + let m_X = AOBasis.ortho basis in + (* Project other vectors in the current basis *) + let m_C = + gemm m_S @@ mo_coef other + in + (* Append dummy vectors to the input vectors *) + let result = + let vecs = Mat.to_col_vecs m_X in + Array.iteri (fun i v -> if (i < Array.length vecs) then vecs.(i) <- v) + (Mat.to_col_vecs m_C) ; + Mat.of_col_vecs vecs + in + (* Gram-Schmidt Orthonormalization *) + gemm m_X @@ (Util.qr_ortho @@ gemm ~transa:`T m_X result) + |> Util.remove_epsilons + |> Conventions.rephase + in + + let mo_occupation = + let occ = mo_occupation other in + Vec.init (Mat.dim2 mo_coef) (fun i -> + if (i <= Vec.dim occ) then occ.{i} + 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 + in + let finish = + match finish with + | 0 -> cols + | x -> x + in + + let rec aux first = + + if (first > finish) then () + else + begin + Format.fprintf ppf "@[@[@[%s@;" "Eigenvalues:"; + + 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); + + Format.fprintf ppf "@]@;"; + Format.fprintf ppf "@[%a@]" + (Lacaml.Io.pp_lfmat + ~row_labels: + (Array.init rows (fun i -> Printf.sprintf "%d " (i + 1))) + ~col_labels: + (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)) ; + Format.fprintf ppf "@]@;@;@]"; + (aux [@tailcall]) (first+5) + end + in + aux start diff --git a/mo/lib/basis.mli b/mo/lib/basis.mli new file mode 100644 index 0000000..b243765 --- /dev/null +++ b/mo/lib/basis.mli @@ -0,0 +1,97 @@ +(** Data structure to represent the molecular orbitals. + + The MO indices start from 1. + +*) + +open Linear_algebra +open Common + +type mo_type = + | RHF | ROHF | UHF | CASSCF | Projected + | Natural of string + | Localized of string + +type t + +(** {1 Accessors} *) + +val simulation : t -> Simulation.t +(** Simulation which produced the MOs *) + +val mo_type : t -> mo_type +(** Kind of MOs (RHF, CASSCF, Localized...) *) + +val ao_basis : t -> Ao.Basis.t +(** Matrix of the MO coefficients in the AO basis *) + +val mo_occupation : t -> t Vector.t +(** Occupation numbers *) + +val mo_coef : t -> (Ao.Basis.t, t) Matrix.t +(** Molecular orbitcal coefficients *) + +val eN_ints : t -> (t,t) Matrix.t +(** Electron-nucleus potential integrals *) + +val ee_ints : t -> t Four_idx_storage.t +(** Electron-electron repulsion integrals *) + +val kin_ints : t -> (t,t) Matrix.t +(** Kinetic energy integrals *) + +val one_e_ints : t -> (t,t) Matrix.t +(** One-electron integrals {% $\hat{T} + V$ %} *) + +val two_e_ints : t -> t Four_idx_storage.t +(** Electron-electron repulsion integrals *) + +(* TODO +val f12_ints : t -> F12.t +(** F12 integrals *) +*) + +val size : t -> int +(** Number of molecular orbitals in the basis *) + +val mo_energies : t -> t Vector.t +(** Fock MO energies *) + +val values : t -> Coordinate.t -> t 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 -> + unit -> t +(** Function to build a data structure representing the molecular orbitals. *) + +val of_hartree_fock : HartreeFock.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 +(** 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 +(** Build a matrix in AO basis from a matrix in MO basis. *) + +(** {1 Printers} *) + +val pp : ?start:int -> ?finish:int -> Format.formatter -> t -> unit + + + diff --git a/mo/lib/class.ml b/mo/lib/class.ml new file mode 100644 index 0000000..265fe4b --- /dev/null +++ b/mo/lib/class.ml @@ -0,0 +1,141 @@ +open Particles + +type mo_class = + | Core of int (* Always doubly occupied *) + | Inactive of int (* With 0,1 or 2 holes *) + | Active of int (* With 0,1 or 2 holes or particles *) + | Virtual of int (* With 0,1 or 2 particles *) + | Deleted of int (* Always unoccupied *) + | Auxiliary of int (* Auxiliary basis function *) + +type t = mo_class list + + +let pp_mo_class ppf = function + | Core i -> Format.fprintf ppf "@[Core %d@]" i + | Inactive i -> Format.fprintf ppf "@[Inactive %d@]" i + | Active i -> Format.fprintf ppf "@[Active %d@]" i + | Virtual i -> Format.fprintf ppf "@[Virtual %d@]" i + | Deleted i -> Format.fprintf ppf "@[Deleted %d@]" i + | Auxiliary i -> Format.fprintf ppf "@[Auxiliary %d@]" i + +let pp ppf t = + Format.fprintf ppf "@[[@,"; + let rec aux = function + | [] -> Format.fprintf ppf "]@]" + | x :: [] -> Format.fprintf ppf "%a@,]@]" pp_mo_class x + | x :: rest -> ( Format.fprintf ppf "%a@,;@," pp_mo_class x; aux rest ) + in + aux t + + +let of_list t = t + +let to_list t = t + + +let core_mos t = + List.filter_map (fun x -> + match x with + | Core i -> Some i + | _ -> None) t + + +let inactive_mos t = + List.filter_map (fun x -> + match x with + | Inactive i -> Some i + | _ -> None ) t + + +let active_mos t = + List.filter_map (fun x -> + match x with + | Active i -> Some i + | _ -> None ) t + + +let virtual_mos t = + List.filter_map (fun x -> + match x with + | Virtual i -> Some i + | _ -> None ) t + + +let deleted_mos t = + List.filter_map (fun x -> + match x with + | Deleted i -> Some i + | _ -> None ) t + + + +let auxiliary_mos t = + List.filter_map (fun x -> + match x with + | Auxiliary i -> Some i + | _ -> None ) t + + +let mo_class_array t = + let sze = List.length t + 1 in + let result = Array.make sze (Deleted 0) in + List.iter (fun c -> + match c with + | Core i -> result.(i) <- Core i + | Inactive i -> result.(i) <- Inactive i + | Active i -> result.(i) <- Active i + | Virtual i -> result.(i) <- Virtual i + | Deleted i -> result.(i) <- Deleted i + | Auxiliary i -> result.(i) <- Auxiliary i + ) t; + result + + +let fci ~frozen_core mo_basis = + let mo_num = Basis.size mo_basis in + let ncore = (Nuclei.small_core @@ Simulation.nuclei @@ Basis.simulation mo_basis) / 2 in + of_list ( + if frozen_core then + List.concat [ + Util.list_range 1 ncore + |> List.map (fun i -> Core i) ; + Util.list_range (ncore+1) mo_num + |> List.map (fun i -> Active i) + ] + else + Util.list_range 1 mo_num + |> List.map (fun i -> Active i) + ) + +let cas_sd mo_basis ~frozen_core n m = + let mo_num = Basis.size mo_basis in + let n_alfa = Basis.simulation mo_basis |> Simulation.electrons |> Electrons.n_alfa in + let n_beta = Basis.simulation mo_basis |> Simulation.electrons |> Electrons.n_beta in + let n_unpaired = n_alfa - n_beta in + let n_alfa_in_cas = (n - n_unpaired)/2 + n_unpaired in + let last_inactive = n_alfa - n_alfa_in_cas in + let last_active = last_inactive + m in + let ncore = + if frozen_core then + (Nuclei.small_core @@ Simulation.nuclei @@ Basis.simulation mo_basis) / 2 + |> min last_inactive + else 0 + in + of_list ( + List.concat [ + if ncore > 0 then + Util.list_range 1 ncore + |> List.map (fun i -> Core i) + else + [] ; + Util.list_range (ncore+1) last_inactive + |> List.map (fun i -> Inactive i) ; + Util.list_range (last_inactive+1) last_active + |> List.map (fun i -> Active i) ; + Util.list_range (last_active+1) mo_num + |> List.map (fun i -> Virtual i) + ] + ) + + diff --git a/mo/lib/class.mli b/mo/lib/class.mli new file mode 100644 index 0000000..0f1be50 --- /dev/null +++ b/mo/lib/class.mli @@ -0,0 +1,58 @@ +(** Classes of MOs : active, inactive, etc *) + +type mo_class = + | Core of int (* Always doubly occupied *) + | Inactive of int (* With 0,1 or 2 holes *) + | Active of int (* With 0,1 or 2 holes or particles *) + | Virtual of int (* With 0,1 or 2 particles *) + | Deleted of int (* Always unoccupied *) + | Auxiliary of int (* Function of the auxiliary basis set *) + +type t + +(** Creation *) +val of_list : mo_class list -> t + +val to_list : t -> mo_class list + +val fci : frozen_core:bool -> Basis.t -> t +(** Creates the MO classes for FCI calculations : all [Active]. The + [n] lowest MOs are [Core] if [frozen_core = true]. +*) + +val cas_sd: Basis.t -> frozen_core:bool -> int -> int -> t +(** [cas_sd mo_basis n m ] creates the MO classes for CAS(n,m) + SD + calculations. lowest MOs are [Core], then all the next MOs are [Inactive], + then [Active], then [Virtual]. +*) + + +val core_mos : t -> int list +(** Returns a list containing the indices of the core MOs. *) + +val active_mos : t -> int list +(** Returns a list containing the indices of the active MOs. *) + +val virtual_mos : t -> int list +(** Returns a list containing the indices of the virtual MOs. *) + +val inactive_mos : t -> int list +(** Returns a list containing the indices of the inactive MOs. *) + +val deleted_mos : t -> int list +(** Returns a list containing the indices of the deleted MOs. *) + +val auxiliary_mos : t -> int list +(** Returns a list containing the indices of the auxiliary MOs. *) + +val mo_class_array : t -> mo_class array +(** Returns an array [a] such that [a.(i)] returns the class of MO [i]. + As the MO indices start from [1], the array has an extra zero entry + that should be ignored. *) + +(** {2 Printers} *) + +val pp_mo_class : Format.formatter -> mo_class -> unit + +val pp : Format.formatter -> t -> unit + diff --git a/mo/lib/dune b/mo/lib/dune new file mode 100644 index 0000000..f53c2bf --- /dev/null +++ b/mo/lib/dune @@ -0,0 +1,9 @@ +; name = name of the supermodule that will wrap all source files as submodules +; public_name = name of the library for ocamlfind and opam +(library + (name mo) + (public_name qcaml.mo) + (libraries + qcaml.simulation + ) + (synopsis "Molecular orbitals.")) diff --git a/simulation/lib/dune b/simulation/lib/dune new file mode 100644 index 0000000..994cfd9 --- /dev/null +++ b/simulation/lib/dune @@ -0,0 +1,14 @@ +; name = name of the supermodule that will wrap all source files as submodules +; public_name = name of the library for ocamlfind and opam +(library + (name simulation) + (public_name qcaml.simulation) + (libraries + qcaml.common + qcaml.particles + qcaml.gaussian_basis + qcaml.gaussian_integrals + qcaml.operators + qcaml.ao + ) + (synopsis "Contains data describing a simulation (AOs, operators, nuclear coordinate,...)")) diff --git a/qcaml/lib/simulation.ml b/simulation/lib/simulation.ml similarity index 100% rename from qcaml/lib/simulation.ml rename to simulation/lib/simulation.ml diff --git a/qcaml/lib/simulation.mli b/simulation/lib/simulation.mli similarity index 100% rename from qcaml/lib/simulation.mli rename to simulation/lib/simulation.mli