Skip to content

Commit e1e006f

Browse files
drobnyjtdrobnyjt
authored
Merge pull request #339 from lcpp-org/new_rcpr
New rcpr
2 parents d42fd1e + db0a527 commit e1e006f

9 files changed

Lines changed: 623 additions & 54 deletions

File tree

‎Cargo.toml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -24,7 +24,7 @@ rayon = "1.10.0"
2424
geo = {version = "0.25", optional = false}
2525
indicatif = {version = "0.18.6", features=["rayon"]}
2626
serde = { version = "1.0.163", features = ["derive"] }
27-
rcpr = { git = "https://github.com/drobnyjt/rcpr", optional = true}
27+
rcpr = { git = "https://github.com/drobnyjt/rcpr", branch="main", optional=true}
2828
ndarray = {version = "0.17.2", features = ["serde"], optional = true}
2929
parry3d-f64 = {optional = true, version="0.2.0"}
3030
pyo3 = {version = "0.29.0", optional=true}

‎examples/f_kporha_et_al_benchmark.py‎

Lines changed: 301 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,301 @@
1+
from libRustBCA import *
2+
import numpy as np
3+
import matplotlib.pyplot as plt
4+
import sys
5+
import os
6+
#This should allow the script to find materials and formulas from anywhere
7+
sys.path.append(os.path.dirname(__file__)+'/../scripts')
8+
sys.path.append('scripts')
9+
import time
10+
from materials import *
11+
from tomlkit import parse, dumps
12+
13+
def input_file(ion, target, incident_energy, angle, number_ions=1000, pot="WW", Es=None):
14+
15+
if ion == tungsten:
16+
ion_interaction_index = 1
17+
else:
18+
ion_interaction_index = 0
19+
20+
mfp = (target["n"]/10**30)**(-1./3.)
21+
22+
if not Es:
23+
Es = target["Es"]
24+
25+
if pot == "WW":
26+
pot_name = "WW"
27+
else:
28+
pot_name = "MORSE"
29+
30+
cpr = {'CPR': {'n0': 2, 'nmax': 32, 'epsilon': 5e-4, 'complex_threshold': 1E-9, 'truncation_threshold': 1E-9, 'far_from_zero': 1e3, 'interval_limit': 1E-4, 'derivative_free': True}}
31+
options = {
32+
'name': f'input_file_{ion["symbol"]}_{target["symbol"]}_{np.round(angle, 1)}_{np.round(incident_energy/1000, 4)}_{pot_name}_{Es}',
33+
'track_trajectories': False, # whether to track trajectories for plotting; memory intensive
34+
'track_recoils': True, # whether to track recoils; must enable for sputtering
35+
'track_recoil_trajectories': False, # whether to track recoil trajectories for plotting
36+
'track_displacements': False, # whether to track collisions with T > Ed for each species
37+
'track_energy_losses': False, # whether to track detailed collision energies; memory intensive
38+
'write_buffer_size': 8192, # how big the buffer is for file writing
39+
'weak_collision_order': 0, # weak collisions at radii (k + 1)*r; enable only when required
40+
'suppress_deep_recoils': False, # suppress recoils too deep to ever sputter
41+
'high_energy_free_flight_paths': False, # SRIM-style high energy free flight distances; use with caution
42+
'num_threads': 6, # number of threads to run in parallel
43+
'num_chunks': 10, # code will write to file every nth chunk; for very large simulations, increase num_chunks
44+
'electronic_stopping_mode': 'LOW_ENERGY_NONLOCAL',
45+
'mean_free_path_model': 'LIQUID', # liquid is amorphous (constant mean free path); gas is exponentially-distributed mean free paths
46+
'interaction_potential': [['ZBL', 'ZBL'],
47+
['ZBL', pot]],
48+
'scattering_integral': [
49+
[{'GAUSS_MEHLER': {'n_points': 5}}, {'GAUSS_MEHLER': {'n_points': 5}}],
50+
[{'GAUSS_MEHLER': {'n_points': 5}}, {'GAUSS_MEHLER': {'n_points': 5}}],
51+
],
52+
53+
'root_finder': [
54+
["DEFAULTNEWTON", "DEFAULTNEWTON"],
55+
["DEFAULTNEWTON", cpr]
56+
],
57+
'seed': 0 # if <0, will generate a seed from thread-local PRNG; if >0, will be used as seed to PRNG
58+
}
59+
60+
# material parameters are per-species
61+
material_parameters = {
62+
'energy_unit': 'EV',
63+
'mass_unit': 'AMU',
64+
# bulk binding energy; typically zero as a model choice
65+
'Eb': [
66+
target["Eb"],
67+
],
68+
# surface binding energy
69+
'Es': [
70+
Es
71+
],
72+
# cutoff energy - particles with E < Ec stop
73+
'Ec': [
74+
target["Ec"],
75+
],
76+
# atomic number
77+
'Z': [
78+
target["Z"],
79+
],
80+
# atomic mass
81+
'm': [
82+
target["m"],
83+
],
84+
# used to pick interaction potential from matrix in [options]
85+
'interaction_index': [1],
86+
'surface_binding_model': {
87+
"PLANAR": {'calculation': "INDIVIDUAL"}
88+
},
89+
'bulk_binding_model': 'INDIVIDUAL'
90+
}
91+
92+
particle_parameters = {
93+
'interaction_index': [ion_interaction_index],
94+
'length_unit': 'ANGSTROM',
95+
'energy_unit': 'EV',
96+
'mass_unit': 'AMU',
97+
# number of computational ions of this species to run at this energy
98+
'N': [number_ions],
99+
# atomic mass
100+
'm': [ion["m"]],
101+
# atomic number
102+
'Z': [ion["Z"]],
103+
# incident energy
104+
'E': [incident_energy],
105+
# cutoff energy - if E < Ec, particle stops
106+
'Ec': [ion["Ec"]],
107+
# surface binding energy
108+
'Es': [ion["Es"]],
109+
# initial position - if Es significant and E low, start (n)^(-1/3) above surface
110+
# otherwise 0, 0, 0 is fine; most geometry modes have surface at x=0 with target x>0
111+
'pos': [[-7.5, 0.0, 0.0]],
112+
# initial direction unit vector; most geometry modes have x-axis into the surface
113+
'dir': [
114+
[
115+
np.cos(angle*np.pi/180.0),
116+
np.sin(angle*np.pi/180.0),
117+
0.0
118+
]
119+
],
120+
}
121+
122+
geometry_0D = {
123+
'length_unit': 'ANGSTROM',
124+
# used to correct nonlocal stopping for known compound discrpancies
125+
'electronic_stopping_correction_factor': 0.0,
126+
# number densities of each species; toml densities with length_unit ANGSTROM are in 1/A^3
127+
'densities': [ target["n"]/10**30 ]
128+
}
129+
130+
input_data = {
131+
'options': options,
132+
'material_parameters': material_parameters,
133+
'particle_parameters': particle_parameters,
134+
'geometry_input': geometry_0D
135+
}
136+
137+
return input_data
138+
139+
140+
def tungsten_sputtering(ion, energy, angle, num_ions=1000, run_sim=False, pot="WW", Es=None):
141+
142+
input = input_file(ion, tungsten, energy, angle, number_ions=num_ions, pot=pot, Es=Es)
143+
if run_sim: rustbca_py(input, geometry_mode="0D")
144+
sputtered = np.atleast_2d(np.genfromtxt(f'{input["options"]["name"]}sputtered.output', delimiter=','))
145+
if np.size(sputtered) > 0:
146+
num_sputtered = np.shape(sputtered)[0]
147+
else:
148+
num_sputtered = 0
149+
return num_sputtered/num_ions
150+
151+
# From F. Kphorha et al.
152+
# picked a representative surface at random
153+
# digitized with WebPlotDigitizer
154+
w_tot= np.array([
155+
[-0.11725293132328218, 0.12903225806451624],
156+
[14.891122278056951, 0.2129032258064516],
157+
[30.01675041876047, 0.5032258064516131],
158+
[44.79061976549414, 0.9483870967741936],
159+
[59.916247906197654, 1.2516129032258065],
160+
[74.69011725293133, 1.1677419354838712],
161+
])
162+
w_sp = np.array([
163+
[-0.23450586264656437, 0.1322580645161291],
164+
[14.891122278056951, 0.19354838709677402],
165+
[29.782244556113906, 0.42258064516129035],
166+
[44.79061976549414, 0.7161290322580646],
167+
[59.916247906197654, 0.7516129032258064],
168+
[74.69011725293133, 0.1870967741935483],
169+
])
170+
ne = np.array([
171+
[-0.11725293132328218, 0.3225806451612905],
172+
[15.008375209380233, 0.3193548387096774],
173+
[30.134003350083756, 0.3322580645161288],
174+
[45.0251256281407, 0.38064516129032255],
175+
[60.03350083752095, 0.2193548387096773],
176+
[75.15912897822446, -0.003225806451613078],
177+
])
178+
ar = np.array([
179+
[-0.4690117252931323, 0.42903225806451606],
180+
[14.773869346733669, 0.44838709677419364],
181+
[29.782244556113906, 0.5064516129032259],
182+
[45.0251256281407, 0.5806451612903225],
183+
[60.26800670016752, 0.15483870967741908],
184+
[75.27638190954774, -0.003225806451613078],
185+
])
186+
energy = 200
187+
188+
ions = [tungsten, neon, argon]
189+
datasets = [w_sp, ne, ar]
190+
colors = []
191+
192+
num_ions = 10000
193+
run_sim = False
194+
195+
tungsten["Es"] = 8.79
196+
tungsten["Eb"] = 0.0
197+
tungsten["Ec"] = 8.79
198+
pot_morse ={'MORSE': {'D': 2.87454*1.602e-19, 'alpha': 13.3682*1e9, 'r0': 0.238631*1e-9}}
199+
200+
angles = np.linspace(0.0, 89.9, 12)
201+
for ion, dataset in zip(ions, datasets):
202+
203+
# keep computational cost to a minimum
204+
# probably suppresses R_N
205+
ion["Ec"] = 8.0
206+
207+
print(f'Running {ion["name"]}...')
208+
209+
Y_RustBCA_WW = [tungsten_sputtering(ion, energy, angle, num_ions, run_sim) for angle in angles]
210+
Y_RustBCA_Morse = [tungsten_sputtering(ion, energy, angle, num_ions, run_sim, pot=pot_morse) for angle in angles]
211+
Y_RustBCA_Morse_1175_eV = [tungsten_sputtering(ion, energy, angle, num_ions, run_sim, pot=pot_morse, Es=11.75) for angle in angles]
212+
Y_RustBCA_KRC = [sputtering_yield(ion, tungsten, energy, angle, num_ions) for angle in angles]
213+
214+
color = plt.plot(angles, Y_RustBCA_WW, label=f'RustBCA {ion["symbol"]} on W - ZBL/WW', linestyle='--')[0].get_color()
215+
colors.append(color)
216+
plt.plot(angles, Y_RustBCA_KRC, label=f'RustBCA {ion["symbol"]} on W - Kr-C', linestyle=':', color=color)
217+
218+
plt.plot(angles, Y_RustBCA_Morse, label=f'RustBCA {ion["symbol"]} on W - ZBL/Morse', linestyle='-.', color=color)
219+
plt.plot(angles, Y_RustBCA_Morse_1175_eV, label=f'RustBCA {ion["symbol"]} on W - ZBL/Morse (Es=11.75 eV)', linestyle='-.', marker='o', color=color)
220+
plt.plot(dataset[:, 0], dataset[:, 1], label=f'MD {ion["symbol"]} on W', color=color)
221+
222+
plt.title('Comparison to F Kporha et al., W Sputtering Yields')
223+
plt.scatter([0.0, 0.0], [0.13, 0.2], label='Exp. Ne on W (Laegreid et al., Stuart et al.)', marker='*', s=100, color=colors[1])
224+
plt.scatter([0.0, 0.0, 0.0], [0.29, 0.6, 0.32], label='Exp. Ar on W (Laegreid et al., Stuart et al., Somogyvári et al.)', marker='*', s=100, color=colors[2])
225+
plt.scatter([0.0,], [0.123], label='Exp. W on W', marker='*', s=100, color=colors[0])
226+
plt.gca().set_ylim([0.0, 1.4])
227+
plt.legend()
228+
plt.xlabel('E [eV]')
229+
plt.ylabel('Y [at/ion]')
230+
231+
# Produced by Somogyvari et al. (2012)
232+
# digitized with WebPlotDigitizer
233+
data_Somogyvari_Ar_W = np.array([
234+
[39.95309128868367, 0.00009332543007969924],
235+
[39.668850177558554, 0.0001621810097358933],
236+
[54.3104899786253, 0.004365158322401661],
237+
[79.29116634226808, 0.046773514128719856],
238+
[104.00501083088497, 0.10964781961431856],
239+
[156.2422374575878, 0.20892961308540398],
240+
[202.03488119514256, 0.32359365692962827],
241+
[255.71248811339458, 0.35481338923357547],
242+
[303.5082793629512, 0.4365158322401659],
243+
])
244+
245+
plt.figure()
246+
plt.title('Comparison to Somogyvari et al.: Ar on W sputtering')
247+
energies = np.logspace(np.log10(30), np.log10(300), 12)
248+
rustbca_ar_w_krc = np.array([sputtering_yield(argon, tungsten, energy, 0.0, num_ions) for energy in energies])
249+
rustbca_ar_w_ww = np.array([tungsten_sputtering(argon, energy, 0.0, num_ions, run_sim) for energy in energies])
250+
rustbca_ar_w_morse = np.array([tungsten_sputtering(argon, energy, 0.0, num_ions, run_sim, pot=pot_morse) for energy in energies])
251+
rustbca_ar_w_morse_1175_eV = np.array([tungsten_sputtering(argon, energy, 0.0, num_ions, run_sim, pot=pot_morse, Es=11.75) for energy in energies])
252+
253+
plt.loglog(data_Somogyvari_Ar_W[:, 0], data_Somogyvari_Ar_W[:, 1], linestyle='', marker='o', label='Exp. Somogyvari et al.')
254+
color = plt.loglog(energies, rustbca_ar_w_ww, linestyle='--', label='RustBCA, ZBL/WW potential')[0].get_color()
255+
plt.loglog(energies, rustbca_ar_w_krc, linestyle=':', label='RustBCA, Kr-C', color=color)
256+
plt.loglog(energies, rustbca_ar_w_morse, linestyle='-.', label='RustBCA, ZBL/Morse', color=color)
257+
plt.loglog(energies, rustbca_ar_w_morse_1175_eV, linestyle='-.', label='RustBCA, ZBL/Morse (Es=11.75eV)', color=color, marker='o')
258+
plt.scatter([100.0, 200.0], [0.1, 0.425], marker='*', label='F Kporha et al., Li Potential')
259+
plt.legend()
260+
plt.xlabel('E [eV]')
261+
plt.ylabel('Y [at/ion]')
262+
263+
# Compiled/produced by Eckstein and Biersack (1986)
264+
# digitized with WebPlotDigitizer from Fig. 1c
265+
data_exp_Saidoh_Sone_1983 = np.array([
266+
[246.28270037954246, 0.20727494578828917],
267+
[495.1576962754319, 0.5491441819199058],
268+
[993.8892543387418, 0.7166462062713239],
269+
[1943.3723882799943, 1.1996619797357773],
270+
[6947.722728427526, 3.106222121244088],
271+
])
272+
data_tridyn = np.array([
273+
[68.14894695911738, 0.0007716158998594478],
274+
[97.04978626836672, 0.008084212159224518],
275+
[143.63389679442, 0.04173038662408122],
276+
[190.63348477243179, 0.09674248249006288],
277+
[301.12742785416896, 0.23046349454010032],
278+
[488.4687632973695, 0.5009442817195837],
279+
[692.5511880621929, 0.7839253399618061],
280+
[981.6598083496964, 1.1045904715883672],
281+
[1919.284303781459, 1.7777507596853939],
282+
[4905.436808264407, 3.1000711459090384],
283+
[9847.180039865756, 4.207988934522391],
284+
])
285+
286+
energies = np.logspace(np.log10(50), np.log10(8000), 10)
287+
rustbca_w_w_krc = np.array([sputtering_yield(tungsten, tungsten, energy, 0.0, num_ions) for energy in energies])
288+
rustbca_w_w_ww = np.array([tungsten_sputtering(tungsten, energy, 0.0, num_ions, run_sim) for energy in energies])
289+
rustbca_w_w_morse = np.array([tungsten_sputtering(tungsten, energy, 0.0, num_ions, run_sim, pot=pot_morse) for energy in energies])
290+
rustbca_w_w_morse_1175_eV = np.array([tungsten_sputtering(tungsten, energy, 0.0, num_ions, run_sim, pot=pot_morse, Es=11.75) for energy in energies])
291+
plt.figure()
292+
plt.title('Comparison to Eckstein and Biersack: W on W')
293+
plt.loglog(energies, rustbca_w_w_krc, label='RustBCA Kr-C')
294+
plt.loglog(energies, rustbca_w_w_ww, label='RustBCA WW/ZBL')
295+
plt.loglog(energies, rustbca_w_w_morse, label='RustBCA ZBL/Morse')
296+
plt.loglog(energies, rustbca_w_w_morse_1175_eV, label='RustBCA ZBL/Morse (Es=11.75 eV)')
297+
plt.scatter(data_exp_Saidoh_Sone_1983[:, 0], data_exp_Saidoh_Sone_1983[:, 1], label='Exp. Saidoh Sone 1983')
298+
plt.loglog(data_tridyn[:, 0], data_tridyn[:, 1], label='TRIDYN, Eckstein and Biersack 1986')
299+
plt.scatter([200.0], [0.2], marker='*', label='F Kporha et al., Li Potential')
300+
plt.legend()
301+
plt.show()

‎examples/test_morse.py‎

Lines changed: 19 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -13,6 +13,10 @@
1313
hydrogen['Ec'] = 0.1
1414
hydrogen['Es'] = 1.5
1515

16+
epsilon = 1e-4
17+
interval_limit = 1e-3
18+
nmax = 32
19+
n0=3
1620
#This function simply contains an entire input file as a multi-line f-string to modify some inputs.
1721
def run_morse_potential(energy, index, num_samples=10000, run_sim=True):
1822

@@ -25,7 +29,7 @@ def run_morse_potential(energy, index, num_samples=10000, run_sim=True):
2529
mean_free_path_model = "LIQUID"
2630
interaction_potential = [[{{"MORSE"={{D=5.4971E-20, r0=2.782E-10, alpha=1.4198E10}}}}]]
2731
scattering_integral = [["GAUSS_LEGENDRE"]]
28-
root_finder = [[{{"CPR"={{n0=3, nmax=100, epsilon=1E-9, complex_threshold=1E-9, truncation_threshold=1E-9, far_from_zero=1E9, interval_limit=1E-13, derivative_free=true}}}}]]
32+
root_finder = [[{{"CPR"={{n0={n0}, nmax={nmax}, epsilon={epsilon}, complex_threshold=1E-9, truncation_threshold=1E-9, far_from_zero=1E22, interval_limit={interval_limit}, derivative_free=true}}}}]]
2933
num_threads = 4
3034
num_chunks = 10
3135
@@ -83,11 +87,11 @@ def run_krc_morse_potential(energy, index, num_samples=10000, run_sim=True):
8387
name = "krc_morse_{index}"
8488
track_recoils = false
8589
weak_collision_order = 0
86-
electronic_stopping_mode = "LOW_ENERGY_NONLOCAL"
90+
electronic_stopping_mode = "INTERPOLATED"
8791
mean_free_path_model = "LIQUID"
88-
interaction_potential = [[{{"KRC_MORSE"={{D=5.4971E-20, r0=2.782E-10, alpha=1.4198E10, k=7E10, x0=0.75E-10}}}}]]
92+
interaction_potential = [[{{"KRC_MORSE"={{D=5.4971E-20, r0=2.782E-10, alpha=1.4198E10, k=8E10, x0=0.75E-10}}}}]]
8993
scattering_integral = [["GAUSS_LEGENDRE"]]
90-
root_finder = [[{{"CPR"={{n0=2, nmax=200, epsilon=1E-9, complex_threshold=1E-9, truncation_threshold=1E-9, far_from_zero=1E9, interval_limit=1E-13, derivative_free=true}}}}]]
94+
root_finder = [[{{"CPR"={{n0=3, nmax={nmax}, epsilon={epsilon}, complex_threshold=1E-9, truncation_threshold=1E-9, far_from_zero=1E22, interval_limit={interval_limit}, derivative_free=true}}}}]]
9195
num_threads = 6
9296
num_chunks = 1
9397
@@ -170,10 +174,10 @@ def run_krc_morse_potential(energy, index, num_samples=10000, run_sim=True):
170174
plt.semilogx(energies, r_benchmark, marker='^', linestyle='', label='Exp.')
171175

172176
#Running and plotting the H-Ni simulations with the Morse potential and updated Es
173-
num_energies = 15
177+
num_energies = 20
174178
energies = np.logspace(-1, 4, num_energies)
175179
run_sim = True
176-
num_samples = 100
180+
num_samples = 1000
177181
R_N = np.zeros(num_energies)
178182
R_E = np.zeros(num_energies)
179183
R_N_2 = np.zeros(num_energies)
@@ -184,26 +188,26 @@ def run_krc_morse_potential(energy, index, num_samples=10000, run_sim=True):
184188
R_N_2[index], R_E_2[index] = run_morse_potential(energy, index, num_samples=num_samples, run_sim=run_sim)
185189

186190
R_N_test = [
187-
0.00, 0.01, 0.28, 0.60, 0.90,
188-
0.95, 0.88, 0.81, 0.70, 0.49,
189-
0.30, 0.24, 0.17, 0.11, 0.10
191+
0.0, 0.028, 0.141, 0.327, 0.65, 0.832, 0.913, 0.926,
192+
0.889, 0.844, 0.783, 0.652, 0.47, 0.38, 0.344, 0.292,
193+
0.253, 0.197, 0.139, 0.084
190194
]
191195

192196
R_N_2_test = [
193-
0.00, 0.01, 0.28, 0.60, 0.90,
194-
0.95, 0.88, 0.81, 0.74, 0.60,
195-
0.46, 0.23, 0.05, 0.00, 0.00
197+
0.0, 0.028, 0.141, 0.327, 0.65, 0.832, 0.913, 0.926,
198+
0.889, 0.844, 0.786, 0.722, 0.622, 0.531, 0.411, 0.254,
199+
0.104, 0.022, 0.003, 0.0
196200
]
197201

198-
np.testing.assert_allclose(R_N, R_N_test)
199-
np.testing.assert_allclose(R_N_2, R_N_2_test)
202+
np.testing.assert_allclose(R_N, R_N_test, atol=0.1)
203+
np.testing.assert_allclose(R_N_2, R_N_2_test, atol=0.1)
200204

201205
plt.semilogx(energies, R_N, label='R_N Morse-Kr-C H-Ni, Es=1.5eV', color='purple')
202206
plt.semilogx(energies, R_N_2, label='R_N Morse H-Ni, Es=1.5eV', color='green')
203207

204208
#Plotting RustBCA data points, using the ergonomic helper function reflection_coefficient().
205209
energies = np.logspace(-1, 4, 50)
206-
r_rustbca = np.array([reflection_coefficient(hydrogen, nickel, energy, 0.0, 10000) for energy in energies])
210+
r_rustbca = np.array([reflection_coefficient(hydrogen, nickel, energy, 0.0, 1000) for energy in energies])
207211
r_n = r_rustbca[:, 0]
208212
r_e = r_rustbca[:, 1]
209213
plt.semilogx(energies, r_n, label='R_N, Default Settings', color='black')

0 commit comments

Comments
 (0)