2to3_changes.log

Disabled circleCI
RefactoringTool: Skipping optional fixer: buffer
RefactoringTool: Skipping optional fixer: idioms
RefactoringTool: Skipping optional fixer: set_literal
RefactoringTool: Skipping optional fixer: ws_comma
RefactoringTool: No changes to ./setup.py
RefactoringTool: No changes to ./doc/create_man.py
RefactoringTool: No changes to ./doc/help.py
RefactoringTool: No changes to ./doc/help_list.py
RefactoringTool: Refactored ./examples/ASE/W2-PIMD-MP2/run-ase.py
RefactoringTool: Refactored ./examples/cp2k/nvt-cl/idtau_plot.py
RefactoringTool: Refactored ./examples/hswfqmc/FixWF_1Bead/intau_plot.py
RefactoringTool: Refactored ./examples/hswfqmc/WFOpt_2Bead/intau_plot.py
RefactoringTool: Refactored ./examples/lammps/h2o-planetary-64/estmod_dipole.py
RefactoringTool: Refactored ./examples/lammps/isof-vapor/process_out.py
RefactoringTool: Refactored ./examples/lammps/isof-water/process_out.py
RefactoringTool: No changes to ./examples/lammps/isofsc-vapor/process_dif.py
RefactoringTool: Refactored ./examples/lammps/isofsc-vapor/process_out.py
RefactoringTool: Refactored ./examples/lammps/isofsc-water/process_out.py
--- ./examples/ASE/W2-PIMD-MP2/run-ase.py 	(original)
+++ ./examples/ASE/W2-PIMD-MP2/run-ase.py 	(refactored)
@@ -14,17 +14,17 @@
     nthreads = int(sys.argv[1])
     dir_calc = sys.argv[2]
 except IndexError:
-    print 'command line arguments: <number of Gaussian threads> <calculator directory name>'
+    print('command line arguments: <number of Gaussian threads> <calculator directory name>')
     sys.exit(1)
 
 
 # print some info
-print 'Gaussian ASE i-PI socket runner'
-print '-------------------------------'
-print
-print '   number of threads: {:d}'.format(nthreads)
-print 'calculator directory: {:s}'.format(dir_calc)
-print
+print('Gaussian ASE i-PI socket runner')
+print('-------------------------------')
+print()
+print('   number of threads: {:d}'.format(nthreads))
+print('calculator directory: {:s}'.format(dir_calc))
+print()
 
 # fast scratch directory for Gaussian
 os.environ['GAUSS_SCRDIR'] = '/dev/shm'
--- ./examples/cp2k/nvt-cl/idtau_plot.py	(original)
+++ ./examples/cp2k/nvt-cl/idtau_plot.py	(refactored)
@@ -10,7 +10,7 @@
 
 nsteps = len(idtdata)
 maxt = nsteps * 0.1
-print("Mean IDTAU for " + str(maxt * 0.5) + "<t<" + str(maxt) + " fs: " + str(mean(idtdata[nsteps // 2:nsteps])))
+print(("Mean IDTAU for " + str(maxt * 0.5) + "<t<" + str(maxt) + " fs: " + str(mean(idtdata[nsteps // 2:nsteps]))))
 
 fig = figure()
 ax = fig.add_subplot(111)
--- ./examples/hswfqmc/FixWF_1Bead/intau_plot.py	(original)
+++ ./examples/hswfqmc/FixWF_1Bead/intau_plot.py	(refactored)
@@ -9,7 +9,7 @@
 
 nsteps = len(invdata)
 maxt = nsteps * 0.1
-print("Mean INTAU for " + str(maxt * 0.5) + "<t<" + str(maxt) + " fs: " + str(mean(invdata[nsteps // 2:nsteps])))
+print(("Mean INTAU for " + str(maxt * 0.5) + "<t<" + str(maxt) + " fs: " + str(mean(invdata[nsteps // 2:nsteps]))))
 
 fig = figure()
 ax = fig.add_subplot(111)
--- ./examples/hswfqmc/WFOpt_2Bead/intau_plot.py	(original)
+++ ./examples/hswfqmc/WFOpt_2Bead/intau_plot.py	(refactored)
@@ -9,7 +9,7 @@
 
 nsteps = len(invdata)
 maxt = nsteps * 0.1
-print("Mean INTAU for " + str(maxt * 0.5) + "<t<" + str(maxt) + " fs: " + str(mean(invdata[nsteps // 2:nsteps])))
+print(("Mean INTAU for " + str(maxt * 0.5) + "<t<" + str(maxt) + " fs: " + str(mean(invdata[nsteps // 2:nsteps]))))
 
 fig = figure()
 ax = fig.add_subplot(111)
--- ./examples/lammps/h2o-planetary-64/estmod_dipole.py	(original)
+++ ./examples/lammps/h2o-planetary-64/estmod_dipole.py	(refactored)
@@ -23,9 +23,9 @@
     nmolec = natoms / 3
     ans = np.zeros((npl, 3))
     x = np.zeros((nmolec, 3, 3))
-    for i in xrange(npl):
+    for i in range(npl):
         x[:] = qsum[:, i].reshape(nmolec, 3, 3)
-        for j in xrange(nmolec):
+        for j in range(nmolec):
             ans[i] += (gam * x[j, 0, :] + 0.5 * (1.0 - gam) * (x[j, 1, :] + x[j, 2, :])) * qm + \
                 x[j, 1, :] * qh + x[j, 2, :] * qh
     return ans
--- ./examples/lammps/isof-vapor/process_out.py	(original)
+++ ./examples/lammps/isof-vapor/process_out.py	(refactored)
@@ -20,7 +20,7 @@
 table = np.loadtxt('simulation.out')
 
 (asize, bsize) = np.shape(table)
-print asize, bsize
+print(asize, bsize)
 
 out = np.zeros((asize - equilibriumstep, 7))
 
--- ./examples/lammps/isof-water/process_out.py	(original)
+++ ./examples/lammps/isof-water/process_out.py	(refactored)
@@ -20,7 +20,7 @@
 table = np.loadtxt('simulation.out')
 
 (asize, bsize) = np.shape(table)
-print asize, bsize
+print(asize, bsize)
 
 out = np.zeros((asize - equilibriumstep, 7))
 
--- ./examples/lammps/isofsc-vapor/process_out.py	(original)
+++ ./examples/lammps/isofsc-vapor/process_out.py	(refactored)
@@ -20,7 +20,7 @@
 table = np.loadtxt('simulation.out')
 
 (asize, bsize) = np.shape(table)
-print asize, bsize
+print(asize, bsize)
 
 out = np.zeros((asize - equilibriumstep - 1, 7))
 
--- ./examples/lammps/isofsc-water/process_out.py	(original)
+++ ./examples/lammps/isofsc-water/process_out.py	(refactored)
@@ -20,7 +20,7 @@
 table = np.loadtxt('simulation.out')
 
 (asize, bsize) = np.shape(table)
-print asize, bsize
+print(asize, bsize)
 
 out = np.zeros((asize - equilibriumstep RefactoringTool: No changes to ./examples/yaff/mil53_ffsocket/run.py
RefactoringTool: Refactored ./examples/yaff/mil53_ffsocket/yaffdriver.py
RefactoringTool: No changes to ./ipi/__init__.py
RefactoringTool: No changes to ./ipi/engine/__init__.py
RefactoringTool: No changes to ./ipi/engine/atoms.py
RefactoringTool: No changes to ./ipi/engine/barostats.py
RefactoringTool: Refactored ./ipi/engine/beads.py
RefactoringTool: No changes to ./ipi/engine/cell.py
RefactoringTool: Refactored ./ipi/engine/ensembles.py
RefactoringTool: Refactored ./ipi/engine/forcefields.py
RefactoringTool: Refactored ./ipi/engine/forces.py
- 1, 7))
 
--- ./examples/yaff/mil53_ffsocket/yaffdriver.py	(original)
+++ ./examples/yaff/mil53_ffsocket/yaffdriver.py	(refactored)
@@ -38,7 +38,7 @@
     def await_header(self):
         data = self.s.recv(HDRLEN * L_CHAR)
         if self.verbose:
-            print 'RECV', (data,)
+            print('RECV', (data,))
         with log.section('DRIVER'):
             log("RECV %s %s" % (data, datetime.now()))
         return data
@@ -49,21 +49,21 @@
             stub = self.s.recv(size - len(data))
             data += stub
             if self.verbose:
-                print 'RECV %i/%i' % (len(data), size)  # , hexlify(stub)
+                print('RECV %i/%i' % (len(data), size))  # , hexlify(stub)
         return data
 
     def send_header(self, header):
         assert len(header) <= 12
         header = header.ljust(12, ' ')
         if self.verbose:
-            print 'SEND', (header,)
+            print('SEND', (header,))
         with log.section('DRIVER'):
             log("SEND %s %s" % (header, datetime.now()))
         self.s.send(header)
 
     def send_data(self, data):
         if self.verbose:
-            print 'SEND', len(data)  # , hexlify(data)
+            print('SEND', len(data))  # , hexlify(data)
         self.s.send(data)
 
 
@@ -155,4 +155,4 @@
                         log("EXIT %s" % datetime.now())
                 break
             else:
-                raise NotImplementedError, "Received unknown message %s" % header
+                raise NotImplementedError("Received unknown message %s" % header)
--- ./ipi/engine/beads.py	(original)
+++ ./ipi/engine/beads.py	(refactored)
@@ -233,7 +233,7 @@
             else:
                 dq = q[b, :] - q[self.nbeads - 1, :]
             epath += np.dot(dq, m * dq)
-        print "WARNING: RETURNS AN INCORRECT RESULT IF OPEN PATHS ARE BEING USED. CALL NM.VSPRING INSTEAD!!"
+        print("WARNING: RETURNS AN INCORRECT RESULT IF OPEN PATHS ARE BEING USED. CALL NM.VSPRING INSTEAD!!")
         return epath * 0.5
 
     def get_fpath(self):
--- ./ipi/engine/ensembles.py	(original)
+++ ./ipi/engine/ensembles.py	(refactored)
@@ -167,7 +167,7 @@
         dself.hweights = depend_array(name="hweights", value=np.asarray(self.hweights))
 
         # we use ScaledForceComponents to replicate the physical forces without (hopefully) them being actually recomputed
-        for ic in xrange(len(self.forces.mforces)):
+        for ic in range(len(self.forces.mforces)):
             sfc = ScaledForceComponent(self.forces.mforces[ic], 1.0)
             self.bias.add_component(self.forces.mbeads[ic], self.forces.mrpc[ic], sfc)
             dd(sfc).scaling._func = lambda i=ic: self.hweights[i] - 1
--- ./ipi/engine/forcefields.py	(original)
+++ ./ipi/engine/forcefields.py	(refactored)
@@ -128,7 +128,7 @@
         par_str = " "
 
         if not self.pars is None:
-            for k, v in self.pars.items():
+            for k, v in list(self.pars.items()):
                 par_str += k + " : " + str(v) + " , "
         else:
             par_str = " "
@@ -216,8 +216,8 @@
                 try:
                     self.requests.remove(request)
                 except ValueError:
-                    print "failed removing request", id(request), ' ',
-                    print [id(r) for r in self.requests], "@", threading.currentThread()
+                    print("failed removing request", id(request), ' ', end=' ')
+                    print([id(r) for r in self.requests], "@", threading.currentThread())
                     raise
 
     def stop(self):
--- ./ipi/engine/forces.py	(original)
+++ ./ipi/engine/forces.py	(refactored)
@@ -736,7 +736,7 @@
         if len(self.mforces) != len(refforce.mforces):
             raise ValueError("Cannot copy forces between objects with different numbers of components")
 
-        for k in xrange(len(self.mforces)):
+        for k in range(len(self.mforces)):
             mreff = refforce.mforces[k]
             mself = self.mforces[k]
             if mreff.nbeads != mself.nbeads:
@@ -753,7 +753,7 @@
             # tainted - it shRefactoringTool: Refactored ./ipi/engine/initializer.py
RefactoringTool: Refactored ./ipi/engine/normalmodes.py
RefactoringTool: Refactored ./ipi/engine/outputs.py
RefactoringTool: Refactored ./ipi/engine/properties.py
RefactoringTool: Refactored ./ipi/engine/simulation.py
ould not be as it's an internal of the force and
             # therefore get copied
             dd(mself.beads).q.set(mreff.beads.q, manual=False)
-            for b in xrange(mself.nbeads):
+            for b in range(mself.nbeads):
                 dfkbref = dd(mreff._forces[b])
                 dfkbself = dd(mself._forces[b])
 
--- ./ipi/engine/initializer.py	(original)
+++ ./ipi/engine/initializer.py	(refactored)
@@ -56,7 +56,7 @@
         self.mode = mode
         self.units = units
 
-        for (o, v) in others.items():
+        for (o, v) in list(others.items()):
             self.__dict__[o] = v
 
 
--- ./ipi/engine/normalmodes.py	(original)
+++ ./ipi/engine/normalmodes.py	(refactored)
@@ -145,7 +145,7 @@
         # if ( (len(self.bosons) > 0) and (len(self.bosons) < self.natoms) ):
         #raise(IOError("@NormalModes : Currently, only full bosonic/distinguishable simulations are allowed"))
         if(len(self.bosons) > self.natoms):
-            raise(IOError("@NormalModes : number of bosons is larger than number of atoms!"))
+            raise IOError
 
         dself = dd(self)
 
@@ -528,8 +528,8 @@
 
         # dynamical masses for the open paths
         for j in self.open_paths:
-            for a in xrange(3 * j, 3 * (j + 1)):
-                for k in xrange(1, self.nbeads):
+            for a in range(3 * j, 3 * (j + 1)):
+                for k in range(1, self.nbeads):
                     dm3[k, a] = self.beads.m3[k, a] * self.o_nm_factor[k]
         return dm3
 
@@ -655,8 +655,8 @@
             # and do open path propagation
             pq = np.zeros(2)
             for j in self.open_paths:
-                for a in xrange(3 * j, 3 * (j + 1)):
-                    for k in xrange(1, self.nbeads):
+                for a in range(3 * j, 3 * (j + 1)):
+                    for k in range(1, self.nbeads):
                         pq[0] = self.pnm[k, a] / sm[k, a]
                         pq[1] = self.qnm[k, a] * sm[k, a]
                         pq = np.dot(o_prop_pq[k], pq)
--- ./ipi/engine/outputs.py	(original)
+++ ./ipi/engine/outputs.py	(refactored)
@@ -163,8 +163,8 @@
         self.system = system
         for what in self.outlist:
             key = getkey(what)
-            if not key in system.properties.property_dict.keys():
-                print "Computable properties list: ", system.properties.property_dict.keys()
+            if not key in list(system.properties.property_dict.keys()):
+                print("Computable properties list: ", list(system.properties.property_dict.keys()))
                 raise KeyError(key + " is not a recognized property")
 
         super(PropertyOutput, self).bind(mode)
@@ -286,8 +286,8 @@
         self.system = system
         # Checks as soon as possible if some asked-for trajs are missing or mispelled
         key = getkey(self.what)
-        if not key in self.system.trajs.traj_dict.keys():
-            print "Computable trajectories list: ", self.system.trajs.traj_dict.keys()
+        if not key in list(self.system.trajs.traj_dict.keys()):
+            print("Computable trajectories list: ", list(self.system.trajs.traj_dict.keys()))
             raise KeyError(key + " is not a recognized output trajectory")
 
         super(TrajectoryOutput, self).bind(mode)
--- ./ipi/engine/properties.py	(original)
+++ ./ipi/engine/properties.py	(refactored)
@@ -835,7 +835,7 @@
         # subtracts centroid
         q = dstrip(self.beads.q).copy()
         qc = dstrip(self.beads.qc)
-        for b in xrange(self.beads.nbeads):
+        for b in range(self.beads.nbeads):
             q[b] -= qc
 
         # zeroes components that are not requested
@@ -933,7 +933,7 @@
                 if len(self.fqref) != 3 * self.beads.natoms:
                     raise ValueError("Atom number mismatch in reference file for virial_fq")
         fq = 0.0
-        for b in xrange(self.beads.nbeads):
+        for b in range(self.beads.nbeads):
             fq += np.dot(self.forces.f[b], self.beads.q[b] - self.fqref)
 
         return fq * 0.5 / self.beads.nbeads
--- ./ipi/engine/simRefactoringTool: Refactored ./ipi/engine/system.py
RefactoringTool: Refactored ./ipi/engine/thermostats.py
RefactoringTool: No changes to ./ipi/engine/motion/__init__.py
RefactoringTool: Refactored ./ipi/engine/motion/al6xxx_kmc.py
ulation.py	(original)
+++ ./ipi/engine/simulation.py	(refactored)
@@ -96,13 +96,13 @@
 
         # echo the input file if verbose enough
         if verbosity.level > 0:
-            print " # i-PI loaded input file: ", fn_input
+            print(" # i-PI loaded input file: ", fn_input)
         if verbosity.level > 1:
-            print " --- begin input file content ---"
+            print(" --- begin input file content ---")
             ifile = open(fn_input, "r")
             for line in ifile.readlines():
-                print line,
-            print " ---  end input file content  ---"
+                print(line, end=' ')
+            print(" ---  end input file content  ---")
             ifile.close()
 
         return simulation
@@ -178,7 +178,7 @@
 
         # start forcefields here so we avoid having a shitload of files printed
         # out only to find the socket is busy or whatever prevented starting the threads
-        for k, f in self.fflist.iteritems():
+        for k, f in self.fflist.items():
             f.start()
 
         # Checks for repeated filenames.
@@ -273,7 +273,7 @@
         #tttime = 0.0
         ttot = 0.0
         # main MD loop
-        for self.step in xrange(self.step, self.tsteps):
+        for self.step in range(self.step, self.tsteps):
             # stores the state before doing a step.
             # this is a bit time-consuming but makes sure that we can honor soft
             # exit requests without screwing the trajectory
--- ./ipi/engine/system.py	(original)
+++ ./ipi/engine/system.py	(refactored)
@@ -93,7 +93,7 @@
         self.simul = simul  # keeps a handle to the parent simulation object
 
         # binds important computation engines
-        print "NOW BINDING THE FORCES.... "
+        print("NOW BINDING THE FORCES.... ")
         self.forces.bind(self.beads, self.cell, self.fcomp, self.simul.fflist, open_paths=self.nm.open_paths)
         self.nm.bind(self.ensemble, self.motion, beads=self.beads, forces=self.forces)
         self.ensemble.bind(self.beads, self.nm, self.cell, self.forces, self.simul.fflist)
--- ./ipi/engine/thermostats.py	(original)
+++ ./ipi/engine/thermostats.py	(refactored)
@@ -994,10 +994,10 @@
 
             if self.intau != 0:
                 if mytemp != 0: self.intau /= (mytemp / self.temp)**(self.dt / self.apat)
-                print("ThermoCL inherent noise time scale: " + str(self.intau))
+                print(("ThermoCL inherent noise time scale: " + str(self.intau)))
             else:
                 self.idtau *= (mytemp / self.temp)**(self.dt / self.apat)
-                print("ThermoCL inherent dissipation time scale: " + str(self.idtau))
+                print(("ThermoCL inherent dissipation time scale: " + str(self.idtau)))
 
         self.idstep = not self.idstep
 
--- ./ipi/engine/motion/al6xxx_kmc.py	(original)
+++ ./ipi/engine/motion/al6xxx_kmc.py	(refactored)
@@ -114,23 +114,23 @@
         self.dcell = Cell()
         self.dcell.h = self.scell*self.ncell
 
-        print "LATTICE PARAM ", self.a0
+        print("LATTICE PARAM ", self.a0)
         # this is the list of lattice sites, in 3D coordinates
-        ix,iy,iz = np.meshgrid(range(self.ncell), range(self.ncell), range(self.ncell), indexing='ij')
+        ix,iy,iz = np.meshgrid(list(range(self.ncell)), list(range(self.ncell)), list(range(self.ncell)), indexing='ij')
         self.sites = np.dot(np.asarray([ix.flatten(),iy.flatten(),iz.flatten()]).T, self.scell.T)
-        print len(self.sites), self.nsites, "###"
+        print(len(self.sites), self.nsites, "###")
         # now we build list of nearest neighbors (fcc-lattice hardcoded!)
         self.neigh = np.zeros((self.nsites,12),int)
         nneigh = np.zeros(self.nsites,int)
         # could be done in a more analytic way but whatever, I'm too lazy
         a02 = 1.01*0.5*self.a0**2                                        # perhaps 1.01 it is not enough, must check!
-        for i in xrange(self.nsites): # determines the connectivity of the lattice
+        for i in range(self.nsites): # determines the connectivity of the lattice
             rij = self.sites.copy().flatten()
-            for j in xrange(self.nsites):
+            for j in range(self.nsites):
                 rij[3*j:3*j+3] -= self.sites[i]
             self.dcell.array_pbc(rij)
             rij.shape = (self.nsites,3)
-            for j in xrange(i):
+            for j in range(i):
                 if np.dot(rij[j],rij[j]) < a02: # found nearest neighbor
                     self.neigh[i,nneigh[i]] = j
                     self.neigh[j,nneigh[j]] = i
@@ -147,7 +147,7 @@
         # geop should not trigger exit if there is early convergence, but just carry on.
         # we hard-code this option to avoid early-termination that would be hard to debug for a user
         geop["exit_on_convergence"] = False
-        for i in xrange(self.neval):
+        for i in range(self.neval):
             # geometry optimizer should not have *any* hystory dependence
             self.geop[i] = GeopMotion(fixcom=fixcom, fixatoms=fixatoms,**geop) #mode="cg", ls_options={"tolerance": 1, "iter": 20,  "step": 1e-3, "adaptive": 0.0}, tolerances={"energy": 1e-7, "force": 1e-2, "position": 1e-4}, ) #!TODO: set the geop parameters properly
 
@@ -161,9 +161,9 @@
             ff = open(self.qcache_file, "rb")
             self.qcache = pickle.load(ff)
             ff.close()
-            print "Loaded %d cached energies" % (len(self.ecache))
+            print("Loaded %d cached energies" % (len(self.ecache)))
         except:
-            print "Couldn't load cache files "+self.ecache_file+","+self.qcache_file+" - resetting"
+            print("Couldn't load cache files "+self.ecache_file+","+self.qcache_file+" - resetting")
             self.ecache = {}
             self.qcache = {}
         self.ncache = len(self.ecache)
@@ -206,7 +206,7 @@
         f_restart = True
         if self.idx is None or len(self.idx)==0:
             f_restart = False
-            idx =np.asarray(range(self.ncell**3), int) # initialize random distribution of atoms
+            idx =np.asarray(list(range(self.ncell**3)), int) # initialize random distribution of atoms
             self.prng.rng.shuffle(idx)
             self.idx = idx
 
@@ -224,20 +224,20 @@
 
         # reverse lookup index [i.e. ridx[i] gives the index of the atoms at site i]
         self.ridx = np.zeros(self.nsites,int)
-        self.ridx[self.idx] = range(self.nsites)
+        self.ridx[self.idx] = list(range(self.nsites))
 
         self.state = np.asarray(list(state))
-        print  "".join(self.state)
-
-
-        print "CHECKING INITIAL ASSIGNMENTS"
-        for i in xrange(self.nsites):
+        print("".join(self.state))
+
+
+        print("CHECKING INITIAL ASSIGNMENTS")
+        for i in range(self.nsites):
             if self.ridx[i]<self.natoms:
-                print self.beads.names[self.ridx[i]], self.state[i]
+                print(self.beads.names[self.ridx[i]], self.state[i])
             else:
-                print "V", self.state[i]
+                print("V", self.state[i])
             if self.idx[self.ridx[i]] != i:
-                print "inconsistent site string for atom ",i, " and site ", self.ridx[i]
+                print("inconsistent site string for atom ",i, " and site ", self.ridx[i])
 
         if not f_restart:
             self.beads.q[0,:] = self.sites[self.idx].flatten() # also sets global q so we can use it further down
@@ -247,7 +247,7 @@
         self.dnm = [None] * self.neval
         self.dens = [None] * self.neval
         self.dbias = [None] * self.neval
-        for i in xrange(self.neval):
+        for i in range(self.neval):
             self.dbeads[i] = beads.copy()
             self.dforces[i] = bforce.copy(self.dbeads[i], self.dcell)
             self.dnm[i] = nm.copy()
@@ -264,7 +264,7 @@
         self.geop[ieval].reset()
         ipot = self.dforces[ieval].pot
 
-        for i in xrange(self.nstep):
+        for i in range(self.nstep):
             # print "geop ", i, self.dforces[ieval].pot
             self.geop[ieval].step(i)
             #if self.geop[ieval].converged[0]: break
@@ -289,8 +289,8 @@
             self.feval[ieval] = 1
 
         with self._threadlock:
-            print "Finished ", nstr
-            print "Energy, initial - TS - final: ", ipot, nevent[-1], newpot
+            print("Finished ", nstr)
+            print("Energy, initial - TS - final: ", ipot, nevent[-1], newpot)
 
     # threaded ts evaluation
     def ts_thread(self, ieval, ostr, nstr, nevent, setfev=1):
@@ -334,7 +334,7 @@
                     break
         with self._threadlock:
             # finds free evaluator
-            for e in xrange(self.neval):
+            for e in range(self.neval):
                 if self.feval[e] == 1:
                     ieval = e
                     self.feval[ieval] = 0
@@ -379,15 +379,15 @@
         uid_1 = self.unique_idx(state_1)
         ru1 = np.zeros(self.nsites,int)
         # this is the reverse map. what is the atom index that sits in a given site?
-        ru1[uid_1] = np.asarray(range(self.nsites),int)
+        ru1[uid_1] = np.asarray(list(range(self.nsites)),int)
 
         uid_2 = self.unique_idx(state_2)
         ru2 = np.zeros(self.nsites,int)
-        ru2[uid_2] = np.asarray(range(self.nsites),int)  # this says which atom is in a given site in u2
+        ru2[uid_2] = np.asarray(list(range(self.nsites)),int)  # this says which atom is in a given site in u2
 
         iu12 = ru2[uid_1]
         iu21 = ru1[uid_2]
-        for i in xrange(self.natoms):
+        for i in range(self.natoms):
             if iu12[i] >= self.natoms:
                 #print "found vacancy swap 1->2", i, u1[i], u2[iu12[i]], iu12[i]
                 i1vac2 = i
@@ -418,7 +418,7 @@
         levents = []
         ethreads = [None] * self.neval
         # loops over the vacancy
-        for ivac in xrange(self.natoms, self.natoms + self.nvac):
+        for ivac in range(self.natoms, self.natoms + self.nvac):
             svac = self.idx[ivac] # lattice site associated with this vacancy
             if self.state[svac] != "V":
                 raise IndexError("Something got screwed and a vacancy state is not a vacancy anymore!")
@@ -455,7 +455,7 @@
                     st.start()
                     ethreads[ieval] = st
                 else:
-                    print "Found state ", nstr, " retrieving cached energy ", self.ecache[nstr]
+                    print("Found state ", nstr, " retrieving cached energy ", self.ecache[nstr])
 
                     # fetch energy from previous calculation
                     nevent = [svac, sneigh, self.ecache[nstr], self.qcache[nstr], 0.0]
@@ -479,17 +479,17 @@
             while not st is None and st.isAlive():
                 st.join(2)
 
-        print "Computed ", len(levents), " possible reactions. Cache len ", len(self.ecache)
+        print("Computed ", len(levents), " possible reactions. Cache len ", len(self.ecache))
 
         # get list of rates
         rates = np.zeros(len(levents), float)
         crates = np.zeros(len(levents), float)
         cdf = 0.0
-        for i in xrange(len(levents)):
+        for i in range(len(levents)):
             #print ("Barrier, naive: %f, static: %f" % (0.5*(ecurr + levents[i][2]) + self.diffusion_barrier_al, levents[i][4]))
 
             ets = 0.5*(ecurr + levents[i][2]) + self.barriers[levents[i][-1]]  #naive heuristic for the ts energy
-            print "Event ", i, levents[i][-1], ecurr, ">>", ets, ">>", levents[i][2]
+            print("Event ", i, levents[i][-1], ecurr, ">>", ets, ">>", levents[i][2])
             rates[i] = self.prefactors[levents[i][-1]] * np.exp(-(ets-ecurr)/kT)
             cdf += rates[i]
             crates[i] = cdf
@@ -500,8 +500,8 @@
         while fpick > crates[isel]:
             isel += 1
         dt = -1.0/cdf*np.log(1.0-self.prng.u)
-        print ("Time spent %12.5e at %s nrg %12.5e" % (dt, ostr,ecurr))
-        print "Selected event ", isel, " with rate ", rates[isel], " / ", cdf
+        print(("Time spent %12.5e at %s nrg %12.5e" % (dt, ostr,ecurr)))
+        print("Selected event ", isel, " with rate ", rates[isel], " / ", cdf)
 
         iev = levRefactoringTool: Refactored ./ipi/engine/motion/alchemy.py
RefactoringTool: Refactored ./ipi/engine/motion/atomswap.py
RefactoringTool: Refactored ./ipi/engine/motion/constrained_dynamics.py
RefactoringTool: Refactored ./ipi/engine/motion/dynamics.py
RefactoringTool: Refactored ./ipi/engine/motion/geop.py
ents[isel] # levents[self.prng.randint(len(levents))]
         svac, sneigh = iev[0], iev[1]
@@ -518,21 +518,21 @@
         self.kmcfile.flush()
         self.tottime += dt
         self.ensemble.time += dt  # updates time counter
-        print  "Finishing step at ", "".join(self.state)
+        print("Finishing step at ", "".join(self.state))
 
         # updates the positions
         self.cell.h = self.dcell.h
 
         uidx = self.unique_idx(self.state)
         ruidx = np.zeros(self.nsites,int)
-        ruidx[uidx] = range(self.nsites)
+        ruidx[uidx] = list(range(self.nsites))
 
         self.sites[self.unique_idx(self.state)]
         oldq = dstrip(self.beads.q[0]).copy()
 
         newq = np.zeros(self.nsites*3, float)
         # we want continuity (modulo PBC jumps, that we'll take care of later...)
-        for i in xrange(self.nsites):
+        for i in range(self.nsites):
             # in which site sits atom i?
             isite = self.idx[i]
             # which atom sits in this site in the unique-mapped structure?
--- ./ipi/engine/motion/alchemy.py	(original)
+++ ./ipi/engine/motion/alchemy.py	(refactored)
@@ -120,7 +120,7 @@
         # this would be double-counting, we already have a bail-out condition above
         # if (1.0/self.nxc < self.prng.u) : return  # tries a round of exhanges with probability 1/nmc
 
-        for x in xrange(ntries):
+        for x in range(ntries):
             i = self.prng.rng.randint(lenlist)
             j = self.prng.rng.randint(lenlist)
             while self.beads.names[axlist[i]] == self.beads.names[axlist[j]]:
--- ./ipi/engine/motion/atomswap.py	(original)
+++ ./ipi/engine/motion/atomswap.py	(refactored)
@@ -103,7 +103,7 @@
         # this would be double-counting, we already have a bail-out condition above
         # if (1.0/self.nxc < self.prng.u) : return  # tries a round of exhanges with probability 1/nmc
         self.dcell.h = self.cell.h  # just in case the cell gets updated in the other motion classes
-        for x in xrange(ntries):
+        for x in range(ntries):
             i = self.prng.rng.randint(lenlist)
             j = self.prng.rng.randint(lenlist)
             while self.beads.names[axlist[i]] == self.beads.names[axlist[j]]:
--- ./ipi/engine/motion/constrained_dynamics.py	(original)
+++ ./ipi/engine/motion/constrained_dynamics.py	(refactored)
@@ -281,7 +281,7 @@
             constr.q = q[ic]
 
             # iterative projection on the manifold
-            for i in xrange(self.maxit):
+            for i in range(self.maxit):
                 g = dstrip(constr.g)
                 # bailout condition
                 if self.tolerance > np.linalg.norm(g, ord=self.norm_order):
@@ -340,13 +340,13 @@
             dself.nsteps_geo = depend_value(name="nsteps_geo", value=1)
         else:
             dself.nsteps_geo = depend_value(name="nsteps_geo", value=motion.nsteps_geo)
-        print self.nsteps_geo
+        print(self.nsteps_geo)
         dself.qdt.add_dependency(dself.nsteps_geo)
         if motion.nsteps_o is None:
             dself.nsteps_o = depend_value(name="nsteps_o", value=1)
         else:
             dself.nsteps_o = depend_value(name="nsteps_o", value=motion.nsteps_o)
-        print self.nsteps_o
+        print(self.nsteps_o)
         dself.tdt.add_dependency(dself.nsteps_o)
         dd(self).csolver = motion.csolver
         dpipe(dself.qdt, dd(self.csolver).dt)
@@ -485,7 +485,7 @@
         (e.g. in GLE) by doing a MTS splitting scheme """
 
         m3 = dstrip(self.beads.m3)
-        for i in xrange(self.nsteps_o):
+        for i in range(self.nsteps_o):
             self.thermostat.step()
 
             # accumulates conserved quantity
--- ./ipi/engine/motion/dynamics.py	(original)
+++ ./ipi/engine/motion/dynamics.py	(refactored)
@@ -218,7 +218,7 @@
 
     def get_pdt(self):
         dtl = 1.0 / self.nmts
-        for i in xrange(1, len(dtl)):
+        for i in range(1, len(dtl)):
             dtl[i] *= dtl[i - 1]
         dtl *= self.dt * 0.5
         return dtl
--- ./ipi/engine/motion/geop.py	(original)
+++ ./ipi/RefactoringTool: No changes to ./ipi/engine/motion/instanton.py
RefactoringTool: No changes to ./ipi/engine/motion/motion.py
RefactoringTool: No changes to ./ipi/engine/motion/multi.py
RefactoringTool: Refactored ./ipi/engine/motion/neb.py
RefactoringTool: Refactored ./ipi/engine/motion/phonons.py
RefactoringTool: Refactored ./ipi/engine/motion/planetary.py
RefactoringTool: No changes to ./ipi/engine/motion/ramp.py
RefactoringTool: No changes to ./ipi/engine/motion/replay.py
RefactoringTool: No changes to ./ipi/engine/smotion/__init__.py
RefactoringTool: No changes to ./ipi/engine/smotion/metad.py
RefactoringTool: No changes to ./ipi/engine/smotion/multi.py
RefactoringTool: Refactored ./ipi/engine/smotion/remd.py
RefactoringTool: No changes to ./ipi/engine/smotion/smotion.py
RefactoringTool: Refactored ./ipi/external/importlib/__init__.py
RefactoringTool: Refactored ./ipi/external/importlib/bundledimportlib.py
RefactoringTool: No changes to ./ipi/inputs/__init__.py
RefactoringTool: No changes to ./ipi/inputs/atoms.py
RefactoringTool: No changes to ./ipi/inputs/barostats.py
RefactoringTool: No changes to ./ipi/inputs/beads.py
RefactoringTool: No changes to ./ipi/inputs/cell.py
RefactoringTool: No changes to ./ipi/inputs/ensembles.py
RefactoringTool: No changes to ./ipi/inputs/forcefields.py
RefactoringTool: No changes to ./ipi/inputs/forces.py
RefactoringTool: No changes to ./ipi/inputs/initializer.py
RefactoringTool: No changes to ./ipi/inputs/interface.py
RefactoringTool: No changes to ./ipi/inputs/normalmodes.py
RefactoringTool: No changes to ./ipi/inputs/outputs.py
RefactoringTool: No changes to ./ipi/inputs/prng.py
RefactoringTool: Refactored ./ipi/inputs/simulation.py
RefactoringTool: Refactored ./ipi/inputs/system.py
engine/motion/geop.py	(refactored)
@@ -553,7 +553,7 @@
 
         if step == 0:
             info(" @GEOP: Initializing L-BFGS", verbosity.debug)
-            print self.d
+            print(self.d)
             self.d += dstrip(self.forces.f) / np.sqrt(np.dot(self.forces.f.flatten(), self.forces.f.flatten()))
 
         self.old_x[:] = self.beads.q
--- ./ipi/engine/motion/neb.py	(original)
+++ ./ipi/engine/motion/neb.py	(refactored)
@@ -252,7 +252,7 @@
                     btau[ii] = d2 * maxpot + d1 * minpot
 
                 else:
-                    print "Error in NEB tangents: Energy of images are equal"
+                    print("Error in NEB tangents: Energy of images are equal")
 
             btau[ii] *= 1.0 / np.linalg.norm(btau)
 
--- ./ipi/engine/motion/phonons.py	(original)
+++ ./ipi/engine/motion/phonons.py	(refactored)
@@ -332,7 +332,7 @@
 
         self.dm.w2, self.dm.U = np.linalg.eigh(self.dm.dynmatrix)
         self.dm.V = self.dm.U.copy()
-        for i in xrange(len(self.dm.V)):
+        for i in range(len(self.dm.V)):
             self.dm.V[:, i] *= self.dm.ism
 
     def step(self, step=None):
--- ./ipi/engine/motion/planetary.py	(original)
+++ ./ipi/engine/motion/planetary.py	(refactored)
@@ -209,7 +209,7 @@
         self.tmtx += time.time()
 
         # sample by constrained-centroid dynamics
-        for istep in xrange(self.nsamples):
+        for istep in range(self.nsamples):
             self.tmc -= time.time()
             self.ccdyn.step(step)
             self.tmc += time.time()
--- ./ipi/engine/smotion/remd.py	(original)
+++ ./ipi/engine/smotion/remd.py	(refactored)
@@ -84,7 +84,7 @@
         super(ReplicaExchange, self).bind(syslist, prng, omaker)
 
         if self.repindex is None or len(self.repindex) == 0:
-            self.repindex = np.asarray(range(len(self.syslist)))
+            self.repindex = np.asarray(list(range(len(self.syslist))))
         else:
             if len(self.syslist) != len(self.repindex):
                 raise ValueError("Size of replica index does not match number of systems replicas")
--- ./ipi/external/importlib/__init__.py	(original)
+++ ./ipi/external/importlib/__init__.py	(refactored)
@@ -2,4 +2,4 @@
     # importlib is a python stdlib from python 2.7 on
     from importlib import *
 except ImportError:
-    from bundledimportlib import *
+    from .bundledimportlib import *
--- ./ipi/external/importlib/bundledimportlib.py	(original)
+++ ./ipi/external/importlib/bundledimportlib.py	(refactored)
@@ -9,7 +9,7 @@
     if not hasattr(package, 'rindex'):
         raise ValueError("'package' not set to a string")
     dot = len(package)
-    for x in xrange(level, 1, -1):
+    for x in range(level, 1, -1):
         try:
             dot = package.rindex('.', 0, dot)
         except ValueError:
--- ./ipi/inputs/simulation.py	(original)
+++ ./ipi/inputs/simulation.py	(refactored)
@@ -141,7 +141,7 @@
 
         self.mode.store(simul.mode)
 
-        _fflist = [v for k, v in sorted(simul.fflist.iteritems())]
+        _fflist = [v for k, v in sorted(simul.fflist.items())]
         if len(self.extra) != len(_fflist) + len(simul.syslist):
             self.extra = [0] * (len(_fflist) + len(simul.syslist))
 
@@ -209,7 +209,7 @@
             elif k == "system_template":
                 syslist += v.fetch()  # this will actually generate automatically a bunch of system objects with the desired properties set automatically to many values
             elif k == "ffsocket" or k == "fflj" or k == 'ffquip' or k == "ffdebye" or k == "ffplumed" or k == "ffsgdml":
-                print "fetching", k
+                print("fetching", k)
                 fflist.append(v.fetch())
             elif k == "ffyaff":
                 fflist.append(v.fetch())
--- ./ipi/inputs/system.py	(original)
+++ ./ipi/inputs/system.py	(refactored)
@@ -84,9 +84,9 @@
                 sys = template
                 if len(labels) != len(ins):
                     raise ValueError("Labels and instance length mismatch")
-                for l in xrange(len(ins)):  # string replacement within tRefactoringTool: No changes to ./ipi/inputs/thermostats.py
RefactoringTool: No changes to ./ipi/inputs/motion/__init__.py
RefactoringTool: Refactored ./ipi/inputs/motion/al6xxx_kmc.py
RefactoringTool: No changes to ./ipi/inputs/motion/alchemy.py
RefactoringTool: No changes to ./ipi/inputs/motion/atomswap.py
RefactoringTool: No changes to ./ipi/inputs/motion/constrained_dynamics.py
RefactoringTool: No changes to ./ipi/inputs/motion/dynamics.py
RefactoringTool: No changes to ./ipi/inputs/motion/geop.py
RefactoringTool: No changes to ./ipi/inputs/motion/instanton.py
RefactoringTool: No changes to ./ipi/inputs/motion/motion.py
RefactoringTool: No changes to ./ipi/inputs/motion/neb.py
RefactoringTool: No changes to ./ipi/inputs/motion/phonons.py
RefactoringTool: No changes to ./ipi/inputs/motion/planetary.py
RefactoringTool: No changes to ./ipi/inputs/motion/ramp.py
RefactoringTool: No changes to ./ipi/inputs/smotion/__init__.py
RefactoringTool: No changes to ./ipi/inputs/smotion/metad.py
RefactoringTool: No changes to ./ipi/inputs/smotion/remd.py
RefactoringTool: No changes to ./ipi/inputs/smotion/smotion.py
RefactoringTool: No changes to ./ipi/interfaces/__init__.py
RefactoringTool: Refactored ./ipi/interfaces/clients.py
he template
+                for l in range(len(ins)):  # string replacement within the template
                     sys = sys.replace(labels[l], ins[l])
-                print "Generating system from template: \n", sys
+                print("Generating system from template: \n", sys)
                 xsys = xml_parse_string(sys)  # parses the string to an XML object
                 isys = InputSystem()
                 isys.parse(xsys.fields[0][1])  # parses the XML object into an InputSystem
--- ./ipi/inputs/motion/al6xxx_kmc.py	(original)
+++ ./ipi/inputs/motion/al6xxx_kmc.py	(refactored)
@@ -146,12 +146,12 @@
         # only stores cache after a decent amount of new structures have been found
         if kmc.ncache_stored*self.STORE_STRIDE<kmc.ncache:
             if kmc.ecache_file != "":
-                print "Storing ECACHE in ", kmc.ecache_file
+                print("Storing ECACHE in ", kmc.ecache_file)
                 ff = open(kmc.ecache_file, "wb")
                 pickle.dump(kmc.ecache, ff)
                 ff.close()
             if kmc.qcache_file != "":
-                print "Storing QCACHE in ", kmc.qcache_file
+                print("Storing QCACHE in ", kmc.qcache_file)
                 ff = open(kmc.qcache_file, "wb")
                 pickle.dump(kmc.qcache, ff)
                 ff.close()
--- ./ipi/interfaces/clients.py	(original)
+++ ./ipi/interfaces/clients.py	(refactored)
@@ -50,7 +50,7 @@
                     _socket = socket.socket(socket.AF_UNIX, socket.SOCK_STREAM)
                     _socket.connect("/tmp/ipi_" + address)
                 except socket.error:
-                    print 'Could not connect to UNIX socket: %s' % ("/tmp/ipi_" + address)
+                    print('Could not connect to UNIX socket: %s' % ("/tmp/ipi_" + address))
                     sys.exit(1)
             else:
                 raise NameError("Interface mode " + mode + " is not implemented (should be unix/inet)")
@@ -97,17 +97,17 @@
         fmt_step = '{0:6d} {1:10.3f} {2:10.3f}'
 
         if t_max is None:
-            print 'Starting communication loop with no maximum run time.'
+            print('Starting communication loop with no maximum run time.')
         else:
-            print 'Starting communication loop with a maximum run time of {0:d} seconds.'.format(t_max)
+            print('Starting communication loop with a maximum run time of {0:d} seconds.'.format(t_max))
             fmt_header += ' {3:>10s}'
             fmt_step += ' {3:10.1f}'
 
         if verbose:
             header = fmt_header.format('step', 'time', 'avg time', 'remaining')
-            print
-            print header
-            print len(header) * '-'
+            print()
+            print(header)
+            print(len(header) * '-')
 
         i_step = 0
         t_step_tot = 0.0
@@ -121,7 +121,7 @@
 
                 # process message and respond
                 if msg == "":
-                    print "Server shut down."
+                    print("Server shut down.")
                     break
                 elif msg == Message("status"):
                     if self.havedata:
@@ -142,7 +142,7 @@
                         t_step_avg = t_step_tot / (i_step + 1)
                         if t_max is not None:
                             t_remain = t_max - (t_now - t0)
-                        print fmt_step.format(i_step, t_step, t_step_avg, t_remain)
+                        print(fmt_step.format(i_step, t_step, t_step_avg, t_remain))
                     self.havedata = True
                     i_step += 1
                 elif msg == Message("getforce"):
@@ -154,25 +154,25 @@
                     self.sendall(np.int32(0), 4)
                     self.havedata = False
                 else:
-                    print >> sys.stderr, "Client could not understand command:", msg
+                    print("Client could not understand command:", msg, file=sys.stderr)
                     break
 
                 # check exit conditions - run time or exit file
                 if t_max is not None and time.time() RefactoringTool: Refactored ./ipi/interfaces/sockets.py
RefactoringTool: No changes to ./ipi/tests/__init__.py
RefactoringTool: Refactored ./ipi/tests/common.py
RefactoringTool: Refactored ./ipi/tests/test_atoms.py
RefactoringTool: Refactored ./ipi/tests/test_contraction.py
RefactoringTool: Refactored ./ipi/tests/test_depend.py
RefactoringTool: No changes to ./ipi/tests/test_interface.py
RefactoringTool: Refactored ./ipi/tests/test_io.py
- t0 > t_max:
-                    print 'Maximum run time of {0:d} seconds exceeded.'.format(t_max)
+                    print('Maximum run time of {0:d} seconds exceeded.'.format(t_max))
                     break
                 if fn_exit is not None and os.path.exists(fn_exit):
-                    print 'Exit file "{0:s}" found. Removing file.'.format(fn_exit)
+                    print('Exit file "{0:s}" found. Removing file.'.format(fn_exit))
                     os.remove(fn_exit)
                     break
 
         except socket.error as e:
-            print 'Error communicating through socket: [{0}] {1}'.format(e.errno, e.strerror)
+            print('Error communicating through socket: [{0}] {1}'.format(e.errno, e.strerror))
         except KeyboardInterrupt:
-            print ' Keyboard interrupt.'
-
-        print 'Communication loop finished.'
-        print
+            print(' Keyboard interrupt.')
+
+        print('Communication loop finished.')
+        print()
 
 
 class ClientASE(Client):
--- ./ipi/interfaces/sockets.py	(original)
+++ ./ipi/interfaces/sockets.py	(refactored)
@@ -334,7 +334,7 @@
                 self.sendall(pos)
                 self.status = Status.Up | Status.Busy
             except:
-                print "Error in sendall, resetting status"
+                print("Error in sendall, resetting status")
                 self.get_status()
                 return
         else:
--- ./ipi/tests/common.py	(original)
+++ ./ipi/tests/common.py	(refactored)
@@ -73,8 +73,8 @@
 
         # Start simulation
         # TODO
-        print subprocess.check_output("ls", shell=True)
-        print subprocess.check_output("pwd", shell=True)
+        print(subprocess.check_output("ls", shell=True))
+        print(subprocess.check_output("pwd", shell=True))
 
         # wait for driver to finish
         p.wait()
--- ./ipi/tests/test_atoms.py	(original)
+++ ./ipi/tests/test_atoms.py	(refactored)
@@ -5,7 +5,7 @@
 # See the "licenses" directory for full license information.
 
 
-from common import local
+from .common import local
 
 from ipi.utils.io import read_file
 
--- ./ipi/tests/test_contraction.py	(original)
+++ ./ipi/tests/test_contraction.py	(refactored)
@@ -21,17 +21,17 @@
     """
 
     rescale = nmtransform.nm_rescale(q.shape[0], n)
-    print "Initial position of the beads:"
-    print q, q.shape, (q.shape[0], n)
+    print("Initial position of the beads:")
+    print(q, q.shape, (q.shape[0], n))
 
     # rescale up to the n beads
     beads_n = rescale.b1tob2(q)
-    print "Upscaled to %d beads:" % n
-    print beads_n, beads_n.shape
+    print("Upscaled to %d beads:" % n)
+    print(beads_n, beads_n.shape)
 
     beads_final = rescale.b2tob1(beads_n)
-    print "Final position of the beads:"
-    print beads_final
+    print("Final position of the beads:")
+    print(beads_final)
 
     assert_almost_equal(q, beads_final)
     return beads_n
@@ -74,7 +74,7 @@
     assert_almost_equal(centroid_q, centroid_big)
 
 
-numbers_to_check = range(10, 56, 9)
+numbers_to_check = list(range(10, 56, 9))
 
 
 qs_to_check = [
--- ./ipi/tests/test_depend.py	(original)
+++ ./ipi/tests/test_depend.py	(refactored)
@@ -18,14 +18,14 @@
 def test_slicing():
     """Depend: Slicing test"""
     c = a[0]
-    print(type(c))
+    print((type(c)))
     assert isinstance(c, dp.depend_array)
 
 
 def test_addition():
     """Depend: Addition test"""
     c = a + b
-    print(type(c))
+    print((type(c)))
     assert isinstance(c, np.ndarray)
 
 
@@ -33,14 +33,14 @@
     """Depend: Increment test"""
     c = np.zeros((2, 2))
     c += a
-    print(type(c))
+    print((type(c)))
     assert isinstance(c, np.ndarray)
 
 
 def test_dot():
     """Depend: Dot test"""
     c = np.dot(a, b)
-    print(type(c))
+    print((type(c)))
     assert isinstance(c, np.ndarray)
 
 
--- ./ipi/tests/test_io.py	(original)
+++ ./ipi/tests/test_io.py	(refactored)
@@ -15,7 +15,7 @@
 import numpy as np
 from numpy.testing import assert_equal
 
-from common import local
+from .common import local
 from ipi.utils.units import uniRefactoringTool: Refactored ./ipi/tests/test_runs.py
RefactoringTool: No changes to ./ipi/tests/test_units.py
RefactoringTool: No changes to ./ipi/utils/__init__.py
RefactoringTool: No changes to ./ipi/utils/constrtools.py
RefactoringTool: No changes to ./ipi/utils/decorators.py
RefactoringTool: Refactored ./ipi/utils/depend.py
RefactoringTool: No changes to ./ipi/utils/exchange.py
RefactoringTool: Refactored ./ipi/utils/hesstools.py
RefactoringTool: Refactored ./ipi/utils/inputvalue.py
t_to_internal
 from ipi.utils.io import iter_file, read_file, print_file
 
--- ./ipi/tests/test_runs.py	(original)
+++ ./ipi/tests/test_runs.py	(refactored)
@@ -5,7 +5,7 @@
 # See the "licenses" directory for full license information.
 
 
-from common import SimulationTest
+from .common import SimulationTest
 
 
 def test_lj_gas():
--- ./ipi/utils/depend.py	(original)
+++ ./ipi/utils/depend.py	(refactored)
@@ -231,7 +231,7 @@
             if (not item()._tainted[0]):
                 item().taint()
         if self._synchro is not None:
-            for v in self._synchro.synced.values():
+            for v in list(self._synchro.synced.values()):
                 if (not v._tainted[0]) and (v is not self):
                     v.taint(taintme=True)
             self._tainted[:] = (taintme and (not self._name == self._synchro.manual))
--- ./ipi/utils/hesstools.py	(original)
+++ ./ipi/utils/hesstools.py	(refactored)
@@ -16,7 +16,7 @@
     # Check dimensions
 
     if h.shape != (m3.shape[1], m3.shape[1] * nbeads):
-        print h.shape, m3.shape
+        print(h.shape, m3.shape)
         raise ValueError("@get_dynmat: The provided hessian hasn't the proper dimension (3*natoms, 3*natoms*nbeads) ")
 
     ism = m3.reshape((1, -1)) ** (-0.5)
@@ -141,8 +141,8 @@
             info(" Warning @Clean hessian: We have deleted %d 'zero' frequencies " % (nzero.size), verbosity.high)
             info(" but the norm is greater than 0.01 cm^-1.  This should not happen.", verbosity.high)
 
-        d = np.delete(d, range(nneg.size, nneg.size + nzero.size))
-        w = np.delete(w, range(nneg.size, nneg.size + nzero.size), axis=1)
+        d = np.delete(d, list(range(nneg.size, nneg.size + nzero.size)))
+        w = np.delete(w, list(range(nneg.size, nneg.size + nzero.size)), axis=1)
 
     if mofi:
         if asr == 'poly':
@@ -187,7 +187,7 @@
         else:
             h[:, :] = b[:, :]
             i0 = i
-            print('We have found a temporary file ( hessian_' + str(i) + '.tmp). ')
+            print(('We have found a temporary file ( hessian_' + str(i) + '.tmp). '))
             if b.shape == h.shape:  # Check that the last temporary file was properly written
                 break
             else:
--- ./ipi/utils/inputvalue.py	(original)
+++ ./ipi/utils/inputvalue.py	(refactored)
@@ -172,10 +172,10 @@
         # merge instencefields with the static class fields
         self.instancefields.update(self.fields)
 
-        for f, v in self.instancefields.iteritems():
+        for f, v in self.instancefields.items():
             self.__dict__[f] = v[0](**v[1])
 
-        for a, v in self.attribs.iteritems():
+        for a, v in self.attribs.items():
             self.__dict__[a] = v[0](**v[1])
 
         self.set_default()
@@ -235,7 +235,7 @@
             newfield = self.dynamic[name][0](**self.dynamic[name][1])
             newfield.parse(xml)
         except:
-            print "Error parsing " + name + " from " + str(xml)
+            print("Error parsing " + name + " from " + str(xml))
             raise
         self.extra.append((name, newfield))
 
@@ -324,7 +324,7 @@
         if xml is None:
             self._text = text
         else:
-            for a, v in xml.attribs.iteritems():
+            for a, v in xml.attribs.items():
                 if a in self.attribs:
                     self.__dict__[a].parse(text=v)
                 elif a == "_text":
@@ -486,7 +486,7 @@
                 rstr += self.__dict__[f].help_latex(name=f, level=level + 1, stop_level=stop_level, standalone=standalone)
 
         if len(self.dynamic) != 0 and level != stop_level:
-            for f, v in self.dynamic.iteritems():
+            for f, v in self.dynamic.items():
                 dummy_obj = v[0](**v[1])
                 rstr += dummy_obj.help_latex(name=f, level=level + 1, stop_level=stop_level, standalone=standalone)
 
@@ -647,7 +647,7 @@
         if show_fields:
             for f in self.instancefields:
                 rstr += self.__dict__[f].help_xml(f, "   " + indent, level + 1, stop_level)
-            for f, v in sRefactoringTool: Refactored ./ipi/utils/instools.py
elf.dynamic.iteritems():
+            for f, v in self.dynamic.items():
                 # we must create the object manually, as dynamic objects are
                 # not automatically added to the input object's dictionary
                 dummy_obj = v[0](**v[1])
@@ -697,7 +697,7 @@
         """Base function for storing data passed as a dictionary"""
 
         self._explicit = True
-        for f, v in value.iteritems():
+        for f, v in value.items():
             self.__dict__[f].store(value[f])
 
         pass
@@ -707,7 +707,7 @@
 
         self.check()
         rdic = {}
-        for f, v in self.instancefields.iteritems():
+        for f, v in self.instancefields.items():
             rdic[f] = self.__dict__[f].fetch()
         return rdic
 
@@ -873,7 +873,7 @@
         super(InputValue, self).fetch()
 
         if self._dimension != "undefined":
-            print "returning ", self.value * unit_to_internal(self._dimension, self.units.fetch(), 1.0)
+            print("returning ", self.value * unit_to_internal(self._dimension, self.units.fetch(), 1.0))
             return self.value * unit_to_internal(self._dimension, self.units.fetch(), 1.0)
         else:
             return self.value
--- ./ipi/utils/instools.py	(original)
+++ ./ipi/utils/instools.py	(refactored)
@@ -93,7 +93,7 @@
     M = A.shape[1]
     newA = np.empty((u + l + 1, M))
     newA[:u + 1] = A
-    for i in xrange(1, l + 1):
+    for i in range(1, l + 1):
         newA[u + i, :M - i] = A[-1 - i, i:]
     return newA
 
@@ -194,9 +194,9 @@
 def print_instanton_geo(prefix, step, nbeads, natoms, names, q, f, pots, cell, shift, output_maker):
 
     outfile = output_maker.get_output(prefix + '_' + str(step) + '.ener', 'w')
-    print >> outfile, ('#Bead    Energy (eV)')
-    for i in range(nbeads):
-        print >> outfile, (str(i) + '     ' + str(units.unit_to_user('energy', "electronvolt", pots[i] - shift)))
+    print(('#Bead    Energy (eV)'), file=outfile)
+    for i in range(nbeads):
+        print((str(i) + '     ' + str(units.unit_to_user('energy', "electronvolt", pots[i] - shift))), file=outfile)
     outfile.close()
 
     # print_file("xyz", pos[0], cell, out, title='positions{angstrom}')
@@ -208,24 +208,24 @@
     outfile = output_maker.get_output(prefix + '_' + str(step) + '.xyz', 'w')
     outfile2 = output_maker.get_output(prefix + '_forces_' + str(step) + '.xyz', 'w')
     for i in range(nbeads):
-        print >> outfile, natoms
-        print >> outfile2, natoms
-
-        print >> outfile, ('CELL(abcABC):  %f %f %f %f %f %f cell{atomic_unit}  Traj: positions{%s}   Bead:       %i' % (a, b, c, alpha, beta, gamma, unit, i))
-        print >> outfile2, ('CELL(abcABC):  %f %f %f %f %f %f cell{atomic_unit}  Traj: positions{%s}   Bead:       %i' % (a, b, c, alpha, beta, gamma, unit2, i))
+        print(natoms, file=outfile)
+        print(natoms, file=outfile2)
+
+        print(('CELL(abcABC):  %f %f %f %f %f %f cell{atomic_unit}  Traj: positions{%s}   Bead:       %i' % (a, b, c, alpha, beta, gamma, unit, i)), file=outfile)
+        print(('CELL(abcABC):  %f %f %f %f %f %f cell{atomic_unit}  Traj: positions{%s}   Bead:       %i' % (a, b, c, alpha, beta, gamma, unit2, i)), file=outfile2)
         # print >> outfile, ('#Potential (eV):   ' + str(units.unit_to_user('energy', "electronvolt", pots[i] - shift)))
 
         for j in range(natoms):
-            print >> outfile, names[j], \
+            print(names[j], \
                 str(units.unit_to_user('length', unit, q[i, 3 * j])), \
                 str(units.unit_to_user('length', unit, q[i, 3 * j + 1])), \
-                str(units.unit_to_user('length', unit, q[i, 3 * j + 2]))
+                str(units.unit_to_user('length', unit, q[i, 3 * j + 2])), file=outfile)
 
         for j in range(natoms):
-            print >> outfile2, names[j], \
+            print(names[j], \
                 str(units.unit_to_user('force', unit2, f[i, 3 * j])), \
                 str(units.unit_to_user('force', unit2, f[i, 3 * j + 1])), \
-                str(units.unit_to_user('force', unit2, fRefactoringTool: No changes to ./ipi/utils/mathtools.py
RefactoringTool: Refactored ./ipi/utils/messages.py
RefactoringTool: Refactored ./ipi/utils/mintools.py
RefactoringTool: No changes to ./ipi/utils/nmtransform.py
RefactoringTool: No changes to ./ipi/utils/prng.py
RefactoringTool: Refactored ./ipi/utils/softexit.py
RefactoringTool: Refactored ./ipi/utils/sparse.py
RefactoringTool: No changes to ./ipi/utils/units.py
RefactoringTool: Refactored ./ipi/utils/io/__init__.py
[i, 3 * j + 2]))
+                str(units.unit_to_user('force', unit2, f[i, 3 * j + 2])), file=outfile2)
     outfile.close()
     outfile2.close()
 
--- ./ipi/utils/messages.py	(original)
+++ ./ipi/utils/messages.py	(refactored)
@@ -147,4 +147,4 @@
         return
     if verbosity.trace:
         traceback.print_stack(file=sys.stdout)
-    print(" !W! " + text)
+    print((" !W! " + text))
--- ./ipi/utils/mintools.py	(original)
+++ ./ipi/utils/mintools.py	(refactored)
@@ -669,7 +669,7 @@
     if (min_d > 0.0):
 
         if(neg != 0):
-            print "problem in 'find'!!!"
+            print("problem in 'find'!!!")
         if (np.linalg.norm(DXE) < tr):
             DX = np.dot(w, DXE)
             DX = DX.reshape(shape)
@@ -694,7 +694,7 @@
             lamb_min = max(lamb, lamb_min)
 
         if dy > 0.0 or lamb_min > lamb_max:
-            print 'Problem in find. II'
+            print('Problem in find. II')
 
         lamb = lamb - y / dy
         if lamb <= lamb_min or lamb >= lamb_max:
@@ -743,7 +743,7 @@
         info(" @MINIMIZE: moving uphill, resetting search direction! ", verbosity.debug)
         d0 = g0 / np.sqrt(np.dot(g0.flatten(), g0.flatten()))
 
-    print "@ GEOP step ", big_step
+    print("@ GEOP step ", big_step)
     # Perform approximate line minimization in direction d0
     x, u, g = min_approx(fdf, x0, fdf0, d0, big_step, tol, itmax)
 
@@ -1134,7 +1134,7 @@
             info(" @MINIMIZE: Scaled step size", verbosity.debug)
 
         x = np.add(x0, d0)
-        print "step size:", np.sqrt(np.dot(d0.flatten(), d0.flatten()))
+        print("step size:", np.sqrt(np.dot(d0.flatten(), d0.flatten())))
         fx, g = fdf(x)
 
     info(" @MINIMIZE: Started L-BFGS", verbosity.debug)
--- ./ipi/utils/softexit.py	(original)
+++ ./ipi/utils/softexit.py	(refactored)
@@ -77,7 +77,7 @@
            message: The message to output to standard output.
         """
 
-        print "SOFTEXIT CALLED FROM THREAD", threading.currentThread(), message
+        print("SOFTEXIT CALLED FROM THREAD", threading.currentThread(), message)
         if not self.triggered:    # avoid double calls from different threads
             self.exiting = True
             self.triggered = True
@@ -90,7 +90,7 @@
                 try:
                     f()
                 except RuntimeError as err:
-                    print "Error running emergency softexit, ", err
+                    print("Error running emergency softexit, ", err)
                     pass
 
             self.exiting = False  # emergency is over, signal we can be relaxed
--- ./ipi/utils/sparse.py	(original)
+++ ./ipi/utils/sparse.py	(refactored)
@@ -108,7 +108,7 @@
     # Convert to dense format
     def toarray(self):
         nparray = np.zeros((self.m, self.n), dtype=self.a.dtype)
-        for row in xrange(self.m):
+        for row in range(self.m):
             nparray[row, self.ja[self.ia[row]:self.ia[row + 1]]] = \
                 self.a[self.ia[row]:self.ia[row + 1]]
         return nparray
@@ -144,7 +144,7 @@
     # Convert to dense format
     def toarray(self):
         nparray = np.zeros((self.m, self.n), dtype=self.a.dtype)
-        for col in xrange(self.n):
+        for col in range(self.n):
             nparray[self.ja[self.ia[col]:self.ia[col + 1]], col] = \
                 self.a[self.ia[col]:self.ia[col + 1]]
         return nparray
--- ./ipi/utils/io/__init__.py	(original)
+++ ./ipi/utils/io/__init__.py	(refactored)
@@ -13,7 +13,7 @@
 
 import sys
 import os
-import cStringIO
+import io
 import numpy as np
 
 from ipi.utils.messages import info, verbosity
@@ -54,13 +54,13 @@
             mode = mode_map[mode]
         module = importlib.import_module("ipi.utils.io.backends.io_%s" % mode)
     except ImportError:
-        print "Error: mode %s is not supported." % mode
+        print("Error: mode %s is not supported." % mode)
         sys.exit()
 
     try:
         func = getattr(module, io_map[io] % mode)
     except KeyError:
-        print "Error: io %s is not supported with mode %s." % (io, mode)
+        print("Error: io %s iRefactoringTool: Refactored ./ipi/utils/io/io_units.py
RefactoringTool: No changes to ./ipi/utils/io/backends/__init__.py
RefactoringTool: No changes to ./ipi/utils/io/backends/io_binary.py
RefactoringTool: No changes to ./ipi/utils/io/backends/io_json.py
RefactoringTool: No changes to ./ipi/utils/io/backends/io_pdb.py
RefactoringTool: Refactored ./ipi/utils/io/backends/io_xyz.py
RefactoringTool: No changes to ./ipi/utils/io/inputs/__init__.py
RefactoringTool: Refactored ./ipi/utils/io/inputs/io_xml.py
RefactoringTool: Refactored ./ipi_tests/pdb_generator.py
s not supported with mode %s." % (io, mode))
         sys.exit()
 
     return func
@@ -332,7 +332,7 @@
 
 
 def netstring_encoded_savez(ofile, compressed=True, **named_objs):
-    output = cStringIO.StringIO()
+    output = io.StringIO()
     if compressed:
         # np.savez_compressed(output,*unnamed_objs,**named_objs)
         np.savez_compressed(output, **named_objs)
@@ -358,7 +358,7 @@
     if not ifile.read(1) == ",":
         raise ValueError("Invalid netstring delimiter")
 
-    istr = cStringIO.StringIO(content)
+    istr = io.StringIO(content)
     npz = np.load(istr)
     rdic = {}
     for a in npz.files:
--- ./ipi/utils/io/io_units.py	(original)
+++ ./ipi/utils/io/io_units.py	(refactored)
@@ -18,7 +18,7 @@
 cell_unit_re = re.compile(r'cell\{([A-Za-z_]*)\}')       # cell unit pattern
 traj_dict = Traj().traj_dict                             # trajectory dictionary
 traj_re = [re.compile('%s%s' % (key, r'\{[A-Za-z_]*\}'))
-           for key in traj_dict.keys()]  # trajectory patterns
+           for key in list(traj_dict.keys())]  # trajectory patterns
 
 
 def auto_units(comment="", dimension="automatic", units="automatic", cell_units="automatic", mode="xyz"):
@@ -38,8 +38,8 @@
     if comment != "":  # but they can be overridden by a special comment line
         # tries to guess units from the input
         # Extracting trajectory units
-        is_comment_useful = filter(None, [key.search(comment.strip())
-                                          for key in traj_re])
+        is_comment_useful = [_f for _f in [key.search(comment.strip())
+                                          for key in traj_re] if _f]
         if len(is_comment_useful) > 0:
             traj = is_comment_useful[0].group()[:-1].split('{')
             auto_dimension, auto_units = traj_dict[traj[0]]['dimension'], traj[1]
--- ./ipi/utils/io/backends/io_xyz.py	(original)
+++ ./ipi/utils/io/backends/io_xyz.py	(refactored)
@@ -87,11 +87,11 @@
     """
 
     try:
-        natoms = int(filedesc.next())
+        natoms = int(next(filedesc))
     except (StopIteration, ValueError):
         raise EOFError
 
-    comment = filedesc.next()
+    comment = next(filedesc)
 
     # Extracting cell
     cell = [key.search(comment) for key in cell_re]
--- ./ipi/utils/io/inputs/io_xml.py	(original)
+++ ./ipi/utils/io/inputs/io_xml.py	(refactored)
@@ -101,7 +101,7 @@
         """
 
         # creates a new node
-        newnode = xml_node(attribs=dict((k, attrs[k]) for k in attrs.keys()), name=name, fields=[])
+        newnode = xml_node(attribs=dict((k, attrs[k]) for k in list(attrs.keys())), name=name, fields=[])
         # adds it to the list of open nodes
         self.open.append(newnode)
         # adds it to the list of fields of the parent tag
@@ -202,7 +202,7 @@
     rstr = ""
     if not name == "":
         rstr = indent + "<" + name;
-        for a, v in xml.attribs.iteritems():
+        for a, v in xml.attribs.items():
             rstr += " " + a + "='" + v + "'"
         rstr += ">"
 
@@ -426,7 +426,7 @@
 
     rdict = {}
     for s in rlist:
-        rtuple = map(mystrip, s.split(key_split))
+        rtuple = list(map(mystrip, s.split(key_split)))
         if not len(rtuple) == 2:
             raise ValueError("Format for a key:value format is wrong for item " + s)
         rdict[rtuple[0]] = rtuple[1]
--- ./ipi_tests/pdb_generator.py	(original)
+++ ./ipi_tests/pdb_generator.py	(refactored)
@@ -6,7 +6,7 @@
 import numpy as np
 from ipi.utils.units import Elements
 from copy import copy
-all_elem = Elements.mass_list.keys()
+all_elem = list(Elements.mass_list.keys())
 at_names = None
 
 
@@ -25,7 +25,7 @@
     if not comment.endswith('\n'):
         comment += '\n'
     output = comment
-    for i in xrange(natoms):
+    for i in range(natoms):
         output += ('%-6s'    # Record name                       1- 6  1
                    '%5i'    # Atom serial number                7-11  2
                    ' '      # Space                               12
@@ -54,7 +54,7 @@
     """ Generate a fake xyz trajectory. Atoms and coordinaRefactoringTool: Refactored ./ipi_tests/xyz_generator.py
RefactoringTool: Refactored ./ipi_tests/engine/test_initializer.py
RefactoringTool: No changes to ./ipi_tests/engine/test_properties.py
RefactoringTool: Refactored ./ipi_tests/utils/test_depend.py
RefactoringTool: Refactored ./ipi_tests/utils/io/backends/test__init__.py
RefactoringTool: Refactored ./ipi_tests/utils/io/backends/test_io_units.py
tes are random.
     """
     output, xyz, all_names = xyz_rand(natoms, comment)
-    for _ in xrange(nframe - 1):
+    for _ in range(nframe - 1):
         raw_output = xyz_rand(natoms, comment, names=True)
         output += raw_output[0]
         xyz = np.concatenate([xyz, raw_output[1]])
@@ -77,4 +77,4 @@
 
     # Fast autocheck... if the test is wrong itself... it is bad ;)
     natoms = 100
-    print xyz_rand(natoms, '')[0]
+    print(xyz_rand(natoms, '')[0])
--- ./ipi_tests/xyz_generator.py	(original)
+++ ./ipi_tests/xyz_generator.py	(refactored)
@@ -6,7 +6,7 @@
 import numpy as np
 from ipi.utils.units import Elements
 from copy import copy
-all_elem = Elements.mass_list.keys()
+all_elem = list(Elements.mass_list.keys())
 at_names = None
 
 
@@ -27,7 +27,7 @@
     if not comment.endswith('\n'):
         comment += '\n'
     output += comment
-    for i in xrange(natoms):
+    for i in range(natoms):
         output += '%8s %12.5e %12.5e %12.5e\n' % \
             (at_names[i], xyz[i, 0], xyz[i, 1], xyz[i, 2])
     return (output, xyz.flatten(), copy(at_names))
@@ -37,7 +37,7 @@
     """ Generate a fake xyz trajectory. Atoms and coordinates are random.
     """
     output, xyz, all_names = xyz_rand(natoms, comment)
-    for _ in xrange(nframe - 1):
+    for _ in range(nframe - 1):
         raw_output = xyz_rand(natoms, comment, names=True)
         output += raw_output[0]
         xyz = np.concatenate([xyz, raw_output[1]])
--- ./ipi_tests/engine/test_initializer.py	(original)
+++ ./ipi_tests/engine/test_initializer.py	(refactored)
@@ -57,7 +57,7 @@
         masses[_ii] = Elements.mass(_at)
 
     ratoms = []
-    for _fr in xrange(frames):
+    for _fr in range(frames):
         ratoms.append(Atoms(natoms))
         ratoms[-1].q = xyz[_fr * natoms * 3:3 * (_fr + 1) * natoms] * units_conv_at
         ratoms[-1].m = masses[_fr * natoms:(_fr + 1) * natoms]
--- ./ipi_tests/utils/test_depend.py	(original)
+++ ./ipi_tests/utils/test_depend.py	(refactored)
@@ -63,22 +63,22 @@
 
 
 def threadA(Aobj):
-    for i in xrange(10000):
+    for i in range(10000):
         Aobj.scalar = np.sqrt(i)
         time.sleep(0.0001)
 
 
 def threadB(Aobj, Bobj):
-    for i in xrange(10000):
+    for i in range(10000):
         Aobj.scalar = i
         time.sleep(0.0001)
         Bobj.scalar = 4 + i
         Bobj.vector = np.sqrt(i)
-        print Bobj.dvector[0]
+        print(Bobj.dvector[0])
 
 
 # myB.dvector should always contain B.vector*B.scalar- A.scalar
-print myB.dvector - mbv * (mbs - mas)
+print(myB.dvector - mbv * (mbs - mas))
 
 # now try with threads
 ta = threading.Thread(target=threadA, name="TA", kwargs={"Aobj": myA})
@@ -90,4 +90,4 @@
 ta.join()
 tb.join()
 
-print myB.dvector - myB.vector * (myB.scalar - myA.scalar)
+print(myB.dvector - myB.vector * (myB.scalar - myA.scalar))
--- ./ipi_tests/utils/io/backends/test__init__.py	(original)
+++ ./ipi_tests/utils/io/backends/test__init__.py	(refactored)
@@ -77,4 +77,4 @@
     npt.assert_almost_equal(expected_q * unit_conv_q, returned_q, 5)
     npt.assert_almost_equal(expected_cell.flatten() * unit_conv_cell, returned_cell, 5)
     if output_type.strip() == 'objects':
-        assert all([expected_names[_ii] == returned_names[_ii] for _ii in xrange(len(expected_names))])
+        assert all([expected_names[_ii] == returned_names[_ii] for _ii in range(len(expected_names))])
--- ./ipi_tests/utils/io/backends/test_io_units.py	(original)
+++ ./ipi_tests/utils/io/backends/test_io_units.py	(refactored)
@@ -40,9 +40,9 @@
 
     masses = []
     _count = 0
-    print 'atom_names', atom_names
+    print('atom_names', atom_names)
     for _at in atom_names:
-        print Elements.mass(_at), _count, _at
+        print(Elements.mass(_at), _count, _at)
         masses.append(Elements.mass(_at))
         _count += 1
 
@@ -50,7 +50,7 @@
     type(masses)
     res = testing.process_units(comment, cell.copy(), xyz.copy(), np.array(atom_names).copy(), np.array(masses).copy(), output=output)
 
-    print xyz, res['data']
+    print(xyz, res['data'])
     npt.assert_array_almost_equal(res['dRefactoringTool: Refactored ./ipi_tests/utils/io/backends/test_io_xyz.py
RefactoringTool: No changes to ./tests/test_docs.py
RefactoringTool: Refactored ./tools/py/Instanton_interpolation.py
ata'], xyz * conver_xyz, 5)
     npt.assert_array_almost_equal(res['masses'], masses, 5)
     npt.assert_array_almost_equal(res['cell'], cell * conver_cell, 5)
@@ -73,7 +73,7 @@
     _count = 0
 
     for _at in atom_names:
-        print Elements.mass(_at), _count, _at
+        print(Elements.mass(_at), _count, _at)
         masses.append(Elements.mass(_at))
         _count += 1
 
--- ./ipi_tests/utils/io/backends/test_io_xyz.py	(original)
+++ ./ipi_tests/utils/io/backends/test_io_xyz.py	(refactored)
@@ -71,7 +71,7 @@
     if check_comment:
         genh = np.array(check_comment.group(1).split()[:9], float).reshape((3, 3))
         invgenh = np.linalg.inv(genh)
-        for _ui in xrange(natoms * frames):
+        for _ui in range(natoms * frames):
             _uu = np.array([xyz[3 * _ui], xyz[3 * _ui + 1], xyz[3 * _ui + 2]])
             _us = np.dot(_uu, invgenh)
             _uu = np.dot(expected_cell, _us)
@@ -87,7 +87,7 @@
     filedesc, xyz, atom_names, \
         natoms, frames, comment, expected_cell, precision = create_random_xyz_traj_to_read
 
-    for _fr in xrange(frames):
+    for _fr in range(frames):
 
         tcomment, tcell, tqatoms, tnames, tmasses = io_xyz.read_xyz(filedesc)
 
@@ -152,7 +152,7 @@
     cell_list = []
     atoms_list = []
 
-    for _fr in xrange(frames):
+    for _fr in range(frames):
         cell = Cell(expected_cell)
         atoms = Atoms(natoms)
         atoms.q[:] = xyz[_fr * natoms * 3:(_fr + 1) * natoms * 3]
--- ./tools/py/Instanton_interpolation.py	(original)
+++ ./tools/py/Instanton_interpolation.py	(refactored)
@@ -70,7 +70,7 @@
         if (os.path.exists(input_geo)):
             ipos = open(input_geo, "r")
         else:
-            print("We can't find".format(input_geo))
+            print(("We can't find".format(input_geo)))
             sys.exit()
 
         pos = list()
@@ -94,7 +94,7 @@
         if (os.path.exists(chk)):
             simulation = Simulation.load_from_xml(chk, custom_verbosity='low', request_banner=False)
         else:
-            print("We can't find".format(chk))
+            print(("We can't find".format(chk)))
             sys.exit()
         cell = simulation.syslist[0].cell
         beads = simulation.syslist[0].motion.beads.copy()
@@ -104,8 +104,8 @@
         atom = beads._blist[0]
 
     print(' ')
-    print('We have a half ring polymer made of {} beads and {} atoms.'.format(nbeads, natoms))
-    print('We will expand the ring polymer to get a half polymer of {} beads.'.format(nbeadsNew))
+    print(('We have a half ring polymer made of {} beads and {} atoms.'.format(nbeads, natoms)))
+    print(('We will expand the ring polymer to get a half polymer of {} beads.'.format(nbeadsNew)))
 
     # Make the rpc step (standar)
     # q2 = np.concatenate((q, np.flipud(q)), axis=0)   # Compose the full ring polymer.
@@ -116,7 +116,7 @@
     # new_q = rpc.b1tob2(q )
 
     # Make the rpc step (open path)
-    rpc = nm_rescale(nbeads, nbeadsNew, np.asarray(range(natoms)))
+    rpc = nm_rescale(nbeads, nbeadsNew, np.asarray(list(range(natoms))))
     new_q = rpc.b1tob2(q)
 
     # Print
@@ -130,7 +130,7 @@
     print('The new Instanton geometry (half polymer) was generated')
     print('Check new_instanton.xyz')
     print('')
-    print("Don't forget to change the number of beads to the new value ({}) in your input file".format(nbeadsNew))
+    print(("Don't forget to change the number of beads to the new value ({}) in your input file".format(nbeadsNew)))
     print('when starting your new simulation with an increased number of beads.')
     print('')
 
@@ -140,7 +140,7 @@
         try:
             hess = open(input_hess, "r")
         except:
-            print("We can't find".format(input_hess))
+            print(("We can't find".format(input_hess)))
             sys.exit()
         h = np.zeros((natoms * 3)**2 * nbeads)
         aux = hess.readline().split()
@@ -161,7 +161,7 @@
             print("We don't have a hessian so there is nothing more to do")
             sys.exit()
 
-    print('The new hessian is {} x {}.'.format(3 * natoms, natoRefactoringTool: Refactored ./tools/py/Instanton_postproc.py
ms * 3 * nbeadsNew))
+    print(('The new hessian is {} x {}.'.format(3 * natoms, natoms * 3 * nbeadsNew)))
     out = open("new_hessian.dat", "w")
 
     print('Creating matrix... ')
@@ -199,7 +199,7 @@
     print('')
     print('Remeber to adapt/add the following line in your input:')
     print('')
-    print(" <hessian mode='file' shape='({}, {})' >hessian.dat</hessian>".format(3 * natoms, natoms * 3 * nbeadsNew))
+    print((" <hessian mode='file' shape='({}, {})' >hessian.dat</hessian>".format(3 * natoms, natoms * 3 * nbeadsNew)))
     print('')
 
 sys.exit()
--- ./tools/py/Instanton_postproc.py	(original)
+++ ./tools/py/Instanton_postproc.py	(refactored)
@@ -180,7 +180,7 @@
 # -----END of some functions-----------------
 
 # -----READ---------------------------------
-print('\nWe are ready to start. Reading {} ... (This can take a while)'.format(inputt))
+print(('\nWe are ready to start. Reading {} ... (This can take a while)'.format(inputt)))
 
 from ipi.utils.io.inputs.io_xml import xml_parse_file
 import ipi.inputs.simulation as isimulation
@@ -203,8 +203,8 @@
         sys.exit()
 
 if case != 'instanton' and nbeads > 1:
-    print('Incompatibility between case and nbeads in {}.'.format(inputt))
-    print('case {} , beads {}'.format(case, nbeads))
+    print(('Incompatibility between case and nbeads in {}.'.format(inputt)))
+    print(('case {} , beads {}'.format(case, nbeads)))
     sys.exit()
 
 
@@ -262,7 +262,7 @@
             print(' You can generate that file using this script in the case reactant.')
             sys.exit()
 
-        print('Our linear polymer has  {}'.format(nbeads))
+        print(('Our linear polymer has  {}'.format(nbeads)))
         pos = beads.q
         m3 = beads.m3
         omega2 = (temp * nbeads * kb / hbar) ** 2
@@ -282,18 +282,18 @@
 beta = 1.0 / (kb * temp)
 betaP = 1.0 / (kb * (nbeads) * temp)
 
-print('\nTemperature: {} K'.format(temp / K2au))
-print('NBEADS: {}'.format(nbeads))
-print('atoms:  {}'.format(natoms))
-print('ASR:    {}'.format(natoms))
-print('1/(betaP*hbar) = {:8.5f}'.format((1 / (betaP * hbar))))
+print(('\nTemperature: {} K'.format(temp / K2au)))
+print(('NBEADS: {}'.format(nbeads)))
+print(('atoms:  {}'.format(natoms)))
+print(('ASR:    {}'.format(natoms)))
+print(('1/(betaP*hbar) = {:8.5f}'.format((1 / (betaP * hbar)))))
 
 if not quiet:
     print('Diagonalization ... \n\n')
     d, w, detI = clean_hessian(h, pos, natoms, nbeads, m, m3, asr, mofi=True)
     print("Final lowest 10 frequencies (cm^-1)")
     d10 = np.array2string(np.sign(d[0:10]) * np.absolute(d[0:10]) ** 0.5 / cm2au, precision=2, max_line_width=100, formatter={'float_kind': lambda x: "%.2f" % x})
-    print('{}'.format(d10))
+    print(('{}'.format(d10)))
 
 if case == 'reactant':
     Qtras = ((np.sum(m)) / (2 * np.pi * beta * hbar**2))**1.5
@@ -318,9 +318,9 @@
     # logQvib    = -np.sum( np.log( 2*np.sinh( (beta*hbar*np.sqrt(d)/2.0) )  ))   #Limit n->inf
     logQvib_rp = -get_rp_freq(d, nbeadsR, temp)
 
-    print('\nWe are done. Reactants. Nbeads {}'.format(nbeadsR))
-    print('{:14s} | {:8s} | {:8s}'.format('Qtras(bohr^-3)', 'Qrot', 'logQvib_rp'))
-    print('{:14.3f} | {:8.3f} |{:8.3f}\n'.format(Qtras, Qrot, logQvib_rp))
+    print(('\nWe are done. Reactants. Nbeads {}'.format(nbeadsR)))
+    print(('{:14s} | {:8s} | {:8s}'.format('Qtras(bohr^-3)', 'Qrot', 'logQvib_rp')))
+    print(('{:14.3f} | {:8.3f} |{:8.3f}\n'.format(Qtras, Qrot, logQvib_rp)))
     print('A file with the frequencies in atomic units was generated\n')
 
 elif case == 'TS':
@@ -336,11 +336,11 @@
     U = (pots.sum() - V0)
 
     print('\nWe are done. TS')
-    print('Partition functions at {} K'.format(temp / K2au))
-    print('\nQtras: {}'.format(Qtras))
-    print('Qrot: {}'.format(Qrot))
-    print('logQvib: {}'.format(logQvib))
-    print('Potential energy at TS:  {} eV, V/kBT {}\n'.format(U / eV2au, U / (kb * temp)))
+    print(('Partition functions at {} K'.format(temp / K2au)))
+    print(('\nQtras: {}'.format(Qtras)))
+    print(('Qrot: {}'.format(Qrot)))
+    print(('logQvib: {}'.format(logQvib)))
+    print(('Potential energy at TS:  {} eV, V/kBT {}\n'.format(U / eV2au, U / (kb * temp))))
 
 elif case == 'instanton':
 
@@ -353,11 +353,11 @@
             Qrot = 1.0
 
         if not quiet:
-            print('Deleted frequency: {:8.3f} cm^-1'.format((np.sign(d[1]) * np.absolute(d[1]) ** 0.5 / cm2au)))
+            print(('Deleted frequency: {:8.3f} cm^-1'.format((np.sign(d[1]) * np.absolute(d[1]) ** 0.5 / cm2au))))
             if asr != 'poly':
                 print('WARNING asr != poly')
                 print('First 10 eigenvalues')
-                print('{}'.format(np.sign(d[0:10]) * np.absolute(d[0:10]) ** 0.5 / cm2au))
+                print(('{}'.format(np.sign(d[0:10]) * np.absolute(d[0:10]) ** 0.5 / cm2au)))
                 print("Please check that this you don't have any unwanted zero frequency")
             logQvib = -np.sum(np.log(betaP * hbar * np.sqrt(np.absolute(np.delete(d, 1))))) + 6 * np.log(nbeads) + np.log(nbeads)
         else:
@@ -368,9 +368,9 @@
         action1 = (2 * pots.sum() * factor - nbeads * V0) * 1. / (temp * nbeads * kb)
         action2 = spring_pot(nbeads, pos, omega2, m3) / (temp * nbeads * kb)
 
-        print('\nWe are done. Instanton rate. Nbeads {} (diff only {})'.format(nbeads, nbeads / 2))
-        print('   {:8s} {:8s}  | {:11s} | {:11s} | {:11s} | {:8s} ( {:8s},{:8s} ) |'.format('BN', '(BN*N)', 'Qt(bohr^-3)', 'Qrot', 'log(Qvib*N)', 'S/hbar', 'S1/hbar', 'S2/hbar'))
-        print('{:8.3f} ( {:8.3f} ) | {:11.3f} | {:11.3f} | {:11.3f} | {:8.3f} ( {:8.3f} {:8.3f} ) |'.format(BN, BN * nbeads, Qtras, Qrot, logQvib, (action1 + action2), action1, action2))
+        print(('\nWe are done. Instanton rate. Nbeads {} (diff only {})'.format(nbeads, nbeads / 2)))
+        print(('   {:8s} {:8s}  | {:11s} | {:11s} | {:11s} | {:8s} ( {:8s},{:8s} ) |'.format('BN', '(BN*N)', 'Qt(bohr^-3)', 'Qrot', 'log(Qvib*N)', 'S/hbar', 'S1/hbar', 'S2/hbar')))
+        print(('{:8.3f} ( {:8.3f} ) | {:11.3f} | {:11.3f} | {:11.3f} | {:8.3f} ( {:8.3f} {:8.3f} ) |'.format(BN, BN * nbeads, Qtras, Qrot, logQvib, (action1 + action2), action1, action2)))
         print('\n\n')
 
     elif mode == 'splitting':
@@ -379,8 +379,8 @@
         d_min = np.zeros(natoms * 3)
         aux = out.readline().split()
         if len(aux) != (natoms * 3):
-            print('We are expecting {} frequencies.'.format((natoms * 3 - 6)))
-            print('instead we have read  {}'.format(len(aux)))
+            print(('We are expecting {} frequencies.'.format((natoms * 3 - 6))))
+            print(('instead we have read  {}'.format(len(aux))))
         for i in range((natoms * 3)):
             d_min[i] = float(aux[i])
         d_min = d_min.reshape((natoms * 3))
@@ -408,19 +408,19 @@
         # cm2au= (2 * np.pi * 3e10 * 2.4188843e-17)
 
         print('\n\nWe are done')
-        print('Nbeads {}, betaP {} a.u.,hbar {} a.u'.format((nbeads, betaP, hbar)))
+        print(('Nbeads {}, betaP {} a.u.,hbar {} a.u'.format((nbeads, betaP, hbar))))
         print('')
-        print('V0  {} eV ( {} Kcal/mol) '.format((V0 / eV2au, V0 / cal2au / 1000)))
-        print('S1/hbar {} ,S2/hbar {} ,S/hbar {}'.format(action1 / hbar, action2 / hbar, action / hbar))
-        print('BN {} a.u.'.format(BN))
-        print('BN/(hbar^2 * betaN)  {}  (should be same as S/hbar) '.format((BN / ((hbar**2) * betaP))))
+        print(('V0  {} eV ( {} Kcal/mol) '.format((V0 / eV2au, V0 / cal2au / 1000))))
+        print(('S1/hbar {} ,S2/hbar {} ,S/hbar {}'.format(action1 / hbar, action2 / hbar, action / hbar)))
+        print(('BN {} a.u.'.format(BN)))
+        print(('BN/(hbar^2 * betaN)  {}  (should be same as S/hbar) '.format((BN / ((hbar**2) * betaP)))))
         print('')
         if quiet:
             print ('phi is not computed because you specified the quiet option')
-            print ('We can provied only Tetaphi which value is {} a.u. '.format(tetaphi))
+            print(('We can provied only Tetaphi which value is {} a.u. '.format(tetaphi)))
         else:
-            print('phi {} a.u.   Teta {} a.u. '.format(phi, tetaphi / phi))
-          RefactoringTool: Refactored ./tools/py/a2b.py
RefactoringTool: No changes to ./tools/py/bin2xyz.py
RefactoringTool: Refactored ./tools/py/contract-trajectory.py
RefactoringTool: Refactored ./tools/py/effective_temperatures.py
RefactoringTool: Refactored ./tools/py/energies_ppi.py
  print('Tunnelling splitting matrix element (h)  {} a.u ({} cm^-1)'.format((h, h / cm2au)))
+            print(('phi {} a.u.   Teta {} a.u. '.format(phi, tetaphi / phi)))
+            print(('Tunnelling splitting matrix element (h)  {} a.u ({} cm^-1)'.format((h, h / cm2au))))
     else:
         print('We can not recongnize the mode.')
         sys.exit()
--- ./tools/py/a2b.py	(original)
+++ ./tools/py/a2b.py	(refactored)
@@ -24,7 +24,7 @@
 
 cell_unit_re = re.compile(r'cell\{([A-Za-z_]*)\}')       # cell unit pattern
 traj_dict = Traj().traj_dict                             # trajectory dictionary
-traj_re = [re.compile('%s%s' % (key, r'\{[A-Za-z_]*\}')) for key in traj_dict.keys()]  # trajectory patterns
+traj_re = [re.compile('%s%s' % (key, r'\{[A-Za-z_]*\}')) for key in list(traj_dict.keys())]  # trajectory patterns
 verbosity.level = "low"
 
 
@@ -56,7 +56,7 @@
 
     is_comment_useful = []
     if comment != "":
-        is_comment_useful = filter(None, [key.search(comment.strip()) for key in traj_re])
+        is_comment_useful = [_f for _f in [key.search(comment.strip()) for key in traj_re] if _f]
         if len(is_comment_useful) > 0:
             traj = is_comment_useful[0].group()[:-1].split('{')
             key, auto_dimension, auto_units = traj[0], traj_dict[traj[0]]['dimension'], traj[1]
--- ./tools/py/contract-trajectory.py	(original)
+++ ./tools/py/contract-trajectory.py	(refactored)
@@ -1,6 +1,6 @@
 #!/usr/bin/env python2
 
-from __future__ import print_function
+
 
 import os
 import sys
@@ -69,7 +69,7 @@
     while True:
         try:
             # Get the frames for all beads.
-            frames = [trj.next() for trj in trjs_in]
+            frames = [next(trj) for trj in trjs_in]
         except StopIteration:
             # Stop when any of the trajectories runs out of frames.
             break
--- ./tools/py/effective_temperatures.py	(original)
+++ ./tools/py/effective_temperatures.py	(refactored)
@@ -45,17 +45,17 @@
     nbeads = len(fns_for)
 
     # print some information
-    print 'temperature = {:f} K'.format(float(temp))
-    print
-    print 'number of beads = {:d}'.format(nbeads)
-    print
-    print 'forces file names:'
+    print('temperature = {:f} K'.format(float(temp)))
+    print()
+    print('number of beads = {:d}'.format(nbeads))
+    print()
+    print('forces file names:')
     for fn_for in fns_for:
-        print '{:s}'.format(fn_for)
-    print
-    print 'output file name:'
-    print fn_out
-    print
+        print('{:s}'.format(fn_for))
+    print()
+    print('output file name:')
+    print(fn_out)
+    print()
 
     # open input and output files
     ifor = [open(fn, "r") for fn in fns_for]
@@ -72,7 +72,7 @@
     while True:  # Reading input files and calculating PPI correction
 
         if ifr % 100 == 0:
-            print '\rProcessing frame {:d}'.format(ifr),
+            print('\rProcessing frame {:d}'.format(ifr), end=' ')
             sys.stdout.flush()
 
         try:
--- ./tools/py/energies_ppi.py	(original)
+++ ./tools/py/energies_ppi.py	(refactored)
@@ -67,24 +67,24 @@
     # check that we found the same number of positions and forces files
     nbeads = len(fns_pos)
     if nbeads != len(fns_for):
-        print fns_pos
-        print fns_for
+        print(fns_pos)
+        print(fns_for)
         raise ValueError("Mismatch between number of input files for forces and positions.")
 
     # print some information
-    print 'temperature = {:f} K'.format(float(temp))
-    print
-    print 'number of beads = {:d}'.format(nbeads)
-    print
-    print 'positions and forces file names:'
+    print('temperature = {:f} K'.format(float(temp)))
+    print()
+    print('number of beads = {:d}'.format(nbeads))
+    print()
+    print('positions and forces file names:')
     for fn_pos, fn_for in zip(fns_pos, fns_for):
-        print '{:s}   {:s}'.format(fn_pos, fn_for)
-    print
-    print 'potential energy file: {:s}'.format(fns_iU)
-    print
-    print 'output file name:'
-    print fn_out_en
-    print
+        print('{:s}   {:s}'.format(fn_pos, fn_for))RefactoringTool: Refactored ./tools/py/energy_ppi.py
RefactoringTool: Refactored ./tools/py/estmod_example.py
RefactoringTool: No changes to ./tools/py/fixcom.py
RefactoringTool: Refactored ./tools/py/get_np.py
RefactoringTool: Refactored ./tools/py/get_np_rad.py

+    print()
+    print('potential energy file: {:s}'.format(fns_iU))
+    print()
+    print('output file name:')
+    print(fn_out_en)
+    print()
 
     # open input and output files
     ipos = [open(fn, "r") for fn in fns_pos]
@@ -121,7 +121,7 @@
     while True:  # Reading input files and calculating PPI correction
 
         if ifr % 100 == 0:
-            print '\rProcessing frame {:d}'.format(ifr),
+            print('\rProcessing frame {:d}'.format(ifr), end=' ')
             sys.stdout.flush()
 
         try:
--- ./tools/py/energy_ppi.py	(original)
+++ ./tools/py/energy_ppi.py	(refactored)
@@ -130,7 +130,7 @@
             f2ePA_av += f2ePA
             eVir_av += eVir
             ifr += 1
-            print(ifr - skipSteps)  # Printing current time frame (excluding thermalization)
+            print((ifr - skipSteps))  # Printing current time frame (excluding thermalization)
 
             dE = (3.0 * Constants.kb * temperature + ePA_av / float(ifr - skipSteps)) * f2_av / float(ifr - skipSteps) - \
                 f2ePA_av / float(ifr - skipSteps)
--- ./tools/py/estmod_example.py	(original)
+++ ./tools/py/estmod_example.py	(refactored)
@@ -66,7 +66,7 @@
     """
     nndof, npl = qsum.shape
     ans = np.zeros((npl,) + Bshape)
-    for i in xrange(npl):
+    for i in range(npl):
         ans[i] = 0.0  # Change this bit
     return ans
 
@@ -82,7 +82,7 @@
     """
     nndof, npl = qsum.shape
     ans = np.zeros((npl,) + Ashape)
-    for i in xrange(npl):
+    for i in range(npl):
         ans[i] = 0.0  # Change this bit
     return ans
 
@@ -95,7 +95,7 @@
     """
     nndof, npl = qsum.shape
     ans = np.zeros((npl,) + Ashape + (ndof,))
-    for i in xrange(npl):
+    for i in range(npl):
         ans[i] = 0.0  # Change this bit
     return ans
 
@@ -108,7 +108,7 @@
     """
     nndof, npl = qsum.shape
     ans = np.zeros((npl,) + Ashape + (ndof, ndof))
-    for i in xrange(npl):
+    for i in range(npl):
         ans[i] = 0.0  # Change this bit
     return ans
 
@@ -120,7 +120,7 @@
     """
     npts, npl = A.shape[:2]
     tcf = np.zeros(npts)
-    for i in xrange(npl):
+    for i in range(npl):
         tcf[:] += 0.0  # Change this bit
     return tcf
 
--- ./tools/py/get_np.py	(original)
+++ ./tools/py/get_np.py	(refactored)
@@ -77,7 +77,7 @@
     step = np.shape(delta)[0]
     n_block = int(step / bsize)
 
-    for x in xrange(n_block):
+    for x in range(n_block):
         dq = delta[x * bsize: (x + 1) * bsize]
         hx = histo(np.concatenate((dq.T[0], -dq.T[0])), dqxgrid, kernel, 0, np.sqrt(T * P * m))
         hy = histo(np.concatenate((dq.T[1], -dq.T[1])), dqygrid, kernel, 0, np.sqrt(T * P * m))
@@ -145,9 +145,9 @@
         psqmedz = psqmedz + np.dot(pzgrid**2, np.asarray(nplistz)[i, :]) / normz
         psqmed2z = psqmed2z + (np.dot(pzgrid**2, np.asarray(nplistz)[i, :]) / normz)**2
 
-    print 'number of blocks', n_block
-    print 'av_px^2', psqmedx / n_block, 'sigmax', np.sqrt((psqmed2x / n_block) - (psqmedx / n_block)**2) / np.sqrt(n_block)
-    print 'av_py^2', psqmedy / n_block, 'sigmay', np.sqrt((psqmed2y / n_block) - (psqmedy / n_block)**2) / np.sqrt(n_block)
-    print 'av_pz^2', psqmedz / n_block, 'sigmaz', np.sqrt((psqmed2z / n_block) - (psqmedz / n_block)**2) / np.sqrt(n_block)
+    print('number of blocks', n_block)
+    print('av_px^2', psqmedx / n_block, 'sigmax', np.sqrt((psqmed2x / n_block) - (psqmedx / n_block)**2) / np.sqrt(n_block))
+    print('av_py^2', psqmedy / n_block, 'sigmay', np.sqrt((psqmed2y / n_block) - (psqmedy / n_block)**2) / np.sqrt(n_block))
+    print('av_pz^2', psqmedz / n_block, 'sigmaz', np.sqrt((psqmed2z / n_block) - (psqmedz / n_block)**2) / np.sqrt(n_block))
 
 if __name__ == '__main__': get_np(*sys.argv[1:])
--- ./tools/py/get_np_rad.py	(original)
+++ ./tools/py/get_np_rad.py	(refactored)
@@ -126,7 +126,7 @@
     mw_P2 = k_param[3]
     P = qpath.shape[1] / 3
 
-    for i in xrange(len(qpath)):
+    for i in range(len(qpath)):
         # Instantaneous position of all the beads of the ring polymer
         q = qpath[i]
 
@@ -233,10 +233,10 @@
     RefactoringTool: Refactored ./tools/py/get_np_vec.py
n_blocks = int(len(qpath) / bsize)
 
     if n_blocks == 0:
-        print "# ERROR: Not enough data to build a block"
+        print("# ERROR: Not enough data to build a block")
         exit()
 
-    for x in xrange(n_blocks):
+    for x in range(n_blocks):
 
         if der == False:
             qpath_block = qpath[x * bsize: (x + 1) * bsize]
@@ -295,7 +295,7 @@
     avg_r2_4pi_h = avg_r2_4pi_h / norm_r2_4pi_h
     err_r2_4pi_h = np.std(np.asarray(r2_4pi_h_list) / (norm_r2_4pi_h / n_blocks), axis=0) / np.sqrt(n_blocks)
     np.savetxt(str(prefix + "4pi_r2_h" + ".data"), np.c_[r, avg_r2_4pi_h, err_r2_4pi_h])
-    print "# Printing the radial distribution of the end-to-end distance :", str(prefix + "4pi_r2_h" + ".data")
+    print("# Printing the radial distribution of the end-to-end distance :", str(prefix + "4pi_r2_h" + ".data"))
 
     # Block averages the radial momentum distributions and estimates errors.
     avg_p2_4pi_np = np.sum(np.asarray(p2_4pi_np_list), axis=0)
@@ -303,12 +303,12 @@
     avg_p2_4pi_np = avg_p2_4pi_np / norm_p2_4pi_np
     err_p2_4pi_np = np.std(np.asarray(p2_4pi_np_list) / (norm_p2_4pi_np / n_blocks), axis=0) / np.sqrt(n_blocks)
     np.savetxt(str(prefix + "4pi_p2_np" + ".data"), np.c_[p, avg_p2_4pi_np, err_p2_4pi_np])
-    print "# Printing the radial distribution of the particle momentum :", str(prefix + "4pi_p2_np" + ".data")
+    print("# Printing the radial distribution of the particle momentum :", str(prefix + "4pi_p2_np" + ".data"))
 
     # Also calulates the average value of p^2.
     avg_avgp2 = np.sum(np.asarray(avgp2_list), axis=0) / norm_p2_4pi_np
     err_avgp2 = np.std(np.asarray(avgp2_list) / (norm_p2_4pi_np / n_blocks), axis=0) / np.sqrt(n_blocks)
-    print "# avg <p2> :", avg_avgp2, "+/-", err_avgp2
+    print("# avg <p2> :", avg_avgp2, "+/-", err_avgp2)
 
 
 if __name__ == '__main__':
--- ./tools/py/get_np_vec.py	(original)
+++ ./tools/py/get_np_vec.py	(refactored)
@@ -10,6 +10,7 @@
 import numpy as np
 import time
 from scipy.interpolate import RegularGridInterpolator
+from functools import reduce
 
 
 def histo3d(qdata, dqxgrid, dqygrid, dqzgrid, ns, cut, invsigma, bsize):
@@ -157,11 +158,11 @@
     if der == False:
         n_block = int(step / bsize)
         if (n_block == 0):
-            print 'not enough data to build a block'
+            print('not enough data to build a block')
             exit()
-        for x in xrange(n_block):
+        for x in range(n_block):
             dq = delta[x * bsize: (x + 1) * bsize]
-            print "# Computing 3D histogram."
+            print("# Computing 3D histogram.")
             h3d = histo3d(np.concatenate((dq, -dq)), dqxgrid, dqygrid, dqzgrid, ns, cut, np.sqrt(T * P * m), bsize)
             h3dlist.append(h3d)
 
@@ -174,7 +175,7 @@
             ygrid = dqygrid
             zgrid = dqzgrid
 
-            print "# NORM OF THE 3D HISTO:", h3d.flatten().sum() * dqxstep**3
+            print("# NORM OF THE 3D HISTO:", h3d.flatten().sum() * dqxstep**3)
 
             # Creates an interpolation function on a 3D grid
             hxyz = RegularGridInterpolator((xgrid, ygrid, zgrid), h3d)
@@ -194,7 +195,7 @@
             fth0y0 = h0y0 * 0.0
             fth00z = h00z * 0.0
 
-            print "# computing FT for block #", x
+            print("# computing FT for block #", x)
             for i in range(len(fthx00)):
                 fthx00[i] = (hx00 * np.cos(pxgrid[i] * dqxgrid)).sum() * dqxstep
                 fth0y0[i] = (h0y0 * np.cos(pygrid[i] * dqygrid)).sum() * dqystep
@@ -226,28 +227,28 @@
         np.savetxt(prefix + "_hy.data", np.c_[xgrid, h0y0])
         np.savetxt(prefix + "_hz.data", np.c_[xgrid, h00z])
 
-        print "# NORM OF npx", fthx00.sum() * pxstep
-        print "# NORM OF npx", fth0y0.sum() * pystep
-        print "# NORM OF npx", fth00z.sum() * pzstep
+        print("# NORM OF npx", fthx00.sum() * pxstep)
+        print("# NORM OF npx", fth0y0.sum() * pystep)
+        print("# NORM OF npx", fth00z.sum() * pzstep)
 
         # Calculates the average values of the second moments.
-        print RefactoringTool: Refactored ./tools/py/get_np_xyz.py
"px^2:", np.sum(np.asarray(px2list), axis=0) * norm, "+/-", np.std(np.asarray(px2list)) * np.sqrt(n_block) * norm
-        print "py^2:", np.sum(np.asarray(py2list), axis=0) * norm, "+/-", np.std(np.asarray(py2list)) * np.sqrt(n_block) * norm
-        print "pz^2:", np.sum(np.asarray(pz2list), axis=0) * norm, "+/-", np.std(np.asarray(pz2list)) * np.sqrt(n_block) * norm
-
-        print "# time taken (s)", time.clock() - start
+        print("px^2:", np.sum(np.asarray(px2list), axis=0) * norm, "+/-", np.std(np.asarray(px2list)) * np.sqrt(n_block) * norm)
+        print("py^2:", np.sum(np.asarray(py2list), axis=0) * norm, "+/-", np.std(np.asarray(py2list)) * np.sqrt(n_block) * norm)
+        print("pz^2:", np.sum(np.asarray(pz2list), axis=0) * norm, "+/-", np.std(np.asarray(pz2list)) * np.sqrt(n_block) * norm)
+
+        print("# time taken (s)", time.clock() - start)
 
     else:
         n_block = int(step / bsize)
         if (n_block == 0):
-            print 'not enough data to build a block'
+            print('not enough data to build a block')
             exit()
-        for x in xrange(n_block):
+        for x in range(n_block):
             dq = delta[x * bsize: (x + 1) * bsize]
             dfx = fx[x * bsize: (x + 1) * bsize]
             dfy = fy[x * bsize: (x + 1) * bsize]
             dfz = fz[x * bsize: (x + 1) * bsize]
-            print "# Computing 3D histogram."
+            print("# Computing 3D histogram.")
             h3d_der = histo3d_der(dq, dfx, dqxgrid, dqygrid, dqzgrid, ns, cut, np.sqrt(T * P * m), bsize, m, P, T)
             h3dlist.append(h3d_der)
 
@@ -268,12 +269,12 @@
         h = np.cumsum((hx00 - hx00[::-1]) * 0.5) * dqxstep
         h = h / (h.sum() * dqxstep) * n_block * bsize
         hx00 = hx00 / (h.sum() * dqxstep) * n_block * bsize
-        print h.sum() * dqxstep
+        print(h.sum() * dqxstep)
         np.savetxt("hx.data", np.c_[xgrid, h])
 
-        print "# px^2 (from the 2nd derivative of the histogram)", (30.0 * avghx[(ns - 1) / 2] - 16.0 * avghx[(ns - 1) / 2 + 1] - 16.0 * avghx[(ns - 1) / 2 - 1] + avghx[(ns - 1) / 2 - 2] + avghx[(ns - 1) / 2 + 2]) / dqxstep**2 / norm_npx / 12.0
-
-        print (30.0 * h[(ns - 1) / 2] - 16.0 * h[(ns - 1) / 2 + 1] - 16.0 * h[(ns - 1) / 2 - 1] + h[(ns - 1) / 2 - 2] + h[(ns - 1) / 2 + 2]) / dqxstep**2 / (n_block * bsize) / 12.0
+        print("# px^2 (from the 2nd derivative of the histogram)", (30.0 * avghx[(ns - 1) / 2] - 16.0 * avghx[(ns - 1) / 2 + 1] - 16.0 * avghx[(ns - 1) / 2 - 1] + avghx[(ns - 1) / 2 - 2] + avghx[(ns - 1) / 2 + 2]) / dqxstep**2 / norm_npx / 12.0)
+
+        print((30.0 * h[(ns - 1) / 2] - 16.0 * h[(ns - 1) / 2 + 1] - 16.0 * h[(ns - 1) / 2 - 1] + h[(ns - 1) / 2 - 2] + h[(ns - 1) / 2 + 2]) / dqxstep**2 / (n_block * bsize) / 12.0)
 
 if __name__ == '__main__':
 
--- ./tools/py/get_np_xyz.py	(original)
+++ ./tools/py/get_np_xyz.py	(refactored)
@@ -103,12 +103,12 @@
     n_block = int(step / bsize)
 
     if (n_block == 0):
-        print 'not enough data to build a block'
+        print('not enough data to build a block')
         exit()
 
     if der == False:
-        for x in xrange(n_block):
-            print "# building the histogram for block $", x + 1
+        for x in range(n_block):
+            print("# building the histogram for block $", x + 1)
             dq = bq[x * bsize: (x + 1) * bsize]
             hx = histo(dq[:, 0] - dq[:, 3 * (P - 1)], dqxgrid, kernel, np.sqrt(T * P * m))
             hx = ((hx + hx[::-1]) * 0.5)
@@ -124,7 +124,7 @@
             npx = hx * 0.0
             npy = hy * 0.0
             npz = hz * 0.0
-            print "# computing FT for block #", x + 1
+            print("# computing FT for block #", x + 1)
             for i in range(len(hx)):
                 npx[i] = (hx * np.cos(pxgrid[i] * dqxgrid)).sum() * dqxstep
                 npy[i] = (hy * np.cos(pygrid[i] * dqygrid)).sum() * dqystep
@@ -134,15 +134,15 @@
             nplisty.append(npy)
             nplistz.append(npz)
     else:
-        for x in xrange(n_block):
-            print "# building the histogram for block $", x + 1
+        for x in range(n_block):
+            print("# building the histogram for block $", x + 1)
             dq = bq[x * bsize: (x + 1) * bsize]
             df = bf[x * bsize: (x + 1) * bsize]
 
             c = 0.5 * np.asarray([-1 + 2 * float(j) / float(P - 1) for j in range(P)])
             mwp2 = m * (P * T)**2 / (P - 1)
             bp = 1.0 / (P * T)
-            for i in xrange(len(dq)):
+            for i in range(len(dq)):
                 q = dq[i]
                 f = df[i]
                 x = q[:3] - q[-3:]
@@ -154,8 +154,8 @@
                 sc = np.asarray([(s[0::3] * c).sum(), (s[1::3] * c).sum(), (s[2::3] * c).sum()])
                 sc *= -mwp2
                 fi = -bp * (fc + sc)
-                print "@*@", x, fi
-                print "@@@", np.sqrt((x * x).sum()), np.dot(x, fi)
+                print("@*@", x, fi)
+                print("@@@", np.sqrt((x * x).sum()), np.dot(x, fi))
 
             hx = histo_der(dq[:, 0::3], df[:, 0::3], dqxgrid, kernel, np.sqrt(T * P * m), m, P, T)
             hx = np.cumsum((hx - hx[::-1]) * 0.5) * dqxstep
@@ -163,7 +163,7 @@
             hy = np.cumsum((hy - hy[::-1]) * 0.5) * dqystep
             hz = histo_der(dq[:, 2::3], df[:, 2::3], dqzgrid, kernel, np.sqrt(T * P * m), m, P, T)
             hz = np.cumsum((hz - hz[::-1]) * 0.5) * dqzstep
-            print hx.sum() * dqxstep, hy.sum() * dqystep, hz.sum() * dqzstep
+            print(hx.sum() * dqxstep, hy.sum() * dqystep, hz.sum() * dqzstep)
             hxlist.append(hx)
             hylist.append(hy)
             hzlist.append(hz)
@@ -172,7 +172,7 @@
             npx = hx * 0.0
             npy = hy * 0.0
             npz = hz * 0.0
-            print "# computing FT for block #", x + 1
+            print("# computing FT for block #", x + 1)
             for i in range(len(hx)):
                 npx[i] = (hx * np.cos(pxgrid[i] * dqxgrid)).sum() * dqxstep
                 npy[i] = (hy * np.cos(pygrid[i] * dqygrid)).sum() * dqystep
@@ -194,13 +194,13 @@
     norm_npy = avghy[(ns - 1) / 2]
     norm_npz = avghz[(ns - 1) / 2]
 
-    print "# Dx^2", np.dot(dqxgrid**2, avghx) * dqxstep / (bsize * n_block)
-    print "# Dy^2", np.dot(dqygrid**2, avghy) * dqystep / (bsize * n_block)
-    print "# Dz^2", np.dot(dqzgrid**2, avghz) * dqzstep / (bsize * n_block)
-
-    print "# px^2 (from the 2nd derivative of the histogram)", (30.0 * avghx[(ns - 1) / 2] - 16.0 * avghx[(ns - 1) / 2 + 1] - 16.0 * avghx[(ns - 1) / 2 - 1] + avghx[(ns - 1) / 2 - 2] + avghx[(ns - 1) / 2 + 2]) / dqxstep**2 / norm_npx / 12.0
-    print "# py^2 (from the 2nd derivative of the histogram)", (30.0 * avghy[(ns - 1) / 2] - 16.0 * avghy[(ns - 1) / 2 + 1] - 16.0 * avghy[(ns - 1) / 2 - 1] + avghy[(ns - 1) / 2 - 2] + avghy[(ns - 1) / 2 + 2]) / dqystep**2 / norm_npy / 12.0
-    print "# pz^2 (from the 2nd derivative of the histogram)", (30.0 * avghz[(ns - 1) / 2] - 16.0 * avghz[(ns - 1) / 2 + 1] - 16.0 * avghz[(ns - 1) / 2 - 1] + avghz[(ns - 1) / 2 - 2] + avghz[(ns - 1) / 2 + 2]) / dqzstep**2 / norm_npz / 12.0
+    print("# Dx^2", np.dot(dqxgrid**2, avghx) * dqxstep / (bsize * n_block))
+    print("# Dy^2", np.dot(dqygrid**2, avghy) * dqystep / (bsize * n_block))
+    print("# Dz^2", np.dot(dqzgrid**2, avghz) * dqzstep / (bsize * n_block))
+
+    print("# px^2 (from the 2nd derivative of the histogram)", (30.0 * avghx[(ns - 1) / 2] - 16.0 * avghx[(ns - 1) / 2 + 1] - 16.0 * avghx[(ns - 1) / 2 - 1] + avghx[(ns - 1) / 2 - 2] + avghx[(ns - 1) / 2 + 2]) / dqxstep**2 / norm_npx / 12.0)
+    print("# py^2 (from the 2nd derivative of the histogram)", (30.0 * avghy[(ns - 1) / 2] - 16.0 * avghy[(ns - 1) / 2 + 1] - 16.0 * avghy[(ns - 1) / 2 - 1] + avghy[(ns - 1) / 2 - 2] + avghy[(ns - 1) / 2 + 2]) / dqystep**2 / norm_npy / 12.0)
+    print("# pz^2 (from the 2nd derivative of the histogram)", (30.0 * avghz[(ns - 1) / 2] - 16.0 * avghz[(ns - 1) / 2 + 1] - 16.0 * avghz[(ns - 1) / 2 - 1] + avghz[(ns - 1) / 2 - 2] + avghz[(ns - 1) / 2 + 2]) / dqzstep**2 / norm_npz / 12.0)
 
     np.savetxt("hxx.data", avghx)
     np.savetxt(str(prefix + "histo.dRefactoringTool: Refactored ./tools/py/getacf.py
RefactoringTool: Refactored ./tools/py/getproperty.py
RefactoringTool: Refactored ./tools/py/gleacf.py
ata"), np.c_[dqxgrid, avghx / (bsize * n_block), errhx, dqygrid, avghy / (bsize * n_block), errhy, dqzgrid, avghz / (bsize * n_block), errhz])
@@ -225,13 +225,13 @@
         py2.append(np.dot(pygrid**2, np.asarray(nplisty)[i, :]) * pystep / norm_npy / 2.0 / np.pi * n_block)
         pz2.append(np.dot(pzgrid**2, np.asarray(nplistz)[i, :]) * pzstep / norm_npz / 2.0 / np.pi * n_block)
 
-    print '# number of blocks', n_block
-    print '# px^2', np.mean(np.asarray(px2)), 'sigmax', np.std(np.asarray(px2)) / np.sqrt(n_block)
-    print '# py^2', np.mean(np.asarray(py2)), 'sigmay', np.std(np.asarray(py2)) / np.sqrt(n_block)
-    print '# pz^2', np.mean(np.asarray(pz2)), 'sigmaz', np.std(np.asarray(pz2)) / np.sqrt(n_block)
-    print '# p^2', np.mean(np.asarray(px2) + np.asarray(py2) + np.asarray(pz2)), 'sigma', np.std(np.asarray(px2) + np.asarray(py2) + np.asarray(pz2)) / np.sqrt(n_block)
-
-    print "# time taken (s)", time.clock() - start
+    print('# number of blocks', n_block)
+    print('# px^2', np.mean(np.asarray(px2)), 'sigmax', np.std(np.asarray(px2)) / np.sqrt(n_block))
+    print('# py^2', np.mean(np.asarray(py2)), 'sigmay', np.std(np.asarray(py2)) / np.sqrt(n_block))
+    print('# pz^2', np.mean(np.asarray(pz2)), 'sigmaz', np.std(np.asarray(pz2)) / np.sqrt(n_block))
+    print('# p^2', np.mean(np.asarray(px2) + np.asarray(py2) + np.asarray(pz2)), 'sigma', np.std(np.asarray(px2) + np.asarray(py2) + np.asarray(pz2)) / np.sqrt(n_block))
+
+    print("# time taken (s)", time.clock() - start)
 
 if __name__ == '__main__':
     parser = argparse.ArgumentParser(description=None)
--- ./tools/py/getacf.py	(original)
+++ ./tools/py/getacf.py	(refactored)
@@ -62,8 +62,8 @@
     nblocks = 0
     dt = unit_to_internal("time", timestep[1], float(timestep[0]))
     data = np.zeros((bsize, nspecies, 3), float)
-    time = np.asarray(range(mlag + 1)) * dt
-    omega = np.asarray(range(2 * (mlag + npad))) / float(2 * (mlag + npad)) * (2 * np.pi / dt)
+    time = np.asarray(list(range(mlag + 1))) * dt
+    omega = np.asarray(list(range(2 * (mlag + npad)))) / float(2 * (mlag + npad)) * (2 * np.pi / dt)
     fvvacf = omega.copy() * 0.0
     fvvacf2 = fvvacf.copy() * 0.0
     vvacf = time.copy() * 0.0
@@ -83,7 +83,7 @@
 
     ff = open(ifile)
     # Skips the first fskip frames
-    for x in xrange(fskip):
+    for x in range(fskip):
         rr = read_file_raw("xyz", ff)
 
     while True:
--- ./tools/py/getproperty.py	(original)
+++ ./tools/py/getproperty.py	(refactored)
@@ -47,7 +47,7 @@
                 warning("Could not find " + propertyname + " in file " + inputfile)
                 raise EOFError
             line = line.split()
-            if (step >= skip): print line[icol]
+            if (step >= skip): print(line[icol])
             step += 1
         except EOFError:
             break
--- ./tools/py/gleacf.py	(original)
+++ ./tools/py/gleacf.py	(refactored)
@@ -77,8 +77,8 @@
     a, O = np.linalg.eig(dAqp)
     O1 = np.linalg.inv(O)
     W = np.dot(np.dot(O1, dDqp), O1.T)
-    for i in xrange(len(W)):
-        for j in xrange(len(W)):
+    for i in range(len(W)):
+        for j in range(len(W)):
             W[i, j] /= (a[i] + a[j])
     nC = np.real(np.dot(O, np.dot(W, O.T)))
 
@@ -96,7 +96,7 @@
     om2list = omlist**2
     y = 0
     if Ap[0, 0] < 2.0 * dw:
-        print "# WARNING: White-noise term is weaker than the spacing of the frequency grid. Will increase automatically to avoid instabilities in the numerical integration."
+        print("# WARNING: White-noise term is weaker than the spacing of the frequency grid. Will increase automatically to avoid instabilities in the numerical integration.")
 
     # outer loop over the physical frequency
     for omega_0 in omlist:
@@ -198,15 +198,15 @@
     iy = ifacf[:, 1]
 
     # computes the facf kernel
-    print "# computing the kernel."
+    print("# computing the kernel.")
     ker = gleKernel(ix, Ap, Dp)
 
     # (de-)convolutes the spectrum
     if(action == "conv"):
-        print "# printing the output spectrum."
+        print("# printing the output speRefactoringTool: Refactored ./tools/py/heat_capacity_ppi.py
RefactoringTool: No changes to ./tools/py/kinetic2tag.py
RefactoringTool: Refactored ./tools/py/kinetic_energy_ppi.py
RefactoringTool: No changes to ./tools/py/mergebeadspdb.py
RefactoringTool: Refactored ./tools/py/mux-positions.py
ctrum.")
         output_facf((ix, np.dot(iy, ker.T)), oprefix, input_facf(path2ifacf, nrows, 1))
     elif(action == "deconv"):
-        print "# deconvoluting the input spectrum."
+        print("# deconvoluting the input spectrum.")
         oy = ISRA(ix, ker, iy, dparam, oprefix)
 
 if __name__ == '__main__':
--- ./tools/py/heat_capacity_ppi.py	(original)
+++ ./tools/py/heat_capacity_ppi.py	(refactored)
@@ -65,24 +65,24 @@
     # check that we found the same number of positions and forces files
     nbeads = len(fns_pos)
     if nbeads != len(fns_for):
-        print fns_pos
-        print fns_for
+        print(fns_pos)
+        print(fns_for)
         raise ValueError("Mismatch between number of input files for forces and positions.")
 
     # print some information
-    print 'temperature = {:f} K'.format(float(temp))
-    print
-    print 'number of beads = {:d}'.format(nbeads)
-    print
-    print 'positions and forces file names:'
+    print('temperature = {:f} K'.format(float(temp)))
+    print()
+    print('number of beads = {:d}'.format(nbeads))
+    print()
+    print('positions and forces file names:')
     for fn_pos, fn_for in zip(fns_pos, fns_for):
-        print '{:s}   {:s}'.format(fn_pos, fn_for)
-    print
-    print 'potential energy file: {:s}'.format(fns_iU)
-    print
-    print 'output file name:'
-    print fn_out_en
-    print
+        print('{:s}   {:s}'.format(fn_pos, fn_for))
+    print()
+    print('potential energy file: {:s}'.format(fns_iU))
+    print()
+    print('output file name:')
+    print(fn_out_en)
+    print()
 
     # open input and output files
     ipos = [open(fn, "r") for fn in fns_pos]
@@ -112,7 +112,7 @@
     while True:  # Reading input files and calculating PPI correction
 
         if ifr % 100 == 0:
-            print '\rProcessing frame {:d}'.format(ifr),
+            print('\rProcessing frame {:d}'.format(ifr), end=' ')
             sys.stdout.flush()
 
         try:
--- ./tools/py/kinetic_energy_ppi.py	(original)
+++ ./tools/py/kinetic_energy_ppi.py	(refactored)
@@ -66,24 +66,24 @@
     # check that we found the same number of positions and forces files
     nbeads = len(fns_pos)
     if nbeads != len(fns_for):
-        print fns_pos
-        print fns_for
+        print(fns_pos)
+        print(fns_for)
         raise ValueError("Mismatch between number of input files for forces and positions.")
 
     # print some information
-    print 'temperature = {:f} K'.format(float(temp))
-    print
-    print 'number of beads = {:d}'.format(nbeads)
-    print
-    print 'positions and forces file names:'
+    print('temperature = {:f} K'.format(float(temp)))
+    print()
+    print('number of beads = {:d}'.format(nbeads))
+    print()
+    print('positions and forces file names:')
     for fn_pos, fn_for in zip(fns_pos, fns_for):
-        print '{:s}   {:s}'.format(fn_pos, fn_for)
-    print
-    print 'potential energy file: {:s}'.format(fns_iU)
-    print
-    print 'output file name:'
-    print fn_out_en
-    print
+        print('{:s}   {:s}'.format(fn_pos, fn_for))
+    print()
+    print('potential energy file: {:s}'.format(fns_iU))
+    print()
+    print('output file name:')
+    print(fn_out_en)
+    print()
 
     # open input and output files
     ipos = [open(fn, "r") for fn in fns_pos]
@@ -116,7 +116,7 @@
     while True:  # Reading input files and calculating PPI correction
 
         if ifr % 100 == 0:
-            print '\rProcessing frame {:d}'.format(ifr),
+            print('\rProcessing frame {:d}'.format(ifr), end=' ')
             sys.stdout.flush()
 
         try:
--- ./tools/py/mux-positions.py	(original)
+++ ./tools/py/mux-positions.py	(refactored)
@@ -1,6 +1,6 @@
 #!/usr/bin/env python2
 
-from __future__ import print_function
+
 
 import os
 import sys
@@ -63,7 +63,7 @@
         try:
             # Get the frames from all trajectories...
             for idx, trj in enumerate(trjs_in):
-                frame = trj.next()
+                frame = next(trj)
 
                 # gets units from first frame
                 dimension, units, cellRefactoringTool: Refactored ./tools/py/paraweights.py
RefactoringTool: Refactored ./tools/py/pepper.py
RefactoringTool: Refactored ./tools/py/planetary.py
_units = auto_units(comment=frame["comment"])
--- ./tools/py/paraweights.py	(original)
+++ ./tools/py/paraweights.py	(refactored)
@@ -56,7 +56,7 @@
     isimul.parse(xmlrestart.fields[0][1])
     verbosity.level = "quiet"
     banner()
-    print "# Printing out temperature re-weighing factors for a parallel tempering simulation"
+    print("# Printing out temperature re-weighing factors for a parallel tempering simulation")
     simul = isimul.fetch()
 
     if simul.mode != "paratemp":
--- ./tools/py/pepper.py	(original)
+++ ./tools/py/pepper.py	(refactored)
@@ -102,28 +102,28 @@
                 python_files = python_files + [filename]
             else:
                 if verbosity != 'silent':
-                    print filename, 'is not a python file, skipping '
+                    print(filename, 'is not a python file, skipping ')
         # execute autopep8 only on python files, without being recursive
         autopep8_args = autopep8_args + python_files
         if not python_files:
-            print 'No python files to process.'
+            print('No python files to process.')
             sys.exit()
         if verbosity != 'silent':
-            print 'Running autopep8 on the files: ', ' '.join(python_files)
+            print('Running autopep8 on the files: ', ' '.join(python_files))
     else:
         # perform recursive search in the given directory
         if os.path.isdir(path):
             os.chdir(path)
             autopep8_args = autopep8_args + ['-r'] + ['.']
             if verbosity != 'silent':
-                print 'Running autopep8 recursively in the directory: ', os.getcwd()
+                print('Running autopep8 recursively in the directory: ', os.getcwd())
         else:
-            print 'The given directory does not exist: ', path
+            print('The given directory does not exist: ', path)
             sys.exit()
 
     if is_in_check_mode:
         # perform check, do not change files
-        print 'Performing check of PEP-8 compliance'
+        print('Performing check of PEP-8 compliance')
         style_check = pycodestyle.StyleGuide(select=styles_to_be_corrected)
         if files is not None:
             for filename in files:
@@ -145,19 +145,19 @@
             # We must strip, otherwise we get double newline
             line = process.stdout.readline().rstrip()
             if re.match('\[Errno.*\]', line):
-                print line
+                print(line)
             else:
                 if verbosity == 'high':
-                    print line
+                    print(line)
                 elif verbosity == 'medium':
                     if re.match('\[file:.*\]', line):
                         # new file is processed. If something left unfixed, print message now
                         if issue_line_buffer is not '':
-                            print issue_line_buffer
+                            print(issue_line_buffer)
                         # Pattern: [file:filename]
                         filename_line = re.search(r'\[file:(\S+)\]', line)
                         # Print only filename
-                        print filename_line.group(1)
+                        print(filename_line.group(1))
                     elif re.match('.*0 issue.*', line):
                         # there are 0 issues left, clean the buffer
                         issue_line_buffer = ''
@@ -181,6 +181,6 @@
         if verbosity == 'medium':
             # end of processing. If something left unfixed, print the message now
             if issue_line_buffer is not '':
-                print issue_line_buffer
+                print(issue_line_buffer)
     if verbosity != 'silent':
-        print 'autopep8 terminated'
+        print('autopep8 terminated')
--- ./tools/py/planetary.py	(original)
+++ ./tools/py/planetary.py	(refactored)
@@ -1,6 +1,6 @@
 #!/usr/bin/env python2
 
-from __future__ import division
+
 
 import os
 import sys
@@ -101,7 +101,7 @@
 
     with open("{}.xc.xyz".format(args.prefix), "r") as f:
         args.natoms = int(f.readline())
-        for i in xrange(args.natoms + 2):
+        for i in range(args.natoms + 2):
             f.readline()
         stride2 = int(f.readline().split("Step:")[1].split("Bead:")[0])
 
@@ -120,12 +120,12 @@
         estimators.append(mod)
 
     args.estmods = ", ".join(args.estmods)
-    print "----------------------"
-    print "DETERMINED INPUT PARAMETERS"
-    print "----------------------"
+    print("----------------------")
+    print("DETERMINED INPUT PARAMETERS")
+    print("----------------------")
     for k, v in sorted(args.__dict__.items()):
-        print "{:<30s}{}".format(k, v)
-    print "----------------------"
+        print("{:<30s}{}".format(k, v))
+    print("----------------------")
     sys.stdout.flush()
 
     return args.prefix, args.natoms, args.iseed, args.npl, args.npts, args.stride, args.dt, args.temperature, estimators
@@ -149,8 +149,8 @@
     npts, npl, ndim = A.shape
     tcf = np.zeros(npts)
     weights = np.arange(npts, 0, -1, dtype=float)
-    for i in xrange(npl):
-        for j in xrange(ndim):
+    for i in range(npl):
+        for j in range(ndim):
             tcf[:] += np.correlate(A[:, i, j], B[:, i, j], mode="full")[npts - 1::-1] / weights
     # for i in xrange(npl):
     #    for j in xrange(ndim):
@@ -343,7 +343,7 @@
         self.read_qcpc()
         self.matrix_setup()
 
-        for j in xrange(self.stride):
+        for j in range(self.stride):
             # Linear interpolation of frequency matrix
             self.omega_interp[:] = (j * self.omega + (self.stride - j) * self.omega_old) / self.stride
             # Velocity verlet integrator
@@ -385,7 +385,7 @@
                 # ( npl,    prod(mod.Ashape),    3*self.natoms)
                 A1q = mod.Afunc1(self.qsum)
 
-                for j in xrange(self.npl):
+                for j in range(self.npl):
                     mod.Atemp[:] = 0.0
 
                     for k, val in enumerate(A1q[j]):
@@ -400,7 +400,7 @@
                 # ( npl,    prod(mod.Ashape),    3*self.natoms,    3*self.natoms )
                 A2q = mod.Afunc2(self.qsum)
 
-                for j in xrange(self.npl):
+                for j in range(self.npl):
                     mod.Atemp[:] = 0.0
 
                     for k, val in enumerate(A2q[j]):
@@ -428,9 +428,9 @@
         for mod in self.estimators:
             if hasattr(mod, "Afunc0"):
                 np.savetxt("{}_ce_{}.dat".format(self.prefix, mod.name), np.transpose([tarr, mod.TCF_c]))
-                print "Saved {} (centroid) to {}_ce_{}.dat".format(mod.name, self.prefix, mod.name)
+                print("Saved {} (centroid) to {}_ce_{}.dat".format(mod.name, self.prefix, mod.name))
             np.savetxt("{}_pl_{}.dat".format(self.prefix, mod.name), np.transpose([tarr, mod.TCF]))
-            print "Saved {} (full) to {}_pl_{}.dat".format(mod.name, self.prefix, mod.name)
+            print("Saved {} (full) to {}_pl_{}.dat".format(mod.name, self.prefix, mod.name))
 
     def shutdown(self):
         """
@@ -466,7 +466,7 @@
                        by calling self.shutdown()
 
         """
-        print "STARTING PLANETARY SIMULATION"
+        print("STARTING PLANETARY SIMULATION")
         self.read_omega2()
         self.read_qcpc()
         self.get_masses()
@@ -478,20 +478,20 @@
         self.psum[:] = self.p + self.pc[:, np.newaxis]
         self.estimate(0)
         self.count = 0
-        for i in xrange(self.npts - 1):
+        for i in range(self.npts - 1):
             self.step()
             self.count += 1  # this might be pointless, obviously self.count = i
             self.estimate(i + 1)
-            print "Completed step {:d} of {:d}".format(self.count, self.npts - 1)
+            print("Completed step {:d} of {:d}".format(self.count, self.npts - 1))
             sys.stdout.flush()
         if correlate:
-            print "CALCULATING TCFS"
+            print("CALCULATING TCFS")
             self.correlate()
         if write:
-            print "SAVING TCFS"
+            print("SAVING TCFS")
             self.write_tcfs()
         if shutdown:
-            print "SIRefactoringTool: Refactored ./tools/py/posforce2kinetic.py
RefactoringTool: Refactored ./tools/py/potential_energy_ppi.py
RefactoringTool: Refactored ./tools/py/rdf_ppi.py
MULATION COMPLETE. SHUTTING DOWN."
+            print("SIMULATION COMPLETE. SHUTTING DOWN.")
             self.shutdown()
         else:
             warning("Reached end of simulation but files may remain open.")
--- ./tools/py/posforce2kinetic.py	(original)
+++ ./tools/py/posforce2kinetic.py	(refactored)
@@ -37,23 +37,23 @@
     # check that we found the same number of positions and forces files
     nbeads = len(fns_pos)
     if nbeads != len(fns_for):
-        print fns_pos
-        print fns_for
+        print(fns_pos)
+        print(fns_for)
         raise ValueError("Mismatch between number of input files for forces and positions.")
 
     # print some information
-    print 'temperature = {:f} K'.format(T)
-    print
-    print 'number of beads = {:d}'.format(nbeads)
-    print
-    print 'positions and forces file names:'
+    print('temperature = {:f} K'.format(T))
+    print()
+    print('number of beads = {:d}'.format(nbeads))
+    print()
+    print('positions and forces file names:')
     for fn_pos, fn_for in zip(fns_pos, fns_for):
-        print '{:s}   {:s}'.format(fn_pos, fn_for)
-    print
-    print 'output file names:'
-    print fn_out_kin
-    print fn_out_kod
-    print
+        print('{:s}   {:s}'.format(fn_pos, fn_for))
+    print()
+    print('output file names:')
+    print(fn_out_kin)
+    print(fn_out_kod)
+    print()
 
     temp = unit_to_internal("energy", "kelvin", T)
 
@@ -69,7 +69,7 @@
 
         # print progress
         if ifr % 100 == 0:
-            print '\rProcessing frame {:d}'.format(ifr),
+            print('\rProcessing frame {:d}'.format(ifr), end=' ')
             sys.stdout.flush()
 
         # load one frame
@@ -116,7 +116,7 @@
 
         ifr += 1
 
-    print '\rProcessed {:d} frames.'.format(ifr)
+    print('\rProcessed {:d} frames.'.format(ifr))
 
     ikin.close()
     ikod.close()
--- ./tools/py/potential_energy_ppi.py	(original)
+++ ./tools/py/potential_energy_ppi.py	(refactored)
@@ -63,19 +63,19 @@
     nbeads = len(fns_for)
 
     # print some information
-    print 'temperature = {:f} K'.format(float(temp))
-    print
-    print 'number of beads = {:d}'.format(nbeads)
-    print
-    print 'forces file names:'
+    print('temperature = {:f} K'.format(float(temp)))
+    print()
+    print('number of beads = {:d}'.format(nbeads))
+    print()
+    print('forces file names:')
     for fn_for in fns_for:
-        print '{:s}'.format(fn_for)
-    print
-    print 'potential energy file: {:s}'.format(fns_iU)
-    print
-    print 'output file name:'
-    print fn_out_en
-    print
+        print('{:s}'.format(fn_for))
+    print()
+    print('potential energy file: {:s}'.format(fns_iU))
+    print()
+    print('output file name:')
+    print(fn_out_en)
+    print()
 
     # open input and output files
     ifor = [open(fn, "r") for fn in fns_for]
@@ -103,7 +103,7 @@
     while True:  # Reading input files and calculating PPI correction
 
         if ifr % 100 == 0:
-            print '\rProcessing frame {:d}'.format(ifr),
+            print('\rProcessing frame {:d}'.format(ifr), end=' ')
             sys.stdout.flush()
 
         try:
--- ./tools/py/rdf_ppi.py	(original)
+++ ./tools/py/rdf_ppi.py	(refactored)
@@ -56,23 +56,23 @@
     # check that we found the same number of positions and forces files
     nbeads = len(fns_pos)
     if nbeads != len(fns_for):
-        print fns_pos
-        print fns_for
+        print(fns_pos)
+        print(fns_for)
         raise ValueError("Mismatch between number of input files for forces and positions.")
 
     # print some information
-    print 'temperature = {:f} K'.format(float(temp))
-    print
-    print 'number of beads = {:d}'.format(nbeads)
-    print
-    print 'positions and forces file names:'
+    print('temperature = {:f} K'.format(float(temp)))
+    print()
+    print('number of beads = {:d}'.format(nbeads))
+    print()
+    print('positions and forces file names:')
     for fn_pos, fn_for in zip(fns_pos, fns_for):
-        print '{:s}   {:s}'.format(fn_pos, fn_for)
-    print
-    print 'output file names:'
-    pRefactoringTool: Refactored ./tools/py/regtest-parallel.py
RefactoringTool: Refactored ./tools/py/regtest.py
rint fn_out_rdf
-    print fn_out_rdf_q
-    print
+        print('{:s}   {:s}'.format(fn_pos, fn_for))
+    print()
+    print('output file names:')
+    print(fn_out_rdf)
+    print(fn_out_rdf_q)
+    print()
 
     # open input and output files
     ipos = [open(fn, "r") for fn in fns_pos]
@@ -109,7 +109,7 @@
     while noteof:  # Reading input files and calculating PPI correction
 
         if ifr % 100 == 0:
-            print '\rProcessing frame {:d}'.format(ifr),
+            print('\rProcessing frame {:d}'.format(ifr), end=' ')
             sys.stdout.flush()
 
         try:
--- ./tools/py/regtest-parallel.py	(original)
+++ ./tools/py/regtest-parallel.py	(refactored)
@@ -91,7 +91,7 @@
 
 import argparse
 import os
-import Queue
+import queue
 import re
 import shutil
 import shlex
@@ -144,8 +144,8 @@
 except KeyError:
     IPI_ROOT_FOLDER = os.path.abspath('../')
 
-QUEUE_ALL = Queue.Queue()  # Queue to access the "run"
-QUEUE_COM = Queue.Queue()  # Queue to access the "compare"
+QUEUE_ALL = queue.Queue()  # Queue to access the "run"
+QUEUE_COM = queue.Queue()  # Queue to access the "compare"
 
 
 def main():
@@ -174,7 +174,7 @@
         QUEUE_ALL.put(test_obj)
         index += 10
 
-    print 'Starting tests'
+    print('Starting tests')
     running_test = []
     running_com = []
     if int(_parser()['nproc']) > 1:
@@ -185,7 +185,7 @@
             if len(running_test) < _parser()['nproc']:
                 try:
                     running_test.append(QUEUE_ALL.get_nowait())
-                except Queue.Empty:
+                except queue.Empty:
                     pass
                 else:
                     running_test[-1].generate_output = True
@@ -203,7 +203,7 @@
             if len(running_com) < 1:
                 try:
                     running_com.append(QUEUE_COM.get_nowait())
-                except Queue.Empty:
+                except queue.Empty:
                     pass
                 else:
                     running_com[-1].copy_reference = _parser()['create_reference']
@@ -225,7 +225,7 @@
         for _thr in running_test:
             _thr.die = True
 
-    print
+    print()
 
 
 def _parser():
@@ -343,9 +343,9 @@
                 msg += ' > ' + str(os.path.split(root)[1]) + '\n'
 
     if len(test_list) < 1:
-        print "**No test found!**"
+        print("**No test found!**")
     else:
-        print msg
+        print(msg)
 
     return test_list
 
@@ -359,7 +359,7 @@
 
     with open(path_to_test) as _file:
         _text = _file.read()
-    print _text[:100]
+    print(_text[:100])
     return len([x.group(1) for x in REGTEST_STRING_RGX.finditer(_text)]) > 0
 
 
@@ -654,7 +654,7 @@
                         remove_file(straj['old_filename'])
                         shutil.copy2(straj['new_filename'],
                                      straj['old_filename'])
-            except IOError, e:
+            except IOError as e:
                 self.test_status = 'ERROR'
                 self.msg += 'Error while copying the new reference!!\n'
                 self.msg += "Unable to copy file. %s" % e
@@ -1000,7 +1000,7 @@
     while answer.lower() not in _yes and answer not in _no:
         sys.stderr.write(msg)
 
-        answer = raw_input()
+        answer = input()
 
         if answer.lower() in _yes:
             return True
--- ./tools/py/regtest.py	(original)
+++ ./tools/py/regtest.py	(refactored)
@@ -72,14 +72,14 @@
     try:
         os.makedirs(run_dir)
     except OSError:
-        print 'The directory %s exists. Do you want to delete its contents and continue? (y/n)' % run_dir
+        print('The directory %s exists. Do you want to delete its contents and continue? (y/n)' % run_dir)
         if answer_is_y():
             shutil.rmtree(run_dir)
             os.makedirs(run_dir)
         else:
             quit('Terminated')
     if is_in_reference_mode:
-        print 'Do you agree to replace references if they exist? (y/n)'
+        print('Do you agree to replace references if they exist? (y/n)')
         if answer_is_y():
             replace_references = True
         else:
@@ -90,20 +90,20 @@
     for test_candidate in list_of_test_candidates:
         try:
             list_of_test_cases.append(TestCase(test_candidate))
-        except TypeError, e:
-            print '--> Could not create test case:', test_candidate.name
-            print str(e)
+        except TypeError as e:
+            print('--> Could not create test case:', test_candidate.name)
+            print(str(e))
             continue
     queue_of_test_instances = deque([])
     counter = Counter()
-    print
-    print 'List of test cases to be executed:'
+    print()
+    print('List of test cases to be executed:')
     for test_case in list_of_test_cases:
         test_instance_dir = os.path.join(run_dir, test_case.name)
         queue_of_test_instances.append(TestInstance(test_case, test_instance_dir, counter))
-        print '-->', test_case.name
-    print
-    print 'Executing following test instances:'
+        print('-->', test_case.name)
+    print()
+    print('Executing following test instances:')
     while queue_of_test_instances:
         test_instance = queue_of_test_instances.popleft()
         try:
@@ -118,13 +118,13 @@
                         if not results.files_are_equal():
                             test_passed = False
                     if test_passed:
-                        print '    PASSED'
+                        print('    PASSED')
                     else:
-                        print '    FAILED'
+                        print('    FAILED')
                         for results in differences:
                             results.print_differences()
                 except ValueError as e:
-                    print >> sys.stderr, str(e)
+                    print(str(e), file=sys.stderr)
             else:
                 try:
                     test_instance.put_output_as_reference()
@@ -132,13 +132,13 @@
                     if replace_references:
                         test_instance.test_case.remove_reference()
                         test_instance.put_output_as_reference()
-                        print '    SUCCESS: References replaced'
+                        print('    SUCCESS: References replaced')
                     else:
-                        print '    SKIPPING: References not copied'
+                        print('    SKIPPING: References not copied')
 
         except (RegtestTimeout, WrongDriverCommand, IPIError, OutputCorrupted) as e:
-            print '    ERROR'
-            print >> sys.stderr, str(e), 'Test instance run directory:', test_instance.run_dir
+            print('    ERROR')
+            print(str(e), 'Test instance run directory:', test_instance.run_dir, file=sys.stderr)
 
             continue
         except WrongIPICommand as e:
@@ -402,10 +402,10 @@
                     ipi_out.write(stderr)
                 raise IPIError('I-PI Error\n' + stderr)
         except KeyboardInterrupt as e:
-            print '\nKeyboard interrupt!'
+            print('\nKeyboard interrupt!')
             traceback.print_exc(file=sys.stdout)
             if ipi_proc.poll() is None:
-                print 'Trying to shut down i-PI cleanly...'
+                print('Trying to shut down i-PI cleanly...')
                 # Let I-PI start before we terminate it, otherwise it can hang up
                 time.sleep(Parameters.ipi_starting_time)
                 ipi_proc.terminate()
@@ -413,7 +413,7 @@
                 if stderr:
                     with open(ipi_output_path, 'a') as ipi_out:
                         ipi_out.write(stderr)
-                print 'i-PI terminated'
+                print('i-PI terminated')
             for prc in driver_prcs:
                 if prc.poll() is None:
                     prc.terminate()
@@ -546,9 +546,9 @@
         """
         Prints differences between file1 and file2
         """
-        print 'Files:', self.file1, self.file2
+        print('Files:', self.file1, self.file2)
         for element in self.list_of_differences:
-            print 'Difference in line', element[0], 'in word', elRefactoringTool: Refactored ./tools/py/remdsort.py
RefactoringTool: Refactored ./tools/py/total_energy_ppi.py
RefactoringTool: No changes to ./tools/py/trimsim.py
RefactoringTool: No changes to ./tools/py/xyz2bin.py
RefactoringTool: No changes to ./tools/py/xyz2pdb.py
RefactoringTool: Files that need to be modified:
RefactoringTool: ./setup.py
RefactoringTool: ./doc/create_man.py
RefactoringTool: ./doc/help.py
RefactoringTool: ./doc/help_list.py
RefactoringTool: ./examples/ASE/W2-PIMD-MP2/run-ase.py
RefactoringTool: ./examples/cp2k/nvt-cl/idtau_plot.py
RefactoringTool: ./examples/hswfqmc/FixWF_1Bead/intau_plot.py
RefactoringTool: ./examples/hswfqmc/WFOpt_2Bead/intau_plot.py
RefactoringTool: ./examples/lammps/h2o-planetary-64/estmod_dipole.py
RefactoringTool: ./examples/lammps/isof-vapor/process_out.py
RefactoringTool: ./examples/lammps/isof-water/process_out.py
RefactoringTool: ./examples/lammps/isofsc-vapor/process_dif.py
RefactoringTool: ./examples/lammps/isofsc-vapor/process_out.py
RefactoringTool: ./examples/lammps/isofsc-water/process_out.py
RefactoringTool: ./examples/yaff/mil53_ffsocket/run.py
RefactoringTool: ./examples/yaff/mil53_ffsocket/yaffdriver.py
RefactoringTool: ./ipi/__init__.py
RefactoringTool: ./ipi/engine/__init__.py
RefactoringTool: ./ipi/engine/atoms.py
RefactoringTool: ./ipi/engine/barostats.py
RefactoringTool: ./ipi/engine/beads.py
RefactoringTool: ./ipi/engine/cell.py
RefactoringTool: ./ipi/engine/ensembles.py
RefactoringTool: ./ipi/engine/forcefields.py
RefactoringTool: ./ipi/engine/forces.py
RefactoringTool: ./ipi/engine/initializer.py
RefactoringTool: ./ipi/engine/normalmodes.py
RefactoringTool: ./ipi/engine/outputs.py
RefactoringTool: ./ipi/engine/properties.py
RefactoringTool: ./ipi/engine/simulation.py
RefactoringTool: ./ipi/engine/system.py
RefactoringTool: ./ipi/engine/thermostats.py
RefactoringTool: ./ipi/engine/motion/__init__.py
RefactoringTool: ./ipi/engine/motion/al6xxx_kmc.py
RefactoringTool: ./ipi/engine/motion/alchemy.py
RefactoringTool: ./ipi/engine/motion/atomswap.py
RefactoringTool: ./ipi/engine/motion/constrained_dynamics.py
RefactoringTool: ./ipi/engine/motion/dynamics.py
RefactoringTool: ./ipi/engine/motion/geop.py
RefactoringTool: ./ipi/engine/motion/instanton.py
RefactoringTool: ./ipi/engine/motion/motion.py
RefactoringTool: ./ipi/engine/motion/multi.py
RefactoringTool: ./ipi/engine/motion/neb.py
RefactoringTool: ./ipi/engine/motion/phonons.py
RefactoringTool: ./ipi/engine/motion/planetary.py
RefactoringTool: ./ipi/engine/motion/ramp.py
RefactoringTool: ./ipi/engine/motion/replay.py
RefactoringTool: ./ipi/engine/smotion/__init__.py
RefactoringTool: ./ipi/engine/smotion/metad.py
RefactoringTool: ./ipi/engine/smotion/multi.py
RefactoringTool: ./ipi/engine/smotion/remd.py
RefactoringTool: ./ipi/engine/smotion/smotion.py
RefactoringTool: ./ipi/external/importlib/__init__.py
RefactoringTool: ./ipi/external/importlib/bundledimportlib.py
RefactoringTool: ./ipi/inputs/__init__.py
RefactoringTool: ./ipi/inputs/atoms.py
RefactoringTool: ./ipi/inputs/barostats.py
RefactoringTool: ./ipi/inputs/beads.py
RefactoringTool: ./ipi/inputs/cell.py
RefactoringTool: ./ipi/inputs/ensembles.py
RefactoringTool: ./ipi/inputs/forcefields.py
RefactoringTool: ./ipi/inputs/forces.py
RefactoringTool: ./ipi/inputs/initializer.py
RefactoringTool: ./ipi/inputs/interface.py
RefactoringTool: ./ipi/inputs/normalmodes.py
RefactoringTool: ./ipi/inputs/outputs.py
RefactoringTool: ./ipi/inputs/prng.py
RefactoringTool: ./ipi/inputs/simulation.py
RefactoringTool: ./ipi/inputs/system.py
RefactoringTool: ./ipi/inputs/thermostats.py
RefactoringTool: ./ipi/inputs/motion/__init__.py
RefactoringTool: ./ipi/inputs/motion/al6xxx_kmc.py
RefactoringTool: ./ipi/inputs/motion/alchemy.py
RefactoringTool: ./ipi/inputs/motion/atomswap.py
RefactoringTool: ./ipi/inputs/motion/constrained_dynamics.py
RefactoringTool: ./ipi/inputs/motion/dynamics.py
RefactoringTool: ./ipi/inputs/motion/geop.py
RefactoringTool: ./ipi/inputs/motion/instanton.py
RefactoringTool: ./ipi/inputs/motion/motion.py
RefactoringTool: ./ipi/inputs/motion/neb.py
RefactoringTool: ./ipi/inputs/motion/phonons.py
RefactoringTool: ./ipi/inputs/motion/planetary.py
RefactoringTool: ./ipi/inputs/motion/ramp.py
RefactoringTool: ./ipi/inputs/smotion/__init__.py
RefactoringTool: ./ipi/inputs/smotion/metad.py
RefactoringTool: ./ipi/inputs/smotion/remd.py
RefactoringTool: ./ipi/inputs/smotion/smotion.py
RefactoringTool: ./ipi/interfaces/__init__.py
RefactoringTool: ./ipi/interfaces/clients.py
RefactoringTool: ./ipi/interfaces/sockets.py
RefactoringTool: ./ipi/tests/__init__.py
RefactoringTool: ./ipi/tests/common.py
RefactoringTool: ./ipi/tests/test_atoms.py
RefactoringTool: ./ipi/tests/test_contraction.py
RefactoringTool: ./ipi/tests/test_depend.py
RefactoringTool: ./ipi/tests/test_interface.py
RefactoringTool: ./ipi/tests/test_io.py
RefactoringTool: ./ipi/tests/test_runs.py
RefactoringTool: ./ipi/tests/test_units.py
RefactoringTool: ./ipi/utils/__init__.py
RefactoringTool: ./ipi/utils/constrtools.py
RefactoringTool: ./ipi/utils/decorators.py
RefactoringTool: ./ipi/utils/depend.py
RefactoringTool: ./ipi/utils/exchange.py
RefactoringTool: ./ipi/utils/hesstools.py
RefactoringTool: ./ipi/utils/inputvalue.py
RefactoringTool: ./ipi/utils/instools.py
RefactoringTool: ./ipi/utils/mathtools.py
RefactoringTool: ./ipi/utils/messages.py
RefactoringTool: ./ipi/utils/mintools.py
RefactoringTool: ./ipi/utils/nmtransform.py
RefactoringTool: ./ipi/utils/prng.py
RefactoringTool: ./ipi/utils/softexit.py
RefactoringTool: ./ipi/utils/sparse.py
RefactoringTool: ./ipi/utils/units.py
RefactoringTool: ./ipi/utils/io/__init__.py
RefactoringTool: ./ipi/utils/io/io_units.py
RefactoringTool: ./ipi/utils/io/backends/__init__.py
RefactoringTool: ./ipi/utils/io/backends/io_binary.py
RefactoringTool: ./ipi/utils/io/backends/io_json.py
RefactoringTool: ./ipi/utils/io/backends/io_pdb.py
RefactoringTool: ./ipi/utils/io/backends/io_xyz.py
RefactoringTool: ./ipi/utils/io/inputs/__init__.py
RefactoringTool: ./ipi/utils/io/inputs/io_xml.py
RefactoringTool: ./ipi_tests/pdb_generator.py
RefactoringTool: ./ipi_tests/xyz_generator.py
RefactoringTool: ./ipi_tests/engine/test_initializer.py
RefactoringTool: ./ipi_tests/engine/test_properties.py
RefactoringTool: ./ipi_tests/utils/test_depend.py
RefactoringTool: ./ipi_tests/utils/io/backends/test__init__.py
RefactoringTool: ./ipi_tests/utils/io/backends/test_io_units.py
RefactoringTool: ./ipi_tests/utils/io/backends/test_io_xyz.py
RefactoringTool: ./tests/test_docs.py
RefactoringTool: ./tools/py/Instanton_interpolation.py
RefactoringTool: ./tools/py/Instanton_postproc.py
RefactoringTool: ./tools/py/a2b.py
RefactoringTool: ./tools/py/bin2xyz.py
RefactoringTool: ./tools/py/contract-trajectory.py
RefactoringTool: ./tools/py/effective_temperatures.py
RefactoringTool: ./tools/py/energies_ppi.py
RefactoringTool: ./tools/py/energy_ppi.py
RefactoringTool: ./tools/py/estmod_example.py
RefactoringTool: ./tools/py/fixcom.py
RefactoringTool: ./tools/py/get_np.py
RefactoringTool: ./tools/py/get_np_rad.py
RefactoringTool: ./tools/py/get_np_vec.py
RefactoringTool: ./tools/py/get_np_xyz.py
RefactoringTool: ./tools/py/getacf.py
RefactoringTool: ./tools/py/getproperty.py
RefactoringTool: ./tools/py/gleacf.py
RefactoringTool: ./tools/py/heat_capacity_ppi.py
RefactoringTool: ./tools/py/kinetic2tag.py
RefactoringTool: ./tools/py/kinetic_energy_ppi.py
RefactoringTool: ./tools/py/mergebeadspdb.py
RefactoringTool: ./tools/py/mux-positions.py
RefactoringTool: ./tools/py/paraweights.py
RefactoringTool: ./tools/py/pepper.py
RefactoringTool: ./tools/py/planetary.py
RefactoringTool: ./tools/py/posforce2kinetic.py
RefactoringTool: ./tools/py/potential_energy_ppi.py
RefactoringTool: ./tools/py/rdf_ppi.py
RefactoringTool: ./tools/py/regtest-parallel.py
RefactoringTool: ./tools/py/regtest.py
RefactoringTool: ./tools/py/remdsort.py
RefactoringTool: ./tools/py/total_energy_ppi.py
RefactoringTool: ./tools/py/trimsim.py
RefactoringTool: ./tools/py/xyz2bin.py
RefactoringTool: ./tools/py/xyz2pdb.py
ement[1]
+            print('Difference in line', element[0], 'in word', element[1])
 
     def files_are_equal(self):
         """
@@ -942,7 +942,7 @@
 
     _yes = ['yes', 'y']
     _no = ['no', 'n']
-    answer = raw_input()
+    answer = input()
     if answer.lower() in _yes:
         return True
     elif answer.lower() in _no:
--- ./tools/py/remdsort.py	(original)
+++ ./tools/py/remdsort.py	(refactored)
@@ -137,7 +137,7 @@
     # now reads files one frame at a time, and re-direct output to the appropriate location
 
     line = ptfile.readline().split()
-    irep = range(nsys)  # Could this be harmful?
+    irep = list(range(nsys))  # Could this be harmful?
     step = 0
     while True:
         # reads one line from index file
--- ./tools/py/total_energy_ppi.py	(original)
+++ ./tools/py/total_energy_ppi.py	(refactored)
@@ -66,24 +66,24 @@
     # check that we found the same number of positions and forces files
     nbeads = len(fns_pos)
     if nbeads != len(fns_for):
-        print fns_pos
-        print fns_for
+        print(fns_pos)
+        print(fns_for)
         raise ValueError("Mismatch between number of input files for forces and positions.")
 
     # print some information
-    print 'temperature = {:f} K'.format(float(temp))
-    print
-    print 'number of beads = {:d}'.format(nbeads)
-    print
-    print 'positions and forces file names:'
+    print('temperature = {:f} K'.format(float(temp)))
+    print()
+    print('number of beads = {:d}'.format(nbeads))
+    print()
+    print('positions and forces file names:')
     for fn_pos, fn_for in zip(fns_pos, fns_for):
-        print '{:s}   {:s}'.format(fn_pos, fn_for)
-    print
-    print 'potential energy file: {:s}'.format(fns_iU)
-    print
-    print 'output file name:'
-    print fn_out_en
-    print
+        print('{:s}   {:s}'.format(fn_pos, fn_for))
+    print()
+    print('potential energy file: {:s}'.format(fns_iU))
+    print()
+    print('output file name:')
+    print(fn_out_en)
+    print()
 
     # open input and output files
     ipos = [open(fn, "r") for fn in fns_pos]
@@ -116,7 +116,7 @@
     while True:  # Reading input files and calculating PPI correction
 
         if ifr % 100 == 0:
-            print '\rProcessing frame {:d}'.format(ifr),
+            print('\rProcessing frame {:d}'.format(ifr), end=' ')
             sys.stdout.flush()
 
         try: