Adding a Glue Type / N-Body Empirical Potential - the Finnis-Sinclair 1984 potential for single element systems

[eschmidt42]

Glue: Finnis-Sinclair Potential

Hi, I'm new to Julia and was very happy to find that this Molly.jl repo exists. Playing with it (great documentation btw! 🙂) and wanting to simulate metals, I noticed potentials for metals seemed to be missing. So, I've written a small piece of code implementing one of the standard empirical potential types for metal, the Finnis-Sinclair empirical potential, for myself. btw thank you for documenting, how to best extend Molly.jl, so even newbies in Julia like myself can get started making custom changes 🙂. Metal potentials are not on your list of features you want to add as far as I can see, but I thought this still might be interesting for some users, hence this pull request.

Main Change

• added the Finnis-Sinclair form (Finnis and Sinclair 1984, A simple empirical N-body potential for transition metals glue type empirical potential as a new interaction: src/interaction/glue_fs.jl

Note: The "glue" (kind of like local electron density, but not really) is a scalar pre-computed for each atom and is a required quantity to compute energies and forces. So simply iterating over atom pairs to compute energy and force contributions is insufficient, without having pre-computed/updaed the glue values. Hence, in a few places, the MD code needed to be adjusted to allow for this extra step.

Minor Changes

• added a notebook fs.ipynb to document the usage of the added interaction type
• modified accelerations in src/forces.jl: to enable glue calculation before pair contributions to the forces
• modified potential_energy in src/analysis.jl to enable the use of a potential_energy function for the glue interaction type
• modified src/loggers.jl by adding new loggers which proved useful for debugging: GlueDensityLogger, VelocityLogger, ForceLogger
• modified Simulation in src/types.jl by adding glue_densities as an array where per-atom glue values are stored for the glue potential
• added test/glue.jl as a new test to run to verify the correctness of the predicitons of potential energies, forces using reference values from the Finnis and Sinclair paper linked above as well as testing the behavior of the glue scalar values.

Note: A peculiarity of the Finnis-Sinclair potential is that the potential energy of an atom depends on the scalar glue value of that atom by taking the square root. Hence, cases in which the glue value, for whatever reason, becomes negative, will crash a simulation, since it is not defined in the model how to handle complex potential energy values. But negative glue values should really never occur, and if they do they indicate that something else is wrong with the simulation (crystal configured incorrectly or so). A test is present in test/glue.jl to verify that this doesn't happen due to algebraic sign errors or so in the implementation of the potential itself for a correctly configured simulation.

Next steps

Since this is a minimal first implementation of one empirical potential type for metals, which is only parameterized for single element systems (Fe, Nb, ...) there are a couple more things that could be done. The more relevant next steps could be to:

• add shifted force / potential version of the glue potential
• add parallelism for the glue computation
• test CUDA compatibility
• add other types of glue potentials, like the Embedded Atom method and
• verify the correctness of the implementation for two, three and more element systems (the implemented Finnis-Sinclair paper only deals with single element systems).

[eschmidt42]

Glad to hear this could be a relevant contribution and thank you for your suggestions 😃. I'll incorporate them them over the next couple of days.

[jgreener64]

Thanks for making this PR, and I'm glad you liked the package and docs. This potential would certainly be welcome in some form. As a first pass, addressing @SebastianM-C's useful comments would be good. In particular we should avoid adding any more dependencies unless they are very light or add something crucial.

More broadly though, I've been thinking about this PR over the last few days and how it backs up something we have discussed previously: Molly needs a "free-form" interaction type to go with SpecificInteraction and GeneralInteraction (which is really a pairwise interaction), where the new type has all atoms available to it. Then you could do the glue calculation with this interaction type and update some array in the interaction struct, or even update the properties of a custom atom type.

That way, you wouldn't need to change the force/energy calculation functions or modify the Simulation type, both of which I am uneasy with. In particular, code like if typeof(inter) <: GlueInteraction would be dealt with via dispatch rather than branching.

One of the principles of the package is that you should be able to implement custom potentials solely by adding your own interaction types. I'd rather add a new class of interactions which has all atoms available to facilitate that, rather than special casing the summation functions.

I might have a chance to implement such an interaction type over Easter, and if so will look to make it work with this example, at which point this PR could be reworked to fit on top of that. All ideas welcome of course.