Skip to content

Commit f32afa2

Browse files
drobnyjtdrobnyjt
authored
Merge pull request #170 from lcpp-org/dev
Merge dev into Main for version 1.2.0
2 parents 99bf81c + 20ad62e commit f32afa2

15 files changed

+921
-60
lines changed

‎.github/workflows/rustbca_compile_check.yml

+15-5
Original file line numberDiff line numberDiff line change
@@ -27,23 +27,33 @@ jobs:
2727
sudo apt-get install rustc
2828
- name: Install pip for Python-3
2929
run: |
30-
sudo apt-get install python3-pip
30+
sudo apt-get install python3-pip python3-dev
3131
- name: Install Python libraries
3232
run: |
3333
python3 -m pip install numpy shapely scipy
3434
- name: Install python TOML library from source
3535
run: |
3636
git clone https://github.com/uiri/toml.git
3737
cd toml
38-
sudo python3 setup.py install
38+
python3 setup.py install --root .
3939
- name: Install HDF5 Libraries
4040
run: |
4141
sudo apt install libhdf5-dev
4242
- name: test Python Bindings
4343
run: |
44-
sudo python3 -m pip install setuptools_rust testresources
45-
sudo python3 setup.py install
46-
python3 -c "from libRustBCA.pybca import *; print(simple_bca_py)"
44+
python3 -m pip install setuptools_rust testresources
45+
python3 -m pip install .
46+
python3 -c "from libRustBCA.pybca import *;"
47+
- name: Test Fortran and C bindings
48+
run : |
49+
cargo build --release
50+
cp examples/test_rustbca.f90 .
51+
gfortran -c rustbca.f90 target/release/liblibRustBCA.so
52+
gfortran test_rustbca.f90 rustbca.f90 target/release/liblibRustBCA.so
53+
./a.out
54+
cp examples/RustBCA.c .
55+
g++ RustBCA.c RustBCA.h target/release/liblibRustBCA.so -Iexamples/ -I/usr/include/python3.8
56+
./a.out
4757
- name: Test RustBCA
4858
run: |
4959
sudo cargo test --features cpr_rootfinder_netlib,hdf5_input,distributions,parry3d

‎Cargo.toml

+3-1
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,7 @@
11
[package]
22
name = "RustBCA"
3-
version = "1.1.1"
3+
version = "1.2.0"
4+
45
authors = ["Jon Drobny <[email protected]>"]
56
edition = "2018"
67

@@ -14,6 +15,7 @@ crate-type = ["cdylib"]
1415

1516
[dependencies]
1617
rand = "0.8.3"
18+
rand_distr = "0.4.2"
1719
toml = "0.5.8"
1820
anyhow = "1.0.38"
1921
itertools = "0.10.0"

‎RustBCA.h

+30
Original file line numberDiff line numberDiff line change
@@ -97,8 +97,38 @@ struct InputCompoundBCA {
9797
double *Eb2;
9898
};
9999

100+
struct OutputTaggedBCA {
101+
uintptr_t len;
102+
double (*particles)[9];
103+
double *weights;
104+
int32_t *tags;
105+
};
106+
107+
struct InputTaggedBCA {
108+
uintptr_t len;
109+
/// x y z
110+
double (*positions)[3];
111+
/// vx, vy, vz
112+
double (*velocities)[3];
113+
double Z1;
114+
double m1;
115+
double Ec1;
116+
double Es1;
117+
uintptr_t num_species_target;
118+
double *Z2;
119+
double *m2;
120+
double *n2;
121+
double *Ec2;
122+
double *Es2;
123+
double *Eb2;
124+
int32_t *tags;
125+
double *weights;
126+
};
127+
100128
extern "C" {
101129

130+
OutputTaggedBCA compound_tagged_bca_list_c(InputTaggedBCA input);
131+
102132
OutputBCA simple_bca_list_c(InputSimpleBCA input);
103133

104134
OutputBCA compound_bca_list_c(InputCompoundBCA input);

‎examples/RustBCA.c

+54
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,54 @@
1+
#include "RustBCA.h"
2+
#include <iostream>
3+
#include <vector>
4+
5+
int main(int argc, char * argv[]) {
6+
OutputTaggedBCA output;
7+
double velocities[2][3] = {{500000.0, 0.1, 0.0}, {500000.0, 0.1, 0.0}};
8+
double positions[2][3] = {{0.0, 0.0, 0.0}, {1.0, 1.0, 1.0}};
9+
int tags[2] = {0, 1};
10+
double weights[2] = {1.0, 1.0};
11+
double Z[3] = {74.0, 74.0};
12+
double m[2] = {184.0, 184.0};
13+
double n[2] = {0.06306, 0.06306};
14+
double Ec[2] = {1.0, 1.0};
15+
double Es[2] = {8.79, 8.79};
16+
double Eb[2] = {0.0, 0.0};
17+
18+
InputTaggedBCA input = {
19+
2,
20+
positions,
21+
velocities,
22+
1.0,
23+
1.0,
24+
1.0,
25+
1.0,
26+
2,
27+
Z,
28+
m,
29+
n,
30+
Ec,
31+
Es,
32+
Eb,
33+
tags,
34+
weights
35+
};
36+
37+
//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);
38+
//output = compound_bca_list_c(input);
39+
output = compound_tagged_bca_list_c(input);
40+
41+
std::cout << "Particle 1 Z: ";
42+
std::cout << output.particles[0][0];
43+
std::cout << std::endl;
44+
std::cout << "Particle 1 E [eV]: ";
45+
std::cout << output.particles[0][2];
46+
std::cout << std::endl;
47+
std::cout << "Particle 2 Z: ";
48+
std::cout << output.particles[1][0];
49+
std::cout << std::endl;
50+
std::cout << "Particle 2 E [eV]: ";
51+
std::cout << output.particles[1][2];
52+
std::cout << std::endl;
53+
return 0;
54+
}

‎examples/test_rustbca.f90

+75
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,75 @@
1+
2+
program test_rustbca
3+
4+
use rustbca
5+
use, intrinsic :: iso_c_binding
6+
7+
integer :: N_ions
8+
real(c_double), allocatable, dimension(:) :: ux, uy, uz, E, Z1, m1, Ec1, Es1
9+
integer(c_int) :: num_species_target, num_emitted_particles
10+
real(c_double), target :: Z2(2), m2(2), Ec2(2), Es2(2), Eb2(2), n2(2)
11+
real(c_double) :: ux1, uy1, uz1, E1
12+
type(c_ptr) :: bca_output_c
13+
real(c_double), pointer, dimension(:,:) :: bca_output_f
14+
real :: start, stop
15+
logical(c_bool) :: track_recoils
16+
17+
!Initial ion conditions
18+
N_ions = 100000
19+
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))
20+
ux(:) = 0.999
21+
uy(:) = sqrt(1.0 - 0.999*0.999)
22+
uz(:) = 0.0
23+
E(:) = 1000.0_8
24+
25+
!Hydrogen
26+
Z1(:) = 1.0_8
27+
m1(:) = 1.008_8
28+
Ec1(:) = 1.0_8
29+
Es1(:) = 1.5_8
30+
31+
!Titanium Hydride
32+
num_species_target = 2
33+
Z2(1) = 22.0_8
34+
m2(1) = 47.867_8
35+
Ec2(1) = 4.84_8
36+
Es2(1) = 4.84_8
37+
Eb2(1) = 3.0_8
38+
n2(1) = 0.04527_8
39+
40+
Z2(2) = 1.0_8
41+
m2(2) = 1.008_8
42+
Ec2(2) = 1.5_8
43+
Es2(2) = 1.5_8
44+
Eb2(2) = 0.0_8
45+
n2(2) = 0.09054_8
46+
47+
track_recoils = .false.
48+
49+
call cpu_time(start)
50+
bca_output_c = compound_bca_list_fortran(N_ions, track_recoils, ux, uy, uz, E, &
51+
Z1, m1, Ec1, Es1, &
52+
num_species_target, Z2, m2, Ec2, Es2, Eb2, n2, &
53+
num_emitted_particles)
54+
call c_f_pointer(bca_output_c, bca_output_f, [num_emitted_particles, 6])
55+
call cpu_time(stop)
56+
57+
write(*,*) "Elapsed time in seconds per ion per eV: ", (stop - start)/N_ions/1000.0
58+
59+
!write(*,*) bca_output_f
60+
61+
call cpu_time(start)
62+
do i = 0, N_ions
63+
!Test reflect_single_ion routine
64+
ux1 = 0.999
65+
uy1 = sqrt(1.0 - 0.999*0.999)
66+
uz1 = 0.0
67+
E1 = 1000.0
68+
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)
69+
end do
70+
call cpu_time(stop)
71+
write(*,*) "Elapsed time in ions per eV per s: ", (stop - start)/N_ions/1000.0
72+
73+
!call exit(1)
74+
75+
end program test_rustbca

‎rustbca.f90

+95
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,95 @@
1+
module rustbca
2+
3+
use, intrinsic :: iso_c_binding
4+
5+
interface
6+
7+
subroutine reflect_single_ion_c(num_species_target, ux, uy, uz, E1, &
8+
Z1, m1, Ec1, Es1, Z2, m2, Ec2, Es2, Eb2, n2) bind(c)
9+
10+
!Runs a single ion BCA trajectory with no recoils
11+
!Args:
12+
! num_species_target (integer): number of species in target
13+
! ux (real(c_double)): x-direction of incident ion, x-direction of reflected ion
14+
! uy (real(c_double)): y-direction of incident ion, y-direction of reflected ion
15+
! uz (real(c_double)): z-direction of incident ion, z-direction of reflected ion
16+
! E1 (real(c_double)): initial energy of incident ion in eV
17+
! Z1 (real(c_double)): atomic number of incident ion
18+
! m1 (real(c_double)): atomic mass of incident ion in eV
19+
! Ec1 (real(c_double)): cutoff energy of incident ion in eV
20+
! Es1 (real(c_double)): surface binding energy of incident ion in eV
21+
! Z2 (real(c_double), dimension(:)): list of atomic numbers of target speciesd
22+
! m2 (real(c_double), dimension(:)): list of atomic masses of target species in amu
23+
! Ec2 (real(c_double), dimension(:)): list of cutoff energies of target species in eV
24+
! Es2 (real(c_double), dimension(:)): list of surface binding energies of target species in eV
25+
! Eb2 (real(c_double), dimension(:)): list of bulk binding energies of target species in eV
26+
! n2 (real(c_double), dimension(:)): list of number densities of target species in 1/angstrom^3
27+
28+
use, intrinsic :: iso_c_binding
29+
real(c_double), intent(inout) :: ux, uy, uz, E1
30+
real(c_double), intent(in) :: Z1, m1, Ec1, Es1
31+
integer(c_int), intent(in) :: num_species_target
32+
real(c_double), intent(in), dimension(*) :: Z2, m2, Ec2, Es2, Eb2, n2
33+
34+
end subroutine reflect_single_ion_c
35+
36+
function compound_bca_list_fortran(num_incident_ions, track_recoils, ux, uy, uz, E1, &
37+
Z1, m1, Ec1, Es1, &
38+
num_species_target, Z2, m2, Ec2, Es2, Eb2, n2, &
39+
num_emitted_particles) bind(c) result(output)
40+
41+
42+
!Runs a homogeneous, flat, compound target BCA with an arbitrary list of ions.
43+
!Args:
44+
! num_incident_ion (integer(c_int)): number of incident ions
45+
! track_recoils (logical(c_bool)): whether to generate recoils (disable to turn off sputtering)
46+
! ux (real(c_double), dimension(:)): x-direction of incident ion, x-direction of reflected ion
47+
! uy (real(c_double), dimension(:)): y-direction of incident ion, y-direction of reflected ion
48+
! uz (real(c_double), dimension(:)): z-direction of incident ion, z-direction of reflected ion
49+
! E1 (real(c_double), dimension(:)): initial energy of incident ion in eV
50+
! Z1 (real(c_double), dimension(:)): atomic number of incident ion
51+
! m1 (real(c_double), dimension(:)): atomic mass of incident ion in eV
52+
! Ec1 (real(c_double), dimension(:)): cutoff energy of incident ion in eV
53+
! num_species_target(integer(c_int)): number of species in target
54+
! Es1 (real(c_double), dimension(:)): surface binding energy of incident ion in eV
55+
! Z2 (real(c_double), dimension(:)): list of atomic numbers of target speciesd
56+
! m2 (real(c_double), dimension(:)): list of atomic masses of target species in amu
57+
! Ec2 (real(c_double), dimension(:)): list of cutoff energies of target species in eV
58+
! Es2 (real(c_double), dimension(:)): list of surface binding energies of target species in eV
59+
! Eb2 (real(c_double), dimension(:)): list of bulk binding energies of target species in eV
60+
! n2 (real(c_double), dimension(:)): list of number densities of target species in 1/angstrom^3
61+
! num_emitted_particles (integer(c_int), intent(out)): NOTE THAT THIS IS INTENT(OUT) number of emitted particles in output
62+
!Returns:
63+
! 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
64+
65+
use, intrinsic :: iso_c_binding
66+
logical(c_bool), intent(in) :: track_recoils
67+
integer(c_int), intent(in) :: num_incident_ions, num_species_target
68+
integer(c_int), intent(out) :: num_emitted_particles
69+
real(c_double), intent(in), dimension(*) :: ux, uy, uz, E1, Z1, m1, Ec1, Es1
70+
real(c_double), intent(in), dimension(*) :: Z2, m2, Ec2, Es2, EB2, n2
71+
type(c_ptr) :: output
72+
end function compound_bca_list_fortran
73+
74+
end interface
75+
76+
contains
77+
78+
subroutine transform_to_local_angle(ux, uy, uz, alpha)
79+
80+
!Rotates a vector in 2D
81+
!Args:
82+
! ux (real(c_double)): x-direction
83+
! uy (real(c_double)): y-direction
84+
! uz (real(c_double)): z-direction
85+
! alpha (real(c_double)): local surface angle measured counter-clockwise from x-axis in radians
86+
87+
real(8), intent(inout) :: ux, uy, uz
88+
real(8), intent(in) :: alpha
89+
90+
ux = ux*cos(alpha) - uy*sin(alpha)
91+
uy = ux*sin(alpha) + uy*cos(alpha)
92+
93+
end subroutine transform_to_local_angle
94+
95+
end module rustbca

‎setup.py

+12-1
Original file line numberDiff line numberDiff line change
@@ -3,7 +3,18 @@
33

44
setup(
55
name="RustBCA",
6-
rust_extensions=[RustExtension("libRustBCA.pybca", binding=Binding.PyO3, features=["python"])],
6+
version="1.1.2",
7+
rust_extensions=[
8+
RustExtension(
9+
"libRustBCA.pybca",
10+
binding=Binding.PyO3,
11+
features=["python"],
12+
#args=["+nightly", "--edition 2018", "-Z unstable-options"],
13+
#optional=True,
14+
rust_version="1.56.1"
15+
16+
)
17+
],
718
# rust extensions are not zip safe, just like C-extensions.
819
zip_safe=False,
920
)

‎src/bca.rs

+16-11
Original file line numberDiff line numberDiff line change
@@ -377,15 +377,19 @@ pub fn choose_collision_partner<T: Geometry>(particle_1: &particle::Particle, ma
377377

378378
//Choose recoil Z, M
379379
let (species_index, Z_recoil, M_recoil, Ec_recoil, Es_recoil, interaction_index) = material.choose(x_recoil, y_recoil, z_recoil);
380-
381-
return (species_index,
382-
particle::Particle::new(
383-
M_recoil, Z_recoil, 0., Ec_recoil, Es_recoil,
384-
x_recoil, y_recoil, z_recoil,
385-
cosx, cosy, cosz,
386-
false, options.track_recoil_trajectories, interaction_index
387-
)
388-
)
380+
let mut new_particle = particle::Particle::new(
381+
M_recoil, Z_recoil, 0., Ec_recoil, Es_recoil,
382+
x_recoil, y_recoil, z_recoil,
383+
cosx, cosy, cosz,
384+
false, options.track_recoil_trajectories, interaction_index
385+
);
386+
387+
// Carry through tag, weight, and tracked vector from originating particle
388+
new_particle.weight = particle_1.weight;
389+
new_particle.tag = particle_1.tag;
390+
new_particle.tracked_vector = particle_1.tracked_vector;
391+
392+
return (species_index, new_particle)
389393
}
390394

391395
/// 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,
443447
recoil_energy: f64, x0: f64, strong_collision_Z: f64, strong_collision_index: usize, options: &Options) {
444448

445449
//If particle energy drops below zero before electronic stopping calcualtion, it produces NaNs
450+
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);
446451
particle_1.E -= recoil_energy;
447452
assert!(!particle_1.E.is_nan(), "Numerical error: particle energy is NaN following collision.");
448453
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
694699
return Ok(x0);
695700
}
696701
}
697-
return Err(anyhow!("Numerical error: exceeded maximum number of Newton-Raphson iterations, {}. E: {}; x0: {}; Error: {}; Tolerance: {}",
698-
max_iterations, E0, x0, err, tolerance));
702+
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",
703+
max_iterations, E0/Q, x0, err, tolerance, Za, Zb, Ma/AMU, Mb/AMU, a, impact_parameter/ANGSTROM));
699704
}
700705

701706
/// Gauss-Mehler quadrature.

0 commit comments

Comments
 (0)