Merge dev into Main for version 1.2.0

.github/workflows/rustbca_compile_check.yml

@@ -27,23 +27,33 @@ jobs:
         sudo apt-get install rustc
     - name: Install pip for Python-3
       run: |
-        sudo apt-get install python3-pip
+        sudo apt-get install python3-pip python3-dev
     - name: Install Python libraries
       run: |
         python3 -m pip install numpy shapely scipy
     - name: Install python TOML library from source
       run: |
         git clone https://github.com/uiri/toml.git
         cd toml
-        sudo python3 setup.py install
+        python3 setup.py install --root .
     - name: Install HDF5 Libraries
       run: |
         sudo apt install libhdf5-dev
     - name: test Python Bindings
       run: |
-        sudo python3 -m pip install setuptools_rust testresources
-        sudo python3 setup.py install
-        python3 -c "from libRustBCA.pybca import *; print(simple_bca_py)"
+        python3 -m pip install setuptools_rust testresources
+        python3 -m pip install .
+        python3 -c "from libRustBCA.pybca import *;"
+    - name: Test Fortran and C bindings
+      run : |
+        cargo build --release
+        cp examples/test_rustbca.f90 .
+        gfortran -c rustbca.f90 target/release/liblibRustBCA.so
+        gfortran test_rustbca.f90 rustbca.f90 target/release/liblibRustBCA.so
+        ./a.out
+        cp examples/RustBCA.c .
+        g++ RustBCA.c RustBCA.h target/release/liblibRustBCA.so -Iexamples/ -I/usr/include/python3.8
+        ./a.out
     - name: Test RustBCA
       run: |
         sudo cargo test --features cpr_rootfinder_netlib,hdf5_input,distributions,parry3d

Cargo.toml

@@ -1,6 +1,7 @@
 [package]
 name = "RustBCA"
-version = "1.1.1"
+version = "1.2.0"
+
 authors = ["Jon Drobny <drobny2@illinois.edu>"]
 edition = "2018"
 
@@ -14,6 +15,7 @@ crate-type = ["cdylib"]
 
 [dependencies]
 rand = "0.8.3"
+rand_distr = "0.4.2"
 toml = "0.5.8"
 anyhow = "1.0.38"
 itertools = "0.10.0"

RustBCA.h

@@ -97,8 +97,38 @@ struct InputCompoundBCA {
   double *Eb2;
 };
 
+struct OutputTaggedBCA {
+  uintptr_t len;
+  double (*particles)[9];
+  double *weights;
+  int32_t *tags;
+};
+
+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);

examples/RustBCA.c

@@ -0,0 +1,54 @@
+#include "RustBCA.h"
+#include <iostream>
+#include <vector>
+
+int main(int argc, char * argv[]) {
+  OutputTaggedBCA output;
+  double velocities[2][3] = {{500000.0, 0.1, 0.0}, {500000.0, 0.1, 0.0}};
+  double positions[2][3] = {{0.0, 0.0, 0.0}, {1.0, 1.0, 1.0}};
+  int tags[2] = {0, 1};
+  double weights[2] = {1.0, 1.0};
+  double Z[3] = {74.0, 74.0};
+  double m[2] = {184.0, 184.0};
+  double n[2] = {0.06306, 0.06306};
+  double Ec[2] = {1.0, 1.0};
+  double Es[2] = {8.79, 8.79};
+  double Eb[2] = {0.0, 0.0};
+
+  InputTaggedBCA input = {
+    2,
+    positions,
+    velocities,
+    1.0,
+    1.0,
+    1.0,
+    1.0,
+    2,
+    Z,
+    m,
+    n,
+    Ec,
+    Es,
+    Eb,
+    tags,
+    weights
+  };
+
+  //output = simple_bca_c(0., 0., 0., 0.5, 0.5, 0.00, 2000.0, 2.0, 4.0, 1.0, 0.0, 74.0, 184.0, 1.0, 8.79, 0.06306, 0.0);
+  //output = compound_bca_list_c(input);
+  output = compound_tagged_bca_list_c(input);
+
+  std::cout << "Particle 1 Z: ";
+  std::cout << output.particles[0][0];
+  std::cout << std::endl;
+  std::cout << "Particle 1 E [eV]: ";
+  std::cout << output.particles[0][2];
+  std::cout << std::endl;
+  std::cout << "Particle 2 Z: ";
+  std::cout << output.particles[1][0];
+  std::cout << std::endl;
+  std::cout << "Particle 2 E [eV]: ";
+  std::cout << output.particles[1][2];
+  std::cout << std::endl;
+  return 0;
+}

examples/test_rustbca.f90

@@ -0,0 +1,75 @@
+
+program test_rustbca
+
+    use rustbca
+    use, intrinsic :: iso_c_binding
+
+    integer :: N_ions
+    real(c_double), allocatable, dimension(:) :: ux, uy, uz, E, Z1, m1, Ec1, Es1
+    integer(c_int) :: num_species_target, num_emitted_particles
+    real(c_double), target :: Z2(2), m2(2), Ec2(2), Es2(2), Eb2(2), n2(2)
+    real(c_double) :: ux1, uy1, uz1, E1
+    type(c_ptr) :: bca_output_c
+    real(c_double), pointer, dimension(:,:) :: bca_output_f
+    real :: start, stop
+    logical(c_bool) :: track_recoils
+
+    !Initial ion conditions
+    N_ions = 100000
+    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)
+    uz(:) = 0.0
+    E(:) = 1000.0_8
+
+    !Hydrogen
+    Z1(:) = 1.0_8
+    m1(:) = 1.008_8
+    Ec1(:) = 1.0_8
+    Es1(:) = 1.5_8
+
+    !Titanium Hydride
+    num_species_target = 2
+    Z2(1) = 22.0_8
+    m2(1) = 47.867_8
+    Ec2(1) = 4.84_8
+    Es2(1) = 4.84_8
+    Eb2(1) = 3.0_8
+    n2(1) = 0.04527_8
+
+    Z2(2) = 1.0_8
+    m2(2) = 1.008_8
+    Ec2(2) = 1.5_8
+    Es2(2) = 1.5_8
+    Eb2(2) = 0.0_8
+    n2(2) = 0.09054_8
+
+    track_recoils = .false.
+
+    call cpu_time(start)
+    bca_output_c = compound_bca_list_fortran(N_ions, track_recoils, ux, uy, uz, E, &
+        Z1, m1, Ec1, Es1, &
+        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
+
+    call cpu_time(start)
+    do i = 0, N_ions
+        !Test reflect_single_ion routine
+        ux1 = 0.999
+        uy1 = sqrt(1.0 - 0.999*0.999)
+        uz1 = 0.0
+        E1 = 1000.0
+        call reflect_single_ion_c(num_species_target, ux1, uy1, uz1, E1, Z1(1), m1(1), Ec1(1), Es1(1), Z2, m2, Ec2, Es2, Eb2, n2)
+    end do
+    call cpu_time(stop)
+    write(*,*) "Elapsed time in ions per eV per s: ", (stop - start)/N_ions/1000.0
+
+    !call exit(1)
+
+end program test_rustbca

rustbca.f90

@@ -0,0 +1,95 @@
+module rustbca
+
+    use, intrinsic :: iso_c_binding
+
+    interface
+
+        subroutine reflect_single_ion_c(num_species_target, ux, uy, uz, E1, &
+            Z1, m1, Ec1, Es1, Z2, m2, Ec2, Es2, Eb2, n2) bind(c)
+
+            !Runs a single ion BCA trajectory with no recoils
+            !Args:
+            !   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
+            !   E1 (real(c_double)): initial energy of incident ion in eV
+            !   Z1 (real(c_double)): atomic number of incident ion
+            !   m1 (real(c_double)): atomic mass of incident ion in eV
+            !   Ec1 (real(c_double)): cutoff energy of incident ion in eV
+            !   Es1 (real(c_double)): surface binding energy of incident ion in eV
+            !   Z2 (real(c_double), dimension(:)): list of atomic numbers of target speciesd
+            !   m2 (real(c_double), dimension(:)): list of atomic masses of target species in amu
+            !   Ec2 (real(c_double), dimension(:)): list of cutoff energies of target species in eV
+            !   Es2 (real(c_double), dimension(:)): list of surface binding energies of target species in eV
+            !   Eb2 (real(c_double), dimension(:)): list of bulk binding energies of target species in eV
+            !   n2 (real(c_double), dimension(:)): list of number densities of target species in 1/angstrom^3
+
+            use, intrinsic :: iso_c_binding
+            real(c_double), intent(inout) :: ux, uy, uz, E1
+            real(c_double), intent(in) :: Z1, m1, Ec1, Es1
+            integer(c_int), intent(in) :: num_species_target
+            real(c_double), intent(in), dimension(*) :: Z2, m2, Ec2, Es2, Eb2, n2
+
+        end subroutine reflect_single_ion_c
+
+        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
+            !   track_recoils (logical(c_bool)): whether to generate recoils (disable to turn off sputtering)
+            !   ux (real(c_double), dimension(:)): x-direction of incident ion, x-direction of reflected ion
+            !   uy (real(c_double), dimension(:)): y-direction of incident ion, y-direction of reflected ion
+            !   uz (real(c_double), dimension(:)): z-direction of incident ion, z-direction of reflected ion
+            !   E1 (real(c_double), dimension(:)): initial energy of incident ion in eV
+            !   Z1 (real(c_double), dimension(:)): atomic number of incident ion
+            !   m1 (real(c_double), dimension(:)): atomic mass of incident ion in eV
+            !   Ec1 (real(c_double), dimension(:)): cutoff energy of incident ion in eV
+            !   num_species_target(integer(c_int)): number of species in target
+            !   Es1 (real(c_double), dimension(:)): surface binding energy of incident ion in eV
+            !   Z2 (real(c_double), dimension(:)): list of atomic numbers of target speciesd
+            !   m2 (real(c_double), dimension(:)): list of atomic masses of target species in amu
+            !   Ec2 (real(c_double), dimension(:)): list of cutoff energies of target species in eV
+            !   Es2 (real(c_double), dimension(:)): list of surface binding energies of target species in eV
+            !   Eb2 (real(c_double), dimension(:)): list of bulk binding energies of target species in eV
+            !   n2 (real(c_double), dimension(:)): list of number densities of target species in 1/angstrom^3
+            !   num_emitted_particles (integer(c_int), intent(out)): NOTE THAT THIS IS INTENT(OUT) number of emitted particles in output
+            !Returns:
+            !   output (type(c_ptr)): a c pointer to a 2D array of size (num_emitted_particles, 6) that consists of Z, m, E, ux, uy, uz
+
+            use, intrinsic :: iso_c_binding
+            logical(c_bool), intent(in) :: track_recoils
+            integer(c_int), intent(in) :: num_incident_ions, num_species_target
+            integer(c_int), intent(out) :: num_emitted_particles
+            real(c_double), intent(in), dimension(*) :: ux, uy, uz, E1, Z1, m1, Ec1, Es1
+            real(c_double), intent(in), dimension(*) :: Z2, m2, Ec2, Es2, EB2, n2
+            type(c_ptr) :: output
+        end function compound_bca_list_fortran
+
+    end interface
+
+    contains
+
+        subroutine transform_to_local_angle(ux, uy, uz, alpha)
+
+            !Rotates a vector in 2D
+            !Args:
+            !   ux (real(c_double)): x-direction
+            !   uy (real(c_double)): y-direction
+            !   uz (real(c_double)): z-direction
+            !   alpha (real(c_double)): local surface angle measured counter-clockwise from x-axis in radians
+
+            real(8), intent(inout) :: ux, uy, uz
+            real(8), intent(in) :: alpha
+
+            ux = ux*cos(alpha) - uy*sin(alpha)
+            uy = ux*sin(alpha) + uy*cos(alpha)
+
+        end subroutine transform_to_local_angle
+
+end module rustbca

setup.py

@@ -3,7 +3,18 @@
 
 setup(
     name="RustBCA",
-    rust_extensions=[RustExtension("libRustBCA.pybca", binding=Binding.PyO3, features=["python"])],
+    version="1.1.2",
+    rust_extensions=[
+        RustExtension(
+            "libRustBCA.pybca",
+            binding=Binding.PyO3,
+            features=["python"],
+            #args=["+nightly", "--edition 2018", "-Z unstable-options"],
+            #optional=True,
+            rust_version="1.56.1"
+
+        )
+    ],
     # rust extensions are not zip safe, just like C-extensions.
     zip_safe=False,
 )

src/bca.rs

@@ -377,15 +377,19 @@ pub fn choose_collision_partner<T: Geometry>(particle_1: &particle::Particle, ma
 
     //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);
-
-    return (species_index,
-        particle::Particle::new(
-            M_recoil, Z_recoil, 0., Ec_recoil, Es_recoil,
-            x_recoil, y_recoil, z_recoil,
-            cosx, cosy, cosz,
-            false, options.track_recoil_trajectories, interaction_index
-        )
-    )
+    let mut new_particle = particle::Particle::new(
+        M_recoil, Z_recoil, 0., Ec_recoil, Es_recoil,
+        x_recoil, y_recoil, z_recoil,
+        cosx, cosy, cosz,
+        false, options.track_recoil_trajectories, interaction_index
+    );
+
+    // Carry through tag, weight, and tracked vector from originating particle
+    new_particle.weight = particle_1.weight;
+    new_particle.tag = particle_1.tag;
+    new_particle.tracked_vector = particle_1.tracked_vector;
+
+    return (species_index, new_particle)
 }
 
 /// Calculate the distance of closest approach of two particles given a particular binary collision geometry.
@@ -443,6 +447,7 @@ pub fn update_particle_energy<T: Geometry>(particle_1: &mut particle::Particle,
     recoil_energy: f64, x0: f64, strong_collision_Z: f64, strong_collision_index: usize, options: &Options) {
 
     //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. {
@@ -694,8 +699,8 @@ pub fn newton_rootfinder(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_par
             return Ok(x0);
         }
     }
-    return Err(anyhow!("Numerical error: exceeded maximum number of Newton-Raphson iterations, {}. E: {}; x0: {}; Error: {}; Tolerance: {}",
-        max_iterations, E0, x0, err, tolerance));
+    return Err(anyhow!("Numerical error: exceeded maximum number of Newton-Raphson iterations, {}. E: {} eV; x0: {}; Error: {}; Tolerance: {}; Za: {}; Zb: {}; Ma: {} amu; Mb: {} amu; a: {}; p: {} A",
+        max_iterations, E0/Q, x0, err, tolerance, Za, Zb, Ma/AMU, Mb/AMU, a, impact_parameter/ANGSTROM));
 }
 
 /// Gauss-Mehler quadrature.