lcpp-org / RustBCA

A free, open-source Binary Collision Approximation (BCA) code for ion-material interactions including sputtering, implantation, and reflection
https://github.com/lcpp-org/RustBCA/wiki
GNU General Public License v3.0
42 stars 15 forks source link

[question] apparent rarely occurring infinite loop at high energies for certain CPR rootfinder parameters with the combined Kr-C-Morse potential #224

Open drobnyjt opened 1 year ago

drobnyjt commented 1 year ago

There seems to be a rare event at high energies (>10 keV) for the Kr-C-Morse combined potential that causes an infinite loop. I suspect a value is going to +/-INF and is not being caught by a NaN check, but more investigation is required.

The occurrence rate is somewhere between 1/250000 ions and 1/10000 ions, but any numerical situation that causes an infinite loop should instead panic. The question is: 1) where is the issue occurring? 2) can it be ameliorated with changes to the CPR rootfinder parameters? 3) what changes should be made to the code?

Example input file that exhibits this issue:

    [options]
    name = "krc_morse_23_default_Es_1_5"
    track_recoils = false
    weak_collision_order = 0
    electronic_stopping_mode = "LOW_ENERGY_NONLOCAL"
    mean_free_path_model = "LIQUID"
    interaction_potential = [[{"KRC_MORSE"={D=6.758838e-20, r0=2.7820000000000003e-10, alpha=14197999999.999998, k=70000000000.0, x0=7.5e-11}}]]
    scattering_integral = [["GAUSS_LEGENDRE"]]
    root_finder = [[{"CPR"={n0=2, nmax=100, epsilon=1E-7, complex_threshold=1E-6, truncation_threshold=1E-9, far_from_zero=1E9, interval_limit=1E-12, derivative_free=true}}]]
    num_threads = 4
    num_chunks = 1

    [particle_parameters]
    length_unit = "ANGSTROM"
    energy_unit = "EV"
    mass_unit = "AMU"
    N = [ 10000 ]
    m = [ 1.008 ]
    Z = [ 1 ]
    E = [ 6189.65818891261 ]
    Ec = [ 0.1 ]
    Es = [ 1.5 ]
    interaction_index = [ 0 ]
    pos = [ [ -4.4, 0.0, 0.0,] ]
    dir = [ [ 0.9999999999984769, 1.7453292519934434e-6, 0.0,] ]

    [geometry_input]
    length_unit = "ANGSTROM"
    electronic_stopping_correction_factor = 1.09
    densities = [ 0.0914 ]

    [material_parameters]
    energy_unit = "EV"
    mass_unit = "AMU"
    Eb = [ 0.0 ]
    Es = [ 5.61 ]
    Ec = [ 3.0 ]
    Z = [ 28 ]
    m = [ 58.69 ]
    interaction_index = [ 0 ]
    surface_binding_model = {"PLANAR"={calculation="INDIVIDUAL"}}
    bulk_binding_model = "AVERAGE"
drobnyjt commented 1 year ago

Potential cause: if the CPR rootfinder truncation threshold is too high, there may be an infinite loop caused by truncated terms leading to identical Chebyshev polynomials each iteration and the error not decreasing as n increases, leading to infinite subdivision - although it should hit either nmax or the subdivision interval limit first. To test this, reducing those two parameters should be tried to see if it's not actually an infinite loop, just a very, very long loop caused by many large matrix inversions.

drobnyjt commented 1 year ago

Interestingly, I cannot reproduce this on a different machine. More investigation is needed.

drobnyjt commented 10 months ago

On the original machine with this issue, it appears to be limited to high energy, low nmax, and very small interval limits.