Update RustBCA to 2.0.0: Improvements to Python bindings, code cleanup, simplify Mesh2D input

.github/workflows/rustbca_compile_check.yml

@@ -44,7 +44,7 @@ jobs:
         python3 -m pip install setuptools_rust testresources setuptools wheel
         python3 -m pip install --upgrade pip setuptools wheel
         python3 -m pip install .
-        python3 -c "from libRustBCA.pybca import *;"
+        python3 -c "from libRustBCA import *;"
         python3 examples/test_rustbca.py
     - name: Test Fortran and C bindings
       run : |
@@ -65,9 +65,9 @@ jobs:
         ./target/release/RustBCA 0D examples/titanium_dioxide_0D.toml
         ./target/release/RustBCA 1D examples/layered_geometry_1D.toml
         cat 2000.0eV_0.0001deg_He_TiO2_Al_Sisummary.output
-        ./target/release/RustBCA examples/boron_nitride.toml
         ./target/release/RustBCA examples/layered_geometry.toml
         cat 2000.0eV_0.0001deg_He_TiO2_Al_Sisummary.output
         ./target/release/RustBCA SPHERE examples/boron_nitride_sphere.toml
         cargo run --release --features parry3d TRIMESH examples/tungsten_twist_trimesh.toml
+        ./target/release/RustBCA examples/boron_nitride_wire.toml
 

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "RustBCA"
-version = "1.3.0"
+version = "2.0.0"
 default-run = "RustBCA"
 authors = ["Jon Drobny <drobny2@illinois.edu>", "Jon Drobny <jdrobny@tae.com>"]
 edition = "2018"
@@ -12,7 +12,7 @@ path = "src/main.rs"
 [lib]
 name = "libRustBCA"
 path = "src/lib.rs"
-crate-type = ["cdylib"]
+crate-type = ["cdylib", "lib"]
 
 [dependencies]
 rand = "0.8.3"
@@ -30,7 +30,7 @@ netlib-src = {version  = "0.8", optional = true}
 intel-mkl-src = {version = "0.6.0", optional = true}
 rcpr = { git = "https://github.com/drobnyjt/rcpr", optional = true}
 ndarray = {version = "0.14.0", features = ["serde"], optional = true}
-parry3d-f64 = {version = "0.2.0", optional = true}
+parry3d-f64 = {optional = true, version="0.2.0"}
 egui = {version = "0.15.0", optional = true}
 
 [dependencies.pyo3]

README.md

@@ -10,7 +10,7 @@ and Fortran.
 By discretizing the collision cascade into a sequence of binary collisions,
 [BCA] codes can accurately and efficiently model the prompt interaction
 between an energetic ion and a target material. This includes reflection,
-implantation, and transmission of the incident ion, s well as sputtering
+implantation, and transmission of the incident ion, as well as sputtering
 and displacement damage of the target. Generally, [BCA] codes can be
 valid for incident ion energies between approximately ~1 eV/nucleon 
 to <1 GeV/nucleon.
@@ -23,16 +23,17 @@ Journal of Open Source Software by clicking the badge below:
 
 ## Getting started
 
-The easiest way to get started is with the ergonomic Python functions currently on the development branch. Follow these steps to install, build, and run simple RustBCA simulations for sputtering yields and reflection coefficients:
+The easiest way to get started is with the ergonomic Python functions.
+Follow these steps to install, build, and run simple RustBCA simulations
+for sputtering yields and reflection coefficients:
 ```
 git clone https://github.com/lcpp-org/rustbca
 cd rustbca
-git checkout dev
 python -m pip install .
 python
 Python 3.9.6 (tags/v3.9.6:db3ff76, Jun 28 2021, 15:26:21) [MSC v.1929 64 bit (AMD64)] on win32
 Type "help", "copyright", "credits" or "license" for more information.
->>> from libRustBCA.pybca import *; from scripts.materials import *
+>>> from libRustBCA import *; from scripts.materials import *
 >>> angle = 0.0 # deg
 >>> energy = 1000.0 # eV
 >>> num_samples = 10000
@@ -44,7 +45,8 @@ Type "help", "copyright", "credits" or "license" for more information.
 ```
 
 For those eager to get started with the standalone code, try running one of the examples in the
-`RustBCA/examples` directory. Note that to automatically manipulate input files and reproduce the plots located on the [Wiki], these require several optional, but common,
+`RustBCA/examples` directory. Note that to automatically manipulate input files and reproduce
+the plots located on the [Wiki], these require several optional, but common,
 [Python] packages (`matplotlib`, `numpy`, `scipy`, `shapely`, and `toml`).
 
 ### H trajectories and collision cascades in a boron nitride dust grain
@@ -94,17 +96,17 @@ plt.show()
 The following features are implemented in `rustBCA`:
 
 * Ion-material interactions for all combinations of incident ion and target species.
-* Infinite, homogeneous targets (Mesh0D), Layered, finite-depth inhomogeneous targets (Mesh1D), arbitrary 2D geometry composition through a triangular mesh (Mesh2D), homogeneous spherical geometry (Sphere) and homogeneous, arbitrary triangular mesh geometry (TriMesh).
+* Infinite, homogeneous targets (Mesh0D), Layered, finite-depth inhomogeneous targets (Mesh1D), arbitrary 2D composition through a triangular mesh (Mesh2D), homogeneous spherical geometry (Sphere) and homogeneous 3D triangular mesh geometry (TriMesh).
 * Amorphous Solid/Liquid targets, Gaseous targets, and targets with both solid/liquid and gaseous elements
 * Low energy (< 25 keV/nucleon) electronic stopping modes including:
   * local (Oen-Robinson),
   * nonlocal (Lindhard-Scharff),
-  * and equipartition forms.
+  * and equipartition
 * Biersack-Varelas interpolation is also included for electronic stopping up to ~1 GeV/nucleon. Note that high energy physics beyond electronic stopping are not included.
 * Optionally, the Biersack-Haggmark treatment of high-energy free-flight paths between collisions can be included to greatly speed up high-energy simulations (i.e., by neglecting very small angle scattering).
 * A wide range of interaction potentials are provided, including:
   * the Kr-C, ZBL, Lenz-Jensen, and Moliere universal, screened-Coulomb potentials.
-  * the Lennard-Jones 12-6, Lennard-Jones 6.5-6, and Morse attractive-repulsive potentials.
+  * the Lennard-Jones 12-6 and Morse attractive-repulsive potentials.
 * Solving the distance-of-closest-approach problem is achieved using:
   * the Newton-Raphson method for simple root-finding,
   * or, for attractive-repulsive potentials, an Adaptive Chebyshev Proxy Rootfinder with Automatic Subdivision algorithm and a Polynomial root-finding algorithm are provided through the [rcpr] crate.
@@ -249,7 +251,7 @@ automatically during the build.
 
 ## Usage
 
-To use `RustBCA`, modify the `input.toml` file, which is used to configure each
+To use `RustBCA`, modify an `input.toml` file, which is used to configure each
 simulation.
 To run a simulation, execute:
 
@@ -279,4 +281,4 @@ Also have a look at the examples on the [Wiki] for writing `.toml` input files.
 [rustup]: https://rustup.rs
 [Rust]: https://en.wikipedia.org/wiki/Rust_(programming_language)
 [TOML]: https://en.wikipedia.org/wiki/TOML
-[Wiki]: https://github.com/lcpp-org/RustBCA/wiki
+[Wiki]: https://github.com/lcpp-org/RustBCA/wiki

examples/benchmark_eam.py

@@ -0,0 +1,60 @@
+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')
+from materials import *
+from formulas import *
+
+#Data digitized from Fig. 1 of https://doi.org/10.1016/0022-3115(84)90433-1
+#Modeled reflection coefficients using EAM at 100 eV and below and experimentally at 100 eV and above
+#Used to show that TRIM breaks down at higher incident energies
+data = np.array(
+    [[0.2962771, 0.18217821],
+    [0.96344626, 0.4871287],
+    [3.0372448, 0.8930693],
+    [9.876638, 0.86534655],
+    [30.184526, 0.7841584],
+    [96.644066, 0.6217822], #note: this is the last EAM data point
+    [102.83226, 0.4950495], #note: this is the first experimental data point
+    [203.52855, 0.37623763],
+    [516.34265, 0.3980198],
+    [809.75903, 0.32277226],
+    [1006.2257, 0.2990099],
+    [2054.3218, 0.17821783],
+    [4129.5522, 0.13069306],
+    [6890.884, 0.0990099]]
+)
+
+#Setting the cutoff low is necessary since the incident energies are also low.
+hydrogen['Ec'] = 0.001
+
+#Plotting the EAM data points.
+energies = data[:6, 0]
+r_benchmark = data[:6, 1]
+plt.semilogx(energies, r_benchmark, marker='o', linestyle='', label='EAM')
+
+#Plotting the experimental data points.
+energies = data[6:, 0]
+r_benchmark = data[6:, 1]
+plt.semilogx(energies, r_benchmark, marker='^', linestyle='', label='Exp.')
+
+#Plotting RustBCA data points, using the ergonomic helper function reflection_coefficient().
+energies = np.logspace(-1, 4, 25)
+r_rustbca = [reflection_coefficient(hydrogen, nickel, energy, 0.0, 1000)[0] for energy in energies]
+plt.semilogx(energies, r_rustbca, label='RustBCA', color='black')
+
+r_thomas = [thomas_reflection(hydrogen, nickel, energy) for energy in energies]
+plt.semilogx(energies, r_thomas, label='Thomas', color='red', linestyle='--')
+
+plt.title('Reflection Coefficients of Normal Incidence H on Ni')
+plt.xlabel('E [eV]')
+plt.ylabel('R')
+plt.gca().set_ylim([0.0, 1.0])
+plt.legend(loc='upper right')
+
+plt.show()
+

examples/benchmark_nnp.py

@@ -0,0 +1,132 @@
+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')
+from materials import *
+from formulas import *
+
+R_T_Be = np.array(
+    [
+    [9.67742,	0.016754618],
+    [20.32258,	0.09511873],
+    [30.0,	0.13153034],
+    [49.677418,	0.102242745],
+    [74.83871,	0.08403694],
+    [100.32258,	0.031002639],
+    ]
+)
+
+R_D_Be = np.array(
+    [
+    [10.32258,	0.11649077],
+    [20.32258,	0.19168866],
+    [30.32258,	0.19722955],
+    [49.677418,	0.17269129],
+    [75.16129,	0.11174142],
+    [100.0,	0.044459105],
+    ]
+)
+
+R_H_Be = np.array(
+    [
+    [9.67742,	0.19485489],
+    [19.67742,	0.32150397],
+    [30.0,	0.33575198],
+    [50.322582,	0.25105542],
+    [75.16129,	0.1378628],
+    [100.645164,	0.044459105],
+    ]
+)
+
+#Plotting the MD data points.
+energies = R_H_Be[:, 0]
+R = R_H_Be[:, 1]
+plt.plot(energies, R, marker='o', linestyle='', label='H on Be MD+NNP')
+
+#Plotting the MD data points.
+energies = R_D_Be[:, 0]
+R = R_D_Be[:, 1]
+plt.plot(energies, R, marker='x', linestyle='', label='D on Be MD+NNP')
+
+#Plotting the MD data points.
+energies = R_T_Be[:, 0]
+R = R_T_Be[:, 1]
+plt.plot(energies, R, marker='^', linestyle='', label='T on Be MD+NNP')
+
+
+for ion in [hydrogen, deuterium, tritium]:
+    energies = np.logspace(1, 2, 25)
+    r_rustbca = [reflection_coefficient(ion, beryllium, energy, 0.0, 10000)[0] for energy in energies]
+    r_thomas = [thomas_reflection(ion, beryllium, energy) for energy in energies]
+    line = plt.plot(energies, r_rustbca, label=f'{ion["symbol"]} on Be RustBCA')[0]
+    plt.plot(energies, r_thomas, label=f'{ion["symbol"]} on Be Thomas', linestyle='--', color=line.get_color())
+
+plt.xlabel('E [eV]')
+plt.ylabel('R_N')
+plt.title('Reflection Coefficients')
+plt.legend()
+plt.savefig('reflection_md_nnp_rustbca.png')
+
+#Data digitized from https://doi.org/10.1088/1741-4326/ac592a
+#MD+NNP Sputtering Yields
+y_H_Be = np.array([
+    [20.004726, 0.0029361406],
+    [29.880848, 0.013341782],
+    [49.791046, 0.024575423],
+    [75.00529,	0.016207783],
+    [99.68844,	0.017530797],
+])
+
+y_D_Be = np.array([
+    [19.9699, 0.005002439],
+    [29.601622, 0.021064475],
+    [49.76866,	0.03461476],
+    [74.71549,	0.030082114],
+    [99.954605, 0.013559438],
+])
+
+y_T_Be = np.array([
+    [20.013433, 0.0025699555],
+    [29.85846,	0.01879205],
+    [49.756844, 0.04147379],
+    [74.707405, 0.034042966],
+    [99.68098, 0.01965117],
+])
+
+
+plt.figure()
+
+#Plotting the MD data points.
+energies = y_H_Be[:, 0]
+y = y_H_Be[:, 1]
+plt.semilogy(energies, y, marker='o', linestyle='', label='H on Be MD+NNP')
+
+#Plotting the MD data points.
+energies = y_D_Be[:, 0]
+y = y_D_Be[:, 1]
+plt.semilogy(energies, y, marker='x', linestyle='', label='D on Be MD+NNP')
+
+#Plotting the MD data points.
+energies = y_T_Be[:, 0]
+y = y_T_Be[:, 1]
+plt.semilogy(energies, y, marker='^', linestyle='', label='T on Be MD+NNP')
+
+for ion in [hydrogen, deuterium, tritium]:
+    energies = np.logspace(1, 2, 25)
+    y_rustbca = [sputtering_yield(ion, beryllium, energy, 0.0, 100000) for energy in energies]
+    plt.semilogy(energies, y_rustbca, label=f'{ion["symbol"]} on Be RustBCA')
+
+yamamura_yield = [yamamura(hydrogen, beryllium, energy) for energy in energies]
+plt.semilogy(energies, yamamura_yield, label='Yamamura H on Be', linestyle='--')
+
+plt.title('Sputtering Yields of Hydrogenic Species on Be')
+plt.xlabel('E [eV]')
+plt.ylabel('Y [at/ion]')
+plt.gca().set_ylim([1e-4, 1e-1])
+plt.legend(loc='lower right')
+
+plt.savefig('sputtering_md_nnp_rustbca.png')