diff --git a/c/vasp/.gitignore b/c/vasp/.gitignore
new file mode 100644
index 00000000..0d20b648
--- /dev/null
+++ b/c/vasp/.gitignore
@@ -0,0 +1 @@
+*.pyc
diff --git a/python/vasp/plotools.py b/python/vasp/plotools.py
index c1c087f1..eb229b7a 100644
--- a/python/vasp/plotools.py
+++ b/python/vasp/plotools.py
@@ -1,6 +1,8 @@
 
 import itertools as it
 import numpy as np
+from proj_group import ProjectorGroup
+from proj_shell import ProjectorShell
 
 np.set_printoptions(suppress=True)
 
@@ -18,37 +20,6 @@ def issue_warning(message):
     print "  !!! WARNING !!!: " + message
     print
 
-class Projector:
-    """
-    Class describing a local-orbital projector.
-    """
-  
-    def __init__(self, matrix, ib1=1, ib2=None, nion=1):
-        self.p_matrix = matrix.astype(np.complex128)
-        self.norb, self.nb = matrix.shape
-        self.nion = nion
-        self.ib1 = ib1 - 1
-        if not ib2 is None:
-            self.ib2 = ib2 - 1
-        else:
-            self.ib2 = self.nb - 1
-  
-    def project_up(self, mat):
-        return np.dot(self.p_matrix.conj().T, np.dot(mat, self.p_matrix))
-  
-    def project_down(self, mat):
-        assert mat.shape == (self.nb, self.nb), "  Matrix must match projector in size"
-        return np.dot(self.p_matrix, np.dot(mat, self.p_matrix.conj().T))
-  
-    def orthogonalize(self):
-        """
-        Orthogonalizes a projector.
-        Returns an overlap matrix and its eigenvalues for initial projectors.
-        """
-        self.p_matrix, overlap, over_eig = orthogonalize_projector(self.p_matrix)
-
-        return (overlap, over_eig)
-        
 ################################################################################
 # check_data_consistency()
 ################################################################################
diff --git a/python/vasp/proj_group.py b/python/vasp/proj_group.py
index 53f98698..e1774488 100644
--- a/python/vasp/proj_group.py
+++ b/python/vasp/proj_group.py
@@ -3,95 +3,6 @@ import numpy as np
 
 np.set_printoptions(suppress=True)
 
-################################################################################
-#
-# orthogonalize_projector_matrix()
-#
-################################################################################
-def orthogonalize_projector_matrix(p_matrix):
-    """
-    Orthogonalizes a projector defined by a rectangular matrix `p_matrix`.
-
-    Parameters
-    ----------
-
-    p_matrix (numpy.array[complex]) : matrix `Nm x Nb`, where `Nm` is
-      the number of orbitals, `Nb` number of bands
-
-    Returns
-    -------
-
-    Orthogonalized projector matrix, initial overlap matrix and its eigenvalues.
-    """
-# Overlap matrix O_{m m'} = \sum_{v} P_{m v} P^{*}_{v m'}
-    overlap = np.dot(p_matrix, p_matrix.conj().T)
-# Calculate [O^{-1/2}]_{m m'}
-    eig, eigv = np.linalg.eigh(overlap)
-    assert np.all(eig > 0.0), ("Negative eigenvalues of the overlap matrix:"
-       "projectors are ill-defined")
-    sqrt_eig = np.diag(1.0 / np.sqrt(eig))
-    shalf = np.dot(eigv, np.dot(sqrt_eig, eigv.conj().T))
-# Apply \tilde{P}_{m v} = \sum_{m'} [O^{-1/2}]_{m m'} P_{m' v}
-    p_ortho = np.dot(shalf, p_matrix)
-
-    return (p_ortho, overlap, eig)
-
-################################################################################
-#
-# select_bands()
-#
-################################################################################
-def select_bands(eigvals, emin, emax):
-    """
-    Select a subset of bands lying within a given energy window.
-    The band energies are assumed to be sorted in an ascending order.
-
-    Parameters
-    ----------
-    
-    eigvals (numpy.array) : all eigenvalues
-    emin, emax (float) : energy window
-
-    Returns
-    -------
-
-    ib_win, nb_min, nb_max : 
-    """
-# Sanity check
-    if emin > eigvals.max() or emax < eigvals.min():
-        raise Exception("Energy window does not overlap with the band structure")
-
-    nk, nband, ns_band = eigvals.shape
-    ib_win = np.zeros((nk, ns_band, 2), dtype=np.int32)
-
-    ib_min = 10000000
-    ib_max = 0
-    for isp in xrange(ns_band):
-        for ik in xrange(nk):
-            for ib in xrange(nband):
-                en = eigvals[ik, ib, isp]
-                if en >= emin:
-                    break
-            ib1 = ib
-            for ib in xrange(ib1, nband):
-                en = eigvals[ik, ib, isp]
-                if en > emax:
-                    break
-            else:
-# If we reached the last band add 1 to get the correct bound
-                ib += 1
-            ib2 = ib - 1
-
-            assert ib1 <= ib2, "No bands inside the window for ik = %s"%(ik)
-
-            ib_win[ik, isp, 0] = ib1
-            ib_win[ik, isp, 1] = ib2
-
-            ib_min = min(ib_min, ib1)
-            ib_max = max(ib_max, ib2)
-
-    return ib_win, ib_min, ib_max
-
 ################################################################################
 ################################################################################
 #
@@ -125,7 +36,7 @@ class ProjectorGroup:
         self.shells = shells
 
 # Determine the minimum and maximum band numbers
-        ib_win, ib_min, ib_max = select_bands(eigvals, self.emin, self.emax)
+        ib_win, ib_min, ib_max = self.select_bands(eigvals)
         self.ib_win = ib_win
         self.ib_min = ib_min
         self.ib_max = ib_max
@@ -254,7 +165,7 @@ class ProjectorGroup:
                         i1, i2 = blocks[ion]
                         p_mat[i1:i2, :nb] = shell.proj_win[ion, isp, ik, :nlm, :nb]
 # Now orthogonalize the obtained block projector
-                p_orth, overl, eig = orthogonalize_projector_matrix(p_mat)
+                p_orth, overl, eig = self.orthogonalize_projector_matrix(p_mat)
 #                print "ik = ", ik
 #                print overl.real
 # Distribute back projectors in the same order
@@ -265,5 +176,93 @@ class ProjectorGroup:
                         i1, i2 = blocks[ion]
                         shell.proj_win[ion, isp, ik, :nlm, :nb] = p_orth[i1:i2, :nb]
 
+################################################################################
+#
+# orthogonalize_projector_matrix()
+#
+################################################################################
+    def orthogonalize_projector_matrix(self, p_matrix):
+        """
+        Orthogonalizes a projector defined by a rectangular matrix `p_matrix`.
+
+        Parameters
+        ----------
+
+        p_matrix (numpy.array[complex]) : matrix `Nm x Nb`, where `Nm` is
+          the number of orbitals, `Nb` number of bands
+
+        Returns
+        -------
+
+        Orthogonalized projector matrix, initial overlap matrix and its eigenvalues.
+        """
+# Overlap matrix O_{m m'} = \sum_{v} P_{m v} P^{*}_{v m'}
+        overlap = np.dot(p_matrix, p_matrix.conj().T)
+# Calculate [O^{-1/2}]_{m m'}
+        eig, eigv = np.linalg.eigh(overlap)
+        assert np.all(eig > 0.0), ("Negative eigenvalues of the overlap matrix:"
+           "projectors are ill-defined")
+        sqrt_eig = np.diag(1.0 / np.sqrt(eig))
+        shalf = np.dot(eigv, np.dot(sqrt_eig, eigv.conj().T))
+# Apply \tilde{P}_{m v} = \sum_{m'} [O^{-1/2}]_{m m'} P_{m' v}
+        p_ortho = np.dot(shalf, p_matrix)
+
+        return (p_ortho, overlap, eig)
+
+################################################################################
+#
+# select_bands()
+#
+################################################################################
+    def select_bands(self, eigvals):
+        """
+        Select a subset of bands lying within a given energy window.
+        The band energies are assumed to be sorted in an ascending order.
+
+        Parameters
+        ----------
+        
+        eigvals (numpy.array) : all eigenvalues
+        emin, emax (float) : energy window
+
+        Returns
+        -------
+
+        ib_win, nb_min, nb_max : 
+        """
+# Sanity check
+        if self.emin > eigvals.max() or self.emax < eigvals.min():
+            raise Exception("Energy window does not overlap with the band structure")
+
+        nk, nband, ns_band = eigvals.shape
+        ib_win = np.zeros((nk, ns_band, 2), dtype=np.int32)
+
+        ib_min = 10000000
+        ib_max = 0
+        for isp in xrange(ns_band):
+            for ik in xrange(nk):
+                for ib in xrange(nband):
+                    en = eigvals[ik, ib, isp]
+                    if en >= self.emin:
+                        break
+                ib1 = ib
+                for ib in xrange(ib1, nband):
+                    en = eigvals[ik, ib, isp]
+                    if en > self.emax:
+                        break
+                else:
+# If we reached the last band add 1 to get the correct bound
+                    ib += 1
+                ib2 = ib - 1
+
+                assert ib1 <= ib2, "No bands inside the window for ik = %s"%(ik)
+
+                ib_win[ik, isp, 0] = ib1
+                ib_win[ik, isp, 1] = ib2
+
+                ib_min = min(ib_min, ib1)
+                ib_max = max(ib_max, ib2)
+
+        return ib_win, ib_min, ib_max