from ase.io import write, read import numpy as np import sys import os from ase.visualize import view from ase import Atom from ase.constraints import FixAtoms import random as rng import math # Seed rng rng.seed(420697418008) # Control panel Ni_conc = [0.02, 0.05, 0.08, 0.1, 0.12, 0.15, 0.2] # List for concentrations with_print = False # Text output? size = 6 # Size of cell compared to the 1*3 cell normally use cluster_size = 3 # Size of clusters, used to calculate how many seeds to remove times = 130#00 # Number of times to run the surface constructor for each molar fraction ViewStructures = False # View structures generated? # Initiating variables #cluster_matrix = np.zeros((1,3)) lenRows = 12.791/4 countArray = np.zeros((len(Ni_conc),times)) if with_print == True: print('\nThis script creates random surface configurations for RuNiO2 and counts the number of RuNi bridge sites \nWritten by Adrian Frandsen\n\n') atoms = read('../Ni/OOO.traj') del atoms[71:] atoms *= (size,size,1) atoms.set_positions(atoms.get_positions() + [0.5*lenRows, 0.775, 0], apply_constraint = False) cellX = atoms.cell[0][0] cellY = atoms.cell[1][1] for i in range(len(Ni_conc)): if True:#with_print == True: print('-------- Ni concentration is: ' + str(Ni_conc[i]) + ' --------\n') for j in range(times): # Do this enough times to get a good sample size, 10? 20? #if with_print == True: print('Round ' + str(j+1) + '\n' ) # Place Ni for k in [atom for atom in atoms if atom.symbol == 'Ru']: if rng.random() < Ni_conc[i]: k.symbol = 'Ni' numberofseeds = len([atom.index for atom in atoms if atom.symbol == 'Ni']) # Number of cluster seeds if with_print == True: print('Number of initial seeds placed: ' + str(numberofseeds)) seedsneeded = 24*(size**2)*Ni_conc[i]/cluster_size # The number of seeds that would be expected following the concentration if with_print == True: print('Allowed number of seeds: ' + str(24*size**2) + ' x ' + str(Ni_conc[i]) + ' / (cluster size of '+ str(cluster_size) +')' + ' = ' + str(seedsneeded)) seedstoremove = round(numberofseeds - seedsneeded) # The number of cluster seeds to remove if with_print == True: print('\nRemoving ' + str(seedstoremove) + ' cluster seeds') for k in range(seedstoremove): temp = rng.randint(0,len([atom for atom in atoms if atom.symbol == 'Ni'])-1) [atom for atom in atoms if atom.symbol == 'Ni'][temp].symbol = 'Ru' if with_print == True: print('There are now ' + str(len([atom for atom in atoms if atom.symbol == 'Ni'])) + ' cluster seeds left.\n\nAfter clustering this constitutes a Ni fraction of : ' + str(3*len([atom for atom in atoms if atom.symbol == 'Ni'])/(24*size*size))) print('Compared to the ideal: ' + str(Ni_conc[i])+ '\n') # Create clusters # -Identify clustering spots for each # -- Choose clustering spot clusters = [atom for atom in atoms if atom.symbol == 'Ni'] if with_print == True: print('The chosen seeds are:', [atom.index for atom in atoms if atom.symbol == 'Ni'], '\nClustering happens in random order.\n\n') rng.shuffle(clusters) for spot in clusters: # Clustering spot-list cspots = [atom.index for atom in atoms if atom.symbol == 'Ru' # Is Ru AND in one of two cases AND fulfills y coordinate requirements and ( ( # Within same layer, diagonal (case 1), x values ((spot.position[0]+0.6*lenRows)%cellX < atom.position[0] < (spot.position[0]+1.4*lenRows)%cellX or (spot.position[0]-0.6*lenRows)%cellX > atom.position[0] > (spot.position[0]-1.5*lenRows)%cellX ) # z values, no wraparound, no modulus and spot.position[2]-0.2 < atom.position[2] < spot.position[2]+0.2 ) or #Within layers above and below, diagonal (case 2) ( # x values spot.position[0]-0.2 < atom.position[0] < spot.position[0]+0.2 and # z values (spot.position[2]-3.5 < atom.position[2] < spot.position[2]-3.1 or spot.position[2]+3.1 < atom.position[2] < spot.position[2]+3.5) ) ) and # y values, these are the same no matter if we go along x or z ((spot.position[1]+1.3) %cellY < atom.position[1] < (spot.position[1] + 1.7)%cellY or (spot.position[1]-1.3) %cellY > atom.position[1] > (spot.position[1] - 1.7)%cellY ) ] # Choose two if with_print == True: print('Atom number ' + str(spot.index) + ' sees the following potential clustering spots:' ) print(cspots) # Choose one nmb1 = rng.randint(0,len(cspots)-1) # Exchange 1 atoms[cspots[nmb1]].symbol = 'Ni' # Find neighbors neighbors = [atom.index for atom in atoms if ( (((atoms[cspots[nmb1]].position[1] + 3)%cellY < atom.position[1] < (atoms[cspots[nmb1]].position[1] + 3.2)%cellY) or ((atoms[cspots[nmb1]].position[1] - 3)%cellY > atom.position[1] > (atoms[cspots[nmb1]].position[1] - 3.2)%cellY)) and ((atoms[cspots[nmb1]].position[2] - 0.2) < atom.position[2] < (atoms[cspots[nmb1]].position[2] + 0.2)) and ((atoms[cspots[nmb1]].position[0] - 0.2) < atom.position[0] < (atoms[cspots[nmb1]].position[0] + 0.2)) )] # Choose second nmb2 = rng.randint(0,len(cspots)-1) # Ask if it is the same number, or if it is considered a neighbor. If any of these are true, reroll # If you want to go to high concentrations you should consider writing an if statement asking whether there is even any available Ru atoms to convert temp = 0 while nmb1 == nmb2 or atoms[cspots[nmb2]].index in neighbors: nmb2 = rng.randint(0,len(cspots)-1) temp = temp + 1 if temp > 9: print('Exiting "while True" loop') break # Exchange 2 atoms[cspots[nmb2]].symbol = 'Ni' if with_print == True: print('Chooses ' + str(cspots[nmb1]) + ' and ' + str(cspots[nmb2]) + '.\n') # Count surface configurations: Am I missing wraparound? count = 0 for l in range(size): for m in range(3*size): temp = [atom.symbol for atom in atoms if l*cellX / size < atom.position[0] < (l+0.5)*cellX / size and ( m*cellY / (3*size) < atom.position[1] < (m+2)*cellY / (3*size) or (m+1)*cellY / (3*size) % cellY < atom.position[1] < (m+2)*cellY/(3*size) % cellY ) and 24.6 < atom.position[2] < 24.8 ] if temp[0] != temp[1]: count = count + 1 # Save sitedensity = round(count/(3*size*size),4) countArray[i,j] = sitedensity if ViewStructures == True and (j == 0 or j == 1 or j == 2): view(atoms) if with_print == True: print('Number of bridge sites: ', 3*size*size) print('Number of LOER active sites: ', count) print('Fraction of LOERsites/bridgesites: ', sitedensity,'\n') # Clean for atom in [atom for atom in atoms if atom.symbol == 'Ni']: atom.symbol = 'Ru' # Save the fraction of LOER to total bridge sites #print('Number of sites (total): ', times*27) #print('Number of LOER active sites: ', count) #print('Fraction LOERsites/Totalsites: ', sitedensity,'\n\n') print('Array of site densities:\n') print(countArray) np.save('3clusterdata.npy', countArray) # Model 2: Ni along (121) planes sizeList = [3,4,5,6,7,8,9,12,20] planeData = np.zeros((len(sizeList), 2)) for i in range(len(sizeList)): atoms = read('../Ni/OOO.traj') del atoms[71:] atoms *= (sizeList[i],sizeList[i],1) cellX = atoms.cell[0][0] cellY = atoms.cell[1][1] atoms.set_positions(atoms.get_positions() + [0.5*lenRows, 0.775, 0], apply_constraint = False) # Define two planes from points normal121 = np.array((-0.5, 0.35355339, 1.06066017)) # What points? upperPoint = np.array((1.63, 2.33, 25.737)) lowerPoint = np.array((1.63, 2.33, 23.737)) # Find d for ax+by+cz+d=0: a = normal121[0] b = normal121[1] c = normal121[2] d_upper = -(a*upperPoint[0] + b*upperPoint[1] + c*upperPoint[2]) d_lower = -(a*lowerPoint[0] + b*lowerPoint[1] + c*lowerPoint[2]) # Find z = - (ax + by + d) / c for each atom and ask if it is below upper and above lower NiPlane = [atom for atom in atoms if atom.symbol == "Ru" and -(a*atom.position[0] + b*atom.position[1] + d_lower)/c < atom.position[2] < - (a*atom.position[0] + b*atom.position[1] + d_upper)/c ] for atom in NiPlane: atom.symbol = "Ni" #if ViewStructures == True: # view(atoms) # Find Ni fraction print("\n------------ Ni 121 plane (size: "+str(sizeList[i])+")-------------", "\n\nTotal Ni fraction is " + str(len(NiPlane)/(sizeList[i]*sizeList[i]*24))) Ni_top_layer = [atom for atom in atoms if atom.symbol == "Ni" and atom.position[2] > 24] Ni_count = len(Ni_top_layer) print('Ni count for top layer: ' + str(Ni_count)) print('Ni fraction for top layer: ' + str(Ni_count/(sizeList[i]*sizeList[i]*6))) print("Ni introduced for 121 plane: " +str(len(NiPlane))) # Count surface active sites count = 0 for l in range(sizeList[i]): for m in range(sizeList[i]*3): temp = [atom.symbol for atom in atoms if l*cellX / sizeList[i] < atom.position[0] < (l+0.5)*cellX / sizeList[i] and ( m*cellY / (sizeList[i]*3) < atom.position[1] < (m+2)*cellY / (sizeList[i]*3) or (m+1)*cellY / (sizeList[i]*3) % cellY < atom.position[1] < (m+2)*cellY/(sizeList[i]*3) % cellY ) and 24.6 < atom.position[2] < 24.8 ] if temp[0] != temp[1]: count = count + 1 print('Number of sites (total): ', sizeList[i]*sizeList[i]*3) print('Number of LOER active sites: ', count) print('Fraction LOERsites/Totalsites: ', round(count/(sizeList[i]*sizeList[i]*3), 4),'\n\n') planeData[i,0] = Ni_count/(sizeList[i]*sizeList[i]*6) # Nickel fraction for top layer (representative layer) planeData[i,1] = round(count/(sizeList[i]*sizeList[i]*3), 4) # active site density. np.save('planeData.npy', planeData)