[question] Grazing incidence angle with 3D TRIMESH option

[Torben-S]

Hi,

I plan on using RustBCA for BCA simulations for rough surfaces. Doing some tests with the code I noticed, that I get very different results for ion impact on a flat surface at grazing incident angles for the 0D and TRIMESH options.

For my trimesh I just created a rectangle, using coordinates that are large enough, such that nothing leaves the surface anywhere but at the top (or is deposited). Its from -50nm to +50nm in x- and y-direction and from 0nm to -100nm in z direction. Particles hit at x,y,z=0,0,0.

Doing that I get very unusual sputter yields for 200 eV Ne impacting a W surface at grazing incidence angles that are fundamentally different from results I got from RustBCA in 0D and also 3D-SDTrimSP. The yields and reflection coefficients from the "baseline" 0D RustBCA and 3D-SDTrimSP simulations look like this

When using my Mesh with 3D RustBCA the yield diverges after 40°, and even increases for the grazing incidence.

At angles of incidence of $\geq87^\circ$ I would expect (and the 0D RustBCA and the 3D-SDTrimSP simulations back this) almost only reflection and no sputtering.

To check what the issue might be, I did a simulation with 0D and 3D where I tracked the full trajectories for 10 incoming ions

In the bottom figures there are the trajectories from the 0D case, all Ne ions are reflected forward, which is expected.
For the 3D case on the top we can see that the particles enter the surface instead of being reflected immediately on the top. Also the Ne ions cause some sputtering, because of the collision cascades within the surface.

Do you have any idea why that might be? In the input file all the parameters are the same. Except for the coordinate system used, because my mesh starts at z=0 and goes to z=-50, so the ions are initiated with
pos = [ [ -19.081173196775577, 0, 1,] ]
dir = [ [ 0.9986295399881645, 0, -0.05233585637999017,] ]
and for the 0D case with the surface at x=0 in positive x-direction with
pos = [ [ -1, -19.081173196775577, 0,] ]
dir = [ [ 0.05233585637999017, 0.9986295399881645, 0,] ]

I tested different resolutions of the mesh, but none worked, and if I understood the workings of RustBCA correctly, for this flat case the mesh resolution shouldn't matter, right?

I would appreciate any help, because I would really like to use RustBCA for my 3D BCA work. If you have any further question what options I used, let me know.

Cheers, Torben

[drobnyjt]

Hey @Torben-S; thanks for making an issue, and including so much information. This is definitely unexpected behavior.

It is possible this is related to Issue #286, in which case I should be able to diagnose the issue fairly easily.

Can you share an example input file with the trimesh, or just the trimesh section of the input file? That'll help me figure out what's going on.

[drobnyjt]

test_trimesh_vs_0d.py

Doing Neon on Tungsten with the script linked above for 200 eV and angles between 0 and 89 degrees, I get this behavior for 0D versus TRIMESH, as expected:

So something is going on in your simulations, we just need to figure out what.

Oh, also, if you're able to share SDTrimSP input files as well I would appreciate it - I don't have access to the code but I am familiar with the various options and ideally (and in much of my prior testing against TRIDYN) the two codes should be a little bit more similar than in your first two plots if the input options correspond exactly.

[Torben-S]

Thanks for answering so quickly.

I tried your script and it worked, I get the same reasonable result you do.

After that, I input your geometry and definition for particle position and direction into my own toml. it also worked fine.

My own geometry input looks like this. I switched from nm to Angstrom, to be closer to your input. Also I dont use the A but 1e-10, because I calculate the density according to this value.

length_unit = "1.0e-10"
electronic_stopping_correction_factor = 1.0
densities = [ 0.06306, ]
vertices = [[-500.0, -500.0, 0.0], [500.0, -500.0, 0.0], [-500.0, 500.0, 0.0], [500.0, 500.0, 0.0], [-500.0, -500.0, -500.0], [500.0, -500.0, -500.0], [-500.0, 500.0, -500.0], [500.0, 500.0, -500.0]]
indices = [[0, 1, 2], [1, 3, 2], [6, 5, 4], [6, 7, 5], [4, 1, 0], [4, 5, 1], [5, 3, 1], [5, 7, 3], [7, 2, 3], [7, 6, 2], [6, 0, 2], [6, 4, 0]]

It seems to me the only remaining difference is the geometry and coordinate system. For the particles I use
pos = [[-float(np.sin(alpha))/float(np.cos(alpha)),0,1]]
dir = [[float(np.sin(alpha)),0,-float(np.cos(alpha))]]

I could try to rewrite my geometry, so that the particles fly in the positive x-direction, but it should work independent of the direction right? Or do you see another problem with my geometry?

My full input toml is the following, I think we used similar settings:
#Use with TRIMESH geometry option only
[options]
name = ""
track_trajectories = false
track_recoils = true
track_recoil_trajectories = false
write_buffer_size = 8000
weak_collision_order = 0
suppress_deep_recoils = false
high_energy_free_flight_paths = false
num_threads = 10
num_chunks = 50
use_hdf5 = false
electronic_stopping_mode = "LOW_ENERGY_NONLOCAL"
mean_free_path_model = "LIQUID"
interaction_potential = [["KR_C"]]
scattering_integral = [['MAGIC']]
root_finder = [[{"NEWTON"={max_iterations=100, tolerance=1E-3}}]]
track_displacements = false
track_energy_losses = false

[material_parameters]
energy_unit = "EV"
mass_unit = "AMU"
Eb = [ 0.0, ]
Es = [ 8.79, ]
Ec = [ 2.93, ]
Z = [ 74.0, ]
m = [ 183.84, ]
interaction_index = [ 0, ]
surface_binding_model = {"PLANAR"={calculation="INDIVIDUAL"}}
bulk_binding_model = "INDIVIDUAL"

[particle_parameters]
length_unit = "1.0e-10"
energy_unit = "EV"
mass_unit = "AMU"
N = [ 5000, ]
m = [ 20.1797, ]
Z = [ 10.0, ]
E = [ 200, ]
Ec = [ 0.5, ]
Es = [ 0.0, ]
pos = [ [ -5.6712851359633385, 0, 1,] ]
dir = [ [ 0.9848077703770209, 0, -0.1736480791861542,] ]
interaction_index = [0]
particle_input_filename = "None"

[geometry_input]
length_unit = "1.0e-10"
electronic_stopping_correction_factor = 1.0
densities = [ 0.06306, ]
vertices = [[-500.0, -500.0, 0.0], [500.0, -500.0, 0.0], [-500.0, 500.0, 0.0], [500.0, 500.0, 0.0], [-500.0, -500.0, -500.0], [500.0, -500.0, -500.0], [-500.0, 500.0, -500.0], [500.0, 500.0, -500.0]]
indices = [[0, 1, 2], [1, 3, 2], [6, 5, 4], [6, 7, 5], [4, 1, 0], [4, 5, 1], [5, 3, 1], [5, 7, 3], [7, 2, 3], [7, 6, 2], [6, 0, 2], [6, 4, 0]]

[drobnyjt]

@Torben-S Thanks! Yeah, it should definitely work for arbitrary directions - that was the problem in the issue mentioned above; that should have been resolved, but it is possible there are remaining edge cases I did not account for.

However, I am actually having trouble replicating the problem you observed on my end - I've added this block to the test script:

    if not run_3d:
        x0 = -5.0
        y0 = 0.0
        z0 = 1.0
        ux = np.cos(angle*np.pi/180.)
        uy = np.sin(angle*np.pi/180.)
        uz = 0.

    else:
        x0 = -float(np.sin(angle*np.pi/180.))/float(np.cos(angle*np.pi/180.))
        y0 = 0.0
        z0 = 1.0
        print(x0, y0, z0)
        ux = np.sin(angle*np.pi/180.)
        uy = 0.0
        uz = -np.cos(angle*np.pi/180.)

I thought it was possible that the initial position you were using could be the problem (I usually just initiate atoms at ~one interatomic distance from the surface - if I recall correctly, this was the recommendation from the original TRIM.SP paper and seems to give good results) but that isn't causing any issues as far as I can tell. If you got too far away though, it might cause issues, but that is not happening here.

I also replaced the entire input file in my test script with the one you shared, but am still getting this result:

which looks okay...

Here's the new script using your input file. I'll keep fiddling with it and see if there's something I'm missing. Just to check - are you using the most recent version of rustbca? Are you using the version from the main or dev branches?

test_trimesh_vs_0d.py

[drobnyjt]

@Torben-S I suspect I've figured it out - it looks like the fixes for issue #286 were never merged from dev into main - if I revert to main and run the python script, I get weird results:

So this is probably a sign that that merge is very much overdue. I'll start that process now.

[Torben-S]

I see, yeah I also tried switching to the dev branch now, there it works fine, even with my script.

Thank you very much once again for your swift help, I will continue on the main branch once the fix is included. I hope my next steps will not cause similar problems, but it's overall a nice code to work with.

I can send you my SDTrimSP input some time tomorrow, but will check once more whether all the energies etc. really match, if you say the yields should be closer together.

[drobnyjt]

Great, glad to hear it, you're very welcome. I've just merged dev onto main and published a new release (v2.12.3) with all of the fixes.

I would also be suspicious of the scattering integral - during my dissertation, I discovered that the coefficients for the MAGIC algorithm in TRIDYN (and presumably TRIM.SP) produced incorrect results for Kr-C compared to quadrature, and I do not know if that was ever fixed for SDTrimSP - SDTrimSP should have quadrature available though, and that will be the best point of comparison to RustBCA.

This figure from my dissertation (available here) shows the effect of the incorrect MAGIC coefficients for Kr-C in TRIDYN on scattering angles, but again, I do not know if the coefficients in SDTrimSP were fixed or not, I think the version of TRIDYN I had access to in grad school dated to the 90s originally.