diff --git a/Cargo.toml b/Cargo.toml
index ee9637a..1183c7b 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -24,7 +24,7 @@ rayon = "1.10.0"
geo = {version = "0.25", optional = false}
indicatif = {version = "0.18.6", features=["rayon"]}
serde = { version = "1.0.163", features = ["derive"] }
-rcpr = { git = "https://github.com/drobnyjt/rcpr", optional = true}
+rcpr = { git = "https://github.com/drobnyjt/rcpr", branch="main", optional=true}
ndarray = {version = "0.17.2", features = ["serde"], optional = true}
parry3d-f64 = {optional = true, version="0.2.0"}
pyo3 = {version = "0.29.0", optional=true}
diff --git a/examples/f_kporha_et_al_benchmark.py b/examples/f_kporha_et_al_benchmark.py
new file mode 100644
index 0000000..d54bdaa
--- /dev/null
+++ b/examples/f_kporha_et_al_benchmark.py
@@ -0,0 +1,301 @@
+from libRustBCA import *
+import numpy as np
+import matplotlib.pyplot as plt
+import sys
+import os
+#This should allow the script to find materials and formulas from anywhere
+sys.path.append(os.path.dirname(__file__)+'/../scripts')
+sys.path.append('scripts')
+import time
+from materials import *
+from tomlkit import parse, dumps
+
+def input_file(ion, target, incident_energy, angle, number_ions=1000, pot="WW", Es=None):
+
+ if ion == tungsten:
+ ion_interaction_index = 1
+ else:
+ ion_interaction_index = 0
+
+ mfp = (target["n"]/10**30)**(-1./3.)
+
+ if not Es:
+ Es = target["Es"]
+
+ if pot == "WW":
+ pot_name = "WW"
+ else:
+ pot_name = "MORSE"
+
+ 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}}
+ options = {
+ 'name': f'input_file_{ion["symbol"]}_{target["symbol"]}_{np.round(angle, 1)}_{np.round(incident_energy/1000, 4)}_{pot_name}_{Es}',
+ 'track_trajectories': False, # whether to track trajectories for plotting; memory intensive
+ 'track_recoils': True, # whether to track recoils; must enable for sputtering
+ 'track_recoil_trajectories': False, # whether to track recoil trajectories for plotting
+ 'track_displacements': False, # whether to track collisions with T > Ed for each species
+ 'track_energy_losses': False, # whether to track detailed collision energies; memory intensive
+ 'write_buffer_size': 8192, # how big the buffer is for file writing
+ 'weak_collision_order': 0, # weak collisions at radii (k + 1)*r; enable only when required
+ 'suppress_deep_recoils': False, # suppress recoils too deep to ever sputter
+ 'high_energy_free_flight_paths': False, # SRIM-style high energy free flight distances; use with caution
+ 'num_threads': 6, # number of threads to run in parallel
+ 'num_chunks': 10, # code will write to file every nth chunk; for very large simulations, increase num_chunks
+ 'electronic_stopping_mode': 'LOW_ENERGY_NONLOCAL',
+ 'mean_free_path_model': 'LIQUID', # liquid is amorphous (constant mean free path); gas is exponentially-distributed mean free paths
+ 'interaction_potential': [['ZBL', 'ZBL'],
+ ['ZBL', pot]],
+ 'scattering_integral': [
+ [{'GAUSS_MEHLER': {'n_points': 5}}, {'GAUSS_MEHLER': {'n_points': 5}}],
+ [{'GAUSS_MEHLER': {'n_points': 5}}, {'GAUSS_MEHLER': {'n_points': 5}}],
+ ],
+
+ 'root_finder': [
+ ["DEFAULTNEWTON", "DEFAULTNEWTON"],
+ ["DEFAULTNEWTON", cpr]
+ ],
+ 'seed': 0 # if <0, will generate a seed from thread-local PRNG; if >0, will be used as seed to PRNG
+ }
+
+ # material parameters are per-species
+ material_parameters = {
+ 'energy_unit': 'EV',
+ 'mass_unit': 'AMU',
+ # bulk binding energy; typically zero as a model choice
+ 'Eb': [
+ target["Eb"],
+ ],
+ # surface binding energy
+ 'Es': [
+ Es
+ ],
+ # cutoff energy - particles with E < Ec stop
+ 'Ec': [
+ target["Ec"],
+ ],
+ # atomic number
+ 'Z': [
+ target["Z"],
+ ],
+ # atomic mass
+ 'm': [
+ target["m"],
+ ],
+ # used to pick interaction potential from matrix in [options]
+ 'interaction_index': [1],
+ 'surface_binding_model': {
+ "PLANAR": {'calculation': "INDIVIDUAL"}
+ },
+ 'bulk_binding_model': 'INDIVIDUAL'
+ }
+
+ particle_parameters = {
+ 'interaction_index': [ion_interaction_index],
+ 'length_unit': 'ANGSTROM',
+ 'energy_unit': 'EV',
+ 'mass_unit': 'AMU',
+ # number of computational ions of this species to run at this energy
+ 'N': [number_ions],
+ # atomic mass
+ 'm': [ion["m"]],
+ # atomic number
+ 'Z': [ion["Z"]],
+ # incident energy
+ 'E': [incident_energy],
+ # cutoff energy - if E < Ec, particle stops
+ 'Ec': [ion["Ec"]],
+ # surface binding energy
+ 'Es': [ion["Es"]],
+ # initial position - if Es significant and E low, start (n)^(-1/3) above surface
+ # otherwise 0, 0, 0 is fine; most geometry modes have surface at x=0 with target x>0
+ 'pos': [[-7.5, 0.0, 0.0]],
+ # initial direction unit vector; most geometry modes have x-axis into the surface
+ 'dir': [
+ [
+ np.cos(angle*np.pi/180.0),
+ np.sin(angle*np.pi/180.0),
+ 0.0
+ ]
+ ],
+ }
+
+ geometry_0D = {
+ 'length_unit': 'ANGSTROM',
+ # used to correct nonlocal stopping for known compound discrpancies
+ 'electronic_stopping_correction_factor': 0.0,
+ # number densities of each species; toml densities with length_unit ANGSTROM are in 1/A^3
+ 'densities': [ target["n"]/10**30 ]
+ }
+
+ input_data = {
+ 'options': options,
+ 'material_parameters': material_parameters,
+ 'particle_parameters': particle_parameters,
+ 'geometry_input': geometry_0D
+ }
+
+ return input_data
+
+
+def tungsten_sputtering(ion, energy, angle, num_ions=1000, run_sim=False, pot="WW", Es=None):
+
+ input = input_file(ion, tungsten, energy, angle, number_ions=num_ions, pot=pot, Es=Es)
+ if run_sim: rustbca_py(input, geometry_mode="0D")
+ sputtered = np.atleast_2d(np.genfromtxt(f'{input["options"]["name"]}sputtered.output', delimiter=','))
+ if np.size(sputtered) > 0:
+ num_sputtered = np.shape(sputtered)[0]
+ else:
+ num_sputtered = 0
+ return num_sputtered/num_ions
+
+# From F. Kphorha et al.
+# picked a representative surface at random
+# digitized with WebPlotDigitizer
+w_tot= np.array([
+ [-0.11725293132328218, 0.12903225806451624],
+ [14.891122278056951, 0.2129032258064516],
+ [30.01675041876047, 0.5032258064516131],
+ [44.79061976549414, 0.9483870967741936],
+ [59.916247906197654, 1.2516129032258065],
+ [74.69011725293133, 1.1677419354838712],
+])
+w_sp = np.array([
+ [-0.23450586264656437, 0.1322580645161291],
+ [14.891122278056951, 0.19354838709677402],
+ [29.782244556113906, 0.42258064516129035],
+ [44.79061976549414, 0.7161290322580646],
+ [59.916247906197654, 0.7516129032258064],
+ [74.69011725293133, 0.1870967741935483],
+])
+ne = np.array([
+ [-0.11725293132328218, 0.3225806451612905],
+ [15.008375209380233, 0.3193548387096774],
+ [30.134003350083756, 0.3322580645161288],
+ [45.0251256281407, 0.38064516129032255],
+ [60.03350083752095, 0.2193548387096773],
+ [75.15912897822446, -0.003225806451613078],
+])
+ar = np.array([
+ [-0.4690117252931323, 0.42903225806451606],
+ [14.773869346733669, 0.44838709677419364],
+ [29.782244556113906, 0.5064516129032259],
+ [45.0251256281407, 0.5806451612903225],
+ [60.26800670016752, 0.15483870967741908],
+ [75.27638190954774, -0.003225806451613078],
+])
+energy = 200
+
+ions = [tungsten, neon, argon]
+datasets = [w_sp, ne, ar]
+colors = []
+
+num_ions = 10000
+run_sim = False
+
+tungsten["Es"] = 8.79
+tungsten["Eb"] = 0.0
+tungsten["Ec"] = 8.79
+pot_morse ={'MORSE': {'D': 2.87454*1.602e-19, 'alpha': 13.3682*1e9, 'r0': 0.238631*1e-9}}
+
+angles = np.linspace(0.0, 89.9, 12)
+for ion, dataset in zip(ions, datasets):
+
+ # keep computational cost to a minimum
+ # probably suppresses R_N
+ ion["Ec"] = 8.0
+
+ print(f'Running {ion["name"]}...')
+
+ Y_RustBCA_WW = [tungsten_sputtering(ion, energy, angle, num_ions, run_sim) for angle in angles]
+ Y_RustBCA_Morse = [tungsten_sputtering(ion, energy, angle, num_ions, run_sim, pot=pot_morse) for angle in angles]
+ Y_RustBCA_Morse_1175_eV = [tungsten_sputtering(ion, energy, angle, num_ions, run_sim, pot=pot_morse, Es=11.75) for angle in angles]
+ Y_RustBCA_KRC = [sputtering_yield(ion, tungsten, energy, angle, num_ions) for angle in angles]
+
+ color = plt.plot(angles, Y_RustBCA_WW, label=f'RustBCA {ion["symbol"]} on W - ZBL/WW', linestyle='--')[0].get_color()
+ colors.append(color)
+ plt.plot(angles, Y_RustBCA_KRC, label=f'RustBCA {ion["symbol"]} on W - Kr-C', linestyle=':', color=color)
+
+ plt.plot(angles, Y_RustBCA_Morse, label=f'RustBCA {ion["symbol"]} on W - ZBL/Morse', linestyle='-.', color=color)
+ 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)
+ plt.plot(dataset[:, 0], dataset[:, 1], label=f'MD {ion["symbol"]} on W', color=color)
+
+plt.title('Comparison to F Kporha et al., W Sputtering Yields')
+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])
+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])
+plt.scatter([0.0,], [0.123], label='Exp. W on W', marker='*', s=100, color=colors[0])
+plt.gca().set_ylim([0.0, 1.4])
+plt.legend()
+plt.xlabel('E [eV]')
+plt.ylabel('Y [at/ion]')
+
+# Produced by Somogyvari et al. (2012)
+# digitized with WebPlotDigitizer
+data_Somogyvari_Ar_W = np.array([
+ [39.95309128868367, 0.00009332543007969924],
+ [39.668850177558554, 0.0001621810097358933],
+ [54.3104899786253, 0.004365158322401661],
+ [79.29116634226808, 0.046773514128719856],
+ [104.00501083088497, 0.10964781961431856],
+ [156.2422374575878, 0.20892961308540398],
+ [202.03488119514256, 0.32359365692962827],
+ [255.71248811339458, 0.35481338923357547],
+ [303.5082793629512, 0.4365158322401659],
+])
+
+plt.figure()
+plt.title('Comparison to Somogyvari et al.: Ar on W sputtering')
+energies = np.logspace(np.log10(30), np.log10(300), 12)
+rustbca_ar_w_krc = np.array([sputtering_yield(argon, tungsten, energy, 0.0, num_ions) for energy in energies])
+rustbca_ar_w_ww = np.array([tungsten_sputtering(argon, energy, 0.0, num_ions, run_sim) for energy in energies])
+rustbca_ar_w_morse = np.array([tungsten_sputtering(argon, energy, 0.0, num_ions, run_sim, pot=pot_morse) for energy in energies])
+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])
+
+plt.loglog(data_Somogyvari_Ar_W[:, 0], data_Somogyvari_Ar_W[:, 1], linestyle='', marker='o', label='Exp. Somogyvari et al.')
+color = plt.loglog(energies, rustbca_ar_w_ww, linestyle='--', label='RustBCA, ZBL/WW potential')[0].get_color()
+plt.loglog(energies, rustbca_ar_w_krc, linestyle=':', label='RustBCA, Kr-C', color=color)
+plt.loglog(energies, rustbca_ar_w_morse, linestyle='-.', label='RustBCA, ZBL/Morse', color=color)
+plt.loglog(energies, rustbca_ar_w_morse_1175_eV, linestyle='-.', label='RustBCA, ZBL/Morse (Es=11.75eV)', color=color, marker='o')
+plt.scatter([100.0, 200.0], [0.1, 0.425], marker='*', label='F Kporha et al., Li Potential')
+plt.legend()
+plt.xlabel('E [eV]')
+plt.ylabel('Y [at/ion]')
+
+# Compiled/produced by Eckstein and Biersack (1986)
+# digitized with WebPlotDigitizer from Fig. 1c
+data_exp_Saidoh_Sone_1983 = np.array([
+ [246.28270037954246, 0.20727494578828917],
+ [495.1576962754319, 0.5491441819199058],
+ [993.8892543387418, 0.7166462062713239],
+ [1943.3723882799943, 1.1996619797357773],
+ [6947.722728427526, 3.106222121244088],
+])
+data_tridyn = np.array([
+ [68.14894695911738, 0.0007716158998594478],
+ [97.04978626836672, 0.008084212159224518],
+ [143.63389679442, 0.04173038662408122],
+ [190.63348477243179, 0.09674248249006288],
+ [301.12742785416896, 0.23046349454010032],
+ [488.4687632973695, 0.5009442817195837],
+ [692.5511880621929, 0.7839253399618061],
+ [981.6598083496964, 1.1045904715883672],
+ [1919.284303781459, 1.7777507596853939],
+ [4905.436808264407, 3.1000711459090384],
+ [9847.180039865756, 4.207988934522391],
+])
+
+energies = np.logspace(np.log10(50), np.log10(8000), 10)
+rustbca_w_w_krc = np.array([sputtering_yield(tungsten, tungsten, energy, 0.0, num_ions) for energy in energies])
+rustbca_w_w_ww = np.array([tungsten_sputtering(tungsten, energy, 0.0, num_ions, run_sim) for energy in energies])
+rustbca_w_w_morse = np.array([tungsten_sputtering(tungsten, energy, 0.0, num_ions, run_sim, pot=pot_morse) for energy in energies])
+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])
+plt.figure()
+plt.title('Comparison to Eckstein and Biersack: W on W')
+plt.loglog(energies, rustbca_w_w_krc, label='RustBCA Kr-C')
+plt.loglog(energies, rustbca_w_w_ww, label='RustBCA WW/ZBL')
+plt.loglog(energies, rustbca_w_w_morse, label='RustBCA ZBL/Morse')
+plt.loglog(energies, rustbca_w_w_morse_1175_eV, label='RustBCA ZBL/Morse (Es=11.75 eV)')
+plt.scatter(data_exp_Saidoh_Sone_1983[:, 0], data_exp_Saidoh_Sone_1983[:, 1], label='Exp. Saidoh Sone 1983')
+plt.loglog(data_tridyn[:, 0], data_tridyn[:, 1], label='TRIDYN, Eckstein and Biersack 1986')
+plt.scatter([200.0], [0.2], marker='*', label='F Kporha et al., Li Potential')
+plt.legend()
+plt.show()
\ No newline at end of file
diff --git a/examples/test_morse.py b/examples/test_morse.py
index 2716eb5..3f17c54 100644
--- a/examples/test_morse.py
+++ b/examples/test_morse.py
@@ -13,6 +13,10 @@
hydrogen['Ec'] = 0.1
hydrogen['Es'] = 1.5
+epsilon = 1e-4
+interval_limit = 1e-3
+nmax = 32
+n0=3
#This function simply contains an entire input file as a multi-line f-string to modify some inputs.
def run_morse_potential(energy, index, num_samples=10000, run_sim=True):
@@ -25,7 +29,7 @@ def run_morse_potential(energy, index, num_samples=10000, run_sim=True):
mean_free_path_model = "LIQUID"
interaction_potential = [[{{"MORSE"={{D=5.4971E-20, r0=2.782E-10, alpha=1.4198E10}}}}]]
scattering_integral = [["GAUSS_LEGENDRE"]]
- 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}}}}]]
+ 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}}}}]]
num_threads = 4
num_chunks = 10
@@ -83,11 +87,11 @@ def run_krc_morse_potential(energy, index, num_samples=10000, run_sim=True):
name = "krc_morse_{index}"
track_recoils = false
weak_collision_order = 0
- electronic_stopping_mode = "LOW_ENERGY_NONLOCAL"
+ electronic_stopping_mode = "INTERPOLATED"
mean_free_path_model = "LIQUID"
- interaction_potential = [[{{"KRC_MORSE"={{D=5.4971E-20, r0=2.782E-10, alpha=1.4198E10, k=7E10, x0=0.75E-10}}}}]]
+ interaction_potential = [[{{"KRC_MORSE"={{D=5.4971E-20, r0=2.782E-10, alpha=1.4198E10, k=8E10, x0=0.75E-10}}}}]]
scattering_integral = [["GAUSS_LEGENDRE"]]
- 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}}}}]]
+ 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}}}}]]
num_threads = 6
num_chunks = 1
@@ -170,10 +174,10 @@ def run_krc_morse_potential(energy, index, num_samples=10000, run_sim=True):
plt.semilogx(energies, r_benchmark, marker='^', linestyle='', label='Exp.')
#Running and plotting the H-Ni simulations with the Morse potential and updated Es
-num_energies = 15
+num_energies = 20
energies = np.logspace(-1, 4, num_energies)
run_sim = True
-num_samples = 100
+num_samples = 1000
R_N = np.zeros(num_energies)
R_E = np.zeros(num_energies)
R_N_2 = np.zeros(num_energies)
@@ -184,26 +188,26 @@ def run_krc_morse_potential(energy, index, num_samples=10000, run_sim=True):
R_N_2[index], R_E_2[index] = run_morse_potential(energy, index, num_samples=num_samples, run_sim=run_sim)
R_N_test = [
- 0.00, 0.01, 0.28, 0.60, 0.90,
- 0.95, 0.88, 0.81, 0.70, 0.49,
- 0.30, 0.24, 0.17, 0.11, 0.10
+ 0.0, 0.028, 0.141, 0.327, 0.65, 0.832, 0.913, 0.926,
+ 0.889, 0.844, 0.783, 0.652, 0.47, 0.38, 0.344, 0.292,
+ 0.253, 0.197, 0.139, 0.084
]
R_N_2_test = [
- 0.00, 0.01, 0.28, 0.60, 0.90,
- 0.95, 0.88, 0.81, 0.74, 0.60,
- 0.46, 0.23, 0.05, 0.00, 0.00
+ 0.0, 0.028, 0.141, 0.327, 0.65, 0.832, 0.913, 0.926,
+ 0.889, 0.844, 0.786, 0.722, 0.622, 0.531, 0.411, 0.254,
+ 0.104, 0.022, 0.003, 0.0
]
-np.testing.assert_allclose(R_N, R_N_test)
-np.testing.assert_allclose(R_N_2, R_N_2_test)
+np.testing.assert_allclose(R_N, R_N_test, atol=0.1)
+np.testing.assert_allclose(R_N_2, R_N_2_test, atol=0.1)
plt.semilogx(energies, R_N, label='R_N Morse-Kr-C H-Ni, Es=1.5eV', color='purple')
plt.semilogx(energies, R_N_2, label='R_N Morse H-Ni, Es=1.5eV', color='green')
#Plotting RustBCA data points, using the ergonomic helper function reflection_coefficient().
energies = np.logspace(-1, 4, 50)
-r_rustbca = np.array([reflection_coefficient(hydrogen, nickel, energy, 0.0, 10000) for energy in energies])
+r_rustbca = np.array([reflection_coefficient(hydrogen, nickel, energy, 0.0, 1000) for energy in energies])
r_n = r_rustbca[:, 0]
r_e = r_rustbca[:, 1]
plt.semilogx(energies, r_n, label='R_N, Default Settings', color='black')
diff --git a/examples/xenon_sputtering_benchmark.py b/examples/xenon_sputtering_benchmark.py
new file mode 100644
index 0000000..007021d
--- /dev/null
+++ b/examples/xenon_sputtering_benchmark.py
@@ -0,0 +1,238 @@
+from libRustBCA import *
+import numpy as np
+import matplotlib.pyplot as plt
+import sys
+import os
+#This should allow the script to find materials and formulas from anywhere
+sys.path.append(os.path.dirname(__file__)+'/../scripts')
+sys.path.append('scripts')
+import time
+from materials import *
+from tomlkit import parse, dumps
+
+def input_file(ion, target, incident_energy, angle, number_ions=1000):
+
+ mfp = (target["n"]/10**30)**(-1./3.)
+
+ cpr = {'CPR': {'n0': 2, 'nmax': 32, 'epsilon': 1e-3, 'complex_threshold': 1E-9, 'truncation_threshold': 1E-12, 'far_from_zero': 1e3, 'interval_limit': 1E-3, 'derivative_free': True}}
+ options = {
+ 'name': f'input_file_{ion["symbol"]}_{target["symbol"]}_{np.round(angle, 1)}_{np.round(incident_energy/1000, 4)}',
+ 'track_trajectories': False, # whether to track trajectories for plotting; memory intensive
+ 'track_recoils': True, # whether to track recoils; must enable for sputtering
+ 'track_recoil_trajectories': False, # whether to track recoil trajectories for plotting
+ 'track_displacements': False, # whether to track collisions with T > Ed for each species
+ 'track_energy_losses': False, # whether to track detailed collision energies; memory intensive
+ 'write_buffer_size': 8192, # how big the buffer is for file writing
+ 'weak_collision_order': 0, # weak collisions at radii (k + 1)*r; enable only when required
+ 'suppress_deep_recoils': False, # suppress recoils too deep to ever sputter
+ 'high_energy_free_flight_paths': False, # SRIM-style high energy free flight distances; use with caution
+ 'num_threads': 6, # number of threads to run in parallel
+ 'num_chunks': 10, # code will write to file every nth chunk; for very large simulations, increase num_chunks
+ 'electronic_stopping_mode': 'INTERPOLATED',
+ 'mean_free_path_model': 'LIQUID', # liquid is amorphous (constant mean free path); gas is exponentially-distributed mean free paths
+ 'interaction_potential': [['KR_C', 'KR_C'],
+ ['KR_C', 'KR_C']],
+ 'scattering_integral': [
+ [{'GAUSS_MEHLER': {'n_points': 5}}, {'GAUSS_MEHLER': {'n_points': 5}}],
+ [{'GAUSS_MEHLER': {'n_points': 5}}, {'GAUSS_MEHLER': {'n_points': 5}}],
+ ],
+
+ 'root_finder': [
+ ["DEFAULTNEWTON", "DEFAULTNEWTON"],
+ ["DEFAULTNEWTON", "DEFAULTNEWTON"]
+ ],
+ 'seed': 0 # if <0, will generate a seed from thread-local PRNG; if >0, will be used as seed to PRNG
+ }
+
+ # material parameters are per-species
+ material_parameters = {
+ 'energy_unit': 'EV',
+ 'mass_unit': 'AMU',
+ # bulk binding energy; typically zero as a model choice
+ 'Eb': [
+ target["Eb"],
+ ],
+ # surface binding energy
+ 'Es': [
+ target["Es"]
+ ],
+ # cutoff energy - particles with E < Ec stop
+ 'Ec': [
+ target["Ec"],
+ ],
+ # atomic number
+ 'Z': [
+ target["Z"],
+ ],
+ # atomic mass
+ 'm': [
+ target["m"],
+ ],
+ # used to pick interaction potential from matrix in [options]
+ 'interaction_index': [1],
+ 'surface_binding_model': {
+ "PLANAR": {'calculation': "INDIVIDUAL"}
+ },
+ 'bulk_binding_model': 'INDIVIDUAL'
+ }
+
+ particle_parameters = {
+ 'interaction_index': [0],
+ 'length_unit': 'ANGSTROM',
+ 'energy_unit': 'EV',
+ 'mass_unit': 'AMU',
+ # number of computational ions of this species to run at this energy
+ 'N': [number_ions],
+ # atomic mass
+ 'm': [ion["m"]],
+ # atomic number
+ 'Z': [ion["Z"]],
+ # incident energy
+ 'E': [incident_energy],
+ # cutoff energy - if E < Ec, particle stops
+ 'Ec': [ion["Ec"]],
+ # surface binding energy
+ 'Es': [ion["Es"]],
+ # initial position - if Es significant and E low, start (n)^(-1/3) above surface
+ # otherwise 0, 0, 0 is fine; most geometry modes have surface at x=0 with target x>0
+ 'pos': [[-2.*mfp, 0.0, 0.0]],
+ # initial direction unit vector; most geometry modes have x-axis into the surface
+ 'dir': [
+ [
+ np.cos(angle*np.pi/180.0),
+ np.sin(angle*np.pi/180.0),
+ 0.0
+ ]
+ ],
+ }
+
+ geometry_0D = {
+ 'length_unit': 'ANGSTROM',
+ # used to correct nonlocal stopping for known compound discrpancies
+ 'electronic_stopping_correction_factor': 1.0,
+ # number densities of each species; toml densities with length_unit ANGSTROM are in 1/A^3
+ 'densities': [ target["n"]/10**30 ]
+ }
+
+ input_data = {
+ 'options': options,
+ 'material_parameters': material_parameters,
+ 'particle_parameters': particle_parameters,
+ 'geometry_input': geometry_0D
+ }
+
+ return input_data
+
+kolasinski = np.array([
+[78.52760736196319, 0.07331378299120273],
+[98.15950920245393, 0.11290322580645173],
+[149.69325153374234, 0.22580645161290347],
+[201.2269938650307, 0.35483870967741926],
+[250.3067484662576, 0.46627565982404695],
+[299.3865030674846, 0.5835777126099708],
+[397.5460122699387, 0.7771260997067451],
+[500.61349693251526, 0.9149560117302054],
+[598.7730061349691, 1.1348973607038124],
+[699.3865030674845, 1.2580645161290325],
+[998.7730061349691, 1.7243401759530792],
+])
+tartz = np.array([
+[74.84662576687117, 0.08797653958944296],
+[98.15950920245393, 0.12023460410557174],
+[228.22085889570548, 0.30498533724340193],
+[426.9938650306747, 0.6598240469208212],
+[638.0368098159508, 0.9618768328445748],
+[840.4907975460121, 1.2023460410557185],
+[1040.490797546012, 1.4428152492668622],
+[1240.4907975460121, 1.656891495601173],
+[1440.4907975460121, 1.8005865102639296],
+])
+doerner = np.array([
+[123.92638036809814, 0.09384164222873892],
+[149.69325153374234, 0.10850439882697938],
+[174.2331288343558, 0.12023460410557196],
+[200, 0.1304985337243405],
+])
+yalin = np.array([
+[200, 0.20674486803519065],
+[250.3067484662576, 0.2903225806451615],
+[348.46625766871165, 0.4545454545454546],
+[500.61349693251526, 0.7360703812316718],
+[748.4662576687115, 0.9237536656891496],
+])
+blandino = np.array([
+[500.61349693251526, 0.43695014662756604],
+[750.920245398773, 0.6774193548387097],
+])
+zalm = np.array([
+[198.7730061349693, 0.8211143695014664],
+[500.61349693251526, 1.598240469208211],
+])
+bhattacharjee = np.array([
+[99.38650306748468, 0.14369501466275691],
+[198.7730061349693, 0.32844574780058666],
+[299.3865030674846, 0.6070381231671553],
+[402.4539877300612, 0.6832844574780059],
+[500.61349693251526, 0.7741935483870968],
+[603.6809815950919, 0.7976539589442815],
+])
+weijsenfeld = np.array([
+[200, 0.23313782991202325],
+[299.38650306748474, 0.4164222873900294],
+[399.9999999999997, 0.6011730205278589],
+[500.61349693251526, 0.7302052785923754],
+[600, 0.9002932551319645],
+[699.3865030674845, 1.0835777126099706],
+[799.9999999999998, 1.2492668621700878],
+[900.6134969325149, 1.4222873900293256],
+[1001.2269938650306, 1.5865102639296187],
+])
+rosenberg = np.array([
+[198.77300613496942, 0.3225806451612898],
+[299.3865030674846, 0.5425219941348978],
+[398.7730061349693, 0.7258064516129035],
+[599.9999999999995, 1.0747800586510257],
+])
+
+num_ions = 10000
+num_energies = 25
+run_sim = True
+energies = np.logspace(np.log10(25), np.log10(1600), num_energies)
+angle = 0.0
+
+datasets = [rosenberg, weijsenfeld, bhattacharjee, zalm, blandino, yalin, doerner, kolasinski, tartz]
+dataset_names = [
+ 'Rosenberg 1962*',
+ 'Weijsenfeld 1967*',
+ 'Bhattacharjee 1997*',
+ 'Zalm 1983*',
+ 'Blandino 1996*',
+ 'Yalin 2007*',
+ 'Doerner 2003*',
+ 'Kolasinski 2005*',
+ 'Tartz 2011',
+]
+
+Y_Xe_Mo = np.zeros(num_energies)
+molybdenum['n'] = 6.452e28
+molybdenum['Eb'] = 0.0
+molybdenum['Ec'] = 3.0
+
+for index, energy in enumerate(energies):
+
+ input_data = input_file(xenon, molybdenum, energy, angle, num_ions)
+ if run_sim: rustbca_py(input_data, geometry_mode="0D")
+ sputtered = np.genfromtxt(f'{input_data["options"]["name"]}sputtered.output', delimiter=',')
+
+ Y_Xe_Mo[index] = np.shape(sputtered)[0]/num_ions
+
+for dataset_name, dataset in zip(dataset_names, datasets):
+ plt.scatter(dataset[:, 0], dataset[:, 1], label=dataset_name)
+
+plt.plot(energies, Y_Xe_Mo, label='RustBCA Default')
+plt.gca().set_xscale('log')
+plt.legend()
+plt.xlabel('E [eV]')
+plt.ylabel('Y [at/ion]')
+plt.title('Xe on Mo Sputtering Yields')
+plt.show()
\ No newline at end of file
diff --git a/src/bca.rs b/src/bca.rs
index eb9a8ab..25082e5 100644
--- a/src/bca.rs
+++ b/src/bca.rs
@@ -1,8 +1,13 @@
use super::*;
use rand::RngExt;
+use anyhow::ensure;
#[cfg(feature = "cpr_rootfinder")]
-use rcpr::chebyshev::*;
+use rcpr::rootfinders::{
+ find_roots,
+ real_polynomial_roots,
+ Config
+};
/// Geometrical quantities of binary collision.
pub struct BinaryCollisionGeometry {
@@ -587,6 +592,16 @@ pub fn polynomial_rootfinder(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact
}
}
+#[cfg(feature = "cpr_rootfinder")]
+fn transform(x: f64) -> f64 {
+ 2.0/(x*PI/2.).tan().powi(2)
+}
+
+#[cfg(feature = "cpr_rootfinder")]
+fn inverse_transform(x: f64) -> f64 {
+ 2./PI*((2.0/x).sqrt()).atan()
+}
+
#[cfg(feature = "cpr_rootfinder")]
/// Computes the distance of closest approach of two particles with atomic numbers `Za`, `Zb` and masses `Ma`, `Mb` for an arbitrary interaction potential (e.g., Morse) for a given impact parameter and incident energy `E0` using the Chebyshev-Proxy Root-Finder method.
///
@@ -606,7 +621,7 @@ pub fn polynomial_rootfinder(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact
/// `derivative_free`: if false, use Newton's method to polish roots from the CPR. If true, use the secant method.
///
/// # Returns
-/// Returns the distance of closest approach or an error if the root-finder failed.
+/// Returns the distance of closest approach (reduced by a) or an error if the root-finder failed.
pub fn cpr_rootfinder(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parameter: f64,
interaction_potential: InteractionPotential, n0: usize, nmax: usize, epsilon: f64,
complex_threshold: f64, truncation_threshold: f64, far_from_zero: f64,
@@ -614,32 +629,35 @@ pub fn cpr_rootfinder(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parame
//Lindhard screening length and reduced energy
let a = interactions::screening_length(Za, Zb, interaction_potential);
- let reduced_energy = LINDHARD_REDUCED_ENERGY_PREFACTOR*a*Mb/(Ma+Mb)/Za/Zb*E0;
let relative_energy = E0*Mb/(Ma + Mb);
- let p = impact_parameter;
- let f = |r: f64| -> f64 {interactions::distance_of_closest_approach_function(r, a, Za, Zb, relative_energy, impact_parameter, interaction_potential)};
- let g = |r: f64| -> f64 {interactions::distance_of_closest_approach_function_singularity_free(r, a, Za, Zb, relative_energy, impact_parameter, interaction_potential)*
- interactions::scaling_function(r, impact_parameter, interaction_potential)};
-
- let upper_bound = impact_parameter + interactions::crossing_point_doca(interaction_potential);
- let lower_bound = impact_parameter / 1000.0;
-
- let roots = match derivative_free {
- true => find_roots_with_secant_polishing(&g, &f, lower_bound, upper_bound,
- n0, epsilon, nmax, complex_threshold,
- truncation_threshold, interval_limit, far_from_zero),
-
- false => {
- let df = |r: f64| -> f64 {interactions::diff_distance_of_closest_approach_function(r, a, Za, Zb, relative_energy, impact_parameter, interaction_potential)};
- find_roots_with_newton_polishing(&g, &f, &df, lower_bound, upper_bound,
- n0, epsilon, nmax, complex_threshold,
- truncation_threshold, interval_limit, far_from_zero)
- }
- }.with_context(|| format!("Numerical error: CPR Rootfinder failed to converge when calculating distance of closest approach for Er = {} eV p = {} A using {}.",
- relative_energy/EV, impact_parameter/ANGSTROM, interaction_potential))?;
+ let g = |r: f64| -> f64 {
+ interactions::distance_of_closest_approach_function_singularity_free(transform(r)*a, a, Za, Zb, relative_energy, impact_parameter, interaction_potential)*
+ interactions::scaling_function(transform(r)*a, a, interaction_potential)
+ };
+
+ ensure!(1.0 - inverse_transform(impact_parameter/a) > interval_limit, "Numerical error: impact parameter {} A smaller than interval limit.", impact_parameter/a);
+
+ let upper_bound = 1.0;
+ let lower_bound = 1e-4;
+
+ let delta = 1e-5;
+ let config = Config::new(
+ epsilon,
+ delta,
+ n0,
+ nmax,
+ complex_threshold,
+ truncation_threshold,
+ far_from_zero,
+ interval_limit
+ );
+
+ let roots = find_roots(&g, vec![(lower_bound, upper_bound)], config)?;
- let max_root = roots.iter().cloned().fold(f64::NAN, f64::max)/a;
+ // Since above the arg to doca is transform(r)*a, this is already scaled as output
+ //let max_root = roots.iter().map(|&x| transform(x)).fold(f64::NAN, f64::max);
+ let max_root = roots.iter().map(|&x| transform(x)).max_by(f64::total_cmp).expect("Numerical error: failed to find maximum root.");
if roots.is_empty() || max_root.is_nan() {
return Err(anyhow!("Numerical error: CPR rootfinder failed to find root. x0: {}, F(a): {}, F(b): {};", max_root, g(0.), g(upper_bound)));
diff --git a/src/interactions.rs b/src/interactions.rs
index 5ba58ee..5ffb125 100644
--- a/src/interactions.rs
+++ b/src/interactions.rs
@@ -179,7 +179,7 @@ pub fn scaling_function(r: f64, a: f64, interaction_potential: InteractionPotent
1./(1. + (r/ANGSTROM).powi(2))
},
InteractionPotential::KRC_MORSE{D, alpha, r0, k, x0} => {
- 1./(1. + (r*alpha).powi(2))
+ 1.
}
InteractionPotential::COULOMB{..} => panic!("Coulombic potential cannot be used with rootfinder.")
}
@@ -282,11 +282,11 @@ pub fn dphi(xi: f64, interaction_potential: InteractionPotential) -> f64 {
pub fn screening_length(Za: f64, Zb: f64, interaction_potential: InteractionPotential) -> f64 {
match interaction_potential {
//ZBL screening length, Eckstein (4.1.8)
- InteractionPotential::ZBL => zbl_screening_length_lookup(Za as u64, Zb as u64),
+ InteractionPotential::ZBL | InteractionPotential::WW => zbl_screening_length_lookup(Za as u64, Zb as u64),
//Lindhard/Firsov screening length, Eckstein (4.1.5)
- InteractionPotential::MOLIERE | InteractionPotential::KR_C | InteractionPotential::LENZ_JENSEN | InteractionPotential::TRIDYN | InteractionPotential::WW => lindhard_screening_length_lookup(Za as u64, Zb as u64),
+ InteractionPotential::MOLIERE | InteractionPotential::KR_C | InteractionPotential::LENZ_JENSEN | InteractionPotential::TRIDYN => lindhard_screening_length_lookup(Za as u64, Zb as u64),
InteractionPotential::LENNARD_JONES_12_6{..} | InteractionPotential::LENNARD_JONES_65_6{..} => lindhard_screening_length_lookup(Za as u64, Zb as u64),
- InteractionPotential::MORSE{D, alpha, r0} => alpha,
+ InteractionPotential::MORSE{D, alpha, r0} => 1./alpha,
InteractionPotential::COULOMB{Za: Z1, Zb: Z2} => zbl_screening_length_lookup(Za as u64, Zb as u64),
InteractionPotential::KRC_MORSE{..} => lindhard_screening_length_lookup(Za as u64, Zb as u64),
InteractionPotential::FOUR_EIGHT{..} => lindhard_screening_length_lookup(Za as u64, Zb as u64),
@@ -434,7 +434,13 @@ pub fn doca_morse(r: f64, impact_parameter: f64, relative_energy: f64, D: f64, a
/// Distance of closest approach function for Morse potential.
pub fn doca_krc_morse(r: f64, impact_parameter: f64, relative_energy: f64, a: f64, Za: f64, Zb: f64, D: f64, alpha: f64, r0: f64, k: f64, x0: f64) -> f64 {
- (r*alpha).powi(2) - (r*alpha).powi(2)/relative_energy*krc_morse(r, a, Za, Zb, D, alpha, r0, k, x0) - (impact_parameter*alpha).powi(2)
+ let K = coulomb_constant(Za, Zb);
+ let ralpha = r*alpha;
+ let term_1 = (ralpha).powi(2) - (impact_parameter*alpha).powi(2);
+ let term_2 = -(ralpha)*alpha*(K/relative_energy)*phi(r/a, InteractionPotential::KR_C)*smootherstep(r, -k, x0);
+ let term_3 = -(ralpha).powi(2)*(morse(r, D, alpha, r0)/relative_energy)*smootherstep(r, k, x0);
+ let scale = 1./(1. + ralpha).powi(2);
+ term_1*scale + term_2*scale + term_3*scale
}
/// Distance of closest approach function for LJ 6.5-6 potential.
diff --git a/src/lib.rs b/src/lib.rs
index 4ebb80f..6b0f24e 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -54,6 +54,7 @@ use pyo3::exceptions::{PyValueError, PyRuntimeError};
//Load internal modules
pub mod material;
pub mod particle;
+#[cfg(test)]
pub mod tests;
pub mod interactions;
pub mod bca;
diff --git a/src/main.rs b/src/main.rs
index 9422bd3..f9f0509 100644
--- a/src/main.rs
+++ b/src/main.rs
@@ -11,10 +11,6 @@ use anyhow::{Result, Context, anyhow};
//Serializing/Deserializing crate
use serde::*;
-//Parallelization
-//use rayon::prelude::*;
-//use rayon::ThreadPoolBuilder;
-
//I/O
use std::fs::OpenOptions;
use std::io::prelude::*;
@@ -45,7 +41,6 @@ pub mod structs;
pub mod sphere;
pub mod physics;
pub mod math;
-
#[cfg(feature = "parry3d")]
pub mod parry;
@@ -67,7 +62,6 @@ macro_rules! main_loop {
($geometry_type:ident, $input_file:expr) => {
{
let (particle_input_array, material, options, output_units) = input::input::<$geometry_type>($input_file);
- //Initialize threads with rayon
println!("Processing {} ions...", particle_input_array.len());
println!("Initializing with {} threads...", options.num_threads);
let _ = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build_global();
@@ -98,7 +92,7 @@ fn main() {
_ => panic!("Too many command line arguments. RustBCA accepts 0 (use 'input.toml') 1 () or 2 ( )"),
};
- // This invokes the above macro that expands into the physics loop invocation for each type
+ // This invokes the above macro that expands into the physics loop invocation for each geometry type
match geometry_type {
GeometryType::MESH0D => main_loop!(Mesh0D, input_file),
GeometryType::MESH1D => main_loop!(Mesh1D, input_file),
diff --git a/src/tests.rs b/src/tests.rs
index afba5de..3c2b4b4 100644
--- a/src/tests.rs
+++ b/src/tests.rs
@@ -4,6 +4,8 @@ use super::*;
use float_cmp::*;
#[cfg(test)]
use rand::RngExt;
+#[cfg(feature = "cpr_rootfinder")]
+use rcpr::rootfinders::*;
#[test]
#[cfg(feature = "cpr_rootfinder")]