import json import numpy as np import argparse import os def energy_to_wavelength(energy): ''' Convert electronvolt to nanometer ''' with np.errstate(divide='ignore'): return 1239.842 / energy def wavelength_to_energy(wavelength): ''' Convert nanometer to electronvolt ''' return 1239.842 / wavelength def get_fitted_values(target_concentrations, concentrations, real_sets, imag_sets, order=6): ''' Fit dielectric functions and evaluate at target_concentrations. ''' fitted_real_sets = [] fitted_imag_sets = [] N = real_sets.shape[1] # Fit real parts for i in range(N): polynomial = np.polyfit(concentrations, real_sets[:, i], order) eps_real = np.poly1d(polynomial)(target_concentrations) fitted_real_sets.append(eps_real) # Fit imag parts for i in range(N): polynomial = np.polyfit(concentrations, imag_sets[:, i], order) eps_imag = np.poly1d(polynomial)(target_concentrations) fitted_imag_sets.append(eps_imag) fitted_real_sets = np.array(fitted_real_sets).transpose() fitted_imag_sets = np.array(fitted_imag_sets).transpose() return fitted_real_sets, fitted_imag_sets def get_dielectric_data(compound, database='dielectric-functions.json'): ''' Get dielectric data, along with energies, from json file ''' compound = '-'.join(compound) concentrations = [] energy_sets = [] real_sets = [] imag_sets = [] with open(database, 'r') as f: dfs = json.load(f) if compound not in dfs: raise ValueError('There are no records for {} in {}'.format(compound, database)) for df in dfs[compound].values(): concentrations.append(df['concentration']) energy_sets.append(df['energies']) real_sets.append(df['real']) imag_sets.append(df['imag']) # If elements were fed in in reverse alphabetical order, # we flip the concentration # Sort according to concentration and return NumPy arrays concentrations, energy_sets, real_sets, imag_sets = \ zip(*sorted(zip(concentrations, energy_sets, real_sets, imag_sets))) return concentrations, energy_sets, real_sets, imag_sets def homogenize_energy_sets(energy_sets, real_sets, imag_sets, tol=1e-6): ''' If the dielectric data correspond to different energies, we need homogenize the data. ''' real_sets_cut = [] imag_sets_cut = [] # Find the largest energy interval that all sets contain maximum_min_e = -1 minimum_max_e = 1e6 for energy_set in energy_sets: maximum_min_e = max(maximum_min_e, min(energy_set)) minimum_max_e = min(minimum_max_e, max(energy_set)) # Pick out the interval that each data set contains. Will only succeed if # the energies match (otherwise one would need to interpolate, which is # not yet implemented) energies = None for energy_set, real_set, imag_set in zip(energy_sets, real_sets, imag_sets): min_index = np.nonzero( abs(np.array(energy_set) - maximum_min_e) < tol)[0][0] max_index = np.nonzero( abs(np.array(energy_set) - minimum_max_e) < tol)[0][0] energy_new = energy_set[min_index:max_index + 1] if energies is None: energies = energy_new else: assert len(energy_new) == len(energies) assert np.allclose(energy_new, energies) real_sets_cut.append(real_set[min_index:max_index + 1]) imag_sets_cut.append(imag_set[min_index:max_index + 1]) return (np.array(energies), np.array(real_sets_cut), np.array(imag_sets_cut)) def fit_dielectric_function(element1, element2, target_concentrations, order=6, database='dielectric-functions.json'): ''' Fit dielectric function at unknown concentrations to the already calculated. For each energy that all dielectric functions for the given alloy has been calculated, a polynomial fit (default orer: 6) is performed). The dielectric function is then evaluated at the desired concentration. Paramaters ---------- element1 : str chemical species element2 : str chemical specise target_concentrations : list of floats concentrations of element2 at which the dielectric functions should be fitted. order : int order of polynomial to which the dielectric function will be fitted. database : stri filename for the database containing the dielectric functions in json format Returns ------- list of list of complex numbers dielectric function at each of the target concentrations list of floats energies at which dielectric funnctions are evaluated ''' compound = sorted([element1, element2]) if min(target_concentrations) < 0 or max(target_concentrations) > 100: raise ValueError('target_concentrations must all fall be within ' 'the range [0, 100]') if element1 == compound[1]: # If order was swithched target_concentrations = 100 - np.array(target_concentrations) else: target_concentrations = np.array(target_concentrations) target_concentrations = target_concentrations / 100 concentrations, energy_sets, real_sets, imag_sets \ = get_dielectric_data(compound, database=database) energies, real_sets, imag_sets = homogenize_energy_sets( energy_sets, real_sets, imag_sets) real_sets, imag_sets = get_fitted_values(target_concentrations, concentrations, real_sets, imag_sets, order=order) return energies, real_sets, imag_sets if __name__ == '__main__': __doc__ = '''Interpolate dielectric functions.''' parser = argparse.ArgumentParser(description=__doc__) parser.add_argument('element1', help='Element number 1') parser.add_argument('element2', help='Element number 2') parser.add_argument('target_concentrations', nargs='+', action='append', type=float, help='Target concentrations (as a percentage) ' 'of element 2') defval = 'df-' parser.add_argument('-f', '--filename', default=defval, help='Preamble to filename. The filename will always ' 'end with something like -Au63-Pd37.txt.') parser.add_argument('-p', '--plot', action='store_true', help='Skip file-writing and plot instead.') parser.add_argument('-w', '--wavelength', action='store_true', help='Convert to wavelength') parser.add_argument('--xmin', type=float, help='Minimum energy (or wavelength if -w is used).' 'If specified and lower than available data, the ' 'minimum energy/wavelength available will be used. ' 'Default is to use all available data, but note that ' 'it is not expected to be accurate at high energies.') parser.add_argument('--xmax', type=float, help='Maximum energy (or wavelength if -w is used. ' 'If specified and higher than available data, the ' 'maximum energy/wavelength available will be used. ' 'Default is to use all available data, but note that ' 'it is not expected to be accurate at high energies.') defval = 6 parser.add_argument('--order', type=int, default=defval, help='Polynomial order of the fit. ' 'Default is {}.'.format(defval)) defval = 'dielectric_functions.json' parser.add_argument('-d', '--database', default=defval, help='Database in json format containing dielectric ' 'data') args = parser.parse_args() tol = 1e-5 elements = ['Au', 'Ag', 'Cu', 'Pd', 'Pt'] if args.element1 not in elements or args.element2 not in elements: print('Error: element1 and element2 need to be {} or {}.'.format( ', '.join(elements[:-1]), elements[-1])) exit(0) target_concentrations = args.target_concentrations[0] energies, real_sets, imag_sets = \ fit_dielectric_function(args.element1, args.element2, target_concentrations, order=args.order, database=args.database) # Define limits if none are specified if args.xmax: if args.wavelength: emin = wavelength_to_energy(args.xmax) else: emax = args.xmax else: if args.wavelength: emin = min(energies) else: emax = max(energies) if args.xmin: if args.wavelength: emax = wavelength_to_energy(args.xmin) else: emin = args.xmin else: if args.wavelength: emax = max(energies) else: emin = min(energies) # If we are not going to save, prepare for plotting if args.plot: import matplotlib.pyplot as plt fig, ax = plt.subplots() # Write to file for i in range(len(target_concentrations)): conc1 = 100 - target_concentrations[i] conc2 = target_concentrations[i] compound = '{}{}-{}{}'.format(args.element1, int(np.round(conc1)), args.element2, int(np.round(conc2))) if not args.plot: # Construct filename using user input filename = args.filename + compound + '.txt' # If name contains directory, create that directory if os.path.dirname(filename): try: os.makedirs(os.path.dirname(filename)) except OSError: pass # Print to file with open(filename, 'w') as f: f.write('# Interpolated from data in {}\n'.format( args.database)) if args.wavelength: f.write('# Col 1: Wavelength (nm)\n') else: f.write('# Col 1: Energy (eV)\n') f.write('# Col 2: Real part (Re(epsilon))\n') f.write('# Col 3: Imag part (Im(epsilon))\n') for energy, real, imag in zip(energies, real_sets[i], imag_sets[i]): # Eliminate points outside boundaries if energy < emin - tol or energy > emax + tol: continue # Convert to wavelength if we want to if args.wavelength: x = energy_to_wavelength(energy) else: x = energy # Write f.write('{:15.6f}'.format(x)) f.write('{:15.6f}'.format(real)) f.write('{:15.6f}\n'.format(imag)) else: # Let's plot if args.wavelength: ax.plot(energy_to_wavelength(energies), imag_sets[i], label=compound) else: ax.plot(energies, imag_sets[i], label=compound) # Finish plot if args.plot: ax.set_ylim([0, 12]) if args.wavelength: xlim_min = energy_to_wavelength(emax) xlim_max = min(energy_to_wavelength(emin), 1e4) ax.set_xlim([xlim_min, xlim_max]) else: ax.set_xlim([emin, emax]) ax.set_ylabel('Imaginary dielectric function') if args.wavelength: ax.set_xlabel('Wavelength (nm)') else: ax.set_xlabel('Energy (eV)') ax.legend() plt.show()