JuliaMolSim / AtomsBase.jl

A Julian abstract interface for atomic structures.
https://juliamolsim.github.io/AtomsBase.jl/
MIT License
81 stars 16 forks source link

AtomView and system getindex #74

Closed ejmeitz closed 1 year ago

ejmeitz commented 1 year ago

How am I intended to implement getindex(v.system, v.index, x) so that getindex(v::AtomView, x::Symbol) works without calling get()? If I use get() there is a cyclic reference and I don't think you can use getfield() with an AtomView.

Couldn't find a library that uses the AtomView approach to AtomsBase besides Molly. Maybe this is an issue with Molly's implementation, but we do not have arrays for mass, charge, atomic number etc stored as part of the system object so those are not accessible directly via getfield(v.system,x)[i] or something along those lines. I can access specific properties by name (e.g. atomic_mass(v.system, 1)) but not for a generic property passed as a symbol to getindex().

Base.getindex(v::AtomView, x::Symbol) = getindex(v.system, v.index, x)
function Base.get(v::AtomView, x::Symbol, default)
    hasatomkey(v.system, x) ? v[x] : default
end
mfherbst commented 1 year ago

Thanks for your question and sorry for the late reply.

I have also spotted the cyclic reference once and that's probably something that should be resolved. I'm not sure how, though as each implementation is sort of convenient for certain scenarios, e.g. https://github.com/JuliaMolSim/AtomsBase.jl/pull/68 added parts of the cyclic thing to AtomsView.

Regarding examples: AtomView is used in ExtXYZ and the FastSystem implementation, e.g. https://github.com/libAtoms/ExtXYZ.jl/blob/master/src/atoms.jl#L249.

Regarding your question, I think the simplest is you define the position, velocity etc. functions in the way they return the right thing in their one-argument position(::MollySystemType) and/or two-argument position(::MollySystemType, i::Integer) form and then

function Base.getindex(system::MollySystemType, i, prop::Symbol)
  if hasatomkey(system, x)
    getproperty(AtomsBase, prop)(system, i)  # uses the two-argument version
    # or
    getproperty(AtomsBase, prop)(system)[i]  # uses the one-argument version
  else
    # throw keyerror
  end
end

I have to admit that is not particularly nice, but should work. Long-term we should think weather we want to have something like this as the default solution in AtomsBase, such that downstream users do not need to worry so much about this. Feel free to make a PR here if you want, but I would probably first try it in Molly and see it does not cause unexpected complications.

ejmeitz commented 1 year ago

Unfortunately, this does not work for the :charge atomkey as it is not part of the AtomsBase interface. If I remove charge from my AtomsBase.atomkeys() function then all the tests pass but if I add anything that does not have an AtomsBase.<property>function it will break. Basically since there is no AtomsBase.charge(sys, i) function which breaks the function you provided. It works fine for position, velocity, atomic_mass, and atomic_symbol which all have AtomsBase overloads.

Maybe I do not understand the AtomViewAPI but charge is not one of the fields for an AtomView even though its in our system. So should I just remove charge from my atom keys (leaving only the default fields)? Is there some way to get this info into the AtomView or any non-mandatory key for that matter?

mfherbst commented 1 year ago

Actually the error is on my end here. I completely forgot that there are optional atomskeys when I was suggesting my implementation. Charge is of course one of them. In that case I'd suggest to not use atomkeys, but instead explicitly list the mandatory keys (to use the getproperty(AtomsBase, prop)(system, i) trick) and then add further explicit branches to do what there is to be done for the charges etc.

ejmeitz commented 1 year ago

No worries, thanks for all your help!

ejmeitz commented 1 year ago

For documentation, this is what the final solution was:

function Base.getindex(system::Union{System, ReplicaSystem}, i, x::Symbol)
    atomsbase_keys = (:position, :velocity, :atomic_mass, :atomic_number)
    custom_keys = (:charge, )
    if hasatomkey(system, x)
        if x ∈ atomsbase_keys
            return getproperty(AtomsBase, x)(system, i)
        elseif x ∈ custom_keys
            return getproperty(Molly, x)(system, i)
        end
    else
      throw(KeyError("Key $(x) not present in system."))
    end
end