diff --git a/EPW/bin/allocate.py b/EPW/bin/allocate.py
index 2d24d2591..99587ac85 100644
--- a/EPW/bin/allocate.py
+++ b/EPW/bin/allocate.py
@@ -1,13 +1,15 @@
-# 
-# Script to automatically add status and error message to ALLOCATE and DEALLOCATE in 
-# Fortran files. 
+#
+# Script to automatically add status and error message to ALLOCATE and DEALLOCATE in
+# Fortran files.
 #
 # Author: S. Ponce
 # Date: Sept. 2019
 #
 
+from __future__ import print_function
 import numpy as np
 
+
 # File name
 file_name = 'elphel2_shuffle.f90'
 
@@ -18,7 +20,7 @@ with open(file_name,'r') as F:
     #tmp = lines.split()
     if len(tmp_split) < 1:
       continue
-    # 
+    #
     #print tmp_split[0][0:9]
     if tmp_split[0] == 'SUBROUTINE':
       for ii in np.arange(len(tmp)):
@@ -28,7 +30,7 @@ with open(file_name,'r') as F:
           end = ii
 
       name_sub = str(tmp[start+2:end])
-    
+
     if tmp_split[0] == 'FUNCTION':
       for ii in np.arange(len(tmp)):
         if tmp[ii] == 'N' and tmp[ii+1] == ' ':
@@ -47,19 +49,18 @@ with open(file_name,'r') as F:
         if tmp[ii] == ')':
           final = ii
       #print 'start ',start
-      print str(tmp[0:start-9])+'ALLOCATE('+str(tmp[start:final])+', STAT = ierr)'
-      print str(tmp[0:start-9])+'IF (ierr /= 0) CALL errore(\''+str(name_sub)+'\', \'Error allocating '+str(tmp[start:end])+'\', 1)'     
+      print('{}ALLOCATE({}, STAT = ierr)'.format(tmp[0:start-9], tmp[start:final]))
+      print("{}IF (ierr /= 0) CALL errore('{}', 'Error allocating {}', 1)".format(
+          tmp[0:start-9], name_sub, tmp[start:end]
+      ))
     elif tmp_split[0][0:11] == 'DEALLOCATE(':
       for ii in np.arange(len(tmp)):
         if tmp[ii] == '(':
           start = ii
         if tmp[ii] == ')':
           end = ii
-      print str(tmp[0:start-10])+'DEALLOCATE('+str(tmp[start+1:end])+', STAT = ierr)' 
-      print str(tmp[0:start-10])+'IF (ierr /= 0) CALL errore(\''+str(name_sub)+'\', \'Error deallocating '+str(tmp[start+1:end])+'\', 1)' 
+      print('{}DEALLOCATE({}, STAT = ierr)'.format(tmp[0:start-10], tmp[start+1:end]))
+      print("{}IF (ierr /= 0) CALL errore('{}', 'Error deallocating {}', 1)".format(
+          tmp[0:start-10], name_sub, tmp[start+1:end]))
     else:
-      print str(lines),
-
-
-
-
+      print(lines, end='')
diff --git a/EPW/bin/pp-xml-depreciated.py b/EPW/bin/pp-xml-depreciated.py
index 471e6346a..699359594 100644
--- a/EPW/bin/pp-xml-depreciated.py
+++ b/EPW/bin/pp-xml-depreciated.py
@@ -1,11 +1,13 @@
-#!/usr/bin/python
+#!/usr/bin/env python3
 #
 # Post-processing script from of PH data in format used by EPW
 # 14/07/2015 - Creation of the script - Samuel Ponce
 # 14/03/2018 - Automatically reads the number of q-points - Michael Waters
-# 14/03/2018 - Detect if SOC is included in the calculation - Samuel Ponce 
+# 14/03/2018 - Detect if SOC is included in the calculation - Samuel Ponce
 # 13/11/2018 - Write dyn files in xml format for SOC case - Shunhong Zhang (USTC)
-# 
+#
+from __future__ import print_function
+from builtins import input
 import numpy as np
 import os
 from xml.dom import minidom
@@ -14,12 +16,12 @@ from xml.dom import minidom
 def dyn2xml(prefix):
     ndyn=int(os.popen('head -2 {0}.dyn0|tail -1'.format(prefix)).read())
     for idyn in range(1,ndyn+1):
-        print '{0}.dyn{1} to {0}.dyn_q{1}.xml'.format(prefix,idyn)
+        print('{0}.dyn{1} to {0}.dyn_q{1}.xml'.format(prefix, idyn))
         dynmat=dyn(prefix,idyn)
         dynmat._write_xml()
 def get_geom_info():
-    if os.path.isfile('ph.out')==False:
-       print 'cannot extract geometry info from ph.out'
+    if not os.path.isfile('ph.out'):
+       print('cannot extract geometry info from ph.out')
        return 1
     else:
        volm=float(os.popen('grep -a volume ph.out 2>/dev/null|tail -1').readline().split()[-2])
@@ -47,7 +49,7 @@ class dyn(object):
         self._at=np.zeros((3,3),float)
         self._bg=np.zeros((3,3),float)
         try: self._volm,self._at,self._bg = get_geom_info()
-        except: print 'warning: lattice info not found'
+        except Exception: print('warning: lattice info not found')
         for i in range(0, 4):
             f.readline()
         self._species=[];
@@ -215,7 +217,7 @@ def hasSOC(prefix):
   xmldoc = minidom.parse(fname)
   item = xmldoc.getElementsByTagName('spinorbit')[0]
   lSOC = item.childNodes[0].data
-  
+
   return lSOC
 
 # Check if the calculation was done in sequential
@@ -225,11 +227,11 @@ def isSEQ(prefix):
     lseq = True
   else:
     lseq = False
- 
+
   return lseq
-    
+
 # Enter the number of irr. q-points
-user_input = raw_input('Enter the prefix used for PH calculations (e.g. diam)\n')
+user_input = input('Enter the prefix used for PH calculations (e.g. diam)\n')
 prefix = str(user_input)
 
 # Test if SOC
@@ -237,14 +239,14 @@ SOC = hasSOC(prefix)
 
 # If SOC detected, but dyn is not in XML and we want to convert it
 if SOC=='true':
-  user_input = raw_input('Calculation with SOC detected. Do you want to convert dyn in XML format [y/n]?\n')
+  user_input = input('Calculation with SOC detected. Do you want to convert dyn in XML format [y/n]?\n')
   if str(user_input) == 'y':
     dyn2xml(prefix)
     os.system('mv {0}.dyn*.xml save'.format(prefix))
 
-# If no SOC detected, do you want to convert into XML format 
+# If no SOC detected, do you want to convert into XML format
 if SOC=='false':
-  user_input = raw_input('Calculation without SOC detected. Do you want to convert to xml anyway [y/n]?\n')
+  user_input = input('Calculation without SOC detected. Do you want to convert to xml anyway [y/n]?\n')
   if str(user_input) == 'y':
     SOC = 'true'
     dyn2xml(prefix)
@@ -258,7 +260,7 @@ if True: # this gets the nqpt from the outputfiles
 
 else:
   # Enter the number of irr. q-points
-  user_input = raw_input('Enter the number of irreducible q-points\n')
+  user_input = input('Enter the number of irreducible q-points\n')
   nqpt = user_input
   try:
     nqpt = int(user_input)
@@ -293,7 +295,7 @@ for iqpt in np.arange(1,nqpt+1):
       else:
         os.system('cp _ph0/'+prefix+'.q_'+str(iqpt)+'/'+prefix+'.dvscf save/'+prefix+'.dvscf_q'+label)
         os.system('rm _ph0/'+prefix+'.q_'+str(iqpt)+'/*wfc*' )
- 
+
   else:
     # Case with SOC
     if SOC == 'true':
@@ -316,4 +318,3 @@ for iqpt in np.arange(1,nqpt+1):
       else:
         os.system('cp _ph0/'+prefix+'.q_'+str(iqpt)+'/'+prefix+'.dvscf1 save/'+prefix+'.dvscf_q'+label)
         os.system('rm _ph0/'+prefix+'.q_'+str(iqpt)+'/*wfc*' )
-
diff --git a/EPW/bin/pp.py b/EPW/bin/pp.py
index c353ee554..f61d7011d 100644
--- a/EPW/bin/pp.py
+++ b/EPW/bin/pp.py
@@ -1,4 +1,4 @@
-#!/usr/bin/python
+#!/usr/bin/env python3
 #
 # Post-processing script from of PH data in format used by EPW
 # 14/07/2015 - Creation of the script - Samuel Ponce
diff --git a/EPW/examples/diamond/epw/meshes/kgen.py b/EPW/examples/diamond/epw/meshes/kgen.py
index 3df9050b9..179533542 100644
--- a/EPW/examples/diamond/epw/meshes/kgen.py
+++ b/EPW/examples/diamond/epw/meshes/kgen.py
@@ -1,15 +1,16 @@
 #
 # 14/07/2015  Samuel Ponce
 #
+from __future__ import print_function
 import numpy as np
 
-i=0
-for ii in np.arange(0.5,0.0,-1.0/200):
-  print str(ii)+' 0.0 0.0 '+str(1.0/201)
-  i +=1
+i = 0
+for ii in np.arange(0.5, 0.0, -1.0 / 200):
+    print('{0} 0.0 0.0 {1}'.format(ii, 1.0 / 201))
+    i += 1
 
-for ii in np.arange(0.0,0.5+1.0/200,1.0/200):
-  print str(ii)+' '+str(ii)+' 0.0 '+str(1.0/201)
-  i +=1
+for ii in np.arange(0.0, 0.5 + 1.0 / 200, 1.0 / 200):
+    print('{0} {0} 0.0 {1}'.format(ii, 1.0 / 201))
+    i += 1
 
-print i
+print(i)
diff --git a/EPW/examples/diamond/phonons/pp.py b/EPW/examples/diamond/phonons/pp.py
index ca36eadc8..ae85d68c6 100755
--- a/EPW/examples/diamond/phonons/pp.py
+++ b/EPW/examples/diamond/phonons/pp.py
@@ -1,20 +1,19 @@
 #
-# Post-processing script QE --> EPW 
+# Post-processing script QE --> EPW
 # 14/07/2015 - Samuel Ponce
 #
 
+from builtins import input
 import numpy as np
 import os
 
-# Enter the number of irr. q-points 
-user_input = raw_input('Enter the prefix used for PH calculations (e.g. diam)\n')
-prefix = str(user_input)
+# Enter the number of irr. q-points
+prefix = input('Enter the prefix used for PH calculations (e.g. diam)\n')
 
-# Enter the number of irr. q-points 
-user_input = raw_input('Enter the number of irreducible q-points\n')
-nqpt = user_input
+# Enter the number of irr. q-points
+nqpt = input('Enter the number of irreducible q-points\n')
 try:
-  nqpt = int(user_input)
+  nqpt = int(nqpt)
 except ValueError:
   raise Exception('The value you enter is not an integer!')
 
@@ -30,4 +29,3 @@ for iqpt in np.arange(1,nqpt+1):
   else:
     os.system('cp _ph0/'+prefix+'.q_'+str(iqpt)+'/'+prefix+'.dvscf1 save/'+prefix+'.dvscf_q'+label)
     os.system('rm _ph0/'+prefix+'.q_'+str(iqpt)+'/*wfc*' )
-
diff --git a/EPW/examples/gan/phonons/pp.py b/EPW/examples/gan/phonons/pp.py
index ca36eadc8..ae85d68c6 100644
--- a/EPW/examples/gan/phonons/pp.py
+++ b/EPW/examples/gan/phonons/pp.py
@@ -1,20 +1,19 @@
 #
-# Post-processing script QE --> EPW 
+# Post-processing script QE --> EPW
 # 14/07/2015 - Samuel Ponce
 #
 
+from builtins import input
 import numpy as np
 import os
 
-# Enter the number of irr. q-points 
-user_input = raw_input('Enter the prefix used for PH calculations (e.g. diam)\n')
-prefix = str(user_input)
+# Enter the number of irr. q-points
+prefix = input('Enter the prefix used for PH calculations (e.g. diam)\n')
 
-# Enter the number of irr. q-points 
-user_input = raw_input('Enter the number of irreducible q-points\n')
-nqpt = user_input
+# Enter the number of irr. q-points
+nqpt = input('Enter the number of irreducible q-points\n')
 try:
-  nqpt = int(user_input)
+  nqpt = int(nqpt)
 except ValueError:
   raise Exception('The value you enter is not an integer!')
 
@@ -30,4 +29,3 @@ for iqpt in np.arange(1,nqpt+1):
   else:
     os.system('cp _ph0/'+prefix+'.q_'+str(iqpt)+'/'+prefix+'.dvscf1 save/'+prefix+'.dvscf_q'+label)
     os.system('rm _ph0/'+prefix+'.q_'+str(iqpt)+'/*wfc*' )
-
diff --git a/EPW/examples/mgb2/phonons/pp.py b/EPW/examples/mgb2/phonons/pp.py
index ca36eadc8..ae85d68c6 100755
--- a/EPW/examples/mgb2/phonons/pp.py
+++ b/EPW/examples/mgb2/phonons/pp.py
@@ -1,20 +1,19 @@
 #
-# Post-processing script QE --> EPW 
+# Post-processing script QE --> EPW
 # 14/07/2015 - Samuel Ponce
 #
 
+from builtins import input
 import numpy as np
 import os
 
-# Enter the number of irr. q-points 
-user_input = raw_input('Enter the prefix used for PH calculations (e.g. diam)\n')
-prefix = str(user_input)
+# Enter the number of irr. q-points
+prefix = input('Enter the prefix used for PH calculations (e.g. diam)\n')
 
-# Enter the number of irr. q-points 
-user_input = raw_input('Enter the number of irreducible q-points\n')
-nqpt = user_input
+# Enter the number of irr. q-points
+nqpt = input('Enter the number of irreducible q-points\n')
 try:
-  nqpt = int(user_input)
+  nqpt = int(nqpt)
 except ValueError:
   raise Exception('The value you enter is not an integer!')
 
@@ -30,4 +29,3 @@ for iqpt in np.arange(1,nqpt+1):
   else:
     os.system('cp _ph0/'+prefix+'.q_'+str(iqpt)+'/'+prefix+'.dvscf1 save/'+prefix+'.dvscf_q'+label)
     os.system('rm _ph0/'+prefix+'.q_'+str(iqpt)+'/*wfc*' )
-
diff --git a/EPW/examples/pb/wSOC/phonons/pp.py b/EPW/examples/pb/wSOC/phonons/pp.py
index ca36eadc8..ae85d68c6 100644
--- a/EPW/examples/pb/wSOC/phonons/pp.py
+++ b/EPW/examples/pb/wSOC/phonons/pp.py
@@ -1,20 +1,19 @@
 #
-# Post-processing script QE --> EPW 
+# Post-processing script QE --> EPW
 # 14/07/2015 - Samuel Ponce
 #
 
+from builtins import input
 import numpy as np
 import os
 
-# Enter the number of irr. q-points 
-user_input = raw_input('Enter the prefix used for PH calculations (e.g. diam)\n')
-prefix = str(user_input)
+# Enter the number of irr. q-points
+prefix = input('Enter the prefix used for PH calculations (e.g. diam)\n')
 
-# Enter the number of irr. q-points 
-user_input = raw_input('Enter the number of irreducible q-points\n')
-nqpt = user_input
+# Enter the number of irr. q-points
+nqpt = input('Enter the number of irreducible q-points\n')
 try:
-  nqpt = int(user_input)
+  nqpt = int(nqpt)
 except ValueError:
   raise Exception('The value you enter is not an integer!')
 
@@ -30,4 +29,3 @@ for iqpt in np.arange(1,nqpt+1):
   else:
     os.system('cp _ph0/'+prefix+'.q_'+str(iqpt)+'/'+prefix+'.dvscf1 save/'+prefix+'.dvscf_q'+label)
     os.system('rm _ph0/'+prefix+'.q_'+str(iqpt)+'/*wfc*' )
-
diff --git a/EPW/examples/pb/woSOC/phonons/pp.py b/EPW/examples/pb/woSOC/phonons/pp.py
index ca36eadc8..ae85d68c6 100644
--- a/EPW/examples/pb/woSOC/phonons/pp.py
+++ b/EPW/examples/pb/woSOC/phonons/pp.py
@@ -1,20 +1,19 @@
 #
-# Post-processing script QE --> EPW 
+# Post-processing script QE --> EPW
 # 14/07/2015 - Samuel Ponce
 #
 
+from builtins import input
 import numpy as np
 import os
 
-# Enter the number of irr. q-points 
-user_input = raw_input('Enter the prefix used for PH calculations (e.g. diam)\n')
-prefix = str(user_input)
+# Enter the number of irr. q-points
+prefix = input('Enter the prefix used for PH calculations (e.g. diam)\n')
 
-# Enter the number of irr. q-points 
-user_input = raw_input('Enter the number of irreducible q-points\n')
-nqpt = user_input
+# Enter the number of irr. q-points
+nqpt = input('Enter the number of irreducible q-points\n')
 try:
-  nqpt = int(user_input)
+  nqpt = int(nqpt)
 except ValueError:
   raise Exception('The value you enter is not an integer!')
 
@@ -30,4 +29,3 @@ for iqpt in np.arange(1,nqpt+1):
   else:
     os.system('cp _ph0/'+prefix+'.q_'+str(iqpt)+'/'+prefix+'.dvscf1 save/'+prefix+'.dvscf_q'+label)
     os.system('rm _ph0/'+prefix+'.q_'+str(iqpt)+'/*wfc*' )
-
diff --git a/EPW/examples/sic/phonons/pp.py b/EPW/examples/sic/phonons/pp.py
index ca36eadc8..ae85d68c6 100755
--- a/EPW/examples/sic/phonons/pp.py
+++ b/EPW/examples/sic/phonons/pp.py
@@ -1,20 +1,19 @@
 #
-# Post-processing script QE --> EPW 
+# Post-processing script QE --> EPW
 # 14/07/2015 - Samuel Ponce
 #
 
+from builtins import input
 import numpy as np
 import os
 
-# Enter the number of irr. q-points 
-user_input = raw_input('Enter the prefix used for PH calculations (e.g. diam)\n')
-prefix = str(user_input)
+# Enter the number of irr. q-points
+prefix = input('Enter the prefix used for PH calculations (e.g. diam)\n')
 
-# Enter the number of irr. q-points 
-user_input = raw_input('Enter the number of irreducible q-points\n')
-nqpt = user_input
+# Enter the number of irr. q-points
+nqpt = input('Enter the number of irreducible q-points\n')
 try:
-  nqpt = int(user_input)
+  nqpt = int(nqpt)
 except ValueError:
   raise Exception('The value you enter is not an integer!')
 
@@ -30,4 +29,3 @@ for iqpt in np.arange(1,nqpt+1):
   else:
     os.system('cp _ph0/'+prefix+'.q_'+str(iqpt)+'/'+prefix+'.dvscf1 save/'+prefix+'.dvscf_q'+label)
     os.system('rm _ph0/'+prefix+'.q_'+str(iqpt)+'/*wfc*' )
-
diff --git a/FFTXlib/gen_test_params.py b/FFTXlib/gen_test_params.py
index 9f6ff2e2d..6ebffbbfd 100644
--- a/FFTXlib/gen_test_params.py
+++ b/FFTXlib/gen_test_params.py
@@ -1,3 +1,4 @@
+from __future__ import print_function
 import re, sys
 
 info = """
diff --git a/PP/tools/sum_states.py b/PP/tools/sum_states.py
index 2db5528dd..7996def57 100755
--- a/PP/tools/sum_states.py
+++ b/PP/tools/sum_states.py
@@ -1,4 +1,4 @@
-#! /usr/bin/python
+#!/usr/bin/env python3
 
 ######  SUM STATES #######
 # Python script for summing and ploting the data from the Density Of States
@@ -31,17 +31,17 @@
 # ni.pdos_atm#1\(Ni\)_wfc#2\(d\) -t "Example PP/02" -xr -6 2
 #
 #
-# The procedure for obtaining the DOS files is explained 
+# The procedure for obtaining the DOS files is explained
 # i.e. in (espresso-dir)/PP/examples/example02/
-# 
+#
 # Author: Dr. Julen Larrucea
 #         University of Bremen,
 #         Bremen Centre for Computational Materials Science, HMI Group
 #         julenl [at] gmail.com  or larrucea [at] hmi.uni-bremen.de
 #
-#  This file is distributed under the terms of the GNU General Public 
+#  This file is distributed under the terms of the GNU General Public
 #  License. See the file `License'
-#  in the root directory of the present distribution, 
+#  in the root directory of the present distribution,
 #  or http://www.gnu.org/copyleft/gpl.txt .
 #######################
 
@@ -60,7 +60,7 @@ min_y,max_y="",""
 output_file_name="sum_dos.out"
 prt="no"
 
-print " #### sum_states.py version "+str(version)+" #### "  
+print(" #### sum_states.py version {} #### ".format(version))
 
 
 # Check if X11, mathplotlib and gnuplot are available
@@ -89,14 +89,14 @@ if len(sys.argv)>1:
    if option=="t":
     graphtitle= sys.argv[sys.argv.index('-t')+1]
    if option=="xr":
-    min_x,max_x= float(sys.argv[sys.argv.index('-xr')+1]),float(sys.argv[sys.argv.index('-xr')+2]) 
+    min_x,max_x= float(sys.argv[sys.argv.index('-xr')+1]),float(sys.argv[sys.argv.index('-xr')+2])
    if option=="yr":
-    min_y,max_y= float(sys.argv[sys.argv.index('-yr')+1]),float(sys.argv[sys.argv.index('-yr')+2]) 
+    min_y,max_y= float(sys.argv[sys.argv.index('-yr')+1]),float(sys.argv[sys.argv.index('-yr')+2])
    if option=="v":
-    print "sum_dos.py version: "+version
+    print("sum_dos.py version:", version)
     sys.exit()
    if option=="h":
-    print '''
+    print('''
     -o QE output file name (for grepping Fermi E)
     -s Selection of atoms for summing the DOSes. "*" for all, *1*Fe*d* for first Fe atom " (def. "*")
     -p Print output to a file and aditionaly provide an output name (def. no output and "sum_dos.out")
@@ -105,8 +105,8 @@ if len(sys.argv)>1:
     -yr set min and max y value for the axes in the graph
     -h print this help
     -v print version
-       Example: sum_states.py --s sys.pdos_atm#4\(Fe2\)_wfc#2\(d\) -t "Wustite LDA+U single Fe" -xr -9 4 
-   '''
+       Example: sum_states.py --s sys.pdos_atm#4\(Fe2\)_wfc#2\(d\) -t "Wustite LDA+U single Fe" -xr -9 4
+   ''')
     sys.exit()
 
 
@@ -116,13 +116,13 @@ if len(os.popen('echo $DISPLAY').read()) > 1:
   try:
    from pylab import *
    mplplot="yes"
-   print "pylab imported"
+   print("pylab imported")
   except:
-   print "There is no mathplotlib installed. Using gnuplot."
+   print("There is no mathplotlib installed. Using gnuplot.")
    mplplot="no"
    prt="yes"
 else:
-  print "No X11. Trying to plot on terminal"
+  print("No X11. Trying to plot on terminal")
   graphic_plot="no"
   if prog_gnuplot=="no":
    prt="yes"
@@ -132,7 +132,7 @@ else:
 if pwout == "":
  for filen in filter(os.path.isfile, os.listdir('.')):
   if "Program PWSCF" in linecache.getline(filen, 2):
-   print "Using " + filen + " as pw.x output. You can specify another one with the -o option."
+   print("Using " + filen + " as pw.x output. You can specify another one with the -o option.")
    pwout=filen
 
 # Parse Fermi energy from the pw.x output
@@ -140,35 +140,32 @@ if pwout!="":
  try:
   os.popen("grep -a 'the Fermi energy is' "+pwout ).read()
   fermi=float(os.popen("grep -a 'the Fermi energy is' "+pwout ).read().split()[4])
-  print "Fermi energy = ", fermi, "a.u."
+  print("Fermi energy = ", fermi, "a.u.")
  except:
-  print "WARNING: No Fermi energy found. Using 0 e.V. instead"
+  print("WARNING: No Fermi energy found. Using 0 e.V. instead")
   fermi=0
 else:
- print "WARNING: No pw.x output found. Using E Fermi = 0 e.V."
+ print("WARNING: No pw.x output found. Using E Fermi = 0 e.V.")
  fermi=0
 
-# List of all DOS files to add 
+# List of all DOS files to add
 dosfiles=[]
 for dfile in os.listdir('.'):
    if fnmatch.fnmatch(dfile, selat):
-     dosfiles.append(dfile) 
+     dosfiles.append(dfile)
 if len(dosfiles)==0:
- print "ERROR: Provide a (list of) valid DOS file(s)"
+ print("ERROR: Provide a (list of) valid DOS file(s)")
  sys.exit()
 
-print "dosfiles list: ",
-for dosfile in dosfiles:
-  print dosfile,
-print ""
+print("dosfiles list:", " ".join(dosfiles))
 
 # Check wetter we have k-solved DOS
 if open(dosfiles[0],'r').readline().split()[1]=="E":
   ksolved="no"
-  print "no ksolved"
+  print("no ksolved")
 elif open(dosfiles[0],'r').readline().split()[1]=="ik":
   ksolved="yes"
-  print "ksolved"
+  print("ksolved")
 
 # Sum over all k-points and files
 mat=[]  # matrix with total sum of ldos
@@ -191,28 +188,28 @@ for i in range(len(dosfiles)):
       elif ik == k and k > 1:
          oldmat.append([float(line.split()[1]),float(line.split()[2]),float(line.split()[3])])
       elif len(line) < 5 and k > 1:  #if blank line, sum k-frame to the total
-       for j in range(len(oldmat)):  
+       for j in range(len(oldmat)):
          mati[j]=[mati[j][0],mati[j][1]+oldmat[j][1],mati[j][2]+oldmat[j][2]]
 
  if mat == []: # if it is the first dos file, copy total matrix (mat) = the first dos files's data
     mat=mati[:]
  else:
     for j in range(len(mati)): # if it is not the first file, sum values
-        mat[j]=[mat[j][0],mat[j][1]+mati[j][1],mat[j][2]+mati[j][2]]  
+        mat[j]=[mat[j][0],mat[j][1]+mati[j][1],mat[j][2]+mati[j][2]]
 
 
-print "...ploting..."
+print("...ploting...")
 
 
 if prt=="yes":
  out=open(output_file_name,"w")
-x,y1,y2=[],[],[]    
+x,y1,y2=[],[],[]
 for i in mat:
  x.append(i[0]-fermi)
  y1.append(i[1])
  y2.append(-i[2])
  if prt=="yes":  # print to a file
-  print>>out, i[0]-fermi, i[1], i[2]
+  print(i[0] - fermi, i[1], i[2], file=out)
 if prt=="yes":
  out.close()
 
@@ -222,7 +219,7 @@ if graphic_plot=="yes":
   if mplplot=="yes":
     plot(x,y1,linewidth=1.0)
     plot(x,y2,linewidth=1.0)
-    print min(y2),max(y1)
+    print(min(y2), max(y1))
     plt.title(graphtitle)
     plt.xlabel('E (eV)')
     plt.ylabel('States')
@@ -231,15 +228,13 @@ if graphic_plot=="yes":
     plt.fill(x,y1,color='0.8')
     plt.fill(x,y2,color='0.9')
     if min_x and max_x:
-     fromx,tox=min_x,max_x 
+      fromx, tox = min_x, max_x
     plt.axis([fromx, tox, min(y2), max(y1)])
-    show()   
+    show()
   elif mplplot=="no" and prog_gnuplot=="yes":  # If no mathplotlib available, use gnuplot
-     os.system("echo \"plot '"+ output_file_name + "' using ($1-"+str(fermi)+"):2 w l, '' u ($1"+str(fermi)+"):3 w l\" | gnuplot -persist")  
+    os.system("""echo "plot '{0}' using ($1-{1}):2 w l, '' u ($1-{1}):3 w l" | gnuplot -persist""".format(
+        output_file_name, fermi))
 elif graphic_plot=="no":  # If no X forwarding available, show graph in terminal
   if prog_gnuplot=="yes":
-     os.system("echo \"set terminal dumb; plot '"+ output_file_name + "' using ($1-"+str(fermi)+"):2 w l, '' u ($1-"+str(fermi)+"):3 w l\" | gnuplot -persist")
-
-
-
-
+    os.system("""echo "set terminal dumb; plot '{0}' using ($1-{1}):2 w l, '' u ($1-{1}):3 w l" | gnuplot -persist""".format(
+        output_file_name, fermi))
diff --git a/dev-tools/get_device_props.py b/dev-tools/get_device_props.py
index 44a3f28cf..562e7c9cb 100755
--- a/dev-tools/get_device_props.py
+++ b/dev-tools/get_device_props.py
@@ -1,20 +1,25 @@
+from __future__ import print_function
 import sys
 import subprocess as sp
 
 have_yaml = True
 try:
     import yaml
-except:
-    have_yaml = False
-    
-try:
-    from StringIO import StringIO
 except ImportError:
-    from io import BytesIO as StringIO
+    have_yaml = False
+
+try:
+    from io import StringIO
+except ImportError:
+    try:
+        from StringIO import StringIO
+    except ImportError:
+        from io import BytesIO as StringIO
+
 
 print("\n\n This is a helper tool to check the details of your GPU before configuring QE.\n\n")
-print("""Remeber to load CUDA environemt and run this on the COMPUTE NODE 
-if you are not sure that the frontend node shares 
+print("""Remeber to load CUDA environemt and run this on the COMPUTE NODE
+if you are not sure that the frontend node shares
 the same configuration as the backend nodes\n\n""")
 
 
@@ -52,13 +57,11 @@ else:
                 minor = str(5 if int(minor)>=5 else 0)
             if 'major' in line:
                 _, major = line.split(':')
-            
+
             if minor != "" and major != "":
                 print("Compute capabilities for dev {}: {}.{}".format(str(devnum), major.strip(),minor.strip()))
                 conf_cc=major.strip()+minor.strip()
                 minor = ""; major = ""; devnum += 1
-                
+
     print("\n If all compute capabilities match, configure QE with:")
     print("./configure --with-cuda=$CUDA_HOME --with-cuda-cc={} --with-cuda-runtime={}\n".format(conf_cc, conf_rt))
-
-    
diff --git a/dev-tools/mem_analyse.py b/dev-tools/mem_analyse.py
index 71e27c8d1..0b3a5bb2b 100755
--- a/dev-tools/mem_analyse.py
+++ b/dev-tools/mem_analyse.py
@@ -1,8 +1,9 @@
+from __future__ import print_function
 import numpy as np
 
 #1) Run dev-tools/mem_counter inside EPW/src
 #2) compile UtilXlib/mem_counter.f90  with -D__DEBUG flag
-#3) Run EPW 
+#3) Run EPW
 #4) grep ' allocating' epw1.out > alloc.txt
 #5) grep 'deallocating' epw1.out > dealloc.txt
 #6) Run this script after having changed the correct allocation lengths
@@ -37,18 +38,18 @@ with open('dealloc.txt','r') as R:
 deall_found = [ False ] * dealloc_len
 
 for ii in np.arange(alloc_len):
-  print ii,' / ',alloc_len
+  print(ii, ' / ', alloc_len)
   name = alloc_name[ii]
-  
+
   found = False
   for jj in np.arange(dealloc_len):
     if name == dealloc_name[jj]:
-      if (alloc_size[ii] == dealloc_size[jj] and deall_found[jj]==False ) :
+      if alloc_size[ii] == dealloc_size[jj] and not deall_found[jj]:
         # We found the corresponding all/deall pair
         deall_found[jj] = True
         found = True
         break
-  if (found == False):   
+  if not found:
     with open('mem_analyse.out','a') as O:
       O.write('We did not find a maching pair in '+str(alloc_sub[ii])+'\n')
       O.write('Allocate:   '+str(name)+' '+str(alloc_size[ii])+'\n')
@@ -56,6 +57,3 @@ for ii in np.arange(alloc_len):
 #    print 'We did not find a maching pair in ', alloc_sub[ii]
 #    print 'Allocate:   ',name,' ',alloc_size[ii]
 #    print 'Deallocate: ',dealloc_name[jj],' ',dealloc_size[jj]
-
-
-
diff --git a/dev-tools/mem_counter.py b/dev-tools/mem_counter.py
index 855c673b0..2b967adec 100755
--- a/dev-tools/mem_counter.py
+++ b/dev-tools/mem_counter.py
@@ -24,6 +24,7 @@ for file in *.bkp; do cp $file ${file%.bkp}; done
 
 """
 
+from __future__ import print_function
 import os, sys, re
 from shutil import copyfile
 
diff --git a/dev-tools/src-normal.py b/dev-tools/src-normal.py
index 53ea51302..b57ca7ca2 100755
--- a/dev-tools/src-normal.py
+++ b/dev-tools/src-normal.py
@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/env python3
 
 # (C) 2010 Norbert Nemec
 #
diff --git a/test-suite/buildbot/Cineca_farm/check_smtp.py b/test-suite/buildbot/Cineca_farm/check_smtp.py
index e9b4d6c44..8ad7c7fdc 100644
--- a/test-suite/buildbot/Cineca_farm/check_smtp.py
+++ b/test-suite/buildbot/Cineca_farm/check_smtp.py
@@ -1,9 +1,9 @@
-#! /usr/bin/python -tt
+#!/usr/bin/env python3
 
 from __future__ import absolute_import
 from __future__ import division
 from __future__ import print_function
-
+from builtins import input
 from getpass import getpass
 from smtplib import SMTP
 
@@ -14,7 +14,6 @@ name and password would allow to send e-mail through it.
 """
 
 
-
 SMTP_HOST = 'smtp.gmail.com:587'
 LOCAL_HOST = ''
 
@@ -28,12 +27,12 @@ def main():
 
 #    server.ehlo()
 #    print(server.ehlo())
- 
+
     server.starttls()
 
     print(server.ehlo(LOCAL_HOST))
 
-    user = raw_input('user: ')
+    user = input('user: ')
     password = getpass('password: ')
 
     print(server.login(user, password))
@@ -41,17 +40,17 @@ def main():
     fromaddr = 'testfarmqef@gmail.com'
     toaddrs  = 'samuel.pon@gmail.com'
     msg = "\r\n".join([
-      "From: testfarmqef@gmail.com",
-      "To: samuel.pon@gmail.com",
-      "Subject: Buildbot",
-      "",
-      "Why, oh why"
-      ])
+        "From: testfarmqef@gmail.com",
+        "To: samuel.pon@gmail.com",
+        "Subject: Buildbot",
+        "",
+        "Why, oh why"
+    ])
 
-
-    server.sendmail(fromaddr, toaddrs, msg) 
+    server.sendmail(fromaddr, toaddrs, msg)
     server.close()
 
+
 if __name__ == '__main__':
     main()
 
diff --git a/test-suite/epw_mob/pp.py b/test-suite/epw_mob/pp.py
index cc7dd82f9..5cd57ca96 100644
--- a/test-suite/epw_mob/pp.py
+++ b/test-suite/epw_mob/pp.py
@@ -1,20 +1,19 @@
 #
-# Post-processing script QE --> EPW 
+# Post-processing script QE --> EPW
 # 14/07/2015 - Samuel Ponce
 #
 
+from builtins import input
 import numpy as np
 import os
 
-# Enter the number of irr. q-points 
-user_input = raw_input('Enter the prefix used for PH calculations (e.g. diam)\n')
-prefix = str(user_input)
+# Enter the number of irr. q-points
+prefix = input('Enter the prefix used for PH calculations (e.g. diam)\n')
 
-# Enter the number of irr. q-points 
-user_input = raw_input('Enter the number of irreducible q-points\n')
-nqpt = user_input
+# Enter the number of irr. q-points
+nqpt = input('Enter the number of irreducible q-points\n')
 try:
-  nqpt = int(user_input)
+  nqpt = int(nqpt)
 except ValueError:
   raise Exception('The value you enter is not an integer!')
 
@@ -30,5 +29,3 @@ for iqpt in np.arange(1,nqpt+1):
   else:
     os.system('cp _ph0/'+prefix+'.q_'+str(iqpt)+'/'+prefix+'.dvscf1 save/'+prefix+'.dvscf_q'+label)
     os.system('rm -f _ph0/'+prefix+'.q_'+str(iqpt)+'/*wfc*')
-
-
diff --git a/test-suite/not_epw_comp/grid.py b/test-suite/not_epw_comp/grid.py
index ae646f342..7f0916f13 100644
--- a/test-suite/not_epw_comp/grid.py
+++ b/test-suite/not_epw_comp/grid.py
@@ -1,11 +1,12 @@
+from __future__ import print_function
 import numpy as np
 
-tot = 0  
+tot = 0
 
-for ii in np.arange(0,1,0.25):
-  for jj in np.arange(0,1,0.25):
-    for kk in np.arange(0,1,0.25):
-      print ii,' ',jj,' ',kk,' ',1.0/64
+for ii in np.arange(0, 1, 0.25):
+  for jj in np.arange(0, 1, 0.25):
+    for kk in np.arange(0, 1, 0.25):
+      print(ii, ' ', jj, ' ', kk, ' ', 1.0 / 64)
       tot += 1
 
-print tot
+print(tot)
diff --git a/test-suite/not_epw_comp/pp.py b/test-suite/not_epw_comp/pp.py
index 8bcaadbdb..db6322b50 100644
--- a/test-suite/not_epw_comp/pp.py
+++ b/test-suite/not_epw_comp/pp.py
@@ -1,20 +1,19 @@
 #
-# Post-processing script QE --> EPW 
+# Post-processing script QE --> EPW
 # 14/07/2015 - Samuel Ponce
 #
 
+from builtins import input
 import numpy as np
 import os
 
-# Enter the number of irr. q-points 
-user_input = raw_input('Enter the prefix used for PH calculations (e.g. diam)\n')
-prefix = str(user_input)
+# Enter the number of irr. q-points
+prefix = input('Enter the prefix used for PH calculations (e.g. diam)\n')
 
-# Enter the number of irr. q-points 
-user_input = raw_input('Enter the number of irreducible q-points\n')
-nqpt = user_input
+# Enter the number of irr. q-points
+nqpt = input('Enter the number of irreducible q-points\n')
 try:
-  nqpt = int(user_input)
+  nqpt = int(nqpt)
 except ValueError:
   raise Exception('The value you enter is not an integer!')
 
@@ -28,4 +27,3 @@ for iqpt in np.arange(1,nqpt+1):
   else:
     os.system('cp _ph0/'+prefix+'.q_'+str(iqpt)+'/'+prefix+'.dvscf1 save/'+prefix+'.dvscf_q'+label)
     os.system('rm -f _ph0/'+prefix+'.q_'+str(iqpt)+'/*wfc*')
-
diff --git a/test-suite/testcode/bin/testcode.py b/test-suite/testcode/bin/testcode.py
index ed474227f..93517712c 100755
--- a/test-suite/testcode/bin/testcode.py
+++ b/test-suite/testcode/bin/testcode.py
@@ -26,6 +26,7 @@ documentation for further details.'''
 # copyright: (c) 2012 James Spencer
 # license: modified BSD; see LICENSE for more details
 
+from __future__ import print_function
 import glob
 import optparse
 import os
diff --git a/test-suite/testcode/lib/testcode2/__init__.py b/test-suite/testcode/lib/testcode2/__init__.py
index acf422121..b823c058b 100644
--- a/test-suite/testcode/lib/testcode2/__init__.py
+++ b/test-suite/testcode/lib/testcode2/__init__.py
@@ -8,6 +8,7 @@ A framework for regression testing numerical programs.
 :license: modified BSD; see LICENSE for more details.
 '''
 
+from __future__ import print_function
 import glob
 import os
 import pipes
diff --git a/test-suite/testcode/lib/testcode2/config.py b/test-suite/testcode/lib/testcode2/config.py
index b4d0ebfdc..4ac1cfdaa 100644
--- a/test-suite/testcode/lib/testcode2/config.py
+++ b/test-suite/testcode/lib/testcode2/config.py
@@ -8,6 +8,7 @@ Parse jobconfig and userconfig ini files.
 :license: modified BSD; see LICENSE for more details.
 '''
 
+from __future__ import print_function
 import copy
 import glob
 import os
diff --git a/test-suite/testcode/lib/testcode2/exceptions.py b/test-suite/testcode/lib/testcode2/exceptions.py
index 11652169e..ea356dc23 100644
--- a/test-suite/testcode/lib/testcode2/exceptions.py
+++ b/test-suite/testcode/lib/testcode2/exceptions.py
@@ -8,6 +8,7 @@ Custom exceptions.  Initialise signal handler for the interrupt signal.
 :license: modified BSD; see LICENSE for more details.
 '''
 
+from __future__ import print_function
 import signal
 import sys
 
diff --git a/test-suite/testcode/lib/testcode2/validation.py b/test-suite/testcode/lib/testcode2/validation.py
index 08abb4d88..ad9b0c535 100644
--- a/test-suite/testcode/lib/testcode2/validation.py
+++ b/test-suite/testcode/lib/testcode2/validation.py
@@ -8,6 +8,7 @@ Classes and functions for comparing data.
 :license: modified BSD; see LICENSE for more details.
 '''
 
+from __future__ import print_function
 import re
 import sys
 import warnings
diff --git a/upflib/fixfiles.py b/upflib/fixfiles.py
index ec4f8e06b..5a708a83d 100755
--- a/upflib/fixfiles.py
+++ b/upflib/fixfiles.py
@@ -4,6 +4,7 @@ simple script to mark as CDATA all text contained in PP_INPUTFILE element. It
 avoids any issue related to the presence on XML reserved characters in that
 section.  
 """
+from __future__ import print_function
 from subprocess import Popen, PIPE
 from xml.etree import cElementTree as eT
 import os