Commit 36a4cdcd authored by Parish, Chad's avatar Parish, Chad
Browse files

'minor'

parent 9af7ccb0
%% Cell type:markdown id: tags:
# Rule #1: generate figures programmatically
<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.
<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.
%% Cell type:markdown id: tags:
### Imports:
%% Cell type:code id: tags:
``` python
import tkinter
from tkinter import filedialog
import os
import glob
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import gridspec
from warnings import warn
plt.style.use({'font.sans-serif': 'arial', 'font.size': 16, 'font.weight': 'bold', 'svg.fonttype': 'none'})
```
%% Cell type:markdown id: tags:
### Find and read in the needed text files
In this case, PHONON.txt, RANGE.txt, and TDATA.txt are needed
%% Cell type:code id: tags:
``` python
# get the data directory
root = tkinter.Tk()
path_dir = filedialog.askdirectory(
initialdir="/",
title="Find SRIM data:")
print(f'Path: {path_dir}')
root.destroy() # remove the GUI panel
os.chdir(path_dir)
# find PHONON.txt
fn = glob.glob('PHONON.txt')
assert len(fn) == 1, 'Did not find exactly 1 PHONON.txt!'
file_lines_phonon = []
with open(fn[0], mode='r') as fh:
for line in fh:
file_lines_phonon.append(line)
else:
pass
print(f'Lines in PHONON.txt: {len(file_lines_phonon)}')
# find TDATA.txt
fn = glob.glob('TDATA.txt')
assert len(fn) == 1, 'Did not find exactly 1 TDATA.txt!'
file_lines_tdata = []
with open(fn[0], mode='r') as fh:
for line in fh:
file_lines_tdata.append(line)
else:
pass
print(f'Lines in TDATA.txt: {len(file_lines_tdata)}')
# find RANGE.txt
fn = glob.glob('RANGE.txt')
assert len(fn) == 1, 'Did not find exactly 1 RANGE.txt!'
file_lines_range = []
with open(fn[0], mode='r') as fh:
for line in fh:
file_lines_range.append(line)
else:
pass
print(f'Lines in RANGE.txt: {len(file_lines_range)}')
```
%% Output
Path: C:/Users/nmp/Dropbox (ORNL)/Visualization_paper/SRIM_outputs
Path: C:/work/Collaborations/isotope_production/30MeV_He_on_Bi
Lines in PHONON.txt: 124
Lines in TDATA.txt: 39
Lines in RANGE.txt: 137
%% Cell type:markdown id: tags:
### Parse out the important information from PHONON.txt
%% Cell type:code id: tags:
``` python
#Find the important lines:
density_found = 0
for i, l in enumerate(file_lines_phonon):
if 'DEPTH' in l: # line that defines the data columns
i_start = i + 3 # find the first line of data
if 'atoms/cm3' in l: # find the density
i_density = i
density_found += 1
assert density_found <=1, 'Can currently only handle single layers'
# Calculate the length of the data. SRIM is always 100
length = len(file_lines_phonon) - i_start
data_phonon = np.zeros((length, 3), dtype='f4') # hold the data here
# parse out the density (atoms/cm^3)
density_line = file_lines_phonon[i_density]
start_phrase = 'Density ='
start_phrase_length = len(start_phrase)
end_phrase = 'atoms/cm3'
st = density_line.find(start_phrase) + start_phrase_length
en = density_line.find(end_phrase)
density = np.array(density_line[st:en]).astype('f4')
print(f'Density: {density:3.3} atoms/cm^3')
# read the data
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')
```
%% Output
Density: 6.34e+22 atoms/cm^3
Density: 2.82e+22 atoms/cm^3
%% Cell type:markdown id: tags:
### Parse out the important information from TDATA.txt
%% Cell type:code id: tags:
``` python
kp_flag = False # Check for Kinchin-Pease
# loop over the TDATA.txt lines
for i,l in enumerate(file_lines_tdata):
# find the projectile
if 'Ion = ' in l:
st = l.find('Ion = ')
en = l.find('(')
ion_string = l[st+6:en-1]
print(f'Ion string: "{ion_string}"')
#find the displacement energy
if 'Displacement = ' in l:
st = l.find('Displacement = ')
en = l.find(',')
displacement = float(l[st+14:en-2])
print(f'Displacement: {displacement} eV')
# Check for Kincihn-Pease model
if 'Kinchin-Pease Estimates' in l:
kp_flag = True
# Find the ion (projectile) energy
if 'Energy = ' in l:
st = l.find('Energy = ')
en = l.find('keV') # if keV not found, find() returns -1:
if en > -1:
energy = float(l[st+10:en-1])
else:
energy = float(input('Energy not parsable -- input ion energy in keV:'))
print(f'Energy: {energy} keV')
# Find the identity of layer 1
if 'TARGET MATERIAL' in l:
temp_l = file_lines_tdata[i+1]
st = temp_l.find('Layer # 1 - ')
en = temp_l.find('\n')
target = temp_l[st+11:en]
print(f'Target: "{target}"')
# Check for multilayer targets
if 'Layer #2' in l:
warn('Currently this code only parses single layers properly')
# find # ions calculated
if 'Total Ions calculated =' in l:
st = l.find('Total Ions calculated =')
en = l.find('\n')
num_ions = int(l[st+23:en])
print(f'Number of ions: {num_ions}')
# If Kinchin-Pease not found...
if kp_flag is False:
warn('Warning: this code is intended to use Kinchin-Pease data.')
```
%% Output
Ion string: "Cu"
Energy: 500.0 keV
Target: " W"
Displacement: 90.0 eV
Number of ions: 500000
Ion string: "He"
Energy: 30000.0 keV
Target: " Layer 1"
Displacement: 23.0 eV
Number of ions: 200000
%% Cell type:markdown id: tags:
### Parse out the important information from RANGE.txt
%% Cell type:code id: tags:
``` python
#Find the important lines:
for i, l in enumerate(file_lines_range):
if 'DEPTH' in l: # line that defines the data columns
i_start = i + 3 # find the first line of data
# Calculate the length of the data. SRIM is always 100
length = len(file_lines_range) - i_start
data_range = np.zeros((length, 3), dtype='f4') # hold the data here
# read the data
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')
```
%% Cell type:markdown id: tags:
### Perform the calculation
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.
%% Cell type:code id: tags:
``` python
fluence = float(input('Please enter the fluence (ions/cm^2): '))
# Calculate the NRT damage in dpa
depth_in_nm = data_phonon[:, 0] / 10 # always angstroms
ion_contribution = data_phonon[:, 1]
recoil_contribution = data_phonon[:, 2]
P = ion_contribution + recoil_contribution # (P)honons
T_damage = P * fluence / density * 1e8 # 1e8 is angstrom --> cm
NRT = T_damage * 0.8 / 2.0 / displacement # NRT displacements per atom. See citation.
# Calculate the ion implantation in atomic fraction
atom_frac = data_range[:, 1] * fluence / density
```
%% Output
Please enter the fluence (ions/cm^2): 1.5e16
Please enter the fluence (ions/cm^2): 1.3e18
%% Cell type:code id: tags:
``` python
with plt.style.context('seaborn'):
fig = plt.figure(figsize=(12, 7))
gs = gridspec.GridSpec(2, 1, height_ratios=[4, 1])
ax0 = plt.subplot(gs[0])
ax0.plot(depth_in_nm, NRT, color='xkcd:brick red')
ax1 = plt.subplot(gs[1])
ax1.plot(depth_in_nm, atom_frac, color='xkcd:brick red')
```
%% Output
%% Cell type:code id: tags:
``` python
try:
err_fluence = float(input('Uncertainty in fluence (for "+/-20%" input "0.2"): '))
except ValueError:
err_fluence = 0
try:
err_displace = float(input('Uncertainty in displacement energy (for "+/-20%" input "0.2"): '))
except ValueError:
err_displace = 0
err_both = np.sqrt( err_fluence * err_fluence + err_displace * err_displace )
print(f'Combined error: {err_both}')
```
%% Output
Uncertainty in fluence (for "+/-20%" input "0.2"): 0.2
Uncertainty in displacement energy (for "+/-20%" input "0.2"): 0.2
Combined error: 0.28284271247461906
Uncertainty in fluence (for "+/-20%" input "0.2"): 0
Uncertainty in displacement energy (for "+/-20%" input "0.2"): 0
Combined error: 0.0
%% Cell type:markdown id: tags:
### 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.
%% Cell type:code id: tags:
``` python
# the 'with' block sets temporary parameters:
with plt.style.context({'font.sans-serif': 'Arial', 'svg.fonttype': 'none'}):
#Create a figure with unequal-sized subplots
fig = plt.figure(figsize=(12, 7))
gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1])
# create the first subplot
ax0 = plt.subplot(gs[0])
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
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')
ax0.legend(loc=3, prop={'size': 14}) # only show legend if error bounds present
ax0.set_ylabel('NRT damage, dpa', fontsize=20)
# Make annotations
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
fluence_string = '$' + pre + r'\times 10^{' + post + '}$ ions/m$^2$'
disp_string = str(displacement) +' eV $E_D$'
if err_both >0:
flu_err_string = r'Fluence $\pm$ ' + str(err_fluence*100) + '%'
disp_err_string = r'$E_D \pm$ ' + str(err_displace*100) + '%'
else:
flu_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}',
transform=ax0.transAxes, fontsize=14, verticalalignment='top',
horizontalalignment='right')
ax0.tick_params(axis='both', labelsize=20)
ax0.spines['right'].set_visible(False)
ax0.spines['top'].set_visible(False)
# create the smaller plot
ax1 = plt.subplot(gs[1])
if atom_frac.max() < 0.01: # Plot atom fraction if <1%
atom_string = 'Implanted\natom fraction'
atom_mult = 1
else:
atom_string = 'Implanted\natom %'
atom_mult = 100
ax1.plot(depth_in_nm, atom_frac*atom_mult, color='xkcd:brick red', label='Implantation')
if err_both > 0: # if user input error bounds
ax1.fill_between(depth_in_nm, y1=(atom_frac * (1-err_fluence))*atom_mult,
y2=(atom_frac * (1+err_fluence))*atom_mult,
color='xkcd:merlot', alpha=0.15, label='Error bounds')
ax1.legend(prop={'size': 14}) # only show legend if error bounds present
ax1.set_ylabel(atom_string, fontsize=20)
ax1.set_xlabel('Depth, nm', fontsize=20)
ax1.tick_params(axis='both', labelsize=20)
ax1.spines['right'].set_visible(False)
ax1.spines['top'].set_visible(False)
fn = input('File name (enter to skip):')
if len(fn) > 0:
fn = '..\\' + fn # go up one directory
fig.savefig(fn + '.png', dpi=300)
fig.savefig(fn + '.svg', dpi=300)
fig.savefig(fn + '.eps', dpi=300)
```
%% Output
File name (enter to skip):
File name (enter to skip):5hours
%% Cell type:code id: tags:
``` python
NRT[50]
```
%% Output
0.45283368
%% Cell type:code id: tags:
``` 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