diff --git a/crystal_mec.py b/crystal_mec.py
index 42c7fc8..b305993 100644
--- a/crystal_mec.py
+++ b/crystal_mec.py
@@ -31,7 +31,7 @@ if __name__=='__main__':
     inputFile = sys.argv[1]
 
     # Reads all the parameters from the input file
-    rB , rPP, center, X, Y, Z, symmetry, outputFile, pattern, npattern , atoms, dist, a, b, c, alpha, beta, gamma, showBath, evjen, showFrag, notInPseudo, notInFrag, symGenerator, generator, translation = read_input(inputFile)
+    rB , rPP, center, xOy, xOz, yOz, X, Y, Z, symmetry, outputFile, pattern, npattern , atoms, dist, a, b, c, alpha, beta, gamma, showBath, evjen, showFrag, notInPseudo, notInFrag, symGenerator, generator, translation = read_input(inputFile)
 
     if verbose > 0:
         out_input_param(rB , rPP, center, X, Y, Z, symmetry, outputFile, pattern, npattern , atoms, dist, a, b, c, alpha, beta, gamma, showBath, evjen, showFrag, notInPseudo, notInFrag, symGenerator, generator, translation)
@@ -79,25 +79,55 @@ if __name__=='__main__':
 
     
     # Orienting the big cell
+    if xOy != []:
+        a = find_center(xOy[0], coordinates)
+        b = a
+        w = [a]
+        while np.absolute(np.absolute(np.dot(a/np.linalg.norm(a),b/np.linalg.norm(b)))-1) < 1e-6:
+            b = find_center(xOy[1], coordinates, without=w)
+            w.append(b)
+        c = np.cross(a,b)
+        M = rotation_matrix(c, [0,0,1])
+        coordinates = rotate(M, coordinates)
+    if xOz != []:
+        a = find_center(xOz[0], coordinates)
+        b = a
+        w = [a]
+        while np.absolute(np.absolute(np.dot(a/np.linalg.norm(a),b/np.linalg.norm(b)))-1) < 1e-6:
+            b = find_center(xOz[1], coordinates, without=w)
+            w.append(b)
+        c = np.cross(a,b)
+        M = rotation_matrix(c, [0,1,0])
+        coordinates = rotate(M, coordinates)
+    if yOz != []:
+        a = find_center(yOz[0], coordinates)
+        b = a
+        w = [a]
+        while np.absolute(np.absolute(np.dot(a/np.linalg.norm(a),b/np.linalg.norm(b)))-1) < 1e-6:
+            b = find_center(yOz[1], coordinates, without=w)
+            w.append(b)
+        c = np.cross(a,b)
+        M = rotation_matrix(c, [1,0,0])
+        coordinates = rotate(M, coordinates)
     if X != []:
         k = [1,0,0]
 
         xVec = find_center(X,coordinates)
-        M = rotation_matrix(k,xVec)
+        M = rotation_matrix(xVec, k)
 
         coordinates = rotate(M, coordinates)
     if Y != []:
         k = [0,1,0]
 
         yVec = find_center(Y,coordinates)
-        M = rotation_matrix(k,yVec)
+        M = rotation_matrix(yVec, k)
 
         coordinates = rotate(M, coordinates)
     if Z != []:
         k = [0,0,1]
 
         zVec = find_center(Z,coordinates)
-        M = rotation_matrix(k,zVec)
+        M = rotation_matrix(zVec, k)
 
         coordinates = rotate(M, coordinates)
 
diff --git a/src/read_input.py b/src/read_input.py
index 0e0bb37..a8c204a 100644
--- a/src/read_input.py
+++ b/src/read_input.py
@@ -9,6 +9,9 @@ def read_input(inputFile):
     X = []
     Y = []
     Z = []
+    xOy = []
+    xOz = []
+    yOz = []
     symmetry = []
     outputFile = ""
     pattern = []
@@ -67,6 +70,15 @@ def read_input(inputFile):
         elif ls[0].casefold() == 'z_axis':
             ls.pop(0)
             Z = [i for i in ls]
+        elif ls[0].casefold() == "xoy":
+            xOy.append( f.readline().split() )
+            xOy.append( f.readline().split() )
+        elif ls[0].casefold() == "xoz":
+            xOz.append( f.readline().split() )
+            xOz.append( f.readline().split() )
+        elif ls[0].casefold() == "yoz":
+            yOz.append( f.readline().split() )
+            yOz.append( f.readline().split() )
         elif ls[0].casefold() == 'symmetry':
             ls.pop(0)
             symmetry = [i for i in ls]
@@ -192,4 +204,4 @@ def read_input(inputFile):
             print("Bad input : missing the keyword -- %s --"%t)
             sys.exit()
 
-    return rB , rPP, center, X, Y, Z, symmetry, outputFile, pattern, npattern , atoms, dist, a, b, c, alpha, beta, gamma, showBath, evjen, showFrag, notInPseudo, notInFrag, symGenerator, generator, translate
+    return rB , rPP, center, xOy, xOz, yOz, X, Y, Z, symmetry, outputFile, pattern, npattern , atoms, dist, a, b, c, alpha, beta, gamma, showBath, evjen, showFrag, notInPseudo, notInFrag, symGenerator, generator, translate
diff --git a/src/utils.py b/src/utils.py
index 512f712..870a3a9 100644
--- a/src/utils.py
+++ b/src/utils.py
@@ -119,7 +119,7 @@ def translate(v,coordinates):
 
 # Finds the point at the center of the given atoms that are the 
 # closest to the origin
-def find_center(centerList, coordinates):
+def find_center(centerList, coordinates, without=[]):
 
     centers = []
     for i in range(len(centerList)):
@@ -129,6 +129,13 @@ def find_center(centerList, coordinates):
         c.append(distance(c,[0,0,0]))  # Computing the distance to the origin
 
     for at in coordinates:
+        w = True
+        for a in without:
+            if distance(at, a) < 1e-6:
+                w = False
+                break
+        if not w:
+            continue
         if at[3] in centerList:
             centers = sorted(centers, key=operator.itemgetter(3)) # Sorting the list with respect to the distance to the origin
             d = distance(at,[0,0,0])