From c5a52818100fc8598879c30f856d8035dd76d932 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?L=C3=A9o=20Gaspard?= Date: Fri, 26 Feb 2021 15:48:23 +0100 Subject: [PATCH] All the code at once --- main.py | 163 ++++++++++++++++++++ src/out.py | 157 +++++++++++++++++++ src/read_input.py | 195 ++++++++++++++++++++++++ src/utils.py | 373 ++++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 888 insertions(+) create mode 100644 main.py create mode 100644 src/out.py create mode 100644 src/read_input.py create mode 100644 src/utils.py diff --git a/main.py b/main.py new file mode 100644 index 0000000..f833494 --- /dev/null +++ b/main.py @@ -0,0 +1,163 @@ +import sys +import numpy as np + +from src.read_input import * +from src.out import * +from src.utils import * + +if __name__=='__main__': + + verbose = 2 + + if len(sys.argv) == 1: + print("Please provide the input file") + sys.exit() + elif len(sys.argv) == 3 : + if sys.argv[2] == 's': + verbose = 0 + elif sys.argv[2] == 'v': + verbose = 1 + elif sys.argv[2] == 'vv': + verbose = 2 + elif sys.argv[2] == 'vvv': + verbose = 3 + else: + print("Unknown argument -%s-"%sys.argv[3]) + elif len(sys.argv) > 3: + print("Too much arguments, please provide an input file and a verbose level (v, vv, vvv)") + sys.exit() + + 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) + + 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) + + # Converting the angles to radian + alpha = alpha * np.pi / 180.0 + beta = beta * np.pi / 180.0 + gamma = gamma * np.pi / 180.0 + + + # Computing the number of replications needed in each directions using the interreticular distance + # for the planes 100, 010 and 001. + # + # The condition is : d_{hkl} >= 2*bath_radius + |translation_vector| + + fac = 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)) + + nA = int(np.ceil(np.sin(alpha)*(2*rB+np.linalg.norm(translation))/(a*fac)))+1 + nB = int(np.ceil(np.sin(beta)*(2*rB+np.linalg.norm(translation))/(b*fac)))+1 + nC = int(np.ceil(np.sin(gamma)*(2*rB+np.linalg.norm(translation))/(c*fac)))+1 + + if verbose > 1: + print("The big cell will be of dimensions %2ix%2ix%2i\n"%(nA,nB,nC)) + + coordinates = big_cell(generator,symGenerator,a,b,c,alpha,beta,gamma,nA,nB,nC) + + # Computing the translation vector by addition of the user translation vector and a translation + # vector that puts the origin at the center of the big cell + t = [-0.5,-0.5,-0.5] + M = get_cell_matrix(nA*a,nB*b,nC*c,alpha,beta,gamma) + t = np.matmul(M,t) + t = [t[0]+translation[0],t[1]+translation[1],t[2]+translation[2]] + + # Translating the coordinates + coordinates = translate(t, coordinates) + + # Finding the center and translating the coordinates + # If this vector creates a displacment bigger than a, b or c + # in any of the abc directions, this might result in an incomplete + # sphere later, the user should provide a translation vector + # to correct this + if center != []: + c = find_center(center,coordinates) + coordinates = translate(-c,coordinates) + + + # Orienting the big cell + if X != []: + k = [1,0,0] + + xVec = find_center(X,coordinates) + M = rotation_matrix(k,xVec) + + rotate(M, coordinates) + if Y != []: + k = [0,1,0] + + yVec = find_center(Y,coordinates) + M = rotation_matrix(k,yVec) + + rotate(M, coordinates) + if Z != []: + k = [0,0,1] + + zVec = find_center(Z,coordinates) + M = rotation_matrix(k,zVec) + + rotate(M, coordinates) + + if verbose > 2: + print("The big cell contains %5i atoms and will be printed in the file big_cell.xyz\n"%len(coordinates)) + write_coordinates(coordinates,'big_cell.xyz',3) + + # Cutting the sphere in the big cell + coordinates = cut_sphere(coordinates,rB) + + if verbose > 2: + print("The sphere contains %5i atoms and will be printed in the file sphere.xyz\n"%len(coordinates)) + write_coordinates(coordinates,'sphere.xyz',3) + + # Finding the fragment + + nAt, coordinates = find_fragment(coordinates,pattern,npattern,notInFrag) + + if verbose > 2 or showFrag: + print("The fragment contains %3i atoms and will be printed in the file fragment.xyz\n"%nAt) + write_coordinates(coordinates,'fragment.xyz',4,'O') + + coordinates = find_pseudo(coordinates,rPP,notInPseudo) + + if verbose > 2 or showBath: + print("The bath will be printed in the file bath.xyz\n") + write_coordinates(coordinates,'bath.xyz',3) + print("The bath sorted with the fragment/pseudo/charge will be printed in the file bath_coloured.xyz\n") + write_coordinates(coordinates,'bath_coloured.xyz',3,color='yes') + + if evjen: + charges = evjen_charges(coordinates,atoms) + else: + charges = [] + for i in range(len(coordinates)): + li = coordinates[i][3] + ii = np.where(atoms=li)[0] + charges.append(atoms[ii+1]) + + if verbose > 1: + print("The total charge fragment+pseudopotential+bath is : % 8.6f\n"%np.sum(charges)) + + if symmetry != []: + nuc1 = nuclear_repulsion(coordinates,charges) + if verbose > 1: + print("Nuclear repulsion before the symmetry : % 8.6f\n"%nuc1) + + coordinates,charges,indexList = compute_symmetry(coordinates,charges,symmetry) + + nuc2 = nuclear_repulsion(coordinates,charges) + if verbose > 1: + print("Nuclear repulsion after the symmetry : % 8.6f\n"%nuc2) + print("The total charge fragment+pseudopotential+bath after symmetry is : % 8.6f\n"%np.sum(charges)) + if verbose > 2: + print("The symmetrized coordinates contain %5i atoms \n"%len(indexList)) + + else: + indexList = [i for i in range(len(coordinates))] + + write_output(outputFile,coordinates,charges,indexList) + if verbose > 2: + print("The output has been written to %s \n"%outputFile) + out_interatomic_distances(coordinates) + diff --git a/src/out.py b/src/out.py new file mode 100644 index 0000000..86e86a3 --- /dev/null +++ b/src/out.py @@ -0,0 +1,157 @@ +import numpy as np +import operator + +from src.utils import distance + +def 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): + print("Output file is : %s\n"%outputFile) + print("Bath radius : % 16.10f\n"%rB) + print("Pseudo-potential width : % 16.10f\n"%rPP) + print("Lattice parameters :") + print(" a = % 8.5f"%a) + print(" b = % 8.5f"%b) + print(" c = % 8.5f"%c) + print(" alpha = % 6.5f"%alpha) + print(" beta = % 6.5f"%beta) + print(" gamma = % 6.5f\n"%gamma) + print("Group symmetry operations :") + for i in symGenerator: + print(" %10s %10s %10s"%(i[0],i[1],i[2])) + print("\nAtomic positions :") + for i in generator: + print(" %2s % 8.5f % 8.5f % 8.5f"%(i[0],i[1],i[2],i[3])) + print("\nThere %3s %2i %8s :"%('is'*(len(pattern)==1)+'are'*(len(pattern)>1),len(pattern),'pattern'+'s'*(len(pattern)>1))) + for i in range(len(pattern)): + print("The following pattern will appear %2i %5s"%(npattern[i],'time'+'s'*(npattern[i]>1))) + for j in range(0,len(pattern[i]),2): + print(" %2i %2s\n"%(pattern[i][j],pattern[i][j+1]),end='') + print("\nThere %3s %2i %5s"%('is'*(len(atoms)==1)+'are'*(len(atoms)>1), len(atoms), 'atom'+'s'*(len(atoms)>1))) + print("Label Charge Number of neighbours Maximum bond length") + for i in atoms: + print(" %2s % 5.2f %2i % 6.4f\n"%(i[0],i[1],i[2],i[3])) + if(len(center)>0): + print("Centering between the following atom(s) :") + for i in center: + print(" %2s \n"%i,end='') + print() + if(len(X) > 0): + print("Aligning the X axis with the following atom(s) :") + for i in X: + print(" %2s "%i,end='') + print() + if(len(Y) > 0): + print("Aligning the Y axis with the following atom(s) :") + for i in Y: + print(" %2s "%i,end='') + print() + if(len(Z) > 0): + print("Aligning the Z axis with the following atom(s) :") + for i in Z: + print(" %2s "%i,end='') + print() + if translation != [0.0,0.0,0.0]: + print("The following translation will be applied : % 5.3f % 5.3f % 5.3f"%(translation[0],translation[1],translation[2])) + if(len(symmetry) > 0): + print("Treating the following symmetry operations :") + for i in symmetry: + print(" %3s "%i,end='') + print() + if(evjen): + print("The program will reequilibrate the charges at the limits of the spheres using Evjen method\n") + if(showFrag): + print("The program will print the fragment coordinates in the file fragment.xyz\n") + if(showFrag): + print("The program will print the bath coordinates in the files bath.xyz and bath_coloured.xyz\n") + if len(notInPseudo) > 0: + print("The following %5s will not be considered in the pseudopotential shell :"%('atom'+'s'*(len(notInPseudo)>1))) + for i in notInPseudo: + print(" %2s "%i) + print() + if len(notInFrag) > 0: + print("The following %5s will be excluded from the fragment :"%('atom'+'s'*(len(notInFrag)>1))) + for i in notInFrag: + print(" %2s % 8.6f % 8.6f % 8.6f"%(i[0],i[1],i[2],i[3])) + print() + +# Write the coordinates to a file using the string +# at the position 'label' as atom name +# if which is specified, it must be an array of labels +# color = 'no' means that the label will be printed +# color = 'yes' means that the color in the bath will be printed +def write_coordinates(coordinates, fileName, label, which='all',color='no'): + + f = open(fileName,'w') + l = [i[label] for i in coordinates] + + if color=='no': + lab = 3 + else: + lab = 4 + + unique, count = np.unique(l,return_counts=True) + labelDic = dict(zip(unique,count)) + + if which=='all': + f.write("%i \n\n"%len(coordinates)) + for i in coordinates: + f.write("%2s % 10.6f % 10.6f % 10.6f\n"%(i[lab],i[0],i[1],i[2])) + else: + count = 0 + for la in which: + count += labelDic[la] + f.write("%i \n\n"%count) + for i in coordinates: + if i[label] in which: + f.write("%2s % 10.6f % 10.6f % 10.6f\n"%(i[lab],i[0],i[1],i[2])) + + f.close() + +# Writes the output file +def write_output(outputFile, coordinates, charges, indexList): + + f = open(outputFile,'w') + + d = {'O':'FRAGMENT','Cl':'PSEUDOPOTENTIAL','C':'BATH'} + + count = 0 + + for t in d: + f.write("%s\n"%d[t]) + + f.write(" Label x y z Charge\n") + for i in indexList: + if coordinates[i][4] == t: + count += 1 + l = coordinates[i][3]+str(count) + f.write("%6s % 16.10f % 16.10f % 16.10f % 8.5f\n"%(l,coordinates[i][0],coordinates[i][1],coordinates[i][2],charges[i])) + f.write("\n") + + +# Writes the interatomic distances inferior to 4 A +def out_interatomic_distances(coordinates): + distanceList = [] + + for i in range(len(coordinates)-1): + for j in range(i+1,len(coordinates)): + a = coordinates[i] + b = coordinates[j] + d = distance(a,b) + + accept = True + if d > 4: + accept = False + else: + for k in distanceList: + if ((a[3] == k[0] and b[3] == k[1]) or (a[3] == k[1] and b[3] == k[0])) and np.abs(d-k[2]) < 1e-5: + accept = False + + if accept: + distanceList.append([a[3],b[3],d]) + + distanceList = sorted(distanceList,key=operator.itemgetter(2)) + + print("All the interatomic distances < 5 A :") + print("Atom1 Atom2 Distance") + for i in distanceList: + print(" %2s %2s %6.4f"%(i[0],i[1],i[2])) + diff --git a/src/read_input.py b/src/read_input.py new file mode 100644 index 0000000..badc5ea --- /dev/null +++ b/src/read_input.py @@ -0,0 +1,195 @@ +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() + 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() + while line.strip().casefold() != 'end': + symGenerator.append(line.split(',')) + line = f.readline() + 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([ls[0],float(ls[1]),float(ls[2]),float(ls[3])]) + 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 new file mode 100644 index 0000000..a366955 --- /dev/null +++ b/src/utils.py @@ -0,0 +1,373 @@ +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 + r12 = r12 - np.round(r12) + r12 = np.matmul(fracToCart,r12) + d = np.sqrt(r12[0]**2+r12[1]**2+r12[2]**2) + + if(d<1e-5): + # We check if we don't already want to delete this atom + if j not in toDel: + toDel.append(j) + + # 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]: + 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 < 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 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: + 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)): + 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(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