WIP: Stillinger Weber potential as a first prototype

[cortner]

This is the first of two PRs to transfer existing empirical potential implementations from JuLIP. This one is the easier one but still interesting for our special issue plans.

At the moment this doesn't work, I'm just putting it here to start discussing.

[tjjarvinen]

AtomsCalculators interface is now working and I added (very crude) pairpotential as well. The performance is not optimal yet, but will do for now.

I added support to use other than atomic number when defining potential. To use this feature you need to redefine get_id function the default is

@inline get_id(at, V, i) = atomic_number(at, i)

To get this working, I also redefined get_neighbours function to

function get_neighbours(at, V, nlist, i)
   Js, Rs = NeighbourLists.neigs(nlist, i)
   Zs = [ get_id(at, V, j) for j in Js ]
   z0 = get_id(at, V, i) 
   return Js, Rs, Zs, z0 
end

[cortner]

is this ready to review and merge?

[tjjarvinen]

I have not looked for Stillinger Weber part. Site and pair potentials are ready, but need tests.

[cortner]

but need tests

I noticed, please add those for pair at least. For site, the SW is already such a test I think.

[cortner]

ok - finished SW site energy tests. Maybe you can think about how to add total energy and total force tests along similar lines. I should try to get the hessian or at least the preconditioner ready for the next meeting.

[tjjarvinen]

I have added tests for PairPotential and Benchmarks.

I did also additional benchmarking. The main point is that parlist generation is slow - almost two orders of magnitude penalty. Also neigs call is slow and improving it would lead to further performance improvements.

Finally I have been looking about the allocations and while we could use ObjectPools, I would actually prefer that we switch to Bumper. There are couple of reasons for this:

1. Bumper uses core Julia intended method task_local_storage
2. True thread (task) safety
3. There seems to be general movement in to use it
4. Support for StaticCompiler.jl

For these the StaticCompiler support is a point that is very importan for the future use, as it would allow generating statically compiled C-libraries.

So, basically I would like to replace ObjectPools with Bumper.

[cortner]

@tjjarvinen I'm very open-minded about switching to Bumper. ObjectPools was something I wrote because I didn't have anything that could conveniently manage reusing arrays in a way that was source-code localized. But I have no desire to maintain my own garbage collector (which is basically what ObjectPools is).

I would however like to understand how it performs first. Two options:

• we merge what we have as is and remove ObjectPools in a new PR
• we create a new branch now where we remove ObjectPools, confirm the performance, and then merge it into the current PR.

Let me know what you prefer.

There are two ways that I use ObjectPools that are related but distinct and both have to be covered by Bumper to make the switch:

1. Temporary arrays: in some function calculate acquire a temporary array A from the object pool, do some calculators with it, then release!(A) back into the object pool.
2. Reuse outputs: Some function calculate returns an array A, the caller can use A and when it is no longer needed it can be released back into the object pool via release!(A).

My understanding of a first reading of the documentation is that Bumper.jl can do 1 but not 2. Maybe we can come up with an alternative here but my feeling is that whatever we do it will delocalize the memory management again.

[tjjarvinen]

We should first finnish this PR and then look for bumper in separate case. There will be some testing needed to see what we can do with.

Like you said the case 1. is no problem, the case 2. can be also be done at least in some cases, but I need to play around a bit to see what we can do with it in the end.

[cortner]

good, then we are on the same page. Let's both try it out and then discuss.

[tjjarvinen]

I added parameter estimation to pair potentials. It should be now in a shape that we could use it to build and train ACE models.

The AtomsCalculators parameter estimate is not ready. I am using this implementation as a test bed before submitting PR for AtomsCalculators.

In short you can do now

using EmpiricalPotentials
using AtomsIO
using AtomsCalculators

fname = joinpath(pkgdir(EmpiricalPotentials), "data", "TiAl-1024.xyz")
data = load_system(fname)

# energy min, min location, atomnumber1, atomnumber2, cutoff
lj = LennardJones(-0.5u"meV", 2.8u"Å",  13, 22, 6.0u"Å"; parametric=true)

AtomsCalculators.potential_energy(data, lj, lj.parameters)
AtomsCalculators.forces(data, lj, lj.parameters)
AtomsCalculators.virial(data, lj, lj.parameters)

[rkurchin]

This looks generally sensible to me but I would think we might want @jgreener64's thoughts before moving forward, since we would obviously want this to play nice with Molly.

[jgreener64]

Looks good to me for now, I'll need more time to think about how this integrates with Molly but that shouldn't hold this up.

[cortner]

merge now and add hessians in a separate PR?

[tjjarvinen]

I would say that merge now. Helps to keep PRs smaller.

[cortner]

I will wait until noon PST and if I don't hear any dissent until then I'll merge the PR.