diff --git a/crystal_mec.py b/crystal_mec.py
index 42c7fc8..e4fcadb 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, X, Y, Z, xOy, yOz, xOz, 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,59 @@ 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(a) ) ) - 1 ) <= 1e-6:
+            b = find_center(xOy[1], coordinates, without=w)
+            w.append(b)
+        c = np.cross(a, b)
+        k = [0,0,1]
+        M = rotation_matrix(c, k)
+        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(a) ) ) - 1 ) <= 1e-6:
+            b = find_center(xOz[1], coordinates, without=w)
+            w.append(b)
+        c = np.cross(a, b)
+        k = [0,1,0]
+        M = rotation_matrix(c, k)
+        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(a) ) ) - 1 ) <= 1e-6:
+            b = find_center(yOz[1], coordinates, without=w)
+            w.append(b)
+        c = np.cross(a, b)
+        k = [1,0,0]
+        M = rotation_matrix(c, k)
+        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/__pycache__/out.cpython-38.pyc b/src/__pycache__/out.cpython-38.pyc
new file mode 100644
index 0000000..0caae17
Binary files /dev/null and b/src/__pycache__/out.cpython-38.pyc differ
diff --git a/src/__pycache__/out.cpython-39.pyc b/src/__pycache__/out.cpython-39.pyc
new file mode 100644
index 0000000..aeba510
Binary files /dev/null and b/src/__pycache__/out.cpython-39.pyc differ
diff --git a/src/__pycache__/read_input.cpython-38.pyc b/src/__pycache__/read_input.cpython-38.pyc
new file mode 100644
index 0000000..af2794d
Binary files /dev/null and b/src/__pycache__/read_input.cpython-38.pyc differ
diff --git a/src/__pycache__/read_input.cpython-39.pyc b/src/__pycache__/read_input.cpython-39.pyc
new file mode 100644
index 0000000..ac029eb
Binary files /dev/null and b/src/__pycache__/read_input.cpython-39.pyc differ
diff --git a/src/__pycache__/utils.cpython-38.pyc b/src/__pycache__/utils.cpython-38.pyc
new file mode 100644
index 0000000..8fe73bc
Binary files /dev/null and b/src/__pycache__/utils.cpython-38.pyc differ
diff --git a/src/__pycache__/utils.cpython-39.pyc b/src/__pycache__/utils.cpython-39.pyc
new file mode 100644
index 0000000..9300972
Binary files /dev/null and b/src/__pycache__/utils.cpython-39.pyc differ
diff --git a/src/read_input.py b/src/read_input.py
index 0e0bb37..9dd2f91 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,27 @@ def read_input(inputFile):
         elif ls[0].casefold() == 'z_axis':
             ls.pop(0)
             Z = [i for i in ls]
+        elif ls[0].casefold() == "xoy":
+            line = f.readline()
+            ls = line.split()
+            xOy.append(ls)
+            line = f.readline()
+            ls = line.split()
+            xOy.append(ls)
+        elif ls[0].casefold() == "xoz":
+            line = f.readline()
+            ls = line.split()
+            xOz.append(ls)
+            line = f.readline()
+            ls = line.split()
+            xOz.append(ls)
+        elif ls[0].casefold() == "yoz":
+            line = f.readline()
+            ls = line.split()
+            yOz.append(ls)
+            line = f.readline()
+            ls = line.split()
+            yOz.append(ls)
         elif ls[0].casefold() == 'symmetry':
             ls.pop(0)
             symmetry = [i for i in ls]
@@ -192,4 +216,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, X, Y, Z, xOy, yOz, xOz, symmetry, outputFile, pattern, npattern , atoms, dist, a, b, c, alpha, beta, gamma, showBath, evjen, showFrag, notInPseudo, notInFrag, symGenerator, generator, translate
diff --git a/src/read_input.py~ b/src/read_input.py~
new file mode 100644
index 0000000..1682874
--- /dev/null
+++ b/src/read_input.py~
@@ -0,0 +1,196 @@
+import sys
+
+# Parses the input file 
+def read_input(inputFile):
+
+    rB = 0.0
+    rPP = 0.0
+    center = []
+    X = []
+    Y = []
+    Z = []
+    symmetry = []
+    outputFile = ""
+    pattern = []
+    npattern = []
+    atoms = []
+    dist = []
+    a = 0.0
+    b = 0.0
+    c = 0.0
+    alpha = 90.0
+    beta = 90.0
+    gamma = 90.0
+    showBath = False
+    evjen = False
+    showFrag = False
+    notInPseudo = []
+    notInFrag = []
+    symGenerator = []
+    generator = []
+    translate = [0.0,0.0,0.0]
+
+    checkInput = {'bath':False,'pseudo':False,'output':False,'pattern':False,'npattern':False,'a':False,'b':False,'c':False,'atoms':False,'symmetry_generator':False,'generator':False}
+
+    f = open(inputFile,'r')
+
+    line = 'x'
+
+    while line.casefold() != 'end_of_input':
+        line = f.readline().strip()
+        ls = line.split()
+        print(ls)
+        if ls == []:
+            continue
+        if ls[0].casefold() in checkInput:
+            checkInput[ls[0].casefold()] = True
+        if ls[0].casefold() == 'bath':
+            try:
+                rB = float(ls[1])
+            except ValueError:
+                print("Error while parsing the input file : %s is not a valid bath radius value"%(ls[1]))
+                sys.exit()
+        elif ls[0].casefold() == 'pseudo':
+            try:
+                rPP = float(ls[1])
+            except ValueError:
+                print("Error while parsing the input file : %s is not a valid pseudopotential radius value"%(ls[1]))
+                sys.exit()
+        elif ls[0].casefold() == 'center':
+            ls.pop(0)
+            center = [i for i in ls]
+        elif ls[0].casefold() == 'x_axis':
+            ls.pop(0)
+            X = [i for i in ls]
+        elif ls[0].casefold() == 'y_axis':
+            ls.pop(0)
+            Y = [i for i in ls]
+        elif ls[0].casefold() == 'z_axis':
+            ls.pop(0)
+            Z = [i for i in ls]
+        elif ls[0].casefold() == 'symmetry':
+            ls.pop(0)
+            symmetry = [i for i in ls]
+        elif ls[0].casefold() ==  'output':
+            outputFile = ls[1]
+        elif ls[0].casefold() == 'pattern':
+            line = f.readline()
+            while line.strip().casefold() != 'end':
+                pattern.append([])
+                ls = line.split()
+                for i in range(len(ls)):
+                    if i%2 == 0:
+                        pattern[-1].append(int(ls[i]))
+                    else:
+                        pattern[-1].append(ls[i])
+                line = f.readline()
+        elif ls[0].casefold() == 'npattern':
+            ls.pop(0)
+            try:
+                npattern = [int(i) for i in ls]
+            except ValueError:
+                print("Error while parsing the input file : the number of patterns is not valid %s"%(line))
+                sys.exit()
+        elif ls[0].casefold() == 'lattice':
+            line = f.readline()
+            while line.strip().casefold() != 'end':
+                ls = line.split()
+                if ls[0].casefold() in checkInput:
+                    checkInput[ls[0].casefold()] = True
+                if ls[0].casefold() == 'a':
+                    try:
+                        a = float(ls[1])
+                    except ValueError:
+                        print("Error while parsing the input file : bad value for the lattice parameter %s"%ls[1])
+                        sys.exit()
+                elif ls[0].casefold() == 'b':
+                    try:
+                        b = float(ls[1])
+                    except ValueError:
+                        print("Error while parsing the input file : bad value for the lattice parameter %s"%ls[1])
+                        sys.exit()
+                elif ls[0].casefold() == 'c':
+                    try:
+                        c = float(ls[1])
+                    except ValueError:
+                        print("Error while parsing the input file : bad value for the lattice parameter %s"%ls[1])
+                        sys.exit()
+                elif ls[0].casefold() == 'alpha':
+                    try:
+                        alpha = float(ls[1])
+                    except ValueError:
+                        print("Error while parsing the input file : bad value for the lattice parameter %s"%ls[1])
+                        sys.exit()
+                elif ls[0].casefold() == 'beta':
+                    try:
+                        beta = float(ls[1])
+                    except ValueError:
+                        print("Error while parsing the input file : bad value for the lattice parameter %s"%ls[1])
+                        sys.exit()
+                elif ls[0].casefold() == 'gamma':
+                    try:
+                        gamma = float(ls[1])
+                    except ValueError:
+                        print("Error while parsing the input file : bad value for the lattice parameter %s"%ls[1])
+                        sys.exit()
+                line = f.readline()
+        elif ls[0].casefold() == 'atoms':
+            line = f.readline()
+            while line.strip().casefold() != 'end':
+                ls = line.split()
+                if(len(ls)) != 4:
+                        print("Error while parsing the input file : not enough values given for the atom in line %s"%line)
+                        sys.exit()
+                try:
+                    atoms.append([ls[0], float(ls[1]), int(ls[2]), float(ls[3])])
+                except ValueError:
+                    print("Error while parsing the input file : bad value for the atom %s"%line)
+                line = f.readline()
+        elif ls[0].casefold() == 'show_bath':
+            showBath = True
+        elif ls[0].casefold() == 'translate':
+            ls.pop(0)
+            try:
+                translate = [float(ls[0]),float(ls[1]),float(ls[2])]
+            except ValueError:
+                print("Error while parsing the input file : the translation vector is not valid %s"%line)
+                sys.exit()
+        elif ls[0].casefold() == 'not_in_pseudo':
+            ls.pop(0)
+            notInPseudo = [i for i in ls]
+        elif ls[0].casefold() == 'show_frag':
+            showFrag = True
+        elif ls[0].casefold() == 'evjen':
+            evjen = True
+        elif ls[0].casefold() == 'symmetry_generator':
+            line = f.readline().replace("'","")
+            while line.strip().casefold() != 'end':
+                symGenerator.append(line.split(','))
+                line = f.readline().replace("'","")
+        elif ls[0].casefold() == 'generator':
+            line = f.readline()
+            while line.strip().casefold() != 'end':
+                ls = line.split()
+                try:
+                    generator.append([ls[0],float(ls[1]),float(ls[2]),float(ls[3])])
+                except ValueError:
+                    print("Error while parsing the input file : bad value for the generator atom %s"%line)
+                    sys.exit()
+                line = f.readline()
+        elif ls[0].casefold() == 'not_in_frag':
+            line = f.readline()
+            while line.strip().casefold() != 'end':
+                ls = line.split()
+                try:
+                    notInFrag.append([float(ls[0]),float(ls[1]),float(ls[2])])
+                except ValueError:
+                    print("Error while parsing the input file : bad value for the atom %s"%line)
+                    sys.exit()
+                line = f.readline()
+    f.close()
+    for t in checkInput:
+        if checkInput[t] == False:
+            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
diff --git a/src/utils.py b/src/utils.py
index 512f712..0a89a33 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,14 @@ 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 i in without:
+            d = distance(at, i)
+            if d < 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])
diff --git a/src/utils.py~ b/src/utils.py~
new file mode 100644
index 0000000..348e604
--- /dev/null
+++ b/src/utils.py~
@@ -0,0 +1,390 @@
+import numpy as np
+import operator
+import sys
+
+def distance(a,b):
+    # Returns the 3D distance between a and b where 
+    # a and b are array where x, y and z are at the
+    # position 0, 1 and 2
+
+    x = a[0]-b[0]
+    y = a[1]-b[1]
+    z = a[2]-b[2]
+
+    return np.sqrt(x**2+y**2+z**2)
+
+def get_cell_matrix(a,b,c,alpha,beta,gamma):
+    # Computing the volume of the primitive cell
+    omega = a*b*c*np.sqrt(1-np.cos(alpha)**2-np.cos(beta)**2-np.cos(gamma)**2+2*np.cos(alpha)*np.cos(beta)*np.cos(gamma))
+
+    # Computing the matrix 
+    M = [
+            [a,b*np.cos(gamma),c*np.cos(beta)],
+            [0,b*np.sin(gamma),c*(np.cos(alpha)-np.cos(beta)*np.cos(gamma))/(np.sin(gamma))],
+            [0,0,omega/(a*b*np.sin(gamma))]
+        ]
+    return M
+
+def big_cell(generator,symGenerator,a,b,c,alpha,beta,gamma,nA,nB,nC):
+    coords = []
+    
+    # Computing the matrix converting fractional to cartesian
+    fracToCart = get_cell_matrix(a,b,c,alpha,beta,gamma)
+
+    for gen in generator:
+        x = gen[1]
+        y = gen[2]
+        z = gen[3]
+
+        for sym in symGenerator:
+            u = eval(sym[0])
+            v = eval(sym[1])
+            w = eval(sym[2])
+
+            # Making sure the value is within the range [0,1]
+            u = u + 1*(u<0) - 1*(u>1)
+            v = v + 1*(v<0) - 1*(v>1)
+            w = w + 1*(w<0) - 1*(w>1)
+            
+            coords.append([u,v,w,gen[0]])
+
+    # Deleting the redundant atoms
+    toDel = []
+    for i in range(len(coords)-1):
+            for j in range(i+1,len(coords)):
+                # Computing the distance using the minimum image convention
+                # as described in Appendix B equation 9 of 
+                # "Statistical Mechanics : Theory and Molecular Simulations
+                # Mark E. Tuckerman"
+
+                r1 = np.array(coords[i][:3])
+                r2 = np.array(coords[j][:3])
+                r12 = r1-r2
+                da = np.sqrt(r12[0]**2+r12[1]**2+r12[2]**2)
+                r12 = r12 - np.round(r12)
+                db = da - np.sqrt(r12[0]**2+r12[1]**2+r12[2]**2)
+                r12 = np.matmul(fracToCart,r12)
+                d = np.sqrt(r12[0]**2+r12[1]**2+r12[2]**2)
+
+                if(d<1e-2):
+                    # We check if we don't already want to delete this atom
+                    if j not in toDel:
+                        toDel.append(j)
+
+    toDel = sorted(toDel)
+
+    # We delete the atoms in the list
+    for i in range(len(toDel)):
+        coords.pop(toDel[i]-i)
+
+    newCoords = []
+
+    # We replicate the cell nA, nB, nC times 
+    for at in coords:
+        newCoords.append([at[0],at[1],at[2],at[3]])
+        for a in range(1,nA):
+            newCoords.append([at[0]+a,at[1],at[2],at[3]])
+            for b in range(1,nB):
+                newCoords.append([at[0]+a,at[1]+b,at[2],at[3]])
+                for c in range(1,nC):
+                    newCoords.append([at[0]+a,at[1]+b,at[2]+c,at[3]])
+            for c in range(1,nC):
+                newCoords.append([at[0]+a,at[1],at[2]+c,at[3]])
+        for b in range(1,nB):
+            newCoords.append([at[0],at[1]+b,at[2],at[3]])
+            for c in range(1,nC):
+                newCoords.append([at[0],at[1]+b,at[2]+c,at[3]])
+        for c in range(1,nC):
+            newCoords.append([at[0],at[1],at[2]+c,at[3]])
+
+    # Now we convert the fractionnal coordinates to cartesian coordinates
+    coords = []
+
+    for at in newCoords:
+        r = [at[0],at[1],at[2]]
+        rxyz = np.matmul(fracToCart,r)
+        coords.append([rxyz[0],rxyz[1],rxyz[2],at[3],'C'])
+
+
+    # Returns the list of the atoms [x,y,z,label,second_label]
+    return coords
+
+# Translates all the coordinates with the vector v
+def translate(v,coordinates):
+    for c in coordinates:
+        c[0] += v[0]
+        c[1] += v[1]
+        c[2] += v[2]
+    return coordinates
+
+# Finds the point at the center of the given atoms that are the 
+# closest to the origin
+def find_center(centerList, coordinates):
+
+    centers = []
+    for i in range(len(centerList)):
+        centers.append([100,100,100])  # Setting a large value for each center
+
+    for c in centers:
+        c.append(distance(c,[0,0,0]))  # Computing the distance to the origin
+
+    for at in coordinates:
+        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])
+            if d <= centers[-1][-1] and d > 0.0:
+                centers[-1] = [at[0],at[1],at[2],d]
+
+    center = np.mean(centers,axis=0)[:3] # Computing the barycenter
+
+    return center
+
+# Defines a rotation matrix that will put r1 at the position r2
+def rotation_matrix(r1,r2):
+
+    r1 = np.array(r1)/np.linalg.norm(r1)
+    r2 = np.array(r2)/np.linalg.norm(r2)
+
+    # Computing the cross product which is the vector around which
+    # the rotation is done
+    crossProduct = np.cross(r1,r2)
+    crossProduct = crossProduct/np.linalg.norm(crossProduct)
+
+    # Computing the angle of the rotation
+    a = np.arccos(np.dot(r1,r2))
+
+    c = np.cos(a)
+    s = np.sin(a)
+    x = crossProduct[0]
+    y = crossProduct[1]
+    z = crossProduct[2]
+
+    M = [
+            [x**2*(1-c)+c,x*y*(1-c)-z*s,x*z*(1-c)+y*s],
+            [x*y*(1-c)+z*s,y**2*(1-c)+c,y*z*(1-c)-x*s],
+            [x*z*(1-c)-y*s,y*z*(1-c)+x*s,z**2*(1-c)+c]
+            ]
+
+    return M
+
+# Rotates all the coordinates using the rotation matric M
+def rotate(M,coordinates):
+    for i in range(len(coordinates)):
+        r = [coordinates[i][0],coordinates[i][1],coordinates[i][2]]
+        rV = np.matmul(M,r)
+        coordinates[i][0] = rV[0]
+        coordinates[i][1] = rV[1]
+        coordinates[i][2] = rV[2]
+
+    return coordinates
+
+# Cuts a sphere centered on the origin in the coordinates
+def cut_sphere(coordinates,r):
+    sphere = []
+    for i in range(len(coordinates)):
+        if distance(coordinates[i],[0,0,0]) <= r:
+            sphere.append(coordinates[i])
+
+    return sphere
+
+# Finds the fragment in the coordinates
+def find_fragment(coordinates, patterns, npatterns,notInFrag):
+
+    inFrag = []
+
+    for n in range(len(patterns)):
+        pattern = patterns[n]
+        npattern = npatterns[n]
+        for i in range(npattern):
+            c = [100,100,100]
+            dc = distance([0,0,0],c)
+
+            inPattern = []
+            # Finding the closest atom of the first type in the pattern
+            for at in coordinates:
+                if at[3] == pattern[1]:
+                    d = distance([0,0,0],at)
+                    if d > 10:
+                        break
+                    if d < dc :
+                        accept = True
+                        for exc in notInFrag:
+                            d = distance(exc,at)
+                            if d < 1e-5:
+                                accept = False
+                        if accept and coordinates.index(at) not in inFrag:
+                            c = [at[0],at[1],at[2],0.0, coordinates.index(at)]
+                            dc = distance([0,0,0],c)
+            # Finding the rest of the pattern around the atom previously found
+            atIn = []
+            for j in range(0,len(pattern),2):
+                d = distance(c,[100,100,100])
+                # Initializing the atoms 
+                for k in range(pattern[j]):
+                    atIn.append([100,100,100,d])
+
+                for at in coordinates:
+                    if distance(at,[0,0,0]) > 10:
+                        break
+                    if at[3] == pattern[j+1]:
+                        atIn = sorted(atIn,key=operator.itemgetter(3))
+                        d = distance(at,c)
+                        trial = [at[0],at[1],at[2],d,coordinates.index(at)]
+                        if d < atIn[-1][3] and trial not in atIn:
+                            accept = True
+                            for exc in notInFrag:
+                                d = distance(exc,trial)
+                                if d < 1e-5:
+                                    accept = False
+                            if accept:
+                                atIn[-1] = trial
+            for at in atIn:
+                inPattern.append(at[4])
+
+            for at in inPattern:
+                if at not in inFrag:
+                    inFrag.append(at)
+
+    for at in inFrag:
+        coordinates[at][4] = 'O'
+
+    return len(inFrag), coordinates
+
+
+# Finds the pseudopotential layer around
+# the fragment
+def find_pseudo(coordinates, rPP, notInPseudo):
+
+    for at in coordinates:
+        if at[4] != 'O':
+            continue
+        for i in range(len(coordinates)):
+            if coordinates[i][4] != 'C':
+                continue
+            d = distance(at,coordinates[i])
+            if d < rPP:
+                coordinates[i][4] = 'Cl'
+
+    return coordinates
+
+# Creates lists containing the neighbours of each
+# atom
+def find_neighbours(coordinates, atoms):
+    neighbourList = [[] for i in range(len(coordinates))]
+
+    atoms = np.array(atoms).flatten()
+
+    for i in range(len(coordinates)-1):
+        for j in range(i+1,len(coordinates)):
+            li = coordinates[i][3] # Label of the atom i
+            lj = coordinates[j][3] # Label of the atom j
+
+            ii = np.where(atoms==li)[0]
+            jj = np.where(atoms==lj)[0]
+
+            ci = float(atoms[ii+1]) # Charge of the atom i
+            cj = float(atoms[jj+1]) # Charge of the atom j
+
+            if ci*cj < 0: # Checking if the charges have opposite signs
+                d = distance(coordinates[i],coordinates[j])
+
+                if d < float(atoms[ii+3]) and d < float(atoms[jj+3]):
+                    neighbourList[i].append(j)
+                    neighbourList[j].append(i)
+    return neighbourList
+
+# For each atom, finds if it has the correct number of neighbours,
+# if not, modify its charge 
+def evjen_charges(coordinates,atoms):
+    neighbourList = find_neighbours(coordinates,atoms)
+
+    atoms = np.array(atoms).flatten()
+
+    charges = []
+
+    for i in range(len(coordinates)):
+        li = coordinates[i][3]
+        ii = np.where(atoms==li)[0]
+
+        nr = len(neighbourList[i])
+        nt = int(atoms[ii+2])
+        ci = float(atoms[ii+1])
+        
+        if nr > nt:
+            print("Error : too much neighbours for atom n°%i, count %i neighbours where it should have a maximum of %i"%(i,nr,nt))
+            sys.exit()
+        charges.append(ci*nr/nt)
+
+    return charges
+
+# Computes the nuclear repulsion 
+def nuclear_repulsion(coordinates,charges):
+
+    rep = 0.0
+    
+    for i in range(len(coordinates)-1):
+        for j in range(i+1,len(coordinates)):
+            rij = distance(coordinates[i],coordinates[j])
+            ci = charges[i]
+            cj = charges[j]
+
+            if(rij < 1):
+                print(i,j,"\n",coordinates[i],"\n",coordinates[j],"\n",rij)
+
+            rep += (ci*cj)/rij
+    return rep
+
+# Computes the symmetry in the whole system
+def compute_symmetry(coordinates,charges,symmetry):
+    symmetrizedCoordinates = []
+    symmetrizedCharges = []
+    uniqueIndexList = [] # The list containing the indexes of the unique atoms
+
+    treated = [] # Will store the index of the atoms already treated
+
+    symOp = []
+
+    # Storing all the symmetry operations
+    for s in symmetry:
+        if s == 'C2x':
+            symOp.append(np.array([1,-1,-1]))
+        elif s == 'C2y':
+            symOp.append(np.array([-1,1,-1]))
+        elif s == 'C2z':
+            symOp.append(np.array([-1,-1,1]))
+        elif s == 'xOy':
+            symOp.append(np.array([1,1,-1]))
+        elif s == 'xOz':
+            symOp.append(np.array([1,-1,1]))
+        elif s == 'yOz':
+            symOp.append(np.array([-1,1,1]))
+        elif s == 'i':
+            symOp.append(np.array([-1,-1,-1]))
+
+    for i in range(len(coordinates)):
+        print(i)
+        if i in treated:
+            continue
+        
+        treated.append(i)
+        at1 = np.array(coordinates[i][:3])
+        symmetrizedCoordinates.append(coordinates[i])
+        symmetrizedCharges.append(charges[i])
+        uniqueIndexList.append(len(symmetrizedCoordinates)-1)
+
+        for j in range(len(coordinates)):
+            if j in treated or coordinates[i][3] != coordinates[j][3]:
+                continue
+
+            at2 = np.array(coordinates[j][:3])
+
+            for s in symOp:
+                if distance(at2, at1*s) < 5:
+                if distance(at1,at1*s) > 1e-4 and distance(at2,at1*s) < 1e-4: # Checking if op.at1 != at1 and that op.at2 = at1
+                    p = at1*s
+                    treated.append(j)
+                    symmetrizedCoordinates.append([p[0],p[1],p[2],coordinates[i][3],coordinates[i][4]])
+                    symmetrizedCharges.append(charges[i])
+                    break
+
+    return symmetrizedCoordinates,symmetrizedCharges,uniqueIndexList