quantise_diatomic fails for some positions/velocities

[Alexsp32]

When using InitialConditions.QuantisedDiatomic.quantise_diatomic to determine vibrational/rotational quanta of diatomic species, the following error will occur for some combinations of positions and velocities:

ERROR: LoadError: ArgumentError: The interval [a,b] is not a bracketing interval.
You need f(a) and f(b) to have different signs (f(a) * f(b) < 0).
Consider a different bracket or try fzero(f, c) with an initial guess c.

This appears to be related to the calculation of intersection points between two functions with a bisection algorithm to determine multiple points and occurs more often at low velocities.

Stacktrace:
  [1] assert_bracket
    @ ~/.julia/packages/Roots/92AED/src/Bracketing/bracketing.jl:52 [inlined]
  [2] #init_state#59
    @ ~/.julia/packages/Roots/92AED/src/Bracketing/bisection.jl:50 [inlined]
  [3] init_state(::Roots.Bisection, F::Roots.Callable_Function{Val{1}, Val{false}, NQCDynamics.InitialConditions.QuantisedDiatomic.var"#energy_difference#20"{Float64, NQCDynamics.InitialConditions.QuantisedDiatomic.EffectivePotential{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Interpolations.ScaledInterpolation{Float64, 1, Interpolations.BSplineInterpolation{Float64, 1, OffsetArrays.OffsetVector{Float64, Vector{Float64}}, Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}}, Tuple{Base.OneTo{Int64}}}, Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}}, Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}}}, Nothing}, x₀::Float64, x₁::Float64, fx₀::Float64, fx₁::Float64)
    @ Roots ~/.julia/packages/Roots/92AED/src/Bracketing/bisection.jl:34
  [4] init_state(M::Roots.Bisection, F::Roots.Callable_Function{Val{1}, Val{false}, NQCDynamics.InitialConditions.QuantisedDiatomic.var"#energy_difference#20"{Float64, NQCDynamics.InitialConditions.QuantisedDiatomic.EffectivePotential{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Interpolations.ScaledInterpolation{Float64, 1, Interpolations.BSplineInterpolation{Float64, 1, OffsetArrays.OffsetVector{Float64, Vector{Float64}}, Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}}, Tuple{Base.OneTo{Int64}}}, Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}}, Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}}}, Nothing}, x::Tuple{Float64, Float64})
    @ Roots ~/.julia/packages/Roots/92AED/src/Bracketing/bracketing.jl:6
  [5] #init#42
    @ ~/.julia/packages/Roots/92AED/src/find_zero.jl:299 [inlined]
  [6] solve(𝑭𝑿::Roots.ZeroProblem{NQCDynamics.InitialConditions.QuantisedDiatomic.var"#energy_difference#20"{Float64, NQCDynamics.InitialConditions.QuantisedDiatomic.EffectivePotential{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Interpolations.ScaledInterpolation{Float64, 1, Interpolations.BSplineInterpolation{Float64, 1, OffsetArrays.OffsetVector{Float64, Vector{Float64}}, Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}}, Tuple{Base.OneTo{Int64}}}, Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}}, Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}}}, Tuple{Float64, Float64}}, M::Roots.Bisection, p::Nothing; verbose::Bool, kwargs::Base.Pairs{Symbol, Roots.NullTracks, Tuple{Symbol}, NamedTuple{(:tracks,), Tuple{Roots.NullTracks}}})
    @ Roots ~/.julia/packages/Roots/92AED/src/find_zero.jl:491
  [7] #find_zero#39
    @ ~/.julia/packages/Roots/92AED/src/find_zero.jl:220 [inlined]
  [8] find_zero (repeats 2 times)
    @ ~/.julia/packages/Roots/92AED/src/find_zero.jl:210 [inlined]
  [9] #find_zero#40
    @ ~/.julia/packages/Roots/92AED/src/find_zero.jl:243 [inlined]
 [10] find_zero
    @ ~/.julia/packages/Roots/92AED/src/find_zero.jl:243 [inlined]
 [11] macro expansion
    @ ~/.julia/packages/TimerOutputs/LHjFw/src/TimerOutput.jl:237 [inlined]
 [12] find_integral_bounds(total_energy::Float64, V::NQCDynamics.InitialConditions.QuantisedDiatomic.EffectivePotential{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Interpolations.ScaledInterpolation{Float64, 1, Interpolations.BSplineInterpolation{Float64, 1, OffsetArrays.OffsetVector{Float64, Vector{Float64}}, Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}}, Tuple{Base.OneTo{Int64}}}, Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}}, Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}})
    @ NQCDynamics.InitialConditions.QuantisedDiatomic ~/.julia/packages/NQCDynamics/f5aO0/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl:227
 [13] macro expansion
    @ ~/.julia/packages/TimerOutputs/LHjFw/src/TimerOutput.jl:237 [inlined]
 [14] quantise_diatomic(sim::Simulation{Classical, NQCDynamics.Calculators.FrictionCalculator{Float64, H2AgModel{Float64}}, Float64, Quantity{Int64, 𝚯, Unitful.FreeUnits{(K,), 𝚯, nothing}}, PeriodicCell{Float64}}, v::Matrix{Float64}, r::Matrix{Float64}; height::Float64, surface_normal::Vector{Float64}, atom_indices::Vector{Int64}, show_timer::Bool, reset_timer::Bool, bond_lengths::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64})
    @ NQCDynamics.InitialConditions.QuantisedDiatomic ~/.julia/packages/NQCDynamics/f5aO0/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl:272
 [15] run_MDEF_and_sample(atoms::Atoms{Float64}, cell::PeriodicCell{Float64}, model::H2AgModel{Float64}, temperature::Quantity{Int64, 𝚯, Unitful.FreeUnits{(K,), 𝚯, nothing}}, sim_time::Float64, distribution::String, number_trajectories::Int64, ensemble_type::EnsembleDistributed, electronic_friction::Bool)
    @ Main H2scatter.jl:98
 [16] md_wrapper(input_values::Tuple{Int64, Quantity{Float64, 𝐋^2 𝐌 𝐓^-2, Unitful.FreeUnits{(eV,), 𝐋^2 𝐌 𝐓^-2, nothing}}})
    @ Main multiprocess_scatter.jl:83
 [17] #4
    @ ./none:0 [inlined]
 [18] iterate
    @ ./generator.jl:47 [inlined]
 [19] collect(itr::Base.Generator{Vector{Quantity{Float64, 𝐋^2 𝐌 𝐓^-2, Unitful.FreeUnits{(eV,), 𝐋^2 𝐌 𝐓^-2, nothing}}}, var"#4#6"{Int64}})
    @ Base ./array.jl:787
 [20] (::var"#3#5")(i::Int64)
    @ Main ./none:0
 [21] iterate
    @ ./generator.jl:47 [inlined]
 [22] collect(itr::Base.Generator{UnitRange{Int64}, var"#3#5"})
    @ Base ./array.jl:787
 [23] macro expansion
    @ ./timing.jl:262 [inlined]
 [24] top-level scope
    @ multiprocess_scatter.jl:112

[Alexsp32]

Somewhat related: quantise_diatomic doesn't seem to use the correct velocities or atom masses when provided with a structure where the diatomic is not in atom indices 1 and 2

[jamesgardner1421]

@Alexsp32 Shall we close this now?