Adding a Glue Type / N-Body Empirical Potential - the Finnis-Sinclair 1984 potential for single element systems #32
Conversation
…on of glue forces using the general_inters loop
…ivates with FordwardDiff.derivate to re-use the non-derivative forms also used to successfully test ground state energies
…s-Sinclair 1984 empirical model
…o the core components
Merged small changes removing commented out lines
There was a problem hiding this comment.
Hello! Thank you for contributing.
I personally was not aware of this kind of potentials, but it seems like an interesting thing to add.
I worked last summer (during Julia Summer of Code) on optimizing Molly for performance and I have a few suggestions that could help your implementation to be faster.
Regarding the documentation, I would suggest to add it to the docs folder in the format of Documenter pages since that way it would be more integrated with the rest of the package documentation.
Project.toml
Outdated
| KernelDensity = "5ab0869b-81aa-558d-bb23-cbf5423bbe9b" | ||
| LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" | ||
| NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce" | ||
| Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" |
There was a problem hiding this comment.
Plots should not be added as a direct dependency. Instead you can add it in the docs/Project.toml if you need to use the package for documentation purposes. In that case, you should also have the documentation prepared for Documenter.jl accordingly.
src/interactions/glue_fs.jl
Outdated
| end | ||
|
|
||
|
|
||
| kb = 8.617333262145e-5 # eV/K |
There was a problem hiding this comment.
For performance reasons you should avoid global variables. You could make kb a parameter in the interaction as is done in the case of electrostatic interaction for example. See https://docs.julialang.org/en/v1/manual/performance-tips/#Avoid-global-variables
Project.toml
Outdated
| @@ -8,10 +8,14 @@ BioStructures = "de9282ab-8554-53be-b2d6-f6c222edabfc" | |||
| CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" | |||
| Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" | |||
| Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" | |||
| Crystal = "3c6eccdf-2a89-4c24-a1d4-ff210daa8476" | |||
There was a problem hiding this comment.
If you need some packages for tests, you should add them as test dependencies in the [extras] section, not as direct package dependencies.
src/interactions/glue_fs.jl
Outdated
| element_pairings = [string(el, el) for el in elements] | ||
| element_pair_map = Dict(pair => i for (i,pair) in enumerate(element_pairings)) | ||
|
|
||
| df = DataFrame( |
There was a problem hiding this comment.
I'm not sure if adding the DataFrames.jl dependency if worth it just for this convenience functionality.
src/interactions/glue_fs.jl
Outdated
|
|
||
| Derivative of the glue density function. | ||
| """ | ||
| ∂glue_∂r(r, β, d) = ForwardDiff.derivative(r -> glue(r,β,d), r) |
There was a problem hiding this comment.
In this case I think that the derivative could be done analytically since the formula is not very complicated.
Project.toml
Outdated
| Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" | ||
| ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" |
There was a problem hiding this comment.
I don't think there is a need to use ForwardDiff for the derivatives introduced by this change. The derivatives expressions could be derived analytically and used directly. This might also be slightly more performant than doing them numerically.
src/loggers.jl
Outdated
|
|
||
| function ForceLogger(T, n_steps::Integer; dims::Integer=3) | ||
| return ForceLogger(n_steps, | ||
| Array{SArray{Tuple{dims}, T, 1, dims}, 1}[]) |
There was a problem hiding this comment.
You could just use Vector{SVector{dims}} here for readability.
| Array{SArray{Tuple{dims}, T, 1, dims}, 1}[]) | |
| Vector{SVector{dims}}[]) |
src/interactions/glue_fs.jl
Outdated
| element_pair_map::Dict | ||
| params::DataFrame |
There was a problem hiding this comment.
These are not a concrete types and will cause performance problems. I would recommend using a (parametrized) NamedTuple or another (concretely typed) struct instead. See https://docs.julialang.org/en/v1/manual/performance-tips/#Avoid-fields-with-abstract-type
src/interactions/glue_fs.jl
Outdated
|
|
||
| Energy derivative given glue density. | ||
| """ | ||
| ∂Uglue_∂ρ(ρ,A) = ForwardDiff.derivative(ρ -> Uglue(ρ,A), ρ) |
There was a problem hiding this comment.
This could be simply derived analytically
| ∂Uglue_∂ρ(ρ,A) = ForwardDiff.derivative(ρ -> Uglue(ρ,A), ρ) | |
| ∂Uglue_∂ρ(ρ,A) = -A /(2 * √ρ) |
src/forces.jl
Outdated
| end | ||
|
|
||
There was a problem hiding this comment.
Remove whitespace
|
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. |
|
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 That way, you wouldn't need to change the force/energy calculation functions or modify the 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. |
…b as an attribute of FinnisSinclair
…d testsets for legitibility of the test summary in case of failures
… and changes in the package constructing the bcc crystal
test/runtests.jl
Outdated
| @@ -22,6 +22,8 @@ end | |||
|
|
|||
| CUDA.allowscalar(false) # Check that we never do scalar indexing on the GPU | |||
|
|
|||
| include("glue.jl") | |||
There was a problem hiding this comment.
You could also wrap the glue potential tests in a testset, so that it is more clear to identify them.
@jgreener64 I think it would also be a good idea to separate more of the tests in separate files and maybe use SafeTestsets.jl to ensure that we don't have variable definitions leaking from one test to another.
There was a problem hiding this comment.
Yes probably a good idea.
| return Upair(r, pij.c, pij.c₀, pij.c₁, pij.c₂) | ||
| end | ||
| end |
Project.toml
Outdated
| ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" | ||
| Reexport = "189a3867-3050-52da-a836-e630ba90ab69" | ||
| Requires = "ae029012-a4dd-5104-9daa-d747884805df" | ||
| SingleCrystal = "3c6eccdf-2a89-4c24-a1d4-ff210daa8476" |
There was a problem hiding this comment.
Why is SingleCrystal added as a dependency of the package? Should it not be just a test dependency?
There was a problem hiding this comment.
Yep, thanks, should be fixed, as soon as I push my changes.
src/interactions/glue_fs.jl
Outdated
| params::DataFrame | ||
| pairs::Dict{String,FinnisSinclairPair} | ||
| singles::Dict{String,FinnisSinclairSingle} | ||
| kb::Real |
There was a problem hiding this comment.
You should either specify the type or parametrize it (which would be more general). Using abstract types in the fields of a struct can lead to performance problems. See https://docs.julialang.org/en/v1/manual/performance-tips/#Avoid-fields-with-abstract-containers
There was a problem hiding this comment.
Thank you for the link :). Have read it but not sure my changes are correct. Current version better?^^
src/interactions/glue_fs.jl
Outdated
| A::Real | ||
| β::Real | ||
| d::Real |
There was a problem hiding this comment.
Using abstract types in the fields of a struct can lead to performance problems. You could parametrize the struct to solve this.
src/interactions/glue_fs.jl
Outdated
| c::Real | ||
| c₀::Real | ||
| c₁::Real | ||
| c₂::Real |
There was a problem hiding this comment.
Real is an abstract type and thus your structs will not be concrete and that will lead to performance problems.
| Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" | ||
| Molly = "aa0f7f06-fcc0-5ec4-a7f3-a573f33f9c4c" | ||
| Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" |
There was a problem hiding this comment.
The plots in the documentation could be generated during the docs build, but it would be ideal to use the same package for all the plots (preferably one of the Makie backends). @jgreener64 One idea would be to use CairoMakie, but it might also be possible to use WGLMakie (it would be a bit more complex, but we get interactive plots). For this to work, the visualize function should be available with AbstractPlotting. A more elegant solution would be to use recipes, but AbstractPlotting is not a light dependency right now.
There was a problem hiding this comment.
It would indeed be better to use the same plotting package throughout, though at the minute the code in the docs is not run during the docs build and Makie plots are pasted as images into the doc files.
It would be good if the docs notebook appeared in the docs - perhaps it could go in the examples section as markdown. Then these docs dependencies can be removed until we deal with doctests later, as the examples can assume that users have the package installed.
There was a problem hiding this comment.
Tried using Makie before but had some problem I can't remember. Plots is mostly required for the docs/fs.ipynb notebook. Can try to replace Plots there though if required.
docs/Project.toml
Outdated
| @@ -1,6 +1,10 @@ | |||
| [deps] | |||
| Crystal = "3c6eccdf-2a89-4c24-a1d4-ff210daa8476" | |||
There was a problem hiding this comment.
Shouldn't this be SingleCrystal?
| Crystal = "3c6eccdf-2a89-4c24-a1d4-ff210daa8476" | |
| SingleCrystal = "3c6eccdf-2a89-4c24-a1d4-ff210daa8476" |
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
src/interaction/glue_fs.jlNote: 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
fs.ipynbto document the usage of the added interaction typeaccelerationsinsrc/forces.jl: to enable glue calculation before pair contributions to the forcespotential_energyinsrc/analysis.jlto enable the use of apotential_energyfunction for the glue interaction typesrc/loggers.jlby adding new loggers which proved useful for debugging:GlueDensityLogger,VelocityLogger,ForceLoggerSimulationinsrc/types.jlby addingglue_densitiesas an array where per-atom glue values are stored for the glue potentialtest/glue.jlas 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.jlto 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: