From 0b6fe42e55da8ac0d55c855fd824a2efe767b055 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sat, 30 Jan 2021 19:07:59 +0100 Subject: [PATCH] Localization --- docs/linear_algebra.html | 10 +- docs/mo.html | 194 ++++++++++++-- docs/top.html | 10 +- linear_algebra/lib/matrix.ml | 12 +- linear_algebra/lib/matrix.mli | 4 +- mo/lib/hartree_fock.ml | 4 +- mo/lib/localization.ml | 331 +++++++++++++++++++++++ mo/lib/localization.mli | 54 ++++ mo/localization.org | 477 ++++++++++++++++++++++++++++++++++ top/lib/install_printers.ml | 1 + 10 files changed, 1056 insertions(+), 41 deletions(-) create mode 100644 mo/lib/localization.ml create mode 100644 mo/lib/localization.mli create mode 100644 mo/localization.org diff --git a/docs/linear_algebra.html b/docs/linear_algebra.html index 1258139..20ba805 100644 --- a/docs/linear_algebra.html +++ b/docs/linear_algebra.html @@ -3,7 +3,7 @@ "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd"> - + Linear Algebra @@ -250,18 +250,18 @@ org_html_manager.setup(); // activate after the parameters are set

Table of Contents

-
-

1 Summmary

+
+

1 Summmary

Author: Anthony Scemama

-

Created: 2021-01-27 Wed 23:44

+

Created: 2021-01-29 Fri 15:35

Validate

diff --git a/docs/mo.html b/docs/mo.html index e82119f..a4de8c7 100644 --- a/docs/mo.html +++ b/docs/mo.html @@ -3,7 +3,7 @@ "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd"> - + Molecular orbitals @@ -238,6 +238,28 @@ org_html_manager.setup(); // activate after the parameters are set /*]]>*///--> // @license-end + +
@@ -250,36 +272,46 @@ org_html_manager.setup(); // activate after the parameters are set

Table of Contents

-
-

1 Summmary

+
+

1 Summmary

-
-

2 Frozen core

+
+

2 Frozen core

Defines how the core electrons are frozen, for each atom.

-
-

2.1 Type

+
+

2.1 Type

-
type kind =
+
type kind =
   | All_electron
   | Small
   | Large
@@ -293,8 +325,8 @@ Defines how the core electrons are frozen, for each atom.
-
-

2.2 Creation

+
+

2.2 Creation

val make : kind -> Particles.Nuclei.t -> t
@@ -330,7 +362,7 @@ Defines how the core electrons are frozen, for each atom.
 
 
 
-
+
 let f = Frozen_core.(make Small nuclei) ;;
 val f : Frozen_core.t = [|0; 2; 2; 0|]
 
@@ -340,8 +372,8 @@ val f : Frozen_core.t = [|0; 2; 2; 0|]
-
-

2.3 Access

+
+

2.3 Access

val num_elec : t -> int
@@ -370,7 +402,7 @@ val f : Frozen_core.t = [|0; 2; 2; 0|]
 
 
 
-
+
 Frozen_core.num_elec f ;;
 - : int = 4
 
@@ -380,8 +412,8 @@ Frozen_core.num_mos f ;;
-
-

2.4 Printers

+
+

2.4 Printers

val pp : Format.formatter -> t -> unit
@@ -390,10 +422,126 @@ Frozen_core.num_mos f ;;
+ +
+

3 Orbital localization

+
+

+Molecular orbital localization function. +

+ +

+Boys: +

+ +

+Edmiston-Rudenberg: +

+
+ + +
+

3.1 Type

+
+
+
open Linear_algebra
+
+type localization_kind =
+  | Edmiston
+  | Boys
+
+type mo = Mo_dim.t
+type ao = Ao.Ao_dim.t
+type loc
+
+
+ +
+
type localization_data 
+type t 
+
+
+
+
+ +
+

3.2 Edmiston-Rudenberg

+
+ +
+

3.3 Boys

+
+ +
+

3.4 Access

+
+
+
val kind         : t -> localization_kind
+val simulation   : t -> Simulation.t
+val selected_mos : t -> int list
+
+val kappa :
+  kind:localization_kind ->
+  Basis.t ->
+  ( ao,loc) Matrix.t ->
+  (loc,loc) Matrix.t * float
+
+val make  :
+  kind:localization_kind ->
+  ?max_iter:int -> 
+  ?convergence:float -> 
+  Basis.t ->
+  int list ->
+  t
+
+val to_basis : t -> Basis.t
+
+
+
+ + + + +++ ++ + + + + + + + + + + + +
kappaReturns the \(\kappa\) antisymmetric matrix used for the rotation matrix and the value of the localization function
makePerforms the orbital localization
+
+
+ +
+

3.5 Printers

+
+
+
(*
+  val pp : Format.formatter -> t -> unit
+  *)
+
+
+
+
+ +
+

3.6 Tests

+
+

Author: Anthony Scemama

-

Created: 2021-01-27 Wed 23:44

+

Created: 2021-01-30 Sat 19:07

Validate

diff --git a/docs/top.html b/docs/top.html index b9aa28d..5547b12 100644 --- a/docs/top.html +++ b/docs/top.html @@ -3,7 +3,7 @@ "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd"> - + Top-level @@ -250,18 +250,18 @@ org_html_manager.setup(); // activate after the parameters are set

Table of Contents

-
-

1 Summmary

+
+

1 Summmary

Author: Anthony Scemama

-

Created: 2021-01-28 Thu 00:32

+

Created: 2021-01-30 Sat 19:07

Validate

diff --git a/linear_algebra/lib/matrix.ml b/linear_algebra/lib/matrix.ml index 568e968..b1c58c3 100644 --- a/linear_algebra/lib/matrix.ml +++ b/linear_algebra/lib/matrix.ml @@ -393,7 +393,7 @@ let qr a = q, r -let exponential_iterative a = +let exponential a = assert (dim1 a = dim2 a); let rec loop result accu n = let b = scale (1./.n) a in @@ -413,7 +413,7 @@ let exponential_iterative a = loop id id 1. -let exponential a = +let exponential_antisymmetric a = let n = dim1 a in assert (n = dim2 a); @@ -433,8 +433,12 @@ let exponential a = |> amax |> abs_float in - assert (test < Constants.epsilon); - result + + (* If the exponential failed, fall back to the iterative exponential *) + if test < 10. *. Constants.epsilon then + result + else + exponential a let to_file ~filename ?(sym=false) ?(cutoff=0.) t = diff --git a/linear_algebra/lib/matrix.mli b/linear_algebra/lib/matrix.mli index f9e4b44..079e245 100644 --- a/linear_algebra/lib/matrix.mli +++ b/linear_algebra/lib/matrix.mli @@ -306,8 +306,8 @@ val diagonalize_symm : ('a,'a) t -> ('a,'a) t * 'a Vector.t val exponential : ('a,'a) t -> ('a,'a) t (** Computes the exponential of a square matrix. *) -val exponential_iterative : ('a,'a) t -> ('a,'a) t -(** Computes the exponential of a square matrix with an iteratve algorithm. *) +val exponential_antisymmetric: ('a,'a) t -> ('a,'a) t +(** Computes the exponential of an antisymmetric square matrix. *) val xt_o_x : o:('a,'a) t -> x:('a,'b) t -> ('b,'b) t (** Computes {% $\mathbf{X^\dag\, O\, X}$ %} *) diff --git a/mo/lib/hartree_fock.ml b/mo/lib/hartree_fock.ml index 8e307a2..3ac4b5d 100644 --- a/mo/lib/hartree_fock.ml +++ b/mo/lib/hartree_fock.ml @@ -7,7 +7,7 @@ type mo = Mo_dim.t type hartree_fock_data = { - iteration : int ; + iteration : int ; coefficients : (ao, mo) Matrix.t option ; eigenvalues : mo Vector.t option ; error : float option ; @@ -77,7 +77,7 @@ let n_iterations t = let last_iteration t = - Util.of_some @@ Lazy.force (t.data.(n_iterations t - 1)) + Util.of_some @@ Lazy.force t.data.(n_iterations t - 1) let eigenvectors t = let data = last_iteration t in diff --git a/mo/lib/localization.ml b/mo/lib/localization.ml new file mode 100644 index 0000000..c87763d --- /dev/null +++ b/mo/lib/localization.ml @@ -0,0 +1,331 @@ +(* [[file:~/QCaml/mo/localization.org::*Type][Type:3]] *) +open Linear_algebra + +type localization_kind = + | Edmiston + | Boys + +type mo = Mo_dim.t +type ao = Ao.Ao_dim.t +type loc + +type localization_data = + { + coefficients : (ao, loc) Matrix.t ; + kappa : (loc, loc) Matrix.t ; + scaling : float ; + loc_value : float ; + iteration : int ; + } + +type t = + { + kind : localization_kind ; + mo_basis : Basis.t ; + data : localization_data option lazy_t array ; + selected_mos : int list ; + } + +open Common +(* Type:3 ends here *) + +(* Edmiston-Rudenberg *) + + +(* [[file:~/QCaml/mo/localization.org::*Edmiston-Rudenberg][Edmiston-Rudenberg:1]] *) +let kappa_edmiston in_basis m_C = + let ao_basis = + Basis.simulation in_basis + |> Simulation.ao_basis + in + + let ee_ints = Ao.Basis.ee_ints ao_basis in + let n_ao = Ao.Basis.size ao_basis in + + let n_mo = Matrix.dim2 m_C in + + (* Temp arrays for integral transformation *) + let m_pqr = + Bigarray.(Array3.create Float64 fortran_layout n_ao n_ao n_ao) + in + let m_qr_i = Matrix.create (n_ao*n_ao) n_mo in + let m_ri_j = Matrix.create (n_ao*n_mo) n_mo in + let m_ij_k = Matrix.create (n_mo*n_mo) n_mo in + + + let m_a12 = Bigarray.(Array2.create Float64 fortran_layout n_mo n_mo) in + let m_b12 = Bigarray.(Array2.create Float64 fortran_layout n_mo n_mo) in + Bigarray.Array2.fill m_b12 0.; + Bigarray.Array2.fill m_a12 0.; + let v_d = + Vector.init n_mo (fun _ -> 0.) + |> Vector.to_bigarray_inplace + in + + Array.iter (fun s -> + + Array.iter (fun r -> + Array.iter (fun q -> + Array.iter (fun p -> + m_pqr.{p,q,r} <- Four_idx_storage.get_phys ee_ints p q r s + ) (Util.array_range 1 n_ao) + ) (Util.array_range 1 n_ao) + ) (Util.array_range 1 n_ao); + + let m_p_qr = + Bigarray.reshape (Bigarray.genarray_of_array3 m_pqr) [| n_ao ; n_ao*n_ao |] + |> Bigarray.array2_of_genarray + |> Matrix.of_bigarray_inplace + in + + (* (qr,i) = = \sum_p

C_{pi} *) + Matrix.gemm_tn_inplace ~c:m_qr_i m_p_qr m_C; + + let m_q_ri = + let x = Matrix.to_bigarray_inplace m_qr_i |> Bigarray.genarray_of_array2 in + Bigarray.reshape_2 x n_ao (n_ao*n_mo) |> Matrix.of_bigarray_inplace + in + + (* (ri,j) = = \sum_q C_{qj} *) + Matrix.gemm_tn_inplace ~c:m_ri_j m_q_ri m_C; + + let m_r_ij = + let x = Matrix.to_bigarray_inplace m_ri_j |> Bigarray.genarray_of_array2 in + Bigarray.reshape_2 x n_ao (n_mo*n_mo) |> Matrix.of_bigarray_inplace + in + + (* (ij,k) = = \sum_r C_{rk} *) + Matrix.gemm_tn_inplace ~c:m_ij_k m_r_ij m_C; + + let m_ijk = + let x = Matrix.to_bigarray_inplace m_ij_k |> Bigarray.genarray_of_array2 in + Bigarray.reshape x [| n_mo ; n_mo ; n_mo |] + |> Bigarray.array3_of_genarray + in + + let m_Cx = Matrix.to_bigarray_inplace m_C in + Array.iter (fun j -> + Array.iter (fun i -> + m_b12.{i,j} <- m_b12.{i,j} +. m_Cx.{s,j} *. (m_ijk.{i,i,i} -. m_ijk.{j,i,j}); + m_a12.{i,j} <- m_a12.{i,j} +. m_ijk.{i,i,j} *. m_Cx.{s,j} -. + 0.25 *. ( (m_ijk.{i,i,i} -. m_ijk.{j,i,j}) *. m_Cx.{s,i} +. + (m_ijk.{j,j,j} -. m_ijk.{i,j,i}) *. m_Cx.{s,j}) + ) (Util.array_range 1 n_mo); + v_d.{j} <- v_d.{j} +. m_ijk.{j,j,j} *. m_Cx.{s,j} + ) (Util.array_range 1 n_mo) + + ) (Util.array_range 1 n_ao); + + let f i j = + if i=j then + 0. + else + begin + let x = 1./. sqrt (m_b12.{i,j} *. m_b12.{i,j} +. m_a12.{i,j} *. m_a12.{i,j} ) in + if asin (m_b12.{i,j} /. x) > 0. then + 0.25 *. acos( -. m_a12.{i,j} *. x) + else + -. 0.25 *. acos( -. m_a12.{i,j} *. x ) + end + in + ( + Matrix.init_cols n_mo n_mo ( fun i j -> if i<=j then f i j else -. (f j i) ), + Vector.sum (Vector.of_bigarray_inplace v_d) + ) +(* Edmiston-Rudenberg:1 ends here *) + +(* Boys *) + + +(* [[file:~/QCaml/mo/localization.org::*Boys][Boys:1]] *) +let kappa_boys in_basis = + let ao_basis = + Basis.simulation in_basis + |> Simulation.ao_basis + in + let multipole = Ao.Basis.multipole ao_basis in + fun m_C -> + let n_mo = Matrix.dim2 m_C in + + let phi_x_phi = Matrix.xt_o_x ~x:m_C ~o:(multipole "x") in + let phi_y_phi = Matrix.xt_o_x ~x:m_C ~o:(multipole "y") in + let phi_z_phi = Matrix.xt_o_x ~x:m_C ~o:(multipole "z") in + + let m_b12 = + let g x i j = + let x_ii = x%:(i,i) in + let x_jj = x%:(j,j) in + let x_ij = x%:(i,j) in + (x_ii -. x_jj) *. x_ij + in + Matrix.init_cols n_mo n_mo (fun i j -> + let x = + (g phi_x_phi i j) +. (g phi_y_phi i j) +. (g phi_z_phi i j) + in + if (abs_float x > 1.e-15) then x else 0. + ) + in + + let m_a12 = + let g x i j = + let x_ii = x%:(i,i) in + let x_jj = x%:(j,j) in + let x_ij = x%:(i,j) in + let y = x_ii -. x_jj in + (x_ij *. x_ij) -. 0.25 *. y *. y + in + Matrix.init_cols n_mo n_mo (fun i j -> + let x = + (g phi_x_phi i j) +. (g phi_y_phi i j) +. (g phi_z_phi i j) + in + if (abs_float x > 1.e-15) then x else 0. + ) + in + + let f i j = + let pi = Constants.pi in + if i=j then + 0. + else + let x = atan2 (m_b12%:(i,j)) (m_a12%:(i,j)) in + if x >= 0. then + 0.25 *. (pi -. x) + else + -. 0.25 *. ( pi +. x) + in + ( + Matrix.init_cols n_mo n_mo ( fun i j -> if i>j then f i j else -. (f j i) ), + let r2 x y z = x*.x +. y*.y +. z*.z in + Vector.init n_mo ( fun i -> + r2 (phi_x_phi%:(i,i)) (phi_y_phi%:(i,i)) (phi_z_phi%:(i,i))) + |> Vector.sum + ) +(* Boys:1 ends here *) + + + +(* | ~kappa~ | Returns the $\kappa$ antisymmetric matrix used for the rotation matrix and the value of the localization function | + * | ~make~ | Performs the orbital localization | *) + + +(* [[file:~/QCaml/mo/localization.org::*Access][Access:2]] *) +let kind t = t.kind +let simulation t = Basis.simulation t.mo_basis +let selected_mos t = t.selected_mos + +let kappa ~kind = + match kind with + | Edmiston -> kappa_edmiston + | Boys -> kappa_boys + + +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 ao_basis t = Simulation.ao_basis (simulation t) +*) + + + +let make ~kind ?(max_iter=500) ?(convergence=1.e-8) mo_basis selected_mos = + + let kappa_loc = kappa ~kind mo_basis in + + let mo_array = Matrix.to_col_vecs (Basis.mo_coef mo_basis) in + let mos_vec_list = List.map (fun i -> mo_array.(i-1)) selected_mos in + let x: (ao,loc) Matrix.t = + Matrix.of_col_vecs_list mos_vec_list + in + + let next_coef kappa = + let r = Matrix.exponential_antisymmetric kappa in + let m_C = Matrix.gemm_nt x r in + m_C + in + + let data_initial = + let iteration = 0 + and kappa, loc_value = kappa_loc x + and scaling = 1.0 + in + let coefficients = next_coef kappa in + { coefficients ; kappa ; scaling ; loc_value ; iteration } + in + + let iteration data = + let x = data.coefficients in + let iteration = data.iteration + 1 + and kappa, loc_value = kappa_loc x + and scaling = data.scaling + in + let coefficients = next_coef kappa in + Printf.printf "%4d %20f\n" iteration loc_value ; + { coefficients ; kappa ; scaling ; loc_value ; iteration } + in + + let array_data = + + let storage = + Array.make max_iter None + in + + let rec loop = function + | 0 -> Some (iteration data_initial) + | n -> begin + match storage.(n) with + | Some result -> Some result + | None -> begin + let data = loop (n-1) in + match data with + | None -> None + | Some data -> begin + (* Check convergence *) + let converged = + abs_float data.loc_value < convergence + in + if converged then + None + else + begin + storage.(n-1) <- Some data ; + storage.(n) <- Some (iteration data); + storage.(n) + end + end + end + end + in + Array.init max_iter (fun i -> lazy (loop i)) + in + { kind ; mo_basis ; data = array_data ; selected_mos } + + + +let to_basis t = + let mo_basis = t.mo_basis in + let simulation = Basis.simulation mo_basis in + let mo_occupation = Basis.mo_occupation mo_basis in + + let data = last_iteration t in + + let mo_coef_array = Matrix.to_col_vecs (Basis.mo_coef mo_basis) in + let new_mos = + Matrix.to_col_vecs data.coefficients + in + selected_mos t + |> List.iteri (fun i j -> mo_coef_array.(j) <- new_mos.(i)) ; + let mo_coef = Matrix.of_col_vecs mo_coef_array in + Basis.make ~simulation ~mo_type:(Localized "Boys") ~mo_occupation ~mo_coef () +(* Access:2 ends here *) + +(* [[file:~/QCaml/mo/localization.org::*Printers][Printers:2]] *) + +(* Printers:2 ends here *) diff --git a/mo/lib/localization.mli b/mo/lib/localization.mli new file mode 100644 index 0000000..db52f29 --- /dev/null +++ b/mo/lib/localization.mli @@ -0,0 +1,54 @@ +(* Type + * + * #+NAME: types *) + +(* [[file:~/QCaml/mo/localization.org::types][types]] *) +open Linear_algebra + +type localization_kind = + | Edmiston + | Boys + +type mo = Mo_dim.t +type ao = Ao.Ao_dim.t +type loc +(* types ends here *) + +(* [[file:~/QCaml/mo/localization.org::*Type][Type:2]] *) +type localization_data +type t +(* Type:2 ends here *) + +(* Access *) + + +(* [[file:~/QCaml/mo/localization.org::*Access][Access:1]] *) +val kind : t -> localization_kind +val simulation : t -> Simulation.t +val selected_mos : t -> int list + +val kappa : + kind:localization_kind -> + Basis.t -> + ( ao,loc) Matrix.t -> + (loc,loc) Matrix.t * float + +val make : + kind:localization_kind -> + ?max_iter:int -> + ?convergence:float -> + Basis.t -> + int list -> + t + +val to_basis : t -> Basis.t +(* Access:1 ends here *) + +(* Printers *) + + +(* [[file:~/QCaml/mo/localization.org::*Printers][Printers:1]] *) +(* + val pp : Format.formatter -> t -> unit + *) +(* Printers:1 ends here *) diff --git a/mo/localization.org b/mo/localization.org new file mode 100644 index 0000000..fa334b4 --- /dev/null +++ b/mo/localization.org @@ -0,0 +1,477 @@ +#+begin_src elisp tangle: no :results none :exports none +(setq pwd (file-name-directory buffer-file-name)) +(setq name (file-name-nondirectory (substring buffer-file-name 0 -4))) +(setq lib (concat pwd "lib/")) +(setq testdir (concat pwd "test/")) +(setq mli (concat lib name ".mli")) +(setq ml (concat lib name ".ml")) +(setq test-ml (concat testdir name ".ml")) +(org-babel-tangle) +#+end_src + +* Orbital localization + :PROPERTIES: + :header-args: :noweb yes :comments both + :END: + + Molecular orbital localization function. + + Boys: + + Edmiston-Rudenberg: + + +** Type + + #+NAME: types + #+begin_src ocaml :tangle (eval mli) +open Linear_algebra + +type localization_kind = + | Edmiston + | Boys + +type mo = Mo_dim.t +type ao = Ao.Ao_dim.t +type loc + #+end_src + + #+begin_src ocaml :tangle (eval mli) +type localization_data +type t + #+end_src + + #+begin_src ocaml :tangle (eval ml) :exports none +<> + +type localization_data = + { + coefficients : (ao, loc) Matrix.t ; + kappa : (loc, loc) Matrix.t ; + scaling : float ; + loc_value : float ; + iteration : int ; + } + +type t = + { + kind : localization_kind ; + mo_basis : Basis.t ; + data : localization_data option lazy_t array ; + selected_mos : int list ; + } + +open Common + #+end_src + +** Edmiston-Rudenberg + + #+begin_src ocaml :tangle (eval ml) :exports none +let kappa_edmiston in_basis m_C = + let ao_basis = + Basis.simulation in_basis + |> Simulation.ao_basis + in + + let ee_ints = Ao.Basis.ee_ints ao_basis in + let n_ao = Ao.Basis.size ao_basis in + + let n_mo = Matrix.dim2 m_C in + + (* Temp arrays for integral transformation *) + let m_pqr = + Bigarray.(Array3.create Float64 fortran_layout n_ao n_ao n_ao) + in + let m_qr_i = Matrix.create (n_ao*n_ao) n_mo in + let m_ri_j = Matrix.create (n_ao*n_mo) n_mo in + let m_ij_k = Matrix.create (n_mo*n_mo) n_mo in + + + let m_a12 = Bigarray.(Array2.create Float64 fortran_layout n_mo n_mo) in + let m_b12 = Bigarray.(Array2.create Float64 fortran_layout n_mo n_mo) in + Bigarray.Array2.fill m_b12 0.; + Bigarray.Array2.fill m_a12 0.; + let v_d = + Vector.init n_mo (fun _ -> 0.) + |> Vector.to_bigarray_inplace + in + + Array.iter (fun s -> + + Array.iter (fun r -> + Array.iter (fun q -> + Array.iter (fun p -> + m_pqr.{p,q,r} <- Four_idx_storage.get_phys ee_ints p q r s + ) (Util.array_range 1 n_ao) + ) (Util.array_range 1 n_ao) + ) (Util.array_range 1 n_ao); + + let m_p_qr = + Bigarray.reshape (Bigarray.genarray_of_array3 m_pqr) [| n_ao ; n_ao*n_ao |] + |> Bigarray.array2_of_genarray + |> Matrix.of_bigarray_inplace + in + + (* (qr,i) = = \sum_p

C_{pi} *) + Matrix.gemm_tn_inplace ~c:m_qr_i m_p_qr m_C; + + let m_q_ri = + let x = Matrix.to_bigarray_inplace m_qr_i |> Bigarray.genarray_of_array2 in + Bigarray.reshape_2 x n_ao (n_ao*n_mo) |> Matrix.of_bigarray_inplace + in + + (* (ri,j) = = \sum_q C_{qj} *) + Matrix.gemm_tn_inplace ~c:m_ri_j m_q_ri m_C; + + let m_r_ij = + let x = Matrix.to_bigarray_inplace m_ri_j |> Bigarray.genarray_of_array2 in + Bigarray.reshape_2 x n_ao (n_mo*n_mo) |> Matrix.of_bigarray_inplace + in + + (* (ij,k) = = \sum_r C_{rk} *) + Matrix.gemm_tn_inplace ~c:m_ij_k m_r_ij m_C; + + let m_ijk = + let x = Matrix.to_bigarray_inplace m_ij_k |> Bigarray.genarray_of_array2 in + Bigarray.reshape x [| n_mo ; n_mo ; n_mo |] + |> Bigarray.array3_of_genarray + in + + let m_Cx = Matrix.to_bigarray_inplace m_C in + Array.iter (fun j -> + Array.iter (fun i -> + m_b12.{i,j} <- m_b12.{i,j} +. m_Cx.{s,j} *. (m_ijk.{i,i,i} -. m_ijk.{j,i,j}); + m_a12.{i,j} <- m_a12.{i,j} +. m_ijk.{i,i,j} *. m_Cx.{s,j} -. + 0.25 *. ( (m_ijk.{i,i,i} -. m_ijk.{j,i,j}) *. m_Cx.{s,i} +. + (m_ijk.{j,j,j} -. m_ijk.{i,j,i}) *. m_Cx.{s,j}) + ) (Util.array_range 1 n_mo); + v_d.{j} <- v_d.{j} +. m_ijk.{j,j,j} *. m_Cx.{s,j} + ) (Util.array_range 1 n_mo) + + ) (Util.array_range 1 n_ao); + + let f i j = + if i=j then + 0. + else + begin + let x = 1./. sqrt (m_b12.{i,j} *. m_b12.{i,j} +. m_a12.{i,j} *. m_a12.{i,j} ) in + if asin (m_b12.{i,j} /. x) > 0. then + 0.25 *. acos( -. m_a12.{i,j} *. x) + else + -. 0.25 *. acos( -. m_a12.{i,j} *. x ) + end + in + ( + Matrix.init_cols n_mo n_mo ( fun i j -> if i<=j then f i j else -. (f j i) ), + Vector.sum (Vector.of_bigarray_inplace v_d) + ) + + + + #+end_src + +** Boys + + #+begin_src ocaml :tangle (eval ml) :exports none +let kappa_boys in_basis = + let ao_basis = + Basis.simulation in_basis + |> Simulation.ao_basis + in + let multipole = Ao.Basis.multipole ao_basis in + fun m_C -> + let n_mo = Matrix.dim2 m_C in + + let phi_x_phi = Matrix.xt_o_x ~x:m_C ~o:(multipole "x") in + let phi_y_phi = Matrix.xt_o_x ~x:m_C ~o:(multipole "y") in + let phi_z_phi = Matrix.xt_o_x ~x:m_C ~o:(multipole "z") in + + let m_b12 = + let g x i j = + let x_ii = x%:(i,i) in + let x_jj = x%:(j,j) in + let x_ij = x%:(i,j) in + (x_ii -. x_jj) *. x_ij + in + Matrix.init_cols n_mo n_mo (fun i j -> + let x = + (g phi_x_phi i j) +. (g phi_y_phi i j) +. (g phi_z_phi i j) + in + if (abs_float x > 1.e-15) then x else 0. + ) + in + + let m_a12 = + let g x i j = + let x_ii = x%:(i,i) in + let x_jj = x%:(j,j) in + let x_ij = x%:(i,j) in + let y = x_ii -. x_jj in + (x_ij *. x_ij) -. 0.25 *. y *. y + in + Matrix.init_cols n_mo n_mo (fun i j -> + let x = + (g phi_x_phi i j) +. (g phi_y_phi i j) +. (g phi_z_phi i j) + in + if (abs_float x > 1.e-15) then x else 0. + ) + in + + let f i j = + let pi = Constants.pi in + if i=j then + 0. + else + let x = atan2 (m_b12%:(i,j)) (m_a12%:(i,j)) in + if x >= 0. then + 0.25 *. (pi -. x) + else + -. 0.25 *. ( pi +. x) + in + ( + Matrix.init_cols n_mo n_mo ( fun i j -> if i>j then f i j else -. (f j i) ), + let r2 x y z = x*.x +. y*.y +. z*.z in + Vector.init n_mo ( fun i -> + r2 (phi_x_phi%:(i,i)) (phi_y_phi%:(i,i)) (phi_z_phi%:(i,i))) + |> Vector.sum + ) + + #+end_src + + +** Access + + #+begin_src ocaml :tangle (eval mli) +val kind : t -> localization_kind +val simulation : t -> Simulation.t +val selected_mos : t -> int list + +val kappa : + kind:localization_kind -> + Basis.t -> + ( ao,loc) Matrix.t -> + (loc,loc) Matrix.t * float + +val make : + kind:localization_kind -> + ?max_iter:int -> + ?convergence:float -> + Basis.t -> + int list -> + t + +val to_basis : t -> Basis.t + + #+end_src + + | ~kappa~ | Returns the $\kappa$ antisymmetric matrix used for the rotation matrix and the value of the localization function | + | ~make~ | Performs the orbital localization | + + #+begin_src ocaml :tangle (eval ml) :exports none +let kind t = t.kind +let simulation t = Basis.simulation t.mo_basis +let selected_mos t = t.selected_mos + +let kappa ~kind = + match kind with + | Edmiston -> kappa_edmiston + | Boys -> kappa_boys + + +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 ao_basis t = Simulation.ao_basis (simulation t) +*) + + + +let make ~kind ?(max_iter=500) ?(convergence=1.e-8) mo_basis selected_mos = + + let kappa_loc = kappa ~kind mo_basis in + + let mo_array = Matrix.to_col_vecs (Basis.mo_coef mo_basis) in + let mos_vec_list = List.map (fun i -> mo_array.(i-1)) selected_mos in + let x: (ao,loc) Matrix.t = + Matrix.of_col_vecs_list mos_vec_list + in + + let next_coef kappa = + let r = Matrix.exponential_antisymmetric kappa in + let m_C = Matrix.gemm_nt x r in + m_C + in + + let data_initial = + let iteration = 0 + and kappa, loc_value = kappa_loc x + and scaling = 1.0 + in + let coefficients = next_coef kappa in + { coefficients ; kappa ; scaling ; loc_value ; iteration } + in + + let iteration data = + let x = data.coefficients in + let iteration = data.iteration + 1 + and kappa, loc_value = kappa_loc x + and scaling = data.scaling + in + let coefficients = next_coef kappa in + Printf.printf "%4d %20f\n" iteration loc_value ; + { coefficients ; kappa ; scaling ; loc_value ; iteration } + in + + let array_data = + + let storage = + Array.make max_iter None + in + + let rec loop = function + | 0 -> Some (iteration data_initial) + | n -> begin + match storage.(n) with + | Some result -> Some result + | None -> begin + let data = loop (n-1) in + match data with + | None -> None + | Some data -> begin + (* Check convergence *) + let converged = + abs_float data.loc_value < convergence + in + if converged then + None + else + begin + storage.(n-1) <- Some data ; + storage.(n) <- Some (iteration data); + storage.(n) + end + end + end + end + in + Array.init max_iter (fun i -> lazy (loop i)) + in + { kind ; mo_basis ; data = array_data ; selected_mos } + + + +let to_basis t = + let mo_basis = t.mo_basis in + let simulation = Basis.simulation mo_basis in + let mo_occupation = Basis.mo_occupation mo_basis in + + let data = last_iteration t in + + let mo_coef_array = Matrix.to_col_vecs (Basis.mo_coef mo_basis) in + let new_mos = + Matrix.to_col_vecs data.coefficients + in + selected_mos t + |> List.iteri (fun i j -> mo_coef_array.(j) <- new_mos.(i)) ; + let mo_coef = Matrix.of_col_vecs mo_coef_array in + Basis.make ~simulation ~mo_type:(Localized "Boys") ~mo_occupation ~mo_coef () + + #+end_src + +** Printers + + #+begin_src ocaml :tangle (eval mli) +(* + val pp : Format.formatter -> t -> unit + *) + #+end_src + + #+begin_src ocaml :tangle (eval ml) :exports none + + #+end_src + +** Tests + + #+begin_src ocaml :tangle (eval test-ml) :exports none + +let test_localization = + + let nuclei = + Particles.Nuclei.of_xyz_string + " 10 +Hydrogen chain, d=1.8 Angstrom +H -4.286335 0.000000 0.000000 +H -3.333816 0.000000 0.000000 +H -2.381297 0.000000 0.000000 +H -1.428778 0.000000 0.000000 +H -0.476259 0.000000 0.000000 +H 0.476259 0.000000 0.000000 +H 1.428778 0.000000 0.000000 +H 2.381297 0.000000 0.000000 +H 3.333816 0.000000 0.000000 +H 4.286335 0.000000 0.000000 +" in + + let basis_file = "/home/scemama/qp2/data/basis/sto-6g" in + + let ao_basis = + Ao.Basis.of_nuclei_and_basis_filename ~nuclei basis_file + in + + + let charge = 0 in + + let multiplicity = 1 in + + let simulation = + Simulation.make ~charge ~multiplicity ~nuclei ao_basis + in + + let hf = + Mo.Hartree_fock.make ~guess:`Hcore simulation + in + + let mo_basis = + Mo.Basis.of_hartree_fock hf + in + + let localized_mo_basis = + Mo.Localization.make + ~kind:Mo.Localization.Boys + mo_basis + [4;5;6;7;8] + |> Mo.Localization.to_basis + in + + Format.printf "%a" (Mo.Basis.pp ~start:1 ~finish:10) localized_mo_basis + + +(* +open Common +open Alcotest + let wd = Qcaml.root ^ Filename.dir_sep ^ "test" in + +let test_xyz molecule length repulsion charge core = + let xyz = Nuclei.of_xyz_file (wd^Filename.dir_sep^molecule^".xyz") in + check int "length" length (Array.length xyz); + check (float 1.e-4) "repulsion" repulsion (Nuclei.repulsion xyz); + check int "charge" charge (Charge.to_int @@ Nuclei.charge xyz); + check int "small_core" core (Nuclei.small_core xyz); + () + +let tests = [ + "caffeine", `Quick, (fun () -> test_xyz "caffeine" 24 917.0684 102 28); + "water", `Quick, (fun () -> test_xyz "water" 3 9.19497 10 2); +] + ,*) + + #+end_src diff --git a/top/lib/install_printers.ml b/top/lib/install_printers.ml index f8fb0b5..3573de7 100644 --- a/top/lib/install_printers.ml +++ b/top/lib/install_printers.ml @@ -15,6 +15,7 @@ let printers = "Gaussian.Atomic_shell_pair.pp" ; "Gaussian.Atomic_shell_pair_couple.pp" ; "Mo.Frozen_core.pp" ; + "Mo.Localization.pp" ; "Particles.Electrons.pp" ; "Particles.Element.pp" ; "Particles.Nuclei.pp" ;