Redesign methods and structures for multi-element and IAP combinations

[dallasfoster]

Some of the more realistic use cases in the CESMIX universe are systems with multiple atom types (i.e., HfO2, HfB, Hydrogen/Palladium cluster). In such systems there are often 2 modeling choices that occur: (1) the forces between atoms is modeled as a combination of potentials (i.e., an empirical potential like ZBL and a basis potential like SNAP) and (2) the function that is used to calculate the force between two atoms (atom_i, atom_j) might depend on the elements of atom_i and atom_j (i.e., in the GaN example Ga-Ga, Ga-N, and N-N are all modeled differently).

In the hope of providing flexibility and ease of use for users of InteratomicPotentials.jl, I think there should be some mechanism available to easily describe the modeling choices of the system. An example could be something like what is implemented in Turing.jl:

## GaN <: CompiledModel
@model GaN(q_Ga, q_N, rcut_Ga_Ga, rcut_N_N, rcut_Ga_N, ε_Ga, σ_Ga, ε_N, σ_N, a, ρ) = begin
 # Define Galium-Galium Interaction
   :Ga ~ :Ga -> Coulomb(q_Ga, q_Ga, e0, rcut_Ga_Ga) + LennardJones(ε_Ga, σ_Ga)

 # Define Nickel-Nickel Interaction
   :N ~ :N -> Coulomb(q_N, q_N, e0, rcut_N_N) + LennardJones(ε_N, σ_N)


 # Define Galium-Nickel Interaction
  :Ga ~ :N -> Coulomb(q_Ga, q_N, e0, rcut_Ga_N) + BornMayer(a, ρ)
end 

With such a model, you could compute the energy and force like so:

potential_energy(A, GaN) # A <: AbstractSystem
force(A, GaN) 

There are some challenges with this approach. First, I have no idea how to implement such a thing in Julia (help?). Second, if a CompiledModel contained a BasisPotential that needed to be fitted (i.e., SNAP). It would be helpful to extract the descriptors associated with that potential. Normally you would invoke the code

descriptors = evaluate_basis(A, snap)

But this could become more complicated take the following example:

## Suppose we are modeling a cluster system with Hydrogen and Palladium. Suppose we want Hydrogen-Hydrogen,  
## Hydrogen-Palladium, and Palladium-Palladium atoms, to interact with the same SNAP potential but we also want
## Hydrogen-Hydrogen atoms to have an extra Bessel potential.
@model Cluster(β_SNAP, β_Bessel, SNAPParams, BesselParams) = begin
 # Define Hydrogen and Palladium Interaction
   :P ~ :P -> SNAP(β_SNAP, SNAPParams)
   :P ~ :H -> SNAP(β_SNAP, SNAPParams)
  :H ~ :H -> SNAP(β_SNAP, SNAPParams) + Bessel(β_Bessel, BesselParams)
end 

This model only needs the single SNAP potential (not three separate SNAP potentials) and single Bessel potential. If I call the descriptors for this model, I want to make sure that they are separated into the SNAP descriptors and Bessel descriptors. Furthermore, for solving linear systems (those generated by SNAP, ACE, etc.), it is important to separate out the so-called 'reference' energies and forces provided by the empirical potentials from the contributions of the BasisPotentials.

[jrdegreeff]

I don't fully understand the BasisPotentials, but my first thought at least for EmpericalPotentials, would be something like this:

struct HybridPotential <: ArbitraryPotential
    interactions::Dict{Set{Symbol}, EmpericalPotential}
    rcutoff::AbstractFloat
end

get_potential(p::HybridPotential, A::AbstractSystem, i::Integer, j::Integer)
    p.interactions[Set([atomic_species(A, ii), atomic_species(A, jj)])]
end

function force(A::AbstractSystem, p::HybridPotential)
    f = fill(zeros(3), length(A))
    nnlist = neighborlist(A, p.rcutoff)

    for ii in 1:length(A)
        j = nnlist.j[ii]
        r = nnlist.r[ii]
        R = nnlist.R[ii]
        for (jj, rj, Rj) in zip(j, r, R)
            fo = force(Rj, rj, get_potential(p, A, ii, jj))
            f[ii] += fo
            f[jj] -= fo
        end
    end
    SVector{3}.(f)
end

function potential_energy(A::AbstractSystem, p::HybridPotential)
    e = 0.0
    nnlist = neighborlist(A, p.rcutoff)

    for ii in 1:length(A)
         j = nnlist.j[ii]
        R = nnlist.R[ii]
        for (jj, Rj) in zip(j, R)
            e += potential_energy(Rj, get_potential(p, A, ii, jj))
        end
    end
    e
end

function virial(A::AbstractSystem, p::HybridPotential)
    v = virial_stress(A, p)
    v[1] + v[2] + v[3]
end

function virial_stress(A::AbstractSystem, p::HybridPotential)
    v = zeros(6)
    nnlist = neighborlist(A, p.rcutoff)

    for ii in 1:length(A)
        r = nnlist.r[ii]
        R = nnlist.R[ii]
        for (rj, Rj) in zip(r, R)
            fo = force(Rj, rj, get_potential(p, A, ii, jj))
            vi = rj * fo'
            v += [vi[1, 1], vi[2, 2], vi[3, 3], vi[3, 2], vi[3, 1], vi[2, 1]]
        end
    end
    SVector{6}(v)
end

This is just off the top of my head, so the names could probably use some work. I could also come up with a way to dispatch things to avoid having the duplicated code (basically something along the lines of make this the default implementation for all EmpericalPotentials, make the default implementation of get_potential just fall through for normal EmpericalPotentials, and have HybridPotential subtype EmpericalPotential), but this seemed like a good format to get the idea across.

So then the usage would be something like this:

GaN = HybridPotential(
    Dict(
        Set([:Ga, :Ga]) => _____,
        Set([:N, :N]) => _____,
        Set([:Ga, :N]) => _____
    ),
    rcut
)

If you don't want all the initialization syntax, we could probably add some sort of macro, though I'll admit that isn't something I have experience with. Does this design accomplish the main things you are looking for? The one thing I can see it doesn't do is handle different rcut for each potential. If that is important we would have to change the order of steps (i.e. have separate ball trees for each pair of atom types (3 in the case of this example) -- this feels doable but would certainly add some more overhead).

You can also pretty easily add a method of Base.+ so that the sum of two potentials yields a potential that dispatches the force, potential_energy, viral_stress, and virial functions in roughly 10 lines. Something like this:

struct CombinedPotential{A<:ArbitraryPotential, B<:ArbitraryPotential} <: ArbitraryPotential
    p1:A
    p2:B
end

Base.:+(p1::A, p2::B) where {A <: ArbitraryPotential, B <: ArbitraryPotential} = CombinedPotential{A, B}(p1, p2)

force(R::AbstractFloat, r::SVector{3,<:AbstractFloat}, p::CombinedPotential) = force(R, r, p.p1) + force(R, r, p.p2)
potential_energy(R::AbstractFloat, p::CombinedPotential) = potential_energy(R, p.p1) + potential_energy(R, p.p2)

That's the gist of it anyway. I'd have to look a bit more how you are doing things for types other than the EmpericalPotential side of the tree.

If you want me to implement any of this stuff, just let me know and I'll gladly code it up in a PR.

Reading your summary of the SNAP side of things, I have some thoughts about how you could go about it, but there is enough of a gap in my understanding that I'll hold off for now. If you want my (half-formed) thoughts on that, it will be easier to convey live when I can get some real time answers to fill in those gaps.

[emmanuellujan]

Regarding the modeling of problems such as GaN, an alternative based on multiple dispatch could be the following:

lj_Ga_Ga = LennardJones(ε_Ga, σ_Ga)​; ​lj_N_N = LennardJones(ε_N, σ_N)
c_Ga_Ga = Coulomb(q_Ga, q_Ga, e0)​; ​c_N_N = Coulomb(q_Ga, q_N, e0)​;
c_Ga_N = Coulomb(q_Ga, q_N, e0)​; ​bm_Ga_N = BornMayer(a, p)


function potential_energy(particle1::Ga, particle2::Ga)
  r = position(particle1) - position(particle2)
  return  potential_energy(r,c_Ga_Ga) + potential_energy(r,lj_Ga)
end

function potential_energy(particle1::N, particle2::N)
    r = position(particle1) - position(particle2)
    return  potential_energy(r,c_N_N) + potential_energy(r,lj_N)
end

function potential_energy(particle1::Ga, particle2::N)
    r = position(particle1) - position(particle2)
    return  potential_energy(r,c_Ga_N) + potential_energy(r,bm_Ga_N)
end

function potential_energy(particle1::N, particle2::Ga)
   return potential_energy(particle2, particle1)
end

This would avoid directly modeling the GaN potential.

[emmanuellujan]

If we need more flexibility, we can add one more parameter, something like this:

struct GaN <: InteratomicPotential end

function potential_energy(particle1::Ga, particle2::Ga, p::GaN)
  r = position(particle1) - position(particle2)
  return  potential_energy(r,c_Ga_Ga) + potential_energy(r,lj_Ga)
end

function potential_energy(particle1::Ga, particle2::Ga, p)
  r = position(particle1) - position(particle2)
  return  potential_energy(r, p)
end

If we need to compute SNAP plus ZBL:

s = SNAP(...); z = ZBL(...)
p = potential_energy(particle1, particle2, s) + potential_energy(particle1, particle2, z)

or

p = sum( potential_energy.(particle1, particle2, [s, z]) )

That's what comes to my mind at the moment :)

[mfherbst]

From having read this, I have to say I like @emmanuellujan's approach best.

But I think before designing the interface we should think about a few crazy things one would want to be able to do and try to ensure the interface is able to do it, e.g.

• Switch potentials based on a distance cutoff or something (angle ? dihedrals? local environment?)
• Easily switch or tweak parameters (e.g. update them or differentiate with respect to them)
• What about linear combinations, not just addition?

Being able to have full flexibility in defining potentials is quite important in my opinion for our toolkit and as these examples illustrate combinatorial complexity quickly explodes, so I would not start writing such high-level macros or convenience functions until setting up the problem becomes a pain. Undoing these abstractions later is much harder than adding them.

I often found (similar to what @emmanuellujan suggested) plain functions to be a good starting point. Why not just provide a good set of functions that are easy-to-use building blocks and leave the composition completely to the user (and ask him to define his types for that purpose)? The part I don't like about @emmanuellujan's approach that I am not so sure about is to use types to encode atoms. For me this is just "data" and not so much worth a type. Also I guess that leads to potentially a lot a lot of types (think about 4-body interactions for example).

So perhaps just provide a Functor for each kind of interaction as basic building blocks (holding the state of the parameters) and then leave the combination of it all completely to the user. Then for standard tasks (i.e. one type of potential for each atom type and interaction) provide a small set of convenience assemblers to do that combination. As time passes we will notice what are frequent patterns we need and we can add more convenience assemblers for these. We also keep full flexibility: If people need more fancy things, they can use those building blocks to do their own stuff. In the end as long as the final evaluation function they produces matches the requirements they are good to go.