[bug] Yields differ for surfaces with x directed normal and z directed normal

[dmonah]

Description
Calculated yield results are different for surfaces with x directed normal and z directed normal.

I noticed this when comparing 2 models using the parry3d option. If I orient the surface with a +x normal and bombard with -x directed ions I get one value for yield. If I reorientated the 3d model so that I have a +z normal and bombard with -z directed ions, my yield drops by about 5% to 25%.

I hacking the 1d code in the Mesh0D Geometry implementation (modifying the methods: inside, inside_simulation_boundary, inside_energy_barrier, and closest_point) so that I can select the normal direction, I see the exact same behavior in the 1d code with the yield values predicted exactly matching the equivalent 3d cases.

To Reproduce
I can't upload any toml files, images or patches as uploading is blocked by my work firewall. Therefore to reproduce, you will need to reproduce cases from the description, sorry. I can share the hack I made to the Mesh0D class to reproduce in 1d:

fn inside(&self, x: f64, y: f64, z: f64) -> bool {
        let v: f64;

        #[cfg(not(feature = "mesh0d_z_normal"))]
        {
            v = x;
        }

        #[cfg(feature = "mesh0d_z_normal")]
        {
            v = z;
        }

        #[cfg(not(feature = "mesh0d_positive_normal"))]
        {
            v > 0.0
        }

        #[cfg(feature = "mesh0d_positive_normal")]
        {
            v < 0.0
        }
    }

    fn inside_simulation_boundary(&self, x: f64, y: f64, z: f64) -> bool {
        let v: f64;

        #[cfg(not(feature = "mesh0d_z_normal"))]
        {
            v = x;
        }

        #[cfg(feature = "mesh0d_z_normal")]
        {
            v = z;
        }

        #[cfg(not(feature = "mesh0d_positive_normal"))]
        {
            v > -10.*self.energy_barrier_thickness
        }

        #[cfg(feature = "mesh0d_positive_normal")]
        {
            v < 10.*self.energy_barrier_thickness
        }
    }

    fn inside_energy_barrier(&self, x: f64, y: f64, z: f64) -> bool {
        let v: f64;

        #[cfg(not(feature = "mesh0d_z_normal"))]
        {
            v = x;
        }

        #[cfg(feature = "mesh0d_z_normal")]
        {
            v = z;
        }

        #[cfg(not(feature = "mesh0d_positive_normal"))]
        {
            v > -self.energy_barrier_thickness
        }

        #[cfg(feature = "mesh0d_positive_normal")]
        {
            v < self.energy_barrier_thickness
        }
    }

    fn closest_point(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) {
        #[cfg(not(feature = "mesh0d_z_normal"))]
        {
            (0., y, z)
        }

        #[cfg(feature = "mesh0d_z_normal")]
        {
            (x, y, 0.)
        }
    }
}

Expected behavior
Yield should be independent of normal direction in 3d (yes?).

Error messages, output files, or figures
Figures attached.

System (please complete the following information):
OS:

Additional context
I've chased the source of the issue down to the following code in choose_collision_partner function:

    //Find recoil location
    let x_recoil: f64 = x + mfp*cosx - impact_parameter*cosphi*sinx;
    let y_recoil: f64 = y + mfp*cosy - impact_parameter*(sinphi*cosz - cosphi*cosy*cosx)/sinx;
    let z_recoil: f64 = z + mfp*cosz + impact_parameter*(sinphi*cosy - cosphi*cosx*cosz)/sinx;

This code attaches significance to the x direction and does not treat all 3 directions equally. I tried applying a fix which does treat all axes the same, and it removes the asymmetry (see bug fix figure), but i do not get expected behavior at glancing angles (yield approaching 0), so I think the specifics of my fix are incorrect.

[drobnyjt]

@dmonah Thank you very much for looking into this; I will investigate further. I encountered some strangeness previously with incident directions of y and z, but believed those had been resolved by some changes to how RustBCA handles the first collision. Unfortunately it seems those fixes were not sufficient.

Can you clarify for your example cases, are you sending ions in exactly the +/- x or z directions, or with some small deviation? My first suspect would be a singularity in the recoil generation routine that you point out (which comes from TRIM originally). I also suspect that this is likely a first collision effect, since the details of the first collision tend to dominate the sputtering yield values. If this is the case, the fix might be straightforward.

[drobnyjt]

Probably the easiest test case to spin up would be a sphere where each direction can be tested all in one input file.

[dmonah]

Hi, Thanks so much for jumping on this so quickly! My minimum incident angle is 0.001 degrees to avoid the singularity, as recommended in the docs. However, I run a number of anlgles and see that the yield calculated with a +/- z normal is consistently lower than the yield calculated with +/- x normal. This is not surprising when I look at the source code as there is a preferential direction in the code identified in the bug report. Let me see if I can attach some figures to this email. Default behavior: [cid:21fa7d26-9180-453b-b54e-c58b32e82697] Modified "choose_collision_partner" code: [cid:0c99a8aa-4be4-41a9-86d1-74b9f97a07ee] If you see figures attached to this reply, I can try sending some toml cases. Derek Confidential – Limited Access and Use

…

________________________________ From: Jon Drobny ***@***.***> Sent: Tuesday, July 8, 2025 14:28 To: lcpp-org/RustBCA ***@***.***> Cc: Monahan, Derek ***@***.***>; Mention ***@***.***> Subject: Re: [lcpp-org/RustBCA] [bug] Yields differ for surfaces with x directed normal and z directed normal (Issue #286) You don't often get email from ***@***.*** Learn why this is important<https://aka.ms/LearnAboutSenderIdentification> External Email: Do NOT reply, click on links, or open attachments unless you recognize the sender and know the content is safe. If you believe this email may be unsafe, please click on the “Report Phishing” button on the top right of Outlook. [https://avatars.githubusercontent.com/u/37962344?s=20&v=4]drobnyjt left a comment (lcpp-org/RustBCA#286)<#286 (comment)> @dmonah<https://github.com/dmonah> Thank you very much for looking into this; I will investigate further. I encountered some strangeness previously with incident directions of y and z, but believed those had been resolved by some changes to how RustBCA handles the first collision. Can you clarify for your example cases, are you sending ions in exactly the +/- x or z directions, or with some small deviation? My first suspect would be a singularity in the recoil generation routine that you point out (which comes from TRIM originally). I also suspect that this is likely a first collision effect, since the details of the first collision tend to dominate the sputtering yield values. If this is the case, the fix might be straightforward. — Reply to this email directly, view it on GitHub<#286 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/BOKFAK3XK3VDYOT4EMHSMQL3HPBP3AVCNFSM6AAAAACA6YFCACVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZTANBYHE3TIOBZGI>. You are receiving this because you were mentioned. LAM RESEARCH CONFIDENTIALITY NOTICE: This e-mail transmission, and any documents, files, or previous e-mail messages attached to it, (collectively, "E-mail Transmission") may be subject to one or more of the following based on the associated sensitivity level: E-mail Transmission (i) contains confidential information, (ii) is prohibited from distribution outside of Lam, and/or (iii) is intended solely for and restricted to the specified recipient(s). If you are not the intended recipient, or a person responsible for delivering it to the intended recipient, you are hereby notified that any disclosure, copying, distribution or use of any of the information contained in or attached to this message is STRICTLY PROHIBITED. If you have received this transmission in error, please immediately notify the sender and destroy the original transmission and its attachments without reading them or saving them to disk. Thank you.

[drobnyjt]

Hi @dmonah - thanks for the reply; unfortunately the figures did not attach correctly.

I should have time to dig deeply into this tomorrow, including finding the original derivation of those formulas and running some tests.

Without the figures it's not clear what effects you're seeing - my previous understanding was that the apparent preference for the x-direction in those three formulas was due to the choice of a particular axis for a 180-degree rotation and that it should be arbitrary with no effect on the results except for which direction(s) the singularity exists, but that could very well be wrong, or there could be a mistake in the formula or my implementation of it.

If it were totally wrong, I would expect it to have been noticed already in the angular distributions of emitted particles at oblique exit angles, but those results so far have compared favorably to previous codes.

[dmonah]

Again, thank you for being so quick to jump on this. Sorry I can't get any figures, or source code to you. I guess I can manually paste TOML into an email if you struggle to reproduce on your end. Thanks, Derek

[drobnyjt]

@dmonah some updates on this - first, I am able to reproduce the effect (albeit smaller than the 5-25% reported at ~1% in the reflection coefficients - it's possible the sputtering yields will be worse, but they would take a lot longer to run) with this input file using the SPHERE geometry type:

#Use with SPHERE geometry option only
[options]
name = "boron_nitride_"
track_trajectories = false
track_recoils = false
track_recoil_trajectories = false
weak_collision_order = 0
num_threads = 1
num_chunks = 1

[material_parameters]
energy_unit = "EV"
mass_unit = "AMU"
Eb = [ 0.0, 0.0,]
Es = [ 0.0, 0.0,]
Ec = [ 5.76, 1.0,]
Z = [ 5, 7,]
m = [ 10.811, 14.0,]
interaction_index = [0, 0]
surface_binding_model = {"ISOTROPIC"={calculation="TARGET"}}
bulk_binding_model = "AVERAGE"

[particle_parameters]
length_unit = "ANGSTROM"
energy_unit = "EV"
mass_unit = "AMU"
N = [ 1000000, 1000000, 1000000, 1000000 ]
m = [ 1.008, 1.009, 1.010, 1.011 ]
Z = [ 1, 1, 1, 1 ]
E = [ 50.0, 50.0, 50.0, 50.0 ]
Ec = [ 1.0, 1.0, 1.0, 1.0 ]
Es = [ 0.0, 0.0, 0.0, 0.0 ]
pos = [ [ -100.0, 0.0, 0.0,], [ 100.0, 0.0, 0.0,] , [ 0.0, 0.0, 100.0,], [ 0.0, 0.0, -100.0,]]
dir = [ [ 0.999, 0.001, 0.0,], [-0.999, 0.001, 0.0], [0.0, 0.001, -0.999], [0.0, 0.0, 0.999 ] ]
interaction_index = [ 0, 0, 0, 0 ]
particle_input_filename=""

[geometry_input]
length_unit = "ANGSTROM"
radius = 100.0
electronic_stopping_correction_factor = 1.0
densities = [0.065, 0.065]

Notably, the effect persists when surface binding is set to isotropic and Es is set to zero, so it's not a surface binding effect.

When investigating the formula you highlighted, I am not able to reproduce a significant difference between x, y, and z directions when sweeping through impact parameters and azimuthal angles and projecting the result onto the local y-z plane:

EDIT: actually, I misrepresented the impact parameter sweep - I'm sweeping through the functional form that samples from p-space, which is sqrt(-ln(R))*pmax where R is (0, 1).

although I admit there may be differences I am not detecting with this test.

For next steps, I am currently working on rederiving the formulas you highlighted from scratch to see if simply replacing them with an alternate version (that should be numerically equivalent) and will see if replacing them fixes the issue.

My next suspicion would be in the into/out of surface routines - there have been issues with these in complex geometries before, and there's a long-standing issue to make it more robust. It is possible that in certain orientations some collisions near the surface are being erroneously skipped.

[drobnyjt]

I've derived some new formulas for the recoil generation:

x_r = x0 + mfp*cosx + p*(cosz*cosphi+ cosy*sinphi)
y_r = y0 + mfp*cosy + (p*(-cosy*cosz*cosphi+(1 + cosx-cosy**2)*sinphi))/(1 + cosx)
z_r = z0 + mfp*cosz + (p*((1 + cosx - cosz**2)*cosphi- cosy*cosz*sinphi))/(1 + cosx)

Unfortunately I don't have the original derivation from TRIM on hand, which appears to yield simpler formulas. I found these by first finding a rotation matrix R that corresponds to transformations from r1 = (1, 0, 0), the local atom's frame, to r1' = (rx, ry, rz), the lab frame (using this method which I've used previously in RustBCA):

R = [
[ rx, ry, rz ],
[ -ry, 1-ry^2/(1+rx), -((ry rz)/(1+rx))],
[ -rz, -((ry rz)/(1+rx)), 1-rz^2/(1+rx) ],
]

and transforming the local recoil location vector, r2= p*(0, cos(phi), sin(phi)), into r2' = R.r2:

r2 = [p (rz Cos[phi] + ry Sin[phi]), (
 p (-ry rz Cos[phi] + (1 + rx - ry^2) Sin[phi]))/(1 + rx), (
 p ((1 + rx - rz^2) Cos[phi] - ry rz Sin[phi]))/(1 + rx)]

Then it's merely a matter of adding the lab frame vector onto the lab frame location of the next collision:

r_recoil = r0' + mfp*r1' + r2'

which yields the formulas above (substituting cosx for rx etc.).

I do need to check that this version of the local recoil location vector r2 maintains a consistent direction for phi compared to RustBCA, but in practice it doesn't matter because phi is randomly picked from a uniform distribution.

Edit: it technically does matter, because if phi is inconsistent than the generated recoil will be in a different azimuthal position than the deflected/deflecting directions post-collision, so this is a must-check even if it doesn't matter for this test.

They're a little less efficient, and obviously derived a different way, but the results of the linear sweep are the same:

After a little more testing I'll pop these into RustBCA and see if it resolves the observed issue with non-x incident directions - if not, I'll move to looking at the inside/out of surface routines next.

[drobnyjt]

@dmonah would you be able to provide some more detail on the exact debugging process that led you to conclude that those formulas were the source of the error? If you can share the input files that used parry3d that would be ideal, particularly one where you're seeing 25% differences in the yield.

Alternatively, you could try swapping out the current recoil generation formulas with the newly derived ones and reporting back if that resolves the discrepancy.

I thought I had been able to reproduce the effect locally, but the more I test the less convinced I am that I'm observing a significant effect. At most I can get 1-3% and it's possible that's just due to floating point error around the singularity and not reflective of an actual mistake in the implementation.

[dmonah]

Okay, so I still can't upload any files / images so I will copy some TOML files:

1d example:

[options]
name = ""
track_trajectories = false
track_recoils = true
track_recoil_trajectories = false
track_displacements = false
track_energy_losses = false
write_buffer_size = 8000
weak_collision_order = 0
suppress_deep_recoils = false
high_energy_free_flight_paths = false
num_threads = 8
num_chunks = 10
use_hdf5 = false
electronic_stopping_mode = "LOW_ENERGY_NONLOCAL"
mean_free_path_model = "LIQUID"
interaction_potential = [ [ "KR_C",],]
scattering_integral = [ [ "MENDENHALL_WELLER",],]
energy_min = 0.0
energy_max = 1000.0
energy_num = 100
angle_min = 0.0
angle_max = 90.0
angle_num = 45
x_min = 0.0
x_max = 100.0
x_num = 50
y_min = 0.0
y_max = 100.0
y_num = 20
z_min = -1000.0
z_max = 0.0
z_num = 50

[particle_parameters]
length_unit = "ANGSTROM"
energy_unit = "EV"
mass_unit = "AMU"
N = [ 25000,]
m = [ 39.948,]
Z = [ 18.0,]
E = [ 1000.0,]
Ec = [ 0.5,]
Es = [ 0.0,]
interaction_index = [ 0,]
pos = [ [ -10.0, 0.0, 0.0,],]
dir = [ [ 0.9999999998476913, 1.7453292519057202e-5, 0.0,],]

[material_parameters]
energy_unit = "EV"
mass_unit = "AMU"
Eb = [ 2.0,]
Es = [ 4.72,]
Ec = [ 1.573,]
Ed = [ 0,]
Z = [ 14.0,]
m = [ 28.08553,]
interaction_index = [ 0,]
surface_binding_model = "INDIVIDUAL"
bulk_binding_model = "AVERAGE"

[geometry_input]
length_unit = "ANGSTROM"
electronic_stopping_correction_factor = 1.0
densities = [ 0.04994,]

Before running, I make the code changes to the Mesh0D class shown in the first post. Then

Run 1:

• Run command: cargo run --release 0D input.toml
• Repeat for incident angles 50 degrees (dir = [ [ 0.6427876096865394, 0.766044443118978, 0.0],]), 60 degrees (dir = [ [ 0.5000000000000001, 0.8660254037844386, 0.0],]), 70 degrees (dir = [ [ 0.3420201433256688, 0.9396926207859083, 0.0],]), 85 degrees (dir = [ [ 0.08715574274765812, 0.9961946980917455, 0.0,],])

Run 2:

• Change input: dir = [ [ 0.0, 1.7453292519057202e-5, 0.9999999998476913,],]
• Run command: cargo run --release --features mesh0d_z_normal 0D input.toml
• Repeat for incident angles 50 degrees (dir = [ [ 0.0, 0.766044443118978, 0.6427876096865394],]), 60 degrees (dir = [ [ 0.0, 0.8660254037844386, 0.5000000000000001],]), 70 degrees (dir = [ [ 0.0, 0.9396926207859083, 0.3420201433256688],]), 85 degrees (dir = [ [ 0.0, 0.9961946980917455, 0.08715574274765812,],])

Output:

• For all incident angles, sputter yields calculated in Run 1 is > the sputter yield calculated in Run 2 (% difference greatest at 0 degrees incident).
• Reflection coefficients to not differ significantly.

[dmonah]

My parry 3d example needs no code modification, but due to the large number of points, I can't paste the mesh here (page crashes). Here is the TOML with mesh data removed:

[options]
name = ""
track_trajectories = false
track_recoils = true
track_recoil_trajectories = false
track_displacements = false
track_energy_losses = false
write_buffer_size = 8000
weak_collision_order = 0
suppress_deep_recoils = false
high_energy_free_flight_paths = false
num_threads = 8
num_chunks = 10
use_hdf5 = false
electronic_stopping_mode = "LOW_ENERGY_NONLOCAL"
mean_free_path_model = "LIQUID"
interaction_potential = [ [ "KR_C",],]
scattering_integral = [ [ "MENDENHALL_WELLER",],]
energy_min = 0.0
energy_max = 1000.0
energy_num = 100
angle_min = 0.0
angle_max = 90.0
angle_num = 45
x_min = 0.0
x_max = 100.0
x_num = 50
y_min = 0.0
y_max = 100.0
y_num = 20
z_min = -1000.0
z_max = 0.0
z_num = 50

[particle_parameters]
length_unit = "ANGSTROM"
energy_unit = "EV"
mass_unit = "AMU"
N = [ 25000,]
m = [ 39.948,]
Z = [ 18.0,]
E = [ 1000.0,]
Ec = [ 0.5,]
Es = [ 0.0,]
interaction_index = [ 0,]
pos = [ [ 10.0, 0.0, 0.0,],]
dir = [ [ -0.9999999998476914, 1.7453292519057202e-5, 0.0,],]

[material_parameters]
energy_unit = "EV"
mass_unit = "AMU"
Eb = [ 2.0,]
Es = [ 4.72,]
Ec = [ 1.573,]
Ed = [ 0,]
Z = [ 14.0,]
m = [ 28.08553,]
interaction_index = [ 0,]
surface_binding_model = "INDIVIDUAL"
bulk_binding_model = "AVERAGE"

[geometry_input]
length_unit = "ANGSTROM"
electronic_stopping_correction_factor = 1.0
densities = [ 0.04994,]
vertices = [ ... ]
indices = [ ... ]

To run, you can create your own meshes which are large enough to ensure no recoil particles cross the boundaries and have a flat surface with:

Run 1:

• normal: (1, 0, 0)
• plane: x = 0

Run 2:

• normal: (0, 0, 1)
• plane: z = 0
• pos = [ [ 0.0, 0.0, 10.0,],]
• dir = [ [ 0.0, 1.7453292519057202e-5, -0.9999999998476914,],]

Don't bother with different incident angles here as it only adds to the issue due to a lack of periodic boundaries.

Output:

• Again, sputter yields calculated in Run 1 is > the sputter yield calculated in Run 2.
• Sputter yields from Run 1 match those calculated in Run 1 for the Mesh0D case and sputter yield from Run 2 match those calculated in Run 2 for the Mesh 0D case.

[dmonah]

Finally, if I edit the choose_collision_partner so that the 3rd term is ignored, e.g.:

        x_recoil = x + mfp*cosx;
        y_recoil = y + mfp*cosy;
        z_recoil = z + mfp*cosz;

All results 1d, 3d, z directed normal, x directed normal match almost exactly.

[dmonah]

I suspected that the issue might be related to the particles being rotated after the collision and not before, resulting the z directed normal particles incorrectly calculating the collision partner location on the first collision. However, I'm not so sure that is the case now as I can set dir[0] == 0 for my initial direction and I don't get any issues.