import os import numpy as np import re import glob import pandas as pd import linecache import sys import gc import atexit def read_data(): length_ions = {} data_ions = {} inc_angles = [] yields = [] sim_E = [] masses =[] sweepdata = [] # contains all SDtrimSweep data, organised [Angle] main_input = [] # contains all standard simulation parameters repo_info = [] # contains strings with info regarding sputteryield sweep file and emission list sweep file spraypath = os.getcwd() # load inp. file inpfile = "SPRAY.inp" ifile = open(inpfile, "r+") ifile_lines = ifile.readlines() for line in ifile_lines: if line[0] == "#": # comment skip continue else: splitline = re.split(r'\t+', line) # print(splitline) if splitline[0] == "pi": main_input.append(int(splitline[1])) if splitline[0] == "ps": main_input.append(int(splitline[1])) if splitline[0] == "ss": main_input.append(int(splitline[1])) if splitline[0] == "np": main_input.append(int(splitline[1])) if splitline[0] == "sa": main_input.append(int(splitline[1])) if splitline[0] == "ba": main_input.append(int(splitline[1])) if splitline[0] == "ia": main_input.append(int(splitline[1])) if splitline[0] == "et": splitelem = re.split(r',+', splitline[1]) main_input.append(splitelem) if splitline[0] == "di": main_input.append(splitline[1]) if splitline[0] == "as": main_input.append(float(splitline[1])) if splitline[0] == "ri": main_input.append(int(splitline[1])) if splitline[0] == "sf": main_input.append(splitline[1]) if splitline[0] == "aa": main_input.append(splitline[1]) if splitline[0] == "ds": main_input.append(splitline[1]) element_array = main_input[7] angles = np.arange(main_input[4], main_input[5] + 1, main_input[6]) # to be simulated print("\nGeneral Simulation parameters") print("Number of primary ions:\t\t\t\t\t", main_input[0]) print("Number of sputtered atoms (virtually):\t\t\t", main_input[1]) print("Number of secondary sputtered atoms (virtually):\t", main_input[2]) print("Number of reflected ions:\t", main_input[11]) print("List of ion incidence angles respective surface normal:\t\t\t\t", angles) print("Azimuthal angle mode:\t", main_input[12]) print("Number of processors for parallel computation:\t\t", main_input[3]) print("List of elements (first is ion):\t\t\t", element_array) print("path of repository data:\t", main_input[8]) # -> needs to be down below .../Data/ level print("central lateral surface area selected for irradiation :\t", main_input[9]) print("type of the surface file:\t", main_input[10], "\n") print("display setting:\t", main_input[13], "\n") # declare the repository data path: path = "" if main_input[8] == "loc": path = os.getcwd() + "/Data" else: path = main_input[8] + "/Data" # file containing the data of the flat sputter Yields for angels 0-89 degrees (1 degree steps) # at 50-2000 eV (50 eV steps) #UPDATE V3_beta: doesn't have to be 1 degree steps anymore; a delta_alpha is calculated in the main program os.chdir(path) sweepfile = [] list_of_files_data = glob.glob("*.dat") for filename in list_of_files_data: if filename[-20:-4] == "AngleEnergySweep": sweepfile = filename repo_info.append(sweepfile) data = pd.read_csv(sweepfile, sep="\t") if len(sweepfile) == 0: print("\n\n\n\n\nWarning, no AngleEnergySweep file found. Check path and name of your sweep file\n sweep file name needs to end with _AngleEnergySweep.dat\n\n") # saving the values of the flat sputter yield for 2 keV to the variable yields (inc_angles contain the # correcponding angles 0-89 degrees) #Generates a list of files generated by either SDTrimSP or Tri3Dyn. #Finds the angles for which the simulations were conducted. #Angles have to be an integer that is contained in the filename! #Also reads the first 8 characters of the first input file to check whether it was created by SDTrimSP or Tri3Dyn regex = re.compile(r'\d+') angles = [] # list_of_files = glob.glob("Data/Ions/*.dat") list_of_files = glob.glob("Ions/*.dat") #finds integers in the names of the data files and interprets them as angles for filename in list_of_files: for x in regex.findall(filename): angles.append(int(x)) angles.sort() with open(list_of_files[0]) as testfile: simulation_info = testfile.read(8) #read first 8 characters of first data file to find out what BCA code was used os.chdir(spraypath) table = path + "/table1.txt" #table for atomic masses; adapted from standard SDTrimSP table elemental_information = pd.read_csv(table, sep="\s+", usecols=['mass'], skiprows=10) #Checks whether the input files were generated by SDTrimSP or Tri3Dyn and reads them in in the corresponding format. #Nota bene: As of now, SDTrimSP files contain a header, whereas Tri3Dyn files do not. In the future the if-conditional #could be modified in a way to explicitly check also the origin of Tri3Dyn files, should they also contain a header. #Files from SDTrimSP are imported using pandas dataframes, such that it does not matter whether or not they have been #compressed to contain only relevant information; using the names of the columns, also original format files should be #read correctly. if 'SDTrimSP' in simulation_info: os.chdir(path) print("SDTrimSP repository data found. Loading...") repo_info.append("SDtrimSP") # finds array of composition of the target, which is included in the second line teststring = linecache.getline(list_of_files[0], 2) os.chdir(spraypath) beginn = teststring.find("[") ende = teststring.find("]") #split string of composition into a python list element_array = teststring[beginn+1:ende].split(", ") #get rid of quotation marks, such that only the element symbols remain element_array[:] = [x[2:-2] for x in element_array] #write elemental masses for element in element_array: masses.append(elemental_information.loc[element, 'mass']) #Initialisation of Arrays of dictionaries/arrays data_atoms = [] length_atoms = [] #append as many dictionaries as necessary for i in range(len(element_array)-1): data_atoms.append({}) length_atoms.append({}) yields.append([]) #inc_angles.append([]) for i in range(0, len(data), 1): #go through each line of the data file if data.loc[data.index[i], "Angle"] == 0.0: #interested in the lines where angle = 0.0 sim_E.append(float(data.loc[data.index[i], "Energy"])) # sim_E should contain all simulated Energies #get the yields for a flat surface for all the components of the target for i in range(0, len(data), 1): #go through each line of the data file if data.loc[data.index[i], "Energy"] == max(sim_E): #interested in the lines where E=max inc_angles.append(float(data.loc[data.index[i], "Angle"])) for j in range(len(element_array)-1): #inc_angles[j].append(float(data.loc[data.index[i], "Angle"])) yields[j].append(float(data.loc[data.index[i], "Y"+element_array[j+1]])) #get the yields of all components #print(sim_E) for E in range(len(sim_E)): sweepdata.append([]) #list for every energy for angle in range(len(inc_angles)): sweepdata[E].append([]) #list for every angle sweepdata[E][angle].append([]) # list for every angle [E][Angle][0] yields, [E][Angle][1]=refl. sweepdata[E][angle].append([]) for E in range(len(sim_E)): for angle in range(len(inc_angles)): for i in range(0, len(data), 1): # go through each line of the data file if data.loc[data.index[i], "Energy"] == sim_E[E]: # interested in the lines where E=certain value if data.loc[data.index[i], "Angle"] == inc_angles[angle]: # interested in the lines where E=certain value for j in range(len(element_array) - 1): # inc_angles[j].append(float(data.loc[data.index[i], "Angle"])) sweepdata[E][angle][0].append(float(data.loc[data.index[i], "Y" + element_array[j + 1]])) # get the yields of all components sweepdata[E][angle][1].append(float(data.loc[data.index[i], "R" + element_array[0]])) # get the refl coeff of all components spacing_E = sim_E[1]-sim_E[0] spacing_Angle = inc_angles[1]-inc_angles[0] #print(sweepdata) #print(sweepdata[0]) #print(sweepdata[0][0]) #print(sweepdata[0][0][0]) for j in range(len(element_array)-1): #Loop over target components; relevant for yields, sputtered atoms data for i in angles: #file_atom = path +"/Data/Atoms/partic_back_r_" + str(i) + ".dat" file_atom = path + "/Atoms/partic_back_r_" + str(i) + ".dat" intermediate_data = pd.read_csv(file_atom, skiprows=4, sep="\s+") #transform SDTrimSP coordinate system into Tri3Dyn coordinate system #append trajectories with mirrored trajectories #(SDTrimSP gives information in the form of a cosine, which is not bijective in the alpha range [0:2pi]) #In the SDTrimSP files, trajectories are included in the form of the cosine of a polar and acimuthal angle, #which are named cosp and cosa, respectively. #converting to cartesian coordinates, this gives: #x = -cosp #y = sinp * cosa #z= sinp * sina #to transform from SDTrimSP to Tri3Dyn, a rotation about the -x axis is necessary such that #x' = -cosp #y' = isnp * sina #z' = -sinp * cosa # coordinate_matrix = np.array((-1.*intermediate_data.loc[intermediate_data['species'] == j+2]["cosp"], \ # np.sin(np.arccos(intermediate_data.loc[intermediate_data['species'] == j+2]["cosp"]))*\ # np.sin(np.arccos(intermediate_data.loc[intermediate_data['species'] == j+2]["cosa"])),\ # -1.0 * np.sin(np.arccos(intermediate_data.loc[intermediate_data['species'] == j+2]["cosp"]))\ # *intermediate_data.loc[intermediate_data['species'] == j+2]["cosa"])).transpose() coordinate_matrix = np.array( (-1. * intermediate_data.loc[intermediate_data['species'] == j + 2]["cosp"], \ np.sin(np.arccos(intermediate_data.loc[intermediate_data['species'] == j + 2]["cosp"])) * \ np.sin(np.arccos(intermediate_data.loc[intermediate_data['species'] == j + 2]["cosa"])), \ -1.0 * np.sin(np.arccos(intermediate_data.loc[intermediate_data['species'] == j + 2]["cosp"])) \ * intermediate_data.loc[intermediate_data['species'] == j + 2]["cosa"], intermediate_data.loc[intermediate_data['species'] == j + 2]["end-energy"])).transpose() mirrored_coordinates = np.array((coordinate_matrix[:,0], -1.*coordinate_matrix[:,1], coordinate_matrix[:,2], coordinate_matrix[:,3])).transpose() # add energy parameter to coordinate matrix, at last index, to include energy to data_Atoms #write trajectories in the array location corresponding to the target component, write length of trajectory arrays data_atoms[j].update({i: np.append(coordinate_matrix, mirrored_coordinates, axis=0)}) length_atoms[j].update({i: 2*len(coordinate_matrix)}) #data_atoms[j].update({i: np.genfromtxt(file_atom)}) for i in angles: # For ions, only the beam is relevant, no looping over composition is necessary # file_ion = path + "/Data/Ions/partic_back_p_" + str(i) +".dat" file_ion = path + "/Ions/partic_back_p_" + str(i) + ".dat" #same coordinate transformation as before intermediate_data = pd.read_csv(file_ion, skiprows=4, sep="\s+") coordinate_matrix = np.array((intermediate_data["end-energy"], -1.*intermediate_data["cosp"], \ np.sin(np.arccos(intermediate_data["cosp"]))*np.sin(np.arccos(intermediate_data["cosa"])), \ -1.0 * np.sin(np.arccos(intermediate_data["cosp"]))*intermediate_data["cosa"])).transpose() mirrored_coordinates = np.array((coordinate_matrix[:,0], coordinate_matrix[:,1],-1.*coordinate_matrix[:,2], coordinate_matrix[:,3])).transpose() data_ions.update({i: np.append(coordinate_matrix, mirrored_coordinates, axis=0)}) length_ions.update({i: 2*len(coordinate_matrix)}) else: #Tri3dyn input files (simulated via beta = 90deg; ions impinging from positive z Axis in TRI3DYN coordinate system) print("Tri3dyn repository data found. Loading...") repo_info.append("TRI3DYN") # IMPORTANT: Declaration array needs to be defined here: (because TRI3DYN has no header!) #element_array = ["Ar", "W"] # here for pure target of W, first is the projectile ''' teststring = linecache.getline(list_of_files[0], 2) beginn = teststring.find("[") ende = teststring.find("]") # split string of composition into a python list element_array = teststring[beginn + 1:ende].split(", ") # get rid of quotation marks, such that only the element symbols remain element_array[:] = [x[2:-2] for x in element_array] ''' # write elemental masses for element in element_array: masses.append(elemental_information.loc[element, 'mass']) # Initialisation of Arrays of dictionaries/arrays data_atoms = [] length_atoms = [] # append as many dictionaries as necessary for i in range(len(element_array) - 1): data_atoms.append({}) length_atoms.append({}) yields.append([]) # inc_angles.append([]) for i in range(0, len(data), 1): #go through each line of the data file if data.loc[data.index[i], "Angle"] == 0.0: #interested in the lines where angle = 0.0 sim_E.append(float(data.loc[data.index[i], "Energy"])) # sim_E should contain all simulated Energies # get the yields for a flat surface for all the components of the target for i in range(0, len(data), 1): # go through each line of the data file if data.loc[data.index[i], "Energy"] == max(sim_E): # interested in the lines where E=2keV inc_angles.append(float(data.loc[data.index[i], "Angle"])) for j in range(len(element_array) - 1): # inc_angles[j].append(float(data.loc[data.index[i], "Angle"])) yields[j].append( float(data.loc[data.index[i], "Y" + element_array[j + 1]])) # get the yields of all components for E in range(len(sim_E)): sweepdata.append([]) #list for every energy for angle in range(len(inc_angles)): sweepdata[E].append([]) #list for every angle sweepdata[E][angle].append([]) # list for every angle [E][Angle][0] yields, [E][Angle][1]=refl. sweepdata[E][angle].append([]) #print(sim_E) for E in range(len(sim_E)): for angle in range(len(inc_angles)): for i in range(0, len(data), 1): # go through each line of the data file if data.loc[data.index[i], "Energy"] == sim_E[E]: # interested in the lines where E=certain value if data.loc[data.index[i], "Angle"] == inc_angles[angle]: # interested in the lines where E=certain value for j in range(len(element_array) - 1): # inc_angles[j].append(float(data.loc[data.index[i], "Angle"])) sweepdata[E][angle][0].append(float(data.loc[data.index[i], "Y" + element_array[j + 1]])) # get the yields of all components sweepdata[E][angle][1].append(float(data.loc[data.index[i], "R" + element_array[0]])) # get the refl coeff of all components spacing_E = sim_E[1]-sim_E[0] spacing_Angle = inc_angles[1]-inc_angles[0] for i in range(len(element_array) - 1): for j in range(0, 90, 5): # fill dict with total results files (with all species mixed) # file_atom = path + "/Data/Atoms/Ar_Wflt_2k_red_" + str(j) + "deg/Ar_Wflt_2k_splst_red_" + str(j) + "deg_elem_" + str(float(i+2.0)) + ".dat" file_atom = path + "/Atoms/Ar_Wflt_2k_red_" + str(j) + "deg/Ar_Wflt_2k_splst_red_" + str(j) + "deg_elem_" + str(float(i + 2.0)) + ".dat" data_atoms[i].update({j: np.genfromtxt(file_atom)}) length_atoms[i].update({j: int(len(data_atoms[i][j]))}) for j in range(0, 90, 5): # For ions, only the beam is relevant, no looping over composition is necessary. For TRI3dYN always this 5deg from 0 to 85deg interval is needed!!! #file_ion = path + "/Data/Ions/Ar_Wflt_2k_bslst_red_" + str(j) + "deg.dat" file_ion = path + "/Ions/Ar_Wflt_2k_bslst_red_" + str(j) + "deg.dat" data_ions.update({j: np.genfromtxt(file_ion)}) length_ions.update({j: int(len(data_ions[j]))}) return data_ions, length_ions, data_atoms, length_atoms, inc_angles, yields, element_array, masses, sweepdata, spacing_Angle, sim_E, spacing_E, main_input, sim_E, repo_info