From e909bb9abbd218071db9a80d0e70b4958b4b664d Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 2 Jan 2018 22:33:17 +0100 Subject: [PATCH] Added many files for nuclei --- Nuclei/Charge.ml | 6 +- Nuclei/{Mass.mli => Mass.ml} | 0 Nuclei/Nuclei.ml | 2 +- Nuclei/Nuclei_parser.mly | 30 +- Nuclei/{Radius.mli => Radius.ml} | 0 Nuclei/Zmatrix.ml | 313 ++++++++++++++++++ Positive_float.ml => Utils/Positive_float.ml | 2 +- .../Positive_float.mli | 0 Utils/Units.ml | 10 + Utils/Units.mli | 7 + _tags | 2 +- 11 files changed, 359 insertions(+), 13 deletions(-) rename Nuclei/{Mass.mli => Mass.ml} (100%) rename Nuclei/{Radius.mli => Radius.ml} (100%) create mode 100644 Nuclei/Zmatrix.ml rename Positive_float.ml => Utils/Positive_float.ml (91%) rename Positive_float.mli => Utils/Positive_float.mli (100%) create mode 100644 Utils/Units.ml create mode 100644 Utils/Units.mli diff --git a/Nuclei/Charge.ml b/Nuclei/Charge.ml index 8896bef..d4f1738 100644 --- a/Nuclei/Charge.ml +++ b/Nuclei/Charge.ml @@ -1,11 +1,11 @@ type t = float let of_float x = x -let of_int i = Float.of_int i -let of_string s = Float.of_string s +let of_int i = float_of_int i +let of_string s = float_of_string s let to_float x = x -let to_int x = Float.to_int x +let to_int x = int_of_float x let to_string x = if x >= 0. then Printf.sprintf "+%f" x diff --git a/Nuclei/Mass.mli b/Nuclei/Mass.ml similarity index 100% rename from Nuclei/Mass.mli rename to Nuclei/Mass.ml diff --git a/Nuclei/Nuclei.ml b/Nuclei/Nuclei.ml index 9684842..cf030de 100644 --- a/Nuclei/Nuclei.ml +++ b/Nuclei/Nuclei.ml @@ -1 +1 @@ -type xyz_input = string * ( (string * (float array)) array) +type xyz_input = string * ( (Element.t * (float array)) array) diff --git a/Nuclei/Nuclei_parser.mly b/Nuclei/Nuclei_parser.mly index 7c7d888..b43857e 100644 --- a/Nuclei/Nuclei_parser.mly +++ b/Nuclei/Nuclei_parser.mly @@ -11,13 +11,13 @@ exception InputError of string %token FLOAT %token EOF -%start input -%type input +%start input_xyz +%type input_xyz %% /* Grammar rules and actions follow */ -input: - | integer title atoms { +input_xyz: + | integer title atoms_xyz { let len = List.length $3 in if len <> $1 then let error_msg = Printf.sprintf "%d atoms entered, expected %d" len $1 in @@ -43,13 +43,29 @@ title_list: | { "" } | title_list text { $1 ^ $2 } -atoms: +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 { ($2, [| $4;$6;$8 |]) :: $1 } - | atoms_list WORD SPACE FLOAT SPACE FLOAT SPACE FLOAT SPACE EOL { ($2, [| $4;$6;$8 |]) :: $1 } + | atoms_list WORD SPACE FLOAT SPACE FLOAT SPACE FLOAT EOL { (Element.of_string $2, [| $4;$6;$8 |]) :: $1 } + | atoms_list WORD SPACE FLOAT SPACE FLOAT SPACE FLOAT SPACE EOL { (Element.of_string $2, [| $4;$6;$8 |]) :: $1 } + +label: + | FLOAT { Zmatrix.Value $1 } + | WORD { Zmatrix.Label $1 } + +first_line: + | WORD { $1 } + +second_line: + | WORD INTEGER label { ($1,$2,$3) } + +third_line: + | WORD INTEGER label INTEGER label { ($1,$2,$3,$4,$5) } + +nth_line: + | WORD INTEGER label INTEGER label INTEGER label { ($1,$2,$3,$4,$5,$6,$7) } diff --git a/Nuclei/Radius.mli b/Nuclei/Radius.ml similarity index 100% rename from Nuclei/Radius.mli rename to Nuclei/Radius.ml diff --git a/Nuclei/Zmatrix.ml b/Nuclei/Zmatrix.ml new file mode 100644 index 0000000..3eb4577 --- /dev/null +++ b/Nuclei/Zmatrix.ml @@ -0,0 +1,313 @@ +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 + +let pi = acos (-1.) +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 + + +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) + + +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 :: i :: 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) + + +type t = (line array * float StringMap.t) + +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 center_of_mass l = +let (x,y,z) = + let sum_mass, com = + Array.fold_left (fun (s,com) (e,x,y,z) -> + let mass = + Positive_float.to_float @@ Element.mass e + in + (s +. mass, ( mass |. (x,y,z) ) |+ com) ) + (0., (0.,0.,0.)) l + in + (1. /. sum_mass) |. com +in +Printf.printf "%f %f %f\n" x y z ; (x,y,z) + +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 + Array.to_list result + + +let to_xyz_string (l,map) = + String.concat "\n" + ( to_xyz (l,map) + |> List.map (fun (e,x,y,z) -> + Printf.sprintf "%s %f %f %f\n" (Element.to_string e) x y z) ) + + + diff --git a/Positive_float.ml b/Utils/Positive_float.ml similarity index 91% rename from Positive_float.ml rename to Utils/Positive_float.ml index 134c224..6f3faea 100644 --- a/Positive_float.ml +++ b/Utils/Positive_float.ml @@ -1,4 +1,4 @@ -type t +type t = float let of_float x = assert ( x >= 0. ); diff --git a/Positive_float.mli b/Utils/Positive_float.mli similarity index 100% rename from Positive_float.mli rename to Utils/Positive_float.mli diff --git a/Utils/Units.ml b/Utils/Units.ml new file mode 100644 index 0000000..ae8b21c --- /dev/null +++ b/Utils/Units.ml @@ -0,0 +1,10 @@ +type units = +| Bohr +| Angstrom +;; + +let angstrom_to_bohr = 1. /. 0.52917721092 +let bohr_to_angstrom = 0.52917721092 +;; + + diff --git a/Utils/Units.mli b/Utils/Units.mli new file mode 100644 index 0000000..eff3fb2 --- /dev/null +++ b/Utils/Units.mli @@ -0,0 +1,7 @@ +type units = Bohr | Angstrom + +(** Conversion functions *) +val angstrom_to_bohr : float +val bohr_to_angstrom : float + + diff --git a/_tags b/_tags index b8f7c84..633df42 100644 --- a/_tags +++ b/_tags @@ -1 +1 @@ -true: package(re) +true: package(str)