Merge branch 'cleaning_dft' of https://github.com/QuantumPackage/qp2 into cleaning_dft

README.md

@@ -36,8 +36,8 @@ https://arxiv.org/abs/1902.08154
 
 # Build status
 
-* Master [![master build status](https://travis-ci.org/QuantumPackage/qp2.svg?branch=master)](https://travis-ci.org/QuantumPackage/qp2)
+* Master [![master build status](https://travis-ci.com/QuantumPackage/qp2.svg?branch=master)](https://travis-ci.org/QuantumPackage/qp2)
-* Development [![dev build status](https://travis-ci.org/QuantumPackage/qp2.svg?branch=dev)](https://travis-ci.org/QuantumPackage/qp2)
+* Development [![dev build status](https://travis-ci.com/QuantumPackage/qp2.svg?branch=dev)](https://travis-ci.org/QuantumPackage/qp2)
 * Documentation [![Documentation Status](https://readthedocs.org/projects/quantum-package/badge/?version=master)](https://quantum-package.readthedocs.io/en/master/?badge=master)
 
 

RELEASE_NOTES.org

@@ -29,6 +29,7 @@
   - Disk-based Davidson when too much memory is required
   - Fixed bug in DIIS
   - Fixed bug in molden (Au -> Angs)
+  - Fixed bug with non-contiguous MOs in active space and deleter MOs
 
 *** User interface
 
@@ -52,6 +53,10 @@
     - Added ~print_hamiltonian~
     - Added input for two body RDM
     - Added keyword ~save_wf_after_selection~
+    - Added a ~restore_symm~ flag to enforce the restoration of
+      symmetry in matrices
+    - qp_export_as_tgz exports also plugin codes
+    - Added a basis module containing basis set information
 
 *** Code
 
@@ -76,6 +81,8 @@
     - Added ~h_core_guess~ routine
     - Fixed Laplacians in real space (indices)
     - Added LIB file to add extra libs in plugin
+    - Using Intel IPP for sorting when using Intel compiler
+    - Removed parallelism in sorting
 
     ao_one_e_integral_zero
     banned_excitations

VERSION

@@ -0,0 +1 @@
+2.2.1

bin/qp_export_as_tgz

@@ -99,7 +99,9 @@ function find_libs () {
 }
 
 function find_exec () {
-  find ${QP_ROOT}/$1 -perm /u+x -type f
+  for i in $@ ; do
+     find ${QP_ROOT}/$i -perm /u+x -type f
+  done
 }
 
 
@@ -119,7 +121,7 @@ fi
 echo "Copying binary files"
 #     --------------------
 
-FORTRAN_EXEC=$(find_exec src)
+FORTRAN_EXEC=$(find_exec src/*/)
 if [[ -z $FORTRAN_EXEC ]] ; then
   error 'No Fortran binaries found.'
   exit 1

config/ifort.cfg

@@ -7,9 +7,9 @@
 #
 [COMMON]
 FC           : ifort -fpic
-LAPACK_LIB   : -mkl=parallel
+LAPACK_LIB   : -mkl=parallel -lirc -lsvml -limf -lipps
 IRPF90       : irpf90
-IRPF90_FLAGS : --ninja --align=32 
+IRPF90_FLAGS : --ninja --align=32  -DINTEL
 
 # Global options
 ################

config/ifort_avx.cfg

@@ -7,9 +7,9 @@
 #
 [COMMON]
 FC           : ifort -fpic
-LAPACK_LIB   : -mkl=parallel
+LAPACK_LIB   : -mkl=parallel -lirc -lsvml -limf -lipps
 IRPF90       : irpf90
-IRPF90_FLAGS : --ninja --align=32
+IRPF90_FLAGS : --ninja --align=32 -DINTEL
 
 # Global options
 ################

config/ifort_avx_mpi.cfg

@@ -7,9 +7,9 @@
 #
 [COMMON]
 FC           : mpiifort -fpic
-LAPACK_LIB   : -mkl=parallel
+LAPACK_LIB   : -mkl=parallel -lirc -lsvml -limf -lipps
 IRPF90       : irpf90
-IRPF90_FLAGS : --ninja --align=32 -DMPI
+IRPF90_FLAGS : --ninja --align=32 -DMPI -DINTEL
 
 # Global options
 ################

config/ifort_debug.cfg

@@ -7,9 +7,9 @@
 #
 [COMMON]
 FC           : ifort -fpic
-LAPACK_LIB   : -mkl=parallel
+LAPACK_LIB   : -mkl=parallel -lirc -lsvml -limf -lipps
 IRPF90       : irpf90
-IRPF90_FLAGS : --ninja --align=32  --assert
+IRPF90_FLAGS : --ninja --align=32  --assert -DINTEL
 
 # Global options
 ################

config/ifort_mpi.cfg

@@ -7,9 +7,9 @@
 #
 [COMMON]
 FC           : mpiifort -fpic
-LAPACK_LIB   : -mkl=parallel
+LAPACK_LIB   : -mkl=parallel -lirc -lsvml -limf -lipps
 IRPF90       : irpf90
-IRPF90_FLAGS : --ninja --align=32 -DMPI
+IRPF90_FLAGS : --ninja --align=32 -DMPI -DINTEL
 
 # Global options
 ################

config/ifort_rome.cfg

@@ -7,9 +7,9 @@
 #
 [COMMON]
 FC           : ifort -fpic
-LAPACK_LIB   : -mkl=parallel
+LAPACK_LIB   : -mkl=parallel -lirc -lsvml -limf -lipps
 IRPF90       : irpf90
-IRPF90_FLAGS : --ninja --align=32 
+IRPF90_FLAGS : --ninja --align=32  -DINTEL
 
 # Global options
 ################
@@ -31,7 +31,7 @@ OPENMP  : 1          ; Append OpenMP flags
 # -ftz                       : Flushes denormal results to zero
 #
 [OPT]
-FC       : -traceback
+FC       : -traceback -shared-intel
 FCFLAGS  : -O2 -ip -g  -march=core-avx2 -align array64byte -fma -ftz -fomit-frame-pointer 
 
 # Profiling flags

config/ifort_xHost.cfg

@@ -7,9 +7,9 @@
 #
 [COMMON]
 FC           : ifort -fpic
-LAPACK_LIB   : -mkl=parallel
+LAPACK_LIB   : -mkl=parallel -lirc -lsvml -limf -lipps
 IRPF90       : irpf90
-IRPF90_FLAGS : --ninja --align=64
+IRPF90_FLAGS : --ninja --align=64 -DINTEL
 
 # Global options
 ################

etc/local.rc

@@ -19,4 +19,3 @@
 # export QP_NIC=lo
 # export QP_NIC=ib0
 
-

external/ezfio

@@ -1 +1 @@
-Subproject commit ccee52d00c2cde1d628b0d34f4a247143747bf36
+Subproject commit ed1df9f3c1f51752656ca98da5693a4119add05c

external/irpf90

@@ -1 +1 @@
-Subproject commit 132a4a1661c9878d21dcbf0ac14f7fe9a3b110d0
+Subproject commit 33ca5e1018f3bbb5e695e6ee558f5dac0753b271

ocaml/Angmom.ml

@@ -26,7 +26,7 @@ let of_string = function
   | "J" | "j" -> J
   | "K" | "k" -> K
   | "L" | "l" -> L
-  | x -> raise (Failure ("Symmetry should be S|P|D|F|G|H|I|J|K|L,
+  | x -> raise (Failure ("Angmom should be S|P|D|F|G|H|I|J|K|L,
 not "^x^"."))
 
 let of_char = function
@@ -40,7 +40,7 @@ let of_char = function
   | 'J' | 'j' -> J
   | 'K' | 'k' -> K
   | 'L' | 'l' -> L
-  | x -> raise (Failure ("Symmetry should be S|P|D|F|G|H|I|J|K|L"))
+  | x -> raise (Failure ("Angmom should be S|P|D|F|G|H|I|J|K|L"))
 
 let to_l = function
   | S -> Positive_int.of_int 0
@@ -68,7 +68,7 @@ let of_l i =
   | 7 -> J
   | 8 -> K
   | 9 -> L
-  | x -> raise (Failure ("Symmetry should be S|P|D|F|G|H|I|J|K|L"))
+  | x -> raise (Failure ("Angmom should be S|P|D|F|G|H|I|J|K|L"))
 
 
 type st = t

ocaml/GaussianPrimitive.ml

@@ -2,14 +2,14 @@ open Qptypes
 open Sexplib.Std
 
 type t =
-  { sym  : Symmetry.t ;
+  { sym  : Angmom.t ;
     expo : AO_expo.t ;
   } [@@deriving sexp]
 
 let to_string p =
   let { sym = s ; expo = e } = p in
   Printf.sprintf "(%s, %22e)"
-    (Symmetry.to_string s)
+    (Angmom.to_string s)
     (AO_expo.to_float e)
 
 

ocaml/Gto.ml

@@ -9,7 +9,7 @@ type fmt =
 | Gaussian
 
 type t =
-{ sym  : Symmetry.t ;
+{ sym  : Angmom.t ;
   lc   : ((GaussianPrimitive.t * AO_coef.t) list)
 } [@@deriving sexp]
 
@@ -47,7 +47,7 @@ let read_one in_channel =
   in
   let sym_str = String.sub buffer 0 2 in
   let   n_str = String.sub buffer 2 ((String.length buffer)-2) in
-  let sym = Symmetry.of_string (String_ext.strip sym_str) in
+  let sym = Angmom.of_string (String_ext.strip sym_str) in
   let n = int_of_string (String_ext.strip n_str) in
   (* Read all the primitives *)
   let rec read_lines result = function
@@ -82,7 +82,7 @@ let read_one in_channel =
 (** Write the GTO in Gamess format *)
 let to_string_gamess { sym = sym ; lc = lc } =
   let result =
-    Printf.sprintf "%s %3d" (Symmetry.to_string sym) (List.length lc)
+    Printf.sprintf "%s %3d" (Angmom.to_string sym) (List.length lc)
   in
   let rec do_work accu i = function
   | [] -> List.rev accu
@@ -102,7 +102,7 @@ let to_string_gamess { sym = sym ; lc = lc } =
 (** Write the GTO in Gaussian format *)
 let to_string_gaussian { sym = sym ; lc = lc } =
   let result =
-    Printf.sprintf "%s %3d   1.00" (Symmetry.to_string sym) (List.length lc)
+    Printf.sprintf "%s %3d   1.00" (Angmom.to_string sym) (List.length lc)
   in
   let rec do_work accu i = function
   | [] -> List.rev accu

ocaml/Gto.mli

@@ -5,7 +5,7 @@ type fmt =
 | Gaussian
 
 type t =
-  { sym : Symmetry.t ;
+  { sym : Angmom.t ;
     lc  : (GaussianPrimitive.t * Qptypes.AO_coef.t) list;
   } [@@deriving sexp]
 

ocaml/Input_ao_basis.ml

@@ -9,7 +9,7 @@ module Ao_basis : sig
       ao_prim_num     : AO_prim_number.t array;
       ao_prim_num_max : AO_prim_number.t;
       ao_nucl         : Nucl_number.t array;
-      ao_power        : Symmetry.Xyz.t array;
+      ao_power        : Angmom.Xyz.t array;
       ao_coef         : AO_coef.t array;
       ao_expo         : AO_expo.t array;
       ao_cartesian    : bool;
@@ -32,7 +32,7 @@ end = struct
       ao_prim_num     : AO_prim_number.t array;
       ao_prim_num_max : AO_prim_number.t;
       ao_nucl         : Nucl_number.t array;
-      ao_power        : Symmetry.Xyz.t array;
+      ao_power        : Angmom.Xyz.t array;
       ao_coef         : AO_coef.t array;
       ao_expo         : AO_expo.t array;
       ao_cartesian    : bool;
@@ -87,7 +87,7 @@ end = struct
       if (data.(2*dim+i-1) > 0) then
         result.(i-1) <- result.(i-1)^"z"^(string_of_int data.(2*dim+i-1));
     done;
-    Array.map Symmetry.Xyz.of_string result
+    Array.map Angmom.Xyz.of_string result
   ;;
 
   let read_ao_coef () =
@@ -133,7 +133,7 @@ end = struct
     let ao_num = AO_number.to_int b.ao_num in
     let gto_array = Array.init (AO_number.to_int b.ao_num)
       (fun i ->
-        let s = Symmetry.Xyz.to_symmetry b.ao_power.(i) in
+        let s = Angmom.Xyz.to_symmetry b.ao_power.(i) in
         let ao_prim_num = AO_prim_number.to_int b.ao_prim_num.(i) in
         let prims = List.init ao_prim_num (fun j ->
           let prim = { GaussianPrimitive.sym  = s ;
@@ -217,9 +217,9 @@ end = struct
      let ao_power =
        let l = Array.to_list ao_power in 
        List.concat [
-         (list_map (fun a -> Positive_int.to_int a.Symmetry.Xyz.x) l) ;
+         (list_map (fun a -> Positive_int.to_int a.Angmom.Xyz.x) l) ;
-         (list_map (fun a -> Positive_int.to_int a.Symmetry.Xyz.y) l) ;
+         (list_map (fun a -> Positive_int.to_int a.Angmom.Xyz.y) l) ;
-         (list_map (fun a -> Positive_int.to_int a.Symmetry.Xyz.z) l) ]
+         (list_map (fun a -> Positive_int.to_int a.Angmom.Xyz.z) l) ]
      in
      Ezfio.set_ao_basis_ao_power(Ezfio.ezfio_array_of_list
      ~rank:2 ~dim:[| ao_num ; 3 |] ~data:ao_power) ;
@@ -409,7 +409,7 @@ end = struct
       | [] -> []
       | (i,n,x)::tail  ->
           (Printf.sprintf " %5d  %6d     %-8s\n" i (Nucl_number.to_int n)
-            (Symmetry.Xyz.to_string x)
+            (Angmom.Xyz.to_string x)
           )::(do_work tail)
       in do_work l
       |> String.concat ""
@@ -496,7 +496,7 @@ md5                     = %s
     (b.ao_nucl |> Array.to_list |> list_map Nucl_number.to_string |>
       String.concat ", ")
     (b.ao_power |> Array.to_list |> list_map (fun x->
-      "("^(Symmetry.Xyz.to_string x)^")" )|> String.concat ", ")
+      "("^(Angmom.Xyz.to_string x)^")" )|> String.concat ", ")
     (b.ao_coef  |> Array.to_list |> list_map AO_coef.to_string
       |> String.concat ", ")
     (b.ao_expo  |> Array.to_list |> list_map AO_expo.to_string

ocaml/Long_basis.ml

@@ -2,7 +2,7 @@ open Qptypes
 open Qputils
 open Sexplib.Std
 
-type t = (Symmetry.Xyz.t * Gto.t * Nucl_number.t ) list [@@deriving sexp]
+type t = (Angmom.Xyz.t * Gto.t * Nucl_number.t ) list [@@deriving sexp]
 
 let of_basis b =
   let rec do_work accu = function
@@ -10,7 +10,7 @@ let of_basis b =
     | (g,n)::tail ->
         begin
           let new_accu =
-            Symmetry.Xyz.of_symmetry g.Gto.sym
+            Angmom.Xyz.of_symmetry g.Gto.sym
             |> List.rev_map (fun x-> (x,g,n))
           in
           do_work (new_accu@accu) tail
@@ -25,7 +25,7 @@ let to_basis b =
   | [] -> List.rev accu
   | (s,g,n)::tail ->
     let first_sym =
-      Symmetry.Xyz.of_symmetry g.Gto.sym
+      Angmom.Xyz.of_symmetry g.Gto.sym
       |> List.hd
     in
     let new_accu =
@@ -42,7 +42,7 @@ let to_basis b =
 let to_string b =
   let middle = list_map (fun (x,y,z) ->
      "( "^((string_of_int (Nucl_number.to_int z)))^", "^
-     (Symmetry.Xyz.to_string x)^", "^(Gto.to_string y)
+     (Angmom.Xyz.to_string x)^", "^(Gto.to_string y)
      ^" )"
   ) b
   |> String.concat ",\n"

ocaml/Long_basis.mli

@@ -5,16 +5,16 @@ open Qptypes;;
   * all the D orbitals are converted to xx, xy, xz, yy, yx
   * etc
 *)
-type t = (Symmetry.Xyz.t * Gto.t * Nucl_number.t) list [@@deriving sexp]
+type t = (Angmom.Xyz.t * Gto.t * Nucl_number.t) list [@@deriving sexp]
 
 (** Transform a basis to a long basis *)
 val of_basis :
-  (Gto.t * Nucl_number.t) list -> (Symmetry.Xyz.t * Gto.t * Nucl_number.t) list
+  (Gto.t * Nucl_number.t) list -> (Angmom.Xyz.t * Gto.t * Nucl_number.t) list
 
 (** Transform a long basis to a basis *)
 val to_basis :
-  (Symmetry.Xyz.t * Gto.t * Nucl_number.t) list -> (Gto.t * Nucl_number.t) list
+  (Angmom.Xyz.t * Gto.t * Nucl_number.t) list -> (Gto.t * Nucl_number.t) list
 
 (** Convert the basis into its string representation *)
 val to_string :
-  (Symmetry.Xyz.t * Gto.t * Nucl_number.t) list -> string
+  (Angmom.Xyz.t * Gto.t * Nucl_number.t) list -> string

ocaml/Primitive.mli

@@ -1,5 +1,5 @@
 type t =
-{ sym : Symmetry.t;
+{ sym : Angmom.t;
   expo : Qptypes.AO_expo.t;
 } [@@deriving sexp]
 
@@ -7,5 +7,5 @@ type t =
 val to_string : t -> string
 
 (** Creation *)
-val of_sym_expo : Symmetry.t -> Qptypes.AO_expo.t -> t
+val of_sym_expo : Angmom.t -> Qptypes.AO_expo.t -> t
 

ocaml/qp_create_ezfio.ml

@@ -512,13 +512,14 @@ let run ?o b au c d m p cart xyz_file =
       let ao_num = List.length long_basis in
       Ezfio.set_ao_basis_ao_num ao_num;
       Ezfio.set_ao_basis_ao_basis b;
+      Ezfio.set_basis_basis b;
       let ao_prim_num = list_map (fun (_,g,_) -> List.length g.Gto.lc) long_basis 
       and ao_nucl = list_map (fun (_,_,n) -> Nucl_number.to_int n) long_basis 
       and ao_power=
         let l = list_map (fun (x,_,_) -> x) long_basis in
-        (list_map (fun t -> Positive_int.to_int Symmetry.Xyz.(t.x)) l)@
+        (list_map (fun t -> Positive_int.to_int Angmom.Xyz.(t.x)) l)@
-        (list_map (fun t -> Positive_int.to_int Symmetry.Xyz.(t.y)) l)@
+        (list_map (fun t -> Positive_int.to_int Angmom.Xyz.(t.y)) l)@
-        (list_map (fun t -> Positive_int.to_int Symmetry.Xyz.(t.z)) l)
+        (list_map (fun t -> Positive_int.to_int Angmom.Xyz.(t.z)) l)
       in
       let ao_prim_num_max = List.fold_left (fun s x ->
         if x > s then x
@@ -561,6 +562,78 @@ let run ?o b au c d m p cart xyz_file =
       and ao_expo = create_expo_coef `Expos
       in
       let () =
+        let shell_num = List.length basis in
+        let lc : (GaussianPrimitive.t * Qptypes.AO_coef.t) list list =
+             list_map ( fun (g,_) -> g.Gto.lc ) basis
+        in
+        let ang_mom =
+          list_map (fun (l : (GaussianPrimitive.t * Qptypes.AO_coef.t) list)  -> 
+            let x, _ = List.hd l in
+            Angmom.to_l x.GaussianPrimitive.sym |> Qptypes.Positive_int.to_int
+             ) lc
+        in
+        let expo =
+          list_map (fun l -> list_map (fun (x,_) -> Qptypes.AO_expo.to_float x.GaussianPrimitive.expo) l ) lc
+          |> List.concat
+        in
+        let coef =
+          list_map (fun l -> 
+            list_map (fun (_,x) -> Qptypes.AO_coef.to_float x) l
+             ) lc
+          |> List.concat
+        in
+        let shell_prim_num =
+          list_map List.length  lc
+        in
+        let shell_prim_idx =
+          let rec aux count accu = function
+          | [] -> List.rev accu
+          | l::rest ->
+            let newcount = count+(List.length l) in
+            aux newcount (count::accu) rest
+          in
+          aux 1 [] lc
+        in
+        let prim_num = List.length coef in
+        Ezfio.set_basis_typ "Gaussian";
+        Ezfio.set_basis_shell_num shell_num;
+        Ezfio.set_basis_prim_num prim_num ;
+        Ezfio.set_basis_shell_prim_num  (Ezfio.ezfio_array_of_list
+          ~rank:1 ~dim:[| shell_num |] ~data:shell_prim_num);
+        Ezfio.set_basis_shell_ang_mom (Ezfio.ezfio_array_of_list
+          ~rank:1 ~dim:[| shell_num |] ~data:ang_mom ) ;
+        Ezfio.set_basis_shell_prim_index (Ezfio.ezfio_array_of_list
+          ~rank:1 ~dim:[| shell_num |] ~data:shell_prim_idx) ;
+        Ezfio.set_basis_basis_nucleus_index (Ezfio.ezfio_array_of_list
+          ~rank:1 ~dim:[| nucl_num |]
+          ~data:(
+          list_map (fun (_,n) -> Nucl_number.to_int n) basis
+          |> List.fold_left (fun accu i -> 
+            match accu with 
+            | [] -> []
+            | (h,j) :: rest -> if j == i then ((h+1,j)::rest) else ((h+1,i)::(h+1,j)::rest)
+          ) [(0,0)]
+          |> List.rev
+          |> List.map fst
+          )) ;
+        Ezfio.set_basis_nucleus_shell_num(Ezfio.ezfio_array_of_list
+          ~rank:1 ~dim:[| nucl_num |]
+          ~data:(
+          list_map (fun (_,n) -> Nucl_number.to_int n) basis
+          |> List.fold_left (fun accu i -> 
+            match accu with 
+            | [] -> [(1,i)]
+            | (h,j) :: rest -> if j == i then ((h+1,j)::rest) else ((1,i)::(h,j)::rest)
+          ) []
+          |> List.rev
+          |> List.map fst
+          )) ;
+        Ezfio.set_basis_prim_coef (Ezfio.ezfio_array_of_list
+          ~rank:1 ~dim:[| prim_num |] ~data:coef) ;
+        Ezfio.set_basis_prim_expo (Ezfio.ezfio_array_of_list
+          ~rank:1 ~dim:[| prim_num |] ~data:expo) ;
+
+
         Ezfio.set_ao_basis_ao_prim_num (Ezfio.ezfio_array_of_list
           ~rank:1 ~dim:[| ao_num |] ~data:ao_prim_num) ;
         Ezfio.set_ao_basis_ao_nucl(Ezfio.ezfio_array_of_list

src/ao_basis/NEED

@@ -1 +1,2 @@
 nuclei
+basis

src/ao_basis/aos.irp.f

@@ -6,6 +6,23 @@ BEGIN_PROVIDER [ integer, ao_prim_num_max ]
  ao_prim_num_max = maxval(ao_prim_num)
 END_PROVIDER
 
+BEGIN_PROVIDER [ integer, ao_shell, (ao_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Index of the shell to which the AO corresponds
+ END_DOC
+ integer :: i, j, k, n
+ k=0
+ do i=1,shell_num
+   n = shell_ang_mom(i)+1
+   do j=1,(n*(n+1))/2
+     k = k+1
+     ao_shell(k) = i
+   enddo
+ enddo
+
+END_PROVIDER
+
  BEGIN_PROVIDER [ double precision, ao_coef_normalized, (ao_num,ao_prim_num_max) ]
 &BEGIN_PROVIDER [ double precision, ao_coef_normalization_factor, (ao_num) ]
   implicit none
@@ -23,14 +40,15 @@ END_PROVIDER
 
   do i=1,ao_num
 
-    powA(1) = ao_power(i,1)
+    powA(1) = ao_power(i,1) +  ao_power(i,2) +  ao_power(i,3)
-    powA(2) = ao_power(i,2)
+    powA(2) = 0
-    powA(3) = ao_power(i,3)
+    powA(3) = 0
 
     ! Normalization of the primitives
     if (primitives_normalized) then
       do j=1,ao_prim_num(i)
-        call overlap_gaussian_xyz(C_A,C_A,ao_expo(i,j),ao_expo(i,j),powA,powA,overlap_x,overlap_y,overlap_z,norm,nz)
+        call overlap_gaussian_xyz(C_A,C_A,ao_expo(i,j),ao_expo(i,j), &
+           powA,powA,overlap_x,overlap_y,overlap_z,norm,nz)
         ao_coef_normalized(i,j) = ao_coef(i,j)/sqrt(norm)
       enddo
     else
@@ -39,6 +57,10 @@ END_PROVIDER
       enddo
     endif
 
+    powA(1) = ao_power(i,1)
+    powA(2) = ao_power(i,2)
+    powA(3) = ao_power(i,3)
+
     ! Normalization of the contracted basis functions
     if (ao_normalized) then
       norm = 0.d0
@@ -56,39 +78,6 @@ END_PROVIDER
 
 END_PROVIDER
 
-BEGIN_PROVIDER [ double precision, ao_coef_normalization_libint_factor, (ao_num) ]
-  implicit none
-  BEGIN_DOC
-  ! |AO| normalization for interfacing with libint
-  END_DOC
-  double precision               :: norm,overlap_x,overlap_y,overlap_z,C_A(3), c
-  integer                        :: l, powA(3), nz
-  integer                        :: i,j,k
-  nz=100
-  C_A(1) = 0.d0
-  C_A(2) = 0.d0
-  C_A(3) = 0.d0
-
-  do i=1,ao_num
-    powA(1) = ao_l(i)
-    powA(2) = 0
-    powA(3) = 0
-
-    ! Normalization of the contracted basis functions
-    norm = 0.d0
-    do j=1,ao_prim_num(i)
-     do k=1,ao_prim_num(i)
-      call overlap_gaussian_xyz(C_A,C_A,ao_expo(i,j),ao_expo(i,k),powA,powA,overlap_x,overlap_y,overlap_z,c,nz)
-      norm = norm+c*ao_coef_normalized(i,j)*ao_coef_normalized(i,k)
-     enddo
-    enddo
-    ao_coef_normalization_libint_factor(i) = ao_coef_normalization_factor(i) * sqrt(norm)
-
-  enddo
-
-END_PROVIDER
-
-
  BEGIN_PROVIDER [ double precision, ao_coef_normalized_ordered, (ao_num,ao_prim_num_max) ]
 &BEGIN_PROVIDER [ double precision, ao_expo_ordered, (ao_num,ao_prim_num_max) ]
   implicit none
@@ -231,38 +220,11 @@ END_PROVIDER
  END_DOC
  do i = 1, nucl_num
   Nucl_num_shell_Aos(i) = 0
-
   do j = 1, Nucl_N_Aos(i)
-   if(ao_l(Nucl_Aos(i,j))==0)then
+    if (ao_power(Nucl_Aos(i,j),1) == ao_l(Nucl_Aos(i,j))) then
-   ! S type function
-   Nucl_num_shell_Aos(i)+=1
-   Nucl_list_shell_Aos(i,Nucl_num_shell_Aos(i))=Nucl_Aos(i,j)
-   elseif(ao_l(Nucl_Aos(i,j))==1)then
-   ! P type function
-    if(ao_power(Nucl_Aos(i,j),1)==1)then
      Nucl_num_shell_Aos(i)+=1
      Nucl_list_shell_Aos(i,Nucl_num_shell_Aos(i))=Nucl_Aos(i,j)
     endif
-   elseif(ao_l(Nucl_Aos(i,j))==2)then
-   ! D type function
-    if(ao_power(Nucl_Aos(i,j),1)==2)then
-     Nucl_num_shell_Aos(i)+=1
-     Nucl_list_shell_Aos(i,Nucl_num_shell_Aos(i))=Nucl_Aos(i,j)
-    endif
-   elseif(ao_l(Nucl_Aos(i,j))==3)then
-   ! F type function
-    if(ao_power(Nucl_Aos(i,j),1)==3)then
-     Nucl_num_shell_Aos(i)+=1
-     Nucl_list_shell_Aos(i,Nucl_num_shell_Aos(i))=Nucl_Aos(i,j)
-    endif
-   elseif(ao_l(Nucl_Aos(i,j))==4)then
-   ! G type function
-    if(ao_power(Nucl_Aos(i,j),1)==4)then
-     Nucl_num_shell_Aos(i)+=1
-     Nucl_list_shell_Aos(i,Nucl_num_shell_Aos(i))=Nucl_Aos(i,j)
-    endif
-   endif
-
   enddo
  enddo
 

src/ao_one_e_ints/pseudopot.f90

@@ -96,8 +96,12 @@ end
 ! x=cos(theta)
 
         double precision function ylm_real(l,m,x,phi)
-        implicit double precision (a-h,o-z)
+        implicit none
-        DIMENSION PM(0:100,0:100)
+        integer :: MM, iabs_m, m, l
+        double precision :: pi, fourpi, factor, x, phi, coef
+        double precision :: xchap, ychap, zchap
+        double precision, external :: fact
+        double precision :: PM(0:100,0:100), plm
         MM=100
         pi=dacos(-1.d0)
         fourpi=4.d0*pi
@@ -1150,8 +1154,10 @@ end
 !       Output:  PM(m,n) --- Pmn(x)
 !       =====================================================
 !
-        IMPLICIT DOUBLE PRECISION (P,X)
+        implicit none
-        DIMENSION PM(0:MM,0:(N+1))
+!        IMPLICIT DOUBLE PRECISION (P,X)
+        integer :: MM, N, I, J, M
+        double precision :: PM(0:MM,0:(N+1)), X, XQ, XS
         DOUBLE PRECISION, SAVE ::  INVERSE(100) = 0.D0
         DOUBLE PRECISION       :: LS, II, JJ
         IF (INVERSE(1) == 0.d0) THEN
@@ -1202,8 +1208,9 @@ end
 !  P_l^|m|(cos(theta))  exp(i m phi)
 
       subroutine erreur(x,n,rmoy,error)
-      implicit double precision(a-h,o-z)
+      implicit none
-      dimension x(n)
+      integer :: i, n
+      double precision :: x(n), rn, rn1, error, rmoy
 ! calcul de la moyenne
       rmoy=0.d0
       do i=1,n

src/basis/EZFIO.cfg

@@ -0,0 +1,75 @@
+[basis]
+type: character*(256)
+doc: Name of the |AO| basis set
+interface: ezfio, provider
+
+[typ]
+type: character*(32)
+doc: Type of basis set. Only 'Gaussian' is supported
+interface: ezfio, provider
+
+[shell_num]
+type: integer
+doc: Number of shells
+interface: ezfio, provider
+
+[nucleus_shell_num]
+type: integer
+doc: Number of shells per nucleus
+size: (nuclei.nucl_num)
+interface: ezfio, provider
+
+[shell_normalization_factor]
+type: double precision
+doc: Normalization factor applied to the whole shell, ex $1/\sqrt{ <d_{z^2}|d_{z^2}>}$
+size: (basis.shell_num)
+interface: ezfio
+
+[shell_ang_mom]
+type: integer
+doc: Angular momentum of each shell
+size: (basis.shell_num)
+interface: ezfio, provider
+
+[shell_prim_num]
+type: integer
+doc: Number of primitives in a shell
+size: (basis.shell_num)
+interface: ezfio, provider
+
+[shell_prim_index]
+type: integer
+doc: Max number of primitives in a shell
+size: (basis.shell_num)
+interface: ezfio, provider
+
+[basis_nucleus_index]
+type: integer
+doc: Index of the nucleus on which the shell is centered
+size: (nuclei.nucl_num)
+interface: ezfio, provider
+
+[prim_normalization_factor]
+type: double precision
+doc: Normalization factor applied to each primitive
+size: (basis.prim_num)
+interface: ezfio
+
+[prim_num]
+type: integer
+doc: Total number of primitives
+interface: ezfio, provider
+
+[prim_coef]
+type: double precision
+doc: Primitive coefficients
+size: (basis.prim_num)
+interface: ezfio, provider
+
+[prim_expo]
+type: double precision
+doc: Exponents in the shell
+size: (basis.prim_num)
+interface: ezfio, provider
+
+

src/basis/NEED

@@ -0,0 +1 @@
+nuclei

src/basis/README.rst

@@ -0,0 +1,8 @@
+======
+basis
+======
+
+This module contains the basis set information, which will then be used to build the atomic orbitals.
+
+
+

src/basis/basis.irp.f

@@ -0,0 +1,119 @@
+BEGIN_PROVIDER [ double precision, shell_normalization_factor , (shell_num) ]
+  implicit none
+  BEGIN_DOC
+  ! Number of primitives per |AO|
+  END_DOC
+
+  logical                        :: has
+  PROVIDE ezfio_filename
+  if (mpi_master) then
+    if (size(shell_normalization_factor) == 0) return
+
+    call ezfio_has_basis_shell_normalization_factor(has)
+    if (has) then
+      write(6,'(A)') '.. >>>>> [ IO READ: shell_normalization_factor ] <<<<< ..'
+      call ezfio_get_basis_shell_normalization_factor(shell_normalization_factor)
+    else
+
+      double precision               :: norm,overlap_x,overlap_y,overlap_z,C_A(3), c
+      integer                        :: l, powA(3), nz
+      integer                        :: i,j,k
+      nz=100
+      C_A(1) = 0.d0
+      C_A(2) = 0.d0
+      C_A(3) = 0.d0
+
+      do i=1,shell_num
+
+        powA(1) = shell_ang_mom(i)
+        powA(2) = 0
+        powA(3) = 0
+
+        norm = 0.d0
+        do k=shell_prim_index(i),shell_prim_index(i)+shell_prim_num(i)-1
+          do j=shell_prim_index(i),shell_prim_index(i)+shell_prim_num(i)-1
+            call overlap_gaussian_xyz(C_A,C_A,prim_expo(j),prim_expo(k), &
+              powA,powA,overlap_x,overlap_y,overlap_z,c,nz)
+            norm = norm+c*prim_coef(j)*prim_coef(k) * prim_normalization_factor(j) * prim_normalization_factor(k)
+          enddo
+        enddo
+        shell_normalization_factor(i) = 1.d0/dsqrt(norm)
+      enddo
+
+    endif
+  endif
+  IRP_IF MPI_DEBUG
+  print *,  irp_here, mpi_rank
+  call MPI_BARRIER(MPI_COMM_WORLD, ierr)
+  IRP_ENDIF
+  IRP_IF MPI
+  include 'mpif.h'
+  integer                        :: ierr
+  call MPI_BCAST( shell_normalization_factor, (shell_num), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
+  if (ierr /= MPI_SUCCESS) then
+    stop 'Unable to read shell_normalization_factor with MPI'
+  endif
+  IRP_ENDIF
+
+  call write_time(6)
+
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, prim_normalization_factor , (prim_num) ]
+  implicit none
+  BEGIN_DOC
+  ! Number of primitives per |AO|
+  END_DOC
+
+  logical                        :: has
+  PROVIDE ezfio_filename
+  if (mpi_master) then
+    if (size(prim_normalization_factor) == 0) return
+
+    call ezfio_has_basis_prim_normalization_factor(has)
+    if (has) then
+      write(6,'(A)') '.. >>>>> [ IO READ: prim_normalization_factor ] <<<<< ..'
+      call ezfio_get_basis_prim_normalization_factor(prim_normalization_factor)
+    else
+
+      double precision               :: norm,overlap_x,overlap_y,overlap_z,C_A(3), c
+      integer                        :: l, powA(3), nz
+      integer                        :: i,j,k
+      nz=100
+      C_A(1) = 0.d0
+      C_A(2) = 0.d0
+      C_A(3) = 0.d0
+
+      do i=1,shell_num
+
+        powA(1) = shell_ang_mom(i)
+        powA(2) = 0
+        powA(3) = 0
+
+        do k=shell_prim_index(i),shell_prim_index(i)+shell_prim_num(i)-1
+            call overlap_gaussian_xyz(C_A,C_A,prim_expo(k),prim_expo(k), &
+              powA,powA,overlap_x,overlap_y,overlap_z,norm,nz)
+            prim_normalization_factor(k) = 1.d0/dsqrt(norm)
+        enddo
+      enddo
+
+
+    endif
+  endif
+  IRP_IF MPI_DEBUG
+  print *,  irp_here, mpi_rank
+  call MPI_BARRIER(MPI_COMM_WORLD, ierr)
+  IRP_ENDIF
+  IRP_IF MPI
+  include 'mpif.h'
+  integer                        :: ierr
+  call MPI_BCAST( prim_normalization_factor, (prim_num), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
+  if (ierr /= MPI_SUCCESS) then
+    stop 'Unable to read prim_normalization_factor with MPI'
+  endif
+  IRP_ENDIF
+
+  call write_time(6)
+
+END_PROVIDER

src/bitmask/bitmasks.irp.f

@@ -27,9 +27,7 @@ BEGIN_PROVIDER [ integer(bit_kind), full_ijkl_bitmask, (N_int) ]
     full_ijkl_bitmask(j) = 0_bit_kind
     do i=0,bit_kind_size-1
       k=k+1
-      if (mo_class(k) /= 'Deleted') then
       full_ijkl_bitmask(j) = ibset(full_ijkl_bitmask(j),i)
-      endif
       if (k == mo_num) exit
     enddo
   enddo

src/cipsi/pt2_stoch_routines.irp.f

@@ -286,7 +286,7 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in)
       call write_int(6,nproc_target,'Number of threads for PT2')
       call write_double(6,mem,'Memory (Gb)')
 
-      call omp_set_nested(.false.)
+      call omp_set_max_active_levels(1)
 
 
       print '(A)', '========== ======================= ===================== ===================== ==========='
@@ -313,6 +313,7 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in)
       endif
       !$OMP END PARALLEL
       call end_parallel_job(zmq_to_qp_run_socket, zmq_socket_pull, 'pt2')
+      call omp_set_max_active_levels(8)
 
       print '(A)', '========== ======================= ===================== ===================== ==========='
 

src/cipsi/selection.irp.f

@@ -253,12 +253,7 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
   deallocate(exc_degree)
   nmax=k-1
 
-  allocate(iorder(nmax))
+  call isort_noidx(indices,nmax)
-  do i=1,nmax
-    iorder(i) = i
-  enddo
-  call isort(indices,iorder,nmax)
-  deallocate(iorder)
 
   ! Start with 32 elements. Size will double along with the filtering.
   allocate(preinteresting(0:32), prefullinteresting(0:32),     &
@@ -749,7 +744,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
 
       double precision :: eigvalues(N_states+1)
       double precision :: work(1+6*(N_states+1)+2*(N_states+1)**2)
-      integer :: iwork(3+5*(N_states+1)), info, k 
+      integer :: info, k , iwork(N_states+1)
 
       if (do_diag) then
         double precision :: pt2_matrix(N_states+1,N_states+1)
@@ -761,8 +756,8 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
           pt2_matrix(N_states+1,istate) = mat(istate,p1,p2)
         enddo
 
-        call DSYEVD( 'V', 'U', N_states+1, pt2_matrix, N_states+1, eigvalues, &
+        call DSYEV( 'V', 'U', N_states+1, pt2_matrix, N_states+1, eigvalues, &
-                     work, size(work), iwork, size(iwork), info )
+                     work, size(work), info )
         if (info /= 0) then
           print *, 'error in '//irp_here
           stop -1
@@ -770,7 +765,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
         pt2_matrix = dabs(pt2_matrix)
         iwork(1:N_states+1) = maxloc(pt2_matrix,DIM=1)
         do k=1,N_states
-          e_pert(iwork(k)) = eigvalues(k) - E0(iwork(k))
+          e_pert(k) = eigvalues(iwork(k)) - E0(k)
         enddo
       endif
 

src/csf/sigma_vector.irp.f

@@ -1,3 +1,9 @@
+ real*8 function logabsgamma(x)
+ implicit none
+ real*8, intent(in) :: x
+ logabsgamma = log(abs(gamma(x)))
+ end function logabsgamma
+  
   BEGIN_PROVIDER [ integer, NSOMOMax]
  &BEGIN_PROVIDER [ integer, NCSFMax]
  &BEGIN_PROVIDER [ integer*8, NMO]
@@ -22,33 +28,65 @@
   integer NSOMO
   integer dimcsfpercfg
   integer detDimperBF
-  real*8 :: coeff
+  real*8 :: coeff, binom1, binom2
   integer MS
   integer ncfgpersomo
+  real*8, external :: logabsgamma
   detDimperBF = 0
   MS = elec_alpha_num-elec_beta_num
   ! number of cfgs = number of dets for 0 somos
-  n_CSF = cfg_seniority_index(0)-1
+  n_CSF = 0
   ncfgprev = cfg_seniority_index(0)
-  do i = 0-iand(MS,1)+2, NSOMOMax,2
+  ncfgpersomo = ncfgprev
-     if(cfg_seniority_index(i) .EQ. -1)then
+  do i = iand(MS,1), NSOMOMax-2,2
+    if(cfg_seniority_index(i) .EQ. -1) then
+      cycle
+    endif
+    if(cfg_seniority_index(i+2) .EQ. -1) then
       ncfgpersomo = N_configuration + 1
     else
-        ncfgpersomo = cfg_seniority_index(i)
+      if(cfg_seniority_index(i+2) > ncfgpersomo) then
+          ncfgpersomo = cfg_seniority_index(i+2)
+      else
+        k = 0
+        do while(cfg_seniority_index(i+2+k) < ncfgpersomo)
+          k = k + 2
+          ncfgpersomo = cfg_seniority_index(i+2+k)
+        enddo
+      endif
     endif
     ncfg = ncfgpersomo - ncfgprev
-  !detDimperBF = max(1,nint((binom(i,(i+1)/2))))
+    if(iand(MS,1) .EQ. 0) then
-  if (i > 2) then
+      !dimcsfpercfg = max(1,nint((binom(i,i/2)-binom(i,i/2+1))))
-    dimcsfpercfg = max(1,nint((binom(i-2,(i-2+1)/2)-binom(i-2,((i-2+1)/2)+1))))
+      binom1 = dexp(logabsgamma(1.0d0*(i+1))                            &
+                  - logabsgamma(1.0d0*((i/2)+1))                        &
+                  - logabsgamma(1.0d0*(i-((i/2))+1)));
+      binom2 = dexp(logabsgamma(1.0d0*(i+1))                            &
+                  - logabsgamma(1.0d0*(((i/2)+1)+1))                    &
+                  - logabsgamma(1.0d0*(i-((i/2)+1)+1)));
+      dimcsfpercfg = max(1,nint(binom1 - binom2))
     else
-    dimcsfpercfg = 1
+      !dimcsfpercfg = max(1,nint((binom(i,(i+1)/2)-binom(i,(i+3)/2))))
+      binom1 = dexp(logabsgamma(1.0d0*(i+1))                            &
+                  - logabsgamma(1.0d0*(((i+1)/2)+1))                    &
+                  - logabsgamma(1.0d0*(i-(((i+1)/2))+1)));
+      binom2 = dexp(logabsgamma(1.0d0*(i+1))                            &
+                  - logabsgamma(1.0d0*((((i+3)/2)+1)+1))                &
+                  - logabsgamma(1.0d0*(i-(((i+3)/2)+1)+1)));
+      dimcsfpercfg = max(1,nint(binom1 - binom2))
     endif
     n_CSF += ncfg * dimcsfpercfg
-  !if(cfg_seniority_index(i+2) == -1) EXIT
+    if(cfg_seniority_index(i+2) > ncfgprev) then
-  !if(detDimperBF > maxDetDimPerBF) maxDetDimPerBF = detDimperBF
+      ncfgprev = cfg_seniority_index(i+2)
-  ncfgprev = cfg_seniority_index(i)
+    else
+      k = 0
+      do while(cfg_seniority_index(i+2+k) < ncfgprev)
+        k = k + 2
+        ncfgprev = cfg_seniority_index(i+2+k)
       enddo
-  END_PROVIDER
+    endif
+  enddo
+END_PROVIDER
 
 
 subroutine get_phase_qp_to_cfg(Ialpha, Ibeta, phaseout)

src/davidson/diagonalization_hcsf_dressed.irp.f

@@ -197,6 +197,7 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N
   call write_int(6,N_st,'Number of states')
   call write_int(6,N_st_diag,'Number of states in diagonalization')
   call write_int(6,sze,'Number of determinants')
+  call write_int(6,sze_csf,'Number of CSFs')
   call write_int(6,nproc_target,'Number of threads for diagonalization')
   call write_double(6, r1, 'Memory(Gb)')
   if (disk_based) then

src/dressing/run_dress_slave.irp.f

@@ -72,7 +72,7 @@ subroutine run_dress_slave(thread,iproce,energy)
   provide psi_energy
   ending = dress_N_cp+1
   ntask_tbd = 0
-  call omp_set_nested(.true.)
+  call omp_set_max_active_levels(8)
 
   !$OMP PARALLEL DEFAULT(SHARED) &
   !$OMP PRIVATE(interesting, breve_delta_m, task_id) &
@@ -84,7 +84,7 @@ subroutine run_dress_slave(thread,iproce,energy)
   zmq_socket_push      = new_zmq_push_socket(thread)
   integer, external :: connect_to_taskserver
   !$OMP CRITICAL
-  call omp_set_nested(.false.)
+  call omp_set_max_active_levels(1)
   if (connect_to_taskserver(zmq_to_qp_run_socket,worker_id,thread) == -1) then
     print *,  irp_here, ': Unable to connect to task server'
     stop -1
@@ -296,7 +296,7 @@ subroutine run_dress_slave(thread,iproce,energy)
   !$OMP END CRITICAL
 
   !$OMP END PARALLEL
-  call omp_set_nested(.false.)
+  call omp_set_max_active_levels(1)
 !  do i=0,dress_N_cp+1
 !    call omp_destroy_lock(lck_sto(i))
 !  end do

src/mo_one_e_ints/EZFIO.cfg

@@ -47,10 +47,3 @@ type: Disk_access
 doc: Read/Write |MO| one-electron integrals from/to disk [ Write | Read | None ]
 interface: ezfio,provider,ocaml
 default: None
-
-
-[restore_symm]
-type: logical
-doc: If true, try to find symmetry in the MO coefficient matrices
-interface: ezfio,provider,ocaml
-default: True

src/two_body_rdm/two_e_dm_mo.irp.f

@@ -5,18 +5,14 @@ BEGIN_PROVIDER [double precision, two_e_dm_mo, (mo_num,mo_num,mo_num,mo_num)]
    !
    ! <Psi| a^{\dagger}_{i \beta} a^{\dagger}_{j \beta} a_{l \beta} a_{k \beta} |Psi>
    !
-   ! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO ALL OCCUPIED ORBITALS : core, inactive and active
+   ! where the indices (i,j,k,l) belong to all MOs.
    !
-   ! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{elec} * (N_{elec} - 1)/2
+   ! The normalization (i.e. sum of diagonal elements) is set to $N_{elec} * (N_{elec} - 1)/2$
    !
-   !  !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
+   !  !!!!! WARNING !!!!! IF "no_core_density" then all elements involving at least one CORE MO are set to zero
+   ! The state-averaged two-electron energy :
    !
-   !  !!!!! WARNING !!!!! IF "no_core_density" then all elements involving at least one CORE MO is set to zero
+   !   \sum_{i,j,k,l = 1, mo_num} two_e_dm_mo(i,j,k,l) * < ii jj | kk ll >
-   ! The two-electron energy of each state can be computed as:
-   !
-   !   \sum_{i,j,k,l = 1, n_core_inact_act_orb} two_e_dm_mo(i,j,k,l,istate) * < ii jj | kk ll >
-   !
-   !   with ii = list_core_inact_act(i), jj = list_core_inact_act(j), kk = list_core_inact_act(k), ll = list_core_inact_act(l)
    END_DOC
    two_e_dm_mo = 0.d0
    integer                        :: i,j,k,l,iorb,jorb,korb,lorb,istate

src/utils/EZFIO.cfg

@@ -0,0 +1,5 @@
+[restore_symm]
+type: logical
+doc: If true, try to find symmetry in the MO coefficient matrices
+interface: ezfio,provider,ocaml
+default: False

src/utils/intel.f90

@@ -0,0 +1,173 @@
+module intel
+  use, intrinsic :: iso_c_binding
+  interface
+    subroutine ippsSortAscend_32s_I(pSrc, len) bind(C, name='ippsSortAscend_32s_I')
+     use  iso_c_binding
+     integer, intent(in), value :: len
+     integer, intent(inout) :: pSrc(len)
+    end
+  end interface
+  interface
+    subroutine ippsSortAscend_32f_I(pSrc, len) bind(C, name='ippsSortAscend_32f_I')
+     use  iso_c_binding
+     integer, intent(in), value :: len
+     real, intent(inout) :: pSrc(len)
+    end
+  end interface
+  interface
+    subroutine ippsSortAscend_64s_I(pSrc, len) bind(C, name='ippsSortAscend_64s_I')
+     use  iso_c_binding
+     integer, intent(in), value :: len
+     integer*8, intent(inout) :: pSrc(len)
+    end
+  end interface
+  interface
+    subroutine ippsSortAscend_64f_I(pSrc, len) bind(C, name='ippsSortAscend_64f_I')
+     use  iso_c_binding
+     integer, intent(in), value :: len
+     double precision, intent(inout) :: pSrc(len)
+    end
+  end interface
+
+  interface
+    subroutine ippsSortRadixIndexGetBufferSize(len, dataType, pBufSize) bind(C, name='ippsSortRadixIndexGetBufferSize')
+     use  iso_c_binding
+     integer, intent(in), value :: len
+     integer, intent(in), value :: dataType
+     integer, intent(out) :: pBufSize
+    end
+  end interface
+
+  interface
+    subroutine ippsSortRadixAscend_16s_I(pSrc, len, pTmp) bind(C, name='ippsSortRadixAscend_16s_I')
+     use  iso_c_binding
+     integer, intent(in), value :: len
+     integer*2, intent(inout) :: pSrc(len)
+     character, intent(inout) :: pTmp(len)
+    end
+  end interface
+  interface
+    subroutine ippsSortRadixAscend_32s_I(pSrc, len, pTmp) bind(C, name='ippsSortRadixAscend_32s_I')
+     use  iso_c_binding
+     integer, intent(in), value :: len
+     integer, intent(inout) :: pSrc(len)
+     character, intent(inout) :: pTmp(len)
+    end
+  end interface
+  interface
+    subroutine ippsSortRadixAscend_32f_I(pSrc, len, pTmp) bind(C, name='ippsSortRadixAscend_32f_I')
+     use  iso_c_binding
+     integer, intent(in), value :: len
+     real, intent(inout) :: pSrc(len)
+     character, intent(inout) :: pTmp(len)
+    end
+  end interface
+  interface
+    subroutine ippsSortRadixAscend_64s_I(pSrc, len, pTmp) bind(C, name='ippsSortRadixAscend_64s_I')
+     use  iso_c_binding
+     integer, intent(in), value :: len
+     integer*8, intent(inout) :: pSrc(len)
+     character, intent(inout) :: pTmp(len)
+    end
+  end interface
+  interface
+    subroutine ippsSortRadixAscend_64f_I(pSrc, len, pTmp) bind(C, name='ippsSortRadixAscend_64f_I')
+     use  iso_c_binding
+     integer, intent(in), value :: len
+     double precision, intent(inout) :: pSrc(len)
+     character, intent(inout) :: pTmp(len)
+    end
+  end interface
+
+  interface
+    subroutine ippsSortRadixIndexAscend_16s(pSrc, srcStrideBytes, pDstIndx, len, pTmpIndx) bind(C, name='ippsSortRadixIndexAscend_16s')
+     use  iso_c_binding
+     integer, intent(in), value :: len
+     integer*2, intent(inout) :: pSrc(len)
+     integer, intent(in), value :: srcStrideBytes
+     integer, intent(inout) :: pDstIndx(len)
+     character, intent(inout) :: pTmpIndx(len)
+    end
+  end interface
+  interface
+    subroutine ippsSortRadixIndexAscend_32s(pSrc, srcStrideBytes, pDstIndx, len, pTmpIndx) bind(C, name='ippsSortRadixIndexAscend_32s')
+     use  iso_c_binding
+     integer, intent(in), value :: len
+     integer, intent(inout) :: pSrc(len)
+     integer, intent(in), value :: srcStrideBytes
+     integer, intent(inout) :: pDstIndx(len)
+     character, intent(inout) :: pTmpIndx(len)
+    end
+  end interface
+  interface
+    subroutine ippsSortRadixIndexAscend_32f(pSrc, srcStrideBytes, pDstIndx, len, pTmpIndx) bind(C,name='ippsSortRadixIndexAscend_32f')
+     use  iso_c_binding
+     integer, intent(in), value :: len
+     real   , intent(inout) :: pSrc(len)
+     integer, intent(in), value :: srcStrideBytes
+     integer, intent(inout) :: pDstIndx(len)
+     character, intent(inout) :: pTmpIndx(len)
+    end
+  end interface
+  interface
+    subroutine ippsSortRadixIndexAscend_64s(pSrc, srcStrideBytes, pDstIndx, len, pTmpIndx) bind(C, name='ippsSortRadixIndexAscend_64s')
+     use  iso_c_binding
+     integer, intent(in), value :: len
+     integer*8, intent(inout) :: pSrc(len)
+     integer, intent(in), value :: srcStrideBytes
+     integer, intent(inout) :: pDstIndx(len)
+     character, intent(inout) :: pTmpIndx(len)
+    end
+  end interface
+  interface
+    subroutine ippsSortRadixIndexAscend_64f(pSrc, srcStrideBytes, pDstIndx, len, pTmpIndx) bind(C,name='ippsSortRadixIndexAscend_64f')
+     use  iso_c_binding
+     integer, intent(in), value :: len
+     real*8 , intent(inout) :: pSrc(len)
+     integer, intent(in), value :: srcStrideBytes
+     integer, intent(inout) :: pDstIndx(len)
+     character, intent(inout) :: pTmpIndx(len)
+    end
+  end interface
+
+  interface
+    subroutine ippsSortIndexAscend_32f_I(pSrcDst, pDstIndx, len) bind(C,name='ippsSortIndexAscend_32f_I')
+     use  iso_c_binding
+     real(4), intent(in) :: pSrcDst(*)
+     integer(4), intent(inout) :: pDstIndx(*)
+     integer(4), intent(in), value :: len
+    end
+  end interface
+  interface
+    subroutine ippsSortIndexAscend_32s_I(pSrcDst, pDstIndx, len) bind(C,name='ippsSortIndexAscend_32s_I')
+     use  iso_c_binding
+     integer(4), intent(in) :: pSrcDst(*)
+     integer(4), intent(inout) :: pDstIndx(*)
+     integer(4), intent(in), value :: len
+    end
+  end interface
+  interface
+    subroutine ippsSortIndexAscend_64f_I(pSrcDst, pDstIndx, len) bind(C,name='ippsSortIndexAscend_64f_I')
+     use  iso_c_binding
+     real(8), intent(in) :: pSrcDst(*)
+     integer(4), intent(inout) :: pDstIndx(*)
+     integer(4), intent(in), value :: len
+    end
+  end interface
+  interface
+    subroutine ippsSortIndexAscend_64s_I(pSrcDst, pDstIndx, len) bind(C,name='ippsSortIndexAscend_64s_I')
+     use  iso_c_binding
+     integer(8), intent(in) :: pSrcDst(*)
+     integer(4), intent(inout) :: pDstIndx(*)
+     integer(4), intent(in), value :: len
+    end
+  end interface
+  interface
+    subroutine ippsSortIndexAscend_16s_I(pSrcDst, pDstIndx, len) bind(C,name='ippsSortIndexAscend_16s_I')
+     use  iso_c_binding
+     integer(2), intent(in) :: pSrcDst(*)
+     integer(4), intent(inout) :: pDstIndx(*)
+     integer(4), intent(in), value :: len
+    end
+  end interface
+end module

src/utils/sort.irp.f

@@ -38,15 +38,7 @@ BEGIN_TEMPLATE
   $type,intent(inout)            :: x(isize)
   integer,intent(inout)          :: iorder(isize)
   integer, external              :: omp_get_num_threads
-  if (omp_get_num_threads() == 1) then
-    !$OMP PARALLEL DEFAULT(SHARED)
-    !$OMP SINGLE
   call rec_$X_quicksort(x,iorder,isize,1,isize,nproc)
-    !$OMP END SINGLE
-    !$OMP END PARALLEL
-  else
-    call rec_$X_quicksort(x,iorder,isize,1,isize,nproc)
-  endif
  end
 
  recursive subroutine rec_$X_quicksort(x, iorder, isize, first, last, level)
@@ -89,16 +81,11 @@ BEGIN_TEMPLATE
     endif
   else
     if (first < i-1) then
-      !$OMP TASK DEFAULT(SHARED) FIRSTPRIVATE(isize,first,i,level)
       call rec_$X_quicksort(x, iorder, isize, first, i-1,level/2)
-      !$OMP END TASK
     endif
     if (j+1 < last) then
-      !$OMP TASK DEFAULT(SHARED) FIRSTPRIVATE(isize,last,j,level)
       call rec_$X_quicksort(x, iorder, isize, j+1, last,level/2)
-      !$OMP END TASK
     endif
-    !$OMP TASKWAIT
   endif
  end
 
@@ -262,7 +249,60 @@ SUBST [ X, type ]
  i2 ; integer*2 ;;
 END_TEMPLATE
 
+
+!---------------------- INTEL
+IRP_IF INTEL
+
 BEGIN_TEMPLATE
+ subroutine $Xsort(x,iorder,isize)
+  use intel
+  implicit none
+  BEGIN_DOC
+  ! Sort array x(isize).
+  ! iorder in input should be (1,2,3,...,isize), and in output
+  ! contains the new order of the elements.
+  END_DOC
+  integer,intent(in)             :: isize
+  $type,intent(inout)            :: x(isize)
+  integer,intent(inout)          :: iorder(isize)
+  integer                        :: n
+  character, allocatable         :: tmp(:)
+  if (isize < 2) return
+  call ippsSortRadixIndexGetBufferSize(isize, $ippsz, n)
+  allocate(tmp(n))
+  call ippsSortRadixIndexAscend_$ityp(x, $n, iorder, isize, tmp)
+  deallocate(tmp)
+  iorder(1:isize) = iorder(1:isize)+1
+  call $Xset_order(x,iorder,isize)
+ end
+
+ subroutine $Xsort_noidx(x,isize)
+  use intel
+  implicit none
+  BEGIN_DOC
+  ! Sort array x(isize).
+  ! iorder in input should be (1,2,3,...,isize), and in output
+  ! contains the new order of the elements.
+  END_DOC
+  integer,intent(in)             :: isize
+  $type,intent(inout)            :: x(isize)
+  integer                        :: n
+  character, allocatable         :: tmp(:)
+  if (isize < 2) return
+  call ippsSortRadixIndexGetBufferSize(isize, $ippsz, n)
+  allocate(tmp(n))
+  call ippsSortRadixAscend_$ityp_I(x, isize, tmp)
+  deallocate(tmp)
+ end
+
+SUBST [ X, type, ityp, n, ippsz ]
+   ; real ; 32f ; 4 ; 13 ;;
+ i ; integer ; 32s ; 4 ; 11 ;;
+ i2 ; integer*2 ; 16s ; 2 ; 7 ;;
+END_TEMPLATE
+
+BEGIN_TEMPLATE
+
  subroutine $Xsort(x,iorder,isize)
   implicit none
   BEGIN_DOC
@@ -289,12 +329,12 @@ BEGIN_TEMPLATE
   endif
  end subroutine $Xsort
 
-SUBST [ X, type, Y ]
+SUBST [ X, type ]
-   ; real ; i ;;
+ d ; double precision ;;
- d ; double precision ; i8 ;;
 END_TEMPLATE
 
 BEGIN_TEMPLATE
+
  subroutine $Xsort(x,iorder,isize)
   implicit none
   BEGIN_DOC
@@ -306,8 +346,112 @@ BEGIN_TEMPLATE
   $type,intent(inout)            :: x(isize)
   integer,intent(inout)          :: iorder(isize)
   integer                        :: n
-!  call $Xradix_sort(x,iorder,isize,-1)
+  if (isize < 2) then
+    return
+  endif
+  call sorted_$Xnumber(x,isize,n)
+  if (isize == n) then
+    return
+  endif
+  if ( isize < 32) then
+    call insertion_$Xsort(x,iorder,isize)
+  else
+    call $Xradix_sort(x,iorder,isize,-1)
+  endif
+ end subroutine $Xsort
+
+SUBST [ X, type ]
+ i8 ; integer*8 ;;
+END_TEMPLATE
+
+!---------------------- END INTEL
+IRP_ELSE
+!---------------------- NON-INTEL
+BEGIN_TEMPLATE
+
+ subroutine $Xsort_noidx(x,isize)
+  implicit none
+  BEGIN_DOC
+  ! Sort array x(isize).
+  END_DOC
+  integer,intent(in)             :: isize
+  $type,intent(inout)            :: x(isize)
+  integer, allocatable           :: iorder(:)
+  integer                        :: i
+  allocate(iorder(isize))
+  do i=1,isize
+   iorder(i)=i
+  enddo
+  call $Xsort(x,iorder,isize)
+  deallocate(iorder)
+ end subroutine $Xsort_noidx
+
+SUBST [ X, type ]
+   ; real ;;
+ d ; double precision ;;
+ i ; integer ;;
+ i8 ; integer*8 ;;
+ i2 ; integer*2 ;;
+END_TEMPLATE
+
+BEGIN_TEMPLATE
+
+ subroutine $Xsort(x,iorder,isize)
+  implicit none
+  BEGIN_DOC
+  ! Sort array x(isize).
+  ! iorder in input should be (1,2,3,...,isize), and in output
+  ! contains the new order of the elements.
+  END_DOC
+  integer,intent(in)             :: isize
+  $type,intent(inout)            :: x(isize)
+  integer,intent(inout)          :: iorder(isize)
+  integer                        :: n
+  if (isize < 2) then
+    return
+  endif
+!  call sorted_$Xnumber(x,isize,n)
+!  if (isize == n) then
+!    return
+!  endif
+  if ( isize < 32) then
+    call insertion_$Xsort(x,iorder,isize)
+  else
+!    call heap_$Xsort(x,iorder,isize)
     call quick_$Xsort(x,iorder,isize)
+  endif
+ end subroutine $Xsort
+
+SUBST [ X, type ]
+   ; real ;;
+ d ; double precision ;;
+END_TEMPLATE
+
+BEGIN_TEMPLATE
+
+ subroutine $Xsort(x,iorder,isize)
+  implicit none
+  BEGIN_DOC
+  ! Sort array x(isize).
+  ! iorder in input should be (1,2,3,...,isize), and in output
+  ! contains the new order of the elements.
+  END_DOC
+  integer,intent(in)             :: isize
+  $type,intent(inout)            :: x(isize)
+  integer,intent(inout)          :: iorder(isize)
+  integer                        :: n
+  if (isize < 2) then
+    return
+  endif
+  call sorted_$Xnumber(x,isize,n)
+  if (isize == n) then
+    return
+  endif
+  if ( isize < 32) then
+    call insertion_$Xsort(x,iorder,isize)
+  else
+    call $Xradix_sort(x,iorder,isize,-1)
+  endif
  end subroutine $Xsort
 
 SUBST [ X, type ]
@@ -316,6 +460,11 @@ SUBST [ X, type ]
  i2 ; integer*2 ;;
 END_TEMPLATE
 
+IRP_ENDIF
+!---------------------- END NON-INTEL
+
+
+
 BEGIN_TEMPLATE
  subroutine $Xset_order(x,iorder,isize)
   implicit none
@@ -413,10 +562,12 @@ SUBST [ X, type ]
  i2; integer*2 ;;
 END_TEMPLATE
 
+
 BEGIN_TEMPLATE
 
- recursive subroutine $Xradix_sort$big(x,iorder,isize,iradix)
+recursive subroutine $Xradix_sort$big(x,iorder,isize,iradix)
   implicit none
+
   BEGIN_DOC
   ! Sort integer array x(isize) using the radix sort algorithm.
   ! iorder in input should be (1,2,3,...,isize), and in output
@@ -552,24 +703,14 @@ BEGIN_TEMPLATE
     endif
 
 
-!   !$OMP PARALLEL DEFAULT(SHARED) if (isize > 1000000)
-!   !$OMP SINGLE
     if (i3>1_$int_type) then
-!     !$OMP TASK FIRSTPRIVATE(iradix_new,i3) SHARED(x,iorder) if(i3 > 1000000)
       call $Xradix_sort$big(x,iorder,i3,iradix_new-1)
-!     !$OMP END TASK
     endif
 
     if (isize-i3>1_$int_type) then
-!     !$OMP TASK FIRSTPRIVATE(iradix_new,i3) SHARED(x,iorder) if(isize-i3 > 1000000)
       call $Xradix_sort$big(x(i3+1_$int_type),iorder(i3+1_$int_type),isize-i3,iradix_new-1)
-!     !$OMP END TASK
     endif
 
-!   !$OMP TASKWAIT
-!   !$OMP END SINGLE
-!   !$OMP END PARALLEL
-
     return
   endif
 
@@ -624,16 +765,11 @@ BEGIN_TEMPLATE
 
 
   if (i1>1_$int_type) then
-    !$OMP TASK FIRSTPRIVATE(i0,iradix,i1) SHARED(x,iorder) if(i1 >1000000)
     call $Xradix_sort$big(x(i0+1_$int_type),iorder(i0+1_$int_type),i1,iradix-1)
-    !$OMP END TASK
   endif
   if (i0>1) then
-    !$OMP TASK FIRSTPRIVATE(i0,iradix) SHARED(x,iorder) if(i0 >1000000)
     call $Xradix_sort$big(x,iorder,i0,iradix-1)
-    !$OMP END TASK
   endif
-  !$OMP TASKWAIT
 
  end
 
@@ -646,3 +782,4 @@ SUBST [ X, type, integer_size, is_big, big, int_type ]
 END_TEMPLATE
 
 
+

src/zmq/utils.irp.f

@@ -127,7 +127,9 @@ function zmq_port(ishift)
   END_DOC
   integer, intent(in)            :: ishift
   character*(8)                  :: zmq_port
+  !$OMP CRITICAL(write)
   write(zmq_port,'(I8)') zmq_port_start+ishift
+  !$OMP END CRITICAL(write)
   zmq_port = adjustl(trim(zmq_port))
 end
 
@@ -518,7 +520,9 @@ subroutine new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,name_in)
 
   zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
   zmq_socket_pull      = new_zmq_pull_socket ()
+  !$OMP CRITICAL(write)
   write(name,'(A,I8.8)') trim(name_in)//'.', icount
+  !$OMP END CRITICAL(write)
   sze = len(trim(name))
   zmq_state = trim(name)
   call lowercase(name,sze)
@@ -582,7 +586,9 @@ subroutine end_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,name_in)
   integer, save                  :: icount=0
 
   icount = icount+1
+  !$OMP CRITICAL(write)
   write(name,'(A,I8.8)') trim(name_in)//'.', icount
+  !$OMP END CRITICAL(write)
   sze = len(trim(name))
   call lowercase(name,sze)
   if (name /= zmq_state) then
@@ -704,7 +710,9 @@ integer function disconnect_from_taskserver_state(zmq_to_qp_run_socket, worker_i
 
   disconnect_from_taskserver_state = -1
 
+  !$OMP CRITICAL(write)
   write(message,*) 'disconnect '//trim(state), worker_id
+  !$OMP END CRITICAL(write)
 
   sze = min(510,len(trim(message)))
   rc = f77_zmq_send(zmq_to_qp_run_socket, trim(message), sze, 0)
@@ -781,7 +789,9 @@ integer function zmq_abort(zmq_to_qp_run_socket)
   character*(512)                :: message
   zmq_abort = 0
 
+  !$OMP CRITICAL(write)
   write(message,*) 'abort '
+  !$OMP END CRITICAL(write)
 
 
   sze = len(trim(message))
@@ -823,7 +833,9 @@ integer function task_done_to_taskserver(zmq_to_qp_run_socket, worker_id, task_i
 
   task_done_to_taskserver = 0
 
+  !$OMP CRITICAL(write)
   write(message,*) 'task_done '//trim(zmq_state), worker_id, task_id
+  !$OMP END CRITICAL(write)
 
   sze = len(trim(message))
   rc = f77_zmq_send(zmq_to_qp_run_socket, trim(message), sze, 0)
@@ -856,9 +868,11 @@ integer function tasks_done_to_taskserver(zmq_to_qp_run_socket, worker_id, task_
 
   tasks_done_to_taskserver = 0
 
+  !$OMP CRITICAL(write)
   allocate(character(LEN=64+n_tasks*12) :: message)
   write(fmt,*) '(A,X,A,I10,X,', n_tasks, '(I11,1X))'
   write(message,*) 'task_done '//trim(zmq_state), worker_id, (task_id(k), k=1,n_tasks)
+  !$OMP END CRITICAL(write)
 
   sze = len(trim(message))
   rc = f77_zmq_send(zmq_to_qp_run_socket, trim(message), sze, 0)
@@ -900,7 +914,9 @@ integer function get_task_from_taskserver(zmq_to_qp_run_socket,worker_id,task_id
 
   get_task_from_taskserver = 0
 
+  !$OMP CRITICAL(write)
   write(message,*) 'get_task '//trim(zmq_state), worker_id
+  !$OMP END CRITICAL(write)
 
   sze = len(trim(message))
   rc = f77_zmq_send(zmq_to_qp_run_socket, message, sze, 0)
@@ -961,7 +977,9 @@ integer function get_tasks_from_taskserver(zmq_to_qp_run_socket,worker_id,task_i
 
   get_tasks_from_taskserver = 0
 
+  !$OMP CRITICAL(write)
   write(message,'(A,A,X,I10,I10)') 'get_tasks ', trim(zmq_state), worker_id, n_tasks
+  !$OMP END CRITICAL(write)
 
   sze = len(trim(message))
   rc = f77_zmq_send(zmq_to_qp_run_socket, message, sze, 0)
@@ -1061,7 +1079,9 @@ integer function zmq_delete_task(zmq_to_qp_run_socket,zmq_socket_pull,task_id,mo
 
   zmq_delete_task = 0
 
+  !$OMP CRITICAL(write)
   write(message,*) 'del_task ', zmq_state, task_id
+  !$OMP END CRITICAL(write)
   rc = f77_zmq_send(zmq_to_qp_run_socket,trim(message),len(trim(message)),0)
   if (rc /= len(trim(message))) then
     zmq_delete_task = -1
@@ -1101,7 +1121,9 @@ integer function zmq_delete_task_async_send(zmq_to_qp_run_socket,task_id,sending
   endif
   zmq_delete_task_async_send = 0
 
+  !$OMP CRITICAL(write)
   write(message,*) 'del_task ', zmq_state, task_id
+  !$OMP END CRITICAL(write)
   rc = f77_zmq_send(zmq_to_qp_run_socket,trim(message),len(trim(message)),0)
   if (rc /= len(trim(message))) then
     zmq_delete_task_async_send = -1
@@ -1159,8 +1181,10 @@ integer function zmq_delete_tasks(zmq_to_qp_run_socket,zmq_socket_pull,task_id,n
 
   allocate(character(LEN=64+n_tasks*12) :: message)
 
+  !$OMP CRITICAL(write)
   write(fmt,*) '(A,1X,A,1X,', n_tasks, '(I11,1X))'
   write(message,*) 'del_task '//trim(zmq_state), (task_id(k), k=1,n_tasks)
+  !$OMP END CRITICAL(write)
 
 
   rc = f77_zmq_send(zmq_to_qp_run_socket,trim(message),len(trim(message)),0)
@@ -1206,8 +1230,10 @@ integer function zmq_delete_tasks_async_send(zmq_to_qp_run_socket,task_id,n_task
 
   allocate(character(LEN=64+n_tasks*12) :: message)
 
+  !$OMP CRITICAL(write)
   write(fmt,*) '(A,1X,A,1X,', n_tasks, '(I11,1X))'
   write(message,*) 'del_task '//trim(zmq_state), (task_id(k), k=1,n_tasks)
+  !$OMP END CRITICAL(write)
 
 
   rc = f77_zmq_send(zmq_to_qp_run_socket,trim(message),len(trim(message)),0)