diff --git a/particles/lib/mass.ml b/particles/lib/mass.ml index 3fec2d3..f7fc492 100644 --- a/particles/lib/mass.ml +++ b/particles/lib/mass.ml @@ -1,3 +1,3 @@ -(** Atomic mass. *) - +(* [[file:../mass.org::*Atomic mass][Atomic mass:2]] *) include Common.Non_negative_float +(* Atomic mass:2 ends here *) diff --git a/particles/lib/mass.mli b/particles/lib/mass.mli index cd230dc..9d91dd3 100644 --- a/particles/lib/mass.mli +++ b/particles/lib/mass.mli @@ -1,3 +1,12 @@ -(** Atomic mass. *) +(* Atomic mass + * :PROPERTIES: + * :header-args: :noweb yes :comments both + * :END: + * + * Atomic mass, a non-negative float. + * + * #+NAME: types *) +(* [[file:../mass.org::types][types]] *) include module type of Common.Non_negative_float +(* types ends here *) diff --git a/particles/lib/zmatrix.ml b/particles/lib/zmatrix.ml index 27431f3..d404a0b 100644 --- a/particles/lib/zmatrix.ml +++ b/particles/lib/zmatrix.ml @@ -1,3 +1,4 @@ +(* [[file:../zmatrix.org::*Type][Type:2]] *) module StringMap = Map.Make(String) type atom_id = int @@ -5,16 +6,34 @@ type angle = Label of string | Value of float type distance = Label of string | Value of float type dihedral = Label of string | Value of float -let pi = acos (-1.) +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) +(* Type:2 ends here *) + + + +(* | ~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 | *) + + +(* [[file:../zmatrix.org::*Conversion][Conversion:2]] *) +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 + 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) @@ -34,47 +53,38 @@ let atom_id_of_string : string -> atom_id = let distance_of_string : string -> distance = fun s -> - try - distance_of_float @@ float_of_string s - with _ -> Label 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 + 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 - + 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 + | 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 + | 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 - - -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) + | Value x -> x + | Label s -> StringMap.find s map let string_of_line map = @@ -83,11 +93,11 @@ let string_of_line 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 + | 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 = @@ -98,28 +108,26 @@ let line_of_string l = match line_clean with | e :: [] -> First (Element.of_string e) | e :: _ :: r :: [] -> Second - (Element.of_string e, - distance_of_string r) + (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) + (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) + (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) -type t = (line array * float StringMap.t) - let of_string t = let l = Str.split (Str.regexp "\n") t @@ -196,10 +204,10 @@ let rotation_matrix axis angle = 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 @@ -296,6 +304,14 @@ let to_xyz_string (l,map) = Printf.sprintf "%s %f %f %f\n" (Element.to_string e) x y z) |> Array.to_list ) +(* Conversion:2 ends here *) - - +(* [[file:../zmatrix.org::*Printers][Printers:2]] *) +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 "@]" +(* Printers:2 ends here *) diff --git a/particles/lib/zmatrix.mli b/particles/lib/zmatrix.mli new file mode 100644 index 0000000..f964ea4 --- /dev/null +++ b/particles/lib/zmatrix.mli @@ -0,0 +1,23 @@ +(* Type + * + * #+NAME: types *) + +(* [[file:../zmatrix.org::types][types]] *) +type t +(* types ends here *) + +(* Conversion *) + + +(* [[file:../zmatrix.org::*Conversion][Conversion:1]] *) +val of_string : string -> t +val to_xyz : t -> (Element.t * float * float * float) array +val to_xyz_string : t -> string +(* Conversion:1 ends here *) + +(* Printers *) + + +(* [[file:../zmatrix.org::*Printers][Printers:1]] *) +val pp : Format.formatter -> t -> unit +(* Printers:1 ends here *) diff --git a/particles/mass.org b/particles/mass.org new file mode 100644 index 0000000..b727702 --- /dev/null +++ b/particles/mass.org @@ -0,0 +1,27 @@ +#+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/zmatrix.org b/particles/zmatrix.org new file mode 100644 index 0000000..a54f98b --- /dev/null +++ b/particles/zmatrix.org @@ -0,0 +1,354 @@ +#+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_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\n" (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/install_printers.org b/top/install_printers.org index 7b9b2af..002ba96 100644 --- a/top/install_printers.org +++ b/top/install_printers.org @@ -46,13 +46,12 @@ let eval_exn str = let rec install_printers = function - | [] -> true + | [] -> eval_exn "#require \"lacaml.top\";;" | printer :: printers -> let cmd = Printf.sprintf "#install_printer %s;;" printer in eval_exn cmd && install_printers printers let () = - eval_exn "#require \"lacaml.top\";;"; if not (install_printers printers) then Format.eprintf "Problem installing QCaml-printers@." diff --git a/top/lib/install_printers.ml b/top/lib/install_printers.ml index 7edd6bb..ae25f4f 100644 --- a/top/lib/install_printers.ml +++ b/top/lib/install_printers.ml @@ -11,6 +11,7 @@ let printers = "Common.Zkey.pp" ; "Particles.Electrons.pp" ; "Particles.Element.pp" ; + "Particles.Zmatrix.pp" ; ] @@ -21,13 +22,12 @@ let eval_exn str = let rec install_printers = function - | [] -> true + | [] -> eval_exn "#require \"lacaml.top\";;" | printer :: printers -> let cmd = Printf.sprintf "#install_printer %s;;" printer in eval_exn cmd && install_printers printers let () = - eval_exn "#require \"lacaml.top\";;"; if not (install_printers printers) then Format.eprintf "Problem installing QCaml-printers@." (* Intall printers:2 ends here *)