JuliaMolSim / Molly.jl

Molecular simulation in Julia
Other
371 stars 51 forks source link

UndefVarError: `σ6` not defined ERROR when simulating the NPT ensemble #169

Closed GDufenshuoo closed 3 months ago

GDufenshuoo commented 3 months ago

When I attempt to simulate the NPT ensemble, I consistently encounter this error. Could there possibly be a missing definition here?

ERROR:

ERROR: UndefVarError: `σ6` not defined
Stacktrace:
 [1] potential_energy
   @ ~/.julia/packages/Molly/B4xv7/src/interactions/fene_bond.jl:59 [inlined]
 [2] potential_energy_pair_spec!(pe_vec::Vector{…}, coords::Vector{…}, atoms::Vector{…}, pairwise_inters_nonl::Tuple{}, pairwise_inters_nl::Tuple{…}, sils_1_atoms::Tuple{}, sils_2_atoms::Tuple{…}, sils_3_atoms::Tuple{…}, sils_4_atoms::Tuple{}, boundary::CubicBoundary{…}, energy_units::Unitful.FreeUnits{…}, neighbors::NeighborList{…}, n_threads::Int64, ::Val{…})
   @ Molly ~/.julia/packages/Molly/B4xv7/src/energy.jl:201
 [3] potential_energy_pair_spec
   @ ~/.julia/packages/Molly/B4xv7/src/energy.jl:100 [inlined]
 [4] potential_energy(sys::System{…}, neighbors::NeighborList{…}; n_threads::Int64)
   @ Molly ~/.julia/packages/Molly/B4xv7/src/energy.jl:84
 [5] potential_energy
   @ ~/.julia/packages/Molly/B4xv7/src/energy.jl:74 [inlined]
 [6] apply_coupling!(sys::System{…}, barostat::MonteCarloBarostat{…}, sim::Langevin{…}, neighbors::NeighborList{…}, step_n::Int64; n_threads::Int64)
   @ Molly ~/.julia/packages/Molly/B4xv7/src/coupling.jl:199
 [7] simulate!(sys::System{…}, sim::Langevin{…}, n_steps::Int64; n_threads::Int64, run_loggers::Bool, rng::Random._GLOBAL_RNG)
   @ Molly ~/.julia/packages/Molly/B4xv7/src/simulators.jl:352
 [8] simulate!(sys::System{…}, sim::Langevin{…}, n_steps::Int64)
   @ Molly ~/.julia/packages/Molly/B4xv7/src/simulators.jl:323
 [9] top-level scope
   @ ~/Documents/Github/BetaRelaxiaton/Glass/20240217.jl:113
Some type information was truncated. Use `show(err)` to see complete types.

And it can be fixed like this

@inline function potential_energy(b::FENEBond, coord_i, coord_j, boundary)
    dr = vector(coord_i, coord_j, boundary)
    r = norm(dr)
    r2 = r^2
    r2inv = inv(r2)
    r6inv = r2inv^3
    σ6 = b.σ^6    # add
    r02 = b.r0^2
    uwca = zero(b.ϵ)
    if r < (b.σ * 2 ^ (1 / 6))
        uwca = 4 * b.ϵ * ((σ6 * r6inv) ^ 2 - σ6 * r6inv) + b.ϵ
    end
    return -(b.k / 2) * r02 * log(1 - r2 / r02) + uwca
end
jgreener64 commented 3 months ago

Thanks for reporting, fixed in d651459.