diff --git a/crystal_mec.py b/crystal_mec.py index e4fcadb..42c7fc8 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, 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) + 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) 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,59 +79,25 @@ 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(xVec, k) + M = rotation_matrix(k,xVec) coordinates = rotate(M, coordinates) - if Y != []: k = [0,1,0] yVec = find_center(Y,coordinates) - M = rotation_matrix(yVec, k) + M = rotation_matrix(k,yVec) coordinates = rotate(M, coordinates) if Z != []: k = [0,0,1] zVec = find_center(Z,coordinates) - M = rotation_matrix(zVec, k) + M = rotation_matrix(k,zVec) coordinates = rotate(M, coordinates) diff --git a/src/__pycache__/out.cpython-38.pyc b/src/__pycache__/out.cpython-38.pyc deleted file mode 100644 index 0caae17..0000000 Binary files a/src/__pycache__/out.cpython-38.pyc and /dev/null differ diff --git a/src/__pycache__/out.cpython-39.pyc b/src/__pycache__/out.cpython-39.pyc deleted file mode 100644 index aeba510..0000000 Binary files a/src/__pycache__/out.cpython-39.pyc and /dev/null differ diff --git a/src/__pycache__/read_input.cpython-38.pyc b/src/__pycache__/read_input.cpython-38.pyc deleted file mode 100644 index af2794d..0000000 Binary files a/src/__pycache__/read_input.cpython-38.pyc and /dev/null differ diff --git a/src/__pycache__/read_input.cpython-39.pyc b/src/__pycache__/read_input.cpython-39.pyc deleted file mode 100644 index ac029eb..0000000 Binary files a/src/__pycache__/read_input.cpython-39.pyc and /dev/null differ diff --git a/src/__pycache__/utils.cpython-38.pyc b/src/__pycache__/utils.cpython-38.pyc deleted file mode 100644 index 8fe73bc..0000000 Binary files a/src/__pycache__/utils.cpython-38.pyc and /dev/null differ diff --git a/src/__pycache__/utils.cpython-39.pyc b/src/__pycache__/utils.cpython-39.pyc deleted file mode 100644 index 9300972..0000000 Binary files a/src/__pycache__/utils.cpython-39.pyc and /dev/null differ diff --git a/src/read_input.py b/src/read_input.py index 9dd2f91..0e0bb37 100644 --- a/src/read_input.py +++ b/src/read_input.py @@ -9,9 +9,6 @@ def read_input(inputFile): X = [] Y = [] Z = [] - xOy = [] - xOz = [] - yOz = [] symmetry = [] outputFile = "" pattern = [] @@ -70,27 +67,6 @@ 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] @@ -216,4 +192,4 @@ def read_input(inputFile): print("Bad input : missing the keyword -- %s --"%t) sys.exit() - 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 + 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/read_input.py~ b/src/read_input.py~ deleted file mode 100644 index 1682874..0000000 --- a/src/read_input.py~ +++ /dev/null @@ -1,196 +0,0 @@ -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 0a89a33..512f712 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, without=[]): +def find_center(centerList, coordinates): centers = [] for i in range(len(centerList)): @@ -129,14 +129,6 @@ def find_center(centerList, coordinates, without=[]): 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~ deleted file mode 100644 index 348e604..0000000 --- a/src/utils.py~ +++ /dev/null @@ -1,390 +0,0 @@ -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