diff --git a/nexus/lib/basisset.py b/nexus/lib/basisset.py
index a2ed76009..11e97aba5 100644
--- a/nexus/lib/basisset.py
+++ b/nexus/lib/basisset.py
@@ -975,7 +975,7 @@ class GaussianBasisSet(DevBase):
                 else:
                     plotter = plot
                 #end if
-                plotter(range(len(lwidths)),lwidths,lstyle,label=l)
+                plotter(list(range(len(lwidths))),lwidths,lstyle,label=l)
             #end if
         #end for
         if label:
diff --git a/nexus/lib/fileio.py b/nexus/lib/fileio.py
index bd1159de9..59c2ff27f 100644
--- a/nexus/lib/fileio.py
+++ b/nexus/lib/fileio.py
@@ -933,7 +933,7 @@ class XsfFile(StandardFile):
         ndim = 3
         permute = dim!=0
         if permute:
-            r = range(0,ndim)
+            r = list(range(0,ndim))
             r.pop(dim)
             permutation = tuple([dim]+r)
             data = data.transpose(permutation)
@@ -948,7 +948,7 @@ class XsfFile(StandardFile):
 
     def line_plot(self,dim,filepath):
         r,d = self.line_data(dim)
-        savetxt(filepath,array(zip(r,d)))
+        savetxt(filepath,array(list(zip(r,d))))
     #end def line_plot
 #end class XsfFile
 
diff --git a/nexus/lib/nexus.py b/nexus/lib/nexus.py
index 46548df8e..f47c2b981 100644
--- a/nexus/lib/nexus.py
+++ b/nexus/lib/nexus.py
@@ -224,7 +224,7 @@ class Settings(NexusCore):
         self.process_noncore_settings(kw)
 
         # transfer select core data to the global namespace
-        nexus_core_noncore.transfer_from(nexus_core,nexus_core_noncore.keys())
+        nexus_core_noncore.transfer_from(nexus_core,list(nexus_core_noncore.keys()))
         nexus_noncore.set(**nexus_core_noncore.copy()) # prevent write to core namespace
 
         # copy final core and noncore settings
diff --git a/nexus/lib/numerics.py b/nexus/lib/numerics.py
index dfbf013c4..ed8a7dca2 100644
--- a/nexus/lib/numerics.py
+++ b/nexus/lib/numerics.py
@@ -797,11 +797,11 @@ def simstats(x,dim=None):
     reshape = ndim>2
     nblocks = shape[dim]
     if permute:
-        r = range(ndim)
+        r = list(range(ndim))
         r.pop(dim)
         r.append(dim)
         permutation = tuple(r)
-        r = range(ndim)
+        r = list(range(ndim))
         r.pop(ndim-1)
         r.insert(dim,ndim-1)
         invperm     = tuple(r)
diff --git a/nexus/lib/qmcpack.py b/nexus/lib/qmcpack.py
index fe871390d..85880bd61 100644
--- a/nexus/lib/qmcpack.py
+++ b/nexus/lib/qmcpack.py
@@ -120,7 +120,7 @@ class Qmcpack(Simulation):
         if self.system is None:
             self.should_twist_average = False
         elif nexus_core.generate_only:
-            twistnums = range(len(self.system.structure.kpoints))
+            twistnums = list(range(len(self.system.structure.kpoints)))
             if self.should_twist_average:
                 self.twist_average(twistnums)
             #end if
@@ -222,7 +222,7 @@ class Qmcpack(Simulation):
                     self.error('wavefunction file not found:\n'+h5file)
                 #end if
 
-                twistnums = range(len(structure.kpoints))
+                twistnums = list(range(len(structure.kpoints)))
                 if self.should_twist_average:
                     self.twist_average(twistnums)
                 elif not has_twist and orb_elem.twistnum is None:
@@ -446,7 +446,7 @@ class Qmcpack(Simulation):
             output_files = []
         else:
             if self.should_twist_average and not isinstance(self.input,TracedQmcpackInput):
-                self.twist_average(range(len(self.system.structure.kpoints)))
+                self.twist_average(list(range(len(self.system.structure.kpoints))))
                 br = self.bundle_request
                 input = self.input.trace(br.quantity,br.values)
                 input.generate_filenames(self.infile)
diff --git a/nexus/lib/qmcpack_quantity_analyzers.py b/nexus/lib/qmcpack_quantity_analyzers.py
index 6c544cea4..80ea52c79 100644
--- a/nexus/lib/qmcpack_quantity_analyzers.py
+++ b/nexus/lib/qmcpack_quantity_analyzers.py
@@ -3104,7 +3104,7 @@ class RectilinearGrid(SpaceGridBase):
                     self.error('ndu is different than len(gmap)')
                 #end if
                 du = 1./self.odu[d]
-                fine_interval_centers[d] = self.umin + .5*du + du*array(range(ndu))
+                fine_interval_centers[d] = self.umin + .5*du + du*array(list(range(ndu)))
                 find_interval_domains[d] = zeros((ndu,))
             #end for
             #checks are done on each source spacegrid to determine interpolation compatibility 
diff --git a/nexus/lib/qmcpack_result_analyzers.py b/nexus/lib/qmcpack_result_analyzers.py
index 31b1be02d..601c5b5b1 100644
--- a/nexus/lib/qmcpack_result_analyzers.py
+++ b/nexus/lib/qmcpack_result_analyzers.py
@@ -294,7 +294,7 @@ class OptimizationAnalyzer(ResultAnalyzer):
 
         #plot energy and variance
         figure()
-        r = range(nopt)
+        r = list(range(nopt))
         subplot(3,1,1)
         errorbar(r,en,enerr,fmt='b')
         ylabel('Energy (Ha)')
diff --git a/nexus/lib/structure.py b/nexus/lib/structure.py
index 2de709d3f..3fae324d4 100644
--- a/nexus/lib/structure.py
+++ b/nexus/lib/structure.py
@@ -379,7 +379,7 @@ def reduce_tilematrix(tiling):
         tilematrix = T.copy()
         del t
         tbar = identity(dim) #basis for shearing
-        dr = range(dim)
+        dr = list(range(dim))
         #dr = [1,0,2]
         other = dim*[0] # other[d] = dimensions other than d
         for d in dr: 
@@ -1967,7 +1967,7 @@ class Structure(Sobj):
             return
         elif magnetization=='':
             magnetization = identifiers
-            indices = range(len(self.elem))
+            indices = list(range(len(self.elem)))
         else:
             indices = self.locate(identifiers)
         #end if
@@ -2040,14 +2040,14 @@ class Structure(Sobj):
         nrem = len(indices)
         nadd = len(pos)
         if nadd<nrem:
-            ar = array(range(0,nadd))
-            rr = array(range(nadd,nrem))
+            ar = array(list(range(0,nadd)))
+            rr = array(list(range(nadd,nrem)))
             self.elem[indices[ar]] = elem[:]
             self.pos[indices[ar]]  = pos[:]
             self.remove(indices[rr])
         elif nadd>nrem:
-            ar = array(range(0,nrem))
-            er = array(range(nrem,nadd))
+            ar = array(list(range(0,nrem)))
+            er = array(list(range(nrem,nadd)))
             self.elem[indices[ar]] = elem[ar]
             self.pos[indices[ar]]  = pos[ar]
             ii = indices[ar[-1]]
@@ -2072,7 +2072,7 @@ class Structure(Sobj):
         np = len(pos)
         nps= len(self.pos)
         d = empty((np,))
-        ip = array(range(np))
+        ip = array(list(range(np)))
         ips= nn.ravel()
         for i in ip:
             j = ips[i]
@@ -2951,7 +2951,7 @@ class Structure(Sobj):
         # get nearest neighbor index pairs in the large cell
         neigh_pairs = voronoi_neighbors(ss.pos)
         # create a mapping from large to small indices
-        large_to_small = 3**d*range(len(self.pos))
+        large_to_small = 3**d*list(range(len(self.pos)))
         # find the neighbor pairs in the small cell
         neighbors = obj()
         small_inds = set(ss.locate(sn.pos))
@@ -3262,7 +3262,7 @@ class Structure(Sobj):
             #end if
         #end for
         npos,dim = pos.shape
-        drange = range(dim)
+        drange = list(range(dim))
         n = len(surface)
         i=0
         while i<n:
@@ -3600,7 +3600,7 @@ class Structure(Sobj):
         kindices = []
         kpoints  = []
         shift = empty((ndim,))
-        kr = range(nkpoints)
+        kr = list(range(nkpoints))
         for k in range(imin[2],imax[2]+1):
             for j in range(imin[1],imax[1]+1):
                 for i in range(imin[0],imax[0]+1):
@@ -4164,7 +4164,7 @@ class Structure(Sobj):
         #end if
         ref = reference.tile((3,3,3))
         ref.recenter(reference.center)
-        rmap = array(3**3*range(len(reference.pos)),dtype=int)
+        rmap = array(3**3*list(range(len(reference.pos)),dtype=int))
         nn = nearest_neighbors(1,ref.pos,self.pos).ravel()
         displacement = self.pos - ref.pos[nn]
         if not map:
@@ -4203,8 +4203,8 @@ class Structure(Sobj):
             rbonds = rbv[ir]
             ib  = empty((neighbors,),dtype=int)
             ibr = empty((neighbors,),dtype=int)
-            r  = range(neighbors)
-            rr = range(neighbors)
+            r  = list(range(neighbors))
+            rr = list(range(neighbors))
             for n in range(neighbors):
                 mindist = 1e99
                 ibmin  = -1
@@ -4322,7 +4322,7 @@ class Structure(Sobj):
         ntile = ncells*natoms
         pos = empty((ntile,dim))
         ind = empty((ntile,),dtype=int)
-        oind = range(natoms)
+        oind = list(range(natoms))
         for nt in range(ncells):
             n=nt*natoms
             ind[n:n+natoms]=oind[:]
@@ -4695,7 +4695,7 @@ class Structure(Sobj):
         self.pos = pos
         if selective_dynamics or spos.shape[1]>3:
             move = array(spos[:,3:6],dtype=str)
-            self.freeze(range(self.size()),directions=move=='F')
+            self.freeze(list(range(self.size())),directions=move=='F')
         #end if
     #end def read_poscar
 
@@ -5631,7 +5631,7 @@ class DefectStructure(Structure):
         nn = nearest_neighbors(1,dlarge,dsmall)
         dist = dsmall.distances(dlarge[nn.ravel()])
         dmatch = dist<dist_cutoff
-        ismall = array(range(len(dsmall.pos)))
+        ismall = array(list(range(len(dsmall.pos))))
         ismall = ismall[dmatch]
         ilarge = nn[ismall]
         if natoms1<natoms2: