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

eschmidt42 · 2021-03-21T09:41:42Z

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).

…stals

… vacancy

…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

…gy tests

…o the core components

Merged small changes removing commented out lines

Glue

SebastianM-C

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.

SebastianM-C · 2021-03-21T10:14:00Z

Project.toml

 KernelDensity = "5ab0869b-81aa-558d-bb23-cbf5423bbe9b"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"


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.

SebastianM-C · 2021-03-21T10:20:32Z

src/interactions/glue_fs.jl

+end
+
+
+kb = 8.617333262145e-5 # eV/K


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

SebastianM-C · 2021-03-21T11:26:42Z

Project.toml

@@ -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"


If you need some packages for tests, you should add them as test dependencies in the [extras] section, not as direct package dependencies.

SebastianM-C · 2021-03-21T11:29:33Z

src/interactions/glue_fs.jl

+    element_pairings = [string(el, el) for el in elements]
+    element_pair_map = Dict(pair => i for (i,pair) in enumerate(element_pairings))
+
+    df = DataFrame(


I'm not sure if adding the DataFrames.jl dependency if worth it just for this convenience functionality.

SebastianM-C · 2021-03-21T11:31:13Z

src/interactions/glue_fs.jl

+
+Derivative of the glue density function.
+"""
+∂glue_∂r(r, β, d) = ForwardDiff.derivative(r -> glue(r,β,d), r)


In this case I think that the derivative could be done analytically since the formula is not very complicated.

SebastianM-C · 2021-03-21T11:50:41Z

Project.toml

 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"


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.

SebastianM-C · 2021-03-21T12:01:25Z

src/loggers.jl

+
+function ForceLogger(T, n_steps::Integer; dims::Integer=3)
+    return ForceLogger(n_steps,
+                            Array{SArray{Tuple{dims}, T, 1, dims}, 1}[])


You could just use Vector{SVector{dims}} here for readability.

Suggested change

Array{SArray{Tuple{dims}, T, 1, dims}, 1}[])

Vector{SVector{dims}}[])

SebastianM-C · 2021-03-21T12:06:50Z

src/interactions/glue_fs.jl

+    element_pair_map::Dict
+    params::DataFrame


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

SebastianM-C · 2021-03-21T12:10:58Z

src/interactions/glue_fs.jl

+
+Energy derivative given glue density.
+"""
+∂Uglue_∂ρ(ρ,A) = ForwardDiff.derivative(ρ -> Uglue(ρ,A), ρ)


This could be simply derived analytically

Suggested change

∂Uglue_∂ρ(ρ,A) = ForwardDiff.derivative(ρ -> Uglue(ρ,A), ρ)

∂Uglue_∂ρ(ρ,A) = -A /(2 * √ρ)

SebastianM-C · 2021-03-21T12:13:11Z

src/forces.jl

    end
-
+    


Remove whitespace

Suggested change

eschmidt42 · 2021-03-22T18:02:26Z

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 · 2021-03-24T12:07:45Z

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.

…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

SebastianM-C · 2021-05-11T08:51:06Z

test/runtests.jl

@@ -22,6 +22,8 @@ end

 CUDA.allowscalar(false) # Check that we never do scalar indexing on the GPU

+include("glue.jl")


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.

Yes probably a good idea.

SebastianM-C · 2021-05-11T08:51:28Z

src/interactions/glue_fs.jl

    return Upair(r, pij.c, pij.c₀, pij.c₁, pij.c₂)
-end
+end


Suggested change

end

end

SebastianM-C · 2021-05-11T08:56:18Z

Project.toml

 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
+SingleCrystal = "3c6eccdf-2a89-4c24-a1d4-ff210daa8476"


Why is SingleCrystal added as a dependency of the package? Should it not be just a test dependency?

Yep, thanks, should be fixed, as soon as I push my changes.

SebastianM-C · 2021-05-11T09:03:20Z

src/interactions/glue_fs.jl

-    params::DataFrame
+    pairs::Dict{String,FinnisSinclairPair}
+    singles::Dict{String,FinnisSinclairSingle}
+    kb::Real


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

Thank you for the link :). Have read it but not sure my changes are correct. Current version better?^^

SebastianM-C · 2021-05-11T09:06:09Z

src/interactions/glue_fs.jl

+    A::Real
+    β::Real
+    d::Real


Using abstract types in the fields of a struct can lead to performance problems. You could parametrize the struct to solve this.

SebastianM-C · 2021-05-11T09:08:35Z

src/interactions/glue_fs.jl

+    c::Real
+    c₀::Real
+    c₁::Real
+    c₂::Real


Real is an abstract type and thus your structs will not be concrete and that will lead to performance problems.

SebastianM-C · 2021-05-11T09:18:20Z

docs/Project.toml

 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 Molly = "aa0f7f06-fcc0-5ec4-a7f3-a573f33f9c4c"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"


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.

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.

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.

SebastianM-C · 2021-05-11T09:19:09Z

docs/Project.toml

@@ -1,6 +1,10 @@
 [deps]
+Crystal = "3c6eccdf-2a89-4c24-a1d4-ff210daa8476"


Shouldn't this be SingleCrystal?

Suggested change

Crystal = "3c6eccdf-2a89-4c24-a1d4-ff210daa8476"

SingleCrystal = "3c6eccdf-2a89-4c24-a1d4-ff210daa8476"

eschmidt42 added 27 commits February 2, 2021 19:34

added packages

19f43d6

added the finnis-sinclair variant of glue type potentials

a3f4b6a

added documentation

5b805ff

expanded potential energy and force tests to all vacancy free bcc cry…

8942490

…stals

added potential energy and force checks for bcc crystal with a single…

87e7d82

… vacancy

added tests for the finnis-sinclair type of glue potentials

52f46b8

added section to run a longer simulation for the documentation

dd0881c

added fields for forces, velocities and glue densities in Simulation

221b690

refactored glue forces to work as general_inters

eaacf33

modified potential_energy to pre-compute glue densities, if necessary

a9d7ef1

modified accelerations to allow for glue density update and calculati…

fd93544

…on of glue forces using the general_inters loop

added GlueDensityLogger, VelocityLogger and ForcesLogger

95b9635

fixed copy-paste error

36457e0

added ForwardDiff

1808489

refactored reducing redundancy and replacing manually implemented der…

057b9a6

…ivates with FordwardDiff.derivate to re-use the non-derivative forms also used to successfully test ground state energies

updated algebraic signs

acf7ce8

added convenience function to easily generate variables for the Finni…

b99aca3

…s-Sinclair 1984 empirical model

refactored the ground state energy, forces and vacancy formation ener…

fd7cb6f

…gy tests

refactored fs.ipynb: removed superfluous tests and reduced the math t…

f80004e

…o the core components

removed superfluous enumerate

275c296

added docstrings

b136903

cleaned

931f214

removed cell

4e2c73b

removing comments

f88c5b2

Merge branch 'hotfix-removing-comments'

3982887

Merged small changes removing commented out lines

Merge pull request #1 from eschmidt42/glue

b472fd9

Glue

resolved merge conflict

ee4d487

SebastianM-C suggested changes Mar 21, 2021

View reviewed changes

eschmidt42 added 19 commits March 27, 2021 19:34

replaced ForwardDiff.derivate with analytical derivatives and added k…

5a81002

…b as an attribute of FinnisSinclair

removed f comment in docstrings

c0d37fa

replaced dataframes with dicts

a71dc51

replaced another derivative forgotten earlier

4b55ae4

updated tests

b332f60

minor change to improve readability of ForceLogger

825f1ad

removed DataFrames reference

d646f9b

removed DataFrames and Plots reference

945ddb0

relocated packages from Molly.jl to docs

115adf2

made notebook compatible with changes in in src/interactions/glue_fs.jl

30a54b2

removed because relocated to ./docs

331cc70

removed superfluous cells

d87b4c6

adjusted glue tests to changes in Crystal.jl and introduced more name…

bbc47f1

…d testsets for legitibility of the test summary in case of failures

updated test: switched from Crystal to SingleCrystal package

d073ef8

switched toml from Crystal to SingleCrystal

fbf2629

switched to toml SingleCrystal.jl

d95cfbe

updated notebook to work with the new data structure of the potential…

b276156

… and changes in the package constructing the bcc crystal

removed ForwardDiff

71522fe

removed empty spaces and updated doc strings

b9ece8d

SebastianM-C reviewed May 11, 2021

View reviewed changes

eschmidt42 added 7 commits May 14, 2021 19:48

rm SingleCrystal package as direct dependency

d1f113a

replaced abstract types in structs

05d1990

wrapped glue tests in a testset

12c5a31

kind of resolved conflict

996ed6e

welp...

3b9ce51

adjusted to renaming of Simnulation struct

076ad3d

updated test to new Simulation

d2b3375

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

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

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

eschmidt42 commented Mar 21, 2021

SebastianM-C left a comment

SebastianM-C Mar 21, 2021

SebastianM-C Mar 21, 2021

SebastianM-C Mar 21, 2021

SebastianM-C Mar 21, 2021

SebastianM-C Mar 21, 2021

SebastianM-C Mar 21, 2021

SebastianM-C Mar 21, 2021

SebastianM-C Mar 21, 2021

SebastianM-C Mar 21, 2021

SebastianM-C Mar 21, 2021

eschmidt42 commented Mar 22, 2021

jgreener64 commented Mar 24, 2021

SebastianM-C May 11, 2021

jgreener64 May 12, 2021

SebastianM-C May 11, 2021

SebastianM-C May 11, 2021

eschmidt42 May 14, 2021 •

edited

Loading

SebastianM-C May 11, 2021

eschmidt42 May 14, 2021

SebastianM-C May 11, 2021

SebastianM-C May 11, 2021

SebastianM-C May 11, 2021

jgreener64 May 12, 2021

eschmidt42 May 14, 2021

SebastianM-C May 11, 2021

		Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
		ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"

	Array{SArray{Tuple{dims}, T, 1, dims}, 1}[])
	Vector{SVector{dims}}[])

	∂Uglue_∂ρ(ρ,A) = ForwardDiff.derivative(ρ -> Uglue(ρ,A), ρ)
	∂Uglue_∂ρ(ρ,A) = -A /(2 * √ρ)

		@@ -22,6 +22,8 @@ end

		CUDA.allowscalar(false) # Check that we never do scalar indexing on the GPU

		include("glue.jl")

		@@ -1,6 +1,10 @@
		[deps]
		Crystal = "3c6eccdf-2a89-4c24-a1d4-ff210daa8476"

	Crystal = "3c6eccdf-2a89-4c24-a1d4-ff210daa8476"
	SingleCrystal = "3c6eccdf-2a89-4c24-a1d4-ff210daa8476"

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

Are you sure you want to change the base?

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

Conversation

eschmidt42 commented Mar 21, 2021

Glue: Finnis-Sinclair Potential

Main Change

Minor Changes

Next steps

SebastianM-C left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

eschmidt42 commented Mar 22, 2021

jgreener64 commented Mar 24, 2021

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

eschmidt42 May 14, 2021 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

eschmidt42 May 14, 2021 •

edited

Loading