diff --git a/.github/workflows/rustbca_compile_check.yml b/.github/workflows/rustbca_compile_check.yml index 0ab499f..584476c 100644 --- a/.github/workflows/rustbca_compile_check.yml +++ b/.github/workflows/rustbca_compile_check.yml @@ -23,7 +23,7 @@ jobs: sudo apt-get install curl - name: Install rust run: | - curl --proto '=https' --tlsv1.2 -sSf -y https://sh.rustup.rs | sh + curl --proto '=https' --tlsv1.2 -sSf -y https://sh.rustup.rs | sh sudo apt-get install rustc cargo - name: Install pip for Python-3 run: | @@ -43,7 +43,7 @@ jobs: python3 examples/test_rustbca.py - name: Test Fortran and C bindings run : | - cargo build --release + cargo build --release --lib --features parry3d cp examples/test_rustbca.f90 . gfortran -c rustbca.f90 target/release/liblibRustBCA.so gfortran test_rustbca.f90 rustbca.f90 target/release/liblibRustBCA.so @@ -53,7 +53,7 @@ jobs: ./a.out - name: Test RustBCA run: | - cargo test --features cpr_rootfinder_netlib,hdf5_input,distributions,parry3d + cargo test --features cpr_rootfinder,hdf5_input,distributions,parry3d - name: Run Examples run: | cargo run --release 0D examples/boron_nitride_0D.toml @@ -65,4 +65,6 @@ jobs: ./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 - + cat boron_nitride_summary.output + ./target/release/RustBCA HOMOGENEOUS2D examples/boron_nitride_wire_homogeneous.toml + cat boron_nitride_h_summary.output diff --git a/Cargo.toml b/Cargo.toml index 6ec7a91..f4738f1 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,9 +1,9 @@ [package] name = "RustBCA" -version = "2.0.1" +version = "2.8.4" default-run = "RustBCA" authors = ["Jon Drobny ", "Jon Drobny "] -edition = "2018" +edition = "2021" [[bin]] name = "RustBCA" @@ -15,26 +15,22 @@ path = "src/lib.rs" crate-type = ["cdylib", "lib"] [dependencies] -rand = "0.8.3" -rand_distr = "0.4.2" -toml = "0.5.8" -anyhow = "1.0.38" -itertools = "0.10.0" -rayon = "1.5.0" -geo = {version = "0.17.1", optional = false} +rand = "0.8.5" +rand_distr = "0.4.3" +toml = "0.7.4" +anyhow = "1.0.71" +itertools = "0.10.5" +rayon = "1.7.0" +geo = {version = "0.25", optional = false} indicatif = {version = "0.15.0", features=["rayon"]} -serde = { version = "1.0.123", features = ["derive"] } +serde = { version = "1.0.163", features = ["derive"] } hdf5 = {version = "0.7.1", optional = true} -openblas-src = {version = "0.9", optional = true} -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 = {optional = true, version="0.2.0"} -egui = {version = "0.15.0", optional = true} [dependencies.pyo3] -version = "0.13.2" +version = "0.19.0" features = ["extension-module"] optional = true @@ -49,9 +45,7 @@ debug = false [features] hdf5_input = ["hdf5"] -cpr_rootfinder_openblas = ["rcpr", "openblas-src"] -cpr_rootfinder_netlib = ["rcpr", "netlib-src"] -cpr_rootfinder_intel_mkl = ["rcpr", "intel-mkl-src"] +cpr_rootfinder = ["rcpr"] distributions = ["ndarray"] no_list_output = [] parry3d = ["parry3d-f64"] diff --git a/README.md b/README.md index 98ed2be..c9f2e2c 100644 --- a/README.md +++ b/README.md @@ -3,7 +3,7 @@ `RustBCA` is a general-purpose, high-performance code for simulating ion-material interactions including sputtering, reflection, and implantation using the binary collision approximation ([BCA]), written in [Rust]! -RustBCA consists of a standalone code and libraries for including +RustBCA includes a standalone version and libraries for including ion-material interactions in simulations written in C/C++, Python, and Fortran. @@ -13,7 +13,9 @@ between an energetic ion and a target material. This includes reflection, 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. +to <1 GeV/nucleon. Improvements to RustBCA have expanded the regime +of validity for some quantities, such as reflection coefficients, below +1 eV/nucleon. Check out the `RustBCA` [Wiki] for detailed information, installation instructions, use cases, examples, and more. See the RustBCA paper at the @@ -27,6 +29,7 @@ Selected citations of RustBCA as of 5/24/23: * [hPIC2: A hardware-accelerated, hybrid particle-in-cell code for dynamic plasma-material interactions](https://doi.org/10.1016/j.cpc.2022.108569), LT Meredith et al. (2023) * [Global sensitivity analysis of a coupled multiphysics model to predict surface evolution in fusion plasma–surface interactions](https://doi.org/10.1016/j.commatsci.2023.112229), P. Robbe et al. (2023) * [Modeling the effect of nitrogen recycling on the erosion and leakage of tungsten impurities from the SAS-VW divertor in DIII-D during nitrogen gas injection](https://doi.org/10.1016/j.nme.2022.101254), MS Parsons et al. (2023) +* [Enabling attractive-repulsive potentials in binary-collision-approximation monte-carlo codes for ion-surface interactions](https://doi.org/10.1088/2053-1591/ad1262), J Drobny and D Curreli (2023) ## Getting started @@ -53,7 +56,7 @@ 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, +the plots located on the [Wiki], these may require some optional [Python] packages (`matplotlib`, `numpy`, `scipy`, `shapely`, and `toml`). ### H trajectories and collision cascades in a boron nitride dust grain @@ -103,32 +106,33 @@ 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 composition through a triangular mesh (Mesh2D), homogeneous spherical geometry (Sphere) and homogeneous 3D triangular mesh geometry (TriMesh). +* Infinite, homogeneous targets (Mesh0D), Layered, finite-depth inhomogeneous targets (Mesh1D), arbitrary 2D composition through a triangular mesh (Mesh2D), fast homogeneous 2D geometry (Homogeneous2D), 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 -* 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). +* 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, and that Biersack-Varelas may not be as accurate as other methods. +* 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 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. + * the Newton-Raphson method for purely repulsive potentials, + * or, for attractive-repulsive potentials, an Adaptive Chebyshev Proxy Rootfinder with Automatic Subdivision algorithm and a polynomial root-finding algorithm are provided through [rcpr]. * Multiple interaction potentials can be used in a single simulation for any number of potentials/species. * For example, the He-W interaction can be specified using a Lennard-Jones 12-6 potential, while the W-W interaction can be defined using a Kr-C potential. * The scattering integral can be calculated using: * Gauss-Mehler quadrature, * Gauss-Legendre quadrature, * Mendenall-Weller quadrature, - * or the MAGIC algorithm. + * or the MAGIC algorithm (for certain screened Coulomb potentials only). * Input files use the [TOML] format, making them both human-readable and easily parsable. * RustBCA generates user-friendly, context-providing error messages, which help pinpoint the cause of errors and provide suggested fixes to the user. * The simulation results are comma-delimited (`csv` format) and include: * the energies and directions of emitted particles (reflected ions and sputtered atoms), * the final positions of implanted ions, + * displacements, * full trajectory tracking for both the incident ions and target atoms, * and many other parameters such as position of origin of sputtered particles and energy loss along trajectories. * Optionally, the code can produce energy-angle and implantation distributions when built with the `--features distributions` flag and disable space-intensive particle list output with `--features no_list_output`. @@ -140,11 +144,6 @@ Without optional features, `RustBCA` should compile with `cargo` alone on Windows, MacOS, and Linux systems. [HDF5] for particle list input has been tested on Windows, but version 1.10.6 must be used. -[rcpr], the adaptive Chebyshev Proxy Rootfinder with automatic subdivision and -polynomial rootfinder package for [Rust], has not yet been successfully compiled -on Windows. -However, it can be compiled on the Windows Subsystem for Linux (WSL) and, likely, -on Ubuntu for Windows or Cygwin. #### Manual Dependences @@ -157,7 +156,7 @@ on Ubuntu for Windows or Cygwin. #### Optional Dependencies * [HDF5] libraries -* [rcpr], a CPR and polynomial rootfinder, required for using attractive-repulsive interaction potentials such as Lennard-Jones or Morse. It may require additional software (see below). +* [rcpr], a CPR and polynomial rootfinder, required for using attractive-repulsive interaction potentials such as Lennard-Jones or Morse. * For manipulating input files and running associated scripts, the following are required: * [Python] 3.6+ * The [Python] libraries: `numpy`, `matplotlib`, `toml` (must build from source), `shapely`, and `scipy`. @@ -191,17 +190,17 @@ python3 setup.py install ```bash sudo apt-get install gcc gfortran build-essential cmake liblapack-dev libblas-dev liblapacke-dev ``` -8. Install `cargo`: +8. (Optional - should come with rustup) Install `cargo`: ```bash sudo apt-get install cargo ``` 9. Build `RustBCA`: ```bash git clone https://github.com/lcpp-org/rustBCA -cd rustBCA +cd RustBCA cargo build --release ``` -10. (Optional) Build `rustBCA` with optional dependencies, `hdf5` and/or `rcpr` (with your choice of backend: `openblas`, `netlib`, or `intel-mkl`): +10. (Optional) Build `RustBCA` with optional dependencies, `hdf5` and/or `rcpr` (with your choice of backend: `openblas`, `netlib`, or `intel-mkl`): ```bash cargo build --release --features cpr_rootfinder_netlib,hdf5_input cargo build --release --features cpr_rootfinder_openblas,hdf5_input @@ -221,7 +220,7 @@ cargo test --features cpr_rootfinder_intel_mkl ### Detailed instructions for Fedora 33 -Most of the ingredients for building `rustBCA` and running the [Python] helper +Most of the ingredients for building `RustBCA` and running the [Python] helper scripts are available natively in the Fedora software repository, so the setup is relatively painless. @@ -282,7 +281,7 @@ Additionally, `RustBCA` accepts an input file type (one of: `0D`, `1D`, `2D`, `T For further details, have a look at [Usage](https://github.com/lcpp-org/RustBCA/wiki/Usage,-Input-File,-and-Output-Files) -on the `rustBCA` [Wiki] for usage instructions. +on the `RustBCA` [Wiki] for usage instructions. Also have a look at the examples on the [Wiki] for writing `.toml` input files. [BCA]: https://en.wikipedia.org/wiki/Binary_collision_approximation diff --git a/RustBCA.h b/RustBCA.h index c1219a4..8d60a7b 100644 --- a/RustBCA.h +++ b/RustBCA.h @@ -1,155 +1,203 @@ -#include -#include -#include -#include -#include -#include - -static const double PI = M_PI; - -static const double FRAC_2_SQRT_PI = M_2_PI; - -static const double SQRT_2 = M_SQRT2; - -///Fundamental charge in Coulombs. -static const double Q = 1.602176634e-19; - -/// One electron-volt in Joules. -static const double EV = Q; - -/// One atomic mass unit in kilograms. -static const double AMU = 1.66053906660e-27; - -/// One Angstrom in meters. -static const double ANGSTROM = 1e-10; - -/// One micron in meters. -static const double MICRON = 1e-6; - -/// One nanometer in meters. -static const double NM = 1e-9; - -/// One centimeter in meters. -static const double CM = 1e-2; - -/// Vacuum permitivity in Farads/meter. -static const double EPS0 = 8.8541878128e-12; - -/// Bohr radius in meters. -static const double A0 = 5.29177210903e-11; - -/// Electron mass in kilograms. -static const double ME = 9.1093837015e-31; - -/// sqrt(pi). -static const double SQRTPI = (2. / FRAC_2_SQRT_PI); - -/// sqrt(2 * pi). -static const double SQRT2PI = ((2. * SQRT_2) / FRAC_2_SQRT_PI); - -/// Speed of light in meters/second. -static const double C = 299792458.; - -/// Bethe-Bloch electronic stopping prefactor, in SI units. -static const double BETHE_BLOCH_PREFACTOR = ((((((4. * PI) * ((Q * Q) / ((4. * PI) * EPS0))) * ((Q * Q) / ((4. * PI) * EPS0))) / ME) / C) / C); - -/// Lindhard-Scharff electronic stopping prefactor, in SI units. -static const double LINDHARD_SCHARFF_PREFACTOR = (((1.212 * ANGSTROM) * ANGSTROM) * Q); - -/// Lindhard reduced energy prefactor, in SI units. -static const double LINDHARD_REDUCED_ENERGY_PREFACTOR = ((((4. * PI) * EPS0) / Q) / Q); - -struct OutputBCA { - uintptr_t len; - double (*particles)[9]; -}; - -struct InputSimpleBCA { - uintptr_t len; - /// vx, vy, vz - double (*velocities)[3]; - double Z1; - double m1; - double Ec1; - double Es1; - double Z2; - double m2; - double n2; - double Ec2; - double Es2; - double Eb2; -}; - -struct InputCompoundBCA { - uintptr_t len; - /// vx, vy, vz - double (*velocities)[3]; - double Z1; - double m1; - double Ec1; - double Es1; - uintptr_t num_species_target; - double *Z2; - double *m2; - double *n2; - double *Ec2; - double *Es2; - double *Eb2; -}; - -struct OutputTaggedBCA { - uintptr_t len; - double (*particles)[9]; - double *weights; - int32_t *tags; - bool *incident; -}; - -struct InputTaggedBCA { - uintptr_t len; - /// x y z - double (*positions)[3]; - /// vx, vy, vz - double (*velocities)[3]; - double Z1; - double m1; - double Ec1; - double Es1; - uintptr_t num_species_target; - double *Z2; - double *m2; - double *n2; - double *Ec2; - double *Es2; - double *Eb2; - int32_t *tags; - double *weights; -}; - -extern "C" { - -OutputTaggedBCA compound_tagged_bca_list_c(InputTaggedBCA input); - -OutputBCA simple_bca_list_c(InputSimpleBCA input); - -OutputBCA compound_bca_list_c(InputCompoundBCA input); - -OutputBCA simple_bca_c(double x, - double y, - double z, - double ux, - double uy, - double uz, - double E1, - double Z1, - double m1, - double Ec1, - double Es1, - double Z2, - double m2, - double Ec2, - double Es2, - double n2, - double Eb2); - -} // extern "C" +#include +#include +#include +#include +#include +#include + +static const double PI = M_PI; + +static const double FRAC_2_SQRT_PI = M_2_PI; + +static const double SQRT_2 = M_SQRT2; + +///Fundamental charge in Coulombs. +static const double Q = 1.602176634e-19; + +/// One atomic mass unit in kilograms. +static const double AMU = 1.66053906660e-27; + +/// One Angstrom in meters. +static const double ANGSTROM = 1e-10; + +/// One micron in meters. +static const double MICRON = 1e-6; + +/// One nanometer in meters. +static const double NM = 1e-9; + +/// One centimeter in meters. +static const double CM = 1e-2; + +/// One milimter in meters. +static const double MM = 1e-3; + +/// Vacuum permitivity in Farads/meter. +static const double EPS0 = 8.8541878128e-12; + +/// Bohr radius in meters. +static const double A0 = 5.29177210903e-11; + +/// Electron mass in kilograms. +static const double ME = 9.1093837015e-31; + +/// sqrt(pi). +static const double SQRTPI = (2. / FRAC_2_SQRT_PI); + +/// sqrt(2 * pi). +static const double SQRT2PI = ((2. * SQRT_2) / FRAC_2_SQRT_PI); + +/// Speed of light in meters/second. +static const double C = 299792458.; + +/// Bethe-Bloch electronic stopping prefactor, in SI units. +static const double BETHE_BLOCH_PREFACTOR = ((((((4. * PI) * ((Q * Q) / ((4. * PI) * EPS0))) * ((Q * Q) / ((4. * PI) * EPS0))) / ME) / C) / C); + +/// Lindhard-Scharff electronic stopping prefactor, in SI units. +static const double LINDHARD_SCHARFF_PREFACTOR = (((1.212 * ANGSTROM) * ANGSTROM) * Q); + +/// Lindhard reduced energy prefactor, in SI units. +static const double LINDHARD_REDUCED_ENERGY_PREFACTOR = ((((4. * PI) * EPS0) / Q) / Q); + +struct OutputTaggedBCA { + uintptr_t len; + double (*particles)[9]; + double *weights; + int32_t *tags; + bool *incident; +}; + +struct InputTaggedBCA { + uintptr_t len; + /// x y z + double (*positions)[3]; + /// vx, vy, vz + double (*velocities)[3]; + double Z1; + double m1; + double Ec1; + double Es1; + uintptr_t num_species_target; + double *Z2; + double *m2; + double *n2; + double *Ec2; + double *Es2; + double *Eb2; + int32_t *tags; + double *weights; +}; + +struct OutputBCA { + uintptr_t len; + double (*particles)[9]; +}; + +struct InputSimpleBCA { + uintptr_t len; + /// vx, vy, vz + double (*velocities)[3]; + double Z1; + double m1; + double Ec1; + double Es1; + double Z2; + double m2; + double n2; + double Ec2; + double Es2; + double Eb2; +}; + +struct InputCompoundBCA { + uintptr_t len; + /// vx, vy, vz + double (*velocities)[3]; + double Z1; + double m1; + double Ec1; + double Es1; + uintptr_t num_species_target; + double *Z2; + double *m2; + double *n2; + double *Ec2; + double *Es2; + double *Eb2; +}; + +extern "C" { + +OutputTaggedBCA compound_tagged_bca_list_c(InputTaggedBCA input); + +void drop_output_tagged_bca(OutputTaggedBCA output); + +void drop_output_bca(OutputBCA output); + +void reflect_single_ion_c(int *num_species_target, + double *ux, + double *uy, + double *uz, + double *E1, + double *Z1, + double *m1, + double *Ec1, + double *Es1, + double *Z2, + double *m2, + double *Ec2, + double *Es2, + double *Eb2, + double *n2); + +OutputBCA simple_bca_list_c(InputSimpleBCA input); + +OutputBCA compound_bca_list_c(InputCompoundBCA input); + +const double (*compound_bca_list_fortran(int *num_incident_ions, + bool *track_recoils, + double *ux, + double *uy, + double *uz, + double *E1, + double *Z1, + double *m1, + double *Ec1, + double *Es1, + int *num_species_target, + double *Z2, + double *m2, + double *Ec2, + double *Es2, + double *Eb2, + double *n2, + int *num_emitted_particles))[6]; + +OutputBCA simple_bca_c(double x, + double y, + double z, + double ux, + double uy, + double uz, + double E1, + double Z1, + double m1, + double Ec1, + double Es1, + double Z2, + double m2, + double Ec2, + double Es2, + double n2, + double Eb2); + +void rotate_given_surface_normal(double nx, + double ny, + double nz, + double *ux, + double *uy, + double *uz); + +void rotate_back(double nx, double ny, double nz, double *ux, double *uy, double *uz); + +} // extern "C" diff --git a/examples/RustBCA.c b/examples/RustBCA.c index a218531..9f88ef9 100644 --- a/examples/RustBCA.c +++ b/examples/RustBCA.c @@ -50,5 +50,44 @@ int main(int argc, char * argv[]) { std::cout << "Particle 2 E [eV]: "; std::cout << output.particles[1][2]; std::cout << std::endl; + + drop_output_tagged_bca(output); + + double nx = -0.707106; + double ny = -0.707106; + double nz = 0.0; + + double ux = 1.0; + double uy = 0.0; + double uz = 0.0; + + std::cout << "Before rotation: "; + std::cout << ux; + std::cout << ", "; + std::cout << uy; + std::cout << ", "; + std::cout << uz; + std::cout << std::endl; + + rotate_given_surface_normal(nx, ny, nz, &ux, &uy, &uz); + + std::cout << "After rotation: "; + std::cout << ux; + std::cout << ", "; + std::cout << uy; + std::cout << ", "; + std::cout << uz; + std::cout << std::endl; + + rotate_back(nx, ny, nz, &ux, &uy, &uz); + + std::cout << "After rotation back: "; + std::cout << ux; + std::cout << ", "; + std::cout << uy; + std::cout << ", "; + std::cout << uz; + std::cout << std::endl; + return 0; } diff --git a/examples/boron_nitride_0D.toml b/examples/boron_nitride_0D.toml index b6723c4..f6f64f1 100644 --- a/examples/boron_nitride_0D.toml +++ b/examples/boron_nitride_0D.toml @@ -2,12 +2,14 @@ name = "boron_nitride_" weak_collision_order = 0 num_threads = 4 +track_displacements = true [material_parameters] energy_unit = "EV" mass_unit = "AMU" Eb = [ 0.0, 0.0,] Es = [ 5.76, 0.0,] +Ed = [25.0, 25.0] Ec = [ 1.0, 1.0,] Z = [ 5, 7,] m = [ 10.811, 14,] diff --git a/examples/boron_nitride_wire.toml b/examples/boron_nitride_wire.toml index a0007b9..8f7aa84 100644 --- a/examples/boron_nitride_wire.toml +++ b/examples/boron_nitride_wire.toml @@ -13,13 +13,13 @@ Z = [ 5, 7,] m = [ 10.811, 14,] [particle_parameters] -length_unit = "MICRON" +length_unit = "ANGSTROM" energy_unit = "EV" mass_unit = "AMU" N = [ 10000 ] m = [ 1.008 ] Z = [ 1 ] -E = [ 1000.0 ] +E = [ 5000.0 ] Ec = [ 1.0 ] Es = [ 10.0 ] pos = [ [ 0.0, 0.0, 0.0,] ] diff --git a/examples/boron_nitride_wire_homogeneous.toml b/examples/boron_nitride_wire_homogeneous.toml new file mode 100644 index 0000000..c7f40a1 --- /dev/null +++ b/examples/boron_nitride_wire_homogeneous.toml @@ -0,0 +1,33 @@ +[options] +name = "boron_nitride_h_" +weak_collision_order = 0 +num_threads = 4 + +[material_parameters] +energy_unit = "EV" +mass_unit = "AMU" +Eb = [ 0.0, 0.0,] +Es = [ 5.76, 0.0,] +Ec = [ 1.0, 1.0,] +Z = [ 5, 7,] +m = [ 10.811, 14,] + +[particle_parameters] +length_unit = "ANGSTROM" +energy_unit = "EV" +mass_unit = "AMU" +N = [ 10000 ] +m = [ 1.008 ] +Z = [ 1 ] +E = [ 5000.0 ] +Ec = [ 1.0 ] +Es = [ 10.0 ] +pos = [ [ 0.0, 0.0, 0.0,] ] +dir = [ [ 0.9999999999984769, 1.7453292519934434e-6, 0.0,] ] + +[geometry_input] +length_unit = "ANGSTROM" +points = [ [ 0.0, 0.0,], [ 1000.0, 0.0,], [ 998.0267284282716, 62.79051952931337,], [ 992.1147013144779, 125.33323356430427,], [ 982.2872507286887, 187.38131458572462,], [ 968.5831611286311, 248.68988716485478,], [ 951.0565162951535, 309.0169943749474,], [ 929.7764858882514, 368.12455268467795,], [ 904.8270524660195, 425.7792915650727,], [ 876.3066800438636, 481.7536741017153,], [ 844.3279255020151, 535.8267949789966,], [ 809.0169943749474, 587.7852522924732,], [ 770.5132427757892, 637.4239897486898,], [ 728.9686274214115, 684.5471059286888,], [ 684.5471059286887, 728.9686274214115,], [ 637.4239897486897, 770.5132427757893,], [ 587.785252292473, 809.0169943749474,], [ 535.8267949789965, 844.3279255020151,], [ 481.75367410171515, 876.3066800438637,], [ 425.77929156507264, 904.8270524660196,], [ 368.12455268467784, 929.7764858882515,], [ 309.01699437494744, 951.0565162951535,], [ 248.68988716485475, 968.5831611286311,], [ 187.3813145857245, 982.2872507286887,], [ 125.33323356430427, 992.1147013144779,], [ 62.79051952931331, 998.0267284282716,], [ -1.6081226496766366e-13, 1000.0,], [ -62.7905195293134, 998.0267284282716,], [ -125.33323356430436, 992.1147013144778,], [ -187.38131458572482, 982.2872507286886,], [ -248.68988716485484, 968.5831611286311,], [ -309.01699437494756, 951.0565162951535,], [ -368.12455268467795, 929.7764858882515,], [ -425.7792915650727, 904.8270524660195,], [ -481.75367410171543, 876.3066800438635,], [ -535.8267949789969, 844.327925502015,], [ -587.785252292473, 809.0169943749474,], [ -637.4239897486898, 770.5132427757893,], [ -684.5471059286888, 728.9686274214115,], [ -728.9686274214116, 684.5471059286886,], [ -770.5132427757893, 637.4239897486896,], [ -809.0169943749473, 587.7852522924733,], [ -844.3279255020151, 535.8267949789966,], [ -876.3066800438636, 481.7536741017152,], [ -904.8270524660196, 425.77929156507247,], [ -929.7764858882515, 368.1245526846778,], [ -951.0565162951536, 309.0169943749471,], [ -968.5831611286311, 248.6898871648548,], [ -982.2872507286887, 187.38131458572457,], [ -992.1147013144779, 125.3332335643041,], [ -998.0267284282716, 62.79051952931314,], [ -1000.0, -3.216245299353273e-13,], [ -998.0267284282716, -62.79051952931334,], [ -992.1147013144779, -125.33323356430428,], [ -982.2872507286886, -187.38131458572477,], [ -968.5831611286311, -248.689887164855,], [ -951.0565162951535, -309.0169943749477,], [ -929.7764858882513, -368.1245526846783,], [ -904.8270524660195, -425.77929156507264,], [ -876.3066800438636, -481.7536741017154,], [ -844.327925502015, -535.8267949789968,], [ -809.0169943749472, -587.7852522924734,], [ -770.5132427757891, -637.42398974869,], [ -728.9686274214115, -684.5471059286888,], [ -684.5471059286887, -728.9686274214115,], [ -637.4239897486896, -770.5132427757893,], [ -587.7852522924733, -809.0169943749473,], [ -535.8267949789963, -844.3279255020153,], [ -481.75367410171526, -876.3066800438636,], [ -425.7792915650722, -904.8270524660198,], [ -368.1245526846778, -929.7764858882515,], [ -309.01699437494756, -951.0565162951535,], [ -248.68988716485444, -968.5831611286312,], [ -187.38131458572462, -982.2872507286887,], [ -125.33323356430373, -992.1147013144779,], [ -62.79051952931321, -998.0267284282716,], [ -1.8369701987210297e-13, -1000.0,], [ 62.79051952931372, -998.0267284282716,], [ 125.33323356430424, -992.1147013144779,], [ 187.38131458572514, -982.2872507286886,], [ 248.68988716485492, -968.5831611286311,], [ 309.0169943749472, -951.0565162951536,], [ 368.12455268467824, -929.7764858882513,], [ 425.7792915650726, -904.8270524660196,], [ 481.7536741017157, -876.3066800438634,], [ 535.8267949789968, -844.327925502015,], [ 587.7852522924737, -809.016994374947,], [ 637.42398974869, -770.5132427757891,], [ 684.5471059286887, -728.9686274214115,], [ 728.9686274214118, -684.5471059286883,], [ 770.5132427757893, -637.4239897486897,], [ 809.0169943749478, -587.7852522924726,], [ 844.3279255020152, -535.8267949789963,], [ 876.3066800438636, -481.7536741017153,], [ 904.8270524660197, -425.77929156507224,], [ 929.7764858882515, -368.12455268467784,], [ 951.0565162951538, -309.01699437494676,], [ 968.5831611286312, -248.6898871648545,], [ 982.2872507286887, -187.38131458572468,], [ 992.1147013144779, -125.33323356430378,], [ 998.0267284282716, -62.790519529313265,], [ 1000.0, -2.4492935982947065e-13,],] +simulation_boundary_points = [ [ -1100.0, -1100.0,], [ -1100.0, 1100.0,], [ 1100.0, 1100.0,], [ 1100.0, -1100.0,],] +densities = [ 0.065, 0.065,] +electronic_stopping_correction_factor = 1.0 diff --git a/examples/test_morse.py b/examples/test_morse.py index 7aa025f..6de9119 100644 --- a/examples/test_morse.py +++ b/examples/test_morse.py @@ -10,8 +10,8 @@ from formulas import * #Corrections to hydrogen to enable low-E simualtions and to match literature Es for H on Ni -hydrogen['Ec'] = 0.001 -hydrogen['Es'] = 2.5 +hydrogen['Ec'] = 0.1 +hydrogen['Es'] = 1.5 #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): @@ -19,23 +19,15 @@ def run_morse_potential(energy, index, num_samples=10000, run_sim=True): input_file = f''' [options] name = "morse_{index}" - track_trajectories = false track_recoils = false - track_recoil_trajectories = false - write_buffer_size = 8000 weak_collision_order = 0 - suppress_deep_recoils = false - high_energy_free_flight_paths = false electronic_stopping_mode = "LOW_ENERGY_NONLOCAL" mean_free_path_model = "LIQUID" interaction_potential = [[{{"MORSE"={{D=5.4971E-20, r0=2.782E-10, alpha=1.4198E10}}}}]] - scattering_integral = [[{{"GAUSS_MEHLER"={{n_points = 100}}}}]] - root_finder = [[{{"CPR"={{n0=2, nmax=100, epsilon=1E-6, complex_threshold=1E-12, truncation_threshold=1E-12, far_from_zero=1E6, interval_limit=1E-16, derivative_free=false}}}}]] - use_hdf5 = false + scattering_integral = [["GAUSS_LEGENDRE"]] + root_finder = [[{{"CPR"={{n0=2, nmax=100, epsilon=1E-9, complex_threshold=1E-3, truncation_threshold=1E-9, far_from_zero=1E9, interval_limit=1E-12, derivative_free=true}}}}]] num_threads = 4 num_chunks = 10 - track_displacements = false - track_energy_losses = false [particle_parameters] length_unit = "ANGSTROM" @@ -45,8 +37,8 @@ def run_morse_potential(energy, index, num_samples=10000, run_sim=True): m = [ 1.008 ] Z = [ 1 ] E = [ {energy} ] - Ec = [ 0.001 ] - Es = [ 2.5 ] + Ec = [ 0.1 ] + Es = [ 1.5 ] interaction_index = [ 0 ] pos = [ [ -4.4, 0.0, 0.0,] ] dir = [ [ 0.9999999999984769, 1.7453292519934434e-6, 0.0,] ] @@ -72,7 +64,7 @@ def run_morse_potential(energy, index, num_samples=10000, run_sim=True): with open(f'morse_{index}.toml', 'w') as file: file.write(input_file) - if run_sim: os.system(f'cargo run --release --features cpr_rootfinder_netlib 0D morse_{index}.toml') + if run_sim: os.system(f'cargo run --release --features cpr_rootfinder 0D morse_{index}.toml') reflected_list = np.atleast_2d(np.genfromtxt(f'morse_{index}reflected.output', delimiter=',')) if np.size(reflected_list) > 0: @@ -84,10 +76,72 @@ def run_morse_potential(energy, index, num_samples=10000, run_sim=True): return num_reflected/num_samples, energy_reflected/energy/num_samples +def run_krc_morse_potential(energy, index, num_samples=10000, run_sim=True): + + input_file = f''' + [options] + name = "krc_morse_{index}" + track_recoils = false + weak_collision_order = 0 + electronic_stopping_mode = "LOW_ENERGY_NONLOCAL" + 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}}}}]] + scattering_integral = [["GAUSS_LEGENDRE"]] + root_finder = [[{{"CPR"={{n0=2, nmax=100, epsilon=1E-9, complex_threshold=1E-3, truncation_threshold=1E-9, far_from_zero=1E9, interval_limit=1E-12, derivative_free=true}}}}]] + num_threads = 4 + num_chunks = 10 + + [particle_parameters] + length_unit = "ANGSTROM" + energy_unit = "EV" + mass_unit = "AMU" + N = [ {num_samples} ] + m = [ 1.008 ] + Z = [ 1 ] + E = [ {energy} ] + Ec = [ 0.1 ] + Es = [ 1.5 ] + interaction_index = [ 0 ] + pos = [ [ -4.4, 0.0, 0.0,] ] + dir = [ [ 0.9999999999984769, 1.7453292519934434e-6, 0.0,] ] + + [geometry_input] + length_unit = "ANGSTROM" + electronic_stopping_correction_factor = 1.09 + densities = [ 0.0914 ] + + [material_parameters] + energy_unit = "EV" + mass_unit = "AMU" + Eb = [ 0.0 ] + Es = [ 5.61 ] + Ec = [ 3.0 ] + Z = [ 28 ] + m = [ 58.69 ] + interaction_index = [ 0 ] + surface_binding_model = {{"PLANAR"={{calculation="INDIVIDUAL"}}}} + bulk_binding_model = "AVERAGE" + ''' + + with open(f'krc_morse_{index}.toml', 'w') as file: + file.write(input_file) + + if run_sim: os.system(f'cargo run --release --features cpr_rootfinder 0D krc_morse_{index}.toml') + + reflected_list = np.atleast_2d(np.genfromtxt(f'krc_morse_{index}reflected.output', delimiter=',')) + if np.size(reflected_list) > 0: + num_reflected = np.shape(reflected_list)[0] + energy_reflected = np.sum(reflected_list[:, 2]) + else: + num_reflected = 0 + energy_reflected = 0. + + return num_reflected/num_samples, energy_reflected/energy/num_samples + #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 +#Used to show that TRIM breaks down at lower incident energies data = np.array( [[0.2962771, 0.18217821], [0.96344626, 0.4871287], @@ -117,26 +171,28 @@ def run_morse_potential(energy, index, num_samples=10000, run_sim=True): #Running and plotting the H-Ni simulations with the Morse potential and updated Es num_energies = 15 -energies = np.logspace(-1, 3, num_energies) +energies = np.logspace(-1, 4, num_energies) run_sim = True -num_samples = 1000 +num_samples = 10000 R_N = np.zeros(num_energies) R_E = np.zeros(num_energies) +R_N_2 = np.zeros(num_energies) +R_E_2 = np.zeros(num_energies) for index, energy in enumerate(energies): - R_N[index], R_E[index] = run_morse_potential(energy, index, num_samples=num_samples, run_sim=run_sim) - -plt.semilogx(energies, R_N, label='R_N Morse H-Ni, Es=2.5eV', color='purple') + R_N[index], R_E[index] = run_krc_morse_potential(energy, index, num_samples=num_samples, run_sim=False) + R_N_2[index], R_E_2[index] = run_morse_potential(energy, index, num_samples=num_samples, run_sim=False) +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, 40) +energies = np.logspace(-1, 4, 50) r_rustbca = np.array([reflection_coefficient(hydrogen, nickel, energy, 0.0, 10000) 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') - r_thomas = [thomas_reflection(hydrogen, nickel, energy) for energy in energies] plt.semilogx(energies, r_thomas, label='Thomas', color='red', linestyle='--') @@ -144,5 +200,5 @@ def run_morse_potential(energy, index, num_samples=10000, run_sim=True): plt.title('Reflection Coefficients H+ on Ni') plt.xlabel('E [eV]') plt.ylabel('R') - -plt.savefig('Morse.png') \ No newline at end of file +plt.show() +plt.savefig('Morse.png') diff --git a/examples/test_rustbca.f90 b/examples/test_rustbca.f90 index c0b8cae..991184b 100644 --- a/examples/test_rustbca.f90 +++ b/examples/test_rustbca.f90 @@ -15,7 +15,7 @@ program test_rustbca logical(c_bool) :: track_recoils !Initial ion conditions - N_ions = 100000 + N_ions = 1000 allocate(ux(N_ions), uy(N_ions), uz(N_ions), E(N_ions), Z1(N_ions), m1(N_ions), Ec1(N_ions), Es1(N_ions)) ux(:) = 0.999 uy(:) = sqrt(1.0 - 0.999*0.999) @@ -52,11 +52,12 @@ program test_rustbca num_species_target, Z2, m2, Ec2, Es2, Eb2, n2, & num_emitted_particles) call c_f_pointer(bca_output_c, bca_output_f, [num_emitted_particles, 6]) + call cpu_time(stop) write(*,*) "Elapsed time in seconds per ion per eV: ", (stop - start)/N_ions/1000.0 - !write(*,*) bca_output_f + deallocate(bca_output_f) call cpu_time(start) do i = 0, N_ions @@ -70,6 +71,6 @@ program test_rustbca call cpu_time(stop) write(*,*) "Elapsed time in ions per eV per s: ", (stop - start)/N_ions/1000.0 - !call exit(1) - + deallocate(ux, uy, uz, E, Z1, m1, Ec1, Es1) + end program test_rustbca diff --git a/examples/test_rustbca.py b/examples/test_rustbca.py index dd8b58e..3a8ca74 100644 --- a/examples/test_rustbca.py +++ b/examples/test_rustbca.py @@ -11,8 +11,33 @@ import time def main(): + + #test rotation to and from RustBCA coordinates + + #nx, ny, nz is the normal vector (out of surface) + nx = -np.sqrt(2)/2 + ny = -np.sqrt(2)/2 + nz = 0.0 + + #ux, uy, uz is the particle direction (simulation coordinates) + ux = 1.0 + uy = 0.0 + uz = 0.0 + + print(f'Before rotation: ({ux}, {uy}, {uz})') + ux, uy, uz = rotate_given_surface_normal_py(nx, ny, nz, ux, uy, uz) + print(f'After rotation: ({ux}, {uy}, {uz})') + ux, uy, uz = rotate_back_py(nx, ny, nz, ux, uy, uz) + print(f'After rotation back: ({ux}, {uy}, {uz})') + + #After rotating and rotating back, effect should be where you started (minus fp error) + assert(abs(ux - 1.0) < 1e-6) + assert(abs(uy) < 1e-6) + assert(abs(uz) < 1e-6) + #scripts/materials.py has a number of potential ions and targets ion = helium + ion['Eb'] = 0.0 target = tungsten #tungsten Eb = 0 is a 'worst-case' - accepted literature value is 3 eV @@ -23,16 +48,27 @@ def main(): num_samples = 1000 #For automatic control, libRustBCA offers two helpful functions for the sputtering yield and reflection coefficients Y = sputtering_yield(ion, target, energy, angle, num_samples) - + print(f'Sputtering yield for {ion["symbol"]} on {target["symbol"]} at {energy} eV is {Y} at/ion. Yamamura predicts { np.round(yamamura(ion, target, energy),3)} at/ion.') R_N, R_E = reflection_coefficient(ion, target, energy, angle, num_samples) print(f'Particle reflection coefficient for {ion["symbol"]} on {target["symbol"]} at {energy} eV is {R_N}. Thomas predicts {np.round(thomas_reflection(ion, target, energy), 3)}.') print(f'Energy reflection coefficient for {ion["symbol"]} on {target["symbol"]} at {energy} eV is {R_E}') + R_N, R_E = compound_reflection_coefficient(ion, [target, ion], [target['n'], 0.1*target['n']], energy, angle, num_samples) + print(f'Particle reflection coefficient for {ion["symbol"]} on {ion["symbol"]}x{target["symbol"]} where x=0.1 at {energy} eV is {R_N}. Thomas predicts {np.round(thomas_reflection(ion, target, energy), 3)}.') + print(f'Energy reflection coefficient for {ion["symbol"]}x{target["symbol"]} where x=0.1 at {energy} eV is {R_E}') + + vx0 = 1e5 + vy0 = 1e5 + vz0 = 0.0 + vx1, vy1, vz1 = reflect_single_ion_py(ion, target, vx0, vy0, vz0) + + print(f'(vx, vy, vz) before reflection: ({vx0}, {vy0}, {vz0})') + print(f'(vx, vy, vz) after reflection: ({vx1}, {vy1}, {vz1})') #For smooth distributions and good statistics, you should use at least 10k ions - number_ions = 100000 + number_ions = 10000 #1 keV is above the He on W sputtering threshold of ~150 eV energies_eV = 1000.0*np.ones(number_ions) @@ -136,7 +172,7 @@ def main(): heights, _, _ = plt.hist(x[np.logical_and(incident, stopped)], bins=100, density=True, histtype='step') plt.plot([50.0, 50.0], [0.0, np.max(heights)*1.1]) plt.gca().set_ylim([0.0, np.max(heights)*1.1]) - + plt.show() if __name__ == '__main__': diff --git a/rustbca.f90 b/rustbca.f90 index 681a78c..e1a2c0d 100644 --- a/rustbca.f90 +++ b/rustbca.f90 @@ -1,7 +1,7 @@ module rustbca use, intrinsic :: iso_c_binding - + interface subroutine reflect_single_ion_c(num_species_target, ux, uy, uz, E1, & @@ -9,7 +9,7 @@ subroutine reflect_single_ion_c(num_species_target, ux, uy, uz, E1, & !Runs a single ion BCA trajectory with no recoils !Args: - ! num_species_target (integer): number of species in target + ! num_species_target (integer): number of species in target ! ux (real(c_double)): x-direction of incident ion, x-direction of reflected ion ! uy (real(c_double)): y-direction of incident ion, y-direction of reflected ion ! uz (real(c_double)): z-direction of incident ion, z-direction of reflected ion @@ -33,12 +33,28 @@ subroutine reflect_single_ion_c(num_species_target, ux, uy, uz, E1, & end subroutine reflect_single_ion_c + subroutine rotate_given_surface_normal(nx, ny, nz, ux, uy, uz) bind(c) + + use, intrinsic :: iso_c_binding + real(c_double), intent(in) :: nx, ny, nz + real(c_double), intent(inout) :: ux, uy, uz + + end subroutine rotate_given_surface_normal + + subroutine rotate_back(nx, ny, nz, ux, uy, uz) bind(c) + + use, intrinsic :: iso_c_binding + real(c_double), intent(in) :: nx, ny, nz + real(c_double), intent(inout) :: ux, uy, uz + + end subroutine rotate_back + function compound_bca_list_fortran(num_incident_ions, track_recoils, ux, uy, uz, E1, & Z1, m1, Ec1, Es1, & num_species_target, Z2, m2, Ec2, Es2, Eb2, n2, & num_emitted_particles) bind(c) result(output) - + !Runs a homogeneous, flat, compound target BCA with an arbitrary list of ions. !Args: ! num_incident_ion (integer(c_int)): number of incident ions @@ -92,4 +108,4 @@ subroutine transform_to_local_angle(ux, uy, uz, alpha) end subroutine transform_to_local_angle -end module rustbca \ No newline at end of file +end module rustbca diff --git a/scripts/materials.py b/scripts/materials.py index ad13552..b8522fe 100644 --- a/scripts/materials.py +++ b/scripts/materials.py @@ -12,6 +12,29 @@ 's': 2.5 } +carbon = { + 'symbol': 'C', + 'name': 'carbon', + 'Z': 6, + 'm': 12.011, + 'n': 1.1331E29, + 'Es': 7.37, + 'Ec': 1.0, + 'Eb': 3.0, +} + +chromium = { + 'symbol': 'Cr', + 'name': 'chromium', + 'Z': 24, + 'm': 51.9961, + 'Es': 4.10, + 'Ed': 28.0, + 'Ec': 3.0, + 'Eb': 3.0, + 'n': 8.327E28 +} + hydrogen = { 'symbol': 'H', 'name': 'hydrogen', diff --git a/setup.py b/setup.py index e33e419..21ff36d 100644 --- a/setup.py +++ b/setup.py @@ -3,13 +3,13 @@ setup( name="RustBCA", - version="2.0.1", + version="2.8.4", rust_extensions=[ RustExtension( "libRustBCA", binding=Binding.PyO3, features=["python", "parry3d"], - #args=["+nightly", "--edition 2018", "-Z unstable-options"], + #args=["+nightly", "--edition 2021", "-Z unstable-options"], #optional=True, #rust_version="1.57.0" ) diff --git a/src/bca.rs b/src/bca.rs index 32150af..7d3f013 100644 --- a/src/bca.rs +++ b/src/bca.rs @@ -1,6 +1,6 @@ use super::*; -#[cfg(any(feature = "cpr_rootfinder_openblas", feature = "cpr_rootfinder_netlib", feature = "cpr_rootfinder_intel_mkl"))] +#[cfg(feature = "cpr_rootfinder")] use rcpr::chebyshev::*; /// Geometrical quantities of binary collision. @@ -89,7 +89,7 @@ pub fn single_ion_bca(particle: particle::Particle, material: &mate 0. }; - let mut total_energy_loss = 0.; + let mut total_energy_lost_to_recoils = 0.; let mut total_asymptotic_deflection = 0.; let mut normalized_distance_of_closest_approach = 0.; let mut strong_collision_Z = 0.; @@ -123,7 +123,8 @@ pub fn single_ion_bca(particle: particle::Particle, material: &mate particle_2.energy_origin = particle_2.E; //Accumulate energy losses and asymptotic deflections for primary particle - total_energy_loss += binary_collision_result.recoil_energy; + particle_1.E += -binary_collision_result.recoil_energy; + total_energy_lost_to_recoils += binary_collision_result.recoil_energy; total_asymptotic_deflection += binary_collision_result.asymptotic_deflection; particle_1.rotate(binary_collision_result.psi, @@ -184,9 +185,10 @@ pub fn single_ion_bca(particle: particle::Particle, material: &mate binary_collision_geometries[0].mfp + distance_to_target - material.geometry.get_energy_barrier_thickness(), total_asymptotic_deflection); //Subtract total energy from all simultaneous collisions and electronic stopping - bca::update_particle_energy(&mut particle_1, &material, distance_traveled, - total_energy_loss, normalized_distance_of_closest_approach, strong_collision_Z, + let energy_lost_to_electronic_stopping = bca::subtract_electronic_stopping_energy(&mut particle_1, &material, distance_traveled, + normalized_distance_of_closest_approach, strong_collision_Z, strong_collision_index, &options); + particle_1.update_energy_loss_tracker(&options, total_energy_lost_to_recoils, energy_lost_to_electronic_stopping); //Check boundary conditions on leaving and stopping material::boundary_condition_planar(&mut particle_1, &material); @@ -364,9 +366,9 @@ pub fn choose_collision_partner(particle_1: &particle::Particle, ma let z_recoil: f64 = z + mfp*cosz + impact_parameter*(sinphi*cosy - cosphi*cosx*cosz)/sinx; //Choose recoil Z, M - let (species_index, Z_recoil, M_recoil, Ec_recoil, Es_recoil, interaction_index) = material.choose(x_recoil, y_recoil, z_recoil); + let (species_index, Z_recoil, M_recoil, Ec_recoil, Es_recoil, Ed_recoil, interaction_index) = material.choose(x_recoil, y_recoil, z_recoil); let mut new_particle = particle::Particle::new( - M_recoil, Z_recoil, 0., Ec_recoil, Es_recoil, + M_recoil, Z_recoil, 0., Ec_recoil, Es_recoil, Ed_recoil, x_recoil, y_recoil, z_recoil, cosx, cosy, cosz, false, options.track_recoil_trajectories, interaction_index @@ -401,7 +403,7 @@ fn distance_of_closest_approach(particle_1: &particle::Particle, particle_2: &pa options.root_finder[particle_1.interaction_index][particle_2.interaction_index] } else {Rootfinder::NEWTON{max_iterations: 100, tolerance: 1E-6}}; - #[cfg(any(feature = "cpr_rootfinder_openblas", feature = "cpr_rootfinder_netlib", feature = "cpr_rootfinder_intel_mkl"))] + #[cfg(feature = "cpr_rootfinder")] match root_finder { Rootfinder::POLYNOMIAL{complex_threshold} => polynomial_rootfinder(Za, Zb, Ma, Mb, E0, p, interaction_potential, complex_threshold) .with_context(|| format!("Numerical error: Polynomial rootfinder failed for {} at {} eV with p = {} A.", interaction_potential, E0/EV, p/ANGSTROM)) @@ -418,7 +420,7 @@ fn distance_of_closest_approach(particle_1: &particle::Particle, particle_2: &pa .unwrap(), } - #[cfg(not(any(feature = "cpr_rootfinder_openblas", feature = "cpr_rootfinder_netlib", feature = "cpr_rootfinder_intel_mkl")))] + #[cfg(not(feature = "cpr_rootfinder"))] match root_finder { Rootfinder::NEWTON{max_iterations, tolerance} => newton_rootfinder(Za, Zb, Ma, Mb, E0, p, interaction_potential, max_iterations, tolerance) .with_context(|| format!("Numerical error: Newton rootfinder failed for {} at {} eV with p = {} A.", interaction_potential, E0/EV, p/ANGSTROM)) @@ -431,16 +433,10 @@ fn distance_of_closest_approach(particle_1: &particle::Particle, particle_2: &pa } /// Update a particle's energy following nuclear and electronic interactions of a single BCA step. -pub fn update_particle_energy(particle_1: &mut particle::Particle, material: &material::Material, distance_traveled: f64, - recoil_energy: f64, x0: f64, strong_collision_Z: f64, strong_collision_index: usize, options: &Options) { +pub fn subtract_electronic_stopping_energy(particle_1: &mut particle::Particle, material: &material::Material, distance_traveled: f64, + x0: f64, strong_collision_Z: f64, strong_collision_index: usize, options: &Options) -> f64 { - //If particle energy drops below zero before electronic stopping calcualtion, it produces NaNs - assert!(!recoil_energy.is_nan(), "Numerical error: recoil Energy is NaN. E0 = {} Za = {} Ma = {} x0 = {} Zb = {} delta-x = {}", particle_1.E, particle_1.Z, particle_1.m, x0, strong_collision_Z, distance_traveled); - particle_1.E -= recoil_energy; assert!(!particle_1.E.is_nan(), "Numerical error: particle energy is NaN following collision."); - if particle_1.E < 0. { - particle_1.E = 0.; - } let x = particle_1.pos.x; let y = particle_1.pos.y; @@ -452,7 +448,7 @@ pub fn update_particle_energy(particle_1: &mut particle::Particle, let electronic_stopping_powers = material.electronic_stopping_cross_sections(particle_1, options.electronic_stopping_mode); let n = material.number_densities(x, y, z); - let delta_energy = match options.electronic_stopping_mode { + let delta_energy_electronic = match options.electronic_stopping_mode { ElectronicStoppingMode::INTERPOLATED => electronic_stopping_powers.iter().zip(n).map(|(se, number_density)| se*number_density).collect::>().iter().sum::()*distance_traveled, ElectronicStoppingMode::LOW_ENERGY_NONLOCAL => electronic_stopping_powers.iter().zip(n).map(|(se, number_density)| se*number_density).collect::>().iter().sum::()*distance_traveled, ElectronicStoppingMode::LOW_ENERGY_LOCAL => oen_robinson_loss(particle_1.Z, strong_collision_Z, electronic_stopping_powers[strong_collision_index], x0, interaction_potential), @@ -465,16 +461,15 @@ pub fn update_particle_energy(particle_1: &mut particle::Particle, }, }; - particle_1.E += -delta_energy; + particle_1.E += -delta_energy_electronic; //Make sure particle energy doesn't become negative again assert!(!particle_1.E.is_nan(), "Numerical error: particle energy is NaN following electronic stopping."); if particle_1.E < 0. { particle_1.E = 0.; } - - particle_1.energy_loss(&options, recoil_energy, delta_energy); - } else if recoil_energy > 0. { - particle_1.energy_loss(&options, recoil_energy, 0.); + delta_energy_electronic + } else { + 0.0 } } @@ -569,7 +564,7 @@ fn scattering_integral_gauss_legendre(impact_parameter: f64, relative_energy: f6 .unwrap()).sum::() } -#[cfg(any(feature = "cpr_rootfinder_openblas", feature = "cpr_rootfinder_netlib", feature = "cpr_rootfinder_intel_mkl"))] +#[cfg(feature = "cpr_rootfinder")] /// Computes the distance of closest approach of two particles with atomic numbers `Za`, `Zb` and masses `Ma`, `Mb` for an inverse-polynomial interaction potential (e.g., Lennard-Jones) for a given impact parameter and incident energy `E0`. /// Slightly complex roots with an imaginary part smaller than `polynom_complex_threshold` are considered real roots. pub fn polynomial_rootfinder(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parameter: f64, @@ -581,7 +576,21 @@ pub fn polynomial_rootfinder(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact let coefficients = interactions::polynomial_coefficients(relative_energy, impact_parameter, interaction_potential); let roots = real_polynomial_roots(coefficients.clone(), polynom_complex_threshold).unwrap(); + + /* + println!("p={}", impact_parameter/ANGSTROM); + + for coefficient in &coefficients { + println!("{}", {coefficient}); + } + + for root in &roots { + println!("{} A", {root}); + } + */ + let max_root = roots.iter().cloned().fold(f64::NAN, f64::max); + let inverse_transformed_root = interactions::inverse_transform(max_root, interaction_potential); if roots.is_empty() || inverse_transformed_root.is_nan() { @@ -592,7 +601,7 @@ pub fn polynomial_rootfinder(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact } } -#[cfg(any(feature = "cpr_rootfinder_openblas", feature = "cpr_rootfinder_netlib", feature = "cpr_rootfinder_intel_mkl"))] +#[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. /// /// # Args: @@ -632,13 +641,13 @@ pub fn cpr_rootfinder(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parame let upper_bound = impact_parameter + interactions::crossing_point_doca(interaction_potential); let roots = match derivative_free { - true => find_roots_with_secant_polishing(&g, &f, 0., upper_bound, + true => find_roots_with_secant_polishing(&g, &f, 1e-15, 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, 0., upper_bound, + find_roots_with_newton_polishing(&g, &f, &df, 1e-15, upper_bound, n0, epsilon, nmax, complex_threshold, truncation_threshold, interval_limit, far_from_zero) } @@ -715,7 +724,7 @@ pub fn magic(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parameter: f64, //Since this is legacy code I don't think I will clean this up let C_ = match interaction_potential { InteractionPotential::MOLIERE => vec![ 0.6743, 0.009611, 0.005175, 6.314, 10.0 ], - InteractionPotential::KR_C => vec![ 0.7887, 0.01166, 00.006913, 17.16, 10.79 ], + InteractionPotential::KR_C => vec![ 0.7887, 0.01166, 0.006913, 17.16, 10.79 ], InteractionPotential::ZBL => vec![ 0.99229, 0.011615, 0.0071222, 9.3066, 14.813 ], InteractionPotential::TRIDYN => vec![1.0144, 0.235809, 0.126, 69350., 83550.], //Undocumented Tridyn constants _ => panic!("Input error: unimplemented interaction potential {} for MAGIC algorithm. Use a screened Coulomb potential.", interaction_potential) diff --git a/src/enums.rs b/src/enums.rs index 272e923..908c040 100644 --- a/src/enums.rs +++ b/src/enums.rs @@ -8,7 +8,8 @@ pub enum MaterialType { #[cfg(feature = "parry3d")] BALL(material::Material), #[cfg(feature = "parry3d")] - TRIMESH(material::Material) + TRIMESH(material::Material), + HOMOGENEOUS2D(material::Material) } #[derive(Deserialize, PartialEq, Clone, Copy)] @@ -29,6 +30,7 @@ pub enum GeometryType { BALL, #[cfg(feature = "parry3d")] TRIMESH, + HOMOGENEOUS2D, } /// Mode of electronic stopping to use. @@ -169,7 +171,11 @@ pub enum InteractionPotential { /// Tungsten-tungsten cubic spline potential (Following Bjorkas et al.) WW, /// Unscreened Coulombic interatomic potential between ions with charges Za and Zb. - COULOMB{Za: f64, Zb: f64} + COULOMB{Za: f64, Zb: f64}, + /// Morse potential that smoothly interpolates to Kr-C at short ranges. + KRC_MORSE{D: f64, alpha: f64, r0: f64, k: f64, x0: f64}, + /// Born repulsive potential. + FOUR_EIGHT{alpha: f64, beta: f64} } impl fmt::Display for InteractionPotential { @@ -184,7 +190,9 @@ impl fmt::Display for InteractionPotential { InteractionPotential::LENNARD_JONES_65_6{sigma, epsilon} => write!(f, "Lennard-Jones 6.5-6 Potential with sigma = {} A, epsilon = {} eV", sigma/ANGSTROM, epsilon/EV), InteractionPotential::MORSE{D, alpha, r0} => write!(f, "Morse potential with D = {} eV, alpha = {} 1/A, and r0 = {} A", D/EV, alpha*ANGSTROM, r0/ANGSTROM), InteractionPotential::WW => write!(f, "W-W cubic spline interaction potential."), - InteractionPotential::COULOMB{Za, Zb} => write!(f, "Coulombic interaction with Za = {} and Zb = {}", Za, Zb) + InteractionPotential::COULOMB{Za, Zb} => write!(f, "Coulombic interaction with Za = {} and Zb = {}", Za, Zb), + InteractionPotential::KRC_MORSE{D, alpha, r0, k, x0} => write!(f, "Morse potential with D = {} eV, alpha = {} 1/A, and r0 = {} A with short-range Kr-C with interpolation width k = {} 1/A and transition point x0 = {} A", D/EV, alpha*ANGSTROM, r0/ANGSTROM, k*ANGSTROM, x0/ANGSTROM), + InteractionPotential::FOUR_EIGHT{alpha, beta} => write!(f, "Inverse polynomial potential with beta/r^8 (beta = {} A) Born repulsive term and alpha/r^4 (alpha = {} A) attractive term. Used to describe hydrogen interactions with transition metals.", alpha/ANGSTROM, beta/ANGSTROM), } } } diff --git a/src/geometry.rs b/src/geometry.rs index 2a11cbd..8e3d31a 100644 --- a/src/geometry.rs +++ b/src/geometry.rs @@ -265,6 +265,120 @@ impl Geometry for Mesh1D { } } +#[derive(Deserialize, Clone)] +pub struct HomogeneousMesh2DInput { + pub length_unit: String, + pub points: Vec<(f64, f64)>, + pub simulation_boundary_points: Vec<(f64, f64)>, + pub densities: Vec, + pub electronic_stopping_correction_factor: f64 +} + +#[derive(Clone)] +pub struct HomogeneousMesh2D { + pub boundary: Polygon, + pub simulation_boundary: Polygon, + pub energy_barrier_thickness: f64, + pub densities: Vec, + pub electronic_stopping_correction_factor: f64, + pub concentrations: Vec +} + +impl Geometry for HomogeneousMesh2D { + + type InputFileFormat = InputHomogeneous2D; + + fn new(input: &<::InputFileFormat as GeometryInput>::GeometryInput) -> Self { + let length_unit: f64 = match input.length_unit.as_str() { + "MICRON" => MICRON, + "CM" => CM, + "MM" => MM, + "ANGSTROM" => ANGSTROM, + "NM" => NM, + "M" => 1., + _ => input.length_unit.parse() + .expect(format!( + "Input errror: could nor parse length unit {}. Use a valid float or one of ANGSTROM, NM, MICRON, CM, MM, M", + &input.length_unit.as_str() + ).as_str()), + }; + + let boundary_points_converted: Vec<(f64, f64)> = input.points.iter().map(|(x, y)| (x*length_unit, y*length_unit)).collect(); + + let simulation_boundary_points_converted: Vec<(f64, f64)> = input.simulation_boundary_points.iter().map(|(x, y)| (x*length_unit, y*length_unit)).collect(); + + let electronic_stopping_correction_factor = input.electronic_stopping_correction_factor; + + let densities: Vec = input.densities.iter().map(|element| element/(length_unit).powi(3)).collect(); + + let total_density: f64 = densities.iter().sum(); + + let energy_barrier_thickness = total_density.powf(-1./3.)/SQRTPI*2.; + + let concentrations: Vec = densities.iter().map(|&density| density/total_density).collect::>(); + + let boundary = Polygon::new(LineString::from(boundary_points_converted), vec![]); + + let simulation_boundary = Polygon::new(LineString::from(simulation_boundary_points_converted), vec![]); + + HomogeneousMesh2D { + densities, + simulation_boundary, + boundary, + electronic_stopping_correction_factor, + energy_barrier_thickness, + concentrations + } + } + + fn get_densities(&self, x: f64, y: f64, z: f64) -> &Vec { + &self.densities + } + + fn get_ck(&self, x: f64, y: f64, z: f64) -> f64 { + self.electronic_stopping_correction_factor + } + + fn get_total_density(&self, x: f64, y: f64, z: f64) -> f64 { + self.densities.iter().sum::() + } + + fn get_concentrations(&self, x: f64, y: f64, z: f64) -> &Vec { + &self.concentrations + } + + fn inside(&self, x: f64, y: f64, z: f64) -> bool { + self.boundary.contains(&point!(x: x, y: y)) + } + + fn inside_simulation_boundary(&self, x: f64, y: f64, z: f64) -> bool { + self.simulation_boundary.contains(&point!(x: x, y: y)) + } + fn inside_energy_barrier(&self, x: f64, y: f64, z: f64) -> bool { + if self.inside(x, y, z) { + true + } else { + if let Closest::SinglePoint(p) = self.boundary.closest_point(&point!(x: x, y: y)) { + let distance = ((x - p.x()).powf(2.) + (y - p.y()).powf(2.)).sqrt(); + distance < self.energy_barrier_thickness + } else if let Closest::Intersection(p) = self.boundary.closest_point(&point!(x: x, y: y)) { + true + } else { + panic!("Geometry error: closest point routine failed to find single closest point to ({}, {}, {}).", x, y, z); + } + } + } + fn closest_point(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) { + if let Closest::SinglePoint(p) = self.boundary.closest_point(&point!(x: x, y: y)) { + (p.x(), p.y(), z) + } else if let Closest::Intersection(p) = self.boundary.closest_point(&point!(x: x, y: y)) { + return (p.x(), p.y(), z) + } else { + panic!("Geometry error: closest point routine failed to find single closest point to ({}, {}, {}).", x, y, z); + } + } +} + /// Object that contains raw mesh input data. #[derive(Deserialize, Clone)] pub struct Mesh2DInput { @@ -434,6 +548,8 @@ impl Geometry for Mesh2D { fn closest_point(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) { if let Closest::SinglePoint(p) = self.boundary.closest_point(&point!(x: x, y: y)) { (p.x(), p.y(), z) + } else if let Closest::Intersection(p) = self.boundary.closest_point(&point!(x: x, y: y)) { + (p.x(), p.y(), z) } else { panic!("Geometry error: closest point routine failed to find single closest point to ({}, {}, {}).", x, y, z); } diff --git a/src/input.rs b/src/input.rs index e6d5ad9..705c333 100644 --- a/src/input.rs +++ b/src/input.rs @@ -49,6 +49,42 @@ impl InputFile for Input2D { } } +/// Rustbca's internal representation of an input file for homogeneous 2D meshes. +#[derive(Deserialize, Clone)] +pub struct InputHomogeneous2D { + pub options: Options, + pub material_parameters: material::MaterialParameters, + pub particle_parameters: particle::ParticleParameters, + pub geometry_input: geometry::HomogeneousMesh2DInput, +} + +impl GeometryInput for InputHomogeneous2D { + type GeometryInput = geometry::HomogeneousMesh2DInput; +} + +impl InputFile for InputHomogeneous2D { + + fn new(string: &str) -> InputHomogeneous2D { + toml::from_str(string).context( + "Could not parse TOML file. Be sure you are using the correct input file mode (e.g., + ./RustBCA SPHERE sphere.toml or RustBCA.exe 0D mesh_0d.toml)." + ).unwrap() + } + + fn get_options(&self) -> &Options{ + &self.options + } + fn get_material_parameters(&self) -> &material::MaterialParameters{ + &self.material_parameters + } + fn get_particle_parameters(&self) -> &particle::ParticleParameters{ + &self.particle_parameters + } + fn get_geometry_input(&self) -> &Self::GeometryInput{ + &self.geometry_input + } +} + #[derive(Deserialize, Clone)] pub struct Input0D { pub options: Options, @@ -372,6 +408,10 @@ where ::InputFileFormat: Deserialize<'static> + 'static { material.interaction_index = vec![0; material.m.len()]; } + if material.Ed.len() <= 1 { + material.Ed = vec![material.Ed[0]; material.m.len()]; + } + //Check that incompatible options are not on simultaneously if options.high_energy_free_flight_paths { assert!(options.electronic_stopping_mode == ElectronicStoppingMode::INTERPOLATED, @@ -399,7 +439,7 @@ where ::InputFileFormat: Deserialize<'static> + 'static { for ((interaction_potentials, scattering_integrals), root_finders) in options.interaction_potential.iter().zip(options.scattering_integral.clone()).zip(options.root_finder.clone()) { for ((interaction_potential, scattering_integral), root_finder) in interaction_potentials.iter().zip(scattering_integrals).zip(root_finders) { - if cfg!(not(any(feature="cpr_rootfinder_openblas", feature="cpr_rootfinder_netlib", feature="cpr_rootfinder_intel_mkl",))) { + if cfg!(not(feature="cpr_rootfinder")) { assert!( match root_finder { Rootfinder::POLYNOMIAL{..} => false, Rootfinder::CPR{..} => false, @@ -410,6 +450,7 @@ where ::InputFileFormat: Deserialize<'static> + 'static { assert!( match (interaction_potential, root_finder) { + (InteractionPotential::FOUR_EIGHT{..}, Rootfinder::POLYNOMIAL{..}) => true, (InteractionPotential::LENNARD_JONES_12_6{..}, Rootfinder::CPR{..}) => true, (InteractionPotential::LENNARD_JONES_12_6{..}, Rootfinder::POLYNOMIAL{..}) => true, (InteractionPotential::LENNARD_JONES_12_6{..}, _) => false, @@ -523,7 +564,7 @@ where ::InputFileFormat: Deserialize<'static> + 'static { E: match E { Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rand::thread_rng())*energy_unit}, Distributions::UNIFORM{min, max} => {let uniform = Uniform::from(min..max); uniform.sample(&mut rand::thread_rng())*energy_unit}, - Distributions::POINT(x) => x*energy_unit, + Distributions::POINT(E) => E*energy_unit, }, Ec: Ec*energy_unit, Es: Es*energy_unit, @@ -535,27 +576,27 @@ where ::InputFileFormat: Deserialize<'static> + 'static { y: match y { Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rand::thread_rng())*length_unit}, Distributions::UNIFORM{min, max} => {let uniform = Uniform::from(min..max); uniform.sample(&mut rand::thread_rng())*length_unit}, - Distributions::POINT(x) => x*length_unit, + Distributions::POINT(y) => y*length_unit, }, z: match z { Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rand::thread_rng())*length_unit}, Distributions::UNIFORM{min, max} => {let uniform = Uniform::from(min..max); uniform.sample(&mut rand::thread_rng())*length_unit}, - Distributions::POINT(x) => x*length_unit, + Distributions::POINT(z) => z*length_unit, }, ux: match cosx { Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rand::thread_rng())*length_unit}, Distributions::UNIFORM{min, max} => {let uniform = Uniform::from(min..max); uniform.sample(&mut rand::thread_rng())*length_unit}, - Distributions::POINT(x) => {assert!(x != 1.0, "ux cannot equal exactly 1.0 to prevent gimbal lock"); x*length_unit} + Distributions::POINT(ux) => {assert!(ux != 1.0, "ux cannot equal exactly 1.0 to prevent gimbal lock"); ux} }, uy: match cosy { Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rand::thread_rng())*length_unit}, Distributions::UNIFORM{min, max} => {let uniform = Uniform::from(min..max); uniform.sample(&mut rand::thread_rng())*length_unit}, - Distributions::POINT(x) => x*length_unit, + Distributions::POINT(uy) => uy, }, uz: match cosz { Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rand::thread_rng())*length_unit}, Distributions::UNIFORM{min, max} => {let uniform = Uniform::from(min..max); uniform.sample(&mut rand::thread_rng())*length_unit}, - Distributions::POINT(x) => x*length_unit, + Distributions::POINT(uz) => uz, }, interaction_index: interaction_index, tag: 0, diff --git a/src/interactions.rs b/src/interactions.rs index db965a9..b050294 100644 --- a/src/interactions.rs +++ b/src/interactions.rs @@ -12,6 +12,10 @@ pub fn crossing_point_doca(interaction_potential: InteractionPotential) -> f64 { } +fn sigmoid(x: f64, k: f64, x0: f64) -> f64 { + 1./(1. + (-k*(x - x0)).exp()) +} + /// Interaction potential between two particles a and b at a distance `r`. pub fn interaction_potential(r: f64, a: f64, Za: f64, Zb: f64, interaction_potential: InteractionPotential) -> f64 { match interaction_potential { @@ -23,16 +27,23 @@ pub fn interaction_potential(r: f64, a: f64, Za: f64, Zb: f64, interaction_poten }, InteractionPotential::LENNARD_JONES_65_6{sigma, epsilon} => { lennard_jones_65_6(r, sigma, epsilon) - } + }, InteractionPotential::MORSE{D, alpha, r0} => { morse(r, D, alpha, r0) - } + }, InteractionPotential::WW => { tungsten_tungsten_cubic_spline(r) - } + }, InteractionPotential::COULOMB{Za, Zb} => { coulomb(r, Za, Zb) + }, + InteractionPotential::KRC_MORSE{D, alpha, r0, k, x0} => { + krc_morse(r, a, Za, Zb, D, alpha, r0, k, x0) + } + InteractionPotential::FOUR_EIGHT{alpha, beta} => { + four_eight(r, alpha, beta) } + } } @@ -40,7 +51,8 @@ pub fn interaction_potential(r: f64, a: f64, Za: f64, Zb: f64, interaction_poten pub fn energy_threshold_single_root(interaction_potential: InteractionPotential) -> f64 { match interaction_potential{ InteractionPotential::LENNARD_JONES_12_6{..} | InteractionPotential::LENNARD_JONES_65_6{..} => f64::INFINITY, - InteractionPotential::MORSE{..} => f64::INFINITY, + InteractionPotential::MORSE{..} | InteractionPotential::KRC_MORSE{..} => f64::INFINITY, + InteractionPotential::FOUR_EIGHT{..} => f64::INFINITY, InteractionPotential::WW => f64::INFINITY, InteractionPotential::COULOMB{..} => f64::INFINITY, InteractionPotential::MOLIERE | InteractionPotential::KR_C | InteractionPotential::LENZ_JENSEN | InteractionPotential::ZBL | InteractionPotential::TRIDYN => 0., @@ -93,6 +105,12 @@ pub fn distance_of_closest_approach_function(r: f64, a: f64, Za: f64, Zb: f64, r InteractionPotential::WW => { doca_tungsten_tungsten_cubic_spline(r, impact_parameter, relative_energy) }, + InteractionPotential::KRC_MORSE{D, alpha, r0, k, x0} => { + doca_krc_morse(r, impact_parameter, relative_energy, a, Za, Zb, D, alpha, r0, k, x0) + }, + InteractionPotential::FOUR_EIGHT{alpha, beta} => { + doca_four_eight(r, impact_parameter, relative_energy, alpha, beta) + } InteractionPotential::COULOMB{..} => panic!("Coulombic potential cannot be used with rootfinder.") } } @@ -118,6 +136,12 @@ pub fn distance_of_closest_approach_function_singularity_free(r: f64, a: f64, Za InteractionPotential::MORSE{D, alpha, r0} => { doca_morse(r, impact_parameter, relative_energy, D, alpha, r0) }, + InteractionPotential::KRC_MORSE{D, alpha, r0, k, x0} => { + doca_krc_morse(r, impact_parameter, relative_energy, a, Za, Zb, D, alpha, r0, k, x0) + }, + InteractionPotential::FOUR_EIGHT{alpha, beta} => { + doca_four_eight(r, impact_parameter, relative_energy, alpha, beta) + } InteractionPotential::WW => { doca_tungsten_tungsten_cubic_spline(r, impact_parameter, relative_energy) }, @@ -139,12 +163,19 @@ pub fn scaling_function(r: f64, a: f64, interaction_potential: InteractionPotent let n = 6.; 1./(1. + (r/sigma).powf(n)) }, + InteractionPotential::FOUR_EIGHT{alpha, beta} => { + let n = 8.; + 1./(1. + r.powi(8)/beta) + } InteractionPotential::MORSE{D, alpha, r0} => { 1./(1. + (r*alpha).powi(2)) } InteractionPotential::WW => { 1. }, + InteractionPotential::KRC_MORSE{D, alpha, r0, k, x0} => { + 1./(1. + (r*alpha).powi(2)) + } InteractionPotential::COULOMB{..} => panic!("Coulombic potential cannot be used with rootfinder.") } } @@ -239,35 +270,67 @@ pub fn screening_length(Za: f64, Zb: f64, interaction_potential: InteractionPote InteractionPotential::LENNARD_JONES_12_6{..} | InteractionPotential::LENNARD_JONES_65_6{..} => 0.8853*A0*(Za.sqrt() + Zb.sqrt()).powf(-2./3.), InteractionPotential::MORSE{D, alpha, r0} => alpha, InteractionPotential::COULOMB{Za: Z1, Zb: Z2} => 0.88534*A0/(Z1.powf(0.23) + Z2.powf(0.23)), + InteractionPotential::KRC_MORSE{..} => 0.8853*A0*(Za.sqrt() + Zb.sqrt()).powf(-2./3.), + InteractionPotential::FOUR_EIGHT{..} =>0.8853*A0*(Za.sqrt() + Zb.sqrt()).powf(-2./3.), } } /// Coefficients of inverse-polynomial interaction potentials. +/// Note: these go in order of decreasing orders of r, e.g., [r^n, r^n-1, r^n-2, ... r, 1]. pub fn polynomial_coefficients(relative_energy: f64, impact_parameter: f64, interaction_potential: InteractionPotential) -> Vec { match interaction_potential { InteractionPotential::LENNARD_JONES_12_6{sigma, epsilon} => { - vec![1., 0., -impact_parameter.powi(2), 0., 0., 0., 4.*epsilon*sigma.powf(6.)/relative_energy, 0., 0., 0., 0., 0., -4.*epsilon*sigma.powf(12.)/relative_energy] + let impact_parameter_angstroms = impact_parameter/ANGSTROM; + let epsilon_ev = epsilon/EV; + let sigma_angstroms = sigma/ANGSTROM; + let relative_energy_ev = relative_energy/EV; + //vec![1., 0., -impact_parameter.powi(2), 0., 0., 0., 4.*epsilon_ev*sigma.powf(6.)/relative_energy_ev, 0., 0., 0., 0., 0., -4.*epsilon_ev*sigma.powf(12.)/relative_energy_ev] + vec![1.0, -impact_parameter.powi(2), 0.0, 4.*epsilon_ev*sigma.powf(6.)/relative_energy_ev, 0.0, 0.0, -4.*epsilon_ev*sigma.powf(12.)/relative_energy_ev] }, InteractionPotential::LENNARD_JONES_65_6{sigma, epsilon} => { - vec![1., 0., 0., 0., -impact_parameter.powi(2), 0., 0., 0., 0., 0., 0., 0., 4.*epsilon*sigma.powf(6.)/relative_energy, -4.*epsilon*sigma.powf(6.5)/relative_energy] + let impact_parameter_angstroms = impact_parameter/ANGSTROM; + let epsilon_ev = epsilon/EV; + let sigma_angstroms = sigma/ANGSTROM; + let relative_energy_ev = relative_energy/EV; + vec![1., 0., 0., 0., -impact_parameter.powi(2), 0., 0., 0., 0., 0., 0., 0., 4.*epsilon_ev*sigma.powf(6.)/relative_energy_ev, -4.*epsilon_ev*sigma.powf(6.5)/relative_energy_ev] }, + InteractionPotential::FOUR_EIGHT{alpha, beta} => { + //Note: I've transformed to angstroms here to help the rootfinder with numerical issues. + //The units on alpha [m^4 J] and beta [m^8 J] make them unreasonably small to use numerically. + let alpha_angstroms4_joule = alpha/ANGSTROM.powi(4); + let beta_angstroms8_joule = beta/ANGSTROM.powi(8); + let impact_parameteter_angstroms = impact_parameter/ANGSTROM; + vec![1., -(impact_parameteter_angstroms).powi(2), alpha_angstroms4_joule/relative_energy, 0., -beta_angstroms8_joule/relative_energy,] + } _ => panic!("Input error: non-polynomial interaction potential used with polynomial root-finder.") } } -/// Inverse-polynomial interaction potentials transformation to remove singularities for the CPR root-finder. +/// Inverse-polynomial interaction potentials transformation to reduce numercial issues for the polynomial rootfinder. +/// This is for if, for example, you have non-integer powers and need to multiply by a factor to get to integer powers. +/// E.g., the 6.5-6 potential needs to be squared (x*x) giving you a 13th order polynomial in polynomial_coefficients(). pub fn inverse_transform(x: f64, interaction_potential: InteractionPotential) -> f64 { match interaction_potential { InteractionPotential::LENNARD_JONES_12_6{..} => { - x + x.sqrt() }, InteractionPotential::LENNARD_JONES_65_6{..} => { x*x }, + | InteractionPotential::FOUR_EIGHT{..} => { + x.sqrt()*ANGSTROM + } _ => panic!("Input error: non-polynomial interaction potential used with polynomial root-finder transformation.") } } +/// Four-Eight potential (Born repulsion) +pub fn four_eight(r: f64, alpha: f64, beta: f64) -> f64 { + let a = alpha.powf(1./4.); + let b = beta.powf(1./8.); + -(a/r).powi(4) + (b/r).powi(8) +} + /// Lennard-Jones 12-6 pub fn lennard_jones(r: f64, sigma: f64, epsilon: f64) -> f64 { 4.*epsilon*((sigma/r).powf(12.) - (sigma/r).powf(6.)) @@ -283,11 +346,29 @@ pub fn morse(r: f64, D: f64, alpha: f64, r0: f64) -> f64 { D*((-2.*alpha*(r - r0)).exp() - 2.*(-alpha*(r - r0)).exp()) } +/// Kr-C + Morse potential with sigmoid interpolation at crossing point of two potentials +pub fn krc_morse(r: f64, a: f64, Za: f64, Zb: f64, D: f64, alpha: f64, r0: f64, k: f64, x0: f64) -> f64 { + sigmoid(r, k, x0)*morse(r, D, alpha, r0) + sigmoid(r, -k, x0)*screened_coulomb(r, a, Za, Zb, InteractionPotential::KR_C) +} + +/// Distance of closest approach function for four-eight potential (Born repulsion) +pub fn doca_four_eight(r: f64, impact_parameter: f64, relative_energy: f64, alpha: f64, beta: f64) -> f64 { + let a = alpha.powf(1./4.); + let b = beta.powf(1./8.); + let b4 = beta.sqrt(); + (r/b).powf(8.) - (-(a*r/b/b) + 1.)/relative_energy - (impact_parameter*r.powf(3.)/b4) +} + /// Distance of closest approach function for Morse potential. pub fn doca_morse(r: f64, impact_parameter: f64, relative_energy: f64, D: f64, alpha: f64, r0: f64) -> f64 { (r*alpha).powi(2) - (r*alpha).powi(2)*D/relative_energy*((-2.*alpha*(r - r0)).exp() - 2.*(-alpha*(r - r0)).exp()) - (impact_parameter*alpha).powi(2) } +/// 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) +} + /// First derivative w.r.t. `r` of the distance of closest approach function for Morse potential. pub fn diff_doca_morse(r: f64, impact_parameter: f64, relative_energy: f64, D: f64, alpha: f64, r0: f64) -> f64 { 2.*alpha.powi(2)*r - 2.*alpha.powi(2)*D*r*(-2.*alpha*(r - r0) - 1.).exp()*(alpha*r*(alpha*(r - r0)).exp() - 2.*(alpha*(r - r0)).exp() - r*alpha + 1.) @@ -433,8 +514,9 @@ pub fn first_screening_radius(interaction_potential: InteractionPotential) -> f6 InteractionPotential::TRIDYN => 0.278544, InteractionPotential::LENNARD_JONES_12_6{..} => 1., InteractionPotential::LENNARD_JONES_65_6{..} => 1., - InteractionPotential::MORSE{..} => 1., + InteractionPotential::MORSE{..} | InteractionPotential::KRC_MORSE{..} => 1., InteractionPotential::WW => 1., InteractionPotential::COULOMB{..} => 0.3, + InteractionPotential::FOUR_EIGHT{..} => 1., } } diff --git a/src/lib.rs b/src/lib.rs index f851c31..05817c3 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -2,16 +2,12 @@ #![allow(non_snake_case)] #![allow(non_camel_case_types)] -#[cfg(feature = "cpr_rootfinder_openblas")] -extern crate openblas_src; -#[cfg(feature = "cpr_rootfinder_netlib")] -extern crate netlib_src; -#[cfg(feature = "cpr_rootfinder_intel_mkl")] -extern crate intel_mkl_src; - use std::{fmt}; use std::mem::discriminant; +use std::alloc::{dealloc, Layout}; +use std::mem::align_of; + //Parallelization - currently only used in python library functions #[cfg(feature = "python")] use rayon::prelude::*; @@ -77,7 +73,7 @@ pub mod parry; pub use crate::enums::*; pub use crate::consts::*; pub use crate::structs::*; -pub use crate::input::{Input2D, Input1D, Input0D, Options, InputFile, GeometryInput}; +pub use crate::input::{Input2D, InputHomogeneous2D, Input1D, Input0D, Options, InputFile, GeometryInput}; pub use crate::output::{OutputUnits}; pub use crate::geometry::{Geometry, GeometryElement, Mesh0D, Mesh1D, Mesh2D}; pub use crate::sphere::{Sphere, SphereInput, InputSphere}; @@ -96,18 +92,15 @@ pub fn libRustBCA(py: Python, m: &PyModule) -> PyResult<()> { m.add_function(wrap_pyfunction!(compound_bca_list_1D_py, m)?)?; m.add_function(wrap_pyfunction!(sputtering_yield, m)?)?; m.add_function(wrap_pyfunction!(reflection_coefficient, m)?)?; + m.add_function(wrap_pyfunction!(compound_reflection_coefficient, m)?)?; + m.add_function(wrap_pyfunction!(reflect_single_ion_py, m)?)?; #[cfg(feature = "parry3d")] m.add_function(wrap_pyfunction!(rotate_given_surface_normal_py, m)?)?; + #[cfg(feature = "parry3d")] m.add_function(wrap_pyfunction!(rotate_back_py, m)?)?; Ok(()) } -#[repr(C)] -pub struct OutputBCA { - pub len: usize, - pub particles: *mut [f64; 9], -} - #[derive(Debug)] #[repr(C)] pub struct InputSimpleBCA { @@ -168,6 +161,12 @@ pub struct InputTaggedBCA { pub weights: *mut f64, } +#[repr(C)] +pub struct OutputBCA { + pub len: usize, + pub particles: *mut [f64; 9], +} + #[derive(Debug)] #[repr(C)] pub struct OutputTaggedBCA { @@ -178,6 +177,39 @@ pub struct OutputTaggedBCA { pub incident: *mut bool, } +#[no_mangle] +pub extern "C" fn drop_output_tagged_bca(output: OutputTaggedBCA) { + let length = output.len; + + if length > 0 { + + let particles_layout = Layout::from_size_align(length, align_of::<[f64; 9]>()).unwrap(); + let weights_layout = Layout::from_size_align(length, align_of::()).unwrap(); + let tags_layout = Layout::from_size_align(length, align_of::()).unwrap(); + let incident_layout = Layout::from_size_align(length, align_of::()).unwrap(); + + unsafe { + dealloc(output.particles as *mut u8, particles_layout); + dealloc(output.weights as *mut u8, weights_layout); + dealloc(output.tags as *mut u8, tags_layout); + dealloc(output.incident as *mut u8, incident_layout); + }; + } +} + +#[no_mangle] +pub extern "C" fn drop_output_bca(output: OutputBCA) { + let length = output.len; + + if length > 0 { + let particles_layout = Layout::from_size_align(length, align_of::<[f64; 9]>()).unwrap(); + + unsafe { + dealloc(output.particles as *mut u8, particles_layout); + }; + } +} + #[no_mangle] pub extern "C" fn compound_tagged_bca_list_c(input: InputTaggedBCA) -> OutputTaggedBCA { @@ -208,6 +240,7 @@ pub extern "C" fn compound_tagged_bca_list_c(input: InputTaggedBCA) -> OutputTag Eb: Eb2, Es: Es2, Ec: Ec2, + Ed: vec![0.0; input.num_species_target], Z: Z2, m: m2, interaction_index: vec![0; input.num_species_target], @@ -245,6 +278,7 @@ pub extern "C" fn compound_tagged_bca_list_c(input: InputTaggedBCA) -> OutputTag E: E1, Ec: input.Ec1*EV, Es: input.Es1*EV, + Ed: 0.0, pos: Vector::new(x, y, z), dir: Vector::new(ux, uy, uz), pos_origin: Vector::new(x, y, z), @@ -302,6 +336,9 @@ pub extern "C" fn compound_tagged_bca_list_c(input: InputTaggedBCA) -> OutputTag let incident_ptr = output_incident.as_mut_ptr(); std::mem::forget(total_output); + std::mem::forget(output_tags); + std::mem::forget(output_weights); + std::mem::forget(output_incident); OutputTaggedBCA { len, @@ -312,7 +349,6 @@ pub extern "C" fn compound_tagged_bca_list_c(input: InputTaggedBCA) -> OutputTag } } - #[no_mangle] pub extern "C" fn reflect_single_ion_c(num_species_target: &mut c_int, ux: &mut f64, uy: &mut f64, uz: &mut f64, E1: &mut f64, Z1: &mut f64, m1: &mut f64, Ec1: &mut f64, Es1: &mut f64, Z2: *mut f64, m2: *mut f64, Ec2: *mut f64, Es2: *mut f64, Eb2: *mut f64, n2: *mut f64) { @@ -337,6 +373,7 @@ pub extern "C" fn reflect_single_ion_c(num_species_target: &mut c_int, ux: &mut Eb: Eb2, Es: Es2, Ec: Ec2, + Ed: vec![0.0; *num_species_target as usize], Z: Z2, m: m2, interaction_index: vec![0; *num_species_target as usize], @@ -358,6 +395,7 @@ pub extern "C" fn reflect_single_ion_c(num_species_target: &mut c_int, ux: &mut E: *E1*EV, Ec: *Ec1*EV, Es: *Es1*EV, + Ed: 0.0, pos: Vector::new(x, y, z), dir: Vector::new(*ux, *uy, *uz), pos_origin: Vector::new(x, y, z), @@ -409,6 +447,7 @@ pub extern "C" fn simple_bca_list_c(input: InputSimpleBCA) -> OutputBCA { Eb: vec![input.Eb2], Es: vec![input.Es2], Ec: vec![input.Ec2], + Ed: vec![0.0; input.len], Z: vec![input.Z2], m: vec![input.m2], interaction_index: vec![0], @@ -446,6 +485,7 @@ pub extern "C" fn simple_bca_list_c(input: InputSimpleBCA) -> OutputBCA { E: E1, Ec: input.Ec1*EV, Es: input.Es1*EV, + Ed: 0.0, pos: Vector::new(x, y, z), dir: Vector::new(ux, uy, uz), pos_origin: Vector::new(x, y, z), @@ -525,6 +565,7 @@ pub extern "C" fn compound_bca_list_c(input: InputCompoundBCA) -> OutputBCA { Eb: Eb2, Es: Es2, Ec: Ec2, + Ed: vec![0.0; input.num_species_target], Z: Z2, m: m2, interaction_index: vec![0; input.num_species_target], @@ -562,6 +603,7 @@ pub extern "C" fn compound_bca_list_c(input: InputCompoundBCA) -> OutputBCA { E: E1, Ec: input.Ec1*EV, Es: input.Es1*EV, + Ed: 0.0, pos: Vector::new(x, y, z), dir: Vector::new(ux, uy, uz), pos_origin: Vector::new(x, y, z), @@ -617,7 +659,6 @@ pub extern "C" fn compound_bca_list_c(input: InputCompoundBCA) -> OutputBCA { } } - #[no_mangle] pub extern "C" fn compound_bca_list_fortran(num_incident_ions: &mut c_int, track_recoils: &mut bool, ux: *mut f64, uy: *mut f64, uz: *mut f64, E1: *mut f64, @@ -663,6 +704,7 @@ pub extern "C" fn compound_bca_list_fortran(num_incident_ions: &mut c_int, track Eb: Eb2, Es: Es2, Ec: Ec2, + Ed: vec![0.0; *num_species_target as usize], Z: Z2, m: m2, interaction_index: vec![0; *num_species_target as usize], @@ -686,6 +728,7 @@ pub extern "C" fn compound_bca_list_fortran(num_incident_ions: &mut c_int, track E: *E1_*EV, Ec: Ec1_*EV, Es: Es1_*EV, + Ed: 0.0, pos: Vector::new(x, y, z), dir: Vector::new(ux_, uy_, uz_), pos_origin: Vector::new(x, y, z), @@ -806,6 +849,7 @@ pub fn compound_bca_list_py(energies: Vec, ux: Vec, uy: Vec, uz: Eb: Eb2, Es: Es2, Ec: Ec2, + Ed: vec![0.0; num_species_target], Z: Z2, m: m2, interaction_index: vec![0; num_species_target], @@ -870,6 +914,96 @@ pub fn compound_bca_list_py(energies: Vec, ux: Vec, uy: Vec, uz: (total_output, incident) } +#[cfg(feature = "python")] +///reflect_single_ion_py(ion, target, vx, vy, vz) +///Performs a single BCA ion trajectory in target material with specified incident velocity. +///Args: +/// ion (dict): dictionary that defines ion parameters; examples can be found in scripts/materials.py. +/// target (dict): dictionary that defines target parameterrs; examples can be found in scripts/materials.py. +/// vx, vy, vz (float): initial x, y, and z velocity in m/s. +///Returns: +/// vx, vy, vz (float): final x, y, and z velocity in m/s. When ion implants in material, vx, vy, and vz will all be zero. +#[pyfunction] +pub fn reflect_single_ion_py(ion: &PyDict, target: &PyDict, vx: f64, vy: f64, vz: f64) -> (f64, f64, f64){ + + let Z1 = unpack(ion.get_item("Z").expect("Cannot get ion Z from dictionary. Ensure ion['Z'] exists.")); + let m1 = unpack(ion.get_item("m").expect("Cannot get ion mass from dictionary. Ensure ion['m'] exists.")); + let Ec1 = unpack(ion.get_item("Ec").expect("Cannot get ion cutoff energy from dictionary. Ensure ion['Ec'] exists.")); + let Es1 = unpack(ion.get_item("Es").expect("Cannot get ion surface binding energy from dictionary. Ensure ion['Es'] exists.")); + + let Z2 = unpack(target.get_item("Z").expect("Cannot get target Z from dictionary. Ensure target['Z'] exists.")); + let m2 = unpack(target.get_item("m").expect("Cannot get target mass from dictionary. Ensure target['m'] exists.")); + let Ec2 = unpack(target.get_item("Ec").expect("Cannot get target cutoff energy from dictionary. Ensure target['Ec'] exists.")); + let Es2 = unpack(target.get_item("Es").expect("Cannot get target surface binding energy from dictionary. Ensure target['Es'] exists.")); + let Eb2 = unpack(target.get_item("Eb").expect("Cannot get target bulk binding energy from dictionary. Ensure target['Eb'] exists.")); + let n2 = unpack(target.get_item("n").expect("Cannot get target density from dictionary. Ensure target['n'] exists.")); + + assert!(vx > 0.0, "Input error: vx must be greater than zero for incident particles to hit surface at x=0."); + + let options = Options::default_options(false); + + let velocity2 = vx.powf(2.) + vy.powf(2.) + vz.powf(2.); //m/s + let energy_eV = 0.5*m1*AMU*velocity2/EV; //EV + + let ux = vx/velocity2.sqrt(); + let uy = vy/velocity2.sqrt(); + let uz = vz/velocity2.sqrt(); + + let material_parameters = material::MaterialParameters { + energy_unit: "EV".to_string(), + mass_unit: "AMU".to_string(), + Eb: vec![Eb2], + Es: vec![Es2], + Ec: vec![Ec2], + Ed: vec![0.0], + Z: vec![Z2], + m: vec![m2], + interaction_index: vec![0], + surface_binding_model: SurfaceBindingModel::AVERAGE, + bulk_binding_model: BulkBindingModel::INDIVIDUAL, + }; + + let geometry_input = geometry::Mesh0DInput { + length_unit: "M".to_string(), + densities: vec![n2], + electronic_stopping_correction_factor: 1.0 + }; + + let m = material::Material::::new(&material_parameters, &geometry_input); + + let x = -m.geometry.energy_barrier_thickness; + let y = 0.0; + let z = 0.0; + + let p = particle::Particle::default_incident( + m1, + Z1, + energy_eV, + Ec1, + Es1, + x, + ux, + uy, + uz + ); + + let output = bca::single_ion_bca(p, &m, &options); + + let reflected_energy = output[0].E; //Joules + + let reflected_velocity = (2.*reflected_energy/(m1*AMU)).sqrt(); //m/s + + let vx2 = output[0].dir.x*reflected_velocity; + let vy2 = output[0].dir.y*reflected_velocity; + let vz2 = output[0].dir.z*reflected_velocity; + + if output[0].E > 0.0 && output[0].dir.x < 0.0 && output[0].left && output[0].incident { + (vx2, vy2, vz2) + } else { + (0.0, 0.0, 0.0) + } +} + #[cfg(feature = "python")] ///compound_bca_list_1D_py(ux, uy, uz, energies, Z1, m1, Ec1, Es1, Z2, m2, Ec2, Es2, Eb2 n2, dx) /// runs a BCA simulation for a list of particles and outputs a list of sputtered, reflected, and implanted particles. @@ -931,6 +1065,7 @@ pub fn compound_bca_list_1D_py(ux: Vec, uy: Vec, uz: Vec, energie Eb: Eb2, Es: Es2, Ec: Ec2, + Ed: vec![0.0; num_species], Z: Z2, m: m2, interaction_index: vec![0; num_species], @@ -1091,6 +1226,7 @@ pub fn simple_bca(x: f64, y: f64, z: f64, ux: f64, uy: f64, uz: f64, E1: f64, Z1 E: E1*EV, Ec: Ec1*EV, Es: Es1*EV, + Ed: 0.0, pos: Vector::new(x, y, z), dir: Vector::new(ux, uy, uz), pos_origin: Vector::new(x, y, z), @@ -1119,6 +1255,7 @@ pub fn simple_bca(x: f64, y: f64, z: f64, ux: f64, uy: f64, uz: f64, E1: f64, Z1 Eb: vec![Eb2], Es: vec![Es2], Ec: vec![Ec2], + Ed: vec![0.0], Z: vec![Z2], m: vec![m2], interaction_index: vec![0], @@ -1165,6 +1302,7 @@ pub fn simple_compound_bca(x: f64, y: f64, z: f64, ux: f64, uy: f64, uz: f64, E1 E: E1*EV, Ec: Ec1*EV, Es: Es1*EV, + Ed: 0.0, pos: Vector::new(x, y, z), dir: Vector::new(ux, uy, uz), pos_origin: Vector::new(x, y, z), @@ -1193,6 +1331,7 @@ pub fn simple_compound_bca(x: f64, y: f64, z: f64, ux: f64, uy: f64, uz: f64, E1 Eb: Eb2, Es: Es2, Ec: Ec2, + Ed: vec![0.0; Z2.len()], Z: Z2, m: m2, interaction_index: vec![0], @@ -1250,7 +1389,7 @@ pub extern "C" fn rotate_given_surface_normal(nx: f64, ny: f64, nz: f64, ux: &mu //around a non-x axis; y is chosen arbitrarily Rotation3::from_axis_angle(&Vector3::y_axis(), PI).into() }; - + let incident = rotation_matrix*direction; // ux must not be exactly 1.0 to avoid gimbal lock in RustBCA @@ -1275,7 +1414,7 @@ pub extern "C" fn rotate_given_surface_normal(nx: f64, ny: f64, nz: f64, ux: &mu /// rotate_given_surface_normal_py(nx, ny, nz, ux, uy, uz) /// -- /// -/// This function runs a 0D Binary Collision Approximation simulation for the given incident ions and material. +/// This function takes a particle direction and a normal vector and rotates from simulation to RustBCA coordinates. /// Args: /// nx (f64): surface normal in global frame x-component. /// ny (f64): surface normal in global frame y-component. @@ -1286,7 +1425,6 @@ pub extern "C" fn rotate_given_surface_normal(nx: f64, ny: f64, nz: f64, ux: &mu /// Returns: /// direction (f64, f64, f64): direction vector of particle in RustBCA coordinates. pub fn rotate_given_surface_normal_py(nx: f64, ny: f64, nz: f64, ux: f64, uy: f64, uz: f64) -> (f64, f64, f64) { - let mut ux = ux; let mut uy = uy; let mut uz = uz; @@ -1294,26 +1432,13 @@ pub fn rotate_given_surface_normal_py(nx: f64, ny: f64, nz: f64, ux: f64, uy: f6 (ux, uy, uz) } -#[cfg(all(feature = "python", feature = "parry3d"))] -#[pyfunction] -/// rotate_given_surface_normal_py(nx, ny, nz, ux, uy, uz) -/// -- -/// -/// This function runs a 0D Binary Collision Approximation simulation for the given incident ions and material. -/// Args: -/// nx (f64): surface normal in global frame x-component. -/// ny (f64): surface normal in global frame y-component. -/// nz (f64): surface normal in global frame z-component. -/// ux (f64): particle direction in RustBCA frame x-component. -/// uy (f64): particle direction in RustBCA frame normal y-component. -/// uz (f64): particle direction in RustBCA frame normal z-component. -/// Returns: -/// direction (f64, f64, f64): direction vector of particle in global coordinates. -pub fn rotate_back_py(nx: f64, ny: f64, nz: f64, ux: f64, uy: f64, uz: f64) -> (f64, f64, f64) { +#[cfg(feature = "parry3d")] +#[no_mangle] +pub extern "C" fn rotate_back(nx: f64, ny: f64, nz: f64, ux: &mut f64, uy: &mut f64, uz: &mut f64) { let RUSTBCA_DIRECTION: Vector3:: = Vector3::::new(1.0, 0.0, 0.0); let into_surface = Vector3::new(-nx, -ny, -nz); - let direction = Vector3::new(ux, uy, uz); + let direction = Vector3::new(*ux, *uy, *uz); //Rotation to local RustBCA coordinates from global //Here's how this works: a rotation matrix is found that maps the rustbca @@ -1334,7 +1459,32 @@ pub fn rotate_back_py(nx: f64, ny: f64, nz: f64, ux: f64, uy: f64, uz: f64) -> ( let u = rotation_matrix.transpose()*direction; - (u.x, u.y, u.z) + *ux = u.x; + *uy = u.y; + *uz = u.z; +} + +#[cfg(all(feature = "python", feature = "parry3d"))] +#[pyfunction] +/// rotate_back_py(nx, ny, nz, ux, uy, uz) +/// -- +/// +/// This function takes a particle direction and a normal vector and rotates from RustBCA to simulation coordinates. +/// Args: +/// nx (f64): surface normal in global frame x-component. +/// ny (f64): surface normal in global frame y-component. +/// nz (f64): surface normal in global frame z-component. +/// ux (f64): particle direction in RustBCA frame x-component. +/// uy (f64): particle direction in RustBCA frame normal y-component. +/// uz (f64): particle direction in RustBCA frame normal z-component. +/// Returns: +/// direction (f64, f64, f64): direction vector of particle in global coordinates. +pub fn rotate_back_py(nx: f64, ny: f64, nz: f64, ux: f64, uy: f64, uz: f64) -> (f64, f64, f64) { + let mut ux = ux; + let mut uy = uy; + let mut uz = uz; + rotate_back(nx, ny, nz, &mut ux, &mut uy, &mut uz); + (ux, uy, uz) } #[cfg(feature = "python")] @@ -1355,7 +1505,9 @@ fn unpack(python_float: &PyAny) -> f64 { /// num_samples: number of ion trajectories to run; precision will go as 1/sqrt(N) pub fn sputtering_yield(ion: &PyDict, target: &PyDict, energy: f64, angle: f64, num_samples: usize) -> f64 { - const DELTA: f64 = 1e-6; + assert!(angle.abs() <= 90.0, "Incident angle w.r.t. surface normal, {}, cannot exceed 90 degrees.", angle); + + const DELTA: f64 = 1e-6; let Z1 = unpack(ion.get_item("Z").expect("Cannot get ion Z from dictionary. Ensure ion['Z'] exists.")); let m1 = unpack(ion.get_item("m").expect("Cannot get ion mass from dictionary. Ensure ion['m'] exists.")); @@ -1374,8 +1526,8 @@ pub fn sputtering_yield(ion: &PyDict, target: &PyDict, energy: f64, angle: f64, let y = 0.0; let z = 0.0; - let ux = (angle/180.0*PI).cos() - DELTA; - let uy = (angle/180.0*PI).sin() + DELTA; + let ux = (angle/180.0*PI).cos() + DELTA; + let uy = (angle/180.0*PI).sin() - DELTA; let uz = 0.0; let material_parameters = material::MaterialParameters { @@ -1384,6 +1536,7 @@ pub fn sputtering_yield(ion: &PyDict, target: &PyDict, energy: f64, angle: f64, Eb: vec![Eb2], Es: vec![Es2], Ec: vec![Ec2], + Ed: vec![0.0], Z: vec![Z2], m: vec![m2], interaction_index: vec![0], @@ -1427,7 +1580,7 @@ pub fn sputtering_yield(ion: &PyDict, target: &PyDict, energy: f64, angle: f64, } }); let num_sputtered = *num_sputtered.lock().unwrap(); - num_sputtered as f64 / num_samples as f64 + num_sputtered as f64 / num_samples as f64 } #[cfg(feature = "python")] @@ -1445,7 +1598,9 @@ pub fn sputtering_yield(ion: &PyDict, target: &PyDict, energy: f64, angle: f64, /// R_E (f64): energy reflection coefficient (sum of reflected particle energies / total incident energy) pub fn reflection_coefficient(ion: &PyDict, target: &PyDict, energy: f64, angle: f64, num_samples: usize) -> (f64, f64) { - const DELTA: f64 = 1e-6; + const DELTA: f64 = 1e-6; + + assert!(angle.abs() <= 90.0, "Incident angle w.r.t. surface normal, {}, cannot exceed 90 degrees.", angle); let Z1 = unpack(ion.get_item("Z").expect("Cannot get ion Z from dictionary. Ensure ion['Z'] exists.")); let m1 = unpack(ion.get_item("m").expect("Cannot get ion mass from dictionary. Ensure ion['m'] exists.")); @@ -1464,8 +1619,8 @@ pub fn reflection_coefficient(ion: &PyDict, target: &PyDict, energy: f64, angle: let y = 0.0; let z = 0.0; - let ux = (angle/180.0*PI).cos() - DELTA; - let uy = (angle/180.0*PI).sin() + DELTA; + let ux = (angle/180.0*PI).cos() + DELTA; + let uy = (angle/180.0*PI).sin() - DELTA; let uz = 0.0; let mut direction = Vector::new(ux, uy, uz); @@ -1477,6 +1632,7 @@ pub fn reflection_coefficient(ion: &PyDict, target: &PyDict, energy: f64, angle: Eb: vec![Eb2], Es: vec![Es2], Ec: vec![Ec2], + Ed: vec![0.0], Z: vec![Z2], m: vec![m2], interaction_index: vec![0], @@ -1510,7 +1666,110 @@ pub fn reflection_coefficient(ion: &PyDict, target: &PyDict, energy: f64, angle: uy, uz ); - + + let output = bca::single_ion_bca(p, &m, &options); + + for particle in output { + if particle.E > 0.0 && particle.dir.x < 0.0 && particle.left && particle.incident { + let mut num_reflected = num_reflected.lock().unwrap(); + *num_reflected += 1; + let mut energy_reflected = energy_reflected.lock().unwrap(); + *energy_reflected += particle.E; + } + } + }); + let num_reflected = *num_reflected.lock().unwrap(); + let energy_reflected = *energy_reflected.lock().unwrap(); + + (num_reflected as f64 / num_samples as f64, energy_reflected / EV / energy / num_samples as f64) +} + +#[cfg(feature = "python")] +#[pyfunction] +/// compound_reflection_coefficient(ion, target_species, target_number_densities, energy, angle, num_samples) +/// A routine the calculates the reflection coefficient of energetic ions incident upon materials using RustBCA. +/// Args: +/// ion: a dictionary with the keys Z (atomic number), m (atomic mass in AMU), Ec (cutoff energy in eV), Es (surface binding energy in eV) +/// target_species: a list of dictionaries with the keys Z, m, Ec, Es, Eb (bulk binding energy in eV), n2 (number density in 1/m3) +/// target_number_densities (list(f64)): number density of each target species in the compound +/// energy: the incident energy of the ion in eV +/// angle: incident angle of the ion in degrees from surface normal +/// num_samples: number of ion trajectories to run; precision will go as 1/sqrt(N) +/// Returns: +/// R_N (f64): reflection coefficient (number of particles reflected / number of incident particles) +/// R_E (f64): energy reflection coefficient (sum of reflected particle energies / total incident energy) +pub fn compound_reflection_coefficient(ion: &PyDict, targets: Vec<&PyDict>, target_number_densities: Vec, energy: f64, angle: f64, num_samples: usize) -> (f64, f64) { + + const DELTA: f64 = 1e-6; + + assert!(angle.abs() <= 90.0, "Incident angle w.r.t. surface normal, {}, cannot exceed 90 degrees.", angle); + + let Z1 = unpack(ion.get_item("Z").expect("Cannot get ion Z from dictionary. Ensure ion['Z'] exists.")); + let m1 = unpack(ion.get_item("m").expect("Cannot get ion mass from dictionary. Ensure ion['m'] exists.")); + let Ec1 = unpack(ion.get_item("Ec").expect("Cannot get ion cutoff energy from dictionary. Ensure ion['Ec'] exists.")); + let Es1 = unpack(ion.get_item("Es").expect("Cannot get ion surface binding energy from dictionary. Ensure ion['Es'] exists.")); + + let Z2: Vec = targets.iter().enumerate().map( |(index, item)| unpack(item.get_item("Z").unwrap_or_else(|| panic!("Cannot get target Z from dictionary at index {}. Ensure target['Z'] exists.", index)))).collect(); + let m2: Vec = targets.iter().enumerate().map( |(index, item)| unpack(item.get_item("m").unwrap_or_else(|| panic!("Cannot get target m from dictionary at index {}. Ensure target['m'] exists.", index)))).collect(); + let Ec2: Vec = targets.iter().enumerate().map( |(index, item)| unpack(item.get_item("Ec").unwrap_or_else(|| panic!("Cannot get target Ec from dictionary at index {}. Ensure target['Ec'] exists.", index)))).collect(); + let Es2: Vec = targets.iter().enumerate().map( |(index, item)| unpack(item.get_item("Es").unwrap_or_else(|| panic!("Cannot get target Es from dictionary at index {}. Ensure target['Es'] exists.", index)))).collect(); + let Eb2: Vec = targets.iter().enumerate().map( |(index, item)| unpack(item.get_item("Eb").unwrap_or_else(|| panic!("Cannot get target Eb from dictionary at index {}. Ensure target['Eb'] exists.", index)))).collect(); + + let number_target_species = Z2.len(); + + let options = Options::default_options(false); + + let y = 0.0; + let z = 0.0; + + let ux = (angle/180.0*PI).cos() + DELTA; + let uy = (angle/180.0*PI).sin() - DELTA; + let uz = 0.0; + + let mut direction = Vector::new(ux, uy, uz); + direction.normalize(); + + let material_parameters = material::MaterialParameters { + energy_unit: "EV".to_string(), + mass_unit: "AMU".to_string(), + Eb: Eb2, + Es: Es2, + Ec: Ec2, + Ed: vec![0.0; number_target_species], + Z: Z2, + m: m2, + interaction_index: vec![0; number_target_species], + surface_binding_model: SurfaceBindingModel::AVERAGE, + bulk_binding_model: BulkBindingModel::INDIVIDUAL, + }; + + let geometry_input = geometry::Mesh0DInput { + length_unit: "M".to_string(), + densities: target_number_densities, + electronic_stopping_correction_factor: 1.0 + }; + + let m = material::Material::::new(&material_parameters, &geometry_input); + + let x = -m.geometry.energy_barrier_thickness; + + let num_reflected = Mutex::new(0); + let energy_reflected = Mutex::new(0.0); + + (0..num_samples as u64).into_par_iter().for_each( |index| { + + let p = particle::Particle::default_incident( + m1, + Z1, + energy, + Ec1, + Es1, + x, + ux, + uy, + uz + ); + let output = bca::single_ion_bca(p, &m, &options); for particle in output { @@ -1526,4 +1785,4 @@ pub fn reflection_coefficient(ion: &PyDict, target: &PyDict, energy: f64, angle: let energy_reflected = *energy_reflected.lock().unwrap(); (num_reflected as f64 / num_samples as f64, energy_reflected / EV / energy / num_samples as f64) -} \ No newline at end of file +} diff --git a/src/main.rs b/src/main.rs index 49b2a10..7da24c6 100644 --- a/src/main.rs +++ b/src/main.rs @@ -2,13 +2,6 @@ #![allow(non_snake_case)] #![allow(non_camel_case_types)] -#[cfg(feature = "cpr_rootfinder_openblas")] -extern crate openblas_src; -#[cfg(feature = "cpr_rootfinder_netlib")] -extern crate netlib_src; -#[cfg(feature = "cpr_rootfinder_intel_mkl")] -extern crate intel_mkl_src; - use std::{env, fmt}; use std::mem::discriminant; @@ -63,9 +56,9 @@ pub mod parry; pub use crate::enums::*; pub use crate::consts::*; pub use crate::structs::*; -pub use crate::input::{Input2D, Input1D, Input0D, Options, InputFile, GeometryInput}; +pub use crate::input::{Input2D, InputHomogeneous2D, Input1D, Input0D, Options, InputFile, GeometryInput}; pub use crate::output::{OutputUnits}; -pub use crate::geometry::{Geometry, GeometryElement, Mesh0D, Mesh1D, Mesh2D}; +pub use crate::geometry::{Geometry, GeometryElement, Mesh0D, Mesh1D, Mesh2D, HomogeneousMesh2D}; pub use crate::sphere::{Sphere, SphereInput, InputSphere}; pub use crate::physics::{physics_loop}; @@ -88,6 +81,7 @@ fn main() { "BALL" => GeometryType::BALL, #[cfg(feature = "parry3d")] "TRIMESH" => GeometryType::TRIMESH, + "HOMOGENEOUS2D" => GeometryType::HOMOGENEOUS2D, _ => panic!("Unimplemented geometry {}.", args[1].clone()) }), _ => panic!("Too many command line arguments. RustBCA accepts 0 (use 'input.toml') 1 () or 2 ( )"), @@ -120,5 +114,9 @@ fn main() { let (particle_input_array, material, options, output_units) = input::input::(input_file); physics_loop::(particle_input_array, material, options, output_units); } + GeometryType::HOMOGENEOUS2D => { + let (particle_input_array, material, options, output_units) = input::input::(input_file); + physics_loop::(particle_input_array, material, options, output_units); + } } } diff --git a/src/material.rs b/src/material.rs index 9b3eb86..e6931fa 100644 --- a/src/material.rs +++ b/src/material.rs @@ -13,6 +13,10 @@ fn default_vec_zero() -> Vec { vec![0] } +fn default_vec_float_zero() -> Vec { + vec![0.0] +} + /// Holds material input parameters from [material_params]. #[derive(Deserialize, Clone)] pub struct MaterialParameters { @@ -21,6 +25,8 @@ pub struct MaterialParameters { pub Eb: Vec, pub Es: Vec, pub Ec: Vec, + #[serde(default = "default_vec_float_zero")] + pub Ed: Vec, pub Z: Vec, pub m: Vec, #[serde(default = "default_vec_zero")] @@ -38,6 +44,7 @@ pub struct Material { pub Eb: Vec, pub Es: Vec, pub Ec: Vec, + pub Ed: Vec, pub interaction_index: Vec, pub geometry: Box, pub surface_binding_model: SurfaceBindingModel, @@ -73,6 +80,7 @@ impl Material { Eb: material_parameters.Eb.iter().map(|&i| i*energy_unit).collect(), Es: material_parameters.Es.iter().map(|&i| i*energy_unit).collect(), Ec: material_parameters.Ec.iter().map(|&i| i*energy_unit).collect(), + Ed: material_parameters.Ed.iter().map(|&i| i*energy_unit).collect(), interaction_index: material_parameters.interaction_index.clone(), surface_binding_model: material_parameters.surface_binding_model, bulk_binding_model: material_parameters.bulk_binding_model, @@ -235,13 +243,13 @@ impl Material { } ///Choose the parameters of a target atom as a concentration-weighted random draw from the species in the triangle that contains or is nearest to (x, y). - pub fn choose(&self, x: f64, y: f64, z: f64) -> (usize, f64, f64, f64, f64, usize) { + pub fn choose(&self, x: f64, y: f64, z: f64) -> (usize, f64, f64, f64, f64, f64, usize) { let random_number: f64 = rand::random::(); let cumulative_concentrations = self.get_cumulative_concentrations(x, y, z); for (component_index, cumulative_concentration) in cumulative_concentrations.iter().enumerate() { if random_number < *cumulative_concentration { - return (component_index, self.Z[component_index], self.m[component_index], self.Ec[component_index], self.Es[component_index], self.interaction_index[component_index]); + return (component_index, self.Z[component_index], self.m[component_index], self.Ec[component_index], self.Es[component_index], self.Ed[component_index], self.interaction_index[component_index]); } } panic!("Input error: method choose() operation failed to choose a valid species. Check densities."); @@ -395,12 +403,12 @@ pub fn boundary_condition_planar(particle_1: &mut particle::Particl if !material.inside_simulation_boundary(x, y, z) { particle_1.left = true; - particle_1.add_trajectory(); + particle_1.update_trajectory_tracker(); } if (E < Ec) & !particle_1.left & material.inside_energy_barrier(x, y, z) { particle_1.stopped = true; - particle_1.add_trajectory(); + particle_1.update_trajectory_tracker(); } if !particle_1.stopped & !particle_1.left { diff --git a/src/output.rs b/src/output.rs index da0c27c..b28401f 100644 --- a/src/output.rs +++ b/src/output.rs @@ -340,7 +340,7 @@ pub fn output_lists(output_list_streams: &mut OutputListStreams, particle: parti let energy_unit = output_units.energy_unit; let mass_unit = output_units.mass_unit; - if !particle.incident & options.track_displacements { + if !particle.incident & options.track_displacements & (particle.energy_origin > particle.Ed) { writeln!( output_list_streams.displacements_file_stream, "{},{},{},{},{},{},{},{},{}", particle.m/mass_unit, particle.Z, particle.energy_origin/energy_unit, diff --git a/src/particle.rs b/src/particle.rs index 4951b0d..9de3daa 100644 --- a/src/particle.rs +++ b/src/particle.rs @@ -72,6 +72,7 @@ pub struct Particle { pub E: f64, pub Ec: f64, pub Es: f64, + pub Ed: f64, pub pos: Vector, pub dir: Vector, pub pos_old: Vector, @@ -111,6 +112,7 @@ impl Particle { E: input.E, Ec: input.Ec, Es: input.Es, + Ed: 0.0, pos: Vector::new(input.x, input.y, input.z), dir: Vector::new(dirx/dir_mag, diry/dir_mag, dirz/dir_mag), pos_old: Vector::new(input.x, input.y, input.z), @@ -135,7 +137,7 @@ impl Particle { } /// Particle constructor from raw inputs. - pub fn new(m: f64, Z: f64, E: f64, Ec: f64, Es: f64, x: f64, y: f64, z: f64, dirx: f64, diry: f64, dirz: f64, incident: bool, track_trajectories: bool, interaction_index: usize) -> Particle { + pub fn new(m: f64, Z: f64, E: f64, Ec: f64, Es: f64, Ed: f64, x: f64, y: f64, z: f64, dirx: f64, diry: f64, dirz: f64, incident: bool, track_trajectories: bool, interaction_index: usize) -> Particle { let dir_mag = (dirx*dirx + diry*diry + dirz*dirz).sqrt(); Particle { @@ -144,6 +146,7 @@ impl Particle { E, Ec, Es, + Ed, pos: Vector::new(x, y, z), dir: Vector::new(dirx/dir_mag, diry/dir_mag, dirz/dir_mag), pos_old: Vector::new(x, y, z), @@ -166,7 +169,7 @@ impl Particle { tracked_vector: Vector::new(0.0, 0.0, 0.0), } } - + /// Particle constructor from simplified input. pub fn default_incident(m_amu: f64, Z: f64, E_eV: f64, Ec_eV: f64, Es_eV: f64, x: f64, dirx: f64, diry: f64, dirz: f64) -> Particle { let dir_mag = (dirx*dirx + diry*diry + dirz*dirz).sqrt(); @@ -183,6 +186,7 @@ impl Particle { E: E_eV*EV, Ec: Ec_eV*EV, Es: Es_eV*EV, + Ed: 0.0, pos: Vector::new(x, y, z), dir: Vector::new(dirx/dir_mag, diry/dir_mag, dirz/dir_mag), pos_old: Vector::new(x, y, z), @@ -207,14 +211,14 @@ impl Particle { } /// If `track_trajectories`, add the current (E, x, y, z) to the trajectory. - pub fn add_trajectory(&mut self) { + pub fn update_trajectory_tracker(&mut self) { if self.track_trajectories { self.trajectory.push(TrajectoryElement {E: self.E, x: self.pos.x, y: self.pos.y, z: self.pos.z}); } } /// If `track_energy_losses`, add the most recent electronic and nuclear energy loss terms and (x, y, z) to the energy loss tracker. - pub fn energy_loss(&mut self, options: &Options, En: f64, Ee: f64) { + pub fn update_energy_loss_tracker(&mut self, options: &Options, En: f64, Ee: f64) { if self.incident & options.track_energy_losses { self.energies.push(EnergyLoss {Ee, En, x: self.pos.x, y: self.pos.y, z: self.pos.z}); } @@ -256,7 +260,7 @@ impl Particle { pub fn advance(&mut self, mfp: f64, asymptotic_deflection: f64) -> f64 { if self.E > self.Ec { - self.add_trajectory(); + self.update_trajectory_tracker(); } //Update previous position diff --git a/src/physics.rs b/src/physics.rs index 2278b17..9dec666 100644 --- a/src/physics.rs +++ b/src/physics.rs @@ -1,81 +1,81 @@ -use super::*; - -pub fn physics_loop(particle_input_array: Vec, material: material::Material, options: Options, output_units: OutputUnits) { - - println!("Processing {} ions...", particle_input_array.len()); - - let total_count: u64 = particle_input_array.len() as u64; - assert!(total_count/options.num_chunks > 0, "Input error: chunk size == 0 - reduce num_chunks or increase particle count."); - - #[cfg(not(feature = "no_list_output"))] - let mut output_list_streams = output::open_output_lists(&options); - - let mut summary = output::SummaryPerSpecies::new(&options); - - #[cfg(feature = "distributions")] - let mut distributions = output::Distributions::new(&options); - - //Initialize threads with rayon - println!("Initializing with {} threads...", options.num_threads); - if options.num_threads > 1 {let pool = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build_global().unwrap();}; - - //Create and configure progress bar - let bar: ProgressBar = ProgressBar::new(total_count); - bar.set_style(ProgressStyle::default_bar() - .template("[{elapsed_precise}][{bar:40.cyan/blue}][{eta_precise}] {percent}%") - .progress_chars("#>-")); - - //Main loop - for (chunk_index, particle_input_chunk) in particle_input_array.chunks((total_count/options.num_chunks) as usize).enumerate() { - - let mut finished_particles: Vec = Vec::new(); - - if options.num_threads > 1 { - // BCA loop is implemented as parallelized extension of a per-chunk initially empty - // finished particle array via map from particle -> finished particles via BCA - finished_particles.par_extend( - particle_input_chunk.into_par_iter() - .map(|particle_input| { - bar.tick(); - bar.inc(1); - bca::single_ion_bca(particle::Particle::from_input(*particle_input, &options), &material, &options) - }).flatten() - ); - } else { - finished_particles.extend( - particle_input_chunk.iter() - .map(|particle_input| { - bar.tick(); - bar.inc(1); - bca::single_ion_bca(particle::Particle::from_input(*particle_input, &options), &material, &options) - }).flatten() - ); - } - - // Process this chunk of finished particles for output - for particle in finished_particles { - - summary.update(&particle); - - #[cfg(feature = "distributions")] - distributions.update(&particle, &output_units, &options, total_count as usize); - - #[cfg(not(feature = "no_list_output"))] - output::output_lists(&mut output_list_streams, particle, &options, &output_units); - - } - //Flush all file streams before dropping to ensure all data is written - #[cfg(not(feature = "no_list_output"))] - output::output_list_flush(&mut output_list_streams); - } - - summary.print(&options, &output_units); - - //Write distributions to file - #[cfg(feature = "distributions")] - distributions.print(&options); - - bar.finish(); - println!("Finished!"); - -} +use super::*; + +pub fn physics_loop(particle_input_array: Vec, material: material::Material, options: Options, output_units: OutputUnits) { + + println!("Processing {} ions...", particle_input_array.len()); + + let total_count: u64 = particle_input_array.len() as u64; + assert!(total_count/options.num_chunks > 0, "Input error: chunk size == 0 - reduce num_chunks or increase particle count."); + + #[cfg(not(feature = "no_list_output"))] + let mut output_list_streams = output::open_output_lists(&options); + + let mut summary = output::SummaryPerSpecies::new(&options); + + #[cfg(feature = "distributions")] + let mut distributions = output::Distributions::new(&options); + + //Initialize threads with rayon + println!("Initializing with {} threads...", options.num_threads); + if options.num_threads > 1 {let pool = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build_global().unwrap();}; + + //Create and configure progress bar + let bar: ProgressBar = ProgressBar::new(total_count); + bar.set_style(ProgressStyle::default_bar() + .template("[{elapsed_precise}][{bar:40.cyan/blue}][{eta_precise}] {percent}%") + .progress_chars("#>-")); + + //Main loop + for (chunk_index, particle_input_chunk) in particle_input_array.chunks((total_count/options.num_chunks) as usize).enumerate() { + + let mut finished_particles: Vec = Vec::new(); + + if options.num_threads > 1 { + // BCA loop is implemented as parallelized extension of a per-chunk initially empty + // finished particle array via map from particle -> finished particles via BCA + finished_particles.par_extend( + particle_input_chunk.into_par_iter() + .map(|particle_input| { + bar.tick(); + bar.inc(1); + bca::single_ion_bca(particle::Particle::from_input(*particle_input, &options), &material, &options) + }).flatten() + ); + } else { + finished_particles.extend( + particle_input_chunk.iter() + .map(|particle_input| { + bar.tick(); + bar.inc(1); + bca::single_ion_bca(particle::Particle::from_input(*particle_input, &options), &material, &options) + }).flatten() + ); + } + + // Process this chunk of finished particles for output + for particle in finished_particles { + + summary.update(&particle); + + #[cfg(feature = "distributions")] + distributions.update(&particle, &output_units, &options, total_count as usize); + + #[cfg(not(feature = "no_list_output"))] + output::output_lists(&mut output_list_streams, particle, &options, &output_units); + + } + //Flush all file streams before dropping to ensure all data is written + #[cfg(not(feature = "no_list_output"))] + output::output_list_flush(&mut output_list_streams); + } + + summary.print(&options, &output_units); + + //Write distributions to file + #[cfg(feature = "distributions")] + distributions.print(&options); + + bar.finish(); + println!("Finished!"); + +} diff --git a/src/tests.rs b/src/tests.rs index d9805a7..fd263f0 100644 --- a/src/tests.rs +++ b/src/tests.rs @@ -3,6 +3,20 @@ use super::*; #[cfg(test)] use float_cmp::*; + +#[test] +#[cfg(feature = "cpr_rootfinder")] +fn test_polynom() { + use rcpr::chebyshev::*; + let interaction_potential = InteractionPotential::FOUR_EIGHT{alpha: 1., beta: 1.}; + let coefficients = interactions::polynomial_coefficients(1., 1., interaction_potential); + let roots = real_polynomial_roots(coefficients.clone(), 1e-9).unwrap(); + + let max_root = roots.iter().cloned().fold(f64::NAN, f64::max); + println!("{}", max_root); + let inverse_transformed_root = interactions::inverse_transform(max_root, interaction_potential); +} + #[test] #[cfg(feature = "parry3d")] fn test_parry_cuboid() { @@ -17,7 +31,7 @@ fn test_parry_cuboid() { let cosx = 1./(2.0_f64).sqrt(); let cosy = 1./(2.0_f64).sqrt(); let cosz = 0.; - let particle_1 = particle::Particle::new(mass, Z, E, Ec, Es, x, y, z, cosx, cosy, cosz, false, false, 0); + let particle_1 = particle::Particle::new(mass, Z, E, Ec, Es, 0.0, x, y, z, cosx, cosy, cosz, false, false, 0); let material_parameters = material::MaterialParameters{ energy_unit: "EV".to_string(), @@ -25,6 +39,7 @@ fn test_parry_cuboid() { Eb: vec![0.0, 0.0], Es: vec![2.0, 4.0], Ec: vec![1.0, 1.0], + Ed: vec![0.0, 0.0], Z: vec![29., 1.], m: vec![63.54, 1.0008], interaction_index: vec![0, 0], @@ -105,7 +120,7 @@ fn test_parry_sphere() { let cosx = 1./(2.0_f64).sqrt(); let cosy = 1./(2.0_f64).sqrt(); let cosz = 0.; - let mut particle_1 = particle::Particle::new(mass, Z, E, Ec, Es, x, y, z, cosx, cosy, cosz, false, false, 0); + let mut particle_1 = particle::Particle::new(mass, Z, E, Ec, Es, 0.0, x, y, z, cosx, cosy, cosz, false, false, 0); let material_parameters = material::MaterialParameters{ energy_unit: "EV".to_string(), @@ -113,6 +128,7 @@ fn test_parry_sphere() { Eb: vec![0.0, 0.0], Es: vec![2.0, 4.0], Ec: vec![1.0, 1.0], + Ed: vec![0.0, 0.0], Z: vec![29., 1.], m: vec![63.54, 1.0008], interaction_index: vec![0, 0], @@ -247,7 +263,7 @@ fn test_distributions() { let cosx = 0.0; let cosy = 0.0; let cosz = 0.0; - let mut particle = particle::Particle::new(mass, Z, E, Ec, Es, x, y, z, cosx, cosy, cosz, false, false, 0); + let mut particle = particle::Particle::new(mass, Z, E, Ec, Es, 0.0, x, y, z, cosx, cosy, cosz, false, false, 0); let mut distributions = output::Distributions::new(&options); assert_eq!(distributions.x_range[0], 0.); @@ -336,7 +352,7 @@ fn test_spherical_geometry() { let cosx = 1./(2.0_f64).sqrt(); let cosy = 1./(2.0_f64).sqrt(); let cosz = 0.; - let mut particle_1 = particle::Particle::new(mass, Z, E, Ec, Es, x, y, z, cosx, cosy, cosz, false, false, 0); + let mut particle_1 = particle::Particle::new(mass, Z, E, Ec, Es, 0.0, x, y, z, cosx, cosy, cosz, false, false, 0); let material_parameters = material::MaterialParameters{ energy_unit: "EV".to_string(), @@ -344,6 +360,7 @@ fn test_spherical_geometry() { Eb: vec![0.0, 0.0], Es: vec![2.0, 4.0], Ec: vec![1.0, 1.0], + Ed: vec![0.0, 0.0], Z: vec![29., 1.], m: vec![63.54, 1.0008], interaction_index: vec![0, 0], @@ -435,7 +452,7 @@ fn test_geometry() { let cosx = 1./(2.0_f64).sqrt(); let cosy = 1./(2.0_f64).sqrt(); let cosz = 0.; - let mut particle_1 = particle::Particle::new(mass, Z, E, Ec, Es, x, y, z, cosx, cosy, cosz, false, false, 0); + let mut particle_1 = particle::Particle::new(mass, Z, E, Ec, Es, 0.0, x, y, z, cosx, cosy, cosz, false, false, 0); let material_parameters = material::MaterialParameters{ energy_unit: "EV".to_string(), @@ -443,6 +460,7 @@ fn test_geometry() { Eb: vec![0.0, 0.0], Es: vec![2.0, 4.0], Ec: vec![1.0, 1.0], + Ed: vec![0.0, 0.0], Z: vec![29., 1.], m: vec![63.54, 1.0008], interaction_index: vec![0, 0], @@ -464,6 +482,22 @@ fn test_geometry() { electronic_stopping_correction_factors: vec![1.0, 1.0], }; + let geometry_input_homogeneous_2D = geometry::HomogeneousMesh2DInput { + length_unit: "ANGSTROM".to_string(), + points: vec![(0., -thickness/2.), (depth, -thickness/2.), (depth, thickness/2.), (0., thickness/2.)], + densities: vec![0.03, 0.03], + simulation_boundary_points: vec![(0., 1.1*thickness/2.), (depth, 1.1*thickness/2.), (depth, -1.1*thickness/2.), (0., -1.1*thickness/2.), (0., 1.1*thickness/2.)], + electronic_stopping_correction_factor: 1.0, + }; + + let geometry_input_1D = geometry::Mesh1DInput { + length_unit: "ANGSTROM".to_string(), + densities: vec![vec![0.03, 0.03], vec![0.03, 0.03]], + layer_thicknesses: vec![thickness/2., thickness/2.], + electronic_stopping_correction_factors: vec![1.0, 1.0, 1.0] + + }; + let geometry_input_0D = geometry::Mesh0DInput { length_unit: "ANGSTROM".to_string(), densities: vec![0.03, 0.03], @@ -471,6 +505,8 @@ fn test_geometry() { }; let material_2D: material::Material = material::Material::::new(&material_parameters, &geometry_input_2D); + let material_homogeneous_2D: material::Material = material::Material::::new(&material_parameters, &geometry_input_homogeneous_2D); + let material_1D: material::Material = material::Material::::new(&material_parameters, &geometry_input_1D); let mut material_0D: material::Material = material::Material::::new(&material_parameters, &geometry_input_0D); material_0D.geometry.energy_barrier_thickness = 10.*ANGSTROM; @@ -480,6 +516,15 @@ fn test_geometry() { particle_1.pos_old.x = -500.*ANGSTROM; particle_1.pos_old.y = 0.; + assert!(material_2D.inside(particle_1.pos.x, particle_1.pos.y, particle_1.pos.z)); + assert!(!material_2D.inside(particle_1.pos_old.x, particle_1.pos_old.y, particle_1.pos_old.z)); + + assert!(material_homogeneous_2D.inside(particle_1.pos.x, particle_1.pos.y, particle_1.pos.z)); + assert!(!material_homogeneous_2D.inside(particle_1.pos_old.x, particle_1.pos_old.y, particle_1.pos_old.z)); + + assert!(material_0D.inside(particle_1.pos.x, particle_1.pos.y, particle_1.pos.z)); + assert!(!material_0D.inside(particle_1.pos_old.x, particle_1.pos_old.y, particle_1.pos_old.z)); + assert_eq!(material_2D.geometry.get_ck(0., 0., 0.), material_0D.geometry.get_ck(0., 0., 0.)); assert_eq!(material_2D.geometry.get_densities(0., 0., 0.), material_0D.geometry.get_densities(0., 0., 0.)); assert_eq!(material_2D.geometry.get_total_density(0., 0., 0.), material_0D.geometry.get_total_density(0., 0., 0.)); @@ -487,6 +532,14 @@ fn test_geometry() { assert_eq!(material_2D.geometry.closest_point(-10., 0., 5.), material_0D.geometry.closest_point(-10., 0., 5.)); assert_eq!(material_2D.geometry.get_densities(-10., 0., 5.), material_0D.geometry.get_densities(-10., 0., 5.)); assert_eq!(material_2D.geometry.get_ck(-10., 0., 5.), material_0D.geometry.get_ck(-10., 0., 5.)); + + assert_eq!(material_2D.geometry.get_ck(0., 0., 0.), material_homogeneous_2D.geometry.get_ck(0., 0., 0.)); + assert_eq!(material_2D.geometry.get_densities(0., 0., 0.), material_homogeneous_2D.geometry.get_densities(0., 0., 0.)); + assert_eq!(material_2D.geometry.get_total_density(0., 0., 0.), material_homogeneous_2D.geometry.get_total_density(0., 0., 0.)); + assert_eq!(material_2D.geometry.get_concentrations(0., 0., 0.), material_homogeneous_2D.geometry.get_concentrations(0., 0., 0.)); + assert_eq!(material_2D.geometry.closest_point(-10., 0., 5.), material_homogeneous_2D.geometry.closest_point(-10., 0., 5.)); + assert_eq!(material_2D.geometry.get_densities(-10., 0., 5.), material_homogeneous_2D.geometry.get_densities(-10., 0., 5.)); + assert_eq!(material_2D.geometry.get_ck(-10., 0., 5.), material_homogeneous_2D.geometry.get_ck(-10., 0., 5.)); } #[test] @@ -502,7 +555,7 @@ fn test_surface_binding_energy_barrier() { let cosx = 1./(2.0_f64).sqrt(); let cosy = 1./(2.0_f64).sqrt(); let cosz = 0.; - let mut particle_1 = particle::Particle::new(mass, Z, E, Ec, Es, x, y, z, cosx, cosy, cosz, false, false, 0); + let mut particle_1 = particle::Particle::new(mass, Z, E, Ec, Es, 0.0, x, y, z, cosx, cosy, cosz, false, false, 0); let material_parameters = material::MaterialParameters{ energy_unit: "EV".to_string(), @@ -510,6 +563,7 @@ fn test_surface_binding_energy_barrier() { Eb: vec![0.0, 0.0], Es: vec![2.0, 4.0], Ec: vec![1.0, 1.0], + Ed: vec![0.0, 0.0], Z: vec![29., 1.], m: vec![63.54, 1.0008], interaction_index: vec![0, 0], @@ -647,7 +701,7 @@ fn test_surface_refraction() { let cosx = 1./(2.0_f64).sqrt(); let cosy = 1./(2.0_f64).sqrt(); let cosz = 0.; - let mut particle_1 = particle::Particle::new(mass, Z, E, Ec, Es, x, y, z, cosx, cosy, cosz, false, false, 0); + let mut particle_1 = particle::Particle::new(mass, Z, E, Ec, Es, 0.0, x, y, z, cosx, cosy, cosz, false, false, 0); //Test particle entering material and gaining energy @@ -736,6 +790,7 @@ fn test_momentum_conservation() { Eb: vec![0.0], Es: vec![Es2], Ec: vec![Ec2], + Ed: vec![0.0], Z: vec![Z2], m: vec![m2], interaction_index: vec![0], @@ -765,7 +820,7 @@ fn test_momentum_conservation() { println!("Case: {} {} {} {}", energy_eV, high_energy_free_flight_paths, potential, scattering_integral); - let mut particle_1 = particle::Particle::new(m1, Z1, E1, Ec1, Es1, x1, y1, z1, cosx, cosy, cosz, false, false, 0); + let mut particle_1 = particle::Particle::new(m1, Z1, E1, Ec1, Es1, 0.0, x1, y1, z1, cosx, cosy, cosz, false, false, 0); #[cfg(not(feature = "distributions"))] let options = Options { @@ -860,8 +915,9 @@ fn test_momentum_conservation() { binary_collision_geometries[0].phi_azimuthal); //Subtract total energy from all simultaneous collisions and electronic stopping - bca::update_particle_energy(&mut particle_1, &material_1, 0., - binary_collision_result.recoil_energy, 0., particle_2.Z, species_index, &options); + particle_1.E += -binary_collision_result.recoil_energy; + bca::subtract_electronic_stopping_energy(&mut particle_1, &material_1, 0., + 0., particle_2.Z, species_index, &options); let mom1_1 = particle_1.get_momentum(); let mom2_1 = particle_2.get_momentum(); @@ -897,6 +953,7 @@ fn test_rotate() { let E = 1.; let Ec = 1.; let Es = 1.; + let Ed = 0.0; let x = 0.; let y = 0.; let z = 0.; @@ -906,7 +963,7 @@ fn test_rotate() { let psi = -PI/4.; let phi = 0.; - let mut particle = particle::Particle::new(mass, Z, E, Ec, Es, x, y, z, cosx, cosy, cosz, false, false, 0); + let mut particle = particle::Particle::new(mass, Z, E, Ec, Es, Ed, x, y, z, cosx, cosy, cosz, false, false, 0); //Check that rotation in 2D works particle.rotate(psi, phi); @@ -936,6 +993,7 @@ fn test_particle_advance() { let E = 1.; let Ec = 1.; let Es = 1.; + let Ed = 0.0; let x = 0.; let y = 0.; let z = 0.; @@ -945,7 +1003,7 @@ fn test_particle_advance() { let mfp = 1.; let asymptotic_deflection = 0.5; - let mut particle = particle::Particle::new(mass, Z, E, Ec, Es, x, y, z, cosx, cosy, cosz, false, false, 0); + let mut particle = particle::Particle::new(mass, Z, E, Ec, Es, Ed, x, y, z, cosx, cosy, cosz, false, false, 0); let distance_traveled = particle.advance(mfp, asymptotic_deflection); @@ -1027,8 +1085,8 @@ fn test_quadrature() { let x0_newton = bca::newton_rootfinder(Za, Zb, Ma, Mb, E0, p, InteractionPotential::KR_C, 100, 1E-12).unwrap(); //If cpr_rootfinder is enabled, compare Newton to CPR - they should be nearly identical - #[cfg(any(feature = "cpr_rootfinder_openblas", feature = "cpr_rootfinder_netlib", feature = "cpr_rootfinder_intel_mkl"))] - if let Ok(x0_cpr) = bca::cpr_rootfinder(Za, Zb, Ma, Mb, E0, p, InteractionPotential::KR_C, 2, 1000, 1E-13, 1E-16, 1E-18, 1E6, 1E-18, false) { + #[cfg(feature = "cpr_rootfinder")] + if let Ok(x0_cpr) = bca::cpr_rootfinder(Za, Zb, Ma, Mb, E0, p, InteractionPotential::KR_C, 2, 10000, 1E-6, 1E-6, 1E-9, 1E9, 1E-13, true) { println!("CPR: {} Newton: {}", x0_cpr, x0_newton); assert!(approx_eq!(f64, x0_newton, x0_cpr, epsilon=1E-3)); }; @@ -1046,4 +1104,4 @@ fn test_quadrature() { println!("Gauss-Mehler: {} Gauss-Legendre: {} Mendenhall-Weller: {} MAGIC: {}", theta_gm, theta_gl, theta_mw, theta_magic); -} +} \ No newline at end of file