From d56c9dcb0f09099febb8c63b37ca3448689fd5a2 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 17 Nov 2015 22:25:26 +0100 Subject: [PATCH] Option to add automatically dummy basis functions --- ocaml/Atom.ml | 12 +- ocaml/Atom.mli | 4 +- ocaml/Element.ml | 218 ++++++++++++++++++----- ocaml/Element.mli | 4 +- ocaml/Input_nuclei.ml | 4 +- ocaml/Molecule.ml | 96 +++++++--- ocaml/Molecule.mli | 4 + ocaml/Point3d.ml | 25 ++- ocaml/Point3d.mli | 7 +- ocaml/qp_create_ezfio_from_xyz.ml | 97 ++++++++-- ocaml/test_molecule.ml | 9 + scripts/get_basis.sh | 1 - scripts/pseudo/elts_num_ele.py | 3 +- scripts/pseudo/put_pseudo_in_ezfio.py | 32 +++- src/Integrals_Bielec/map_integrals.irp.f | 9 +- 15 files changed, 408 insertions(+), 117 deletions(-) diff --git a/ocaml/Atom.ml b/ocaml/Atom.ml index 2aafe644..832cfa5b 100644 --- a/ocaml/Atom.ml +++ b/ocaml/Atom.ml @@ -8,8 +8,8 @@ type t = coord : Point3d.t ; } with sexp -(** Read xyz coordinates of the atom with unit u *) -let of_string u s = +(** Read xyz coordinates of the atom *) +let of_string ~units s = let buffer = s |> String.split ~on:' ' |> List.filter ~f:(fun x -> x <> "") @@ -18,21 +18,21 @@ let of_string u s = | [ name; charge; x; y; z ] -> { element = Element.of_string name ; charge = Charge.of_string charge ; - coord = Point3d.of_string u (String.concat [x; y; z] ~sep:" ") + coord = Point3d.of_string ~units (String.concat [x; y; z] ~sep:" ") } | [ name; x; y; z ] -> let e = Element.of_string name in { element = e ; charge = Element.to_charge e; - coord = Point3d.of_string u (String.concat [x; y; z] ~sep:" ") + coord = Point3d.of_string ~units (String.concat [x; y; z] ~sep:" ") } | _ -> raise (AtomError s) ;; -let to_string u a = +let to_string ~units a = [ Element.to_string a.element ; Charge.to_string a.charge ; - Point3d.to_string u a.coord ] + Point3d.to_string ~units a.coord ] |> String.concat ~sep:" " ;; diff --git a/ocaml/Atom.mli b/ocaml/Atom.mli index 93f7d885..28915993 100644 --- a/ocaml/Atom.mli +++ b/ocaml/Atom.mli @@ -5,5 +5,5 @@ type t = { element : Element.t; charge : Charge.t; coord : Point3d.t; } val t_of_sexp : Sexplib.Sexp.t -> t val sexp_of_t : t -> Sexplib.Sexp.t -val of_string : Units.units -> string -> t -val to_string : Units.units -> t -> string +val of_string : units:Units.units -> string -> t +val to_string : units:Units.units -> t -> string diff --git a/ocaml/Element.ml b/ocaml/Element.ml index d282340b..6bc2de4e 100644 --- a/ocaml/Element.ml +++ b/ocaml/Element.ml @@ -1,4 +1,5 @@ -open Core.Std;; +open Core.Std +open Qptypes exception ElementError of string @@ -8,49 +9,49 @@ type t = |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 -with sexp;; +with sexp let of_string x = match (String.capitalize (String.lowercase 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 +| "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 | x -> raise (ElementError ("Element "^x^" unknown")) -;; + let to_string = function | X -> "X" @@ -90,7 +91,7 @@ let to_string = function | Se -> "Se" | Br -> "Br" | Kr -> "Kr" -;; + let to_long_string = function | X -> "Dummy" @@ -130,7 +131,7 @@ let to_long_string = function | Se -> "Selenium" | Br -> "Bromine" | Kr -> "Krypton" -;; + let to_charge c = let result = match c with @@ -172,7 +173,7 @@ let to_charge c = | Br -> 35 | Kr -> 36 in Charge.of_int result -;; + let of_charge c = match (Charge.to_int c) with | 0 -> X @@ -213,5 +214,134 @@ let of_charge c = match (Charge.to_int c) with | 35 -> Br | 36 -> Kr | x -> raise (ElementError ("Element of charge "^(string_of_int x)^" unknown")) -;; + + +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 + in + Units.angstrom_to_bohr *. (result x) + |> Positive_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 + in + Units.angstrom_to_bohr *. (result x) + |> Positive_float.of_float + +let mass x = + let result = function + | 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 + in + result x + |> Positive_float.of_float diff --git a/ocaml/Element.mli b/ocaml/Element.mli index 48bd3c61..8d9862c9 100644 --- a/ocaml/Element.mli +++ b/ocaml/Element.mli @@ -13,6 +13,8 @@ val of_string : string -> t val to_string : t -> string val to_long_string : t -> string -(** get the positive charge *) +(** Properties *) 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 diff --git a/ocaml/Input_nuclei.ml b/ocaml/Input_nuclei.ml index cbcb9f46..d050ded9 100644 --- a/ocaml/Input_nuclei.ml +++ b/ocaml/Input_nuclei.ml @@ -147,7 +147,7 @@ nucl_coord = %s (b.nucl_charge |> Array.to_list |> List.map ~f:(Charge.to_string) |> String.concat ~sep:", " ) (b.nucl_coord |> Array.to_list |> List.map - ~f:(Point3d.to_string Units.Bohr) |> String.concat ~sep:"\n" ) + ~f:(Point3d.to_string ~units:Units.Bohr) |> String.concat ~sep:"\n" ) ;; @@ -161,7 +161,7 @@ nucl_coord = %s Printf.sprintf " %-3s %d %s" (b.nucl_label.(i) |> Element.to_string) (b.nucl_charge.(i) |> Charge.to_int ) - (b.nucl_coord.(i) |> Point3d.to_string Units.Angstrom) ) + (b.nucl_coord.(i) |> Point3d.to_string ~units:Units.Angstrom) ) ) |> String.concat ~sep:"\n" in Printf.sprintf " diff --git a/ocaml/Molecule.ml b/ocaml/Molecule.ml index f295420b..f0800f7f 100644 --- a/ocaml/Molecule.ml +++ b/ocaml/Molecule.ml @@ -10,27 +10,36 @@ type t = { } with sexp let get_charge { nuclei ; elec_alpha ; elec_beta } = - let result = (Elec_alpha_number.to_int elec_alpha) + - (Elec_beta_number.to_int elec_beta) in + let result = + (Elec_alpha_number.to_int elec_alpha) + + (Elec_beta_number.to_int elec_beta) + in let rec nucl_charge = function | a::rest -> (Charge.to_float a.Atom.charge) +. nucl_charge rest | [] -> 0. in Charge.of_float (nucl_charge nuclei -. (Float.of_int result)) -;; + let get_multiplicity m = - let elec_alpha = m.elec_alpha in + let elec_alpha = + m.elec_alpha + in Multiplicity.of_alpha_beta elec_alpha m.elec_beta -;; + let get_nucl_num m = - let nmax = (List.length m.nuclei) in + let nmax = + List.length m.nuclei + in Nucl_number.of_int nmax ~max:nmax -;; + let name m = - let cm = Charge.to_int (get_charge m) in + let cm = + get_charge m + |> Charge.to_int + in let c = match cm with | 0 -> "" @@ -39,8 +48,12 @@ let name m = | i when i>1 -> Printf.sprintf " (%d+)" i | i -> Printf.sprintf " (%d-)" (-i) in - let mult = Multiplicity.to_string (get_multiplicity m) in - let { nuclei ; elec_alpha ; elec_beta } = m in + let mult = + get_multiplicity m + |> Multiplicity.to_string + in + let { nuclei ; elec_alpha ; elec_beta } = m + in let rec build_list accu = function | a::rest -> begin @@ -53,7 +66,9 @@ let name m = in let rec build_name accu = function | (a, n)::rest -> - let a = Element.to_string a in + let a = + Element.to_string a + in begin match n with | 1 -> build_name (a::accu) rest @@ -64,19 +79,25 @@ let name m = end | [] -> accu in - let result = build_list [] nuclei |> build_name [c ; ", " ; mult] + let result = + build_list [] nuclei |> build_name [c ; ", " ; mult] in String.concat (result) -;; + let to_string m = - let { nuclei ; elec_alpha ; elec_beta } = m in - let n = List.length nuclei in - let title = name m in - [ Int.to_string n ; title ] @ (List.map ~f:(fun x -> Atom.to_string - Units.Angstrom x) nuclei) + let { nuclei ; elec_alpha ; elec_beta } = m + in + let n = + List.length nuclei + in + let title = + name m + in + [ Int.to_string n ; title ] @ + (List.map ~f:(fun x -> Atom.to_string Units.Angstrom x) nuclei) |> String.concat ~sep:"\n" -;; + let of_xyz_string ?(charge=(Charge.of_int 0)) ?(multiplicity=(Multiplicity.of_int 1)) @@ -94,7 +115,9 @@ let of_xyz_string ) + 1 - (Charge.to_int charge) |> Elec_number.of_int in - let (na,nb) = Multiplicity.to_alpha_beta ne multiplicity in + let (na,nb) = + Multiplicity.to_alpha_beta ne multiplicity + in let result = { nuclei = l ; elec_alpha = na ; @@ -109,7 +132,7 @@ let of_xyz_string raise (MultiplicityError msg); else () ; result -;; + let of_xyz_file @@ -121,8 +144,33 @@ let of_xyz_file let (_,buffer) = String.lsplit2_exn buffer ~on:'\n' in of_xyz_string ~charge:charge ~multiplicity:multiplicity ~units:units buffer -;; -include To_md5;; + + +let distance_matrix molecule = + let coord = + molecule.nuclei + |> List.map ~f:(fun x -> x.Atom.coord) + |> Array.of_list + in + let n = + Array.length coord + in + let result = + Array.make_matrix ~dimx:n ~dimy:n 0. + in + for i = 0 to (n-1) + do + for j = 0 to (n-1) + do + result.(i).(j) <- Point3d.distance coord.(i) coord.(j) + done; + done; + result + + + + +include To_md5 let to_md5 = to_md5 sexp_of_t -;; + diff --git a/ocaml/Molecule.mli b/ocaml/Molecule.mli index 176441d4..1a3d9715 100644 --- a/ocaml/Molecule.mli +++ b/ocaml/Molecule.mli @@ -34,5 +34,9 @@ val of_xyz_string : ?multiplicity:Multiplicity.t -> ?units:Units.units -> string -> t +(** Creates the distance matrix between all the atoms *) +val distance_matrix : + t -> (float array) array + (** Computes the MD5 hash *) val to_md5 : t -> Qptypes.MD5.t diff --git a/ocaml/Point3d.ml b/ocaml/Point3d.ml index 18f091f1..5717ed39 100644 --- a/ocaml/Point3d.ml +++ b/ocaml/Point3d.ml @@ -7,9 +7,16 @@ type t = { z : float ; } with sexp +let of_tuple ~units (x,y,z) = + let f = match units with + | Units.Bohr -> 1. + | Units.Angstrom -> Units.angstrom_to_bohr + in + { x = x *. f ; y = y *. f ; z = z *. f } + (** Read x y z coordinates in string s with units u *) -let of_string u s = - let f = match u with +let of_string ~units s = + let f = match units with | Units.Bohr -> 1. | Units.Angstrom -> Units.angstrom_to_bohr in @@ -22,7 +29,6 @@ let of_string u s = { x = l.(0) *. f ; y = l.(1) *. f ; z = l.(2) *. f } -;; let distance2 p1 p2 = @@ -30,17 +36,18 @@ let distance2 p1 p2 = and { x=x2 ; y=y2 ; z=z2 } = p2 in (x2-.x1)*.(x2-.x1) +. (y2-.y1)*.(y2-.y1) +. (z2-.z1)*.(z2-.z1) |> Positive_float.of_float -;; -let distance p1 p2 = sqrt (Positive_float.to_float (distance2 p1 p2)) -;; -let to_string u p = - let f = match u with +let distance p1 p2 = + sqrt (Positive_float.to_float (distance2 p1 p2)) + + +let to_string ~units p = + let f = match units with | Units.Bohr -> 1. | Units.Angstrom -> Units.bohr_to_angstrom in let { x=x ; y=y ; z=z } = p in Printf.sprintf "%16.8f %16.8f %16.8f" (x*.f) (y*.f) (z*.f) -;; + diff --git a/ocaml/Point3d.mli b/ocaml/Point3d.mli index 066a4514..6d7428ec 100644 --- a/ocaml/Point3d.mli +++ b/ocaml/Point3d.mli @@ -4,11 +4,14 @@ type t = z : float; } with sexp +(** Create from a tuple of floats *) +val of_tuple : units:Units.units -> float*float*float -> t + (** Create from an xyz string *) -val of_string : Units.units -> string -> t +val of_string : units:Units.units -> string -> t (** Convert to a string for printing *) -val to_string : Units.units -> t -> string +val to_string : units:Units.units -> t -> string (** Computes the squared distance between 2 points *) val distance2 : t -> t -> Qptypes.Positive_float.t diff --git a/ocaml/qp_create_ezfio_from_xyz.ml b/ocaml/qp_create_ezfio_from_xyz.ml index 544e6e09..fb10ff2f 100644 --- a/ocaml/qp_create_ezfio_from_xyz.ml +++ b/ocaml/qp_create_ezfio_from_xyz.ml @@ -11,26 +11,92 @@ let spec = ~doc:"string Name of basis set." +> flag "c" (optional_with_default 0 int) ~doc:"int Total charge of the molecule. Default is 0." + +> flag "d" (optional_with_default 0. float) + ~doc:"float Add dummy atoms. x * (covalent radii of the atoms)" +> flag "m" (optional_with_default 1 int) ~doc:"int Spin multiplicity (2S+1) of the molecule. Default is 1." +> flag "p" no_arg ~doc:" Using pseudopotentials" +> anon ("xyz_file" %: string) -;; -let run ?o b c m p xyz_file = + +let dummy_centers ~threshold ~molecule ~nuclei = + let d = + Molecule.distance_matrix molecule + in + let n = + Array.length d + in + let nuclei = + Array.of_list nuclei + in + let rec aux accu = function + | (-1,_) -> accu + | (i,-1) -> aux accu (i-1,i-1) + | (i,j) when (i>j) -> + let new_accu = + let x,y = + Element.covalent_radius (nuclei.(i)).Atom.element |> Positive_float.to_float, + Element.covalent_radius (nuclei.(j)).Atom.element |> Positive_float.to_float + in + let r = + ( x +. y ) *. threshold + in + if d.(i).(j) < r then + (i,x,j,y,d.(i).(j)) :: accu + else + accu + in aux new_accu (i,j-1) + | (i,j) when (i=j) -> aux accu (i,j-1) + | _ -> assert false + in + aux [] (n-1,n-1) + |> List.map ~f:(fun (i,x,j,y,r) -> + let f = + x /. (x +. y) + in + let u = + Point3d.of_tuple ~units:Units.Bohr + ( nuclei.(i).Atom.coord.Point3d.x +. + (nuclei.(j).Atom.coord.Point3d.x -. nuclei.(i).Atom.coord.Point3d.x) *. f, + nuclei.(i).Atom.coord.Point3d.y +. + (nuclei.(j).Atom.coord.Point3d.y -. nuclei.(i).Atom.coord.Point3d.y) *. f, + nuclei.(i).Atom.coord.Point3d.z +. + (nuclei.(j).Atom.coord.Point3d.z -. nuclei.(i).Atom.coord.Point3d.z) *. f) + in + Atom.{ element = Element.X ; charge = Charge.of_int 0 ; coord = u } + ) + + + + + +let run ?o b c d m p xyz_file = (* Read molecule *) let molecule = (Molecule.of_xyz_file xyz_file ~charge:(Charge.of_int c) ~multiplicity:(Multiplicity.of_int m) ) in - let nuclei = molecule.Molecule.nuclei in + let dummy = + dummy_centers ~threshold:d ~molecule ~nuclei:molecule.Molecule.nuclei + in +(* + List.iter dummy ~f:(fun x -> + Printf.printf "%s\n" (Atom.to_string ~units:Units.Angstrom x) + ); +*) + let nuclei = + molecule.Molecule.nuclei @ dummy + in + let basis_table = Hashtbl.Poly.create () in (* Open basis set channels *) let basis_channel element = - let key = Element.to_string element in + let key = + Element.to_string element + in match Hashtbl.find basis_table key with | Some in_channel -> in_channel @@ -80,7 +146,8 @@ let run ?o b c m p xyz_file = in Unix.unlink filename; List.iter nuclei ~f:(fun elem-> - let key = Element.to_string elem.Atom.element + let key = + Element.to_string elem.Atom.element in match Hashtbl.add basis_table ~key:key ~data:new_channel with | `Ok -> () @@ -92,7 +159,8 @@ let run ?o b c m p xyz_file = let elem = Element.of_string key and basis = String.lowercase basis in - let key = Element.to_string elem + let key = + Element.to_string elem in let command = Qpackage.root ^ "/scripts/get_basis.sh \"" ^ temp_filename ^ @@ -175,7 +243,12 @@ let run ?o b c m p xyz_file = |> List.rev |> List.map ~f:(fun (x,i) -> try - Basis.read_element (basis_channel x.Atom.element) i x.Atom.element + let e = + match x.Atom.element with + | Element.X -> Element.H + | e -> e + in + Basis.read_element (basis_channel x.Atom.element) i e with | End_of_file -> begin @@ -264,7 +337,7 @@ let run ?o b c m p xyz_file = | None -> failwith "Error in basis" | Some x -> Input.Ao_basis.write x -;; + let command = Command.basic @@ -279,13 +352,13 @@ elements can be defined as follows: ") spec - (fun o b c m p xyz_file () -> - run ?o b c m p xyz_file ) -;; + (fun o b c d m p xyz_file () -> + run ?o b c d m p xyz_file ) + let () = Command.run command -;; + diff --git a/ocaml/test_molecule.ml b/ocaml/test_molecule.ml index 3abd7e9a..5d2935ed 100644 --- a/ocaml/test_molecule.ml +++ b/ocaml/test_molecule.ml @@ -39,6 +39,15 @@ H 0.54386314 0.00000000 -0.92559535 let m = Molecule.of_xyz_file "c2h6.xyz" in print_string (Molecule.to_string m); + print_string "\nDistance matrix\n"; + print_string "---------------\n"; + let d = + Molecule.distance_matrix m + in + Array.iter d ~f:(fun x -> + Array.iter x ~f:(fun y -> Printf.printf "%12.8f " y); + print_newline (); + ) ;; test_molecule ();; diff --git a/scripts/get_basis.sh b/scripts/get_basis.sh index 1e969d5c..4d1fc61d 100755 --- a/scripts/get_basis.sh +++ b/scripts/get_basis.sh @@ -9,7 +9,6 @@ # - if [[ -z ${QP_ROOT} ]] then print "The QP_ROOT environment variable is not set." diff --git a/scripts/pseudo/elts_num_ele.py b/scripts/pseudo/elts_num_ele.py index 3c4ad09f..8f31f4f7 100644 --- a/scripts/pseudo/elts_num_ele.py +++ b/scripts/pseudo/elts_num_ele.py @@ -1,4 +1,5 @@ -name_to_elec = {"H": 1, +name_to_elec = {"X": 0, + "H": 1, "He": 2, "Li": 3, "Be": 4, diff --git a/scripts/pseudo/put_pseudo_in_ezfio.py b/scripts/pseudo/put_pseudo_in_ezfio.py index f44fb097..6ad69f10 100755 --- a/scripts/pseudo/put_pseudo_in_ezfio.py +++ b/scripts/pseudo/put_pseudo_in_ezfio.py @@ -58,17 +58,35 @@ def get_pseudo_str(l_atom): str_ = "" for a in l_atom: - l_cmd_atom = ["--atom", a] - l_cmd_head = [EMSL_path, "get_basis_data", - "--db_path", db_path, - "--basis", "BFD-Pseudo"] + if a is not 'X': + l_cmd_atom = ["--atom", a] - process = Popen(l_cmd_head + l_cmd_atom, stdout=PIPE, stderr=PIPE) + l_cmd_head = [EMSL_path, "get_basis_data", + "--db_path", db_path, + "--basis", "BFD-Pseudo"] - stdout, _ = process.communicate() - str_ += stdout.strip() + "\n" + process = Popen(l_cmd_head + l_cmd_atom, stdout=PIPE, stderr=PIPE) + stdout, _ = process.communicate() + str_ += stdout.strip() + "\n" + + else: # Dummy atoms + str_ += """Element Symbol: X +Number of replaced protons: 0 +Number of projectors: 0 + +Pseudopotential data: + +Local component: +Coeff. r^n Exp. +0.0 -1 0. +0.0 1 0. +0.0 0 0. + +Non-local component: +Coeff. r^n Exp. Proj. +""" return str_ diff --git a/src/Integrals_Bielec/map_integrals.irp.f b/src/Integrals_Bielec/map_integrals.irp.f index 181075bb..23deec02 100644 --- a/src/Integrals_Bielec/map_integrals.irp.f +++ b/src/Integrals_Bielec/map_integrals.irp.f @@ -525,7 +525,8 @@ integer function load_$ao_integrals(filename) integer*8 :: i integer(cache_key_kind), pointer :: key(:) real(integral_kind), pointer :: val(:) - integer :: iknd, kknd, n, j + integer :: iknd, kknd + integer*8 :: n, j load_$ao_integrals = 1 open(unit=66,file=filename,FORM='unformatted',STATUS='UNKNOWN') read(66,err=98,end=98) iknd, kknd @@ -555,12 +556,8 @@ integer function load_$ao_integrals(filename) return 99 continue call map_deinit($ao_integrals_map) - FREE $ao_integrals_map - if (.True.) then - PROVIDE $ao_integrals_map - endif - stop 'Problem reading $ao_integrals_map file in work/' 98 continue + stop 'Problem reading $ao_integrals_map file in work/' end