Do we really need separate `AbstractParticle` and `AbstractElement` types?

[Gregstrq]

Right now we have two types describing our particles: AbstractParticle and AbstractElement. AbstractElement contains the physical information about the element, besides its coordinates and possibly velocities.

At the same time AbstractParticle by definition is simply a container which stores AbstractElement, vector of coordinates and possibly vector of velocities. Since AbstractParticle is simply a container, it does not add anything new in terms of functionality. Moreover, since it may or may not contain velocities, we need to create at least two concrete subtypes: one without velocities and one with them, which again leads to the redundancy in the interface (we need to implement he functions that deal with AbstractElement and coordinates for both of these subtypes).

In these regards, Rachel mentioned this dicsussion of the indexing in Xtals package.

So, I propose that we remove AbstractParticle type. We can rename AbstractElement to AbstractParticle although, but it would still correspond to element data without coordinates and velocities. And getindex would then simply return a tuple (Element, coordinate) or (Element, coordinate, velocity).

[rkurchin]

I think removing this extra type is possibly a good idea, though if we do, I think getindex should just return another AbstractSystem object. This would allow the return type to be the same regardless of whether getindex was passed a single integer or a slice. It would also be consistent with the way Xtals.jl does things currently.

But I'm curious for @mfherbst's thoughts on this, as I feel like he had a reason originally for there to be a type for this...does it have to do with things that aren't spherically symmetrical needing to know about orientation or something like that perhaps?

[Gregstrq]

I am not sure about returning subsystem for single index. For example, base julia arrays act differently if you use a single index or a range.

In that regards, this is a really good question what kind of behaviour we want when we index with range: do we want to copy data and create a totally new system, or do we want to create some kind of view into existing system?

[rkurchin]

To summarize some recent discussions we've had on this (in case this issue stays open when we start publicizing the repo for broader feedback):

The (mostly-consensus) understanding is that AbstractElement types would be containers for global, static properties of particles in a simulation (for the case of a chemical element, things like atomic number, atomic mass, etc.) while AbstractParticle is "everything else." This could include quantities that might change during a simulation (e.g. position, velocity), but also things that might not change but could be specific to a simulation (e.g. ionization state perhaps).

Pragmatically, the utility in separating these is that in the context of running a simulation, AbstractElement can act as an "index" into a database of simulation parameters (pseudopotentials for DFT, interatomic potential parameters for MD, etc.).

@Gregstrq still has some concerns around how these types of parameters could get "mixed together" in inappropriate ways. I'm working on a sample implementation for an SoA version to complement the existing AoS one, which will hopefully help clarify (and/or reify) some of these issues.

[cortner]

It sounds to me rather that you want to have the properties of a particle divided into constants and variables. Isn't that overthinking it a bit at this level of abstraction? One user might want the mass to be a constant, the next a variable?

Maybe there could be an implementation of AbstractParticle that does makes such a distinction but why should that be enforced on everybody?

[rkurchin]

To my mind, the main use is the thing I said towards the end of the above comment, that you can use the information in an AbstractElement to look up parameters that you might need to actually do an iteration of your simulation without needing to store many copies of them. To be concrete for a moment, if I'm doing a DFT simulation of NaCl in a unit cell with 4 sodiums and 4 chlorines, I can use the Na and Cl AbstractElements to look up pseudopotentials, while the AbstractParticles store these AbstractElements as identifiers as well as the specific spatial coordinates of those atoms.

Conceptually, this splitting into a sort of "generic identity" (in my mind this is usually a chemical element but of course it needn't always be) and a "specific instantiation" (e.g. an atom with a set of coordinates) makes sense, but it is entirely possible I'm too stuck in the DFT reference frames and am missing something broader. Though I think this makes sense also for MD at least?

[cortner]

I understood that, but this seems already like an advanced functionality. What I didnt understand why it needs to be in AtomsBase which is supposed to be lightweight?

Also - why can't you look up the pseudopotentials if the species are stored in the AbstractParticle?

[cortner]

I guess what it comes down to is whether

• you want to enforce this on all users and even somebody who doesn't care now needs a EmptyElement parameter.; or
• you put the functionality into a prototype implementation;
• or you put an extra level of abstraction ConcreteParticle <: AbstractParticleWithElement <: AbstractParticle

and then specify which parts of the interface are supported for which level?

[rkurchin]

These are all fair points! I suppose out of those options, I probably like the adding a layer of abstraction option the best as it seems the most flexible. Keep the highest level the most generic, but still specify the interface for some more involved (and anticipated very common) use cases like what I describe above. The names of everything still TBD, of course...

I guess one question I'd put back to you is: what specific use case do you envision where you wouldn't want some way to specify the types of particles like what this provides?

[Gregstrq]

One user might want the mass to be a constant, the next a variable?

Well, the idea was to use different subtypes of AbstractElement depending on what kind of physical constants you have.

you want to enforce this on all users and even somebody who doesn't care now needs a EmptyElement parameter.

Could you describe a use case with Empy Element? It seems that any kind of particle is bound to have some static physical properties which can be packed into a structure subtyping AbstractElement.

[cortner]

Sure - any context in which there is just one particle species. Often but not always these wil be tou models

[cortner]

But more importantly it isn’t clear to me why this separation is needed /, useful

[Gregstrq]

But more importantly it isn’t clear to me why this separation is needed /, useful

Precisely because we can have single particle species or only a few species. May be the name is misleading, but AbstractElement type corresponds to the concept of species. Its subtypes are supposed to store the information that allows us to distinguish different species. Since you can have many particles of same species but with different coordinates, it is useful to have type that corresponds to species.

Sure - any context in which there is just one particle species

But each particle still has a concrete property: species it belongs to. Note that no one requires you to store a copy of species information for every particle. It is a matter of implementation, not of interface.

[cortner]

But the interface seems to enforce it?

[cortner]

And I'm not against having an interface for species at all. I'm just confused about the need to treat it separately from the concept of a Particle. Why not just have an property species ?

[Gregstrq]

Because the concept of species is different in the different contexts. For example, someone might need the mass and charge of a particle. Someone might want to treat a molecule like a point particle, but it might have additional properties related to the fact that it is molecule. For me, in NMR context, I only need spin value and gyromagnetic ratio.

So, depending on the context, one has a different set of properties which describe what species is. Hence, depending on the context, we have different types (structs) that describe what species is. And since we have to deal with a multitude of types that correspond to the same concept, it is natural to introduce an abstract type that relates to this concept.

[cortner]

That’s all fine and doesn’t contradict what I’m saying. I’m only against enforcing an extra type parameter in both an atom and a structure unless really needed. And I don’t see the need yet.

[Gregstrq]

I'm just confused about the need to treat it separately from the concept of a Particle

We might have a static and dynamic systems consisting of the same species. Although the species are the same, the physical properties of the particles are different, because in the static case there is no velocity.

[Gregstrq]

I’m only against enforcing an extra type parameter in both an atom and a structure unless really needed. And I don’t see the need yet.

What is the problem with that? The type can be automatically computed during construction, user does not have to worry about that.

[cortner]

Still doesn’t prevent you from doing what I suggested

[cortner]

First the user does have to worry because it is enforced on her by the interface. Secondly it makes the entire interface more complex than needed. It smells of over-engineering.

[cortner]

Im really not against it if there is a strong software engineering argument for it but so far I see only downsides.

[rkurchin]

I think it is useful from the interoperability perspective as having a type makes it easy to dispatch on, and this allows checks for various types of compatibility and/or interconversion between two systems expressed in different concrete implementations of the interface, but perhaps with the same AbstractElement type. (BTW, I feel like AbstractSpecies could be a good name for it, if we keep it around...just commented to this effect over in #2 )

Curious for @mfherbst's and @jgreener64's thoughts on all of this as well, though.

[cortner]

Agree on AbstractSpecies, I also like that name better and was going to propose it depending on the outcome of this discussion.

I can see in principle maybe how this might help? At least that's the first point where I'm not sure anymore. Could you give a concrete example to explain what you are envisioning? How might such a conversion work without the authors of the two types manually implementing it? Because once you are back to that, then I again don't see what is gained.

Regarding checks for compatibility ... I don't know. The default in Julia seems to me is, when it breaks, add the methods to make it work. In fact, I am not even 100% certain about the need for an AbstractParticle type. I could implement the "abstract particle" interface for types that are subtypes of something entirely different but I would also like them to behave like particles at the same time.

I used to over-engineer my type hierarchies but eventually deleted almost all abstract types. Nowadays I always start with concrete types and only add abstract type layers when I start duplicating functionality + it is needed for dispatch. The moment you enforce some abstract types on some methods, you can no longer use them for arguments that aren't derived from those and since Julia forbids you to have more than one supertype that's a real problem.

[cortner]

The interface for iterate is a nice example.

[rkurchin]

I will confess upfront that I haven't thought this through in detail, but one hypothetical I keep coming back to for interoperability is imagining trying to implement AIMD in e.g. Molly by calling out to, e.g. DFTK as the force predictor, so I'll try to spin that out in a bit more detail now...

Each simulator could (and very possibly would, since they'd be optimizing different things to maximize performance and might want to store data differently, etc.) use different concrete implementations of the interface. But let's suppose they both use the same ChemicalElement <: AbstractElement type to describe the atoms making up the system. Of course, they're using these elements to look up different specific simulation parameters (pseudos for DFT, potential parameters for MD), but individual particles are still identified by a chemical element - to be obnoxiously explicitly clear, what I mean by that specifically is that apart from things that can change during the simulation like their physical location, two particles (or I guess strictly nuclei) of the same element are identical (up to their physical location in the simulation box). This means that we know unambiguously how to "translate" the MD simulation to a DFT simulation and back again (and the structure is taken care of too since both codes use the interface) and can dispatch on this ChemicalElement type such that if there were some other DFT or MD implementation that we wanted to use instead, provided that they also implement the interface, we can easily swap it in. And if instead of the DFT, we wanted to use some ML code, we could again do it. The ChemicalElement is basically answering the question "what is this system made of?" so that we can assess compatibility.

And eventually, we can also plug this into plotting recipes that will know how to make nice visualizations with standardized colors/labels/etc. for different species again via dispatch on that parameter...

In some sense, the chemical element case is a bit simplistic, since the elemental symbol is sufficient information. I think it gets more interesting in the case where the AbstractElement is something a bit more complicated, perhaps an element including a spin configuration or ionization state, though the more specialized it is, the harder to concoct compelling interoperability cases too...

Maybe this helps?

[cortner]

It helps in that I know what you envision but not in why the complex type hierarchy is needed. Where are you dispatching?

[rkurchin]

To continue the example above, the DFT force computation function would need to take in a structure that would have the type parameter ChemicalElement, but could be any implementation of the structure interface otherwise.

[cortner]

Still not clear why you need to dispatch on that. I get that it seems this way but it really doesn’t. All you need is get_chemical_species

[jgreener64]

As others have mentioned, the difference isn't so much about static and dynamic information as information about a specific particle in the simulation (some AbstractParticle) and information on what species that particle is (some AbstractElement). This arises in many atomistic simulations, for example macromolecular MD where the specific particles might have partial charge/residue number/coordinates (depending on implementation) but there would be multiple particles of the same atom type which might have mass/chemical element information.

I agree though that not all simulations make this distinction and we might not want to impose it. The benefit of imposing it as Rachel says is that people might start to reuse the same species types, i.e.

two systems expressed in different concrete implementations of the interface, but perhaps with the same AbstractElement type

and then people could write functions that dispatch on that element type at the level of the system. The drawback is that some users might have to use empty containers.

I am not sure about the multiple levels of abstraction idea, it replaces complexity with more complexity.

[cortner]

I think we all agree that the vast majority of the simulations we will perform will want some for of species information and having the species "class" as type information is useful. This is not disputed.

I am ~disputing~ questioning only the need for (1) baking this type information into the abstract type hierarchy; and more strongly (2) the need for the abstract type hierarchy to begin with, rather than just specifying an interface.

In the absence of multiple inheritance this really limits how we can use those.

[Gregstrq]

I am disputing questioning only the need for (1) baking this type information into the abstract type hierarchy; and more strongly (2) the need for the abstract type hierarchy to begin with, rather than just specifying an interface.

The way I see it, the interface is both a set of concepts and a set of operations on those concepts. Abstract types correspond to specific concepts, representing them on the level of the code. The functions (methods) correspond to the said operations.

Also, an abstract type hierarchy facilitates the communication about these concepts between the various developers using the same interface. If something is a subtype of an abstract type corresponding to a particular concept, then everyone knows that it implements that particular concept. It also helps in debugging interoperability issues between different libraries using the same interface.

So, the real question here is whether some concept is an integral part of the interface. And it seems that everybody agrees that the concept of species is in fact an important part of the interface.

[cortner]

again - nobody disputes that it is an important part of the interface. I'm only questioning whether it must be baked into quite as "in your face" as it is now.

I appreciate that the second point about the need for an abstract type hierarchy will be more difficult for me to defend, but it is not as obvious as it might seem to you. Much of Base Julia is organized only by convention - again cf the iterator interface. When you look at how some of the most experienced Julia programmers approach package development, they typically don't bother with type hierarchies at all until it becomes necessary for organizing dispatch. Here, I haven't yet seen a concrete example that convinces me that we need it.

So in summary: I agree it is critical to have a well-defined interface for what a particle is and what a species is. But it is not at all obvious to me that a type hierarchy is needed or indeed useful.

[mfherbst]

Interesting discussion. I fully agree that it is always easier to add abstract types later rather than removing them. I have to say that (coming from the C++ world) I also only recently learned to appreciate the "design by convention" aspect.

So if I get it right, what you propose @cortner is to define the interface like define which functions exist and what sort of conceptional entities they return, but not enforce an abstract type hierarchy for these entities. While my feeling is that this might be giving up quite a lot (too much?) of the structure of the problem, it might also make it easier to reach a concensus and push out a first version of the interface (less bike-sheding about names for example :smile:). I'd generally be open for dropping the abstract types, e.g. on the element/species (btw we had species at some point in the past, right?). I'm not so sure about the system, however, because at that level in DFTK we would definitely need dispatching in the setup functions and I assume the same is true for other codes.

Could that be a compromise worth pursuing?

[cortner]

Hi Michael - thanks for considering this option.

I think there are several possible layers:

1. Keep things as they are. I would be moderately unhappy, but will live.
2. Keep AbstractParticle but remove the AbstractSpecies type information i.e. no more ChemicalElement in AbstractParticle{ChemicalElement}. This does NOT prevent us from having an AbstractSpecies type, for this option maybe any option maybe we should. Similarly, I'm unconvinced that AbstractSystem requires the Species and Particle types inside the type information. But for sure it should not be required to have both?
3. Remove all or most abstract types and work only by convention.

I am actually not sure whether 2 or 3 is better. If I'd be designing my own code then I would always go for 3 since I can add abstract types later on at any time. But this is now a community package and once there are 5-8 implementations breaking things will not be good. So I do think that we need to come to a concensus.

I'd also like to re-emphasize a previous point: I'm really not fundamentally against abstract type hierarchies, e.g. I quite like Java style interfaces. But (1) in Julia we don't get multiple inheritance; and (2) such a strict interface requires you to implement everything. Whereas in Julia we have the freedom to only implement the part of the interface that we need. It is really liberating.

I think some thought experiments might help decide? ... coming up ...

[cortner]

(1) Suppose e.g. I do some ML on atoms. Then an atom object ::MyAtom might act at the same time as an Atom and as a node in a graph. Or in my case an ACE.State which contains features (such as position, species, ...). Now should I have MyAtom <: AbstractAtom or MyAtom <: ACE.AbstractState or MyAtom <: GraphNode? I can't have all three.

(2) I keep asking myself why I need to put the species information into the type parameters?

function get_pseudopotential(at::AbstractAtom) 
    sym = chemical symbol(at) 
    return get_this_from some_other_library(sym)
end

(3) Worst-case scenario, if you really need to dispatch on the type of the species, you can do

do_something(at::AbstractAtom) = do_something(at, get _species(at))
do_something(at::AbstractAtom, s::AbstractSpecies) = ...generic
do_something(at::AbstractAtom, s::ChemicalSpecies) = ...special

How many cases can we envision where we need to do this? If many then that would be an indicator that we do need to consider putting the type information in (i.e. option 1). But I would be very surprised. I think it is much more likely that we will just write generic code and then if the concrete species doesn't implement some functionality it will break.

(4) @mfherbst can you explain why you need this in the DFTK setup stage?

[cortner]

Maybe an example from a different domain helps: I want to multiply A * B where A, B contain some objects TA = eltype(A), TB = eltype(B) that are not numbers. If mul would require AbstractMatrix{<: Number} then I'd have to reimplement this from scratch. But it doesn't - so now I just stick in my types and run A*B and it will break, e.g. it will tell me that *(::TA, TB) is not defined. Ok, so I implement it and try again. Next Julia will tell me that I need identity(::TA). Ok, so I implement it and try again ... The point is that matrix multiplication is defined in terms of the operations required to carry it out. And any type that supports those can now be used to perform it.

[Gregstrq]

@cortner, We have decided to go with the third option and drop all the abstract types except AbstractSystem from the interface. However, we will keep the type hierarchy that we came up with in the concrete implementations we are working on right now.

[cortner]

thanks for the update. I'm a bit unsure how this is different from 1? Anyhow ... I will see where it goes and then comment again. Time for me to shut up.

[mfherbst]

To comment on your question Christoph:

(4) @mfherbst can you explain why you need this in the DFTK setup stage?

In DFTK we have a certain flexibility to setup DFT calculations. The most high-level is model_LDA(lattice, atoms) for example and the most low-level Model(lattice; kwargs...) where basically you can bootstrap any operator you want including analytic potentials or potentials generated by the atoms. The point now is that the AbstractSystem structure/convention in AtomsBase sort of supplies both lattice and atoms and potentially a few other things as well. So at the DFTK level this info somewhere needs to be split apart in a transparent way, such that the full flexibility we have is still kept. At that level we would need multiple dispatch on the "system" object and therefore at least an AbstractSystem abstract type would be useful for us.

[cortner]

I don't fully see the issue, but also I'm a lot less worried about the AbstractSystem than the AbstractAtom.

Just to explain why I don't see it: why is it not enough to just define the interface for get_lattice etc.? That's where the dispatch should happen, no?

[mfherbst]

Essentially for compatibility with the interface as we support it now, where the lattice can be just a matrix of numbers or some python object (for ASE compatibility). Especially for the matrix of numbers case it's a bit hackish to define a get_lattice method at all (get_lattice(a::AbstractMatrix) = a is a bit abusive in my opinion).

[rkurchin]

To be clear though, I don't think there's anything stopping you from defining the interface functions on your own type but not actually subtyping AbstractSystem. Some things might not be quite so smooth, but I think it should still work.

And if it wouldn't work, then IMO that would be a problem and it's worth discussing. For example, could what we're after be achieved by some sort of trait instead of actually inserting something into the type hierarchy? I could imagine that other package developers might not always want to make their type a subtype of AbstractSystem, for instance, if they've already chosen to subtype AbstractArray...

[cortner]

I could imagine that other package developers might not always want to make their type a subtype of AbstractSystem, for instance, if they've already chosen to subtype AbstractArray...

yes, that is precisely the point. It is not the end of the world, we just need to start converting then. But if we worked purely with an interface, then no need!

[rkurchin]

I'm going to try gradually removing pieces and only ever commit working versions over at #23, so you can keep an eye out there. So far, I've pushed a version that removes AbstractElement but keeps everything else. I'll plan on going all the way to introducing no abstract types so that we can see what that case might look like as well, which I think will help us decide if that's what we want or if some intermediate case is best.

Thanks @cortner for forcing us to think through all this really carefully; this will undoubtedly be better for it, wherever we land!

[jgreener64]

I can see both sides but personally I would lean towards keeping AbstractSystem for its general use in dispatching, one example being defining plot recipes.

Ironically the use of having something be a subtype of AbstractArray points to the future potential of AbstractSystem. From the Julia interface docs:

If a type is defined as a subtype of AbstractArray, it inherits a very large set of rich behaviors including iteration and multidimensional indexing built on top of single-element access.

I worry we lock ourselves out of this potential if we don't introduce AbstractSystem.

[cortner]

I agree that’s another important point - to inherit behaviour. As I said I don’t have a strong view on AbstractSystem.

That said - one can always “inherit” behaviour via a macro

[rkurchin]

@cortner, curious for your thoughts on the current state of #23 as regards the discussion above (that should get merged later today most likely; I need to add a few more tests)

[cortner]

see my comments in the review. From here on I think I would need to test this.