Commit 68f507c9 authored by Parish, Chad's avatar Parish, Chad
Browse files

Add Agarwal et al. citation.

parent daf3e414
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# Rule #1: generate figures programmatically # Rule #1: generate figures programmatically
<br> <br>
This parses the annoyingly formated .txt files output by SRIM (http://www.srim.org/), uses the parsed data and minimal user input, and generates well-formatted figures of the damage in dpa (displacements per atom) and the ion implantation level in the material identified by parsing the input files. This parses the annoyingly formated .txt files output by SRIM (http://www.srim.org/), uses the parsed data and minimal user input, and generates well-formatted figures of the damage in dpa (displacements per atom) and the ion implantation level in the material identified by parsing the input files.
<br> <br>
<br> <br>
The principle illustrated by this notebook is that a small investment in time can result in an easily reusable code that will make the subsequent generation of figures fast and efficient. The principle illustrated by this notebook is that a small investment in time can result in an easily reusable code that will make the subsequent generation of figures fast and efficient.
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Imports: ### Imports:
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
import tkinter import tkinter
from tkinter import filedialog from tkinter import filedialog
import os import os
import glob import glob
import numpy as np import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from matplotlib import gridspec from matplotlib import gridspec
from warnings import warn from warnings import warn
plt.style.use({'font.sans-serif': 'arial', 'font.size': 16, 'font.weight': 'bold', 'svg.fonttype': 'none'}) plt.style.use({'font.sans-serif': 'arial', 'font.size': 16, 'font.weight': 'bold', 'svg.fonttype': 'none'})
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Find and read in the needed text files ### Find and read in the needed text files
In this case, PHONON.txt, RANGE.txt, and TDATA.txt are needed In this case, PHONON.txt, RANGE.txt, and TDATA.txt are needed
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# get the data directory # get the data directory
root = tkinter.Tk() root = tkinter.Tk()
path_dir = filedialog.askdirectory( path_dir = filedialog.askdirectory(
initialdir="/", initialdir="/",
title="Find SRIM data:") title="Find SRIM data:")
print(f'Path: {path_dir}') print(f'Path: {path_dir}')
root.destroy() # remove the GUI panel root.destroy() # remove the GUI panel
os.chdir(path_dir) os.chdir(path_dir)
# find PHONON.txt # find PHONON.txt
fn = glob.glob('PHONON.txt') fn = glob.glob('PHONON.txt')
assert len(fn) == 1, 'Did not find exactly 1 PHONON.txt!' assert len(fn) == 1, 'Did not find exactly 1 PHONON.txt!'
file_lines_phonon = [] file_lines_phonon = []
with open(fn[0], mode='r') as fh: with open(fn[0], mode='r') as fh:
for line in fh: for line in fh:
file_lines_phonon.append(line) file_lines_phonon.append(line)
else: else:
pass pass
print(f'Lines in PHONON.txt: {len(file_lines_phonon)}') print(f'Lines in PHONON.txt: {len(file_lines_phonon)}')
# find TDATA.txt # find TDATA.txt
fn = glob.glob('TDATA.txt') fn = glob.glob('TDATA.txt')
assert len(fn) == 1, 'Did not find exactly 1 TDATA.txt!' assert len(fn) == 1, 'Did not find exactly 1 TDATA.txt!'
file_lines_tdata = [] file_lines_tdata = []
with open(fn[0], mode='r') as fh: with open(fn[0], mode='r') as fh:
for line in fh: for line in fh:
file_lines_tdata.append(line) file_lines_tdata.append(line)
else: else:
pass pass
print(f'Lines in TDATA.txt: {len(file_lines_tdata)}') print(f'Lines in TDATA.txt: {len(file_lines_tdata)}')
# find RANGE.txt # find RANGE.txt
fn = glob.glob('RANGE.txt') fn = glob.glob('RANGE.txt')
assert len(fn) == 1, 'Did not find exactly 1 RANGE.txt!' assert len(fn) == 1, 'Did not find exactly 1 RANGE.txt!'
file_lines_range = [] file_lines_range = []
with open(fn[0], mode='r') as fh: with open(fn[0], mode='r') as fh:
for line in fh: for line in fh:
file_lines_range.append(line) file_lines_range.append(line)
else: else:
pass pass
print(f'Lines in RANGE.txt: {len(file_lines_range)}') print(f'Lines in RANGE.txt: {len(file_lines_range)}')
``` ```
%% Output %% Output
Path: C:/work/Collaborations/isotope_production/30MeV_He_on_Bi Path: C:/work/Collaborations/isotope_production/30MeV_He_on_Bi
Lines in PHONON.txt: 124 Lines in PHONON.txt: 124
Lines in TDATA.txt: 39 Lines in TDATA.txt: 39
Lines in RANGE.txt: 137 Lines in RANGE.txt: 137
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Parse out the important information from PHONON.txt ### Parse out the important information from PHONON.txt
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
#Find the important lines: #Find the important lines:
density_found = 0 density_found = 0
for i, l in enumerate(file_lines_phonon): for i, l in enumerate(file_lines_phonon):
if 'DEPTH' in l: # line that defines the data columns if 'DEPTH' in l: # line that defines the data columns
i_start = i + 3 # find the first line of data i_start = i + 3 # find the first line of data
if 'atoms/cm3' in l: # find the density if 'atoms/cm3' in l: # find the density
i_density = i i_density = i
density_found += 1 density_found += 1
assert density_found <=1, 'Can currently only handle single layers' assert density_found <=1, 'Can currently only handle single layers'
# Calculate the length of the data. SRIM is always 100 # Calculate the length of the data. SRIM is always 100
length = len(file_lines_phonon) - i_start length = len(file_lines_phonon) - i_start
data_phonon = np.zeros((length, 3), dtype='f4') # hold the data here data_phonon = np.zeros((length, 3), dtype='f4') # hold the data here
# parse out the density (atoms/cm^3) # parse out the density (atoms/cm^3)
density_line = file_lines_phonon[i_density] density_line = file_lines_phonon[i_density]
start_phrase = 'Density =' start_phrase = 'Density ='
start_phrase_length = len(start_phrase) start_phrase_length = len(start_phrase)
end_phrase = 'atoms/cm3' end_phrase = 'atoms/cm3'
st = density_line.find(start_phrase) + start_phrase_length st = density_line.find(start_phrase) + start_phrase_length
en = density_line.find(end_phrase) en = density_line.find(end_phrase)
density = np.array(density_line[st:en]).astype('f4') density = np.array(density_line[st:en]).astype('f4')
print(f'Density: {density:3.3} atoms/cm^3') print(f'Density: {density:3.3} atoms/cm^3')
# read the data # read the data
for i in range(i_start, len(file_lines_phonon),1): for i in range(i_start, len(file_lines_phonon),1):
data_phonon[i-i_start,:] = np.array(file_lines_phonon[i].split()).astype('f4') data_phonon[i-i_start,:] = np.array(file_lines_phonon[i].split()).astype('f4')
``` ```
%% Output %% Output
Density: 2.82e+22 atoms/cm^3 Density: 2.82e+22 atoms/cm^3
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Parse out the important information from TDATA.txt ### Parse out the important information from TDATA.txt
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
kp_flag = False # Check for Kinchin-Pease kp_flag = False # Check for Kinchin-Pease
# loop over the TDATA.txt lines # loop over the TDATA.txt lines
for i,l in enumerate(file_lines_tdata): for i,l in enumerate(file_lines_tdata):
# find the projectile # find the projectile
if 'Ion = ' in l: if 'Ion = ' in l:
st = l.find('Ion = ') st = l.find('Ion = ')
en = l.find('(') en = l.find('(')
ion_string = l[st+6:en-1] ion_string = l[st+6:en-1]
print(f'Ion string: "{ion_string}"') print(f'Ion string: "{ion_string}"')
#find the displacement energy #find the displacement energy
if 'Displacement = ' in l: if 'Displacement = ' in l:
st = l.find('Displacement = ') st = l.find('Displacement = ')
en = l.find(',') en = l.find(',')
displacement = float(l[st+14:en-2]) displacement = float(l[st+14:en-2])
print(f'Displacement: {displacement} eV') print(f'Displacement: {displacement} eV')
# Check for Kincihn-Pease model # Check for Kincihn-Pease model
if 'Kinchin-Pease Estimates' in l: if 'Kinchin-Pease Estimates' in l:
kp_flag = True kp_flag = True
# Find the ion (projectile) energy # Find the ion (projectile) energy
if 'Energy = ' in l: if 'Energy = ' in l:
st = l.find('Energy = ') st = l.find('Energy = ')
en = l.find('keV') # if keV not found, find() returns -1: en = l.find('keV') # if keV not found, find() returns -1:
if en > -1: if en > -1:
energy = float(l[st+10:en-1]) energy = float(l[st+10:en-1])
else: else:
energy = float(input('Energy not parsable -- input ion energy in keV:')) energy = float(input('Energy not parsable -- input ion energy in keV:'))
print(f'Energy: {energy} keV') print(f'Energy: {energy} keV')
# Find the identity of layer 1 # Find the identity of layer 1
if 'TARGET MATERIAL' in l: if 'TARGET MATERIAL' in l:
temp_l = file_lines_tdata[i+1] temp_l = file_lines_tdata[i+1]
st = temp_l.find('Layer # 1 - ') st = temp_l.find('Layer # 1 - ')
en = temp_l.find('\n') en = temp_l.find('\n')
target = temp_l[st+11:en] target = temp_l[st+11:en]
print(f'Target: "{target}"') print(f'Target: "{target}"')
# Check for multilayer targets # Check for multilayer targets
if 'Layer #2' in l: if 'Layer #2' in l:
warn('Currently this code only parses single layers properly') warn('Currently this code only parses single layers properly')
# find # ions calculated # find # ions calculated
if 'Total Ions calculated =' in l: if 'Total Ions calculated =' in l:
st = l.find('Total Ions calculated =') st = l.find('Total Ions calculated =')
en = l.find('\n') en = l.find('\n')
num_ions = int(l[st+23:en]) num_ions = int(l[st+23:en])
print(f'Number of ions: {num_ions}') print(f'Number of ions: {num_ions}')
# If Kinchin-Pease not found... # If Kinchin-Pease not found...
if kp_flag is False: if kp_flag is False:
warn('Warning: this code is intended to use Kinchin-Pease data.') warn('Warning: this code is intended to use Kinchin-Pease data.')
``` ```
%% Output %% Output
Ion string: "He" Ion string: "He"
Energy: 30000.0 keV Energy: 30000.0 keV
Target: " Layer 1" Target: " Layer 1"
Displacement: 23.0 eV Displacement: 23.0 eV
Number of ions: 200000 Number of ions: 200000
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Parse out the important information from RANGE.txt ### Parse out the important information from RANGE.txt
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
#Find the important lines: #Find the important lines:
for i, l in enumerate(file_lines_range): for i, l in enumerate(file_lines_range):
if 'DEPTH' in l: # line that defines the data columns if 'DEPTH' in l: # line that defines the data columns
i_start = i + 3 # find the first line of data i_start = i + 3 # find the first line of data
# Calculate the length of the data. SRIM is always 100 # Calculate the length of the data. SRIM is always 100
length = len(file_lines_range) - i_start length = len(file_lines_range) - i_start
data_range = np.zeros((length, 3), dtype='f4') # hold the data here data_range = np.zeros((length, 3), dtype='f4') # hold the data here
# read the data # read the data
for i in range(i_start, len(file_lines_range),1): for i in range(i_start, len(file_lines_range),1):
data_range[i-i_start,:] = np.array(file_lines_range[i].split()).astype('f4') data_range[i-i_start,:] = np.array(file_lines_range[i].split()).astype('f4')
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Perform the calculation ### Perform the calculation
The details here are from:<br> The details here are from:<br>
Stoller, R.E., Toloczko, M.B., Was, G.S., Certain, A.G., Dwaraknath, S. and Garner, F.A., 2013. On the use of SRIM for computing radiation damage exposure. <i>Nuclear instruments and methods in physics research section B: beam interactions with materials and atoms</i>, <b>310</b>, pp.75-80. Stoller, R.E., Toloczko, M.B., Was, G.S., Certain, A.G., Dwaraknath, S. and Garner, F.A., 2013. On the use of SRIM for computing radiation damage exposure. <i>Nuclear instruments and methods in physics research section B: beam interactions with materials and atoms</i>, <b>310</b>, pp.75-80.
<span style="color:firebrick">Addendum 2021-07-23:</span>
Please see the recent paper S. Agarwal, Y. Lin, C. Li, R. E. Stoller, and S. J. Zinkle, "On the use of SRIM for calculating vacancy production: Quick calculation and full-cascade options," <i>Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms</i>, <b>503</b>, pp.11-29, prior to using results output from this code in order to be sure you understand the proper range of application and internal calculation options of SRIM.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
fluence = float(input('Please enter the fluence (ions/cm^2): ')) fluence = float(input('Please enter the fluence (ions/cm^2): '))
# Calculate the NRT damage in dpa # Calculate the NRT damage in dpa
depth_in_nm = data_phonon[:, 0] / 10 # always angstroms depth_in_nm = data_phonon[:, 0] / 10 # always angstroms
ion_contribution = data_phonon[:, 1] ion_contribution = data_phonon[:, 1]
recoil_contribution = data_phonon[:, 2] recoil_contribution = data_phonon[:, 2]
P = ion_contribution + recoil_contribution # (P)honons P = ion_contribution + recoil_contribution # (P)honons
T_damage = P * fluence / density * 1e8 # 1e8 is angstrom --> cm T_damage = P * fluence / density * 1e8 # 1e8 is angstrom --> cm
NRT = T_damage * 0.8 / 2.0 / displacement # NRT displacements per atom. See citation. NRT = T_damage * 0.8 / 2.0 / displacement # NRT displacements per atom. See citation.
# Calculate the ion implantation in atomic fraction # Calculate the ion implantation in atomic fraction
atom_frac = data_range[:, 1] * fluence / density atom_frac = data_range[:, 1] * fluence / density
``` ```
%% Output %% Output
Please enter the fluence (ions/cm^2): 1.3e18 Please enter the fluence (ions/cm^2): 1.3e18
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
with plt.style.context('seaborn'): with plt.style.context('seaborn'):
fig = plt.figure(figsize=(12, 7)) fig = plt.figure(figsize=(12, 7))
gs = gridspec.GridSpec(2, 1, height_ratios=[4, 1]) gs = gridspec.GridSpec(2, 1, height_ratios=[4, 1])
ax0 = plt.subplot(gs[0]) ax0 = plt.subplot(gs[0])
ax0.plot(depth_in_nm, NRT, color='xkcd:brick red') ax0.plot(depth_in_nm, NRT, color='xkcd:brick red')
ax1 = plt.subplot(gs[1]) ax1 = plt.subplot(gs[1])
ax1.plot(depth_in_nm, atom_frac, color='xkcd:brick red') ax1.plot(depth_in_nm, atom_frac, color='xkcd:brick red')
``` ```
%% Output %% Output
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
try: try:
err_fluence = float(input('Uncertainty in fluence (for "+/-20%" input "0.2"): ')) err_fluence = float(input('Uncertainty in fluence (for "+/-20%" input "0.2"): '))
except ValueError: except ValueError:
err_fluence = 0 err_fluence = 0
try: try:
err_displace = float(input('Uncertainty in displacement energy (for "+/-20%" input "0.2"): ')) err_displace = float(input('Uncertainty in displacement energy (for "+/-20%" input "0.2"): '))
except ValueError: except ValueError:
err_displace = 0 err_displace = 0
err_both = np.sqrt( err_fluence * err_fluence + err_displace * err_displace ) err_both = np.sqrt( err_fluence * err_fluence + err_displace * err_displace )
print(f'Combined error: {err_both}') print(f'Combined error: {err_both}')
``` ```
%% Output %% Output
Uncertainty in fluence (for "+/-20%" input "0.2"): 0 Uncertainty in fluence (for "+/-20%" input "0.2"): 0
Uncertainty in displacement energy (for "+/-20%" input "0.2"): 0 Uncertainty in displacement energy (for "+/-20%" input "0.2"): 0
Combined error: 0.0 Combined error: 0.0
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Now, generate the figure ### Now, generate the figure
This does require some hand-tweaking of parameters on the first attempt, but all subsequent uses will take only a few seconds, where most of the time is spent waiting for the user to click on the new data directory, or to type in the fluence or error. This does require some hand-tweaking of parameters on the first attempt, but all subsequent uses will take only a few seconds, where most of the time is spent waiting for the user to click on the new data directory, or to type in the fluence or error.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# the 'with' block sets temporary parameters: # the 'with' block sets temporary parameters:
with plt.style.context({'font.sans-serif': 'Arial', 'svg.fonttype': 'none'}): with plt.style.context({'font.sans-serif': 'Arial', 'svg.fonttype': 'none'}):
#Create a figure with unequal-sized subplots #Create a figure with unequal-sized subplots
fig = plt.figure(figsize=(12, 7)) fig = plt.figure(figsize=(12, 7))
gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1]) gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1])
# create the first subplot # create the first subplot
ax0 = plt.subplot(gs[0]) ax0 = plt.subplot(gs[0])
ax0.plot(depth_in_nm, NRT, color='xkcd:brick red', label='Damage') # plot the NRT damage ax0.plot(depth_in_nm, NRT, color='xkcd:brick red', label='Damage') # plot the NRT damage
if err_both > 0: # if the user input error bounds, plot the error bounds if err_both > 0: # if the user input error bounds, plot the error bounds
ax0.fill_between(depth_in_nm, y1=(NRT * (1-err_both)), y2=(NRT * (1+err_both)), ax0.fill_between(depth_in_nm, y1=(NRT * (1-err_both)), y2=(NRT * (1+err_both)),
color='xkcd:merlot', alpha=0.15, label='Error bounds') color='xkcd:merlot', alpha=0.15, label='Error bounds')
ax0.legend(loc=3, prop={'size': 14}) # only show legend if error bounds present ax0.legend(loc=3, prop={'size': 14}) # only show legend if error bounds present
ax0.set_ylabel('NRT damage, dpa', fontsize=20) ax0.set_ylabel('NRT damage, dpa', fontsize=20)
# Make annotations # Make annotations
pre, post = str(fluence * 1e4).split('e+') # split fluence into A x 10^B, with conversion from cm^2 to m^2 pre, post = str(fluence * 1e4).split('e+') # split fluence into A x 10^B, with conversion from cm^2 to m^2
ion_val_string = str(num_ions) + ' ' + ion_string + ' ions $\Rightarrow$' + target ion_val_string = str(num_ions) + ' ' + ion_string + ' ions $\Rightarrow$' + target
fluence_string = '$' + pre + r'\times 10^{' + post + '}$ ions/m$^2$' fluence_string = '$' + pre + r'\times 10^{' + post + '}$ ions/m$^2$'
disp_string = str(displacement) +' eV $E_D$' disp_string = str(displacement) +' eV $E_D$'
if err_both >0: if err_both >0:
flu_err_string = r'Fluence $\pm$ ' + str(err_fluence*100) + '%' flu_err_string = r'Fluence $\pm$ ' + str(err_fluence*100) + '%'
disp_err_string = r'$E_D \pm$ ' + str(err_displace*100) + '%' disp_err_string = r'$E_D \pm$ ' + str(err_displace*100) + '%'
else: else:
flu_err_string = '' flu_err_string = ''
disp_err_string = '' disp_err_string = ''
ax0.text(0.98, 0.98, f'{ion_val_string}\n{fluence_string}\n{disp_string}\n{flu_err_string}\n{disp_err_string}', ax0.text(0.98, 0.98, f'{ion_val_string}\n{fluence_string}\n{disp_string}\n{flu_err_string}\n{disp_err_string}',
transform=ax0.transAxes, fontsize=14, verticalalignment='top', transform=ax0.transAxes, fontsize=14, verticalalignment='top',
horizontalalignment='right') horizontalalignment='right')
ax0.tick_params(axis='both', labelsize=20) ax0.tick_params(axis='both', labelsize=20)
ax0.spines['right'].set_visible(False) ax0.spines['right'].set_visible(False)
ax0.spines['top'].set_visible(False) ax0.spines['top'].set_visible(False)
# create the smaller plot # create the smaller plot
ax1 = plt.subplot(gs[1]) ax1 = plt.subplot(gs[1])
if atom_frac.max() < 0.01: # Plot atom fraction if <1% if atom_frac.max() < 0.01: # Plot atom fraction if <1%
atom_string = 'Implanted\natom fraction' atom_string = 'Implanted\natom fraction'
atom_mult = 1 atom_mult = 1
else: else:
atom_string = 'Implanted\natom %' atom_string = 'Implanted\natom %'
atom_mult = 100 atom_mult = 100
ax1.plot(depth_in_nm, atom_frac*atom_mult, color='xkcd:brick red', label='Implantation') ax1.plot(depth_in_nm, atom_frac*atom_mult, color='xkcd:brick red', label='Implantation')
if err_both > 0: # if user input error bounds if err_both > 0: # if user input error bounds
ax1.fill_between(depth_in_nm, y1=(atom_frac * (1-err_fluence))*atom_mult, ax1.fill_between(depth_in_nm, y1=(atom_frac * (1-err_fluence))*atom_mult,
y2=(atom_frac * (1+err_fluence))*atom_mult, y2=(atom_frac * (1+err_fluence))*atom_mult,
color='xkcd:merlot', alpha=0.15, label='Error bounds') color='xkcd:merlot', alpha=0.15, label='Error bounds')
ax1.legend(prop={'size': 14}) # only show legend if error bounds present ax1.legend(prop={'size': 14}) # only show legend if error bounds present
ax1.set_ylabel(atom_string, fontsize=20) ax1.set_ylabel(atom_string, fontsize=20)
ax1.set_xlabel('Depth, nm', fontsize=20) ax1.set_xlabel('Depth, nm', fontsize=20)
ax1.tick_params(axis='both', labelsize=20) ax1.tick_params(axis='both', labelsize=20)
ax1.spines['right'].set_visible(False) ax1.spines['right'].set_visible(False)
ax1.spines['top'].set_visible(False) ax1.spines['top'].set_visible(False)
fn = input('File name (enter to skip):') fn = input('File name (enter to skip):')
if len(fn) > 0: if len(fn) > 0:
fn = '..\\' + fn # go up one directory fn = '..\\' + fn # go up one directory
fig.savefig(fn + '.png', dpi=300) fig.savefig(fn + '.png', dpi=300)
fig.savefig(fn + '.svg', dpi=300) fig.savefig(fn + '.svg', dpi=300)
fig.savefig(fn + '.eps', dpi=300) fig.savefig(fn + '.eps', dpi=300)
``` ```
%% Output %% Output
File name (enter to skip):5hours File name (enter to skip):5hours
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
NRT[50] NRT[50]
``` ```
%% Output %% Output
0.45283368 0.45283368
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
``` ```
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment