Merge pull request #115 from scemama/master

Dummy atoms + Memory bug
2025-01-10 21:18:29 +01:00 · 2015-11-18 00:51:38 +01:00 · 2015-11-18 00:51:38 +01:00 · 4706f21b8d
commit 4706f21b8d
parent d9e73e500a 40cf14a5e5
20 changed files with 618 additions and 310 deletions
--- 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:"   "
 ;;

--- 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
--- 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,7 +9,7 @@ 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
@ -50,7 +51,7 @@ let of_string x =
 |  "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

--- 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
--- 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 "
--- 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
-;;
+
--- 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
--- 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)
-;;
+
  
--- 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
--- 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
-;;
+



--- 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 ();;
--- 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."
--- 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,
--- a/scripts/pseudo/put_pseudo_in_ezfio.py
+++ b/scripts/pseudo/put_pseudo_in_ezfio.py
@ -58,6 +58,8 @@ def get_pseudo_str(l_atom):
    str_ = ""

    for a in l_atom:
+
+        if a is not 'X':
            l_cmd_atom = ["--atom", a]

            l_cmd_head = [EMSL_path, "get_basis_data",
@ -69,6 +71,22 @@ def get_pseudo_str(l_atom):
            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_


--- a/src/Determinants/davidson.irp.f
+++ b/src/Determinants/davidson.irp.f
@ -135,19 +135,19 @@ end subroutine
 subroutine sort_dets_ba_v(key_in, key_out, idx, shortcut, version, N_key, Nint)
  use bitmasks
  implicit none
+  BEGIN_DOC
+! Uncodumented : TODO
+  END_DOC
  integer, intent(in)            :: Nint, N_key
  integer(bit_kind),intent(in)   :: key_in(Nint,2,N_key)
-  integer(bit_kind)              :: key(Nint,2,N_key)
  integer(bit_kind),intent(out)  :: key_out(Nint,N_key)
  integer,intent(out)            :: idx(N_key)
  integer,intent(out)            :: shortcut(0:N_key+1)
  integer(bit_kind),intent(out)  :: version(Nint,N_key+1)
-  integer(bit_kind)              :: tmp(Nint, 2,N_key)
+  integer(bit_kind), allocatable :: key(:,:,:)
  integer                        :: i,ni
  
-  BEGIN_DOC
-! Uncodumented : TODO
-  END_DOC
+  allocate ( key(Nint,2,N_key) )
  do i=1,N_key
    do ni=1,Nint
      key(ni,1,i) = key_in(ni,2,i)
@ -155,8 +155,8 @@ subroutine sort_dets_ba_v(key_in, key_out, idx, shortcut, version, N_key, Nint)
    enddo
  enddo
  
-
  call sort_dets_ab_v(key, key_out, idx, shortcut, version, N_key, Nint)
+  deallocate ( key )
 end subroutine


@ -170,14 +170,15 @@ subroutine sort_dets_ab_v(key_in, key_out, idx, shortcut, version, N_key, Nint)
  END_DOC
  integer, intent(in)                   :: Nint, N_key
  integer(bit_kind),intent(in)          :: key_in(Nint,2,N_key)
-  integer(bit_kind)                     :: key(Nint,2,N_key)
  integer(bit_kind),intent(out)         :: key_out(Nint,N_key)
  integer,intent(out)                   :: idx(N_key)
  integer,intent(out)                   :: shortcut(0:N_key+1)
  integer(bit_kind),intent(out)         :: version(Nint,N_key+1)
+  integer(bit_kind), allocatable        :: key(:,:,:)
  integer(bit_kind)                     :: tmp(Nint, 2)
  integer                               :: tmpidx,i,ni
  
+  allocate (key(Nint,2,N_key))
  do i=1,N_key
    do ni=1,Nint
      key(ni,1,i) = key_in(ni,1,i)
@ -226,6 +227,7 @@ subroutine sort_dets_ab_v(key_in, key_out, idx, shortcut, version, N_key, Nint)
      key_out(ni,i) = key(ni,2,i)
    enddo
  enddo
+  deallocate (key)
 end subroutine


--- a/src/Determinants/s2.irp.f
+++ b/src/Determinants/s2.irp.f
@ -117,11 +117,12 @@ subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2)
  integer                        :: i,j,l,jj,ii
  integer, allocatable           :: idx(:)
  
-  integer                        :: shortcut(0:n+1), sort_idx(n)
-  integer(bit_kind)              :: sorted(N_int,n), version(N_int,n)
+  integer, allocatable           :: shortcut(:), sort_idx(:)
+  integer(bit_kind), allocatable :: sorted(:,:), version(:,:)
  integer                        :: sh, sh2, ni, exa, ext, org_i, org_j, endi, pass
  double precision               :: davidson_threshold_bis
  
+  allocate (shortcut(0:n+1), sort_idx(n), sorted(N_int,n), version(N_int,n))
  s2 = 0.d0
  davidson_threshold_bis = threshold_davidson
  call sort_dets_ab_v(psi_keys_tmp, sorted, sort_idx, shortcut, version, n, N_int)
@ -210,6 +211,7 @@ subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2)
    s2 = s2 + psi_coefs_tmp(i)*psi_coefs_tmp(i)*s2_tmp
  enddo
  s2 = s2 + S_z2_Sz
+  deallocate (shortcut, sort_idx, sorted, version)
 end


--- a/src/Determinants/slater_rules.irp.f
+++ b/src/Determinants/slater_rules.irp.f
@ -1237,11 +1237,10 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint)
  integer                        :: i,j,k,l, jj,ii
  integer                        :: i0, j0
  
-  integer                        :: shortcut(0:n+1), sort_idx(n)
-  integer(bit_kind)              :: sorted(Nint,n), version(Nint,n)
+  integer, allocatable           :: shortcut(:), sort_idx(:)
+  integer(bit_kind), allocatable :: sorted(:,:), version(:,:)
  integer(bit_kind)              :: sorted_i(Nint)
  
-  
  integer                        :: sh, sh2, ni, exa, ext, org_i, org_j, endi
  double precision               :: local_threshold
  
@ -1250,6 +1249,8 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint)
  ASSERT (Nint == N_int)
  ASSERT (n>0)
  PROVIDE ref_bitmask_energy davidson_criterion
+
+  allocate (shortcut(0:n+1), sort_idx(n), sorted(Nint,n), version(Nint,n))
  v_0 = 0.d0

  call sort_dets_ab_v(keys_tmp, sorted, sort_idx, shortcut, version, n, Nint)
@ -1353,5 +1354,6 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint)
  do i=1,n
    v_0(i) += H_jj(i) * u_0(i)
  enddo
+  deallocate (shortcut, sort_idx, sorted, version)
 end

--- 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

--- a/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f
+++ b/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f
@ -5,17 +5,15 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integral, (ao_num_align,ao_num)]
  END_DOC
  if (do_pseudo) then
    ao_pseudo_integral = ao_pseudo_integral_local + ao_pseudo_integral_non_local 
-!    ao_pseudo_integral = ao_pseudo_integral_local 
-!    ao_pseudo_integral = ao_pseudo_integral_non_local 
  else
    ao_pseudo_integral = 0.d0
  endif
 END_PROVIDER

- BEGIN_PROVIDER [ double precision, ao_pseudo_integral_local, (ao_num_align,ao_num)]
+BEGIN_PROVIDER [ double precision, ao_pseudo_integral_local, (ao_num_align,ao_num)]
  implicit none
  BEGIN_DOC
-! Local pseudo-potential
+  ! Local pseudo-potential
  END_DOC
  double precision               :: alpha, beta, gama, delta
  integer                        :: num_A,num_B
@ -26,6 +24,7 @@ END_PROVIDER
  
  double precision               :: cpu_1, cpu_2, wall_1, wall_2, wall_0
  integer                        :: thread_num
+  !$ integer                     :: omp_get_thread_num
  
  ao_pseudo_integral_local = 0.d0
  
@ -35,28 +34,24 @@ END_PROVIDER
  
  allocate(n_k_dump(1:pseudo_klocmax), v_k_dump(1:pseudo_klocmax), dz_k_dump(1:pseudo_klocmax))
  
-
-  !  _                
-  ! /   _. |  _     | 
-  ! \_ (_| | (_ |_| | 
-  !                   
-
-  print*, 'Providing the nuclear electron pseudo integrals '
+  print*, 'Providing the nuclear electron pseudo integrals (local)'
  
  call wall_time(wall_1)
  call cpu_time(cpu_1)
  
+  thread_num = 0
  !$OMP PARALLEL                                                     &
      !$OMP DEFAULT (NONE)                                           &
-  !$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B, &
+      !$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B,&
      !$OMP          num_A,num_B,Z,c,n_pt_in,                        &
      !$OMP          v_k_dump,n_k_dump, dz_k_dump,                   &
      !$OMP          wall_0,wall_2,thread_num)                       &
-  !$OMP SHARED (ao_num,ao_prim_num,ao_expo_ordered_transp,ao_power,ao_nucl,nucl_coord,ao_coef_normalized_ordered_transp, &
+      !$OMP SHARED (ao_num,ao_prim_num,ao_expo_ordered_transp,ao_power,ao_nucl,nucl_coord,ao_coef_normalized_ordered_transp,&
      !$OMP         ao_pseudo_integral_local,nucl_num,nucl_charge,   &
-  !$OMP         pseudo_klocmax,pseudo_lmax,pseudo_kmax,pseudo_v_k,pseudo_n_k, pseudo_dz_k, &
+      !$OMP         pseudo_klocmax,pseudo_lmax,pseudo_kmax,pseudo_v_k,pseudo_n_k, pseudo_dz_k,&
      !$OMP         wall_1)
  
+  !$ thread_num = omp_get_thread_num()
  !$OMP DO SCHEDULE (guided)
  
  do j = 1, ao_num
@ -79,6 +74,10 @@ END_PROVIDER
          double precision               :: c
          c = 0.d0
          
+          if (dabs(ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i))&
+                < 1.d-10) then
+            cycle
+          endif
          do  k = 1, nucl_num
            double precision               :: Z
            Z = nucl_charge(k)
@ -89,11 +88,11 @@ END_PROVIDER
            n_k_dump = pseudo_n_k(k,1:pseudo_klocmax)
            dz_k_dump =  pseudo_dz_k(k,1:pseudo_klocmax)
            
-        c = c + Vloc(pseudo_klocmax, v_k_dump,n_k_dump, dz_k_dump, &
+            c = c + Vloc(pseudo_klocmax, v_k_dump,n_k_dump, dz_k_dump,&
                A_center,power_A,alpha,B_center,power_B,beta,C_center)
            
          enddo
-      ao_pseudo_integral_local(i,j) = ao_pseudo_integral_local(i,j) + &
+          ao_pseudo_integral_local(i,j) = ao_pseudo_integral_local(i,j) +&
              ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i)*c
        enddo
      enddo
@ -121,7 +120,7 @@ END_PROVIDER
 BEGIN_PROVIDER [ double precision, ao_pseudo_integral_non_local, (ao_num_align,ao_num)]
  implicit none
  BEGIN_DOC
-! Local pseudo-potential
+  ! Local pseudo-potential
  END_DOC
  double precision               :: alpha, beta, gama, delta
  integer                        :: num_A,num_B
@ -129,6 +128,7 @@ END_PROVIDER
  integer                        :: power_A(3),power_B(3)
  integer                        :: i,j,k,l,n_pt_in,m
  double precision               :: Vloc, Vpseudo
+  !$ integer                     :: omp_get_thread_num
  
  double precision               :: cpu_1, cpu_2, wall_1, wall_2, wall_0
  integer                        :: thread_num
@ -141,27 +141,23 @@ END_PROVIDER
  
  allocate(n_kl_dump(pseudo_kmax,0:pseudo_lmax), v_kl_dump(pseudo_kmax,0:pseudo_lmax), dz_kl_dump(pseudo_kmax,0:pseudo_lmax))
  
-  !  _                
-  ! /   _. |  _     | 
-  ! \_ (_| | (_ |_| | 
-  !                   
-
-  print*, 'Providing the nuclear electron pseudo integrals '
+  print*, 'Providing the nuclear electron pseudo integrals (non-local)'
  
  call wall_time(wall_1)
  call cpu_time(cpu_1)
-
+  thread_num = 0
  !$OMP PARALLEL                                                     &
      !$OMP DEFAULT (NONE)                                           &
-  !$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B, &
+      !$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B,&
      !$OMP          num_A,num_B,Z,c,n_pt_in,                        &
      !$OMP          n_kl_dump, v_kl_dump, dz_kl_dump,               &
      !$OMP          wall_0,wall_2,thread_num)                       &
-  !$OMP SHARED (ao_num,ao_prim_num,ao_expo_ordered_transp,ao_power,ao_nucl,nucl_coord,ao_coef_normalized_ordered_transp, &
-  !$OMP         ao_pseudo_integral_non_local,nucl_num,nucl_charge, &
-  !$OMP         pseudo_klocmax,pseudo_lmax,pseudo_kmax,pseudo_n_kl, pseudo_v_kl, pseudo_dz_kl, &
+      !$OMP SHARED (ao_num,ao_prim_num,ao_expo_ordered_transp,ao_power,ao_nucl,nucl_coord,ao_coef_normalized_ordered_transp,&
+      !$OMP         ao_pseudo_integral_non_local,nucl_num,nucl_charge,&
+      !$OMP         pseudo_klocmax,pseudo_lmax,pseudo_kmax,pseudo_n_kl, pseudo_v_kl, pseudo_dz_kl,&
      !$OMP         wall_1)
  
+  !$ thread_num = omp_get_thread_num()
  !$OMP DO SCHEDULE (guided)
  
  do j = 1, ao_num
@ -184,6 +180,11 @@ END_PROVIDER
          double precision               :: c
          c = 0.d0
          
+          if (dabs(ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i))&
+                < 1.d-10) then
+            cycle
+          endif
+          
          do  k = 1, nucl_num
            double precision               :: Z
            Z = nucl_charge(k)
@ -197,7 +198,7 @@ END_PROVIDER
            c = c + Vpseudo(pseudo_lmax,pseudo_kmax,v_kl_dump,n_kl_dump,dz_kl_dump,A_center,power_A,alpha,B_center,power_B,beta,C_center)
            
          enddo
-      ao_pseudo_integral_non_local(i,j) = ao_pseudo_integral_non_local(i,j) + &
+          ao_pseudo_integral_non_local(i,j) = ao_pseudo_integral_non_local(i,j) +&
              ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i)*c
        enddo
      enddo
@ -220,7 +221,7 @@ END_PROVIDER
  deallocate(n_kl_dump,v_kl_dump, dz_kl_dump)
  
  
- END_PROVIDER
+END_PROVIDER



--- a/src/Integrals_Monoelec/pseudopot.f90
+++ b/src/Integrals_Monoelec/pseudopot.f90
@ -197,8 +197,8 @@ integer, intent(in) :: n_a(3),n_b(3)
 double precision, intent(in) :: v_kl(kmax,0:lmax),dz_kl(kmax,0:lmax)

 !                     
-! |   _   _  _. |  _  
-! |_ (_) (_ (_| | (/_ 
+! |   _   _  _. | 
+! |_ (_) (_ (_| | 
 !                     

 double precision :: fourpi,f,prod,prodp,binom_func,accu,bigR,bigI,ylm
@ -223,11 +223,6 @@ double precision, allocatable :: array_I_B(:,:,:,:,:)

 double precision :: f1, f2, f3

-!  _                
-! /   _. |  _     | 
-! \_ (_| | (_ |_| | 
-!                   
-
 if (kmax.eq.1.and.lmax.eq.0.and.v_kl(1,0).eq.0.d0) then
  Vpseudo=0.d0
  return
@ -236,7 +231,7 @@ end if
 fourpi=4.d0*dacos(-1.d0)
 ac=dsqrt((a(1)-c(1))**2+(a(2)-c(2))**2+(a(3)-c(3))**2)
 bc=dsqrt((b(1)-c(1))**2+(b(2)-c(2))**2+(b(3)-c(3))**2)
-arg=g_a*ac**2+g_b*bc**2
+arg= g_a*ac*ac + g_b*bc*bc

 if(arg.gt.-dlog(1.d-20))then
  Vpseudo=0.d0
@ -290,6 +285,21 @@ if(ac.eq.0.d0.and.bc.eq.0.d0)then
   enddo
  enddo
 enddo
+! do k=1,kmax
+!   do l=0,lmax
+!     ktot=ntot+n_kl(k,l)
+!     do m=-l,l
+!       prod =bigI(0,0,l,m,n_a(1),n_a(2),n_a(3))*v_kl(k,l)
+!       prodp=bigI(0,0,l,m,n_b(1),n_b(2),n_b(3))*prod
+!       if (dabs (prodp) < 1.d-15) then
+!         cycle
+!       endif
+!       
+!       accu=accu+prodp*int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),0,0,areal,breal,arg)
+!       
+!     enddo
+!   enddo
+! enddo

 !=!=!=!=!
 ! E n d !
@ -625,8 +635,8 @@ double precision, intent(in) :: v_kl(kmax,0:lmax),dz_kl(kmax,0:lmax)
 double precision, intent(in) :: rmax

 !                     
-! |   _   _  _. |  _  
-! |_ (_) (_ (_| | (/_ 
+! |   _   _  _. |
+! |_ (_) (_ (_| |
 !                     

 integer :: l,m,k,kk
@ -1950,6 +1960,7 @@ double precision function int_prod_bessel(l,gam,n,m,a,b,arg)
  double precision :: s_q_0, s_q_k, s_0_0, a_over_b_square
  double precision :: int_prod_bessel_loc
  double precision :: inverses(0:300)
+  double precision :: two_qkmp1, qk

  logical done

@ -2008,19 +2019,18 @@ double precision function int_prod_bessel(l,gam,n,m,a,b,arg)
        stop 'pseudopot.f90 : q > 300'
      endif

+      two_qkmp1 = dble(2*(q+m)+1)
+      qk = dble(q)
      do k=0,q-1
-        s_q_k = ( dble(2*(q-k+m)+1)*dble(q-k)*inverses(k) ) * s_q_k
+        s_q_k = ( two_qkmp1*qk*inverses(k) ) * s_q_k
        sum=sum+s_q_k
+        two_qkmp1 = two_qkmp1-2.d0
+        qk = qk-1.d0 
      enddo
      inverses(q) = a_over_b_square/(dble(2*(q+n)+3) * dble(q+1)) 
 !      do k=0,q
 !        sum=sum+s_q_k
 !        s_q_k = a_over_b_square * ( dble(2*(q-k+m)+1)*dble(q-k)/(dble(2*(k+n)+3) * dble(k+1)) ) * s_q_k
-!      enddo
-      ! Iteration of k
-!      do k=0,q
-!        sum=sum+s_q_k
-!        s_q_k = a_over_b_square * ( dble(2*(q-k+m)+1)*dble(q-k)/(dble(2*(k+n)+3) * dble(k+1)) ) * s_q_k
 !      enddo
    
      int=int+sum