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">
- +Defines how the core electrons are frozen, for each atom.
type kind = +type kind = | All_electron | Small | Large @@ -293,8 +325,8 @@ Defines how the core electrons are frozen, for each atom.
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|]
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 ;;
val pp : Format.formatter -> t -> unit @@ -390,10 +422,126 @@ Frozen_core.num_mos f ;;
+Molecular orbital localization function. +
+ ++Boys: +
+ ++Edmiston-Rudenberg: +
+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 ++
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 + ++
kappa |
+Returns the \(\kappa\) antisymmetric matrix used for the rotation matrix and the value of the localization function | +
make |
+Performs the orbital localization | +
(* + val pp : Format.formatter -> t -> unit + *) ++
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
+< 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" ;