From 21a9b30d2dc975225f5a9fa52f2c95bd88f2f6b4 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 8 Sep 2016 22:37:05 +0200 Subject: [PATCH] Added Zmatrix module --- ocaml/Element.mli | 1 + ocaml/Molecule.ml | 22 +- ocaml/Molecule.mli | 12 + ocaml/Zmatrix.ml | 326 +++++++++++++++++++++ ocaml/_tags | 2 +- ocaml/qp_create_ezfio_from_xyz.ml | 14 +- src/Determinants/slater_rules.irp.f | 16 +- src/Integrals_Bielec/ao_bi_integrals.irp.f | 2 +- 8 files changed, 379 insertions(+), 16 deletions(-) create mode 100644 ocaml/Zmatrix.ml diff --git a/ocaml/Element.mli b/ocaml/Element.mli index 5edfdf31..2c899b3b 100644 --- a/ocaml/Element.mli +++ b/ocaml/Element.mli @@ -19,3 +19,4 @@ val to_charge : t -> Charge.t val of_charge : Charge.t -> t val covalent_radius : t -> Qptypes.Positive_float.t val vdw_radius : t -> Qptypes.Positive_float.t +val mass : t -> Qptypes.Positive_float.t diff --git a/ocaml/Molecule.ml b/ocaml/Molecule.ml index a9d73432..a26e23b5 100644 --- a/ocaml/Molecule.ml +++ b/ocaml/Molecule.ml @@ -147,10 +147,28 @@ let of_xyz_file let (_,buffer) = In_channel.read_all filename |> String.lsplit2_exn ~on:'\n' in let (_,buffer) = String.lsplit2_exn buffer ~on:'\n' in - of_xyz_string ~charge:charge ~multiplicity:multiplicity - ~units:units buffer + of_xyz_string ~charge ~multiplicity ~units buffer +let of_zmt_file + ?(charge=(Charge.of_int 0)) ?(multiplicity=(Multiplicity.of_int 1)) + ?(units=Units.Angstrom) + filename = + In_channel.read_all filename + |> Zmatrix.of_string + |> Zmatrix.to_xyz_string + |> of_xyz_string ~charge ~multiplicity ~units + + +let of_file + ?(charge=(Charge.of_int 0)) ?(multiplicity=(Multiplicity.of_int 1)) + ?(units=Units.Angstrom) + filename = + try + of_xyz_file ~charge ~multiplicity ~units filename + with _ -> + of_zmt_file ~charge ~multiplicity ~units filename + let distance_matrix molecule = let coord = diff --git a/ocaml/Molecule.mli b/ocaml/Molecule.mli index f81f28a3..f6201b18 100644 --- a/ocaml/Molecule.mli +++ b/ocaml/Molecule.mli @@ -29,6 +29,18 @@ val of_xyz_file : ?multiplicity:Multiplicity.t -> ?units:Units.units -> string -> t +(** Creates a molecule from a zmt file *) +val of_zmt_file : + ?charge:Charge.t -> + ?multiplicity:Multiplicity.t -> + ?units:Units.units -> string -> t + +(** Creates a molecule from a file (xyz or zmt) *) +val of_file : + ?charge:Charge.t -> + ?multiplicity:Multiplicity.t -> + ?units:Units.units -> string -> t + (** Creates a molecule from an xyz file in a string *) val of_xyz_string : ?charge:Charge.t -> diff --git a/ocaml/Zmatrix.ml b/ocaml/Zmatrix.ml new file mode 100644 index 00000000..294b1112 --- /dev/null +++ b/ocaml/Zmatrix.ml @@ -0,0 +1,326 @@ +open Qptypes + +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 _ :: [] + | 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)] +(* + [(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 + 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 ) + + + diff --git a/ocaml/_tags b/ocaml/_tags index 3f5cd9b6..0935c0bb 100644 --- a/ocaml/_tags +++ b/ocaml/_tags @@ -1,3 +1,3 @@ -true: package(core,cryptokit,ZMQ,sexplib.syntax) +true: package(core,cryptokit,ZMQ,sexplib.syntax,str) true: thread false: profile diff --git a/ocaml/qp_create_ezfio_from_xyz.ml b/ocaml/qp_create_ezfio_from_xyz.ml index 710523e4..c79bf550 100644 --- a/ocaml/qp_create_ezfio_from_xyz.ml +++ b/ocaml/qp_create_ezfio_from_xyz.ml @@ -19,7 +19,7 @@ let spec = ~doc:"string Name of the pseudopotential" +> flag "cart" no_arg ~doc:" Compute AOs in the Cartesian basis set (6d, 10f, ...)" - +> anon ("xyz_file" %: file ) + +> anon ("(xyz_file|zmt_file)" %: file ) (** Handle dummy atoms placed on bonds *) @@ -93,7 +93,7 @@ let run ?o b c d m p cart xyz_file = (* Read molecule *) let molecule = - (Molecule.of_xyz_file xyz_file ~charge:(Charge.of_int c) + (Molecule.of_file xyz_file ~charge:(Charge.of_int c) ~multiplicity:(Multiplicity.of_int m) ) in let dummy = @@ -309,7 +309,8 @@ let run ?o b c d m p cart xyz_file = | None -> begin match String.rsplit2 ~on:'.' xyz_file with - | Some (x,"xyz") -> x^".ezfio" + | Some (x,"xyz") + | Some (x,"zmt") -> x^".ezfio" | _ -> xyz_file^".ezfio" end in @@ -640,9 +641,10 @@ let command = ============================ -Creates an EZFIO directory from a standard xyz file. The basis set is defined -as a single string if all the atoms are taken from the same basis set, -otherwise specific elements can be defined as follows: +Creates an EZFIO directory from a standard xyz file or from a z-matrix file +in Gaussian format. The basis set is defined as a single string if all the +atoms are taken from the same basis set, otherwise specific elements can be +defined as follows: -b \"cc-pcvdz | H:cc-pvdz | C:6-31g\" diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f index eda2d9b4..21cacaad 100644 --- a/src/Determinants/slater_rules.irp.f +++ b/src/Determinants/slater_rules.irp.f @@ -1730,15 +1730,19 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint) do j=shortcut(sh2,1),endi org_j = sort_idx(j,1) - ext = exa - do ni=1,Nint + ext = exa + popcnt(xor(sorted_i(1), sorted(1,j,1))) + if(ext > 4) then + cycle + endif + do ni=2,Nint ext = ext + popcnt(xor(sorted_i(ni), sorted(ni,j,1))) end do - if(ext <= 4) then - call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij) - vt (org_i) = vt (org_i) + hij*u_0(org_j) - vt (org_j) = vt (org_j) + hij*u_0(org_i) + if(ext > 4) then + cycle endif + call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij) + vt (org_i) = vt (org_i) + hij*u_0(org_j) + vt (org_j) = vt (org_j) + hij*u_0(org_i) enddo enddo enddo diff --git a/src/Integrals_Bielec/ao_bi_integrals.irp.f b/src/Integrals_Bielec/ao_bi_integrals.irp.f index 1dcea81f..2ebb402e 100644 --- a/src/Integrals_Bielec/ao_bi_integrals.irp.f +++ b/src/Integrals_Bielec/ao_bi_integrals.irp.f @@ -1211,7 +1211,7 @@ subroutine compute_ao_integrals_jl(j,l,n_integrals,buffer_i,buffer_value) cycle endif !DIR$ FORCEINLINE - integral = ao_bielec_integral(i,k,j,l) + integral = ao_bielec_integral(i,k,j,l) ! i,k : r1 j,l : r2 if (abs(integral) < thr) then cycle endif