10
1
mirror of https://gitlab.com/scemama/QCaml.git synced 2024-12-23 04:43:32 +01:00
QCaml/particles/lib/zmatrix.ml

314 lines
8.2 KiB
OCaml
Raw Normal View History

2024-02-28 10:34:39 +01:00
(** Type *)
2018-01-02 22:33:17 +01:00
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
2024-02-28 10:34:39 +01:00
type line =
2020-12-29 02:29:43 +01:00
| 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)
2024-02-28 10:34:39 +01:00
(** Conversion *)
2020-12-29 02:29:43 +01:00
let pi = Common.Constants.pi
2018-01-02 22:33:17 +01:00
let to_radian = pi /. 180.
let rec in_range (xmin, xmax) x =
2020-12-29 02:29:43 +01:00
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
2018-01-02 22:33:17 +01:00
2024-02-28 10:34:39 +01:00
let atom_id_of_int : int -> atom_id =
2018-01-02 22:33:17 +01:00
fun x -> ( assert (x>0) ; x)
2024-02-28 10:34:39 +01:00
let distance_of_float : float -> distance =
2018-01-02 22:33:17 +01:00
fun x -> ( assert (x>=0.) ; Value x)
2024-02-28 10:34:39 +01:00
let angle_of_float : float -> angle =
2018-01-02 22:33:17 +01:00
fun x -> Value (in_range (-180., 180.) x)
2024-02-28 10:34:39 +01:00
let dihedral_of_float : float -> dihedral =
2018-01-02 22:33:17 +01:00
fun x -> Value (in_range (-360., 360.) x)
2024-02-28 10:34:39 +01:00
let atom_id_of_string : string -> atom_id =
2018-01-02 22:33:17 +01:00
fun i -> atom_id_of_int @@ int_of_string i
let distance_of_string : string -> distance =
2024-02-28 10:34:39 +01:00
fun s ->
2020-12-29 02:29:43 +01:00
try
distance_of_float @@ float_of_string s
with _ -> Label s
2018-01-02 22:33:17 +01:00
let angle_of_string : string -> angle =
2024-02-28 10:34:39 +01:00
fun s ->
2020-12-29 02:29:43 +01:00
try
angle_of_float @@ float_of_string s
with _ -> Label s
2018-01-02 22:33:17 +01:00
let dihedral_of_string : string -> dihedral =
2024-02-28 10:34:39 +01:00
fun s ->
2020-12-29 02:29:43 +01:00
try
dihedral_of_float @@ float_of_string s
with _ -> Label s
2018-01-02 22:33:17 +01:00
let int_of_atom_id : atom_id -> int = fun x -> x
2024-02-28 10:34:39 +01:00
let float_of_distance : float StringMap.t -> distance -> float =
2018-01-02 22:33:17 +01:00
fun map -> function
2020-12-29 02:29:43 +01:00
| Value x -> x
| Label s -> StringMap.find s map
2018-01-02 22:33:17 +01:00
let float_of_angle : float StringMap.t -> angle -> float =
fun map -> function
2020-12-29 02:29:43 +01:00
| Value x -> x
| Label s -> StringMap.find s map
2018-01-02 22:33:17 +01:00
let float_of_dihedral : float StringMap.t -> dihedral -> float =
fun map -> function
2020-12-29 02:29:43 +01:00
| Value x -> x
| Label s -> StringMap.find s map
2018-01-02 22:33:17 +01:00
2024-02-28 10:34:39 +01:00
let string_of_line map =
2018-01-02 22:33:17 +01:00
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
2020-12-29 02:29:43 +01:00
| 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
2018-01-02 22:33:17 +01:00
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)
2020-09-26 12:02:53 +02:00
| e :: _ :: r :: [] -> Second
2020-12-29 02:29:43 +01:00
(Element.of_string e,
distance_of_string r)
2024-02-28 10:34:39 +01:00
| e :: i :: r :: j :: a :: [] -> Third
2020-12-29 02:29:43 +01:00
(Element.of_string e,
atom_id_of_string i,
distance_of_string r,
atom_id_of_string j,
angle_of_string a)
2018-01-02 22:33:17 +01:00
| e :: i :: r :: j :: a :: k :: d :: [] -> Other
2020-12-29 02:29:43 +01:00
(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)
2018-01-02 22:33:17 +01:00
| 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
2024-02-28 10:34:39 +01:00
let l =
2018-01-02 22:33:17 +01:00
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)
2024-02-28 10:34:39 +01:00
2018-01-02 22:33:17 +01:00
(** 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'))
2024-02-28 10:34:39 +01:00
let rotation_matrix axis angle =
2018-01-02 22:33:17 +01:00
(* Euler-Rodrigues formula for rotation matrix, taken from
https://github.com/jevandezande/zmatrix/blob/master/converter.py
*)
2024-02-28 10:34:39 +01:00
let a =
2018-01-02 22:33:17 +01:00
(cos (angle *. to_radian *. 0.5))
in
2024-02-28 10:34:39 +01:00
let (b, c, d) =
2018-01-02 22:33:17 +01:00
(-. sin (angle *. to_radian *. 0.5)) |. (normalized axis)
in
2024-02-28 10:34:39 +01:00
Array.of_list @@
2018-01-02 22:33:17 +01:00
[(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)]
2024-02-28 10:34:39 +01:00
2018-01-02 22:33:17 +01:00
let apply_rotation_matrix rot u =
(dot rot.(0) u, dot rot.(1) u, dot rot.(2) u)
2024-02-28 10:34:39 +01:00
2020-12-29 02:29:43 +01:00
2018-01-02 22:33:17 +01:00
let to_xyz (z,map) =
let result =
Array.make (Array.length z) None
2024-02-28 10:34:39 +01:00
in
2018-01-02 22:33:17 +01:00
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
2024-02-28 10:34:39 +01:00
| First e ->
2018-01-02 22:33:17 +01:00
result.(i') <- Some (e, 0., 0., 0.)
2024-02-28 10:34:39 +01:00
| Second (e, r) ->
2018-01-02 22:33:17 +01:00
let r =
float_of_distance map r
in
result.(i') <- Some (e, 0., 0., r)
2024-02-28 10:34:39 +01:00
| Third (e, i, r, j, a) ->
2018-01-02 22:33:17 +01:00
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 =
2024-02-28 10:34:39 +01:00
get_cartesian_coord i,
2018-01-02 22:33:17 +01:00
get_cartesian_coord j
in
2024-02-28 10:34:39 +01:00
let u_ij =
2018-01-02 22:33:17 +01:00
(uj |- ui)
in
2024-02-28 10:34:39 +01:00
let rot =
2018-01-02 22:33:17 +01:00
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
2024-02-28 10:34:39 +01:00
| Other (e, i, r, j, a, k, d) ->
2018-01-02 22:33:17 +01:00
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 =
2024-02-28 10:34:39 +01:00
get_cartesian_coord i,
2018-01-02 22:33:17 +01:00
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
2024-02-28 10:34:39 +01:00
let new_vec =
2018-01-02 22:33:17 +01:00
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;
2024-02-28 10:34:39 +01:00
let result =
2018-01-02 22:33:17 +01:00
Array.map (function
| Some x -> x
| None -> failwith "Some atoms were not defined" ) result
in
2018-01-17 15:56:57 +01:00
result
2018-01-02 22:33:17 +01:00
let to_xyz_string (l,map) =
2024-02-28 10:34:39 +01:00
String.concat "\n"
( to_xyz (l,map)
|> Array.map (fun (e,x,y,z) ->
2020-12-29 02:38:12 +01:00
Printf.sprintf "%s %f %f %f" (Element.to_string e) x y z)
2018-01-17 15:56:57 +01:00
|> Array.to_list
)
2020-12-29 02:29:43 +01:00
2024-02-28 10:34:39 +01:00
(** Printers *)
2020-12-29 02:29:43 +01:00
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 "@]"
2024-02-28 10:34:39 +01:00