diff --git a/doc/reference/h5structure.rst b/doc/reference/h5structure.rst
index e38ec01d..e9622a62 100644
--- a/doc/reference/h5structure.rst
+++ b/doc/reference/h5structure.rst
@@ -25,19 +25,19 @@ calculation. The default name of this group is `dft_input`. Its contents are
 =================  ======================================================================  =====================================================================================
 Name               Type                                                                    Meaning
 =================  ======================================================================  =====================================================================================
-energy_unit        numpy.float                                                             Unit of energy used for the calculation.
-n_k                numpy.int                                                               Number of k-points used for the BZ integration.
-k_dep_projection   numpy.int                                                               1 if the dimension of the projection operators depend on the k-point,
+energy_unit              float                                                             Unit of energy used for the calculation.
+n_k                      int                                                               Number of k-points used for the BZ integration.
+k_dep_projection         int                                                               1 if the dimension of the projection operators depend on the k-point,
                                                                                            0 otherwise.
-SP                 numpy.int                                                               1 for spin-polarised Hamiltonian, 0 for paramagnetic Hamiltonian.
-SO                 numpy.int                                                               1 if spin-orbit interaction is included, 0 otherwise.
-charge_below       numpy.float                                                             Number of electrons in the crystal below the correlated orbitals.
+SP                       int                                                               1 for spin-polarised Hamiltonian, 0 for paramagnetic Hamiltonian.
+SO                       int                                                               1 if spin-orbit interaction is included, 0 otherwise.
+charge_below             float                                                             Number of electrons in the crystal below the correlated orbitals.
                                                                                            Note that this is for compatibility with dmftproj, otherwise set to 0
-density_required   numpy.float                                                             Required total electron density. Needed to determine the chemical potential.
+density_required         float                                                             Required total electron density. Needed to determine the chemical potential.
                                                                                            The density in the projection window is then `density_required`-`charge_below`.
-symm_op            numpy.int                                                               1 if symmetry operations are used for the BZ sums,
+symm_op                  int                                                               1 if symmetry operations are used for the BZ sums,
                                                                                            0 if all k-points are directly included in the input.
-n_shells           numpy.int                                                               Number of atomic shells for which post-processing is possible.
+n_shells                 int                                                               Number of atomic shells for which post-processing is possible.
                                                                                            Note: this is `not` the number of correlated orbitals!
                                                                                            If there are two equivalent atoms in the unit cell, `n_shells` is 2.
 shells             list of dict {string:int}, dim n_shells x 4                             Atomic shell information.
@@ -46,17 +46,17 @@ shells             list of dict {string:int}, dim n_shells x 4
                                                                                            'l' is the angular quantum number, 'dim' is the dimension of the atomic shell.
                                                                                            e.g. for two equivalent atoms in the unit cell, `atom` runs from 0 to 1,
                                                                                            but `sort` can take only one value 0.
-n_corr_shells      numpy.int                                                               Number of correlated atomic shells.
+n_corr_shells            int                                                               Number of correlated atomic shells.
                                                                                            If there are two correlated equivalent atoms in the unit cell, `n_corr_shells` is 2.
-n_inequiv_shells   numpy.int                                                               Number of inequivalent atomic shells. Needs to be smaller than `n_corr_shells`.
+n_inequiv_shells         int                                                               Number of inequivalent atomic shells. Needs to be smaller than `n_corr_shells`.
                                                                                            The up / downfolding routines mediate between all correlated shells and the
                                                                                            actual inequivalent shells, by using the self-energy etc. for all equal shells
                                                                                            belonging to the same class of inequivalent shells. The mapping is performed with
                                                                                            information stored in `corr_to_inequiv` and `inequiv_to_corr`.
-corr_to_inequiv    list of numpy.int, dim `n_corr_shells`                                  mapping from correlated shells to inequivalent correlated shells.
+corr_to_inequiv    list of       int, dim `n_corr_shells`                                  mapping from correlated shells to inequivalent correlated shells.
                                                                                            A list of length `n_corr_shells` containing integers, where same numbers mark
                                                                                            equivalent sites.
-inequiv_to_corr    list of numpy.int, dim `n_inequiv_shells`                               A list of length `n_inequiv_shells` containing list indices as integers pointing
+inequiv_to_corr    list of       int, dim `n_inequiv_shells`                               A list of length `n_inequiv_shells` containing list indices as integers pointing
                                                                                            to the corresponding sites in `corr_to_inequiv`.
 corr_shells        list of dict {string:int}, dim n_corr_shells x 6                        Correlated orbital information.
                                                                                            For each correlated shell, have a dict with keys
@@ -64,28 +64,28 @@ corr_shells        list of dict {string:int}, dim n_corr_shells x 6
                                                                                            'atom' is the atom index, 'sort' defines the equivalency of the atoms,
                                                                                            'l' is the angular quantum number, 'dim' is the dimension of the atomic shell.
                                                                                            'SO' is one if spin-orbit is included, 0 otherwise, 'irep' is a dummy integer 0.
-use_rotations      numpy.int                                                               1 if local and global coordinate systems are used, 0 otherwise.
-rot_mat            list of numpy.array.complex,                                            Rotation matrices for correlated shells, if `use_rotations`.
+use_rotations            int                                                               1 if local and global coordinate systems are used, 0 otherwise.
+rot_mat            list of       array.complex,                                            Rotation matrices for correlated shells, if `use_rotations`.
                    dim n_corr_shells x [corr_shells['dim'],corr_shells['dim']]             These rotations are automatically applied for up / downfolding.
                                                                                            Set to the unity matrix if no rotations are used.
-rot_mat_time_inv   list of numpy.int, dim n_corr_shells                                    If `SP` is 1, 1 if the coordinate transformation contains inversion, 0 otherwise.
+rot_mat_time_inv   list of       int, dim n_corr_shells                                    If `SP` is 1, 1 if the coordinate transformation contains inversion, 0 otherwise.
                                                                                            If `use_rotations` or `SP` is 0, give a list of zeros.
-n_reps             numpy.int                                                               Number of irreducible representations of the correlated shell.
+n_reps                   int                                                               Number of irreducible representations of the correlated shell.
                                                                                            e.g. 2 if eg/t2g splitting is used.
-dim_reps           list of numpy.int, dim n_reps                                           Dimension of the representations.
+dim_reps           list of       int, dim n_reps                                           Dimension of the representations.
                                                                                            e.g. [2,3] for eg/t2g subsets.
-T                  list of numpy.array.complex,                                            Transformation matrix from the spherical harmonics to impurity problem basis
+T                  list of       array.complex,                                            Transformation matrix from the spherical harmonics to impurity problem basis
                    dim n_inequiv_corr_shell x                                              normally the real cubic harmonics).
                    [max(corr_shell['dim']),max(corr_shell['dim'])]                         This matrix can be used to calculate the 4-index U matrix, not automatically done.
-n_orbitals         numpy.array.int, dim [n_k,SP+1-SO]                                      Number of Bloch bands included in the projection window for each k-point.
+n_orbitals               array.int, dim [n_k,SP+1-SO]                                      Number of Bloch bands included in the projection window for each k-point.
                                                                                            If SP+1-SO=2, the number of included bands may depend on the spin projection up/down.
-proj_mat           numpy.array.complex,                                                    Projection matrices from Bloch bands to Wannier orbitals.
+proj_mat                 array.complex,                                                    Projection matrices from Bloch bands to Wannier orbitals.
                    dim [n_k,SP+1-SO,n_corr_shells,max(corr_shell['dim']),max(n_orbitals)]  For efficient storage reasons, all matrices must be of the same size
                                                                                            (given by last two indices).
                                                                                            For k-points with fewer bands, only the first entries are used, the rest are zero.
                                                                                            e.g. if number of Bloch bands ranges from 4-6, all matrices are of size 6.
-bz_weights         numpy.array.float, dim n_k                                              Weights of the k-points for the k summation. Soon be replaced by `kpt_weights`
-hopping            numpy.array.complex,                                                    Non-interacting Hamiltonian matrix for each k point.
+bz_weights               array.float, dim n_k                                              Weights of the k-points for the k summation. Soon be replaced by `kpt_weights`
+hopping                  array.complex,                                                    Non-interacting Hamiltonian matrix for each k point.
                    dim [n_k,SP+1-SO,max(n_orbitals),max(n_orbitals)]                       As for `proj_mat`, all matrices have to be of the same size.
 =================  ======================================================================  =====================================================================================
 
@@ -100,18 +100,18 @@ For the Vasp converter:
 =================  ======================================================================  =====================================================================================
 Name               Type                                                                    Meaning
 =================  ======================================================================  =====================================================================================
-kpt_basis          numpy.array.float, dim 3x3                                              Basis for the k-point mesh, reciprocal lattice vectors.
-kpts               numpy.array.float, dim n_k x 3                                          k-points given in reciprocal coordinates.
-kpt_weights        numpy.array.float, dim n_k                                              Weights of the k-points for the k summation.
+kpt_basis                array.float, dim 3x3                                              Basis for the k-point mesh, reciprocal lattice vectors.
+kpts                     array.float, dim n_k x 3                                          k-points given in reciprocal coordinates.
+kpt_weights              array.float, dim n_k                                              Weights of the k-points for the k summation.
 proj_or_hk         string                                                                  Switch determining whether the Vasp converter is running in projection mode `proj`, or
                                                                                            in Hamiltonian mode `hk`. In Hamiltonian mode, the hopping matrix is written in
                                                                                            orbital basis, whereas in projection mode hopping is written in band basis.
-proj_mat_csc       numpy.array.complex,                                                    Projection matrices from Bloch bands to Wannier orbitals for Hamiltonian based `hk`
+proj_mat_csc             array.complex,                                                    Projection matrices from Bloch bands to Wannier orbitals for Hamiltonian based `hk`
                    dim                                                                     approach. No site index is given, since hk is written in orbital basis. The last to
                    [n_k,SP+1-SO, n_corr_shells x max(corr_shell['dim']), max(n_orbitals)]  indices are a square matrix rotating from orbital to band space.
-dft_fermi_weights  numpy.array.float, dim n_k x 1 x max(n_orbitals)                        DFT fermi weights (occupations) of KS eigenstates for each k-point for calculation
+dft_fermi_weights        array.float, dim n_k x 1 x max(n_orbitals)                        DFT fermi weights (occupations) of KS eigenstates for each k-point for calculation
                    (stored in dft_misc_input)                                              of density matrix correction.
-band_window        list of numpy.array.int , dim(SP+1-SO)x n_k x 2                         Band windows as KS band indices in Vasp for each spin channel, and k-point. Needed for
+band_window        list of       array.int , dim(SP+1-SO)x n_k x 2                         Band windows as KS band indices in Vasp for each spin channel, and k-point. Needed for
                    (stored in dft_misc_input)                                              writing out the GAMMA file.
 =================  ======================================================================  =====================================================================================
 
diff --git a/python/triqs_dft_tools/block_structure.py b/python/triqs_dft_tools/block_structure.py
index 840779b2..b949d646 100644
--- a/python/triqs_dft_tools/block_structure.py
+++ b/python/triqs_dft_tools/block_structure.py
@@ -637,7 +637,7 @@ class BlockStructure(object):
 
         return self._create_gf_or_matrix(ish, gf_function, BlockGf, space, **kwargs)
 
-    def create_matrix(self, ish=0, space='solver', dtype=np.complex_):
+    def create_matrix(self, ish=0, space='solver', dtype=complex):
         """ Create a zero matrix having the correct structure.
 
         For ``space='solver'``, the structure is according to
diff --git a/python/triqs_dft_tools/converters/hk.py b/python/triqs_dft_tools/converters/hk.py
index 0984a7b4..16704506 100644
--- a/python/triqs_dft_tools/converters/hk.py
+++ b/python/triqs_dft_tools/converters/hk.py
@@ -131,7 +131,7 @@ class HkConverter(ConverterTools):
 
             use_rotations = 0
             rot_mat = [numpy.identity(
-                corr_shells[icrsh]['dim'], numpy.complex_) for icrsh in range(n_corr_shells)]
+                corr_shells[icrsh]['dim'], complex) for icrsh in range(n_corr_shells)]
             rot_mat_time_inv = [0 for i in range(n_corr_shells)]
 
             # Representative representations are read from file
@@ -149,7 +149,7 @@ class HkConverter(ConverterTools):
                 # Wien2k)
                 ll = 2 * corr_shells[inequiv_to_corr[ish]]['l'] + 1
                 lmax = ll * (corr_shells[inequiv_to_corr[ish]]['SO'] + 1)
-                T.append(numpy.zeros([lmax, lmax], numpy.complex_))
+                T.append(numpy.zeros([lmax, lmax], complex))
 
                 T[ish] = numpy.array([[0.0, 0.0, 1.0, 0.0, 0.0],
                                       [1.0 / sqrt(2.0), 0.0, 0.0,
@@ -167,11 +167,11 @@ class HkConverter(ConverterTools):
             # define the number of n_orbitals for all k points: it is the
             # number of total bands and independent of k!
             n_orbitals = numpy.ones(
-                [n_k, n_spin_blocs], numpy.int) * sum([sh['dim'] for sh in shells])
+                [n_k, n_spin_blocs], int) * sum([sh['dim'] for sh in shells])
 
             # Initialise the projectors:
             proj_mat = numpy.zeros([n_k, n_spin_blocs, n_corr_shells, max(
-                [crsh['dim'] for crsh in corr_shells]), numpy.max(n_orbitals)], numpy.complex_)
+                [crsh['dim'] for crsh in corr_shells]), numpy.max(n_orbitals)], complex)
 
             # Read the projectors from the file:
             for ik in range(n_k):
@@ -193,9 +193,9 @@ class HkConverter(ConverterTools):
 
             # now define the arrays for weights and hopping ...
             # w(k_index),  default normalisation
-            bz_weights = numpy.ones([n_k], numpy.float_) / float(n_k)
+            bz_weights = numpy.ones([n_k], float) / float(n_k)
             hopping = numpy.zeros([n_k, n_spin_blocs, numpy.max(
-                n_orbitals), numpy.max(n_orbitals)], numpy.complex_)
+                n_orbitals), numpy.max(n_orbitals)], complex)
 
             if (weights_in_file):
                 # weights in the file
diff --git a/python/triqs_dft_tools/converters/plovasp/elstruct.py b/python/triqs_dft_tools/converters/plovasp/elstruct.py
index 6418b5ea..728faad0 100644
--- a/python/triqs_dft_tools/converters/plovasp/elstruct.py
+++ b/python/triqs_dft_tools/converters/plovasp/elstruct.py
@@ -152,10 +152,10 @@ class ElectronicStructure:
 # Spin factor
         sp_fac = 2.0 if ns == 1 and not self.nc_flag else 1.0
 
-        den_mat = np.zeros((ns, nproj, nproj), dtype=np.float64)
-        overlap = np.zeros((ns, nproj, nproj), dtype=np.float64)
-#        ov_min = np.ones((ns, nproj, nproj), dtype=np.float64) * 100.0
-#        ov_max = np.zeros((ns, nproj, nproj), dtype=np.float64)
+        den_mat = np.zeros((ns, nproj, nproj), dtype=float)
+        overlap = np.zeros((ns, nproj, nproj), dtype=float)
+#        ov_min = np.ones((ns, nproj, nproj), dtype=float) * 100.0
+#        ov_max = np.zeros((ns, nproj, nproj), dtype=float)
         for ispin in range(ns):
             for ik in range(nk):
                 kweight = self.kmesh['kweights'][ik]
diff --git a/python/triqs_dft_tools/converters/plovasp/inpconf.py b/python/triqs_dft_tools/converters/plovasp/inpconf.py
index 0c2b3a4e..1be71172 100644
--- a/python/triqs_dft_tools/converters/plovasp/inpconf.py
+++ b/python/triqs_dft_tools/converters/plovasp/inpconf.py
@@ -268,7 +268,7 @@ class ConfigParameters:
             err_mess = "Complex matrix must contain 2*M values:\n%s"%(par_str)
             assert 2 * (nm // 2) == nm, err_mess
 
-            tmp = np.array(rows, dtype=np.complex128)
+            tmp = np.array(rows, dtype=complex)
             mat = tmp[:, 0::2] + 1.0j * tmp[:, 1::2]
 
         return mat
diff --git a/python/triqs_dft_tools/converters/plovasp/proj_group.py b/python/triqs_dft_tools/converters/plovasp/proj_group.py
index 0237ae6e..28628a2f 100644
--- a/python/triqs_dft_tools/converters/plovasp/proj_group.py
+++ b/python/triqs_dft_tools/converters/plovasp/proj_group.py
@@ -68,7 +68,7 @@ class ProjectorGroup:
 # Determine the minimum and maximum band numbers
         if 'bands' in gr_pars:
             nk, nband, ns_band = eigvals.shape
-            ib_win = np.zeros((nk, ns_band, 2), dtype=np.int32)
+            ib_win = np.zeros((nk, ns_band, 2), dtype=int)
             ib_win[:,:,0] = gr_pars['bands'][0]-1
             ib_win[:,:,1] = gr_pars['bands'][1]-1
             ib_min = gr_pars['bands'][0] - 1
@@ -152,7 +152,7 @@ class ProjectorGroup:
         block_maps, ndim = self.get_block_matrix_map()
 
         _, ns, nk, _, _ = self.shells[0].proj_win.shape
-        p_mat = np.zeros((ndim, self.nb_max), dtype=np.complex128)
+        p_mat = np.zeros((ndim, self.nb_max), dtype=complex)
 # Note that 'ns' and 'nk' are the same for all shells
         for isp in range(ns):
             for ik in range(nk):
@@ -201,7 +201,7 @@ class ProjectorGroup:
 
         _, ns, nk, _, _ = self.shells[0].proj_win.shape
 
-        self.hk = np.zeros((ns,nk,ndim,ndim), dtype=np.complex128)
+        self.hk = np.zeros((ns,nk,ndim,ndim), dtype=complex)
 # Note that 'ns' and 'nk' are the same for all shells
         for isp in range(ns):
             for ik in range(nk):
@@ -209,7 +209,7 @@ class ProjectorGroup:
                 bmax = self.ib_win[ik, isp, 1]+1
 
                 nb = bmax - bmin
-                p_mat = np.zeros((ndim, nb), dtype=np.complex128)
+                p_mat = np.zeros((ndim, nb), dtype=complex)
                 #print(bmin,bmax,nb)
 # Combine all projectors of the group to one block projector
                 for bl_map in block_maps:
@@ -251,8 +251,8 @@ class ProjectorGroup:
 
         block_maps, ndim = self.get_block_matrix_map()
         _, ns, nk, _, _ = self.shells[0].proj_win.shape
-        p_mat = np.zeros((ndim, self.nb_max), dtype=np.complex128)
-        p_full = np.zeros((1,ns,nk,self.nb_max, self.nb_max), dtype=np.complex128)
+        p_mat = np.zeros((ndim, self.nb_max), dtype=complex)
+        p_full = np.zeros((1,ns,nk,self.nb_max, self.nb_max), dtype=complex)
 
 # Note that 'ns' and 'nk' are the same for all shells
 
@@ -452,7 +452,7 @@ class ProjectorGroup:
             raise Exception("Energy window does not overlap with the band structure")
 
         nk, nband, ns_band = eigvals.shape
-        ib_win = np.zeros((nk, ns_band, 2), dtype=np.int32)
+        ib_win = np.zeros((nk, ns_band, 2), dtype=int)
 
         ib_min = 10000000
         ib_max = 0
diff --git a/python/triqs_dft_tools/converters/plovasp/proj_shell.py b/python/triqs_dft_tools/converters/plovasp/proj_shell.py
index 628357b0..483f16f5 100644
--- a/python/triqs_dft_tools/converters/plovasp/proj_shell.py
+++ b/python/triqs_dft_tools/converters/plovasp/proj_shell.py
@@ -155,7 +155,7 @@ class ProjectorShell:
             assert nr%ns_dim == 0, "Number of rows in TRANSFILE is not compatible with the spin dimension"
             ndim = nr // ns_dim
 
-            self.tmatrices = np.zeros((nion, nr, nm * ns_dim), dtype=np.complex128)
+            self.tmatrices = np.zeros((nion, nr, nm * ns_dim), dtype=complex)
 
             if is_complex:
                 raw_matrices = raw_matrices[:, ::2] + raw_matrices[:, 1::2] * 1j
@@ -187,7 +187,7 @@ class ProjectorShell:
 
             ndim = nrow
 
-            self.tmatrices = np.zeros((nion, nrow, nm), dtype=np.complex128)
+            self.tmatrices = np.zeros((nion, nrow, nm), dtype=complex)
             for io in range(nion):
                 self.tmatrices[io, :, :] = raw_matrix
 
@@ -200,9 +200,9 @@ class ProjectorShell:
         ndim = nm * ns_dim
 
 # We still need the matrices for the output
-        self.tmatrices = np.zeros((nion, ndim, ndim), dtype=np.complex128)
+        self.tmatrices = np.zeros((nion, ndim, ndim), dtype=complex)
         for io in range(nion):
-            self.tmatrices[io, :, :] = np.identity(ndim, dtype=np.complex128)
+            self.tmatrices[io, :, :] = np.identity(ndim, dtype=complex)
 
         return ndim
 
@@ -230,11 +230,11 @@ class ProjectorShell:
 # TODO: implement a non-collinear case
 #       for a non-collinear case 'ndim' is 'ns * nm'
             ndim = self.tmatrices.shape[1]
-            self.proj_arr = np.zeros((nion, ns, nk, ndim, nb), dtype=np.complex128)
+            self.proj_arr = np.zeros((nion, ns, nk, ndim, nb), dtype=complex)
             for ik in range(nk):
                 kp = kmesh['kpoints'][ik]
                 for io, ion in enumerate(self.ion_list):
-                    proj_k = np.zeros((ns, nlm, nb), dtype=np.complex128)
+                    proj_k = np.zeros((ns, nlm, nb), dtype=complex)
                     qcoord = structure['qcoords'][ion]
 #                    kphase = np.exp(-2.0j * np.pi * np.dot(kp, qcoord))
 #                    kphase = 1.0
@@ -249,7 +249,7 @@ class ProjectorShell:
 
         else:
 # No transformation: just copy the projectors as they are
-            self.proj_arr = np.zeros((nion, ns, nk, nlm, nb), dtype=np.complex128)
+            self.proj_arr = np.zeros((nion, ns, nk, nlm, nb), dtype=complex)
             for io, ion in enumerate(self.ion_list):
                 qcoord = structure['qcoords'][ion]
                 for m in range(nlm):
@@ -282,7 +282,7 @@ class ProjectorShell:
 # Set the dimensions of the array
         nion, ns, nk, nlm, nbtot = self.proj_arr.shape
 # !!! Note that the order of the two last indices is different !!!
-        self.proj_win = np.zeros((nion, ns, nk, nlm, nb_max), dtype=np.complex128)
+        self.proj_win = np.zeros((nion, ns, nk, nlm, nb_max), dtype=complex)
 
 # Select projectors for a given energy window
         ns_band = self.ib_win.shape[1]
@@ -310,14 +310,14 @@ class ProjectorShell:
         assert spin_diag, "spin_diag = False is not implemented"
 
         if site_diag:
-            occ_mats = np.zeros((ns, nion, nlm, nlm), dtype=np.float64)
-            overlaps = np.zeros((ns, nion, nlm, nlm), dtype=np.float64)
+            occ_mats = np.zeros((ns, nion, nlm, nlm), dtype=float)
+            overlaps = np.zeros((ns, nion, nlm, nlm), dtype=float)
         else:
             ndim = nion * nlm
-            occ_mats = np.zeros((ns, 1, ndim, ndim), dtype=np.float64)
-            overlaps = np.zeros((ns, 1, ndim, ndim), dtype=np.float64)
+            occ_mats = np.zeros((ns, 1, ndim, ndim), dtype=float)
+            overlaps = np.zeros((ns, 1, ndim, ndim), dtype=float)
 
-#        self.proj_win = np.zeros((nion, ns, nk, nlm, nb_max), dtype=np.complex128)
+#        self.proj_win = np.zeros((nion, ns, nk, nlm, nb_max), dtype=complex)
         kweights = el_struct.kmesh['kweights']
         occnums = el_struct.ferw
         ib1 = self.ib_min
@@ -332,7 +332,7 @@ class ProjectorShell:
                         overlaps[isp, io, :, :] += np.dot(proj_k,
                                                      proj_k.conj().T).real * weight
         else:
-            proj_k = np.zeros((ndim, nbtot), dtype=np.complex128)
+            proj_k = np.zeros((ndim, nbtot), dtype=complex)
             for isp in range(ns):
                 for ik, weight, occ in zip(it.count(), kweights, occnums[isp, :, :]):
                     for io in range(nion):
@@ -363,9 +363,9 @@ class ProjectorShell:
         assert site_diag, "site_diag = False is not implemented"
         assert spin_diag, "spin_diag = False is not implemented"
 
-        loc_ham = np.zeros((ns, nion, nlm, nlm), dtype=np.complex128)
+        loc_ham = np.zeros((ns, nion, nlm, nlm), dtype=complex)
 
-#        self.proj_win = np.zeros((nion, ns, nk, nlm, nb_max), dtype=np.complex128)
+#        self.proj_win = np.zeros((nion, ns, nk, nlm, nb_max), dtype=complex)
         kweights = el_struct.kmesh['kweights']
         occnums = el_struct.ferw
         ib1 = self.ib_min
@@ -403,7 +403,7 @@ class ProjectorShell:
 
         ne = len(emesh)
         dos = np.zeros((ne, ns, nion, nlm))
-        w_k = np.zeros((nk, nb_max, ns, nion, nlm), dtype=np.complex128)
+        w_k = np.zeros((nk, nb_max, ns, nion, nlm), dtype=complex)
         for isp in range(ns):
             for ik in range(nk):
                 is_b = min(isp, ns_band)
diff --git a/python/triqs_dft_tools/converters/plovasp/vaspio.py b/python/triqs_dft_tools/converters/plovasp/vaspio.py
index 322bf47d..0b0876b7 100644
--- a/python/triqs_dft_tools/converters/plovasp/vaspio.py
+++ b/python/triqs_dft_tools/converters/plovasp/vaspio.py
@@ -163,7 +163,7 @@ class Plocar:
             line = f.readline()
             nproj, nspin, nk, nband = list(map(int, line.split()))
 
-        plo = np.zeros((nproj, nspin, nk, nband), dtype=np.complex128)
+        plo = np.zeros((nproj, nspin, nk, nband), dtype=complex)
         proj_params = [{} for i in range(nproj)]
 
         iproj_site = 0
@@ -251,7 +251,7 @@ class Plocar:
             except:
                 print("!!! WARNING !!!: Error reading E-Fermi from LOCPROJ, trying DOSCAR")
 
-            plo = np.zeros((nproj, self.nspin, nk, self.nband), dtype=np.complex128)
+            plo = np.zeros((nproj, self.nspin, nk, self.nband), dtype=complex)
             proj_params = [{} for i in range(nproj)]
 
             iproj_site = 0
@@ -685,7 +685,7 @@ def read_symmcar(vasp_dir, symm_filename='SYMMCAR'):
     print("   {0:>26} {1:d}".format("L_max:", lmax))
 
     rot_mats = np.zeros((nrot, lmax+1, mmax, mmax))
-    rot_map = np.zeros((nrot, ntrans, nion), dtype=np.int32)
+    rot_map = np.zeros((nrot, ntrans, nion), dtype=int)
 
     for irot in range(nrot):
 #   Empty line
diff --git a/python/triqs_dft_tools/converters/vasp.py b/python/triqs_dft_tools/converters/vasp.py
index 6aaf521e..27347def 100644
--- a/python/triqs_dft_tools/converters/vasp.py
+++ b/python/triqs_dft_tools/converters/vasp.py
@@ -256,7 +256,7 @@ class VaspConverter(ConverterTools):
 
 # NB!: these rotation matrices are specific to Wien2K! Set to identity in VASP
             use_rotations = 1
-            rot_mat = [numpy.identity(corr_shells[icrsh]['dim'],numpy.complex_) for icrsh in range(n_corr_shells)]
+            rot_mat = [numpy.identity(corr_shells[icrsh]['dim'],complex) for icrsh in range(n_corr_shells)]
             rot_mat_time_inv = [0 for i in range(n_corr_shells)]
 
 # TODO: implement transformation matrices
@@ -273,16 +273,16 @@ class VaspConverter(ConverterTools):
                 ll = 2 * corr_shells[inequiv_to_corr[ish]]['l']+1
                 lmax = ll * (corr_shells[inequiv_to_corr[ish]]['SO'] + 1)
 # TODO: at the moment put T-matrices to identities
-                T.append(numpy.identity(lmax, numpy.complex_))
+                T.append(numpy.identity(lmax, complex))
 
 #            if nc_flag:
 ## TODO: implement the noncollinear part
 #                raise NotImplementedError("Noncollinear calculations are not implemented")
 #            else:
-            hopping = numpy.zeros([n_k, n_spin_blocs, nb_max, nb_max], numpy.complex_)
-            f_weights = numpy.zeros([n_k, n_spin_blocs, nb_max], numpy.complex_)
+            hopping = numpy.zeros([n_k, n_spin_blocs, nb_max, nb_max], complex)
+            f_weights = numpy.zeros([n_k, n_spin_blocs, nb_max], complex)
             band_window = [numpy.zeros((n_k, 2), dtype=int) for isp in range(n_spin_blocs)]
-            n_orbitals = numpy.zeros([n_k, n_spin_blocs], numpy.int)
+            n_orbitals = numpy.zeros([n_k, n_spin_blocs], int)
 
 
             for isp in range(n_spin_blocs):
@@ -296,7 +296,7 @@ class VaspConverter(ConverterTools):
                         f_weights[ik, isp, ib] = next(rf)
 
             if self.proj_or_hk == 'hk':
-                hopping = numpy.zeros([n_k, n_spin_blocs, n_orbs, n_orbs], numpy.complex_)
+                hopping = numpy.zeros([n_k, n_spin_blocs, n_orbs, n_orbs], complex)
                 # skip header lines
                 hk_file = self.basename + '.hk%i'%(ig + 1)
                 f_hk = open(hk_file, 'rt')
@@ -321,7 +321,7 @@ class VaspConverter(ConverterTools):
 # Projectors
 #            print n_orbitals
 #            print [crsh['dim'] for crsh in corr_shells]
-            proj_mat_csc = numpy.zeros([n_k, n_spin_blocs, sum([sh['dim'] for sh in shells]), numpy.max(n_orbitals)], numpy.complex_)
+            proj_mat_csc = numpy.zeros([n_k, n_spin_blocs, sum([sh['dim'] for sh in shells]), numpy.max(n_orbitals)], complex)
 
 # TODO: implement reading from more than one projector group
 # In 'dmftproj' each ion represents a separate correlated shell.
@@ -348,7 +348,7 @@ class VaspConverter(ConverterTools):
                                     proj_mat_csc[ik, isp, ilm, ib] = complex(pr, pi)
 
 # now save only projectors with flag 'corr' to proj_mat
-            proj_mat = numpy.zeros([n_k, n_spin_blocs, n_corr_shells, max([crsh['dim'] for crsh in corr_shells]), numpy.max(n_orbitals)], numpy.complex_)
+            proj_mat = numpy.zeros([n_k, n_spin_blocs, n_corr_shells, max([crsh['dim'] for crsh in corr_shells]), numpy.max(n_orbitals)], complex)
             if self.proj_or_hk == 'proj':
                 for ish, sh in enumerate(p_shells):
                     if sh['corr']:
diff --git a/python/triqs_dft_tools/converters/wannier90.py b/python/triqs_dft_tools/converters/wannier90.py
index e0d297ee..535dd769 100644
--- a/python/triqs_dft_tools/converters/wannier90.py
+++ b/python/triqs_dft_tools/converters/wannier90.py
@@ -200,7 +200,7 @@ class Wannier90Converter(ConverterTools):
         for ish in range(n_inequiv_shells):
             ll = 2 * corr_shells[inequiv_to_corr[ish]]['l'] + 1
             lmax = ll * (corr_shells[inequiv_to_corr[ish]]['SO'] + 1)
-            T.append(numpy.zeros([lmax, lmax], numpy.complex_))
+            T.append(numpy.zeros([lmax, lmax], complex))
 
         spin_w90name = ['_up', '_down']
         hamr_full = []
@@ -267,7 +267,7 @@ class Wannier90Converter(ConverterTools):
                 # we assume spin up and spin down always have same total number
                 # of WFs
                 n_orbitals = numpy.ones(
-                    [self.n_k, n_spin], numpy.int) * self.nwfs
+                    [self.n_k, n_spin], int) * self.nwfs
 
             else:
                 # consistency check between the _up and _down file contents
@@ -322,7 +322,7 @@ class Wannier90Converter(ConverterTools):
 
         # Third, compute the hoppings in reciprocal space
         hopping = numpy.zeros([self.n_k, n_spin, numpy.max(
-            n_orbitals), numpy.max(n_orbitals)], numpy.complex_)
+            n_orbitals), numpy.max(n_orbitals)], complex)
         for isp in range(n_spin):
             # make Fourier transform H(R) -> H(k) : it can be done one spin at
             # a time
@@ -336,7 +336,7 @@ class Wannier90Converter(ConverterTools):
         # Then, initialise the projectors
         k_dep_projection = 0   # we always have the same number of WFs at each k-point
         proj_mat = numpy.zeros([self.n_k, n_spin, n_corr_shells, max(
-            [crsh['dim'] for crsh in corr_shells]), numpy.max(n_orbitals)], numpy.complex_)
+            [crsh['dim'] for crsh in corr_shells]), numpy.max(n_orbitals)], complex)
         iorb = 0
         # Projectors simply consist in identity matrix blocks selecting those MLWFs that
         # correspond to the specific correlated shell indexed by icrsh.
@@ -346,7 +346,7 @@ class Wannier90Converter(ConverterTools):
         for icrsh in range(n_corr_shells):
             norb = corr_shells[icrsh]['dim']
             proj_mat[:, :, icrsh, 0:norb, iorb:iorb +
-                     norb] = numpy.identity(norb, numpy.complex_)
+                     norb] = numpy.identity(norb, complex)
             iorb += norb
 
         # Finally, save all required data into the HDF archive:
@@ -410,7 +410,7 @@ class Wannier90Converter(ConverterTools):
         # Hamiltonian
         rvec_idx = numpy.zeros((nrpt, 3), dtype=int)
         rvec_deg = numpy.zeros(nrpt, dtype=int)
-        h_of_r = [numpy.zeros((num_wf, num_wf), dtype=numpy.complex_)
+        h_of_r = [numpy.zeros((num_wf, num_wf), dtype=complex)
                   for n in range(nrpt)]
 
         # variable currpos points to the current line in the file
@@ -607,7 +607,7 @@ class Wannier90Converter(ConverterTools):
         """
 
         twopi = 2 * numpy.pi
-        h_of_k = [numpy.zeros((norb, norb), dtype=numpy.complex_)
+        h_of_k = [numpy.zeros((norb, norb), dtype=complex)
                   for ik in range(self.n_k)]
         ridx = numpy.array(list(range(self.nrpt)))
         for ik, ir in product(list(range(self.n_k)), ridx):
diff --git a/python/triqs_dft_tools/converters/wien2k.py b/python/triqs_dft_tools/converters/wien2k.py
index b45ffa72..6ecfaa55 100644
--- a/python/triqs_dft_tools/converters/wien2k.py
+++ b/python/triqs_dft_tools/converters/wien2k.py
@@ -152,7 +152,7 @@ class Wien2kConverter(ConverterTools):
 
             use_rotations = 1
             rot_mat = [numpy.identity(
-                corr_shells[icrsh]['dim'], numpy.complex_) for icrsh in range(n_corr_shells)]
+                corr_shells[icrsh]['dim'], complex) for icrsh in range(n_corr_shells)]
 
             # read the matrices
             rot_mat_time_inv = [0 for i in range(n_corr_shells)]
@@ -183,7 +183,7 @@ class Wien2kConverter(ConverterTools):
                 # is of dimension 2l+1 without SO, and 2*(2l+1) with SO!
                 ll = 2 * corr_shells[inequiv_to_corr[ish]]['l'] + 1
                 lmax = ll * (corr_shells[inequiv_to_corr[ish]]['SO'] + 1)
-                T.append(numpy.zeros([lmax, lmax], numpy.complex_))
+                T.append(numpy.zeros([lmax, lmax], complex))
 
                 # now read it from file:
                 for i in range(lmax):
@@ -197,14 +197,14 @@ class Wien2kConverter(ConverterTools):
             n_spin_blocs = SP + 1 - SO
 
             # read the list of n_orbitals for all k points
-            n_orbitals = numpy.zeros([n_k, n_spin_blocs], numpy.int)
+            n_orbitals = numpy.zeros([n_k, n_spin_blocs], int)
             for isp in range(n_spin_blocs):
                 for ik in range(n_k):
                     n_orbitals[ik, isp] = int(next(R))
 
             # Initialise the projectors:
             proj_mat = numpy.zeros([n_k, n_spin_blocs, n_corr_shells, max(
-                [crsh['dim'] for crsh in corr_shells]), numpy.max(n_orbitals)], numpy.complex_)
+                [crsh['dim'] for crsh in corr_shells]), numpy.max(n_orbitals)], complex)
 
             # Read the projectors from the file:
             for ik in range(n_k):
@@ -224,9 +224,9 @@ class Wien2kConverter(ConverterTools):
 
             # now define the arrays for weights and hopping ...
             # w(k_index),  default normalisation
-            bz_weights = numpy.ones([n_k], numpy.float_) / float(n_k)
+            bz_weights = numpy.ones([n_k], float) / float(n_k)
             hopping = numpy.zeros([n_k, n_spin_blocs, numpy.max(
-                n_orbitals), numpy.max(n_orbitals)], numpy.complex_)
+                n_orbitals), numpy.max(n_orbitals)], complex)
 
             # weights in the file
             for ik in range(n_k):
@@ -301,7 +301,7 @@ class Wien2kConverter(ConverterTools):
 
         mpi.report("Reading input from %s..." % self.parproj_file)
 
-        dens_mat_below = [[numpy.zeros([self.shells[ish]['dim'], self.shells[ish]['dim']], numpy.complex_) for ish in range(self.n_shells)]
+        dens_mat_below = [[numpy.zeros([self.shells[ish]['dim'], self.shells[ish]['dim']], complex) for ish in range(self.n_shells)]
                           for isp in range(self.n_spin_blocs)]
 
         R = ConverterTools.read_fortran_file(
@@ -312,10 +312,10 @@ class Wien2kConverter(ConverterTools):
 
         # Initialise P, here a double list of matrices:
         proj_mat_all = numpy.zeros([self.n_k, self.n_spin_blocs, self.n_shells, max(
-            n_parproj), max([sh['dim'] for sh in self.shells]), numpy.max(self.n_orbitals)], numpy.complex_)
+            n_parproj), max([sh['dim'] for sh in self.shells]), numpy.max(self.n_orbitals)], complex)
 
         rot_mat_all = [numpy.identity(
-            self.shells[ish]['dim'], numpy.complex_) for ish in range(self.n_shells)]
+            self.shells[ish]['dim'], complex) for ish in range(self.n_shells)]
         rot_mat_all_time_inv = [0 for i in range(self.n_shells)]
 
         for ish in range(self.n_shells):
@@ -406,14 +406,14 @@ class Wien2kConverter(ConverterTools):
             n_k = int(next(R))
 
             # read the list of n_orbitals for all k points
-            n_orbitals = numpy.zeros([n_k, self.n_spin_blocs], numpy.int)
+            n_orbitals = numpy.zeros([n_k, self.n_spin_blocs], int)
             for isp in range(self.n_spin_blocs):
                 for ik in range(n_k):
                     n_orbitals[ik, isp] = int(next(R))
 
             # Initialise the projectors:
             proj_mat = numpy.zeros([n_k, self.n_spin_blocs, self.n_corr_shells, max(
-                [crsh['dim'] for crsh in self.corr_shells]), numpy.max(n_orbitals)], numpy.complex_)
+                [crsh['dim'] for crsh in self.corr_shells]), numpy.max(n_orbitals)], complex)
 
             # Read the projectors from the file:
             for ik in range(n_k):
@@ -432,7 +432,7 @@ class Wien2kConverter(ConverterTools):
                                 proj_mat[ik, isp, icrsh, i, j] += 1j * next(R)
 
             hopping = numpy.zeros([n_k, self.n_spin_blocs, numpy.max(
-                n_orbitals), numpy.max(n_orbitals)], numpy.complex_)
+                n_orbitals), numpy.max(n_orbitals)], complex)
 
             # Grab the H
             # we use now the convention of a DIAGONAL Hamiltonian!!!!
@@ -448,7 +448,7 @@ class Wien2kConverter(ConverterTools):
 
             # Initialise P, here a double list of matrices:
             proj_mat_all = numpy.zeros([n_k, self.n_spin_blocs, self.n_shells, max(n_parproj), max(
-                [sh['dim'] for sh in self.shells]), numpy.max(n_orbitals)], numpy.complex_)
+                [sh['dim'] for sh in self.shells]), numpy.max(n_orbitals)], complex)
 
             for ish in range(self.n_shells):
                 for ik in range(n_k):
@@ -751,7 +751,7 @@ class Wien2kConverter(ConverterTools):
             for i_symm in range(n_symm):
 
                 mat.append([numpy.zeros([orbits[orb]['dim'], orbits[orb][
-                           'dim']], numpy.complex_) for orb in range(n_orbits)])
+                           'dim']], complex) for orb in range(n_orbits)])
                 for orb in range(n_orbits):
                     for i in range(orbits[orb]['dim']):
                         for j in range(orbits[orb]['dim']):
@@ -762,7 +762,7 @@ class Wien2kConverter(ConverterTools):
                             mat[i_symm][orb][i, j] += 1j * \
                                 next(R)      # imaginary part
 
-            mat_tinv = [numpy.identity(orbits[orb]['dim'], numpy.complex_)
+            mat_tinv = [numpy.identity(orbits[orb]['dim'], complex)
                         for orb in range(n_orbits)]
 
             if ((SO == 0) and (SP == 0)):
diff --git a/python/triqs_dft_tools/sumk_dft.py b/python/triqs_dft_tools/sumk_dft.py
index 28ff9ac8..427ddfb9 100644
--- a/python/triqs_dft_tools/sumk_dft.py
+++ b/python/triqs_dft_tools/sumk_dft.py
@@ -576,7 +576,7 @@ class SumkDFT(object):
             G_latt << Omega + 1j * broadening
 
         idmat = [numpy.identity(
-            self.n_orbitals[ik, ntoi[sp]], numpy.complex_) for sp in spn]
+            self.n_orbitals[ik, ntoi[sp]], complex) for sp in spn]
         M = copy.deepcopy(idmat)
 
         for ibl in range(self.n_spin_blocks[self.SO]):
@@ -928,7 +928,7 @@ class SumkDFT(object):
             for block, inner in self.gf_struct_solver[ish].items():
                 # get dm for the blocks:
                 dm[block] = numpy.zeros(
-                    [len(inner), len(inner)], numpy.complex_)
+                    [len(inner), len(inner)], complex)
                 for ind1 in inner:
                     for ind2 in inner:
                         block_sumk, ind1_sumk = self.solver_to_sumk[
@@ -1275,7 +1275,7 @@ class SumkDFT(object):
                         """
                         def chi2(y):
                             # reinterpret y as complex number
-                            y = y.view(numpy.complex_)
+                            y = y.view(complex)
                             ret = 0.0
                             for a in range(Z.shape[0]):
                                 for b in range(Z.shape[1]):
@@ -1292,10 +1292,10 @@ class SumkDFT(object):
                         if res.fun > threshold: continue
 
                         # reinterpret the solution as a complex number
-                        y = res.x.view(numpy.complex_)
+                        y = res.x.view(complex)
 
                         # reconstruct the T matrix
-                        T = numpy.zeros(N.shape[:-1], dtype=numpy.complex_)
+                        T = numpy.zeros(N.shape[:-1], dtype=complex)
                         for i in range(len(y)):
                             T += N[:, :, i] * y[i]
 
@@ -1465,7 +1465,7 @@ class SumkDFT(object):
         for icrsh in range(self.n_corr_shells):
             for sp in self.spin_block_names[self.corr_shells[icrsh]['SO']]:
                 dens_mat[icrsh][sp] = numpy.zeros(
-                    [self.corr_shells[icrsh]['dim'], self.corr_shells[icrsh]['dim']], numpy.complex_)
+                    [self.corr_shells[icrsh]['dim'], self.corr_shells[icrsh]['dim']], complex)
 
         ikarray = numpy.array(list(range(self.n_k)))
         for ik in mpi.slice_array(ikarray):
@@ -1483,7 +1483,7 @@ class SumkDFT(object):
                 ntoi = self.spin_names_to_ind[self.SO]
                 spn = self.spin_block_names[self.SO]
                 dims = {sp:self.n_orbitals[ik, ntoi[sp]] for sp in spn}
-                MMat = [numpy.zeros([dims[sp], dims[sp]], numpy.complex_) for sp in spn]
+                MMat = [numpy.zeros([dims[sp], dims[sp]], complex) for sp in spn]
 
                 for isp, sp in enumerate(spn):
                     ind = ntoi[sp]
@@ -1564,7 +1564,7 @@ class SumkDFT(object):
         for ish in range(self.n_inequiv_shells):
             for sp in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]]['SO']]:
                 eff_atlevels[ish][sp] = numpy.identity(
-                    self.corr_shells[self.inequiv_to_corr[ish]]['dim'], numpy.complex_)
+                    self.corr_shells[self.inequiv_to_corr[ish]]['dim'], complex)
                 eff_atlevels[ish][sp] *= -self.chemical_potential
                 eff_atlevels[ish][
                     sp] -= self.dc_imp[self.inequiv_to_corr[ish]][sp]
@@ -1578,13 +1578,13 @@ class SumkDFT(object):
                 dim = self.corr_shells[icrsh]['dim']
                 for sp in self.spin_block_names[self.corr_shells[icrsh]['SO']]:
                     self.Hsumk[icrsh][sp] = numpy.zeros(
-                        [dim, dim], numpy.complex_)
+                        [dim, dim], complex)
                 for isp, sp in enumerate(self.spin_block_names[self.corr_shells[icrsh]['SO']]):
                     ind = self.spin_names_to_ind[
                         self.corr_shells[icrsh]['SO']][sp]
                     for ik in range(self.n_k):
                         n_orb = self.n_orbitals[ik, ind]
-                        MMat = numpy.identity(n_orb, numpy.complex_)
+                        MMat = numpy.identity(n_orb, complex)
                         MMat = self.hopping[
                             ik, ind, 0:n_orb, 0:n_orb] - (1 - 2 * isp) * self.h_field * MMat
                         projmat = self.proj_mat[ik, ind, icrsh, 0:dim, 0:n_orb]
@@ -1626,7 +1626,7 @@ class SumkDFT(object):
             dim = self.corr_shells[icrsh]['dim']
             spn = self.spin_block_names[self.corr_shells[icrsh]['SO']]
             for sp in spn:
-                self.dc_imp[icrsh][sp] = numpy.zeros([dim, dim], numpy.float_)
+                self.dc_imp[icrsh][sp] = numpy.zeros([dim, dim], float)
         self.dc_energ = [0.0 for icrsh in range(self.n_corr_shells)]
 
     def set_dc(self, dc_imp, dc_energ):
@@ -1704,7 +1704,7 @@ class SumkDFT(object):
                 Ncr[bl] += dens_mat[block].real.trace()
             Ncrtot = sum(Ncr.values())
             for sp in spn:
-                self.dc_imp[icrsh][sp] = numpy.identity(dim, numpy.float_)
+                self.dc_imp[icrsh][sp] = numpy.identity(dim, float)
                 if self.SP == 0:  # average the densities if there is no SP:
                     Ncr[sp] = Ncrtot / len(spn)
                 # correction for SO: we have only one block in this case, but
@@ -2025,14 +2025,14 @@ class SumkDFT(object):
 # Convert Fermi weights to a density matrix
             dens_mat_dft = {}
             for sp in spn:
-                dens_mat_dft[sp] = [fermi_weights[ik, ntoi[sp], :].astype(numpy.complex_) for ik in range(self.n_k)]
+                dens_mat_dft[sp] = [fermi_weights[ik, ntoi[sp], :].astype(complex) for ik in range(self.n_k)]
 
 
         # Set up deltaN:
         deltaN = {}
         for sp in spn:
             deltaN[sp] = [numpy.zeros([self.n_orbitals[ik, ntoi[sp]], self.n_orbitals[
-                                      ik, ntoi[sp]]], numpy.complex_) for ik in range(self.n_k)]
+                                      ik, ntoi[sp]]], complex) for ik in range(self.n_k)]
 
         ikarray = numpy.array(list(range(self.n_k)))
         for ik in mpi.slice_array(ikarray):
@@ -2204,7 +2204,7 @@ class SumkDFT(object):
     def check_projectors(self):
         """Calculated the density matrix from projectors (DM = P Pdagger) to check that it is correct and
            specifically that it matches DFT."""
-        dens_mat = [numpy.zeros([self.corr_shells[icrsh]['dim'], self.corr_shells[icrsh]['dim']], numpy.complex_)
+        dens_mat = [numpy.zeros([self.corr_shells[icrsh]['dim'], self.corr_shells[icrsh]['dim']], complex)
                     for icrsh in range(self.n_corr_shells)]
 
         for ik in range(self.n_k):
diff --git a/python/triqs_dft_tools/sumk_dft_tools.py b/python/triqs_dft_tools/sumk_dft_tools.py
index e11f42d7..9bb5bcf0 100644
--- a/python/triqs_dft_tools/sumk_dft_tools.py
+++ b/python/triqs_dft_tools/sumk_dft_tools.py
@@ -100,16 +100,16 @@ class SumkDFTTools(SumkDFT):
         for icrsh in range(self.n_corr_shells):
             G_loc[icrsh].zero()
 
-        DOS = {sp: numpy.zeros([n_om], numpy.float_)
+        DOS = {sp: numpy.zeros([n_om], float)
                for sp in self.spin_block_names[self.SO]}
         DOSproj = [{} for ish in range(self.n_inequiv_shells)]
         DOSproj_orb = [{} for ish in range(self.n_inequiv_shells)]
         for ish in range(self.n_inequiv_shells):
             for sp in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]]['SO']]:
                 dim = self.corr_shells[self.inequiv_to_corr[ish]]['dim']
-                DOSproj[ish][sp] = numpy.zeros([n_om], numpy.float_)
+                DOSproj[ish][sp] = numpy.zeros([n_om], float)
                 DOSproj_orb[ish][sp] = numpy.zeros(
-                    [n_om, dim, dim], numpy.complex_)
+                    [n_om, dim, dim], complex)
 
         ikarray = numpy.array(list(range(self.n_k)))
         for ik in mpi.slice_array(ikarray):
@@ -240,16 +240,16 @@ class SumkDFTTools(SumkDFT):
                      for block, inner in gf_struct_parproj_all[0]]
         G_loc_all = BlockGf(name_list=spn, block_list=glist_all, make_copies=False)
 
-        DOS = {sp: numpy.zeros([n_om], numpy.float_)
+        DOS = {sp: numpy.zeros([n_om], float)
                for sp in self.spin_block_names[self.SO]}
         DOSproj = {}
         DOSproj_orb = {}
 
         for sp in self.spin_block_names[self.SO]:
             dim = n_local_orbs
-            DOSproj[sp] = numpy.zeros([n_om], numpy.float_)
+            DOSproj[sp] = numpy.zeros([n_om], float)
             DOSproj_orb[sp] = numpy.zeros(
-                    [n_om, dim, dim], numpy.complex_)
+                    [n_om, dim, dim], complex)
 
         ikarray = numpy.array(list(range(self.n_k)))
         for ik in mpi.slice_array(ikarray):
@@ -375,16 +375,16 @@ class SumkDFTTools(SumkDFT):
         for ish in range(self.n_shells):
             G_loc[ish].zero()
 
-        DOS = {sp: numpy.zeros([n_om], numpy.float_)
+        DOS = {sp: numpy.zeros([n_om], float)
                for sp in self.spin_block_names[self.SO]}
         DOSproj = [{} for ish in range(self.n_shells)]
         DOSproj_orb = [{} for ish in range(self.n_shells)]
         for ish in range(self.n_shells):
             for sp in self.spin_block_names[self.SO]:
                 dim = self.shells[ish]['dim']
-                DOSproj[ish][sp] = numpy.zeros([n_om], numpy.float_)
+                DOSproj[ish][sp] = numpy.zeros([n_om], float)
                 DOSproj_orb[ish][sp] = numpy.zeros(
-                    [n_om, dim, dim], numpy.complex_)
+                    [n_om, dim, dim], complex)
 
         ikarray = numpy.array(list(range(self.n_k)))
         for ik in mpi.slice_array(ikarray):
@@ -518,11 +518,11 @@ class SumkDFTTools(SumkDFT):
             om_maxplot = plot_range[1]
 
         if ishell is None:
-            Akw = {sp: numpy.zeros([self.n_k, n_om], numpy.float_)
+            Akw = {sp: numpy.zeros([self.n_k, n_om], float)
                    for sp in spn}
         else:
             Akw = {sp: numpy.zeros(
-                [self.shells[ishell]['dim'], self.n_k, n_om], numpy.float_) for sp in spn}
+                [self.shells[ishell]['dim'], self.n_k, n_om], float) for sp in spn}
 
         if not ishell is None:
             gf_struct_parproj = [
@@ -649,7 +649,7 @@ class SumkDFTTools(SumkDFT):
         spn = self.spin_block_names[self.SO]
         ntoi = self.spin_names_to_ind[self.SO]
         # Density matrix in the window
-        self.dens_mat_window = [[numpy.zeros([self.shells[ish]['dim'], self.shells[ish]['dim']], numpy.complex_)
+        self.dens_mat_window = [[numpy.zeros([self.shells[ish]['dim'], self.shells[ish]['dim']], complex)
                                  for ish in range(self.n_shells)]
                                 for isp in range(len(spn))]
         # Set up G_loc
@@ -921,7 +921,7 @@ class SumkDFTTools(SumkDFT):
             print("Omega mesh automatically repined to:  ", self.Om_mesh)
 
         self.Gamma_w = {direction: numpy.zeros(
-            (len(self.Om_mesh), n_om), dtype=numpy.float_) for direction in self.directions}
+            (len(self.Om_mesh), n_om), dtype=float) for direction in self.directions}
 
         # Sum over all k-points
         ikarray = numpy.array(list(range(self.n_k)))
@@ -929,7 +929,7 @@ class SumkDFTTools(SumkDFT):
             # Calculate G_w  for ik and initialize A_kw
             G_w = self.lattice_gf(ik, mu, iw_or_w="w", beta=beta,
                                   broadening=broadening, mesh=mesh, with_Sigma=with_Sigma)
-            A_kw = [numpy.zeros((self.n_orbitals[ik][isp], self.n_orbitals[ik][isp], n_om), dtype=numpy.complex_)
+            A_kw = [numpy.zeros((self.n_orbitals[ik][isp], self.n_orbitals[ik][isp], n_om), dtype=complex)
                     for isp in range(n_inequiv_spin_blocks)]
 
             for isp in range(n_inequiv_spin_blocks):
diff --git a/test/python/analyse_block_structure_from_gf.py b/test/python/analyse_block_structure_from_gf.py
index 328425c4..d90b5c17 100644
--- a/test/python/analyse_block_structure_from_gf.py
+++ b/test/python/analyse_block_structure_from_gf.py
@@ -105,7 +105,7 @@ else:
 
 for conjugate in conjugate_values:
     # construct a random block-diagonal Hloc
-    Hloc = np.zeros((10,10), dtype=np.complex_)
+    Hloc = np.zeros((10,10), dtype=complex)
     # the Hloc of the first three 2x2 blocks is equal
     Hloc0 = get_random_hermitian(2)
     Hloc[:2,:2] = Hloc0
diff --git a/test/python/analyse_block_structure_from_gf2.py b/test/python/analyse_block_structure_from_gf2.py
index 25806815..76e0aa71 100644
--- a/test/python/analyse_block_structure_from_gf2.py
+++ b/test/python/analyse_block_structure_from_gf2.py
@@ -20,7 +20,7 @@ def get_random_transformation(dim):
     return T
 
 # construct a random block-diagonal Hloc
-Hloc = np.zeros((10,10), dtype=np.complex_)
+Hloc = np.zeros((10,10), dtype=complex)
 # the Hloc of the first three 2x2 blocks is equal
 Hloc0 = get_random_hermitian(2)
 Hloc[:2,:2] = Hloc0
@@ -88,7 +88,7 @@ Gt = BlockGf(name_block_generator = [(name,
                 n_points=len(block.mesh),
                 indices=block.indices)) for name, block in G], make_copies=False)
 
-known_moments = np.zeros((2,10,10), dtype=np.complex)
+known_moments = np.zeros((2,10,10), dtype=complex)
 known_moments[1,:] = np.eye(10)
 tail, err = fit_tail(G['ud'], known_moments)
 Gt['ud'].set_from_fourier(G['ud'], tail)
diff --git a/test/python/plovasp/proj_group/test_block_map.py b/test/python/plovasp/proj_group/test_block_map.py
index 2c923226..b6622be8 100644
--- a/test/python/plovasp/proj_group/test_block_map.py
+++ b/test/python/plovasp/proj_group/test_block_map.py
@@ -29,7 +29,7 @@ class TestBlockMap(mytest.MyTestCase):
         self.mock_eigvals = np.zeros((1, 11, 1))
 
         nproj = 16
-        self.mock_plo = np.zeros((nproj, 1, 1, 11), dtype=np.complex128)
+        self.mock_plo = np.zeros((nproj, 1, 1, 11), dtype=complex)
         self.mock_proj_params = [{} for i in range(nproj)]
         ip = 0
 # Mock d-sites
diff --git a/test/python/plovasp/proj_group/test_one_site_compl.py b/test/python/plovasp/proj_group/test_one_site_compl.py
index 2df4d960..f7a71686 100644
--- a/test/python/plovasp/proj_group/test_one_site_compl.py
+++ b/test/python/plovasp/proj_group/test_one_site_compl.py
@@ -73,7 +73,7 @@ class TestProjectorGroupCompl(mytest.MyTestCase):
                 bmax = self.proj_gr.ib_win[ik, isp, 1]+1
 
                 nb = bmax - bmin
-                p_mat = np.zeros((ndim, nb), dtype=np.complex128)
+                p_mat = np.zeros((ndim, nb), dtype=complex)
                 #print(bmin,bmax,nb)
 # Combine all projectors of the group to one block projector
                 for bl_map in block_maps: