lcpp-org/rustbca

A free, open-source Binary Collision Approximation (BCA) code for ion-material interactions including sputtering, implantation, and reflection

rustbca's Issues

Should electronic stopping correction factor be implemented on a per-species basis?

The correction factor is empirical, so I'm not sure if it should be a per-species term, multiplied inside the sum for Bragg's rule, or a material-wide term, multiplied once outside.

[question] bulk binding energy calculation

The bulk binding energy is currently computed as the composition-weighted average of the local bulk binding energies. Obviously this is not strictly correct, but should this be replaced with an option in the style of the surface binding energy? That is, an option to average local or use individual bulk binding energies?

[feature] investigate implementing arbitrary structs for testing suite instead of relying on pre-determined floats

Is your feature request related to a problem? Please describe.
Testing currently relies on so-called "magic numbers" - that is, there are hand-written floats in the testing code. These are intended to cover all normal cases and any edge cases, but this is not necessarily guaranteed.

Proposed solution
Implement the Arbitrary crate's Arbitrary trait on structs such as Particle, Mesh2D, Material, etc.

Alternative solution(s)
Look into automatic testing, property testing, or some other provable testing solution.

Additional context
https://github.com/rust-fuzz/arbitrary

Test issue

[question] What dependencies are necessary to build rcpr, and what is the best way to do so?

See this comment.

rcpr is still early in development, and building it is tricky. What dependencies are actually required to build it as-is? How can those be slimmed down/removed? Is there a standalone eigenvalue algorithm that could be used instead of the heavy numerical backends? Is it worthwhile to write a custom eigenvalue solver in rcpr?

[feature] New output file: summary of yield, reflection coefficient, etc.

Test issue

This is a test of Slack integration of RustBCA.

Unit input is not consistent

If a user puts "CM" as the units, it appears that "energy barrier thickness" is unaffected; as is the densities under 2d_mesh_input.

Avoid large TOML ascii input files at large particle count

Using SD.Trim.SP lookup tables, a frozen oxygen target appears to break the recoil tracking

It's unclear which parameter is causing the code to go haywire; it appears to only be an issue for frozen, gaseous targets.

Add SDTrimSP's "model3" for surface binding energies

Currently, surface binding energies in the code are calculated as follows:

If the particle is an incident particle:
Es = particle.Es

If the particle is a particle from the target:
Es = concentration[0]*Es[0] + concentration[1]*Es[1] + ... + concentration[N]*Es[N]
Where N is the number of species in the target.

This neglects 1) the possibility of the incident particle having different surface binding energies with different pieces of the target (i.e., different triangles in the mesh with different compositions) which will affect the range very slightly and the reflection coefficient at energies ~Es and 2) the possibility of pairwise surface binding energy interactions.

Ensure planar surface binding model is functioning for certain transmitted particles

Currently, the surface binding refraction is calculated using a modified Snell-type law and a single rotation in the x-y plane. However, I have some doubts that, for particles transmitted out of the back or bottom of the target, that this will always refract particles correctly - that is, always ensure that they leave with an angle nearer the surface tangent.

Code needs structure, on path towards dynamic library BCA

https://cheats.rs/#organizing-code

Implement seedable RNG that scales to many cores

Currently, relying on the system RNG means that the code is not repeatable. It should be switched to a threaded RNG that supports seeding.

Best way to cite? [question]

What's the best way to cite RustBCA?

[feature] Architecture.md

Is your feature request related to a problem? Please describe.
An architecture.md file would be really nice to have, especially as the project grows in scope.

Proposed solution
ARCHITECTURE.md

Error when using a new interaction potential

I attempted to add a new interaction potential based on results from a DFT code. The best fit of the data is a modified Lennard-Jones potential:

16.0*( (1.87/xi).powf(6.5) - (1.87/xi).powf(6.0) )

After our conversation yesterday, I attempted to convert this V(r) potential to a dimensionless phi(r) potential by multiplying by r, without worrying about any additional scaling factor.

xi*(16.0*( (1.87/xi).powf(6.5) - (1.87/xi).powf(6.0) ))

This code failed to run, seen in the image below, and I haven't been able to find the backtrace variable to obtain a better backtrace. Let me know if there's any additional information I can give that would be helpful when diagnosing this problem.

Check geometry of this input file

input.txt

Geometry error

[feature] mesh_2d requirement should be replaced with a generic Geometry type

Is your feature request related to a problem? Please describe.
Currently, there is no way to, for example, run rectangular geometry with fewer than 2 triangles. Layered geometry requires many triangle checks which are not technically necessary - a simple bounding box calculation would do.

Proposed solution
The mesh_2D type in Material should be replaced with a generic Geometry type that implements all the necessary traits - determining the composition, whether or not a particle is outside the material, etc.

Additional context
https://doc.rust-lang.org/book/ch10-01-syntax.html

[feature] Add an abstraction layer to simplify calling routines from other codes

Is your feature request related to a problem? Please describe.
Currently, issue #62 concerns adding binning of particle output for coupling to other codes. However, I'm concerned about adding some sort of binning/post-processing routine to RustBCA proper because of future maintenance costs and feature creep.

Proposed solution
However, a good solution would be to add an abstraction layer that includes the input file, so that a separate code can just use RustBCA as a library and call it and do the processing in memory straight out of single_ion_bca(). This would also have the benefit of simplifying main.rs and adding input.rs and potentially output.rs for code readability.

Alternative solution(s)
Python wrapper of Rust code would be simpler, but can streaming data from RustBCA be done succinctly?

Unit parser should be extended to arbitrary values

Instead of restricting to a list of pre-set values, using the parse() method it should be possible to attempt to convert the string to a float and only panic if that fails. See:

https://doc.rust-lang.org/std/string/struct.String.html#method.parse

Firsov/Oen-Robinson stopping needs access to interaction potential to determine correct screening length

[feature] Enhance readability/performance with BCA refactor

Is your feature request related to a problem? Please describe.
Currently, the code has a number of near-hacks to perform the binary collision approximation trajectory steps in the correct order. These include storing previous directions/positions in the Particle object, having the electronic stopping calculation broken up for local/nonlocal forms, potentially redundant boundary condition checks, etc.

Proposed solution
Further separating out certain procedures into functions would be helpful. Additionally, an analysis of the BCA order of operations from a programming perspective should be performed prior to the refactor to ensure correctness and simplicity.

There are many uses of badly-named variables in closures

e.g., in the effective Z calculation, return self.Z.iter().zip(concentrations).map(|(i1, i2)| i1*i2).collect::<Vec<f64>>().iter().sum(); uses "i1" and "i2" - for clarity's sake, these should probably be relevant variable names instead of these placeholders

ZBL potential should use alternative screening length

The ZBL potential requires a different screening length than the implemented Lindhard/Firsov one.

Moliere potential with MAGIC algorithm does not work

I suspect the transcription error lies in Ziegler's book with C4/C5 being switched around, and that is breaking the scattering angle calculation

Add useful documentation comments

https://doc.rust-lang.org/stable/book/ch14-02-publishing-to-crates-io.html?highlight=docs#making-useful-documentation-comments

Consider a doctest suite to ensure documentation stays up to date

Setting E = 0 Sometimes causes code to halt

Consider switching mean free path calculation from sim-wide option to triangle-local option

I've been going back and forth on this one; for a complex mesh, having a massive array of enums ["GASEOUS", "GASEOUS", "GASEOUS", ...] would look really bad - but does it make the input file any less usable than having, say, lists of 100s of floating point numbers (current behavior?)

Oen-Robinson stopping doesn't reference radii in interactions.rs

bca.rs: 355

[feature] Triangle lists should follow GetDP-like structure

Is your feature request related to a problem? Please describe.
Current implementation of triangle lists is unreadable and difficult to manage in code and by hand.

Proposed solution
Switch to:

1. a list of points in the form of x-y pairs: (x, y)
2. triangles become a list of point indexes (0, 1, 2)
3. triangle list becomes a list of three-tuples, each of which defines a triangle with three point indexes

Current implementation of boundary condition does not include particles re-entering material (e.g., for castellated tiles)

Current implementation does not appropriately track particles outside of the energy barrier but inside the simulation boundary - currently, as soon as they cross the energy barrier, they are taken out of the simulation. There needs to be a distinction between the energy barrier and the simulation boundary if macroscopic geometries are to be simulated.

Add multicomponent targets

For homogeneous targets, this should be as simple as completing the Z_eff, number_density, m_eff, Eb, and choose placeholder methods in Material, extending the Material fields n, m, Z, Eb, Es, Ec, ... to Vec, and updating the I/O.

[question] Implementing New Lennard-Jones Interactions

Currently RustBCA has implemented Lennard-Jones potentials with 12-6 and 6.5-6 exponentials. When adding a new potential within interactions.rs, most functions are fairly straightforward, except for the doca_lennard_jones function. Is the generalized form below the correct way to add this function?

4*ε*( (σ/r).powf(m) - (σ/r).powf(n) )
(r/σ).powf(m) - 4*ε/relative_energy*( 1-(r/σ).powf(m-n) ) - (p/σ).powf(2)*(r/σ).powf(m-2)

Thanks for your help!

Multithreading version uses too much memory

Need a solution to handle streaming of chunks to/from disk for very large simulations when using multithreading. Currently, the code can handle input of up to potentially 100M particles, but the output particle struct will be very large and potentially take up too much memory. If rustbca is refactored into purely a library, this is less of an issue, but a standalone version needs a solution.

Order of BCA operations is not correct

Current order of BCA operations is as follows:

Ion start at x0, y0, z0 with direction u0, v0, w0
Determine mfp, impact parameter, azimuthal angle
Choose collision partner at position x_recoil, y_recoil, z_recoil
If x_recoil, y_recoil inside 2D geometry {
Calculate distance of closest approach
Calculate lab-frame scattering angle for ion and recoil
Calculate asymptotic deflection
Rotate ion by lab-frame scattering angle to u1, v1, w1
Rotate recoil by lab-frame recoil angle
Advance ion to x0 + mfp*u1, y0 + mfp*v1, y0 + mfp*v1
If recoil energy greater than cutoff {
Add recoil to particle list
}
}

However, the ion's next position should be:
(x0 + mfp*u0, y0 + mfp*v0, z0 + mfp*w0)

So the order of operations is incorrect. Most succinct solution which preserves the structure of the weak collision loop would be to simply store the particle's previous direction and use that to update the particle location, much as it stores its previous position to determine when it crosses boundaries.

Python toml is unhappy with array of inline tables for multiple potentials

For multiple interaction potentials, I'm using an NXN array of inline tables that is deserialized into a nested NXN Vec of Rootfinder, ScatteringIntegral, and InteractionPotential struct enum variants.

However, arrays of inline tables appear to not actually be supported by the TOML spec - it works in toml-rs, and I can read the files into the Python toml library, but the Atom toml linter is saying it's an error, and an open issue on the Python toml repository suggests that inline tables can't go into arrays in the official TOML spec.

Potential solutions:

• Replace toml with another config file format (yaml, json, ...)
• Use the workaround in rustbca.py - that is, manually append troublesome options to the bottom of the input file

Trajectory tracking doubles up on initial location of ions

[feature] 3D geometries, including mathematically- and mesh-defined

Is your feature request related to a problem? Please describe.
Moving to 3D geometries will enable the study of some unique problems not
feasible in 2D.

Proposed solution

1. Abstract away the input to simplify multiple wrappers for single_ion_bca(). See issue #79.
2. Abstract away the Material using Type Generics so that the traits necessary (distance to surface, local composition, whether or not a point is inside/outside the material). See issue #77
3. Add specific wrapper for mathematically defined 3D geometries
4. Add specific wrapper for mesh-defined 3D geometries

[material_parameters] n is never used

Multiple interaction potentials for different particle pairs

Now that the CPR and polynomial root-finders are working in the cpr_rootfinder branch, the time has come to move away from universal interatomic potentials. Potentials should be specifiable on a species by species basis.

The simplest way to implement this would be the following:

Each species gets a unique integer tag

These tags correspond to indices on an interaction matrix. So for example, if the system is He-W:

He: 0
W:  1
[[ He-He, He-W],
[W-He, W-W]]

The values of the interaction matrix would be integer tags that each correspond to an interaction potential. Additional parameters could be vectors in a same-shaped matrix (e.g., sigma and epsilon for Lennard-Jones).

Implement Bremmstrahlung production / Kinchin-Pease / On-the-fly physics calculations

[bug] progress bar updates with particle chunks instead of per particle

Description
The progress bar updates too slowly to be useful; it should track per particle instead of per chunk.

Setting Z=0 (technically bad input) produces untracked NaNs

Some calculations that may produce NaNs are not checked appropriately - in the above example, the surface refraction code is the first time the NaNs are caught from bad input. I believe the NaNs are produced in the nonlocal electronic stopping section. More robust NaN handling should be added on calculations that, depending on user input, could fail in this way.

Add documentation of Mesh2D object

This issue was brought to my attention by a user - the way that coordinate_sets is used to build the material mesh is not clear. Both changes to nomenclature and additions to the documentation are necessary.

Need test coverage of all functions

[feature]

Remove or suppress warning about metadata from intel-mkl-src in version requirement.

[feature] Python script to create Mesh2d for Input File

Is your feature request related to a problem? Please describe.
For higher precision the geometry will be required to have a large number of points (>100). Creating triangles from a geometry that could include multiple complex geometries is extremely time consuming and not easy for a researcher.

Proposed solution
Add a python script to create Mesh2ds using simple methods that anyone could understand.

What length_unit does energy_barrier_thickness use?

MaterialParameters doesn't include a length unit... how is energy_barrier_thickness calculated from input?

Functions should not rely on tuple arguments/returns

Should create structs for holding function output; e.g., BinaryCollisionResult which would contain doca, theta, psi, etc.