LCPQ/QCaml
10
1
Fork 0
mirror of https://gitlab.com/scemama/QCaml.git synced 2024-11-18 12:03:40 +01:00
Code Releases Activity

Removing org-mode

Browse Source
This commit is contained in:
Anthony Scemama 2024-02-28 10:34:39 +01:00
parent d6ef25d55a
commit feebc5a0e8
21 changed files with 308 additions and 1727 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

86
mo/lib/localization.ml
View File

@ -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 *)

33
mo/lib/localization.mli
View File

@ -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 *)

3
mo/test/localization.ml
View File

@ -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 *)

160
particles/electrons.org
View File

@ -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

328
particles/element.org
View File

@ -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
<<types>>
#+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

30
particles/lib/electrons.ml
View File

@ -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 *)

40
particles/lib/electrons.mli
View File

@ -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 *)

63
particles/lib/element.ml
View File

@ -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 *)

69
particles/lib/element.mli
View File

@ -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 *)

4
particles/lib/mass.ml
View File

@ -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 *)

11
particles/lib/mass.mli
View File

@ -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 *)

40
particles/lib/nuclei.ml
View File

@ -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 *)

56
particles/lib/nuclei.mli
View File

@ -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 *)

8
particles/lib/nuclei_lexer.mll
View File

@ -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 *)

9
particles/lib/xyz_ast.mli
View File

@ -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 *)

129
particles/lib/zmatrix.ml
View File

@ -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 *)

64
particles/lib/zmatrix.mli
View File

@ -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 *)

27
particles/mass.org
View File

@ -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

471
particles/nuclei.org
View File

@ -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 <string> SPACE
%token <string> WORD
%token <int> INTEGER
%token <float> FLOAT
%token EOF
%start input
%type <Xyz_ast.xyz_file> 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

397
particles/zmatrix.org
View File

@ -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

7
top/lib/install_printers.ml
View File

@ -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 *)