From d6ef25d55ac5254102066ff725cf94c915ea6c38 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 26 Jan 2024 11:24:56 +0100 Subject: [PATCH] Removing org-mode --- mo/frozen_core.org | 118 ------------------------------------ mo/lib/frozen_core.ml | 45 +++----------- mo/lib/frozen_core.mli | 54 +++++++++++------ mo/lib/localization.ml | 86 ++++++++++++-------------- mo/lib/localization.mli | 33 +++++----- mo/test/localization.ml | 72 ++++++++++++++++++++++ top/lib/install_printers.ml | 12 +--- 7 files changed, 173 insertions(+), 247 deletions(-) delete mode 100644 mo/frozen_core.org create mode 100644 mo/test/localization.ml diff --git a/mo/frozen_core.org b/mo/frozen_core.org deleted file mode 100644 index 2b8471f..0000000 --- a/mo/frozen_core.org +++ /dev/null @@ -1,118 +0,0 @@ -#+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 - -* Frozen core - :PROPERTIES: - :header-args: :noweb yes :comments both - :END: - - Defines how the core electrons are frozen, for each atom. - -** Type - - #+NAME: types - #+begin_src ocaml :tangle (eval mli) -type kind = - | All_electron - | Small - | Large - #+end_src - - #+begin_src ocaml :tangle (eval mli) -type t - #+end_src - - #+begin_src ocaml :tangle (eval ml) :exports none -<> -type t = int array - #+end_src - -** Creation - - #+begin_src ocaml :tangle (eval mli) -val make : kind -> Particles.Nuclei.t -> t - -val of_int_list : int list -> t -val of_int_array : int array -> t - #+end_src - - | ~make~ | Creates a ~Frozen_core.t~ with the same kind for all atoms | - | ~of_int_array~ | Creates a ~Frozen_core.t~ giving the number of frozen electrons per atom | - | ~of_int_list~ | Creates a ~Frozen_core.t~ giving the number of frozen electrons per atom | - - #+begin_example -let f = Frozen_core.(make Small nuclei) ;; -val f : Frozen_core.t = [|0; 2; 2; 0|] - -let f = Frozen_core.(of_int_list [0; 2; 2; 0]) -val f : Frozen_core.t = [|0; 2; 2; 0|] - #+end_example - - #+begin_src ocaml :tangle (eval ml) :exports none -let make_ae nuclei = - Array.map (fun _ -> 0) nuclei - -let make_small nuclei = - Array.map (fun (e,_) -> Particles.Element.small_core e) nuclei - -let make_large nuclei = - Array.map (fun (e,_) -> Particles.Element.large_core e) nuclei - -let make = function - | All_electron -> make_ae - | Small -> make_small - | Large -> make_large - - -external of_int_array : int array -> t = "%identity" - -let of_int_list = Array.of_list - #+end_src - -** Access - - #+begin_src ocaml :tangle (eval mli) -val num_elec : t -> int -val num_mos : t -> int - #+end_src - - | ~num_elec~ | Number of frozen electrons | - | ~num_mos~ | Number of frozen molecular orbitals | - - #+begin_example -Frozen_core.num_elec f ;; -- : int = 4 - -Frozen_core.num_mos f ;; -- : int = 2 - #+end_example - - #+begin_src ocaml :tangle (eval ml) :exports none -let num_elec t = - Array.fold_left ( + ) 0 t - -let num_mos t = - (num_elec t) / 2 - #+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 -let pp ppf t = - Format.fprintf ppf "@[[|"; - Array.iter (fun x -> Format.fprintf ppf "@,@[%d@]" x) t; - Format.fprintf ppf "|]@]"; - #+end_src - diff --git a/mo/lib/frozen_core.ml b/mo/lib/frozen_core.ml index 19735bf..e17e5a1 100644 --- a/mo/lib/frozen_core.ml +++ b/mo/lib/frozen_core.ml @@ -1,27 +1,14 @@ -(* [[file:~/QCaml/mo/frozen_core.org::*Type][Type:3]] *) +(** Type *) + type kind = | All_electron | Small | Large type t = int array -(* Type:3 ends here *) +(** Creation *) -(* | ~make~ | Creates a ~Frozen_core.t~ with the same kind for all atoms | - * | ~of_int_array~ | Creates a ~Frozen_core.t~ giving the number of frozen electrons per atom | - * | ~of_int_list~ | Creates a ~Frozen_core.t~ giving the number of frozen electrons per atom | - * - * #+begin_example - * let f = Frozen_core.(make Small nuclei) ;; - * val f : Frozen_core.t = [|0; 2; 2; 0|] - * - * let f = Frozen_core.(of_int_list [0; 2; 2; 0]) - * val f : Frozen_core.t = [|0; 2; 2; 0|] - * #+end_example *) - - -(* [[file:~/QCaml/mo/frozen_core.org::*Creation][Creation:2]] *) let make_ae nuclei = Array.map (fun _ -> 0) nuclei @@ -32,41 +19,29 @@ let make_large nuclei = Array.map (fun (e,_) -> Particles.Element.large_core e) nuclei let make = function - | All_electron -> make_ae - | Small -> make_small + | All_electron -> make_ae + | Small -> make_small | Large -> make_large external of_int_array : int array -> t = "%identity" let of_int_list = Array.of_list -(* Creation:2 ends here *) +(** Access *) -(* | ~num_elec~ | Number of frozen electrons | - * | ~num_mos~ | Number of frozen molecular orbitals | - * - * #+begin_example - * Frozen_core.num_elec f ;; - * - : int = 4 - * - * Frozen_core.num_mos f ;; - * - : int = 2 - * #+end_example *) - - -(* [[file:~/QCaml/mo/frozen_core.org::*Access][Access:2]] *) let num_elec t = Array.fold_left ( + ) 0 t let num_mos t = (num_elec t) / 2 -(* Access:2 ends here *) -(* [[file:~/QCaml/mo/frozen_core.org::*Printers][Printers:2]] *) + +(** Printers *) + let pp ppf t = Format.fprintf ppf "@[[|"; Array.iter (fun x -> Format.fprintf ppf "@,@[%d@]" x) t; Format.fprintf ppf "|]@]"; -(* Printers:2 ends here *) + diff --git a/mo/lib/frozen_core.mli b/mo/lib/frozen_core.mli index f6e65bb..cc5330a 100644 --- a/mo/lib/frozen_core.mli +++ b/mo/lib/frozen_core.mli @@ -1,39 +1,55 @@ -(* Type - * - * #+NAME: types *) +(** Type *) + +(** Defines how the core electrons are frozen, for each atom. *) + -(* [[file:~/QCaml/mo/frozen_core.org::types][types]] *) type kind = | All_electron | Small | Large -(* types ends here *) -(* [[file:~/QCaml/mo/frozen_core.org::*Type][Type:2]] *) type t -(* Type:2 ends here *) - -(* Creation *) -(* [[file:~/QCaml/mo/frozen_core.org::*Creation][Creation:1]] *) +(** Creation *) + +(** Example + * + * let f = Frozen_core.(make Small nuclei) ;; + * val f : Frozen_core.t = [|0; 2; 2; 0|] + * + * let f = Frozen_core.(of_int_list [0; 2; 2; 0]) + * val f : Frozen_core.t = [|0; 2; 2; 0|] + * + *) + val make : kind -> Particles.Nuclei.t -> t +(** Creates a ~Frozen_core.t~ with the same kind for all atoms *) val of_int_list : int list -> t +(** Creates a ~Frozen_core.t~ giving the number of frozen electrons per atom *) + val of_int_array : int array -> t -(* Creation:1 ends here *) - -(* Access *) +(** Creates a ~Frozen_core.t~ giving the number of frozen electrons per atom *) -(* [[file:~/QCaml/mo/frozen_core.org::*Access][Access:1]] *) +(** Access *) + +(** Example + * + * Frozen_core.num_elec f ;; + * - : int = 4 + * + * Frozen_core.num_mos f ;; + * + *) + val num_elec : t -> int +(** Number of frozen electrons *) + val num_mos : t -> int -(* Access:1 ends here *) +(** Number of frozen molecular orbitals *) -(* Printers *) +(** Printers *) - -(* [[file:~/QCaml/mo/frozen_core.org::*Printers][Printers:1]] *) val pp : Format.formatter -> t -> unit -(* Printers:1 ends here *) diff --git a/mo/lib/localization.ml b/mo/lib/localization.ml index 31850ce..bd612c1 100644 --- a/mo/lib/localization.ml +++ b/mo/lib/localization.ml @@ -1,6 +1,8 @@ -(* [[file:~/QCaml/mo/localization.org::*Type][Type:3]] *) open Linear_algebra - +open Common + +(** Types *) + type localization_kind = | Edmiston | Boys @@ -18,7 +20,7 @@ type localization_data = convergence : float ; iteration : int ; } - + type t = { kind : localization_kind ; @@ -26,14 +28,10 @@ type 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]] *) +(** Edmiston-Rudenberg *) + let kappa_edmiston in_basis m_C = let ao_basis = Basis.simulation in_basis @@ -117,7 +115,7 @@ let kappa_edmiston in_basis m_C = ) (Util.array_range 1 n_ao); - let f i j = + let f i j = if i=j then 0. else @@ -133,12 +131,10 @@ let kappa_edmiston in_basis m_C = 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]] *) +(** Boys *) + let kappa_boys in_basis = let ao_basis = Basis.simulation in_basis @@ -147,11 +143,11 @@ let kappa_boys in_basis = 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 @@ -160,13 +156,13 @@ let kappa_boys in_basis = (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) + 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 @@ -176,18 +172,18 @@ let kappa_boys in_basis = (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) + 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 f i j = let pi = Constants.pi in if i=j then 0. - else + else let x = atan2 (m_b12%:(i,j)) (m_a12%:(i,j)) in if x >= 0. then 0.25 *. (pi -. x) @@ -205,11 +201,9 @@ let kappa_boys in_basis = -(* | ~kappa~ | Returns the $\kappa$ antisymmetric matrix used for the rotation matrix and the value of the localization function | - * | ~make~ | Performs the orbital localization | *) +(** Access *) -(* [[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 @@ -218,7 +212,7 @@ let kappa ~kind = match kind with | Edmiston -> kappa_edmiston | Boys -> kappa_boys - + let n_iterations t = Array.fold_left (fun accu x -> @@ -238,20 +232,20 @@ let ao_basis t = Simulation.ao_basis (simulation t) let make ~kind ?(max_iter=500) ?(convergence=1.e-6) 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 + Matrix.of_col_vecs_list mos_vec_list in - - let next_coef kappa x = + + let next_coef kappa x = let r = Matrix.exponential_antisymmetric kappa in let m_C = Matrix.gemm_nt x r in m_C in - - let data_initial = + + let data_initial = let iteration = 0 and scaling = 0.1 in @@ -262,7 +256,7 @@ let make ~kind ?(max_iter=500) ?(convergence=1.e-6) mo_basis selected_mos = { coefficients ; kappa ; scaling ; convergence ; loc_value ; iteration } in - let iteration data = + let iteration data = let iteration = data.iteration+1 in let new_kappa, new_loc_value = kappa_loc data.coefficients in let new_convergence = abs_float (Matrix.amax new_kappa) in @@ -318,33 +312,34 @@ let make ~kind ?(max_iter=500) ?(convergence=1.e-6) mo_basis selected_mos = Array.init max_iter (fun i -> lazy (loop i)) in { kind ; mo_basis ; data = array_data ; selected_mos } - - -let to_basis t = + + +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 + Matrix.to_col_vecs data.coefficients in selected_mos t |> List.iteri (fun i j -> mo_coef_array.(j-1) <- 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 *) + 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@]@]@." + Format.fprintf ppf "@[%4s@[%5s@]@,@[%16s@]@,@[%16s@]@,@[%11s@]@]@." "" "#" "Localization " "Convergence" "Scaling"; Format.fprintf ppf "@[%4s%s@]@." "" line; Array.iter (fun data -> @@ -372,4 +367,3 @@ let pp ppf t = ) "MO Localization"; Format.fprintf ppf "@[%s@]@.@." (String.make 70 '='); Format.fprintf ppf "@[%a@]@." pp_iterations t; -(* Printers:2 ends here *) diff --git a/mo/lib/localization.mli b/mo/lib/localization.mli index 12dc2b6..d756fb1 100644 --- a/mo/lib/localization.mli +++ b/mo/lib/localization.mli @@ -1,10 +1,10 @@ -(* Type - * - * #+NAME: types *) +(** Orbital localization *) + + +(** Types *) -(* [[file:~/QCaml/mo/localization.org::types][types]] *) open Linear_algebra - + type localization_kind = | Edmiston | Boys @@ -12,17 +12,13 @@ type localization_kind = 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 localization_data type t -(* Type:2 ends here *) -(* Access *) - -(* [[file:~/QCaml/mo/localization.org::*Access][Access:1]] *) +(** Access *) + val kind : t -> localization_kind val simulation : t -> Simulation.t val selected_mos : t -> int list @@ -32,21 +28,20 @@ val kappa : Basis.t -> ( ao,loc) Matrix.t -> (loc,loc) Matrix.t * float +(** Returns the $\kappa$ antisymmetric matrix used for the rotation matrix and the value of the localization function *) val make : kind:localization_kind -> - ?max_iter:int -> - ?convergence:float -> + ?max_iter:int -> + ?convergence:float -> Basis.t -> int list -> t +(** Performs the orbital localization *) val to_basis : t -> Basis.t -(* Access:1 ends here *) - -(* Printers *) -(* [[file:~/QCaml/mo/localization.org::*Printers][Printers:1]] *) +(** Printers *) + val pp : Format.formatter -> t -> unit -(* Printers:1 ends here *) diff --git a/mo/test/localization.ml b/mo/test/localization.ml new file mode 100644 index 0000000..004bd2b --- /dev/null +++ b/mo/test/localization.ml @@ -0,0 +1,72 @@ +(* Tests *) + +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); +] + *) diff --git a/top/lib/install_printers.ml b/top/lib/install_printers.ml index 4ac491a..8d778c1 100644 --- a/top/lib/install_printers.ml +++ b/top/lib/install_printers.ml @@ -1,13 +1,6 @@ -(* [[file:~/QCaml/top/install_printers.org::*Intall printers][Intall printers:3]] *) +(** Intall printers printers:3]] *) let printers = [ - "Common.Powers.pp" ; - "Common.Range.pp" ; - "Common.Spin.pp" ; - "Common.Zkey.pp" ; - "Gaussian.Atomic_shell.pp" ; - "Gaussian.Atomic_shell_pair.pp" ; - "Gaussian.Atomic_shell_pair_couple.pp" ; "Mo.Frozen_core.pp" ; "Mo.Localization.pp" ; "Particles.Electrons.pp" ; @@ -16,7 +9,7 @@ let printers = "Particles.Zmatrix.pp" ; "Perturbation.Mp2.pp" ; "Simulation.pp" ; - + ] let eval_exn str = @@ -34,4 +27,3 @@ let rec install_printers = function let () = if not (install_printers printers) then Format.eprintf "Problem installing QCaml-printers@." -(* Intall printers:3 ends here *)