diff --git a/mo/lib/localization.ml b/mo/lib/localization.ml index bd612c1..31850ce 100644 --- a/mo/lib/localization.ml +++ b/mo/lib/localization.ml @@ -1,8 +1,6 @@ +(* [[file:~/QCaml/mo/localization.org::*Type][Type:3]] *) open Linear_algebra -open Common - -(** Types *) - + type localization_kind = | Edmiston | Boys @@ -20,7 +18,7 @@ type localization_data = convergence : float ; iteration : int ; } - + type t = { kind : localization_kind ; @@ -28,10 +26,14 @@ type t = data : localization_data option lazy_t array ; selected_mos : int list ; } + +open Common +(* Type:3 ends here *) + +(* Edmiston-Rudenberg *) -(** 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 @@ -115,7 +117,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 @@ -131,10 +133,12 @@ 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 *) -(** Boys *) - +(* [[file:~/QCaml/mo/localization.org::*Boys][Boys:1]] *) let kappa_boys in_basis = let ao_basis = Basis.simulation in_basis @@ -143,11 +147,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 @@ -156,13 +160,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 @@ -172,18 +176,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) @@ -201,9 +205,11 @@ 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 @@ -212,7 +218,7 @@ let kappa ~kind = match kind with | Edmiston -> kappa_edmiston | Boys -> kappa_boys - + let n_iterations t = Array.fold_left (fun accu x -> @@ -232,20 +238,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 @@ -256,7 +262,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 @@ -312,34 +318,33 @@ 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 *) - -(** Printers *) - +(* [[file:~/QCaml/mo/localization.org::*Printers][Printers:2]] *) 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 -> @@ -367,3 +372,4 @@ 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 d756fb1..12dc2b6 100644 --- a/mo/lib/localization.mli +++ b/mo/lib/localization.mli @@ -1,10 +1,10 @@ -(** Orbital localization *) - - -(** Types *) +(* Type + * + * #+NAME: types *) +(* [[file:~/QCaml/mo/localization.org::types][types]] *) open Linear_algebra - + type localization_kind = | Edmiston | Boys @@ -12,13 +12,17 @@ type localization_kind = type mo = Mo_dim.t type ao = Ao.Ao_dim.t type loc +(* types ends here *) -type localization_data +(* [[file:~/QCaml/mo/localization.org::*Type][Type:2]] *) +type localization_data type t +(* Type:2 ends here *) +(* Access *) + -(** 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 @@ -28,20 +32,21 @@ 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 *) -(** Printers *) - +(* [[file:~/QCaml/mo/localization.org::*Printers][Printers:1]] *) val pp : Format.formatter -> t -> unit +(* Printers:1 ends here *) diff --git a/mo/test/localization.ml b/mo/test/localization.ml index 004bd2b..b393568 100644 --- a/mo/test/localization.ml +++ b/mo/test/localization.ml @@ -1,5 +1,7 @@ (* Tests *) + +(* [[file:~/QCaml/mo/localization.org::*Tests][Tests:1]] *) let test_localization = let nuclei = @@ -70,3 +72,4 @@ let tests = [ "water", `Quick, (fun () -> test_xyz "water" 3 9.19497 10 2); ] *) +(* Tests:1 ends here *) diff --git a/particles/electrons.org b/particles/electrons.org deleted file mode 100644 index 6f3f1dc..0000000 --- a/particles/electrons.org +++ /dev/null @@ -1,160 +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 - -* Electrons - :PROPERTIES: - :header-args: :noweb yes :comments both - :END: - - Data structure which contains the number of \alpha and \beta electrons. - -** Type - - #+NAME: types - #+begin_src ocaml :tangle (eval mli) -type t - #+end_src - - #+begin_src ocaml :tangle (eval ml) :exports none -type t = { - n_alfa : int ; - n_beta : int ; -} - #+end_src - - #+begin_src ocaml :tangle (eval test-ml) :exports none -open Common -open Particles -open Alcotest - -let test_all () = - let nuclei = -"3 -Water -O 0. 0. 0. -H -0.756950272703377558 0. -0.585882234512562827 -H 0.756950272703377558 0. -0.585882234512562827 -" - |> Nuclei.of_xyz_string - in - let e = Electrons.of_atoms nuclei in - #+end_src - -** Creation - - #+begin_src ocaml :tangle (eval mli) -open Common - -val make : int -> int -> t - -val of_atoms : ?multiplicity:int -> ?charge:int -> Nuclei.t -> t -(* @param multiplicity default is 1 - @param charge default is 0 - @raise Invalid_argument if the spin multiplicity is not compatible with - the molecule and the total charge. -*) - #+end_src - - | ~make~ | ~make n_alfa n_beta~ | - | ~of_atoms~ | Creates the data relative to electrons for a molecular system described by ~Nuclei.t~ for a given total charge and spin multiplicity. | - - #+begin_src ocaml :tangle (eval ml) :exports none -open Common - -let make n_alfa n_beta = - { n_alfa ; n_beta } - - -let of_atoms ?multiplicity:(multiplicity=1) ?charge:(charge=0) nuclei = - let positive_charges = - Array.fold_left (fun accu (e, _) -> accu + Charge.to_int (Element.to_charge e) ) - 0 nuclei - in - let negative_charges = charge - positive_charges in - let n_elec = - negative_charges in - let n_beta = ((n_elec - multiplicity)+1)/2 in - let n_alfa = n_elec - n_beta in - let result = { n_alfa ; n_beta } in - if multiplicity <> (n_alfa - n_beta)+1 then - invalid_arg (__FILE__^": make"); - result - - #+end_src - - #+begin_src ocaml :tangle (eval test-ml) :exports none -check int "of_atoms alfa" 5 (Electrons.n_alfa e); -check int "of_atoms beta" 5 (Electrons.n_beta e); - #+end_src - -** Access - - #+begin_src ocaml :tangle (eval mli) -val charge : t -> Charge.t -val n_elec : t -> int -val n_alfa : t -> int -val n_beta : t -> int -val multiplicity : t -> int - #+end_src - - | ~charge~ | Sum of the charges of the electrons | - | ~n_elec~ | Number of electrons | - | ~n_alfa~ | Number of alpha electrons | - | ~n_beta~ | Number of beta electrons | - | ~multiplicity~ | Spin multiplicity: $2S+1$ | - - #+begin_src ocaml :tangle (eval ml) :exports none -let charge e = - - (e.n_alfa + e.n_beta) - |> Charge.of_int - -let n_alfa t = t.n_alfa - -let n_beta t = t.n_beta - -let n_elec t = t.n_alfa + t.n_beta - -let multiplicity t = t.n_alfa - t.n_beta + 1 - #+end_src - - #+begin_src ocaml :tangle (eval test-ml) :exports none -check int "charge " (-10) (Charge.to_int @@ Electrons.charge e); -check int "n_elec" 10 (Electrons.n_elec e); -check int "multiplicity" 1 (Electrons.multiplicity e); -check int "of_atoms alfa m3" 6 (Electrons.(of_atoms ~multiplicity:3 nuclei |> n_alfa)); -check int "of_atoms beta m3" 4 (Electrons.(of_atoms ~multiplicity:3 nuclei |> n_beta)); -check int "of_atoms n_elec m3" 10 (Electrons.(of_atoms ~multiplicity:3 nuclei |> n_elec)); -check int "of_atoms alfa m2 c1" 5 (Electrons.(of_atoms ~multiplicity:2 ~charge:1 nuclei |> n_alfa)); -check int "of_atoms beta m2 c1" 4 (Electrons.(of_atoms ~multiplicity:2 ~charge:1 nuclei |> n_beta)); -check int "of_atoms beta m2 c1" 9 (Electrons.(of_atoms ~multiplicity:2 ~charge:1 nuclei |> n_elec)); -check int "of_atoms mult m2 c1" 2 (Electrons.(of_atoms ~multiplicity:2 ~charge:1 nuclei |> multiplicity)); -check bool "make" true (Electrons.make 6 4 = Electrons.(of_atoms ~multiplicity:3 nuclei)); - #+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 "@[n_alfa=%d, n_beta=%d@]" t.n_alfa t.n_beta - #+end_src - -** Tests - - #+begin_src ocaml :tangle (eval test-ml) :exports none - () - -let tests = [ - "all", `Quick, test_all -] - #+end_src diff --git a/particles/element.org b/particles/element.org deleted file mode 100644 index 9852b67..0000000 --- a/particles/element.org +++ /dev/null @@ -1,328 +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 - -* Element - :PROPERTIES: - :header-args: :noweb yes :comments both - :END: - - Chemical elements. - -** Type - - #+NAME: types - #+begin_src ocaml :tangle (eval mli) -type t = - |X - |H |He - |Li|Be |B |C |N |O |F |Ne - |Na|Mg |Al|Si|P |S |Cl|Ar - |K |Ca|Sc|Ti|V |Cr|Mn|Fe|Co|Ni|Cu|Zn|Ga|Ge|As|Se|Br|Kr - |Rb|Sr|Y |Zr|Nb|Mo|Tc|Ru|Rh|Pd|Ag|Cd|In|Sn|Sb|Te|I |Xe - |Pt - -exception ElementError of string - -open Common - #+end_src - - #+begin_src ocaml :tangle (eval ml) :exports none -<> - #+end_src - -** Conversion - - #+begin_src ocaml :tangle (eval mli) -val of_string : string -> t -val to_string : t -> string -val to_long_string : t -> string - -val to_int : t -> int -val of_int : int -> t - -val to_charge : t -> Charge.t -val of_charge : Charge.t -> t - #+end_src - - | ~of_string~ | Creates an ~Element.t~ from a chemical symbol or from the full name of the element (case insensitive) | - | ~to_string~ | Gets the chemical symbol of the ~Element.t~ in a string | - | ~to_long_string~ | Gets the full name of the ~Element.t~ in a string | - | ~to_int~ | Convert to the atomic charge, with ~int~ type | - | ~of_int~ | Create from the atomic charge, with ~int~ type | - | ~to_charge~ | Convert to the atomic charge, with ~Charge.t~ type | - | ~of_charge~ | Create from the atomic charge, with ~Charge.t~ type | - - #+begin_example -Element.of_string "Fe" ;; -- : Element.t = Particles.Element.Fe - -Element.of_string "hydrogen" ;; -- : Element.t = Particles.Element.H - -Element.of_string "Kryptonite" ;; -Exception: Particles.Element.ElementError "Element Kryptonite unknown". - -Element.(to_long_string Fe) ;; -- : string = "Iron" - -Element.(to_string Fe);; -- : string = "Fe" - #+end_example - - #+begin_src ocaml :tangle (eval ml) :exports none -let of_string x = - match (String.capitalize_ascii (String.lowercase_ascii x)) with - | "X" | "Dummy" -> X | "H" | "Hydrogen" -> H - | "He" | "Helium" -> He | "Li" | "Lithium" -> Li - | "Be" | "Beryllium" -> Be | "B" | "Boron" -> B - | "C" | "Carbon" -> C | "N" | "Nitrogen" -> N - | "O" | "Oxygen" -> O | "F" | "Fluorine" -> F - | "Ne" | "Neon" -> Ne | "Na" | "Sodium" -> Na - | "Mg" | "Magnesium" -> Mg | "Al" | "Aluminum" -> Al - | "Si" | "Silicon" -> Si | "P" | "Phosphorus" -> P - | "S" | "Sulfur" -> S | "Cl" | "Chlorine" -> Cl - | "Ar" | "Argon" -> Ar | "K" | "Potassium" -> K - | "Ca" | "Calcium" -> Ca | "Sc" | "Scandium" -> Sc - | "Ti" | "Titanium" -> Ti | "V" | "Vanadium" -> V - | "Cr" | "Chromium" -> Cr | "Mn" | "Manganese" -> Mn - | "Fe" | "Iron" -> Fe | "Co" | "Cobalt" -> Co - | "Ni" | "Nickel" -> Ni | "Cu" | "Copper" -> Cu - | "Zn" | "Zinc" -> Zn | "Ga" | "Gallium" -> Ga - | "Ge" | "Germanium" -> Ge | "As" | "Arsenic" -> As - | "Se" | "Selenium" -> Se | "Br" | "Bromine" -> Br - | "Kr" | "Krypton" -> Kr | "Rb" | "Rubidium" -> Rb - | "Sr" | "Strontium" -> Sr | "Y" | "Yttrium" -> Y - | "Zr" | "Zirconium" -> Zr | "Nb" | "Niobium" -> Nb - | "Mo" | "Molybdenum" -> Mo | "Tc" | "Technetium" -> Tc - | "Ru" | "Ruthenium" -> Ru | "Rh" | "Rhodium" -> Rh - | "Pd" | "Palladium" -> Pd | "Ag" | "Silver" -> Ag - | "Cd" | "Cadmium" -> Cd | "In" | "Indium" -> In - | "Sn" | "Tin" -> Sn | "Sb" | "Antimony" -> Sb - | "Te" | "Tellurium" -> Te | "I" | "Iodine" -> I - | "Xe" | "Xenon" -> Xe | "Pt" | "Platinum" -> Pt - | x -> raise (ElementError ("Element "^x^" unknown")) - - -let to_string = function - | X -> "X" | H -> "H" | He -> "He" | Li -> "Li" - | Be -> "Be" | B -> "B" | C -> "C" | N -> "N" - | O -> "O" | F -> "F" | Ne -> "Ne" | Na -> "Na" - | Mg -> "Mg" | Al -> "Al" | Si -> "Si" | P -> "P" - | S -> "S" | Cl -> "Cl" | Ar -> "Ar" | K -> "K" - | Ca -> "Ca" | Sc -> "Sc" | Ti -> "Ti" | V -> "V" - | Cr -> "Cr" | Mn -> "Mn" | Fe -> "Fe" | Co -> "Co" - | Ni -> "Ni" | Cu -> "Cu" | Zn -> "Zn" | Ga -> "Ga" - | Ge -> "Ge" | As -> "As" | Se -> "Se" | Br -> "Br" - | Kr -> "Kr" | Rb -> "Rb" | Sr -> "Sr" | Y -> "Y" - | Zr -> "Zr" | Nb -> "Nb" | Mo -> "Mo" | Tc -> "Tc" - | Ru -> "Ru" | Rh -> "Rh" | Pd -> "Pd" | Ag -> "Ag" - | Cd -> "Cd" | In -> "In" | Sn -> "Sn" | Sb -> "Sb" - | Te -> "Te" | I -> "I" | Xe -> "Xe" | Pt -> "Pt" - -let to_long_string = function - | X -> "Dummy" | H -> "Hydrogen" | He -> "Helium" - | Li -> "Lithium" | Be -> "Beryllium" | B -> "Boron" - | C -> "Carbon" | N -> "Nitrogen" | O -> "Oxygen" - | F -> "Fluorine" | Ne -> "Neon" | Na -> "Sodium" - | Mg -> "Magnesium" | Al -> "Aluminum" | Si -> "Silicon" - | P -> "Phosphorus" | S -> "Sulfur" | Cl -> "Chlorine" - | Ar -> "Argon" | K -> "Potassium" | Ca -> "Calcium" - | Sc -> "Scandium" | Ti -> "Titanium" | V -> "Vanadium" - | Cr -> "Chromium" | Mn -> "Manganese" | Fe -> "Iron" - | Co -> "Cobalt" | Ni -> "Nickel" | Cu -> "Copper" - | Zn -> "Zinc" | Ga -> "Gallium" | Ge -> "Germanium" - | As -> "Arsenic" | Se -> "Selenium" | Br -> "Bromine" - | Kr -> "Krypton" | Rb -> "Rubidium" | Sr -> "Strontium" - | Y -> "Yttrium" | Zr -> "Zirconium" | Nb -> "Niobium" - | Mo -> "Molybdenum" | Tc -> "Technetium" | Ru -> "Ruthenium" - | Rh -> "Rhodium" | Pd -> "Palladium" | Ag -> "Silver" - | Cd -> "Cadmium" | In -> "Indium" | Sn -> "Tin" - | Sb -> "Antimony" | Te -> "Tellurium" | I -> "Iodine" - | Xe -> "Xenon" | Pt -> "Platinum" - - -let to_int = function - | X -> 0 | H -> 1 | He -> 2 | Li -> 3 - | Be -> 4 | B -> 5 | C -> 6 | N -> 7 - | O -> 8 | F -> 9 | Ne -> 10 | Na -> 11 - | Mg -> 12 | Al -> 13 | Si -> 14 | P -> 15 - | S -> 16 | Cl -> 17 | Ar -> 18 | K -> 19 - | Ca -> 20 | Sc -> 21 | Ti -> 22 | V -> 23 - | Cr -> 24 | Mn -> 25 | Fe -> 26 | Co -> 27 - | Ni -> 28 | Cu -> 29 | Zn -> 30 | Ga -> 31 - | Ge -> 32 | As -> 33 | Se -> 34 | Br -> 35 - | Kr -> 36 | Rb -> 37 | Sr -> 38 | Y -> 39 - | Zr -> 40 | Nb -> 41 | Mo -> 42 | Tc -> 43 - | Ru -> 44 | Rh -> 45 | Pd -> 46 | Ag -> 47 - | Cd -> 48 | In -> 49 | Sn -> 50 | Sb -> 51 - | Te -> 52 | I -> 53 | Xe -> 54 | Pt -> 78 - - -let to_charge c = - to_int c |> Charge.of_int - - -let of_int = function - | 0 -> X | 1 -> H | 2 -> He | 3 -> Li - | 4 -> Be | 5 -> B | 6 -> C | 7 -> N - | 8 -> O | 9 -> F | 10 -> Ne | 11 -> Na - | 12 -> Mg | 13 -> Al | 14 -> Si | 15 -> P - | 16 -> S | 17 -> Cl | 18 -> Ar | 19 -> K - | 20 -> Ca | 21 -> Sc | 22 -> Ti | 23 -> V - | 24 -> Cr | 25 -> Mn | 26 -> Fe | 27 -> Co - | 28 -> Ni | 29 -> Cu | 30 -> Zn | 31 -> Ga - | 32 -> Ge | 33 -> As | 34 -> Se | 35 -> Br - | 36 -> Kr | 37 -> Rb | 38 -> Sr | 39 -> Y - | 40 -> Zr | 41 -> Nb | 42 -> Mo | 43 -> Tc - | 44 -> Ru | 45 -> Rh | 46 -> Pd | 47 -> Ag - | 48 -> Cd | 49 -> In | 50 -> Sn | 51 -> Sb - | 52 -> Te | 53 -> I | 54 -> Xe | 78 -> Pt - | x -> raise (ElementError ("Element of charge "^(string_of_int x)^" unknown")) - - -let of_charge c = - Charge.to_int c |> of_int - - #+end_src - -** Database information - - #+begin_src ocaml :tangle (eval mli) -val covalent_radius : t -> Non_negative_float.t -val vdw_radius : t -> Non_negative_float.t -val mass : t -> Mass.t -val small_core : t -> int -val large_core : t -> int - #+end_src - - | ~covalent_radius~ | Covalent radii of the elements, in atomic units | - | ~vdw_radius~ | Van der Waals radii of the elements, in atomic units | - | ~mass~ | Atomic mass of the elements, in atomic units) | - | ~small_core~ | Number of electrons in the small core model (all except the outermost two shells) | - | ~large_core~ | Number of electrons in the large core model (all except the outermost shell) | - - #+begin_src ocaml :tangle (eval ml) :exports none -let covalent_radius x = - let result = function - | X -> 0. | H -> 0.37 | He -> 0.70 | Li -> 1.23 - | Be -> 0.89 | B -> 0.90 | C -> 0.85 | N -> 0.74 - | O -> 0.74 | F -> 0.72 | Ne -> 0.70 | Na -> 1.00 - | Mg -> 1.36 | Al -> 1.25 | Si -> 1.17 | P -> 1.10 - | S -> 1.10 | Cl -> 0.99 | Ar -> 0.70 | K -> 2.03 - | Ca -> 1.74 | Sc -> 1.44 | Ti -> 1.32 | V -> 1.22 - | Cr -> 0.00 | Mn -> 1.16 | Fe -> 0.00 | Co -> 1.15 - | Ni -> 1.17 | Cu -> 1.25 | Zn -> 1.25 | Ga -> 1.20 - | Ge -> 1.21 | As -> 1.16 | Se -> 0.70 | Br -> 1.24 - | Kr -> 1.91 | Rb -> 2.20 | Sr -> 1.95 | Y -> 1.90 - | Zr -> 1.75 | Nb -> 1.64 | Mo -> 1.54 | Tc -> 1.47 - | Ru -> 1.46 | Rh -> 1.42 | Pd -> 1.39 | Ag -> 1.45 - | Cd -> 1.44 | In -> 1.42 | Sn -> 1.39 | Sb -> 1.39 - | Te -> 1.38 | I -> 1.39 | Xe -> 1.40 | Pt -> 1.30 - in - Constants.a0 *. (result x) - |> Non_negative_float.of_float - - -let vdw_radius x = - let result = function - | X -> 0. | H -> 1.20 | He -> 1.70 | Li -> 1.70 - | Be -> 1.70 | B -> 1.70 | C -> 1.70 | N -> 1.55 - | O -> 1.52 | F -> 1.47 | Ne -> 1.70 | Na -> 1.70 - | Mg -> 1.70 | Al -> 1.94 | Si -> 2.10 | P -> 1.80 - | S -> 1.80 | Cl -> 1.75 | Ar -> 1.70 | K -> 1.70 - | Ca -> 1.70 | Sc -> 1.70 | Ti -> 1.70 | V -> 1.98 - | Cr -> 1.94 | Mn -> 1.93 | Fe -> 1.93 | Co -> 1.92 - | Ni -> 1.70 | Cu -> 1.70 | Zn -> 1.70 | Ga -> 2.02 - | Ge -> 1.70 | As -> 1.96 | Se -> 1.70 | Br -> 2.10 - | Kr -> 1.70 | Rb -> 3.03 | Sr -> 2.49 | Y -> 0. - | Zr -> 0. | Nb -> 0. | Mo -> 0. | Tc -> 0. - | Ru -> 0. | Rh -> 0. | Pd -> 1.63 | Ag -> 1.72 - | Cd -> 1.58 | In -> 1.93 | Sn -> 2.17 | Sb -> 2.06 - | Te -> 2.06 | I -> 1.98 | Xe -> 2.16 | Pt -> 1.75 - in - Constants.a0 *. (result x) - |> Non_negative_float.of_float - - -let mass c = - begin - match c with - | X -> 0. | H -> 1.0079 | He -> 4.00260 | Li -> 6.941 - | Be -> 9.01218 | B -> 10.81 | C -> 12.011 | N -> 14.0067 - | O -> 15.9994 | F -> 18.998403 | Ne -> 20.179 | Na -> 22.98977 - | Mg -> 24.305 | Al -> 26.98154 | Si -> 28.0855 | P -> 30.97376 - | S -> 32.06 | Cl -> 35.453 | Ar -> 39.948 | K -> 39.0983 - | Ca -> 40.08 | Sc -> 44.9559 | Ti -> 47.90 | V -> 50.9415 - | Cr -> 51.996 | Mn -> 54.9380 | Fe -> 55.9332 | Co -> 58.9332 - | Ni -> 58.70 | Cu -> 63.546 | Zn -> 65.38 | Ga -> 69.72 - | Ge -> 72.59 | As -> 74.9216 | Se -> 78.96 | Br -> 79.904 - | Kr -> 83.80 | Rb -> 85.4678 | Sr -> 87.62 | Y -> 88.90584 - | Zr -> 91.224 | Nb -> 92.90637 | Mo -> 95.95 | Tc -> 98. - | Ru -> 101.07 | Rh -> 102.90550 | Pd -> 106.42 | Ag -> 107.8682 - | Cd -> 112.414 | In -> 114.818 | Sn -> 118.710 | Sb -> 121.760 - | Te -> 127.60 | I -> 126.90447 | Xe -> 131.293 | Pt -> 195.084 - end - |> Mass.of_float - - -let noble_gas = - [ He ; Ne ; Ar ; Kr ; Xe ] - - -let large_core t = - let num = to_int t in - let rec loop = function - | gas :: rest -> - if gas < num then - gas - else - loop rest - | [] -> 0 - in - List.rev_map to_int noble_gas - |> loop - - -let small_core t = - let num = to_int t in - let rec loop = function - | large :: small :: rest -> - if large < num then - small - else - loop (small :: rest) - | small :: [] -> - if small < num then - small - else - 0 - | [] -> 0 - in - List.rev_map to_int noble_gas - |> loop - - #+end_src - -** Printers - - #+begin_src ocaml :tangle (eval mli) -val pp : Format.formatter -> t -> unit -val pp_long : Format.formatter -> t -> unit - #+end_src - - #+begin_src ocaml :tangle (eval ml) :exports none -let pp ppf t = - Format.fprintf ppf "@[%s@]" (to_string t) - -let pp_long ppf t = - Format.fprintf ppf "@[%s@]" (to_long_string t) - #+end_src - diff --git a/particles/lib/electrons.ml b/particles/lib/electrons.ml index a3d6784..7ee820b 100644 --- a/particles/lib/electrons.ml +++ b/particles/lib/electrons.ml @@ -1,20 +1,16 @@ -(* [[file:~/QCaml/particles/electrons.org::*Type][Type:2]] *) +(** Type *) type t = { n_alfa : int ; n_beta : int ; } -(* Type:2 ends here *) -(* | ~make~ | ~make n_alfa n_beta~ | - * | ~of_atoms~ | Creates the data relative to electrons for a molecular system described by ~Nuclei.t~ for a given total charge and spin multiplicity. | *) - - -(* [[file:~/QCaml/particles/electrons.org::*Creation][Creation:2]] *) open Common -let make n_alfa n_beta = +(** Creation *) + +let make n_alfa n_beta = { n_alfa ; n_beta } @@ -31,18 +27,9 @@ let of_atoms ?multiplicity:(multiplicity=1) ?charge:(charge=0) nuclei = if multiplicity <> (n_alfa - n_beta)+1 then invalid_arg (__FILE__^": make"); result -(* Creation:2 ends here *) - -(* | ~charge~ | Sum of the charges of the electrons | - * | ~n_elec~ | Number of electrons | - * | ~n_alfa~ | Number of alpha electrons | - * | ~n_beta~ | Number of beta electrons | - * | ~multiplicity~ | Spin multiplicity: $2S+1$ | *) - - -(* [[file:~/QCaml/particles/electrons.org::*Access][Access:2]] *) +(** Access *) let charge e = - (e.n_alfa + e.n_beta) |> Charge.of_int @@ -54,9 +41,10 @@ let n_beta t = t.n_beta let n_elec t = t.n_alfa + t.n_beta let multiplicity t = t.n_alfa - t.n_beta + 1 -(* Access:2 ends here *) -(* [[file:~/QCaml/particles/electrons.org::*Printers][Printers:2]] *) + +(** Printers *) + let pp ppf t = Format.fprintf ppf "@[n_alfa=%d, n_beta=%d@]" t.n_alfa t.n_beta -(* Printers:2 ends here *) + diff --git a/particles/lib/electrons.mli b/particles/lib/electrons.mli index 08609a8..e6bb65d 100644 --- a/particles/lib/electrons.mli +++ b/particles/lib/electrons.mli @@ -1,41 +1,43 @@ -(* Type - * - * #+NAME: types *) +(** Data structure containing the number of alpha and beta electrons. + *) -(* [[file:~/QCaml/particles/electrons.org::types][types]] *) type t -(* types ends here *) -(* Creation *) - - -(* [[file:~/QCaml/particles/electrons.org::*Creation][Creation:1]] *) open Common +(** Creation *) + val make : int -> int -> t +(* make n_alfa n_beta *) val of_atoms : ?multiplicity:int -> ?charge:int -> Nuclei.t -> t -(* @param multiplicity default is 1 +(** Creates the data relative to electrons for a molecular system described by + Nuclei.t for a given total charge and spin multiplicity. + + @param multiplicity default is 1 @param charge default is 0 @raise Invalid_argument if the spin multiplicity is not compatible with the molecule and the total charge. *) -(* Creation:1 ends here *) -(* Access *) +(** Access *) - -(* [[file:~/QCaml/particles/electrons.org::*Access][Access:1]] *) val charge : t -> Charge.t +(** Sum of the charges of the electrons *) + val n_elec : t -> int +(* Number of electrons *) + val n_alfa : t -> int +(* Number of alpha electrons *) + val n_beta : t -> int +(* Number of beta electrons *) + val multiplicity : t -> int -(* Access:1 ends here *) - -(* Printers *) +(* Spin multiplicity: 2S+1 *) -(* [[file:~/QCaml/particles/electrons.org::*Printers][Printers:1]] *) +(** Printers *) + val pp : Format.formatter -> t -> unit -(* Printers:1 ends here *) diff --git a/particles/lib/element.ml b/particles/lib/element.ml index b163d24..00342be 100644 --- a/particles/lib/element.ml +++ b/particles/lib/element.ml @@ -1,4 +1,5 @@ -(* [[file:~/QCaml/particles/element.org::*Type][Type:2]] *) +(** Type *) + type t = |X |H |He @@ -6,43 +7,16 @@ type t = |Na|Mg |Al|Si|P |S |Cl|Ar |K |Ca|Sc|Ti|V |Cr|Mn|Fe|Co|Ni|Cu|Zn|Ga|Ge|As|Se|Br|Kr |Rb|Sr|Y |Zr|Nb|Mo|Tc|Ru|Rh|Pd|Ag|Cd|In|Sn|Sb|Te|I |Xe - |Pt + |Pt exception ElementError of string open Common -(* Type:2 ends here *) +(** Conversion *) -(* | ~of_string~ | Creates an ~Element.t~ from a chemical symbol or from the full name of the element (case insensitive) | - * | ~to_string~ | Gets the chemical symbol of the ~Element.t~ in a string | - * | ~to_long_string~ | Gets the full name of the ~Element.t~ in a string | - * | ~to_int~ | Convert to the atomic charge, with ~int~ type | - * | ~of_int~ | Create from the atomic charge, with ~int~ type | - * | ~to_charge~ | Convert to the atomic charge, with ~Charge.t~ type | - * | ~of_charge~ | Create from the atomic charge, with ~Charge.t~ type | - * - * #+begin_example - * Element.of_string "Fe" ;; - * - : Element.t = Particles.Element.Fe - * - * Element.of_string "hydrogen" ;; - * - : Element.t = Particles.Element.H - * - * Element.of_string "Kryptonite" ;; - * Exception: Particles.Element.ElementError "Element Kryptonite unknown". - * - * Element.(to_long_string Fe) ;; - * - : string = "Iron" - * - * Element.(to_string Fe);; - * - : string = "Fe" - * #+end_example *) - - -(* [[file:~/QCaml/particles/element.org::*Conversion][Conversion:2]] *) -let of_string x = +let of_string x = match (String.capitalize_ascii (String.lowercase_ascii x)) with | "X" | "Dummy" -> X | "H" | "Hydrogen" -> H | "He" | "Helium" -> He | "Li" | "Lithium" -> Li @@ -130,7 +104,7 @@ let to_int = function | Te -> 52 | I -> 53 | Xe -> 54 | Pt -> 78 -let to_charge c = +let to_charge c = to_int c |> Charge.of_int @@ -154,19 +128,11 @@ let of_int = function let of_charge c = Charge.to_int c |> of_int -(* Conversion:2 ends here *) +(** Database information *) -(* | ~covalent_radius~ | Covalent radii of the elements, in atomic units | - * | ~vdw_radius~ | Van der Waals radii of the elements, in atomic units | - * | ~mass~ | Atomic mass of the elements, in atomic units) | - * | ~small_core~ | Number of electrons in the small core model (all except the outermost two shells) | - * | ~large_core~ | Number of electrons in the large core model (all except the outermost shell) | *) - - -(* [[file:~/QCaml/particles/element.org::*Database information][Database information:2]] *) -let covalent_radius x = +let covalent_radius x = let result = function | X -> 0. | H -> 0.37 | He -> 0.70 | Li -> 1.23 | Be -> 0.89 | B -> 0.90 | C -> 0.85 | N -> 0.74 @@ -187,7 +153,7 @@ let covalent_radius x = |> Non_negative_float.of_float -let vdw_radius x = +let vdw_radius x = let result = function | X -> 0. | H -> 1.20 | He -> 1.70 | Li -> 1.70 | Be -> 1.70 | B -> 1.70 | C -> 1.70 | N -> 1.55 @@ -233,7 +199,7 @@ let noble_gas = [ He ; Ne ; Ar ; Kr ; Xe ] -let large_core t = +let large_core t = let num = to_int t in let rec loop = function | gas :: rest -> @@ -247,7 +213,7 @@ let large_core t = |> loop -let small_core t = +let small_core t = let num = to_int t in let rec loop = function | large :: small :: rest -> @@ -264,12 +230,13 @@ let small_core t = in List.rev_map to_int noble_gas |> loop -(* Database information:2 ends here *) -(* [[file:~/QCaml/particles/element.org::*Printers][Printers:2]] *) + +(** Printers *) + let pp ppf t = Format.fprintf ppf "@[%s@]" (to_string t) let pp_long ppf t = Format.fprintf ppf "@[%s@]" (to_long_string t) -(* Printers:2 ends here *) + diff --git a/particles/lib/element.mli b/particles/lib/element.mli index 60dc08f..eee5952 100644 --- a/particles/lib/element.mli +++ b/particles/lib/element.mli @@ -1,8 +1,5 @@ -(* Type - * - * #+NAME: types *) +(** Chemical elements *) -(* [[file:~/QCaml/particles/element.org::types][types]] *) type t = |X |H |He @@ -10,43 +7,73 @@ type t = |Na|Mg |Al|Si|P |S |Cl|Ar |K |Ca|Sc|Ti|V |Cr|Mn|Fe|Co|Ni|Cu|Zn|Ga|Ge|As|Se|Br|Kr |Rb|Sr|Y |Zr|Nb|Mo|Tc|Ru|Rh|Pd|Ag|Cd|In|Sn|Sb|Te|I |Xe - |Pt + |Pt exception ElementError of string open Common -(* types ends here *) - -(* Conversion *) -(* [[file:~/QCaml/particles/element.org::*Conversion][Conversion:1]] *) +(** Conversion *) + +(** + * Element.of_string "Fe" ;; + * - : Element.t = Particles.Element.Fe + * + * Element.of_string "hydrogen" ;; + * - : Element.t = Particles.Element.H + * + * Element.of_string "Kryptonite" ;; + * Exception: Particles.Element.ElementError "Element Kryptonite unknown". + * + * Element.(to_long_string Fe) ;; + * - : string = "Iron" + * + * Element.(to_string Fe);; + * - : string = "Fe" + *) + val of_string : string -> t -val to_string : t -> string -val to_long_string : t -> string +(** Creates an ~Element.t~ from a chemical symbol or from the full name of the element (case insensitive) *) + +val to_string : t -> string +(** Gets the chemical symbol of the ~Element.t~ in a string *) + +val to_long_string : t -> string +(** Gets the full name of the ~Element.t~ in a string *) + +val to_int : t -> int +(** Convert to the atomic charge, with ~int~ type *) -val to_int : t -> int val of_int : int -> t +(** Create from the atomic charge, with ~int~ type *) val to_charge : t -> Charge.t +(** Convert to the atomic charge, with ~Charge.t~ type *) + val of_charge : Charge.t -> t -(* Conversion:1 ends here *) - -(* Database information *) +(** Create from the atomic charge, with ~Charge.t~ type *) -(* [[file:~/QCaml/particles/element.org::*Database information][Database information:1]] *) +(** Database information *) + val covalent_radius : t -> Non_negative_float.t +(** Covalent radii of the elements, in atomic units *) + val vdw_radius : t -> Non_negative_float.t +(** Van der Waals radii of the elements, in atomic units *) + val mass : t -> Mass.t +(** Atomic mass of the elements, in atomic units) *) + val small_core : t -> int +(** Number of electrons in the small core model (all except the outermost two shells) *) + val large_core : t -> int -(* Database information:1 ends here *) - -(* Printers *) +(** Number of electrons in the large core model (all except the outermost shell) *) -(* [[file:~/QCaml/particles/element.org::*Printers][Printers:1]] *) +(** Printers *) + val pp : Format.formatter -> t -> unit val pp_long : Format.formatter -> t -> unit -(* Printers:1 ends here *) diff --git a/particles/lib/mass.ml b/particles/lib/mass.ml index 8c27f97..68c0278 100644 --- a/particles/lib/mass.ml +++ b/particles/lib/mass.ml @@ -1,3 +1,3 @@ -(* [[file:~/QCaml/particles/mass.org::*Atomic mass][Atomic mass:2]] *) +(** Atomic mass *) + include Common.Non_negative_float -(* Atomic mass:2 ends here *) diff --git a/particles/lib/mass.mli b/particles/lib/mass.mli index f6948d0..9da4522 100644 --- a/particles/lib/mass.mli +++ b/particles/lib/mass.mli @@ -1,12 +1,3 @@ -(* Atomic mass - * :PROPERTIES: - * :header-args: :noweb yes :comments both - * :END: - * - * Atomic mass, a non-negative float. - * - * #+NAME: types *) +(** Atomic mass *) -(* [[file:~/QCaml/particles/mass.org::types][types]] *) include module type of Common.Non_negative_float -(* types ends here *) diff --git a/particles/lib/nuclei.ml b/particles/lib/nuclei.ml index 508cfcf..157a28e 100644 --- a/particles/lib/nuclei.ml +++ b/particles/lib/nuclei.ml @@ -1,21 +1,12 @@ -(* [[file:~/QCaml/particles/nuclei.org::*Type][Type:2]] *) +(** Type *) open Common type t = (Element.t * Coordinate.t) array open Xyz_ast -(* Type:2 ends here *) +(** Conversion *) -(* | ~of_xyz_string~ | Create from a string, in xyz format | - * | ~of_xyz_file~ | Create from a file, in xyz format | - * | ~of_zmt_string~ | Create from a string, in z-matrix format | - * | ~of_zmt_file~ | Create from a file, in z-matrix format | - * | ~to_string~ | Transform to a string, for printing | - * | ~of_filename~ | Detects the type of file (xyz, z-matrix) and reads the file | *) - - -(* [[file:~/QCaml/particles/nuclei.org::*Conversion][Conversion:2]] *) let of_xyz_lexbuf lexbuf = let data = Xyz_parser.input Nuclei_lexer.read_all lexbuf @@ -115,18 +106,10 @@ let to_xyz_string t = ) t |> Array.to_list ) |> String.concat "\n" -(* Conversion:2 ends here *) +(** Query *) -(* | ~formula~ | Returns the chemical formula | - * | ~repulsion~ | Nuclear repulsion energy, in atomic units | - * | ~charge~ | Sum of the charges of the nuclei | - * | ~small_core~ | Number of core electrons in the small core model | - * | ~large_core~ | Number of core electrons in the large core model | *) - - -(* [[file:~/QCaml/particles/nuclei.org::*Query][Query:2]] *) let formula t = let dict = Hashtbl.create 67 in Array.iter (fun (e,_) -> @@ -174,9 +157,10 @@ let small_core a = let large_core a = Array.fold_left (fun accu (e,_) -> accu + (Element.large_core e)) 0 a -(* Query:2 ends here *) -(* [[file:~/QCaml/particles/nuclei.org::*Read][Read:2]] *) + +(** Read *) + let of_trexio f = let num = Trexio.read_nucleus_num f in let charge = Trexio.read_nucleus_charge f @@ -188,9 +172,10 @@ let of_trexio f = z = coord.(3*i+2) } in (Element.of_charge charge.(i), Coordinate.make coord) ) -(* Read:2 ends here *) -(* [[file:~/QCaml/particles/nuclei.org::*Write][Write:2]] *) + +(** Write *) + let to_trexio f t = let num = Array.length t in Trexio.write_nucleus_num f num; @@ -210,9 +195,10 @@ let to_trexio f t = repulsion t |> Trexio.write_nucleus_repulsion f -(* Write:2 ends here *) -(* [[file:~/QCaml/particles/nuclei.org::*Printers][Printers:2]] *) + +(** Printers *) + let pp ppf t = Format.fprintf ppf "@[%s@]" (to_string t) -(* Printers:2 ends here *) + diff --git a/particles/lib/nuclei.mli b/particles/lib/nuclei.mli index 775ce49..ec560be 100644 --- a/particles/lib/nuclei.mli +++ b/particles/lib/nuclei.mli @@ -1,58 +1,64 @@ -(* Type - * <<<~Nuclei.t~>>> - * - * #+NAME: types *) +(** Type *) -(* [[file:~/QCaml/particles/nuclei.org::types][types]] *) open Common type t = (Element.t * Coordinate.t) array -(* types ends here *) -(* Conversion *) +(** Conversion *) -(* [[file:~/QCaml/particles/nuclei.org::*Conversion][Conversion:1]] *) val of_xyz_string : string -> t +(** Create from a string, in xyz format *) + val to_xyz_string : t -> string +(** Transform into a string, in xyz format *) + val of_xyz_file : string -> t +(** Create from a file, in xyz format *) val of_zmt_string : string -> t +(** Create from a string, in z-matrix format *) + val of_zmt_file : string -> t +(** Create from a file, in z-matrix format *) val to_string : t -> string +(** Transform to a string, for printing *) val of_filename : string -> t -(* Conversion:1 ends here *) - -(* Query *) +(** Detects the type of file (xyz, z-matrix) and reads the file *) -(* [[file:~/QCaml/particles/nuclei.org::*Query][Query:1]] *) +(** Query *) + val formula : t -> string +(** Returns the chemical formula *) + val repulsion : t -> float +(** Nuclear repulsion energy, in atomic units *) + val charge : t -> Charge.t +(** Sum of the charges of the nuclei *) + val small_core : t -> int +(** Number of core electrons in the small core model *) + val large_core : t -> int -(* Query:1 ends here *) - -(* Read *) +(** Number of core electrons in the large core model | *) -(* [[file:~/QCaml/particles/nuclei.org::*Read][Read:1]] *) +(** Read *) + val of_trexio : Trexio.trexio_file -> t -(* Read:1 ends here *) +(** Read from a file in TREXIO format *) -(* Write *) +(** Write *) - -(* [[file:~/QCaml/particles/nuclei.org::*Write][Write:1]] *) val to_trexio : Trexio.trexio_file -> t -> unit -(* Write:1 ends here *) - -(* Printers *) +(** Write to a file in TREXIO format *) -(* [[file:~/QCaml/particles/nuclei.org::*Printers][Printers:1]] *) +(** Printers *) + val pp : Format.formatter -> t -> unit -(* Printers:1 ends here *) + diff --git a/particles/lib/nuclei_lexer.mll b/particles/lib/nuclei_lexer.mll index ac36dad..3b85da2 100644 --- a/particles/lib/nuclei_lexer.mll +++ b/particles/lib/nuclei_lexer.mll @@ -1,10 +1,5 @@ -(* Lexer - * - * =nuclei_lexer.mll= contains the description of the lexemes used in - * an xyz file. *) +(** Lexer to read xyz files *) - -(* [[file:~/QCaml/particles/nuclei.org::*Lexer][Lexer:1]] *) { open Xyz_parser } @@ -45,4 +40,3 @@ rule read_all = parse done; *) } -(* Lexer:1 ends here *) diff --git a/particles/lib/xyz_ast.mli b/particles/lib/xyz_ast.mli index 6bd214d..256db13 100644 --- a/particles/lib/xyz_ast.mli +++ b/particles/lib/xyz_ast.mli @@ -1,10 +1,9 @@ +(** When an xyz file is read by =xyz_parser.mly=, it is converted into + * an ~xyz_file~ data structure. *) -(* When an xyz file is read by =xyz_parser.mly=, it is converted into - * an ~xyz_file~ data structure. *) +(** Parser *) - -(* [[file:~/QCaml/particles/nuclei.org::*Parser][Parser:2]] *) open Common type nucleus = @@ -19,4 +18,4 @@ type xyz_file = file_title : string ; nuclei : nucleus list ; } -(* Parser:2 ends here *) + diff --git a/particles/lib/zmatrix.ml b/particles/lib/zmatrix.ml index e82fe89..4e1b112 100644 --- a/particles/lib/zmatrix.ml +++ b/particles/lib/zmatrix.ml @@ -1,4 +1,5 @@ -(* [[file:~/QCaml/particles/zmatrix.org::*Type][Type:2]] *) +(** Type *) + module StringMap = Map.Make(String) type atom_id = int @@ -6,7 +7,7 @@ type angle = Label of string | Value of float type distance = Label of string | Value of float type dihedral = Label of string | Value of float -type line = +type line = | First of Element.t | Second of (Element.t * distance) | Third of (Element.t * atom_id * distance * atom_id * angle) @@ -14,59 +15,10 @@ type line = | Coord of (string * float) type t = (line array * float StringMap.t) -(* Type:2 ends here *) +(** Conversion *) -(* | ~of_string~ | Reads a z-matrix from a string | - * | ~to_xyz~ | Converts to xyz format, as in the ~Nuclei~ module | - * | ~to_xyz_string~ | Converts to xyz format, as a string | - * - * #+begin_example - * let zmt = Zmatrix.of_string " - * n - * n 1 nn - * h 1 hn 2 hnn - * h 2 hn 1 hnn 3 dih4 - * h 1 hn 2 hnn 4 dih5 - * h 2 hn 1 hnn 3 dih5 - * - * nn 1.446 - * hn 1.016 - * hnn 106.0 - * dih4 -54.38 - * dih5 54.38 - * " ;; - * - : Zmatrix.t = N - * N 1 1.446000 - * H 1 1.016000 2 106.000000 - * H 2 1.016000 1 106.000000 3 -54.380000 - * H 1 1.016000 2 106.000000 4 54.380000 - * H 2 1.016000 1 106.000000 3 54.380000 - * - * - * Zmatrix.to_xyz zmt ;; - * - : (Element.t * float * float * float) array = - * [|(N, 0., 0., 0.); (N, 0., 0., 1.446); - * (H, -0.976641883073332107, 0., -0.280047553510071046); - * (H, -0.568802835186988709, 0.793909757123734683, 1.726047553510071); - * (H, 0.314092649983635563, 0.924756819385119, -0.280047553510071101); - * (H, -0.568802835186988709, -0.793909757123734683, 1.726047553510071)|] - * - * - * Zmatrix.to_xyz_string zmt ;; - * - : string = - * "N 0.000000 0.000000 0.000000 - * N 0.000000 0.000000 1.446000 - * H -0.976642 0.000000 -0.280048 - * H -0.568803 0.793910 1.726048 - * H 0.314093 0.924757 -0.280048 - * H -0.568803 -0.793910 1.726048" - * #+end_example *) - - - -(* [[file:~/QCaml/particles/zmatrix.org::*Conversion][Conversion:2]] *) let pi = Common.Constants.pi let to_radian = pi /. 180. @@ -78,43 +30,43 @@ let rec in_range (xmin, xmax) x = else x -let atom_id_of_int : int -> atom_id = +let atom_id_of_int : int -> atom_id = fun x -> ( assert (x>0) ; x) -let distance_of_float : float -> distance = +let distance_of_float : float -> distance = fun x -> ( assert (x>=0.) ; Value x) -let angle_of_float : float -> angle = +let angle_of_float : float -> angle = fun x -> Value (in_range (-180., 180.) x) -let dihedral_of_float : float -> dihedral = +let dihedral_of_float : float -> dihedral = fun x -> Value (in_range (-360., 360.) x) -let atom_id_of_string : string -> atom_id = +let atom_id_of_string : string -> atom_id = fun i -> atom_id_of_int @@ int_of_string i let distance_of_string : string -> distance = - fun s -> + fun s -> try distance_of_float @@ float_of_string s with _ -> Label s let angle_of_string : string -> angle = - fun s -> + fun s -> try angle_of_float @@ float_of_string s with _ -> Label s let dihedral_of_string : string -> dihedral = - fun s -> + fun s -> try dihedral_of_float @@ float_of_string s with _ -> Label s let int_of_atom_id : atom_id -> int = fun x -> x -let float_of_distance : float StringMap.t -> distance -> float = +let float_of_distance : float StringMap.t -> distance -> float = fun map -> function | Value x -> x | Label s -> StringMap.find s map @@ -130,7 +82,7 @@ let float_of_dihedral : float StringMap.t -> dihedral -> float = | Label s -> StringMap.find s map -let string_of_line map = +let string_of_line map = let f_r = float_of_distance map and f_a = float_of_angle map and f_d = float_of_dihedral map @@ -153,7 +105,7 @@ let line_of_string l = | e :: _ :: r :: [] -> Second (Element.of_string e, distance_of_string r) - | e :: i :: r :: j :: a :: [] -> Third + | e :: i :: r :: j :: a :: [] -> Third (Element.of_string e, atom_id_of_string i, distance_of_string r, @@ -179,7 +131,7 @@ let of_string t = |> List.map line_of_string in - let l = + let l = match l with | First _ :: Second _ :: Third _ :: _ | First _ :: Second _ :: Coord _ :: [] @@ -200,7 +152,7 @@ let of_string t = work [] (StringMap.empty) l in (Array.of_list l, m) - + (** Linear algebra *) @@ -225,17 +177,17 @@ let normalized u = let cross (x,y,z) (x',y',z') = ((y *. z' -. z *. y'), -. (x *. z' -. z *. x'), (x *. y' -. y *. x')) -let rotation_matrix axis angle = +let rotation_matrix axis angle = (* Euler-Rodrigues formula for rotation matrix, taken from https://github.com/jevandezande/zmatrix/blob/master/converter.py *) - let a = + let a = (cos (angle *. to_radian *. 0.5)) in - let (b, c, d) = + let (b, c, d) = (-. sin (angle *. to_radian *. 0.5)) |. (normalized axis) in - Array.of_list @@ + Array.of_list @@ [(a *. a +. b *. b -. c *. c -. d *. d, 2. *. (b *. c -. a *. d), 2. *. (b *. d +. a *. c)); @@ -245,16 +197,16 @@ let rotation_matrix axis angle = (2. *. (b *. d -. a *. c), 2. *. (c *. d +. a *. b), a *. a +. d *. d -. b *. b -. c *. c)] - + let apply_rotation_matrix rot u = (dot rot.(0) u, dot rot.(1) u, dot rot.(2) u) - + let to_xyz (z,map) = let result = Array.make (Array.length z) None - in + in let get_cartesian_coord i = match result.(i-1) with @@ -265,14 +217,14 @@ let to_xyz (z,map) = let append_line i' = match z.(i') with - | First e -> + | First e -> result.(i') <- Some (e, 0., 0., 0.) - | Second (e, r) -> + | Second (e, r) -> let r = float_of_distance map r in result.(i') <- Some (e, 0., 0., r) - | Third (e, i, r, j, a) -> + | Third (e, i, r, j, a) -> begin let i, r, j, a = int_of_atom_id i, @@ -281,13 +233,13 @@ let to_xyz (z,map) = float_of_angle map a in let ui, uj = - get_cartesian_coord i, + get_cartesian_coord i, get_cartesian_coord j in - let u_ij = + let u_ij = (uj |- ui) in - let rot = + let rot = rotation_matrix (0., 1., 0.) a in let new_vec = @@ -298,7 +250,7 @@ let to_xyz (z,map) = in result.(i') <- Some (e, x, y, z) end - | Other (e, i, r, j, a, k, d) -> + | Other (e, i, r, j, a, k, d) -> begin let i, r, j, a, k, d = int_of_atom_id i, @@ -309,7 +261,7 @@ let to_xyz (z,map) = float_of_dihedral map d in let ui, uj, uk = - get_cartesian_coord i, + get_cartesian_coord i, get_cartesian_coord j, get_cartesian_coord k in @@ -319,7 +271,7 @@ let to_xyz (z,map) = let normal = cross u_ij u_kj in - let new_vec = + let new_vec = r |. (normalized u_ij) |> apply_rotation_matrix (rotation_matrix normal a) |> apply_rotation_matrix (rotation_matrix u_ij d) @@ -332,7 +284,7 @@ let to_xyz (z,map) = | Coord _ -> () in Array.iteri (fun i _ -> append_line i) z; - let result = + let result = Array.map (function | Some x -> x | None -> failwith "Some atoms were not defined" ) result @@ -341,15 +293,16 @@ let to_xyz (z,map) = let to_xyz_string (l,map) = - String.concat "\n" - ( to_xyz (l,map) - |> Array.map (fun (e,x,y,z) -> + String.concat "\n" + ( to_xyz (l,map) + |> Array.map (fun (e,x,y,z) -> Printf.sprintf "%s %f %f %f" (Element.to_string e) x y z) |> Array.to_list ) -(* Conversion:2 ends here *) -(* [[file:~/QCaml/particles/zmatrix.org::*Printers][Printers:2]] *) + +(** Printers *) + let pp ppf (a, map) = let f = string_of_line map in Format.fprintf ppf "@["; @@ -357,4 +310,4 @@ let pp ppf (a, map) = Format.fprintf ppf "%s@." (f line) ) a; Format.fprintf ppf "@]" -(* Printers:2 ends here *) + diff --git a/particles/lib/zmatrix.mli b/particles/lib/zmatrix.mli index 663a790..04758a7 100644 --- a/particles/lib/zmatrix.mli +++ b/particles/lib/zmatrix.mli @@ -1,23 +1,61 @@ -(* Type - * - * #+NAME: types *) +(** Z-matrix representation of nuclear coordinates *) + +(** let zmt = Zmatrix.of_string " + * n + * n 1 nn + * h 1 hn 2 hnn + * h 2 hn 1 hnn 3 dih4 + * h 1 hn 2 hnn 4 dih5 + * h 2 hn 1 hnn 3 dih5 + * + * nn 1.446 + * hn 1.016 + * hnn 106.0 + * dih4 -54.38 + * dih5 54.38 + * " ;; + * - : Zmatrix.t = N + * N 1 1.446000 + * H 1 1.016000 2 106.000000 + * H 2 1.016000 1 106.000000 3 -54.380000 + * H 1 1.016000 2 106.000000 4 54.380000 + * H 2 1.016000 1 106.000000 3 54.380000 + * + * + * Zmatrix.to_xyz zmt ;; + * - : (Element.t * float * float * float) array = + * [|(N, 0., 0., 0.); (N, 0., 0., 1.446); + * (H, -0.976641883073332107, 0., -0.280047553510071046); + * (H, -0.568802835186988709, 0.793909757123734683, 1.726047553510071); + * (H, 0.314092649983635563, 0.924756819385119, -0.280047553510071101); + * (H, -0.568802835186988709, -0.793909757123734683, 1.726047553510071)|] + * + * + * Zmatrix.to_xyz_string zmt ;; + * - : string = + * "N 0.000000 0.000000 0.000000 + * N 0.000000 0.000000 1.446000 + * H -0.976642 0.000000 -0.280048 + * H -0.568803 0.793910 1.726048 + * H 0.314093 0.924757 -0.280048 + * H -0.568803 -0.793910 1.726048" + *) -(* [[file:~/QCaml/particles/zmatrix.org::types][types]] *) type t -(* types ends here *) - -(* Conversion *) -(* [[file:~/QCaml/particles/zmatrix.org::*Conversion][Conversion:1]] *) +(** Conversion *) + val of_string : string -> t +(** Reads a z-matrix from a string *) + val to_xyz : t -> (Element.t * float * float * float) array +(** Converts to xyz format, as in the ~Nuclei~ module *) + val to_xyz_string : t -> string -(* Conversion:1 ends here *) - -(* Printers *) +(** Converts to xyz format, as a string *) -(* [[file:~/QCaml/particles/zmatrix.org::*Printers][Printers:1]] *) +(** Printers *) + val pp : Format.formatter -> t -> unit -(* Printers:1 ends here *) diff --git a/particles/mass.org b/particles/mass.org deleted file mode 100644 index b727702..0000000 --- a/particles/mass.org +++ /dev/null @@ -1,27 +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 - -* Atomic mass - :PROPERTIES: - :header-args: :noweb yes :comments both - :END: - - Atomic mass, a non-negative float. - - #+NAME: types - #+begin_src ocaml :tangle (eval mli) -include module type of Common.Non_negative_float - #+end_src - - #+begin_src ocaml :tangle (eval ml) :exports none -include Common.Non_negative_float - #+end_src - diff --git a/particles/nuclei.org b/particles/nuclei.org deleted file mode 100644 index 669413a..0000000 --- a/particles/nuclei.org +++ /dev/null @@ -1,471 +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 - -* Nuclei - :PROPERTIES: - :header-args: :noweb yes :comments both - :END: - -** Type - <<<~Nuclei.t~>>> - - #+NAME: types - #+begin_src ocaml :tangle (eval mli) -open Common - -type t = (Element.t * Coordinate.t) array - #+end_src - - #+begin_src ocaml :tangle (eval ml) :exports none -open Common - -type t = (Element.t * Coordinate.t) array -open Xyz_ast - #+end_src - -** xyz file lexer/parser - -*** Lexer - - =nuclei_lexer.mll= contains the description of the lexemes used in - an xyz file. - - #+begin_src ocaml :tangle lib/nuclei_lexer.mll :export none -{ -open Xyz_parser -} - -let eol = ['\n'] -let white = [' ' '\t']+ -let word = [^' ' '\t' '\n']+ -let letter = ['A'-'Z' 'a'-'z'] -let integer = ['0'-'9']+ -let real = '-'? (integer '.' integer | integer '.' | '.' integer) (['e' 'E'] ('+'|'-')? integer)? - - -rule read_all = parse - | eof { EOF } - | eol { EOL } - | white as w { SPACE w } - | integer as i { INTEGER (int_of_string i) } - | real as f { FLOAT (float_of_string f) } - | word as w { WORD w } - - -{ -(* DEBUG - let () = - let ic = open_in "h2o.xyz" in - let lexbuf = Lexing.from_channel ic in - while true do - let s = - match read_all lexbuf with - | EOL -> "EOL" - | SPACE w -> "SPACE("^w^")" - | INTEGER i -> "INTEGER("^(string_of_int i)^")" - | FLOAT f -> "FLOAT("^(string_of_float f)^")" - | WORD w -> "WORD("^w^")" - | EOF -> "EOF" - in - print_endline s - done; -*) -} - #+end_src - -*** Parser - - =xyz_parser.mly= parses nuclear coordinates in xyz format. - #+begin_src ocaml :tangle lib/xyz_parser.mly :export none :comments none -%{ -open Common - -let make_angstrom x y z = - Coordinate.(make_angstrom { - x ; y ; z - }) - -let output_of f x y z = - let a = make_angstrom x y z in - fun e -> - { - Xyz_ast. - element = f e; - coord = a ; - } - -let output_of_string = output_of Element.of_string -let output_of_int = output_of Element.of_int - -%} - -%token EOL -%token SPACE -%token WORD -%token INTEGER -%token FLOAT -%token EOF - -%start input -%type input - -%% /* Grammar rules and actions follow */ - -input: - | integer title atoms_xyz { - { - number_of_atoms = $1; - file_title = $2; - nuclei = $3; - } - } -; - - -integer: - | INTEGER EOL { $1 } - | INTEGER SPACE EOL { $1 } - | SPACE INTEGER EOL { $2 } - | SPACE INTEGER SPACE EOL { $2 } -; - -title: - | title_list EOL { $1 } -; - -text: - | WORD { $1 } - | SPACE { $1 } - | FLOAT { (string_of_float $1)} - | INTEGER { (string_of_int $1)} -; - -title_list: - | { "" } - | title_list text { ($1 ^ $2) } -; - -atoms_xyz: - | atoms_list EOL { List.rev $1 } - | atoms_list EOF { List.rev $1 } -; - -atoms_list: - | { [] } - | atoms_list WORD SPACE FLOAT SPACE FLOAT SPACE FLOAT EOL { output_of_string $4 $6 $8 $2 :: $1 } - | atoms_list WORD SPACE FLOAT SPACE FLOAT SPACE FLOAT SPACE EOL { output_of_string $4 $6 $8 $2 :: $1 } - | atoms_list INTEGER SPACE FLOAT SPACE FLOAT SPACE FLOAT EOL { output_of_int $4 $6 $8 $2 :: $1 } - | atoms_list INTEGER SPACE FLOAT SPACE FLOAT SPACE FLOAT SPACE EOL { output_of_int $4 $6 $8 $2 :: $1 } - | atoms_list SPACE WORD SPACE FLOAT SPACE FLOAT SPACE FLOAT EOL { output_of_string $5 $7 $9 $3 :: $1 } - | atoms_list SPACE WORD SPACE FLOAT SPACE FLOAT SPACE FLOAT SPACE EOL { output_of_string $5 $7 $9 $3 :: $1 } - | atoms_list SPACE INTEGER SPACE FLOAT SPACE FLOAT SPACE FLOAT EOL { output_of_int $5 $7 $9 $3 :: $1 } - | atoms_list SPACE INTEGER SPACE FLOAT SPACE FLOAT SPACE FLOAT SPACE EOL { output_of_int $5 $7 $9 $3 :: $1 } -; - #+end_src - - When an xyz file is read by =xyz_parser.mly=, it is converted into - an ~xyz_file~ data structure. - - #+begin_src ocaml :tangle lib/xyz_ast.mli -open Common - -type nucleus = - { - element: Element.t ; - coord : Coordinate.angstrom Coordinate.point; - } - -type xyz_file = - { - number_of_atoms : int ; - file_title : string ; - nuclei : nucleus list ; - } - #+end_src - -** Conversion - - #+begin_src ocaml :tangle (eval mli) -val of_xyz_string : string -> t -val to_xyz_string : t -> string -val of_xyz_file : string -> t - -val of_zmt_string : string -> t -val of_zmt_file : string -> t - -val to_string : t -> string - -val of_filename : string -> t - #+end_src - - | ~of_xyz_string~ | Create from a string, in xyz format | - | ~of_xyz_file~ | Create from a file, in xyz format | - | ~of_zmt_string~ | Create from a string, in z-matrix format | - | ~of_zmt_file~ | Create from a file, in z-matrix format | - | ~to_string~ | Transform to a string, for printing | - | ~of_filename~ | Detects the type of file (xyz, z-matrix) and reads the file | - - #+begin_src ocaml :tangle (eval ml) :exports none -let of_xyz_lexbuf lexbuf = - let data = - Xyz_parser.input Nuclei_lexer.read_all lexbuf - in - - let len = List.length data.nuclei in - if len <> data.number_of_atoms then - Printf.sprintf "Error: expected %d atoms but %d read" - data.number_of_atoms len - |> failwith; - - List.map (fun nucleus -> - nucleus.element, Coordinate.angstrom_to_bohr nucleus.coord - ) data.nuclei - |> Array.of_list - - -let of_xyz_string input_string = - Lexing.from_string input_string - |> of_xyz_lexbuf - - -let of_xyz_file filename = - let ic = open_in filename in - let lexbuf = - Lexing.from_channel ic - in - let result = - of_xyz_lexbuf lexbuf - in - close_in ic; - result - - -let of_zmt_string buffer = - Zmatrix.of_string buffer - |> Zmatrix.to_xyz - |> Array.map (fun (e,x,y,z) -> - (e, Coordinate.(angstrom_to_bohr @@ make_angstrom { x ; y ; z} )) - ) - - -let of_zmt_file filename = - let ic = open_in filename in - let rec aux accu = - try - let line = input_line ic in - aux (line::accu) - with End_of_file -> - close_in ic; - List.rev accu - |> String.concat "\n" - in aux [] - |> of_zmt_string - - -let to_string atoms = - " - Nuclear Coordinates (Angstrom) - ------------------------------ - ------------------------------------------------------------------------ - Center Atomic Element Coordinates (Angstroms) - Number X Y Z ------------------------------------------------------------------------ -" ^ - (Array.mapi (fun i (e, coord) -> - let open Coordinate in - let coord = - bohr_to_angstrom coord - in - Printf.sprintf " %5d %5d %5s %12.6f %12.6f %12.6f" - (i+1) (Element.to_int e) (Element.to_string e) - coord.x coord.y coord.z - ) atoms - |> Array.to_list - |> String.concat "\n" ) ^ - " ------------------------------------------------------------------------ - -" - - -let of_filename filename = - of_xyz_file filename - - -let to_xyz_string t = - [ string_of_int (Array.length t) ; "" ] @ - ( Array.map (fun (e, coord) -> - let open Coordinate in - let coord = - bohr_to_angstrom coord - in - Printf.sprintf " %5s %12.6f %12.6f %12.6f" - (Element.to_string e) coord.x coord.y coord.z - ) t - |> Array.to_list ) - |> String.concat "\n" - - - #+end_src - -** Query - - #+begin_src ocaml :tangle (eval mli) -val formula : t -> string -val repulsion : t -> float -val charge : t -> Charge.t -val small_core : t -> int -val large_core : t -> int - #+end_src - - | ~formula~ | Returns the chemical formula | - | ~repulsion~ | Nuclear repulsion energy, in atomic units | - | ~charge~ | Sum of the charges of the nuclei | - | ~small_core~ | Number of core electrons in the small core model | - | ~large_core~ | Number of core electrons in the large core model | - - #+begin_src ocaml :tangle (eval ml) :exports none -let formula t = - let dict = Hashtbl.create 67 in - Array.iter (fun (e,_) -> - let e = Element.to_string e in - let value = - try (Hashtbl.find dict e) + 1 - with Not_found -> 1 - in - Hashtbl.replace dict e value - ) t; - Hashtbl.to_seq_keys dict - |> List.of_seq - |> List.sort String.compare - |> List.fold_left (fun accu key -> - let x = Hashtbl.find dict key in - accu ^ key ^ "_{" ^ (string_of_int x) ^ "}") "" - - - -let repulsion nuclei = - let get_charge e = - Element.to_charge e - |> Charge.to_float - in - Array.fold_left ( fun accu (e1, coord1) -> - accu +. - Array.fold_left (fun accu (e2, coord2) -> - let r = Coordinate.(norm (coord1 |- coord2)) in - if r > 0. then - accu +. 0.5 *. (get_charge e2) *. (get_charge e1) /. r - else accu - ) 0. nuclei - ) 0. nuclei - - -let charge nuclei = - Array.fold_left (fun accu (e, _) -> accu + Charge.to_int (Element.to_charge e) ) - 0 nuclei - |> Charge.of_int - - -let small_core a = - Array.fold_left (fun accu (e,_) -> accu + (Element.small_core e)) 0 a - - -let large_core a = - Array.fold_left (fun accu (e,_) -> accu + (Element.large_core e)) 0 a - - #+end_src - -** TREXIO - -*** Read - - #+begin_src ocaml :tangle (eval mli) -val of_trexio : Trexio.trexio_file -> t - #+end_src - - #+begin_src ocaml :tangle (eval ml) :exports none -let of_trexio f = - let num = Trexio.read_nucleus_num f in - let charge = Trexio.read_nucleus_charge f - |> Array.map Charge.of_float in - let coord = Trexio.read_nucleus_coord f in - Array.init num (fun i -> - let coord = Coordinate.{ x = coord.(3*i) ; - y = coord.(3*i+1) ; - z = coord.(3*i+2) } in - (Element.of_charge charge.(i), Coordinate.make coord) - ) - #+end_src - -*** Write - - #+begin_src ocaml :tangle (eval mli) -val to_trexio : Trexio.trexio_file -> t -> unit - #+end_src - - #+begin_src ocaml :tangle (eval ml) :exports none -let to_trexio f t = - let num = Array.length t in - Trexio.write_nucleus_num f num; - - Array.map (fun (e, _) -> Element.to_charge e |> Charge.to_float) t - |> Trexio.write_nucleus_charge f; - - Array.map (fun (e, _) -> Element.to_string e) t - |> Trexio.write_nucleus_label f; - - let coord = Array.init (num*3) (fun _ -> 0.) in - Array.iteri (fun i (_, xyz) -> - coord.(3*i) <- Coordinate.(get X xyz) ; - coord.(3*i+1) <- Coordinate.(get Y xyz) ; - coord.(3*i+2) <- Coordinate.(get Z xyz) ) t; - Trexio.write_nucleus_coord f coord; - - repulsion t - |> Trexio.write_nucleus_repulsion f - - #+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 "@[%s@]" (to_string t) - #+end_src - -** Tests - - #+begin_src ocaml :tangle (eval test-ml) :exports none -open Common -open Particles -open Alcotest - -let wd = Common.Qcaml.root ^ Filename.dir_sep ^ "test" - -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/particles/zmatrix.org b/particles/zmatrix.org deleted file mode 100644 index 610da5b..0000000 --- a/particles/zmatrix.org +++ /dev/null @@ -1,397 +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 - -* Z-matrix - :PROPERTIES: - :header-args: :noweb yes :comments both - :END: - - Z-matrix representation of nuclear coordinates. - -** Type - - #+NAME: types - #+begin_src ocaml :tangle (eval mli) -type t - #+end_src - - #+begin_src ocaml :tangle (eval ml) :exports none -module StringMap = Map.Make(String) - -type atom_id = int -type angle = Label of string | Value of float -type distance = Label of string | Value of float -type dihedral = Label of string | Value of float - -type line = -| First of Element.t -| Second of (Element.t * distance) -| Third of (Element.t * atom_id * distance * atom_id * angle) -| Other of (Element.t * atom_id * distance * atom_id * angle * atom_id * dihedral ) -| Coord of (string * float) - -type t = (line array * float StringMap.t) - #+end_src - -** Conversion - - #+begin_src ocaml :tangle (eval mli) -val of_string : string -> t -val to_xyz : t -> (Element.t * float * float * float) array -val to_xyz_string : t -> string - #+end_src - - | ~of_string~ | Reads a z-matrix from a string | - | ~to_xyz~ | Converts to xyz format, as in the ~Nuclei~ module | - | ~to_xyz_string~ | Converts to xyz format, as a string | - - #+begin_example -let zmt = Zmatrix.of_string " - n - n 1 nn - h 1 hn 2 hnn - h 2 hn 1 hnn 3 dih4 - h 1 hn 2 hnn 4 dih5 - h 2 hn 1 hnn 3 dih5 - -nn 1.446 -hn 1.016 -hnn 106.0 -dih4 -54.38 -dih5 54.38 -" ;; -- : Zmatrix.t = N -N 1 1.446000 -H 1 1.016000 2 106.000000 -H 2 1.016000 1 106.000000 3 -54.380000 -H 1 1.016000 2 106.000000 4 54.380000 -H 2 1.016000 1 106.000000 3 54.380000 - - -Zmatrix.to_xyz zmt ;; -- : (Element.t * float * float * float) array = -[|(N, 0., 0., 0.); (N, 0., 0., 1.446); - (H, -0.976641883073332107, 0., -0.280047553510071046); - (H, -0.568802835186988709, 0.793909757123734683, 1.726047553510071); - (H, 0.314092649983635563, 0.924756819385119, -0.280047553510071101); - (H, -0.568802835186988709, -0.793909757123734683, 1.726047553510071)|] - - -Zmatrix.to_xyz_string zmt ;; -- : string = -"N 0.000000 0.000000 0.000000 -N 0.000000 0.000000 1.446000 -H -0.976642 0.000000 -0.280048 -H -0.568803 0.793910 1.726048 -H 0.314093 0.924757 -0.280048 -H -0.568803 -0.793910 1.726048" - #+end_example - - - #+begin_src ocaml :tangle (eval ml) :exports none -let pi = Common.Constants.pi -let to_radian = pi /. 180. - -let rec in_range (xmin, xmax) x = - if (x <= xmin) then - in_range (xmin, xmax) (x -. xmin +. xmax ) - else if (x > xmax) then - in_range (xmin, xmax) (x -. xmax +. xmin ) - else - x - -let atom_id_of_int : int -> atom_id = - fun x -> ( assert (x>0) ; x) - -let distance_of_float : float -> distance = - fun x -> ( assert (x>=0.) ; Value x) - -let angle_of_float : float -> angle = - fun x -> Value (in_range (-180., 180.) x) - -let dihedral_of_float : float -> dihedral = - fun x -> Value (in_range (-360., 360.) x) - - -let atom_id_of_string : string -> atom_id = - fun i -> atom_id_of_int @@ int_of_string i - -let distance_of_string : string -> distance = - fun s -> - try - distance_of_float @@ float_of_string s - with _ -> Label s - -let angle_of_string : string -> angle = - fun s -> - try - angle_of_float @@ float_of_string s - with _ -> Label s - -let dihedral_of_string : string -> dihedral = - fun s -> - try - dihedral_of_float @@ float_of_string s - with _ -> Label s - -let int_of_atom_id : atom_id -> int = fun x -> x - -let float_of_distance : float StringMap.t -> distance -> float = - fun map -> function - | Value x -> x - | Label s -> StringMap.find s map - -let float_of_angle : float StringMap.t -> angle -> float = - fun map -> function - | Value x -> x - | Label s -> StringMap.find s map - -let float_of_dihedral : float StringMap.t -> dihedral -> float = - fun map -> function - | Value x -> x - | Label s -> StringMap.find s map - - -let string_of_line map = - let f_r = float_of_distance map - and f_a = float_of_angle map - and f_d = float_of_dihedral map - and i_i = int_of_atom_id - in function - | First e -> Printf.sprintf "%-3s" (Element.to_string e) - | Second (e, r) -> Printf.sprintf "%-3s %5d %f" (Element.to_string e) 1 (f_r r) - | Third (e, i, r, j, a) -> Printf.sprintf "%-3s %5d %f %5d %f" (Element.to_string e) (i_i i) (f_r r) (i_i j) (f_a a) - | Other (e, i, r, j, a, k, d) -> Printf.sprintf "%-3s %5d %f %5d %f %5d %f" (Element.to_string e) (i_i i) (f_r r) (i_i j) (f_a a) (i_i k) (f_d d) - | Coord (c, f) -> Printf.sprintf "%s %f" c f - - -let line_of_string l = - let line_clean = - Str.split (Str.regexp " ") l - |> List.filter (fun x -> x <> "") - in - match line_clean with - | e :: [] -> First (Element.of_string e) - | e :: _ :: r :: [] -> Second - (Element.of_string e, - distance_of_string r) - | e :: i :: r :: j :: a :: [] -> Third - (Element.of_string e, - atom_id_of_string i, - distance_of_string r, - atom_id_of_string j, - angle_of_string a) - | e :: i :: r :: j :: a :: k :: d :: [] -> Other - (Element.of_string e, - atom_id_of_string i, - distance_of_string r, - atom_id_of_string j, - angle_of_string a, - atom_id_of_string k, - dihedral_of_string d) - | c :: f :: [] -> Coord (c, float_of_string f) - | _ -> failwith ("Syntax error: "^l) - - -let of_string t = - let l = - Str.split (Str.regexp "\n") t - |> List.map String.trim - |> List.filter (fun x -> x <> "") - |> List.map line_of_string - in - - let l = - match l with - | First _ :: Second _ :: Third _ :: _ - | First _ :: Second _ :: Coord _ :: [] - | First _ :: Second _ :: [] - | First _ :: [] -> l - | _ -> failwith "Syntax error" - in - - let (l, m) = - let rec work lst map = function - | (First _ as x) :: rest - | (Second _ as x) :: rest - | (Third _ as x) :: rest - | (Other _ as x) :: rest -> work (x::lst) map rest - | (Coord (c,f)) :: rest -> work lst (StringMap.add c f map) rest - | [] -> (List.rev lst, map) - in - work [] (StringMap.empty) l - in - (Array.of_list l, m) - - -(** Linear algebra *) - -let (|-) (x,y,z) (x',y',z') = - ( x-.x', y-.y', z-.z' ) - -let (|+) (x,y,z) (x',y',z') = - ( x+.x', y+.y', z+.z' ) - -let (|.) s (x,y,z) = - ( s*.x, s*.y, s*.z ) - -let dot (x,y,z) (x',y',z') = - x*.x' +. y*.y' +. z*.z' - -let norm u = - sqrt @@ dot u u - -let normalized u = - 1. /. (norm u) |. u - -let cross (x,y,z) (x',y',z') = - ((y *. z' -. z *. y'), -. (x *. z' -. z *. x'), (x *. y' -. y *. x')) - -let rotation_matrix axis angle = - (* Euler-Rodrigues formula for rotation matrix, taken from - https://github.com/jevandezande/zmatrix/blob/master/converter.py - ,*) - let a = - (cos (angle *. to_radian *. 0.5)) - in - let (b, c, d) = - (-. sin (angle *. to_radian *. 0.5)) |. (normalized axis) - in - Array.of_list @@ - [(a *. a +. b *. b -. c *. c -. d *. d, - 2. *. (b *. c -. a *. d), - 2. *. (b *. d +. a *. c)); - (2. *. (b *. c +. a *. d), - a *. a +. c *. c -.b *. b -. d *. d, - 2. *. (c *. d -. a *. b)); - (2. *. (b *. d -. a *. c), - 2. *. (c *. d +. a *. b), - a *. a +. d *. d -. b *. b -. c *. c)] - - -let apply_rotation_matrix rot u = - (dot rot.(0) u, dot rot.(1) u, dot rot.(2) u) - - -let to_xyz (z,map) = - let result = - Array.make (Array.length z) None - in - - let get_cartesian_coord i = - match result.(i-1) with - | None -> failwith @@ Printf.sprintf "Atom %d is defined in the future" i - | Some (_, x, y, z) -> (x, y, z) - in - - - let append_line i' = - match z.(i') with - | First e -> - result.(i') <- Some (e, 0., 0., 0.) - | Second (e, r) -> - let r = - float_of_distance map r - in - result.(i') <- Some (e, 0., 0., r) - | Third (e, i, r, j, a) -> - begin - let i, r, j, a = - int_of_atom_id i, - float_of_distance map r, - int_of_atom_id j, - float_of_angle map a - in - let ui, uj = - get_cartesian_coord i, - get_cartesian_coord j - in - let u_ij = - (uj |- ui) - in - let rot = - rotation_matrix (0., 1., 0.) a - in - let new_vec = - apply_rotation_matrix rot ( r |. (normalized u_ij)) - in - let (x, y, z) = - new_vec |+ ui - in - result.(i') <- Some (e, x, y, z) - end - | Other (e, i, r, j, a, k, d) -> - begin - let i, r, j, a, k, d = - int_of_atom_id i, - float_of_distance map r, - int_of_atom_id j, - float_of_angle map a, - int_of_atom_id k, - float_of_dihedral map d - in - let ui, uj, uk = - get_cartesian_coord i, - get_cartesian_coord j, - get_cartesian_coord k - in - let u_ij, u_kj = - (uj |- ui) , (uj |- uk) - in - let normal = - cross u_ij u_kj - in - let new_vec = - r |. (normalized u_ij) - |> apply_rotation_matrix (rotation_matrix normal a) - |> apply_rotation_matrix (rotation_matrix u_ij d) - in - let (x, y, z) = - new_vec |+ ui - in - result.(i') <- Some (e, x, y, z) - end - | Coord _ -> () - in - Array.iteri (fun i _ -> append_line i) z; - let result = - Array.map (function - | Some x -> x - | None -> failwith "Some atoms were not defined" ) result - in - result - - -let to_xyz_string (l,map) = - String.concat "\n" - ( to_xyz (l,map) - |> Array.map (fun (e,x,y,z) -> - Printf.sprintf "%s %f %f %f" (Element.to_string e) x y z) - |> Array.to_list - ) - #+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 (a, map) = - let f = string_of_line map in - Format.fprintf ppf "@["; - Array.iter (fun line -> - Format.fprintf ppf "%s@." (f line) - ) a; - Format.fprintf ppf "@]" - #+end_src - diff --git a/top/lib/install_printers.ml b/top/lib/install_printers.ml index 8d778c1..ff95c4a 100644 --- a/top/lib/install_printers.ml +++ b/top/lib/install_printers.ml @@ -1,15 +1,13 @@ -(** Intall printers printers:3]] *) +(* [[file:~/QCaml/top/install_printers.org::*Intall printers][Intall printers:3]] *) let printers = [ - "Mo.Frozen_core.pp" ; "Mo.Localization.pp" ; - "Particles.Electrons.pp" ; "Particles.Element.pp" ; "Particles.Nuclei.pp" ; "Particles.Zmatrix.pp" ; "Perturbation.Mp2.pp" ; "Simulation.pp" ; - + ] let eval_exn str = @@ -27,3 +25,4 @@ let rec install_printers = function let () = if not (install_printers printers) then Format.eprintf "Problem installing QCaml-printers@." +(* Intall printers:3 ends here *)