diff --git a/PKG-INFO b/PKG-INFO
index c0d3897..bf3a165 100644
--- a/PKG-INFO
+++ b/PKG-INFO
@@ -1,6 +1,6 @@
 Metadata-Version: 1.1
 Name: pythtb
-Version: 1.7.1
+Version: 1.7.2
 Summary: Simple solver for tight binding models for use in condensed matter physics and materials science.
 Home-page: http://www.physics.rutgers.edu/pythtb
 Author: Sinisa Coh and David Vanderbilt
diff --git a/examples/0dim.py b/examples/0dim.py
old mode 100755
new mode 100644
index 04caa67..727fec8
--- a/examples/0dim.py
+++ b/examples/0dim.py
@@ -8,7 +8,7 @@
 from __future__ import print_function
 from pythtb import * # import TB model class
 import numpy as np
-import pylab as plt
+import matplotlib.pyplot as plt
 
 # define lattice vectors
 lat=[[1.0,0.0,0.0],[0.0,1.0,0.0],[0.0,0.0,1.0]]
diff --git a/examples/3site_cycle.py b/examples/3site_cycle.py
old mode 100755
new mode 100644
index c1b46d9..4ee56ea
--- a/examples/3site_cycle.py
+++ b/examples/3site_cycle.py
@@ -9,7 +9,7 @@
 from __future__ import print_function
 from pythtb import * # import TB model class
 import numpy as np
-import pylab as plt
+import matplotlib.pyplot as plt
 
 # define lattice vectors
 lat=[[1.0]]
@@ -101,6 +101,7 @@
 ax_band.set_title("Band structure")
 ax_band.set_xlabel("Path in k-vector")
 ax_band.set_ylabel("Band energies")
+ax_band.set_xlim(0.0,1.0)
 fig_band.tight_layout()
 fig_band.savefig("3site_band.pdf")
 # finish plot for Wannier center
diff --git a/examples/3site_cycle_fin.py b/examples/3site_cycle_fin.py
old mode 100755
new mode 100644
index c938d43..aaa054e
--- a/examples/3site_cycle_fin.py
+++ b/examples/3site_cycle_fin.py
@@ -9,7 +9,7 @@
 from __future__ import print_function
 from pythtb import * # import TB model class
 import numpy as np
-import pylab as plt
+import matplotlib.pyplot as plt
 
 # define function to construct model
 def set_model(t,delta,lmbd):
diff --git a/examples/buckled_layer.py b/examples/buckled_layer.py
old mode 100755
new mode 100644
index f47cd88..ecb2fc9
--- a/examples/buckled_layer.py
+++ b/examples/buckled_layer.py
@@ -9,7 +9,7 @@
 from __future__ import print_function
 from pythtb import * # import TB model class
 import numpy as np
-import pylab as plt
+import matplotlib.pyplot as plt
 
 # real space is 3D
 # define lattice vectors
@@ -74,7 +74,7 @@
 ax.set_ylabel("Band energy")
 
 # specify horizontal axis details
-ax.set_xlim([0,k_node[-1]])
+ax.set_xlim(k_node[0],k_node[-1])
 # put tickmarks and labels at node positions
 ax.set_xticks(k_node)
 ax.set_xticklabels(label)
diff --git a/examples/checkerboard.py b/examples/checkerboard.py
old mode 100755
new mode 100644
index 1eb8018..8b47e00
--- a/examples/checkerboard.py
+++ b/examples/checkerboard.py
@@ -8,7 +8,7 @@
 from __future__ import print_function
 from pythtb import * # import TB model class
 import numpy as np
-import pylab as plt
+import matplotlib.pyplot as plt
 
 # define lattice vectors
 lat=[[1.0,0.0],[0.0,1.0]]
@@ -55,7 +55,7 @@
 fig, ax = plt.subplots()
 
 # specify horizontal axis details
-ax.set_xlim([0,k_node[-1]])
+ax.set_xlim(k_node[0],k_node[-1])
 ax.set_xticks(k_node)
 ax.set_xticklabels(label)
 for n in range(len(k_node)):
diff --git a/examples/cone.py b/examples/cone.py
old mode 100755
new mode 100644
index 73acca0..6b3d856
--- a/examples/cone.py
+++ b/examples/cone.py
@@ -9,7 +9,7 @@
 from __future__ import print_function
 from pythtb import * # import TB model class
 import numpy as np
-import pylab as plt
+import matplotlib.pyplot as plt
 
 # define lattice vectors
 lat=[[1.0,0.0],[0.5,np.sqrt(3.0)/2.0]]
diff --git a/examples/edge.py b/examples/edge.py
old mode 100755
new mode 100644
diff --git a/examples/graphene.py b/examples/graphene.py
old mode 100755
new mode 100644
index f99f5f3..a32855b
--- a/examples/graphene.py
+++ b/examples/graphene.py
@@ -8,7 +8,7 @@
 from __future__ import print_function
 from pythtb import * # import TB model class
 import numpy as np
-import pylab as plt
+import matplotlib.pyplot as plt
 
 # define lattice vectors
 lat=[[1.0,0.0],[0.5,np.sqrt(3.0)/2.0]]
@@ -64,7 +64,7 @@
 fig, ax = plt.subplots()
 # specify horizontal axis details
 # set range of horizontal axis
-ax.set_xlim([0,k_node[-1]])
+ax.set_xlim(k_node[0],k_node[-1])
 # put tickmarks and labels at node positions
 ax.set_xticks(k_node)
 ax.set_xticklabels(label)
diff --git a/examples/haldane.py b/examples/haldane.py
old mode 100755
new mode 100644
index c2840ed..79ac4f5
--- a/examples/haldane.py
+++ b/examples/haldane.py
@@ -8,7 +8,7 @@
 from __future__ import print_function
 from pythtb import * # import TB model class
 import numpy as np
-import pylab as plt
+import matplotlib.pyplot as plt
 
 # define lattice vectors
 lat=[[1.0,0.0],[0.5,np.sqrt(3.0)/2.0]]
@@ -66,7 +66,7 @@
 fig, ax = plt.subplots()
 # specify horizontal axis details
 # set range of horizontal axis
-ax.set_xlim([0,k_node[-1]])
+ax.set_xlim(k_node[0],k_node[-1])
 # put tickmarks and labels at node positions
 ax.set_xticks(k_node)
 ax.set_xticklabels(label)
diff --git a/examples/haldane_bp.py b/examples/haldane_bp.py
old mode 100755
new mode 100644
index 3432c22..9874aa3
--- a/examples/haldane_bp.py
+++ b/examples/haldane_bp.py
@@ -9,7 +9,7 @@
 from __future__ import print_function
 from pythtb import * # import TB model class
 import numpy as np
-import pylab as plt
+import matplotlib.pyplot as plt
 
 # define lattice vectors
 lat=[[1.0,0.0],[0.5,np.sqrt(3.0)/2.0]]
@@ -75,6 +75,7 @@
 ax.set_title("Berry phase for lower (red), top (green), both bands (blue)")
 ax.set_xlabel(r"$k_y$")
 ax.set_ylabel(r"Berry phase along $k_x$")
+ax.set_xlim(0.,1.)
 ax.set_ylim(-7.,7.)
 ax.yaxis.set_ticks([-2.*np.pi,-np.pi,0.,np.pi,2.*np.pi])
 ax.set_yticklabels((r'$-2\pi$',r'$-\pi$',r'$0$',r'$\pi$', r'$2\pi$'))
diff --git a/examples/haldane_fin.py b/examples/haldane_fin.py
old mode 100755
new mode 100644
index 221234b..4354f75
--- a/examples/haldane_fin.py
+++ b/examples/haldane_fin.py
@@ -9,7 +9,7 @@
 from __future__ import print_function
 from pythtb import * # import TB model class
 import numpy as np
-import pylab as plt
+import matplotlib.pyplot as plt
 
 # define lattice vectors
 lat=[[1.0,0.0],[0.5,np.sqrt(3.0)/2.0]]
diff --git a/examples/haldane_hwf.py b/examples/haldane_hwf.py
old mode 100755
new mode 100644
index 8550cca..de8a093
--- a/examples/haldane_hwf.py
+++ b/examples/haldane_hwf.py
@@ -12,7 +12,7 @@
 from __future__ import print_function
 from pythtb import * # import TB model class
 import numpy as np
-import pylab as plt
+import matplotlib.pyplot as plt
 
 # set model parameters
 delta=-0.2
@@ -128,7 +128,7 @@
   
 # label both axes
 for ax in [ax1,ax2]:
-  ax.set_xlim([0.,k_node[-1]])
+  ax.set_xlim(k_node[0],k_node[-1])
   ax.set_xticks(k_node)
   ax.set_xticklabels(k_label)
 
diff --git a/examples/kane_mele.py b/examples/kane_mele.py
old mode 100755
new mode 100644
index 2297df3..c139280
--- a/examples/kane_mele.py
+++ b/examples/kane_mele.py
@@ -9,7 +9,7 @@
 from __future__ import print_function
 from pythtb import * # import TB model class
 import numpy as np
-import pylab as plt
+import matplotlib.pyplot as plt
 
 def get_kane_mele(topological):
   "Return a Kane-Mele model in the normal or topological phase."
@@ -93,6 +93,7 @@ def get_kane_mele(topological):
   ax1.set_title("Kane-Mele: "+top_index+" phase")
   ax1.set_xticks(k_node)
   ax1.set_xticklabels(label)
+  ax1.set_xlim(k_node[0],k_node[-1])
   for n in range(len(k_node)):
     ax1.axvline(x=k_node[n],linewidth=0.5, color='k')
   ax1.set_xlabel("k-space")
@@ -127,6 +128,7 @@ def get_kane_mele(topological):
   ax2.set_ylabel('Wannier center along x')
   ax2.set_xlabel(r'$k_y$')
   ax2.set_xticks([0.0,0.5,1.0])
+  ax2.set_xlim(0.0,1.0)
   ax2.set_xticklabels([r"$0$",r"$\pi$", r"$2\pi$"])
   ax2.axvline(x=.5,linewidth=0.5, color='k')
   ax2.set_title("1D Wannier centers: "+top_index+" phase")
diff --git a/examples/simple.py b/examples/simple.py
old mode 100755
new mode 100644
index a47d748..8e11a56
--- a/examples/simple.py
+++ b/examples/simple.py
@@ -6,7 +6,7 @@
 # by Sinisa Coh and David Vanderbilt (see gpl-pythtb.txt)
 
 from pythtb import * # import TB model class
-import pylab as plt
+import matplotlib.pyplot as plt
 
 # specify model
 lat=[[1.0]]
@@ -29,6 +29,7 @@
 ax.set_ylabel("Band energy")
 ax.set_xticks(k_node)
 ax.set_xticklabels(k_label)
+ax.set_xlim(k_node[0],k_node[-1])
 for n in range(len(k_node)):
   ax.axvline(x=k_node[n], linewidth=0.5, color='k')
 fig.tight_layout()
diff --git a/examples/supercell.py b/examples/supercell.py
old mode 100755
new mode 100644
index 8884f8d..be03941
--- a/examples/supercell.py
+++ b/examples/supercell.py
@@ -9,7 +9,7 @@
 from __future__ import print_function
 from pythtb import * # import TB model class
 import numpy as np
-import pylab as plt
+import matplotlib.pyplot as plt
 
 # define lattice vectors
 lat=[[1.0,0.0],[0.5,np.sqrt(3.0)/2.0]]
diff --git a/examples/trestle.py b/examples/trestle.py
old mode 100755
new mode 100644
index be74c18..9a08c77
--- a/examples/trestle.py
+++ b/examples/trestle.py
@@ -8,7 +8,7 @@
 from __future__ import print_function
 from pythtb import * # import TB model class
 import numpy as np
-import pylab as plt
+import matplotlib.pyplot as plt
 
 # define lattice vectors
 lat=[[2.0,0.0],[0.0,1.0]]
@@ -52,7 +52,7 @@
 # First make a figure object
 fig, ax = plt.subplots()
 # specify horizontal axis details
-ax.set_xlim([0,k_node[2]])
+ax.set_xlim(k_node[0],k_node[-1])
 ax.set_xticks(k_node)
 ax.set_xticklabels(k_label)
 ax.axvline(x=k_node[1],linewidth=0.5, color='k')
diff --git a/examples/visualize.py b/examples/visualize.py
old mode 100755
new mode 100644
diff --git a/examples/w90.py b/examples/w90.py
old mode 100755
new mode 100644
index a035775..9be1892
--- a/examples/w90.py
+++ b/examples/w90.py
@@ -5,7 +5,7 @@
 
 from __future__ import print_function
 from pythtb import * # import TB model class
-import pylab as plt
+import matplotlib.pyplot as plt
 
 # read output from Wannier90 that should be in folder named "example_a"
 # see instructions above for how to obtain the example output from Wannier90
diff --git a/examples/w90_quick.py b/examples/w90_quick.py
old mode 100755
new mode 100644
index 9f9b986..d10d767
--- a/examples/w90_quick.py
+++ b/examples/w90_quick.py
@@ -4,7 +4,7 @@
 # by Sinisa Coh and David Vanderbilt (see gpl-pythtb.txt)
 
 from pythtb import * # import TB model class
-import pylab as plt
+import matplotlib.pyplot as plt
 
 # read output from Wannier90 that should be in folder named "example_a"
 #   see instructions above for how to obtain the example output from 
diff --git a/pythtb.py b/pythtb.py
index 3b0ff4e..c8520e2 100644
--- a/pythtb.py
+++ b/pythtb.py
@@ -1,10 +1,10 @@
 from __future__ import print_function
 
 # PythTB python tight binding module.
-# December 22, 2016
-__version__='1.7.1'
+# August 25th, 2017
+__version__='1.7.2'
 
-# Copyright 2010, 2012, 2016 by Sinisa Coh and David Vanderbilt
+# Copyright 2010, 2012, 2016, 2017 by Sinisa Coh and David Vanderbilt
 #
 # This file is part of PythTB.  PythTB is free software: you can
 # redistribute it and/or modify it under the terms of the GNU General
@@ -38,7 +38,9 @@ class tb_model(object):
       real space lattice vectors there are and how many coordinates are
       needed to specify the orbital coordinates.
 
-    .. note:: Parameter *dim_r* can be larger than *dim_k*! For example,
+    .. note::
+
+      Parameter *dim_r* can be larger than *dim_k*! For example,
       a polymer is a three-dimensional molecule (one needs three
       coordinates to specify orbital positions), but it is periodic
       along only one direction. For a polymer, therefore, we should
@@ -375,12 +377,13 @@ def set_hop(self,hop_amp,ind_i,ind_j,ind_R=None,mode="set",allow_conjugate_pair=
           orbital does not have to be in the home unit cell; its unit cell
           position is determined by parameter *ind_R*.
 
-        :param ind_R: Specifies, in reduced coordinates, the shift of
-          the ket orbital. The number of coordinates must equal the
-          dimensionality in real space (*dim_r* parameter) for consistency,
-          but only the periodic directions of ind_R will be considered. If
-          reciprocal space is zero-dimensional (as in a molecule),
-          this parameter does not need to be specified.
+        :param ind_R: Lattice vector (integer array, in reduced
+          coordinates) pointing to the unit cell where the ket
+          orbital is located.  The number of coordinates must equal
+          the dimensionality in real space (*dim_r* parameter) for
+          consistency, but only the periodic directions of ind_R are
+          used. If reciprocal space is zero-dimensional (as in a
+          molecule), this parameter does not need to be specified.
 
         :param mode: Similar to parameter *mode* in function *set_onsite*. 
           Speficies way in which parameter *hop_amp* is
@@ -552,8 +555,7 @@ def _val_to_block(self,val):
                 ret[1,1]+=use_val[0]
                 # sigma_x
                 ret[0,1]+=use_val[1]
-                ret[1,0]+=use_val[1]
-                # sigma_y
+                ret[1,0]+=use_val[1]                # sigma_y
                 ret[0,1]+=use_val[2]*(-1.0j)
                 ret[1,0]+=use_val[2]*( 1.0j)
                 # sigma_z
@@ -623,6 +625,25 @@ def display(self):
                 print(_nice_complex(hopping[0],7,4))
             elif self._nspin==2:
                 print(str(hopping[0]).replace("\n"," "))
+        print('hopping distances:')
+        for i,hopping in enumerate(self._hoppings):
+            print("|  pos(",_nice_int(hopping[1],2),")  - pos(",_nice_int(hopping[2],2), end=' ')
+            if len(hopping)==4:
+                print("+ [", end=' ')
+                for j,v in enumerate(hopping[3]):
+                    print(_nice_int(v,2), end=' ')
+                    if j!=len(hopping[3])-1:
+                        print(",", end=' ')
+                    else:
+                        print("]", end=' ')
+            print(") |  =  ", end=' ')
+            pos_i=np.dot(self._orb[hopping[1]],self._lat)
+            pos_j=np.dot(self._orb[hopping[2]],self._lat)
+            if len(hopping)==4:
+                pos_j+=np.dot(hopping[3],self._lat)
+            dist=np.linalg.norm(pos_j-pos_i)
+            print (_nice_float(dist,7,4))
+
         print()
 
     def visualize(self,dir_first,dir_second=None,eig_dr=None,draw_hoppings=True,ph_color="black"):
@@ -722,7 +743,7 @@ def visualize(self,dir_first,dir_second=None,eig_dr=None,draw_hoppings=True,ph_c
             raise Exception("\n\nNeed to specify index of second coordinate for projection!")
 
         # start a new figure
-        import pylab as plt
+        import matplotlib.pyplot as plt
         fig=plt.figure(figsize=[plt.rcParams["figure.figsize"][0],
                                 plt.rcParams["figure.figsize"][0]])
         ax=fig.add_subplot(111, aspect='equal')
@@ -1517,6 +1538,80 @@ def _shift_to_home(self):
                             self._hoppings[h][3]+=disp_vec
 
 
+    def remove_orb(self,to_remove):
+        r"""
+
+        Returns a model with some orbitals removed.  Note that this
+        will reindex the orbitals with indices higher than those that
+        are removed.  For example.  If model has 6 orbitals and one
+        wants to remove 2nd orbital, then returned model will have 5
+        orbitals indexed as 0,1,2,3,4.  In the returned model orbital
+        indexed as 2 corresponds to the one indexed as 3 in the
+        original model.  Similarly 3 and 4 correspond to 4 and 5.
+        Indices of first two orbitals (0 and 1) are unaffected.
+
+        :param to_remove: List of orbital indices to be removed, or 
+          index of single orbital to be removed
+
+        :returns:
+
+          * **del_tb** -- Object of type :class:`pythtb.tb_model` 
+            representing a model with removed orbitals.
+
+        Example usage::
+        
+          # if original_model has say 10 orbitals then
+          # returned small_model will have only 8 orbitals.
+
+          small_model=original_model.remove_orb([2,5])
+
+        """
+
+        # if a single integer is given, convert to a list with one element
+        if type(to_remove).__name__=='int':
+            orb_index=[to_remove]
+        else:
+            orb_index=copy.deepcopy(to_remove)
+
+        # check range of indices
+        for i,orb_ind in enumerate(orb_index):
+            if orb_ind < 0 or orb_ind > self._norb-1 or type(orb_ind).__name__!='int':
+                raise Exception("\n\nSpecified wrong orbitals to remove!")
+        for i,ind1 in enumerate(orb_index):
+            for ind2 in orb_index[i+1:]:
+                if ind1==ind2:
+                    raise Exception("\n\nSpecified duplicate orbitals to remove!")
+
+        # put the orbitals to be removed in desceding order
+        orb_index = sorted(orb_index,reverse=True)
+
+        # make copy of a model
+        ret=copy.deepcopy(self)
+
+        # adjust some variables in the new model
+        ret._norb-=len(orb_index)
+        ret._nsta-=len(orb_index)*self._nspin
+        # remove indices one by one
+        for i,orb_ind in enumerate(orb_index):
+            # adjust variables
+            ret._orb = np.delete(ret._orb,orb_ind,0)
+            ret._site_energies = np.delete(ret._site_energies,orb_ind,0)
+            ret._site_energies_specified = np.delete(ret._site_energies_specified,orb_ind)
+            # adjust hopping terms (in reverse)
+            for j in range(len(ret._hoppings)-1,-1,-1):
+                h=ret._hoppings[j]
+                # remove all terms that involve this orbital
+                if h[1]==orb_ind or h[2]==orb_ind:
+                    del ret._hoppings[j]
+                else: # otherwise modify term
+                    if h[1]>orb_ind:
+                        ret._hoppings[j][1]-=1
+                    if h[2]>orb_ind:
+                        ret._hoppings[j][2]-=1
+        # return new model
+        return ret
+
+
     def k_uniform_mesh(self,mesh_size):
         r""" 
         Returns a uniform grid of k-points that can be passed to
@@ -2260,7 +2355,7 @@ def impose_pbc(self,mesh_dir,k_dir):
         :math:`\Psi_{n,{\bf k+G}}=\Psi_{n {\bf k}}` not only up to a
         phase, but they are also equal in phase.  It follows that
         the cell-periodic Bloch functions are related by
-        :math:`u_{n,{\bf k+G}}=e^{-i{\bf G}\cdot{\bf r}}\Psi_{n {\bf k}}`.
+        :math:`u_{n,{\bf k+G}}=e^{-i{\bf G}\cdot{\bf r}} u_{n {\bf k}}`.
         See :download:`notes on tight-binding formalism
         <misc/pythtb-formalism.pdf>` section 4.4 and equation 4.18 for
         more detail.  This routine sets the cell-periodic Bloch function
@@ -3354,7 +3449,7 @@ def dist_hop(self):
           (dist,ham)=silicon.dist_hop()
 
           # plot logarithm of the hopping term as a function of distance
-          import pylab as plt
+          import matplotlib.pyplot as plt
           fig, ax = plt.subplots()
           ax.scatter(dist,np.log(np.abs(ham)))
           fig.savefig("localization.pdf")
@@ -3478,7 +3573,7 @@ def w90_bands_consistency(self):
           evals=my_model.solve_all(w90_kpt)
           
           # plot comparison of the two
-          import pylab as plt
+          import matplotlib.pyplot as plt
           fig, ax = plt.subplots() 
           for i in range(evals.shape[0]):
               ax.plot(range(evals.shape[1]),evals[i],"r-",zorder=-50)
diff --git a/setup.py b/setup.py
old mode 100755
new mode 100644
index cdf2753..74f3317
--- a/setup.py
+++ b/setup.py
@@ -3,7 +3,7 @@
 from distutils.core import setup
 
 setup(name='pythtb',
-      version='1.7.1',
+      version='1.7.2',
       author='Sinisa Coh and David Vanderbilt',
       author_email='sinisacoh@gmail.com  dhv@physics.rutgers.edu',
       url='http://www.physics.rutgers.edu/pythtb',
diff --git a/tests/test_cone.py b/tests/test_cone.py
deleted file mode 100644
index f2f570b..0000000
--- a/tests/test_cone.py
+++ /dev/null
@@ -1,10 +0,0 @@
-import sys
-sys.path.append("../examples/")
-sys.path.append("../")
-
-from cone import *
-
-def test_answer():
-    assert np.isclose(w_square.berry_flux([0]),2.17921648013)
-    assert np.isclose(w_square.berry_flux([1]),-2.17921648013)
-    assert np.isclose(w_square.berry_flux([0,1]),0.0)
diff --git a/tests/test_edge.py b/tests/test_edge.py
deleted file mode 100644
index 39c694d..0000000
--- a/tests/test_edge.py
+++ /dev/null
@@ -1,10 +0,0 @@
-import sys
-sys.path.append("../examples/")
-sys.path.append("../")
-
-from edge import *
-
-def test_answer():
-    print(evecs[ed,:])
-    assert np.isclose(evecs[ed,0],2.80366220e-01+0.00000000e+00j)
-    assert np.isclose(evecs[ed,4],1.56548374e-01 +3.48852597e-02j)