From 19f43d66a3c7234a0cb9e9cd05401930dbe28d6a Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Tue, 2 Feb 2021 19:34:47 +0100 Subject: [PATCH 01/49] added packages --- Project.toml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/Project.toml b/Project.toml index 430d0dcf..6eff02af 100644 --- a/Project.toml +++ b/Project.toml @@ -7,9 +7,12 @@ version = "0.2.1" BioStructures = "de9282ab-8554-53be-b2d6-f6c222edabfc" CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" +Crystal = "3c6eccdf-2a89-4c24-a1d4-ff210daa8476" +DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" KernelDensity = "5ab0869b-81aa-558d-bb23-cbf5423bbe9b" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" Requires = "ae029012-a4dd-5104-9daa-d747884805df" From a3f4b6a293aa7107a915371e889784b227634e5f Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Tue, 2 Feb 2021 19:35:52 +0100 Subject: [PATCH 02/49] added the finnis-sinclair variant of glue type potentials --- fs.ipynb | 771 ++++++++++++++++++++++++++++++++++++ src/forces.jl | 10 +- src/interactions/glue_fs.jl | 168 ++++++++ 3 files changed, 947 insertions(+), 2 deletions(-) create mode 100644 fs.ipynb create mode 100644 src/interactions/glue_fs.jl diff --git a/fs.ipynb b/fs.ipynb new file mode 100644 index 00000000..3b4047b3 --- /dev/null +++ b/fs.ipynb @@ -0,0 +1,771 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# The Finnis-Sinclair potential" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Structure:\n", + "1. Computing Finnis-Sinclair forces\n", + "2. Computing Finnis-Sinclair potential energies" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "TODOs:\n", + "* test forces for all potentials using perfect crystals and crystals with a vacancy\n", + "* add tests for potential energy values for each structure\n", + "* other force implementation that does not require altering `function accelerations`?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import Pkg" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "Pkg.activate(\".\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Pkg.add(url=\"https://github.com/eschmidt42/crystal\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "Pkg.add(\"Plots\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "Pkg.status()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "using Molly" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "using DataFrames\n", + "# using Molly\n", + "using Plots\n", + "using Test\n", + "# using LaTeXStrings\n", + "# using LinearAlgebra\n", + "# using SparseArrays\n", + "using Crystal" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Parameterisation of V, Nb, Ta, Cr, Mo, W, Fe" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Parameterisation by Finnis et al. 1984, _A simple empirical N-body potential for transition metals_" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "| element | d | A | $\\beta$ | c | $c_0$ | $c_1$ | $c_2$ | \n", + "| --- | --- | --- | --- | --- | --- | --- | --- |\n", + "| V | 3.692767 | 2.010637 | 0 | 3.8 | -0.8816318 | 1.4907756 | -0.3976370 |\n", + "| Nb | 3.915354 | 3.013789 | 0 | 4.2 | -1.5640104 | 2.0055779 | -0.4663764 |\n", + "| Ta | 4.076980 | 2.591061 | 0 | 4.2 | 1.2157373 | 0.0271471 | -0.1217350 |\n", + "| Cr | 3.915720 | 1.453418 | 1.8 | 2.9 | 29.1429813 | -23.3975027 | 4.7578297 |\n", + "| Mo | 4.114825 | 1.887117 | 0 | 3.25 | 43.4475218 | -31.9332978 | 6.0804249 |\n", + "| W | 4.400224 | 1.896373 | 0 | 3.25 | 47.1346499 | -33.7665655 | 6.2541999 |\n", + "| Fe | 3.699579 | 1.889846 | 1.8 | 3.4 | 1.2110601 | -0.7510840 | 0.1380773 |" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "elements = [\"V\", \"Nb\", \"Ta\", \"Cr\", \"Mo\", \"W\", \"Fe\"]\n", + "element_pairings = [string(el,el) for el in elements]\n", + "element_pair_map = Dict(pair => i for (i,pair) in enumerate(element_pairings))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "df = DataFrame(\n", + " element_pair = element_pairings,\n", + " d = [3.692767, 3.915354, 4.076980, 3.915720, 4.114825, 4.400224, 3.699579],\n", + " A = [2.010637, 3.013789, 2.591061, 1.453418, 1.887117, 1.896373, 1.889846],\n", + " β = [0, 0, 0, 1.8, 0, 0, 1.8],\n", + " c = [3.8, 4.2, 4.2, 2.9, 3.25, 3.25, 3.4],\n", + " c₀ = [-0.8816318, -1.5640104, 1.2157373, 29.1429813, 43.4475218, 47.1346499, 1.2110601],\n", + " c₁ = [1.4907756, 2.0055779, 0.0271471, -23.3975027, -31.9332978, -33.7665655, -0.7510840],\n", + " c₂ = [-0.3976370, -0.4663764, -0.1217350, 4.7578297, 6.0804249, 6.2541999, 0.1380773],\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "masses = Dict(\"V\" => 50.9415, \"Nb\" => 92.9064, \"Ta\" => 180.9479,\n", + " \"Cr\" => 51.996, \"Mo\" => 95.94, \"W\" => 183.85,\n", + " \"Fe\" => 55.847)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Å\n", + "bcc_lattice_constants = Dict(\n", + " \"V\" => 3.0399, \"Nb\" => 3.3008, \n", + " \"Ta\" => 3.3058, \"Cr\" => 2.8845, \"Mo\" => 3.1472, \n", + " \"W\" => 3.1652, \"Fe\" => 2.8665\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Interaction" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Instantiating the interaction" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fs_inter = FinnisSinclair(true, element_pair_map, df)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Glue potential" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The glue potential is the core component which makes the Finnis-Sinclair empirical potential and other similar approaches different to, for example, the Lennard-Jones potential. \n", + "TODO: some more explanation" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\phi(r) = (r-d)^2 + \\beta (r-d)^3/d\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "element_pair, d, A, β, c, c₀, c₁, c₂ = df[1,:] # parameters for Vanadium" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "r = collect(range(0, stop=2*3.3058, length=1000));\n", + "ɸ = Molly.glue_potential.(r, β, d);" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ɸs = [ɸ]\n", + "element_pairs = [element_pair]\n", + "for i in 2:nrow(df)\n", + " element_pair, d, A, β, c, c₀, c₁, c₂ = df[i,:]\n", + " ɸ = Molly.glue_potential.(r, β, d)\n", + " append!(ɸs,[ɸ])\n", + " element_pairs = hcat(element_pairs, string(element_pair))\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plot(r, ɸs, label=element_pairs)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Computing forces" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This is what we need for the simulation." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$\\partial_k = \\partial_{R_k} = $ change of atom position $k$, $r_{ij} = \\|R_{ij}\\|_2$, $R_{ij} = R_j - R_i \\in \\mathbb{R}^3$ \n", + "\n", + "$$\n", + "\\partial_k u = \\partial_k u_\\text{pair} + \\partial_k u_\\text{glue} \n", + "$$\n", + "\n", + "$$\n", + "\\partial_k u_\\text{pair} = \\sum_{i>j} V_{ij}^\\prime(r_{ij})\\partial_k r_{ij}\n", + "$$\n", + "\n", + "$$\n", + "\\partial_k u_\\text{glue} = \\sum_i f_i^\\prime(\\rho_i) \\cdot \\partial_k \\rho_i \n", + "$$\n", + "\n", + "$$\n", + "V_{ij}^\\prime(r_{ij}) = 2(r-c)(c_0 + c_1 r + c_2 r^2) + (r-c)^2(c_1 + 2c_2r)\n", + "$$\n", + "\n", + "$$\n", + "f_i^\\prime(\\rho_i) \\cdot \\partial_k \\rho_i =\n", + "\\begin{cases}\n", + " k = i, & f_k^\\prime(\\rho_k) \\sum_j\\phi_j^\\prime(r_{kj})\\partial_k r_{kj} \\\\\n", + " k \\ne i, & \\sum_{i\\ne k} f_i^\\prime(\\rho_i) \\partial_k \\phi_k^\\prime(r_{ik})\\partial_k r_{ik} \\\\\n", + "\\end{cases}\n", + "$$\n", + "\n", + "$$\n", + "f_i^\\prime = \\frac{1}{2}A_i\\rho_i^{-3/2}\n", + "$$\n", + "\n", + "$$\n", + "\\partial_k\\phi(r) = \\left[2(r-d) + 3\\beta (r-d)^2/d\\right] \\cdot\n", + "\\begin{cases}\n", + " k = i, &\\frac{R_{kj}}{r_{kj}} \\\\ \n", + " k = j, &\\frac{R_{ik}}{r_{ik}} \\\\\n", + "\\end{cases}\n", + "$$\n", + "\n", + "$$\n", + "\\partial_k r_{ij} =\n", + "\\begin{cases}\n", + " k = i, &\\frac{R_{kj}}{r_{kj}} \\\\ \n", + " k = j, &\\frac{R_{ik}}{r_{ik}} \\\\\n", + "\\end{cases}\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "element = \"Fe\"\n", + "a = bcc_lattice_constants[element]\n", + "atoms, coords, box, box_size, box_vectors = Crystal.make_bcc_unitcell(element, a=a)\n", + "sc_atoms, sc_coords, sc_box, sc_box_size = Crystal.make_supercell(atoms, coords, box, box_size, nx=3, ny=3,nz=3)\n", + "n_atoms = length(sc_atoms);" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "specific_inter_list = ((fs_inter,),)\n", + "velocities = [velocity(1., .01, dims=3) for i in 1:n_atoms]\n", + "sim = VelocityVerlet()\n", + "nb_matrix = trues(n_atoms,n_atoms)\n", + "n_steps = 1\n", + "dist_cutoff = 2 * a\n", + "\n", + "nf = DistanceNeighbourFinder(nb_matrix, n_steps, dist_cutoff);\n", + "\n", + "loggers = Dict(\"temperature\" => TemperatureLogger(1))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "s = Simulation(\n", + " simulator=sim, \n", + " atoms=sc_atoms, \n", + " specific_inter_lists=specific_inter_list,\n", + " general_inters=(),\n", + " coords=[SVector{3}(v) for v in sc_coords], \n", + " velocities=velocities,\n", + " temperature=.01, \n", + " box_size=sc_box_size[1,1],\n", + " timestep=.002,\n", + " n_steps=10,\n", + " neighbour_finder=nf,\n", + " loggers=loggers,\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sparse_forces = force.((fs_inter,), (s.coords,), (s,))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Testing that all forces are about 0" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "function test_sparse_forces_zero(sparse_forces, n_atoms; dims=3)\n", + " zeros = [zero(rand(1,3)) for _ in 1:n_atoms]\n", + " forces = getindex.(sparse_forces,2)[1]\n", + " return all(isapprox.(forces, zeros, atol=1e-6))\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "@test test_sparse_forces_zero(sparse_forces, n_atoms)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Modifying `Molly.accelerations` so the forces from the glue interaction are properly used to update the atom forces." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "function test_forces_zero(forces, n_atoms; dims=3)\n", + " zeros = [zero(rand(3)) for _ in 1:n_atoms]\n", + " return all(isapprox.(forces, zeros, atol=1e-6))\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "@test test_forces_zero(accelerations(s, parallel=false), n_atoms)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "simulate!(s, parallel=false)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "s.loggers[\"temperature\"].temperatures" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Computing energies" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This is only really interesting for logging / development of potentials." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Pair energy" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "V_{ij}(r_{ij}) = \n", + "\\begin{cases} \n", + "r \\le c, & (r-c)^2 \\left( c_0 + c_1 r + c_2 r^2 \\right) \\\\\n", + "r > c, & 0 \\\\\n", + "\\end{cases}\n", + "$$\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "element_pair, d, A, β, c, c₀, c₁, c₂ = df[1,:] # parameters for Vanadium" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "V = Molly.pair_potential.(r, c, c₀, c₁, c₂);" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "Vs = [V]\n", + "element_pairs = [element_pair]\n", + "for i in 2:nrow(df)\n", + " element_pair, d, A, β, c, c₀, c₁, c₂ = df[i,:]\n", + " V = Molly.pair_potential.(r, c, c₀, c₁, c₂)\n", + " append!(Vs,[V])\n", + " element_pairs = hcat(element_pairs, string(element_pair))\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plot(r, Vs, label=element_pairs)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Glue energy" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Computing an energy based on local glue values" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "u_\\text{glue} = -A \\cdot \\sqrt{\\rho}\n", + "$$\n", + "\n", + "$$\n", + "\\rho = \\sum_{j \\in \\text{neighborhood}(i)} \\phi(r_{ij})\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ρ = 4. # density that you get summing phi-contributions from neighbours\n", + "Molly.glue_energy(ρ, 1.)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ρ = collect(range(0, stop=50, length=100));\n", + "\n", + "A = df.A[1] # Va\n", + "uₙ = Molly.glue_energy.(ρ, A)\n", + "element_pair = df.element_pair[1]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "uₙs = [uₙ]\n", + "element_pairs = [element_pair]\n", + "for i in 2:nrow(df)\n", + " element_pair, d, A, β, c, c₀, c₁, c₂ = df[i,:]\n", + " uₙ = Molly.glue_energy.(ρ, A)\n", + " append!(uₙs,[uₙ])\n", + " element_pairs = hcat(element_pairs, string(element_pair))\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plot(ρ, uₙs, label=element_pairs)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Pair + glue energy = magic" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$ u_\\text{tot} = u_N + u_P $$\n", + "\n", + "$$ u_P = \\frac{1}{2}\\sum_{i=1,j=1}^{n_\\text{atoms},n_\\text{atoms}} V(r_{ij}) $$\n", + "\n", + "$$ u_N = \\sum_{i=1}^{n_\\text{atoms}} u_\\text{glue}(\\rho_i) $$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "u = Molly.potential_energy(fs_inter, s) / n_atoms" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "@test isapprox(u, -4.28, atol=1e-2)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Running simulation with potential energy logger" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "loggers = Dict(\n", + " \"temperature\" => TemperatureLogger(1),\n", + " \"energy\" => EnergyLogger(1),\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "s = Simulation(\n", + " simulator=sim, \n", + " atoms=sc_atoms, \n", + " specific_inter_lists=specific_inter_list,\n", + " general_inters=(),\n", + " coords=[SVector{3}(v) for v in sc_coords], \n", + " velocities=velocities,\n", + " temperature=.01, \n", + " box_size=sc_box_size[1,1],\n", + " timestep=.002,\n", + " n_steps=10,\n", + " neighbour_finder=nf,\n", + " loggers=loggers,\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "simulate!(s, parallel=false)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "s.loggers[\"energy\"].energies" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 1.5.3", + "language": "julia", + "name": "julia-1.5" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.5.3" + }, + "toc": { + "base_numbering": 1, + "nav_menu": { + "height": "167px", + "width": "309px" + }, + "number_sections": true, + "sideBar": true, + "skip_h1_title": false, + "title_cell": "Table of Contents", + "title_sidebar": "Contents", + "toc_cell": false, + "toc_position": {}, + "toc_section_display": true, + "toc_window_display": false + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/src/forces.jl b/src/forces.jl index bcf09e95..080b5c3f 100644 --- a/src/forces.jl +++ b/src/forces.jl @@ -12,7 +12,8 @@ export Gravity, HarmonicBond, HarmonicAngle, - Torsion + Torsion, + FinnisSinclair """ force(inter, coord_i, coord_j, atom_i, atom_j, box_size) @@ -92,7 +93,11 @@ function accelerations(s::Simulation; parallel::Bool=true) sparse_forces = force.(inter_list, (s.coords,), (s,)) ge1, ge2 = getindex.(sparse_forces, 1), getindex.(sparse_forces, 2) sparse_vec = SparseVector(n_atoms, reduce(vcat, ge1), reduce(vcat, ge2)) - forces += Array(sparse_vec) + if typeof(inter_list[1]) <: GlueInteraction + forces += Array([SVector{3}(v) for v in sparse_vec]) + else + forces += Array(sparse_vec) + end end for i in 1:n_atoms @@ -134,3 +139,4 @@ include("interactions/gravity.jl") include("interactions/harmonic_bond.jl") include("interactions/harmonic_angle.jl") include("interactions/torsion.jl") +include("interactions/glue_fs.jl") diff --git a/src/interactions/glue_fs.jl b/src/interactions/glue_fs.jl new file mode 100644 index 00000000..ac02042a --- /dev/null +++ b/src/interactions/glue_fs.jl @@ -0,0 +1,168 @@ +""" +A glue interaction that will apply to all atom pairs. +Finnis-Sinclair and similar interactions should sub-type this type. +This type should be a GeneralInteraction type. But due to the special +nature of glue interactions and the restriction to pair interactions +of the GeneralInteraction type, glue interactions are for now a sub-type of +SpecificInteraction. +""" +abstract type GlueInteraction <: SpecificInteraction end + +""" + FinnisSinclairInteraction(nl_only,element_pair_map,params) + +The Finnis-Sinclair interaction. +""" +struct FinnisSinclair <: GlueInteraction + nl_only::Bool + element_pair_map::Dict + params::DataFrame +end + +""" + glue_potential(r,β,d) + +The core component of the Finnis-Sinclair type GlueInteraction. Used to calculate contribution to the glue value of an atom, based on its neighbouring atoms. +""" +function glue_potential(r::T, β::T, d::T)::T where T<:Real + return r > d ? 0 : (r-d)^2 + β*(r-d)^3/d +end + +""" +""" +function pair_potential_derivative(r::T, c::T, c₀::T, c₁::T, c₂::T)::T where T<:Real + return (r > c) ? 0 : 2 * (r - c) * (c₀ + c₁*r + c₂*r^2) + (r - c)^2 * (c₁ + 2*c₂*r) +end + +""" +""" +function glue_energy_derivative(ρ::Float64, A::Float64)::Float64 + return A/2 * ρ^-1.5 +end + +""" +""" +function glue_potential_derivative(r::T, β::T, d::T)::T where T<:Real + return r > d ? 0 : 2*(r-d) + 3*β*(r-d)^2/d +end + +function get_pair_params(element1::String, element2::String, inter::FinnisSinclair) + pair = string(sort([element1, element2])...) + return inter.params[inter.element_pair_map[pair],:] +end + +@inline @inbounds function force( + inter::FinnisSinclair, + coords, + s::Simulation + ) + # computing the embedding densities + n_atoms = length(s.coords) + ρs = zeros(n_atoms) + rs = zeros(length(s.neighbours)) + r_vec_norms = zeros(length(s.neighbours),3) + + for (n,(i,j)) in enumerate(s.neighbours) + element_i = s.atoms[i].name + element_j = s.atoms[j].name + element_pair = string(sort([element_i, element_j])...) + pi = get_pair_params(element_i,element_i,inter) # inter.params[inter.element_map[element_i],:] + pj = get_pair_params(element_j,element_j,inter) # inter.params[inter.element_map[element_j],:] + pij = get_pair_params(element_i,element_j,inter) # inter.params[inter.element_map[element_pair],:] + + r_vec = vector(s.coords[i], s.coords[j], s.box_size) + r2 = sum(abs2, r_vec) + r = sqrt(r2) + # storing distance (vectors) so we don't need to recompute + rs[n] = r + r_vec_norms[[n],:] = r_vec / r + # storing glue densities + ρs[i] += glue_potential(r, pj.β, pj.d) + ρs[j] += glue_potential(r, pi.β, pi.d) + end + + fs = [zeros(1,3) for _ in 1:n_atoms] + for (n,(i,j)) in enumerate(s.neighbours) + element_i = s.atoms[i].name + element_j = s.atoms[j].name + element_pair = string(sort([element_i, element_j])...) + pi = get_pair_params(element_i,element_i,inter) # inter.params[inter.element_map[element_i],:] + pj = get_pair_params(element_j,element_j,inter) # inter.params[inter.element_map[element_j],:] + pij = get_pair_params(element_i,element_j,inter) # inter.params[inter.element_map[element_pair],:] + + r = rs[n] + r2 = r^2 + r_vec_norm = r_vec_norms[[n],:] + + # pair contribution + dpairdR_i = r_vec_norm * pair_potential_derivative(r, pij.c, pij.c₀, pij.c₁, pij.c₂) + dpairdR_j = - dpairdR_i + + # glue contribution + dudρ_i = glue_energy_derivative(ρs[i], pi.A) + dudρ_j = glue_energy_derivative(ρs[j], pj.A) + dΦdr_i = glue_potential_derivative(r, pi.β, pi.d) + dΦdr_j = glue_potential_derivative(r, pj.β, pj.d) + + ## density change by moving the current atom + dgluedR_i_curr = r_vec_norm * dudρ_i * dΦdr_j + dgluedR_j_curr = r_vec_norm * dudρ_j * dΦdr_i + ## density change by moving a neighbouring atom + dgluedR_i_neigh = - r_vec_norm * dudρ_j * dΦdr_i + dgluedR_j_neigh = - r_vec_norm * dudρ_i * dΦdr_j + + # updating the forces + f_i = (dpairdR_i + dgluedR_i_curr + dgluedR_i_neigh) + f_j = (dpairdR_j + dgluedR_j_curr + dgluedR_j_neigh) + fs[i] += f_i + fs[j] += f_j + end + + return collect(1:n_atoms), fs +end + +""" +""" +function pair_potential(r::T, c::T, c₀::T, c₁::T, c₂::T)::T where T<:Real + return (r > c) ? 0 : (r - c)^2 * (c₀ + c₁*r + c₂*r^2) +end + +""" +""" +function glue_energy(ρ::Float64, A::Float64)::Float64 + return -A * √ρ +end + +@inline @inbounds function potential_energy(inter::FinnisSinclair, s::Simulation) + #computing the potential energy combining glue and pair components. + + e_pair = 0. + e_glue = 0. + n_atoms = length(s.coords) + ρs = zeros(n_atoms) + for (i,j) in s.neighbours + element_i = s.atoms[i].name + element_j = s.atoms[j].name + element_pair = string(sort([element_i, element_j])...) + pi = get_pair_params(element_i,element_i,inter) # inter.params[inter.element_map[element_i],:] + pj = get_pair_params(element_j,element_j,inter) # inter.params[inter.element_map[element_j],:] + pij = get_pair_params(element_i,element_j,inter) # inter.params[inter.element_map[element_pair],:] + + r_vec = vector(s.coords[i], s.coords[j], s.box_size) + r2 = sum(abs2, r_vec) + r = sqrt(r2) + + e_pair += pair_potential(r, pij.c, pij.c₀, pij.c₁, pij.c₂) + + ρs[i] += glue_potential(r, pj.β, pj.d) + ρs[j] += glue_potential(r, pi.β, pi.d) + end + + es_glue = zeros(n_atoms) + for (i, atom) in enumerate(s.atoms) + A = get_pair_params(atom.name, atom.name, inter).A + es_glue[i] = glue_energy(ρs[i], A) + end + e_glue = sum(es_glue) + return e_pair + e_glue +end \ No newline at end of file From 5b805ffb737b273977eea8aa5f660ded1b79b16c Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Wed, 3 Feb 2021 18:36:00 +0100 Subject: [PATCH 03/49] added documentation --- src/interactions/glue_fs.jl | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/src/interactions/glue_fs.jl b/src/interactions/glue_fs.jl index ac02042a..63b631c9 100644 --- a/src/interactions/glue_fs.jl +++ b/src/interactions/glue_fs.jl @@ -29,23 +29,37 @@ function glue_potential(r::T, β::T, d::T)::T where T<:Real end """ + pair_potential_derivative(r, c, c₀, c₁, c₂) + +Derivative of the pair potential. """ function pair_potential_derivative(r::T, c::T, c₀::T, c₁::T, c₂::T)::T where T<:Real return (r > c) ? 0 : 2 * (r - c) * (c₀ + c₁*r + c₂*r^2) + (r - c)^2 * (c₁ + 2*c₂*r) end """ + glue_energy_derivative(ρ, A) + +Energy derivative given glue density. """ function glue_energy_derivative(ρ::Float64, A::Float64)::Float64 return A/2 * ρ^-1.5 end """ + glue_potential_derivative(r, β, d) + +Derivative of the glue density function. """ function glue_potential_derivative(r::T, β::T, d::T)::T where T<:Real return r > d ? 0 : 2*(r-d) + 3*β*(r-d)^2/d end +""" + get_pair_params(element1, element2, inter) + +Convenience function to generate element pair and return relevant model parameters. +""" function get_pair_params(element1::String, element2::String, inter::FinnisSinclair) pair = string(sort([element1, element2])...) return inter.params[inter.element_pair_map[pair],:] @@ -122,12 +136,18 @@ end end """ + pair_potential(r, c, c₀, c₁, c₂) + +Energy contribution directly from atom pair distances. """ function pair_potential(r::T, c::T, c₀::T, c₁::T, c₂::T)::T where T<:Real return (r > c) ? 0 : (r - c)^2 * (c₀ + c₁*r + c₂*r^2) end """ + glue_energy(ρ, A) + +Energy based on the glue density . """ function glue_energy(ρ::Float64, A::Float64)::Float64 return -A * √ρ From 8942490a05bb2242b11f74d0347e95d542242368 Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Wed, 3 Feb 2021 19:42:44 +0100 Subject: [PATCH 04/49] expanded potential energy and force tests to all vacancy free bcc crystals --- fs.ipynb | 150 +++++++++++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 146 insertions(+), 4 deletions(-) diff --git a/fs.ipynb b/fs.ipynb index 3b4047b3..8591f622 100644 --- a/fs.ipynb +++ b/fs.ipynb @@ -21,8 +21,10 @@ "metadata": {}, "source": [ "TODOs:\n", - "* test forces for all potentials using perfect crystals and crystals with a vacancy\n", - "* add tests for potential energy values for each structure\n", + "* test forces for crystals with a vacancy\n", + "* test energies for crystals with a vacancy\n", + "* copy tests to the appropriate location in the package\n", + "* if moving this nb to docs or so, does `using Molly` from non-root dirs still work?\n", "* other force implementation that does not require altering `function accelerations`?" ] }, @@ -61,7 +63,7 @@ "metadata": {}, "outputs": [], "source": [ - "Pkg.add(\"Plots\")" + "# Pkg.add(\"Plots\")" ] }, { @@ -183,6 +185,18 @@ ")" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "reference_energies = DataFrame(\n", + " element_pair = element_pairings,\n", + " u = [5.31, 7.57, 8.1, 4.1, 6.82, 8.9, 4.28],\n", + ")" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -447,7 +461,7 @@ "source": [ "function test_forces_zero(forces, n_atoms; dims=3)\n", " zeros = [zero(rand(3)) for _ in 1:n_atoms]\n", - " return all(isapprox.(forces, zeros, atol=1e-6))\n", + " return all(isapprox.(forces, zeros, atol=1e-4))\n", "end" ] }, @@ -480,6 +494,67 @@ "s.loggers[\"temperature\"].temperatures" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Testing forces for all elements" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "function test_forces_for_element(element::String, fs_inter; nx::Integer=3, ny::Integer=3, nz::Integer=3, )\n", + " a = bcc_lattice_constants[element]\n", + " atoms, coords, box, box_size, box_vectors = Crystal.make_bcc_unitcell(element, a=a)\n", + " sc_atoms, sc_coords, sc_box, sc_box_size = Crystal.make_supercell(atoms, coords, box, box_size, nx=nx, ny=ny, nz=nz)\n", + " n_atoms = length(sc_atoms)\n", + " \n", + " specific_inter_list = ((fs_inter,),)\n", + " velocities = [velocity(1., .01, dims=3) for i in 1:n_atoms]\n", + " sim = VelocityVerlet()\n", + " nb_matrix = trues(n_atoms,n_atoms)\n", + " n_steps = 1\n", + " dist_cutoff = 2 * a\n", + "\n", + " nf = DistanceNeighbourFinder(nb_matrix, n_steps, dist_cutoff);\n", + "\n", + " loggers = Dict(\"temperature\" => TemperatureLogger(1))\n", + " \n", + " s = Simulation(\n", + " simulator=sim, \n", + " atoms=sc_atoms, \n", + " specific_inter_lists=specific_inter_list,\n", + " general_inters=(),\n", + " coords=[SVector{3}(v) for v in sc_coords], \n", + " velocities=velocities,\n", + " temperature=.01, \n", + " box_size=sc_box_size[1,1],\n", + " timestep=.002,\n", + " n_steps=1,\n", + " neighbour_finder=nf,\n", + " loggers=loggers,\n", + " )\n", + " find_neighbours!(s, s.neighbour_finder, 0)\n", + " forces = accelerations(s, parallel=false)\n", + " return test_forces_zero(forces, n_atoms)\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "for element in elements\n", + " @test test_forces_for_element(element, fs_inter)\n", + "end" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -668,6 +743,73 @@ "@test isapprox(u, -4.28, atol=1e-2)" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Testing potential energies for all elements" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "function test_energies_for_element(element::String, fs_inter, u; nx::Integer=3, ny::Integer=3, nz::Integer=3, )\n", + " a = bcc_lattice_constants[element]\n", + " atoms, coords, box, box_size, box_vectors = Crystal.make_bcc_unitcell(element, a=a)\n", + " sc_atoms, sc_coords, sc_box, sc_box_size = Crystal.make_supercell(atoms, coords, box, box_size, nx=nx, ny=ny, nz=nz)\n", + " n_atoms = length(sc_atoms)\n", + " \n", + " specific_inter_list = ((fs_inter,),)\n", + " velocities = [velocity(1., .01, dims=3) for i in 1:n_atoms]\n", + " sim = VelocityVerlet()\n", + " nb_matrix = trues(n_atoms,n_atoms)\n", + " n_steps = 1\n", + " dist_cutoff = 2 * a\n", + "\n", + " nf = DistanceNeighbourFinder(nb_matrix, n_steps, dist_cutoff);\n", + "\n", + " loggers = Dict(\"temperature\" => TemperatureLogger(1))\n", + " \n", + " s = Simulation(\n", + " simulator=sim, \n", + " atoms=sc_atoms, \n", + " specific_inter_lists=specific_inter_list,\n", + " general_inters=(),\n", + " coords=[SVector{3}(v) for v in sc_coords], \n", + " velocities=velocities,\n", + " temperature=.01, \n", + " box_size=sc_box_size[1,1],\n", + " timestep=.002,\n", + " n_steps=1,\n", + " neighbour_finder=nf,\n", + " loggers=loggers,\n", + " )\n", + " find_neighbours!(s, s.neighbour_finder, 0)\n", + " u_md = Molly.potential_energy(fs_inter, s)/n_atoms\n", + " return isapprox(u_md, u, atol=1e-2)\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "@testset \"potential energies\" begin \n", + " for element in elements\n", + " element_pair = string(element, element)\n", + " row = reference_energies[fs_inter.element_pair_map[element_pair],:]\n", + " @testset \"$element\" begin \n", + " @test test_energies_for_element(element, fs_inter, -row.u) \n", + " end\n", + " end\n", + "end" + ] + }, { "cell_type": "markdown", "metadata": {}, From 87e7d822cd0b79f711eca7420a5262628ee99c00 Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Sun, 7 Feb 2021 17:37:26 +0100 Subject: [PATCH 05/49] added potential energy and force checks for bcc crystal with a single vacancy --- fs.ipynb | 184 +++++++++++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 178 insertions(+), 6 deletions(-) diff --git a/fs.ipynb b/fs.ipynb index 8591f622..dde8672c 100644 --- a/fs.ipynb +++ b/fs.ipynb @@ -21,8 +21,7 @@ "metadata": {}, "source": [ "TODOs:\n", - "* test forces for crystals with a vacancy\n", - "* test energies for crystals with a vacancy\n", + "* create a new section in `docs/src/examples.md` to document finnis-sinclair type glue potential\n", "* copy tests to the appropriate location in the package\n", "* if moving this nb to docs or so, does `using Molly` from non-root dirs still work?\n", "* other force implementation that does not require altering `function accelerations`?" @@ -194,6 +193,7 @@ "reference_energies = DataFrame(\n", " element_pair = element_pairings,\n", " u = [5.31, 7.57, 8.1, 4.1, 6.82, 8.9, 4.28],\n", + " u_vac = [1.92, 2.64, 3.13, 1.97, 2.58, 3.71, 1.77]\n", ")" ] }, @@ -405,7 +405,8 @@ " n_steps=10,\n", " neighbour_finder=nf,\n", " loggers=loggers,\n", - ")" + ")\n", + "find_neighbours!(s, s.neighbour_finder, 0)" ] }, { @@ -550,8 +551,10 @@ "metadata": {}, "outputs": [], "source": [ - "for element in elements\n", - " @test test_forces_for_element(element, fs_inter)\n", + "@testset \"test forces\" begin\n", + " for element in elements\n", + " @test test_forces_for_element(element, fs_inter)\n", + " end\n", "end" ] }, @@ -856,7 +859,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "scrolled": true + }, "outputs": [], "source": [ "simulate!(s, parallel=false)" @@ -871,6 +876,173 @@ "s.loggers[\"energy\"].energies" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Vacancies" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Ensuring that the vacancy formation energies are within .07 eV of the reference." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "function bcc_vacancy_formation_energy(fs_inter;element::String=\"Fe\",nx::Int64=3,ny::Int64=3,nz::Int64=3,)\n", + " a = bcc_lattice_constants[element]\n", + " atoms, coords, box, box_size, box_vectors = Crystal.make_bcc_unitcell(element, a=a)\n", + " sc_atoms, sc_coords, sc_box, sc_box_size = Crystal.make_supercell(atoms, coords, box, box_size, nx=nx, ny=ny, nz=nz)\n", + " n_atoms = length(sc_atoms)\n", + " \n", + " # energy of the vacancy free system\n", + " specific_inter_list = ((fs_inter,),)\n", + " velocities = [velocity(1., .01, dims=3) for i in 1:n_atoms]\n", + " sim = VelocityVerlet()\n", + " nb_matrix = trues(n_atoms,n_atoms)\n", + " n_steps = 1\n", + " dist_cutoff = 2 * a\n", + "\n", + " nf = DistanceNeighbourFinder(nb_matrix, n_steps, dist_cutoff);\n", + "\n", + " loggers = Dict(\"temperature\" => TemperatureLogger(1))\n", + "\n", + " s = Simulation(\n", + " simulator=sim, \n", + " atoms=sc_atoms, \n", + " specific_inter_lists=specific_inter_list,\n", + " general_inters=(),\n", + " coords=[SVector{3}(v) for v in sc_coords], \n", + " velocities=velocities,\n", + " temperature=.01, \n", + " box_size=sc_box_size[1,1],\n", + " timestep=.002,\n", + " n_steps=1,\n", + " neighbour_finder=nf,\n", + " loggers=loggers,\n", + " )\n", + " find_neighbours!(s, s.neighbour_finder, 0)\n", + " u_md = Molly.potential_energy(fs_inter, s)\n", + " \n", + " # introducing a vacancy\n", + " sc_atoms_vac, sc_coords_vac = Crystal.add_vacancies(sc_atoms, sc_coords, ixs=[1])\n", + " n_atoms_vac = length(sc_atoms_vac)\n", + " \n", + " velocities = [velocity(1., .01, dims=3) for i in 1:n_atoms_vac]\n", + " nb_matrix = trues(n_atoms_vac,n_atoms_vac)\n", + "\n", + " nf = DistanceNeighbourFinder(nb_matrix, n_steps, dist_cutoff);\n", + "\n", + " loggers = Dict(\"temperature\" => TemperatureLogger(1))\n", + "\n", + " s = Simulation(\n", + " simulator=sim, \n", + " atoms=sc_atoms_vac, \n", + " specific_inter_lists=specific_inter_list,\n", + " general_inters=(),\n", + " coords=[SVector{3}(v) for v in sc_coords_vac], \n", + " velocities=velocities,\n", + " temperature=.01, \n", + " box_size=sc_box_size[1,1],\n", + " timestep=.002,\n", + " n_steps=1,\n", + " neighbour_finder=nf,\n", + " loggers=loggers,\n", + " )\n", + " find_neighbours!(s, s.neighbour_finder, 0)\n", + " u_md_vac = Molly.potential_energy(fs_inter, s)\n", + " return u_md_vac - n_atoms_vac/n_atoms * u_md\n", + "end\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "@testset \"potential energies: bcc vacancy\" begin \n", + " for element in elements\n", + " element_pair = string(element, element)\n", + " row = reference_energies[fs_inter.element_pair_map[element_pair],:]\n", + " u_vac = bcc_vacancy_formation_energy(fs_inter, element=element)\n", + " @testset \"$element\" begin \n", + " @test isapprox(u_vac, row.u_vac, atol=7e-2) \n", + " end\n", + " end\n", + "end" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Testing that the forces in an unrelaxed crystal with a vacancy are not all zero." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "function test_bcc_vacancy_forces(fs_inter;element::String=\"Fe\",nx::Int64=3,ny::Int64=3,nz::Int64=3,)\n", + " a = bcc_lattice_constants[element]\n", + " atoms, coords, box, box_size, box_vectors = Crystal.make_bcc_unitcell(element, a=a)\n", + " sc_atoms, sc_coords, sc_box, sc_box_size = Crystal.make_supercell(atoms, coords, box, box_size, nx=nx, ny=ny, nz=nz)\n", + " n_atoms = length(sc_atoms)\n", + " \n", + " # introducing a vacancy\n", + " sc_atoms_vac, sc_coords_vac = Crystal.add_vacancies(sc_atoms, sc_coords, ixs=[1])\n", + " n_atoms_vac = length(sc_atoms_vac)\n", + " \n", + " velocities = [velocity(1., .01, dims=3) for i in 1:n_atoms_vac]\n", + " nb_matrix = trues(n_atoms_vac,n_atoms_vac)\n", + "\n", + " nf = DistanceNeighbourFinder(nb_matrix, n_steps, dist_cutoff);\n", + "\n", + " loggers = Dict(\"temperature\" => TemperatureLogger(1))\n", + "\n", + " s = Simulation(\n", + " simulator=sim, \n", + " atoms=sc_atoms_vac, \n", + " specific_inter_lists=specific_inter_list,\n", + " general_inters=(),\n", + " coords=[SVector{3}(v) for v in sc_coords_vac], \n", + " velocities=velocities,\n", + " temperature=.01, \n", + " box_size=sc_box_size[1,1],\n", + " timestep=.002,\n", + " n_steps=1,\n", + " neighbour_finder=nf,\n", + " loggers=loggers,\n", + " )\n", + " find_neighbours!(s, s.neighbour_finder, 0)\n", + " forces = accelerations(s, parallel=false)\n", + " return !test_forces_zero(forces, n_atoms_vac)\n", + "end\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "@testset \"bcc + vacancy: test forces\" begin\n", + " for element in elements\n", + " @test test_bcc_vacancy_forces(fs_inter, element=element)\n", + " end\n", + "end" + ] + }, { "cell_type": "code", "execution_count": null, From 52f46b8c10512549af9687ca38ff2a214f18349d Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Sun, 14 Feb 2021 13:03:36 +0100 Subject: [PATCH 06/49] added tests for the finnis-sinclair type of glue potentials --- test/glue.jl | 256 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 256 insertions(+) create mode 100644 test/glue.jl diff --git a/test/glue.jl b/test/glue.jl new file mode 100644 index 00000000..06c88ae1 --- /dev/null +++ b/test/glue.jl @@ -0,0 +1,256 @@ +using Molly +using Test +using Crystal +using DataFrames + +elements = ["V", "Nb", "Ta", "Cr", "Mo", "W", "Fe"] +element_pairings = [string(el,el) for el in elements] +element_pair_map = Dict(pair => i for (i,pair) in enumerate(element_pairings)) + +df = DataFrame( + element_pair = element_pairings, + d = [3.692767, 3.915354, 4.076980, 3.915720, 4.114825, 4.400224, 3.699579], + A = [2.010637, 3.013789, 2.591061, 1.453418, 1.887117, 1.896373, 1.889846], + β = [0, 0, 0, 1.8, 0, 0, 1.8], + c = [3.8, 4.2, 4.2, 2.9, 3.25, 3.25, 3.4], + c₀ = [-0.8816318, -1.5640104, 1.2157373, 29.1429813, 43.4475218, 47.1346499, 1.2110601], + c₁ = [1.4907756, 2.0055779, 0.0271471, -23.3975027, -31.9332978, -33.7665655, -0.7510840], + c₂ = [-0.3976370, -0.4663764, -0.1217350, 4.7578297, 6.0804249, 6.2541999, 0.1380773], +) + +masses = Dict("V" => 50.9415, "Nb" => 92.9064, "Ta" => 180.9479, + "Cr" => 51.996, "Mo" => 95.94, "W" => 183.85, + "Fe" => 55.847) + +# Å +bcc_lattice_constants = Dict( + "V" => 3.0399, "Nb" => 3.3008, + "Ta" => 3.3058, "Cr" => 2.8845, "Mo" => 3.1472, + "W" => 3.1652, "Fe" => 2.8665 +) + +reference_energies = DataFrame( + element_pair = element_pairings, + u = [5.31, 7.57, 8.1, 4.1, 6.82, 8.9, 4.28], + u_vac = [1.92, 2.64, 3.13, 1.97, 2.58, 3.71, 1.77] +) + +fs_inter = FinnisSinclair(true, element_pair_map, df) + +function test_forces_zero(forces, n_atoms; dims=3) + zeros = [zero(rand(3)) for _ in 1:n_atoms] + return all(isapprox.(forces, zeros, atol=1e-4)) +end + +function test_forces_for_element(element::String, fs_inter; nx::Integer=3, ny::Integer=3, nz::Integer=3, ) + a = bcc_lattice_constants[element] + atoms, coords, box, box_size, box_vectors = Crystal.make_bcc_unitcell(element, a=a) + sc_atoms, sc_coords, sc_box, sc_box_size = Crystal.make_supercell(atoms, coords, box, box_size, nx=nx, ny=ny, nz=nz) + n_atoms = length(sc_atoms) + + specific_inter_list = ((fs_inter,),) + velocities = [velocity(1., .01, dims=3) for i in 1:n_atoms] + sim = VelocityVerlet() + nb_matrix = trues(n_atoms,n_atoms) + n_steps = 1 + dist_cutoff = 2 * a + + nf = DistanceNeighbourFinder(nb_matrix, n_steps, dist_cutoff) + + loggers = Dict("temperature" => TemperatureLogger(1)) + + s = Simulation( + simulator=sim, + atoms=sc_atoms, + specific_inter_lists=specific_inter_list, + general_inters=(), + coords=[SVector{3}(v) for v in sc_coords], + velocities=velocities, + temperature=.01, + box_size=sc_box_size[1,1], + timestep=.002, + n_steps=1, + neighbour_finder=nf, + loggers=loggers, + ) + find_neighbours!(s, s.neighbour_finder, 0) + forces = accelerations(s, parallel=false) + return test_forces_zero(forces, n_atoms) +end + +function test_bcc_vacancy_forces(fs_inter;element::String="Fe",nx::Int64=3,ny::Int64=3,nz::Int64=3,) + a = bcc_lattice_constants[element] + atoms, coords, box, box_size, box_vectors = Crystal.make_bcc_unitcell(element, a=a) + sc_atoms, sc_coords, sc_box, sc_box_size = Crystal.make_supercell(atoms, coords, box, box_size, nx=nx, ny=ny, nz=nz) + n_atoms = length(sc_atoms) + + # introducing a vacancy + sc_atoms_vac, sc_coords_vac = Crystal.add_vacancies(sc_atoms, sc_coords, ixs=[1]) + n_atoms_vac = length(sc_atoms_vac) + + specific_inter_list = ((fs_inter,),) + velocities = [velocity(1., .01, dims=3) for i in 1:n_atoms_vac] + sim = VelocityVerlet() + nb_matrix = trues(n_atoms_vac,n_atoms_vac) + dist_cutoff = 2 * a + n_steps = 1 + nf = DistanceNeighbourFinder(nb_matrix, n_steps, dist_cutoff) + + loggers = Dict("temperature" => TemperatureLogger(1)) + + s = Simulation( + simulator=sim, + atoms=sc_atoms_vac, + specific_inter_lists=specific_inter_list, + general_inters=(), + coords=[SVector{3}(v) for v in sc_coords_vac], + velocities=velocities, + temperature=.01, + box_size=sc_box_size[1,1], + timestep=.002, + n_steps=1, + neighbour_finder=nf, + loggers=loggers, + ) + find_neighbours!(s, s.neighbour_finder, 0) + forces = accelerations(s, parallel=false) + return !test_forces_zero(forces, n_atoms_vac) +end + +function test_energies_for_element(element::String, fs_inter, u; nx::Integer=3, ny::Integer=3, nz::Integer=3, ) + a = bcc_lattice_constants[element] + atoms, coords, box, box_size, box_vectors = Crystal.make_bcc_unitcell(element, a=a) + sc_atoms, sc_coords, sc_box, sc_box_size = Crystal.make_supercell(atoms, coords, box, box_size, nx=nx, ny=ny, nz=nz) + n_atoms = length(sc_atoms) + + specific_inter_list = ((fs_inter,),) + velocities = [velocity(1., .01, dims=3) for i in 1:n_atoms] + sim = VelocityVerlet() + nb_matrix = trues(n_atoms,n_atoms) + n_steps = 1 + dist_cutoff = 2 * a + + nf = DistanceNeighbourFinder(nb_matrix, n_steps, dist_cutoff); + + loggers = Dict("temperature" => TemperatureLogger(1)) + + s = Simulation( + simulator=sim, + atoms=sc_atoms, + specific_inter_lists=specific_inter_list, + general_inters=(), + coords=[SVector{3}(v) for v in sc_coords], + velocities=velocities, + temperature=.01, + box_size=sc_box_size[1,1], + timestep=.002, + n_steps=1, + neighbour_finder=nf, + loggers=loggers, + ) + find_neighbours!(s, s.neighbour_finder, 0) + u_md = Molly.potential_energy(fs_inter, s)/n_atoms + return isapprox(u_md, u, atol=1e-2) +end + +function bcc_vacancy_formation_energy(fs_inter;element::String="Fe",nx::Int64=3,ny::Int64=3,nz::Int64=3,) + a = bcc_lattice_constants[element] + atoms, coords, box, box_size, box_vectors = Crystal.make_bcc_unitcell(element, a=a) + sc_atoms, sc_coords, sc_box, sc_box_size = Crystal.make_supercell(atoms, coords, box, box_size, nx=nx, ny=ny, nz=nz) + n_atoms = length(sc_atoms) + + # energy of the vacancy free system + specific_inter_list = ((fs_inter,),) + velocities = [velocity(1., .01, dims=3) for i in 1:n_atoms] + sim = VelocityVerlet() + nb_matrix = trues(n_atoms,n_atoms) + n_steps = 1 + dist_cutoff = 2 * a + + nf = DistanceNeighbourFinder(nb_matrix, n_steps, dist_cutoff); + + loggers = Dict("temperature" => TemperatureLogger(1)) + + s = Simulation( + simulator=sim, + atoms=sc_atoms, + specific_inter_lists=specific_inter_list, + general_inters=(), + coords=[SVector{3}(v) for v in sc_coords], + velocities=velocities, + temperature=.01, + box_size=sc_box_size[1,1], + timestep=.002, + n_steps=1, + neighbour_finder=nf, + loggers=loggers, + ) + find_neighbours!(s, s.neighbour_finder, 0) + u_md = Molly.potential_energy(fs_inter, s) + + # introducing a vacancy + sc_atoms_vac, sc_coords_vac = Crystal.add_vacancies(sc_atoms, sc_coords, ixs=[1]) + n_atoms_vac = length(sc_atoms_vac) + + velocities = [velocity(1., .01, dims=3) for i in 1:n_atoms_vac] + nb_matrix = trues(n_atoms_vac,n_atoms_vac) + + nf = DistanceNeighbourFinder(nb_matrix, n_steps, dist_cutoff); + + loggers = Dict("temperature" => TemperatureLogger(1)) + + s = Simulation( + simulator=sim, + atoms=sc_atoms_vac, + specific_inter_lists=specific_inter_list, + general_inters=(), + coords=[SVector{3}(v) for v in sc_coords_vac], + velocities=velocities, + temperature=.01, + box_size=sc_box_size[1,1], + timestep=.002, + n_steps=1, + neighbour_finder=nf, + loggers=loggers, + ) + find_neighbours!(s, s.neighbour_finder, 0) + u_md_vac = Molly.potential_energy(fs_inter, s) + return u_md_vac - n_atoms_vac/n_atoms * u_md +end + + +@testset "Finnis-Sinclair" begin + + @testset "test forces" begin + for element in elements + @test test_forces_for_element(element, fs_inter) + end + end + + @testset "bcc + vacancy: test forces" begin + for element in elements + @test test_bcc_vacancy_forces(fs_inter, element=element) + end + end + + @testset "potential energies" begin + for element in elements + element_pair = string(element, element) + row = reference_energies[fs_inter.element_pair_map[element_pair],:] + @testset "$element" begin + @test test_energies_for_element(element, fs_inter, -row.u) + end + end + end + + @testset "potential energies: bcc vacancy" begin + for element in elements + element_pair = string(element, element) + row = reference_energies[fs_inter.element_pair_map[element_pair],:] + u_vac = bcc_vacancy_formation_energy(fs_inter, element=element) + @testset "$element" begin + @test isapprox(u_vac, row.u_vac, atol=7e-2) + end + end + end +end \ No newline at end of file From dd0881cd85174aa239a9c9623118ec220624ce1f Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Sun, 14 Feb 2021 13:40:54 +0100 Subject: [PATCH 07/49] added section to run a longer simulation for the documentation --- fs.ipynb | 102 ++++++++++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 93 insertions(+), 9 deletions(-) diff --git a/fs.ipynb b/fs.ipynb index dde8672c..99857c7c 100644 --- a/fs.ipynb +++ b/fs.ipynb @@ -11,9 +11,11 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Structure:\n", - "1. Computing Finnis-Sinclair forces\n", - "2. Computing Finnis-Sinclair potential energies" + "Content\n", + "1. Finnis-Sinclair forces for perfect bcc crystals\n", + "2. Finnis-Sinclair potential energies for perfect bcc crystals\n", + "3. Finnis-Sinclair forces for bcc crystals with 1 vacancy\n", + "4. Finnis-Sinclair potential energies for bcc crystals with 1 vacancy" ] }, { @@ -21,10 +23,8 @@ "metadata": {}, "source": [ "TODOs:\n", - "* create a new section in `docs/src/examples.md` to document finnis-sinclair type glue potential\n", - "* copy tests to the appropriate location in the package\n", - "* if moving this nb to docs or so, does `using Molly` from non-root dirs still work?\n", - "* other force implementation that does not require altering `function accelerations`?" + "* check if MD trajectories around room temperature are expected (the crystals should stay relatively stable and not explode - currently it seems they evaporate)\n", + "* Find a more elegant solution for the force implementation that does not require altering `function accelerations`" ] }, { @@ -59,10 +59,12 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "scrolled": true + }, "outputs": [], "source": [ - "# Pkg.add(\"Plots\")" + "# Pkg.add(\"Makie\")" ] }, { @@ -1043,6 +1045,88 @@ "end" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "nx = 3\n", + "ny = 3\n", + "nz = 3\n", + "element = \"Fe\"\n", + "T = 300\n", + "n_steps = 1_000\n", + "dt = .002\n", + "\n", + "a = bcc_lattice_constants[element]\n", + "atoms, coords, box, box_size, box_vectors = Crystal.make_bcc_unitcell(element, a=a)\n", + "sc_atoms, sc_coords, sc_box, sc_box_size = Crystal.make_supercell(atoms, coords, box, box_size, nx=nx, ny=ny, nz=nz)\n", + "n_atoms = length(sc_atoms)\n", + "\n", + "# introducing a vacancy\n", + "sc_atoms_vac, sc_coords_vac = Crystal.add_vacancies(sc_atoms, sc_coords, ixs=[1])\n", + "n_atoms_vac = length(sc_atoms_vac)\n", + "\n", + "specific_inter_list = ((fs_inter,),)\n", + "velocities = [velocity(T, sc_atoms_vac[i].mass, dims=3) for i in 1:n_atoms_vac]\n", + "sim = VelocityVerlet()\n", + "nb_matrix = trues(n_atoms_vac,n_atoms_vac)\n", + "dist_cutoff = 2 * a\n", + "nf = DistanceNeighbourFinder(nb_matrix, n_steps, dist_cutoff)\n", + "\n", + "loggers = Dict(\n", + " \"temperature\" => TemperatureLogger(1),\n", + " \"pot\" => EnergyLogger(1),\n", + " \"coords\" => CoordinateLogger(5),\n", + " \"writer\" => StructureWriter(5,\"traj_fs_bcc_vac_fe.pdb\")\n", + ")\n", + "\n", + "s = Simulation(\n", + " simulator=sim, \n", + " atoms=sc_atoms_vac, \n", + " specific_inter_lists=specific_inter_list,\n", + " general_inters=(),\n", + " coords=[SVector{3}(v) for v in sc_coords_vac], \n", + " velocities=velocities,\n", + " temperature=T, \n", + " box_size=sc_box_size[1,1],\n", + " timestep=dt,\n", + " n_steps=n_steps,\n", + " neighbour_finder=nf,\n", + " loggers=loggers,\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "simulate!(s)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "s.loggers[\"temperature\"].temperatures" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "s.loggers[\"pot\"].energies" + ] + }, { "cell_type": "code", "execution_count": null, From 221b690f94f37739668d62218fa97f8bfe1d2101 Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Wed, 10 Mar 2021 21:05:03 +0100 Subject: [PATCH 08/49] added fields for forces, velocities and glue densities in Simulation --- src/types.jl | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/src/types.jl b/src/types.jl index 149a972c..697cccf1 100644 --- a/src/types.jl +++ b/src/types.jl @@ -184,7 +184,7 @@ default values. - `n_steps_made::Vector{Int}=[]`: the number of steps already made during the simulation. This is a `Vector` to allow the `struct` to be immutable. """ -struct Simulation{D, T, A, C, GI, SI} +struct Simulation{D, T, A, C, GI, SI, GL} simulator::Simulator atoms::A specific_inter_lists::SI @@ -200,9 +200,11 @@ struct Simulation{D, T, A, C, GI, SI} timestep::T n_steps::Int n_steps_made::Vector{Int} + glue_densities::GL + forces::C end -Simulation{D}(args...) where {D, T, A, C, GI, SI} = Simulation{D, T, A, C, GI, SI}(args...) +Simulation{D}(args...) where {D, T, A, C, GI, SI, GL} = Simulation{D, T, A, C, GI, SI, GL}(args...) function Simulation(; simulator, @@ -226,10 +228,14 @@ function Simulation(; C = typeof(coords) GI = typeof(general_inters) SI = typeof(specific_inter_lists) - return Simulation{gpu_diff_safe, T, A, C, GI, SI}( + glue_densities = zeros(length(atoms)) + forces = [zero(v) for v in coords] + GL = typeof(glue_densities) + return Simulation{gpu_diff_safe, T, A, C, GI, SI, GL}( simulator, atoms, specific_inter_lists, general_inters, coords, velocities, temperature, box_size, neighbours, neighbour_finder, - thermostat, loggers, timestep, n_steps, n_steps_made) + thermostat, loggers, timestep, n_steps, n_steps_made, glue_densities, + forces) end function Base.show(io::IO, s::Simulation) From eaacf33d6988e4fc5bda3d66fc36a666dd700408 Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Wed, 10 Mar 2021 21:05:50 +0100 Subject: [PATCH 09/49] refactored glue forces to work as general_inters --- src/interactions/glue_fs.jl | 188 ++++++++++++++++++++++++++++++------ 1 file changed, 157 insertions(+), 31 deletions(-) diff --git a/src/interactions/glue_fs.jl b/src/interactions/glue_fs.jl index 63b631c9..c53864d2 100644 --- a/src/interactions/glue_fs.jl +++ b/src/interactions/glue_fs.jl @@ -17,8 +17,12 @@ struct FinnisSinclair <: GlueInteraction nl_only::Bool element_pair_map::Dict params::DataFrame +# step::Int64 +# glue_densities::Dict{Int64,Float64} # place to store glue densities, one per atom end +FinnisSinclair(nl_only, element_pair_map, params) = FinnisSinclair(nl_only, element_pair_map, params, -1, Dict()) + """ glue_potential(r,β,d) @@ -33,8 +37,8 @@ end Derivative of the pair potential. """ -function pair_potential_derivative(r::T, c::T, c₀::T, c₁::T, c₂::T)::T where T<:Real - return (r > c) ? 0 : 2 * (r - c) * (c₀ + c₁*r + c₂*r^2) + (r - c)^2 * (c₁ + 2*c₂*r) +function pair_potential_derivative(r::T, c::T, c₀::T, c₁::T, c₂::T; f::Float64=10.)::T where T<:Real + return (r > c) ? 0 : 2 * ((r - c)*f) * (c₀ + c₁*r + c₂*r^2) + ((r - c)*f)^2 * (c₁ + 2*c₂*r) end """ @@ -43,7 +47,7 @@ end Energy derivative given glue density. """ function glue_energy_derivative(ρ::Float64, A::Float64)::Float64 - return A/2 * ρ^-1.5 + return - A/(2 * √ρ) end """ @@ -65,17 +69,25 @@ function get_pair_params(element1::String, element2::String, inter::FinnisSincla return inter.params[inter.element_pair_map[pair],:] end -@inline @inbounds function force( +""" + Compute glue densities(inter, coords, s, parallel) + +Convenience function to update the densities before the forces are computed in serial/parallel. +""" +function update_glue_densities!( inter::FinnisSinclair, coords, - s::Simulation + s::Simulation; + parallel::Bool=false ) - # computing the embedding densities n_atoms = length(s.coords) - ρs = zeros(n_atoms) - rs = zeros(length(s.neighbours)) - r_vec_norms = zeros(length(s.neighbours),3) + # wiping old glue densities + for i in 1:n_atoms + s.glue_densities[i] = 0 + end + + # updating glue densities for (n,(i,j)) in enumerate(s.neighbours) element_i = s.atoms[i].name element_j = s.atoms[j].name @@ -87,14 +99,86 @@ end r_vec = vector(s.coords[i], s.coords[j], s.box_size) r2 = sum(abs2, r_vec) r = sqrt(r2) - # storing distance (vectors) so we don't need to recompute - rs[n] = r - r_vec_norms[[n],:] = r_vec / r + # storing glue densities - ρs[i] += glue_potential(r, pj.β, pj.d) - ρs[j] += glue_potential(r, pi.β, pi.d) + s.glue_densities[i] += glue_potential(r, pj.β, pj.d) + s.glue_densities[j] += glue_potential(r, pi.β, pi.d) end + return nothing +end + +@inline @inbounds function force!(forces, inter::FinnisSinclair, s::Simulation, i::Integer, j::Integer) + fdr = force(inter, s, i, j) + forces[i] -= fdr + forces[j] += fdr + return nothing +end + +@inline @inbounds function force(inter::FinnisSinclair, s::Simulation, i::Integer, j::Integer) + element_i = s.atoms[i].name + element_j = s.atoms[j].name + element_pair = string(sort([element_i, element_j])...) + pi = get_pair_params(element_i,element_i, inter) # inter.params[inter.element_map[element_i],:] + pj = get_pair_params(element_j, element_j, inter) # inter.params[inter.element_map[element_j],:] + pij = get_pair_params(element_i, element_j, inter) # inter.params[inter.element_map[element_pair],:] + + r_vec = vector(s.coords[i], s.coords[j], s.box_size) + r2 = sum(abs2, r_vec) + r = sqrt(r2) + dr = r_vec / r + + # pair contribution + f_pair = pair_potential_derivative(r, pij.c, pij.c₀, pij.c₁, pij.c₂) +# dpairdR = f_pair * r_vec_norm +# dpairdR_j = - dpairdR_i + + # glue contribution + dudρ_i = glue_energy_derivative(s.glue_densities[i], pi.A) + dudρ_j = glue_energy_derivative(s.glue_densities[j], pj.A) + dΦdr_i = glue_potential_derivative(r, pi.β, pi.d) + dΦdr_j = glue_potential_derivative(r, pj.β, pj.d) + + f_glue = (dudρ_j * dΦdr_i + dudρ_i * dΦdr_j) +# println("f_glue ", f_glue, " f_pair ", f_pair) + f = f_pair + f_glue + return f * dr +end + + +@inline @inbounds function force_old( + inter::FinnisSinclair, + coords, + s::Simulation + ) + # computing the embedding densities + n_atoms = length(s.coords) +# ρs = zeros(n_atoms) +# rs = zeros(length(s.neighbours)) +# r_vec_norms = zeros(length(s.neighbours),3) +# for (n,(i,j)) in enumerate(s.neighbours) +# element_i = s.atoms[i].name +# element_j = s.atoms[j].name +# element_pair = string(sort([element_i, element_j])...) +# pi = get_pair_params(element_i,element_i,inter) # inter.params[inter.element_map[element_i],:] +# pj = get_pair_params(element_j,element_j,inter) # inter.params[inter.element_map[element_j],:] +# pij = get_pair_params(element_i,element_j,inter) # inter.params[inter.element_map[element_pair],:] + +# r_vec = vector(s.coords[i], s.coords[j], s.box_size) +# r2 = sum(abs2, r_vec) +# r = sqrt(r2) +# # storing distance (vectors) so we don't need to recompute +# rs[n] = r +# r_vec_norms[[n],:] = r_vec / r +# # storing glue densities +# ρs[i] += glue_potential(r, pj.β, pj.d) +# ρs[j] += glue_potential(r, pi.β, pi.d) +# end + println("\n\nwuppety") + println("\nmax r ", maximum(rs), " min r ", minimum(rs), " avg r ", sum(rs)/length(rs), " box ", s.box_size) +# println("\ncoords \n", s.coords) + ρs = inter.glue_densities + println("\nrhos \n", ρs) fs = [zeros(1,3) for _ in 1:n_atoms] for (n,(i,j)) in enumerate(s.neighbours) element_i = s.atoms[i].name @@ -104,13 +188,18 @@ end pj = get_pair_params(element_j,element_j,inter) # inter.params[inter.element_map[element_j],:] pij = get_pair_params(element_i,element_j,inter) # inter.params[inter.element_map[element_pair],:] - r = rs[n] + r_vec = vector(s.coords[i], s.coords[j], s.box_size) + r2 = sum(abs2, r_vec) + r = sqrt(r2) +# r = rs[n] r2 = r^2 - r_vec_norm = r_vec_norms[[n],:] + r_vec_norm = r_vec / r +# r_vec_norm = r_vec_norms[[n],:] # pair contribution - dpairdR_i = r_vec_norm * pair_potential_derivative(r, pij.c, pij.c₀, pij.c₁, pij.c₂) - dpairdR_j = - dpairdR_i + f_pair = - pair_potential_derivative(r, pij.c, pij.c₀, pij.c₁, pij.c₂) +# dpairdR = f_pair * r_vec_norm +# dpairdR_j = - dpairdR_i # glue contribution dudρ_i = glue_energy_derivative(ρs[i], pi.A) @@ -118,20 +207,25 @@ end dΦdr_i = glue_potential_derivative(r, pi.β, pi.d) dΦdr_j = glue_potential_derivative(r, pj.β, pj.d) - ## density change by moving the current atom - dgluedR_i_curr = r_vec_norm * dudρ_i * dΦdr_j - dgluedR_j_curr = r_vec_norm * dudρ_j * dΦdr_i - ## density change by moving a neighbouring atom - dgluedR_i_neigh = - r_vec_norm * dudρ_j * dΦdr_i - dgluedR_j_neigh = - r_vec_norm * dudρ_i * dΦdr_j + f_glue = - (dudρ_j * dΦdr_i + dudρ_i * dΦdr_j) +# println("\nr ", r, " f_pair ", f_pair, " f_glue ", f_glue, " f_i ", f_glue+f_pair) +# dgluedR = f_glue * r_vec_norm +# dgluedR_j = - dgluedR_i +# ## density change by moving the current atom +# dgluedR_i_curr = r_vec_norm * dudρ_i * dΦdr_j +# dgluedR_j_curr = r_vec_norm * dudρ_j * dΦdr_i +# ## density change by moving a neighbouring atom +# dgluedR_i_neigh = - r_vec_norm * dudρ_j * dΦdr_i +# dgluedR_j_neigh = - r_vec_norm * dudρ_i * dΦdr_j # updating the forces - f_i = (dpairdR_i + dgluedR_i_curr + dgluedR_i_neigh) - f_j = (dpairdR_j + dgluedR_j_curr + dgluedR_j_neigh) +# f_i = (dpairdR_i + dgluedR_i_curr + dgluedR_i_neigh) +# f_j = (dpairdR_j + dgluedR_j_curr + dgluedR_j_neigh) + f_i = (f_glue + f_pair) * r_vec_norm #dpairdR + dgluedR fs[i] += f_i - fs[j] += f_j + fs[j] -= f_i end - +# println("\nfs \n", fs) return collect(1:n_atoms), fs end @@ -140,8 +234,8 @@ end Energy contribution directly from atom pair distances. """ -function pair_potential(r::T, c::T, c₀::T, c₁::T, c₂::T)::T where T<:Real - return (r > c) ? 0 : (r - c)^2 * (c₀ + c₁*r + c₂*r^2) +function pair_potential(r::T, c::T, c₀::T, c₁::T, c₂::T; f::Float64=10.)::T where T<:Real + return (r > c) ? 0 : ((r - c)*f)^2 * (c₀ + c₁*r + c₂*r^2) end """ @@ -153,8 +247,40 @@ function glue_energy(ρ::Float64, A::Float64)::Float64 return -A * √ρ end +@inline @inbounds function potential_energy(inter::FinnisSinclair, s::Simulation, i) + # logger - general inters - computes the glue energy part only for a single atom + # note: assumes that the glue densities are up to date, only calcutes them when inter.glue_densities is empty + if length(s.glue_densities) == 0 + for (_i,j) in s.neighbours + element_i = s.atoms[_i].name + element_j = s.atoms[j].name + element_pair = string(sort([element_i, element_j])...) + pi = get_pair_params(element_i,element_i,inter) # inter.params[inter.element_map[element_i],:] + pj = get_pair_params(element_j,element_j,inter) + r_vec = vector(s.coords[_i], s.coords[j], s.box_size) + r2 = sum(abs2, r_vec) + r = sqrt(r2) + s.glue_densities[_i] += glue_potential(r, pj.β, pj.d) + s.glue_densities[j] += glue_potential(r, pi.β, pi.d) + end + end + + element = s.atoms[i].name + A = get_pair_params(element, element, inter).A + return glue_energy(s.glue_densities[i], A) +end + +@inline @inbounds function potential_energy(inter::FinnisSinclair, s::Simulation, i, j) + # logger - general inters - computes the pair energy part only for a single atom pair + r_vec = vector(s.coords[i], s.coords[j], s.box_size) + r2 = sum(abs2, r_vec) + r = sqrt(r2) + pij = get_pair_params(s.atoms[i].name, s.atoms[j].name, inter) + return pair_potential(r, pij.c, pij.c₀, pij.c₁, pij.c₂) +end + @inline @inbounds function potential_energy(inter::FinnisSinclair, s::Simulation) - #computing the potential energy combining glue and pair components. + #logger - specific inters - computing the potential energy combining glue and pair components. e_pair = 0. e_glue = 0. From a9d7ef13cc8aed3fe2b05627f4f363e07c672c7a Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Wed, 10 Mar 2021 21:06:51 +0100 Subject: [PATCH 10/49] modified potential_energy to pre-compute glue densities, if necessary --- src/analysis.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/analysis.jl b/src/analysis.jl index 041bd432..c35d0432 100644 --- a/src/analysis.jl +++ b/src/analysis.jl @@ -87,6 +87,11 @@ function potential_energy(s::Simulation) potential = zero(eltype(eltype(s.coords))) for inter in values(s.general_inters) + if typeof(inter) <: GlueInteraction + for i in 1:n_atoms + potential += potential_energy(inter, s, i) + end + end if inter.nl_only neighbours = s.neighbours @inbounds for ni in 1:length(neighbours) From fd935447e224cda076550fd35df22e8ed3ccd5df Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Wed, 10 Mar 2021 21:08:41 +0100 Subject: [PATCH 11/49] modified accelerations to allow for glue density update and calculation of glue forces using the general_inters loop --- src/forces.jl | 23 ++++++++++++++++++----- 1 file changed, 18 insertions(+), 5 deletions(-) diff --git a/src/forces.jl b/src/forces.jl index 080b5c3f..e7f80349 100644 --- a/src/forces.jl +++ b/src/forces.jl @@ -53,6 +53,11 @@ function accelerations(s::Simulation; parallel::Bool=true) # Loop over interactions and calculate the acceleration due to each for inter in values(s.general_inters) + + if typeof(inter) <: GlueInteraction + update_glue_densities!(inter, s.coords, s) + end + if inter.nl_only neighbours = s.neighbours @threads for ni in 1:length(neighbours) @@ -73,6 +78,11 @@ function accelerations(s::Simulation; parallel::Bool=true) forces = zero(s.coords) for inter in values(s.general_inters) + + if typeof(inter) <: GlueInteraction + update_glue_densities!(inter, s.coords, s) + end + if inter.nl_only neighbours = s.neighbours for ni in 1:length(neighbours) @@ -93,16 +103,19 @@ function accelerations(s::Simulation; parallel::Bool=true) sparse_forces = force.(inter_list, (s.coords,), (s,)) ge1, ge2 = getindex.(sparse_forces, 1), getindex.(sparse_forces, 2) sparse_vec = SparseVector(n_atoms, reduce(vcat, ge1), reduce(vcat, ge2)) - if typeof(inter_list[1]) <: GlueInteraction - forces += Array([SVector{3}(v) for v in sparse_vec]) - else - forces += Array(sparse_vec) - end +# if typeof(inter_list[1]) <: GlueInteraction +# sparse_vec = [SVector{3}(v) for v in sparse_vec] +# end +# println("\nsparse_vec \n", sparse_vec) + forces += Array(sparse_vec) end for i in 1:n_atoms forces[i] /= s.atoms[i].mass + s.forces[i] = forces[i] end + +# println("\nforces \n", forces) return forces end From 95b96358fa4c1dd440b8a258505d9e17228895b1 Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Wed, 10 Mar 2021 21:09:12 +0100 Subject: [PATCH 12/49] added GlueDensityLogger, VelocityLogger and ForcesLogger --- src/loggers.jl | 99 +++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 97 insertions(+), 2 deletions(-) diff --git a/src/loggers.jl b/src/loggers.jl index f418f327..0a44dcfd 100644 --- a/src/loggers.jl +++ b/src/loggers.jl @@ -5,8 +5,10 @@ export log_property!, CoordinateLogger, EnergyLogger, - StructureWriter - + StructureWriter, + GlueDensityLogger, + VelocityLogger, + ForceLogger """ TemperatureLogger(n_steps) TemperatureLogger(T, n_steps) @@ -71,6 +73,99 @@ function log_property!(logger::CoordinateLogger, s::Simulation, step_n::Integer) end end +""" + ForceLogger(n_steps; dims=3) + +Log the forces throughout a simulation. +""" +struct ForceLogger{T} <: Logger + n_steps::Int + forces::Vector{Vector{T}} +end + +function ForceLogger(T, n_steps::Integer; dims::Integer=3) + return ForceLogger(n_steps, + Array{SArray{Tuple{dims}, T, 1, dims}, 1}[]) +end + +function ForceLogger(n_steps::Integer; dims::Integer=3) + return ForceLogger(DefaultFloat, n_steps, dims=dims) +end + +function Base.show(io::IO, cl::ForceLogger) + print(io, "ForceLogger{", eltype(eltype(cl.forces)), "} with n_steps ", + cl.n_steps, ", ", length(cl.forces), " frames recorded for ", + length(cl.forces) > 0 ? length(first(cl.forces)) : "?", " atoms") +end + +function log_property!(logger::ForceLogger, s::Simulation, step_n::Integer) + if step_n % logger.n_steps == 0 + push!(logger.forces, deepcopy(s.forces)) + end +end + +""" + VelocityLogger(n_steps; dims=3) + +Log the velocities throughout a simulation. +""" +struct VelocityLogger{T} <: Logger + n_steps::Int + velocities::Vector{Vector{T}} +end + +function VelocityLogger(T, n_steps::Integer; dims::Integer=3) + return VelocityLogger(n_steps, + Array{SArray{Tuple{dims}, T, 1, dims}, 1}[]) +end + +function VelocityLogger(n_steps::Integer; dims::Integer=3) + return VelocityLogger(DefaultFloat, n_steps, dims=dims) +end + +function Base.show(io::IO, cl::VelocityLogger) + print(io, "VelocityLogger{", eltype(eltype(cl.velocities)), "} with n_steps ", + cl.n_steps, ", ", length(cl.velocities), " frames recorded for ", + length(cl.velocities) > 0 ? length(first(cl.velocities)) : "?", " atoms") +end + +function log_property!(logger::VelocityLogger, s::Simulation, step_n::Integer) + if step_n % logger.n_steps == 0 + push!(logger.velocities, deepcopy(s.velocities)) + end +end + +""" + GlueDensityLogger(n_steps) + +Log the glue densities throughout a simulation. +""" +struct GlueDensityLogger{T} <: Logger + n_steps::Int + glue_densities::Vector{Vector{T}} +end + +function GlueDensityLogger(T, n_steps::Integer) + return GlueDensityLogger(n_steps, + Array{T, 1}[]) +end + +function GlueDensityLogger(n_steps::Integer) + return GlueDensityLogger(DefaultFloat, n_steps) +end + +function Base.show(io::IO, cl::CoordinateLogger) + print(io, "GlueDensityLogger{", eltype(eltype(cl.glue_densities)), "} with n_steps ", + cl.n_steps, ", ", length(cl.glue_densities), " frames recorded for ", + length(cl.glue_densities) > 0 ? length(first(cl.glue_densities)) : "?", " atoms") +end + +function log_property!(logger::GlueDensityLogger, s::Simulation, step_n::Integer) + if step_n % logger.n_steps == 0 + push!(logger.glue_densities, deepcopy(s.glue_densities)) + end +end + """ EnergyLogger(n_steps) From 36457e0a56b9b0ff6210dd7e178b315f6fb9075a Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Wed, 10 Mar 2021 21:21:33 +0100 Subject: [PATCH 13/49] fixed copy-paste error --- src/loggers.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/loggers.jl b/src/loggers.jl index 0a44dcfd..17ce7599 100644 --- a/src/loggers.jl +++ b/src/loggers.jl @@ -154,7 +154,7 @@ function GlueDensityLogger(n_steps::Integer) return GlueDensityLogger(DefaultFloat, n_steps) end -function Base.show(io::IO, cl::CoordinateLogger) +function Base.show(io::IO, cl::GlueDensityLogger) print(io, "GlueDensityLogger{", eltype(eltype(cl.glue_densities)), "} with n_steps ", cl.n_steps, ", ", length(cl.glue_densities), " frames recorded for ", length(cl.glue_densities) > 0 ? length(first(cl.glue_densities)) : "?", " atoms") From 180848900896d1ff385fcbb2d7231202d6821920 Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Sat, 13 Mar 2021 00:15:16 +0100 Subject: [PATCH 14/49] added ForwardDiff --- Project.toml | 2 ++ src/Molly.jl | 1 + 2 files changed, 3 insertions(+) diff --git a/Project.toml b/Project.toml index 6eff02af..d8c00e4d 100644 --- a/Project.toml +++ b/Project.toml @@ -10,8 +10,10 @@ Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" Crystal = "3c6eccdf-2a89-4c24-a1d4-ff210daa8476" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" KernelDensity = "5ab0869b-81aa-558d-bb23-cbf5423bbe9b" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" diff --git a/src/Molly.jl b/src/Molly.jl index 25d8851c..b1fa536e 100644 --- a/src/Molly.jl +++ b/src/Molly.jl @@ -14,6 +14,7 @@ using Requires using Base.Threads using LinearAlgebra using SparseArrays +using ForwardDiff include("types.jl") include("cutoffs.jl") From 057b9a6c45bac3d3baa94593c65edfd384fc66ab Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Sat, 13 Mar 2021 00:17:18 +0100 Subject: [PATCH 15/49] refactored reducing redundancy and replacing manually implemented derivates with FordwardDiff.derivate to re-use the non-derivative forms also used to successfully test ground state energies --- src/interactions/glue_fs.jl | 229 ++++++++++++++++++++---------------- 1 file changed, 126 insertions(+), 103 deletions(-) diff --git a/src/interactions/glue_fs.jl b/src/interactions/glue_fs.jl index c53864d2..b4a5a402 100644 --- a/src/interactions/glue_fs.jl +++ b/src/interactions/glue_fs.jl @@ -17,48 +17,72 @@ struct FinnisSinclair <: GlueInteraction nl_only::Bool element_pair_map::Dict params::DataFrame -# step::Int64 -# glue_densities::Dict{Int64,Float64} # place to store glue densities, one per atom + f::Real end -FinnisSinclair(nl_only, element_pair_map, params) = FinnisSinclair(nl_only, element_pair_map, params, -1, Dict()) - """ - glue_potential(r,β,d) + glue(r,β,d; f=10.) -The core component of the Finnis-Sinclair type GlueInteraction. Used to calculate contribution to the glue value of an atom, based on its neighbouring atoms. +The core component of the Finnis-Sinclair type GlueInteraction. +Used to calculate contribution to the glue value of an atom, based on +its neighbouring atoms. +f is a fudge factor to help translate a Å model to a nm model. """ -function glue_potential(r::T, β::T, d::T)::T where T<:Real +function glue(r, β, d) return r > d ? 0 : (r-d)^2 + β*(r-d)^3/d end """ - pair_potential_derivative(r, c, c₀, c₁, c₂) + ∂glue_∂r(r, β, d; f=10.) + +Derivative of the glue density function. +""" +∂glue_∂r(r, β, d) = ForwardDiff.derivative(r -> glue(r,β,d), r) +# function glue_potential_derivative(r::T, β::T, d::T)::T where T<:Real +# return r > d ? 0 : 2*(r-d) + 3*β*(r-d)^2/d +# end +# ∂glue_∂r(r,β,d) = ForwardDiff.derivative(r -> glue_potential(r,β,d), r) -Derivative of the pair potential. """ -function pair_potential_derivative(r::T, c::T, c₀::T, c₁::T, c₂::T; f::Float64=10.)::T where T<:Real - return (r > c) ? 0 : 2 * ((r - c)*f) * (c₀ + c₁*r + c₂*r^2) + ((r - c)*f)^2 * (c₁ + 2*c₂*r) + Uglue(ρ, A) + +Energy based on the glue density . +""" +function Uglue(ρ, A) + return -A * √ρ end + """ - glue_energy_derivative(ρ, A) + ∂Uglue_∂ρ(ρ, A) Energy derivative given glue density. """ -function glue_energy_derivative(ρ::Float64, A::Float64)::Float64 - return - A/(2 * √ρ) -end +∂Uglue_∂ρ(ρ,A) = ForwardDiff.derivative(ρ -> Uglue(ρ,A), ρ) +# function glue_energy_derivative(ρ::Float64, A::Float64)::Float64 +# return - A/(2 * √ρ) +# end """ - glue_potential_derivative(r, β, d) + Upair(r, c, c₀, c₁, c₂; f=10.) -Derivative of the glue density function. +Energy contribution directly from atom pair distances. f is a fudge factor to help translate a Å model to a nm model. """ -function glue_potential_derivative(r::T, β::T, d::T)::T where T<:Real - return r > d ? 0 : 2*(r-d) + 3*β*(r-d)^2/d +function Upair(r, c, c₀, c₁, c₂) + return (r > c) ? 0 : (r - c)^2 * (c₀ + c₁*r + c₂*r^2) end +""" + ∂Upair_∂r(r, c, c₀, c₁, c₂, f=10.) + +Derivative of the pair potential. +f is a fudge factor to help translate a Å model to a nm model. +""" +∂Upair_∂r(r, c, c₀,c₁, c₂) = ForwardDiff.derivative(r -> Upair(r,c,c₀,c₁,c₂), r) +# function pair_potential_derivative(r::T, c::T, c₀::T, c₁::T, c₂::T; f::Float64=10.)::T where T<:Real +# return (r > c) ? 0 : 2 * ((r - c)*f) * (c₀ + c₁*r + c₂*r^2) + ((r - c)*f)^2 * (c₁ + 2*c₂*r) +# end + """ get_pair_params(element1, element2, inter) @@ -70,7 +94,7 @@ function get_pair_params(element1::String, element2::String, inter::FinnisSincla end """ - Compute glue densities(inter, coords, s, parallel) + update_glue_densities!(inter, coords, s, parallel) Convenience function to update the densities before the forces are computed in serial/parallel. """ @@ -89,20 +113,22 @@ function update_glue_densities!( # updating glue densities for (n,(i,j)) in enumerate(s.neighbours) - element_i = s.atoms[i].name - element_j = s.atoms[j].name - element_pair = string(sort([element_i, element_j])...) - pi = get_pair_params(element_i,element_i,inter) # inter.params[inter.element_map[element_i],:] - pj = get_pair_params(element_j,element_j,inter) # inter.params[inter.element_map[element_j],:] - pij = get_pair_params(element_i,element_j,inter) # inter.params[inter.element_map[element_pair],:] + + # collecting parameters + el_i = s.atoms[i].name + el_j = s.atoms[j].name + pi = get_pair_params(el_i,el_i,inter) + pj = get_pair_params(el_j,el_j,inter) + pij = get_pair_params(el_i,el_j,inter) - r_vec = vector(s.coords[i], s.coords[j], s.box_size) + # computing distance + r_vec = vector(s.coords[i], s.coords[j], s.box_size) .* inter.f r2 = sum(abs2, r_vec) r = sqrt(r2) - - # storing glue densities - s.glue_densities[i] += glue_potential(r, pj.β, pj.d) - s.glue_densities[j] += glue_potential(r, pi.β, pi.d) + + # updating glue densities + s.glue_densities[i] += glue(r, pj.β, pj.d) + s.glue_densities[j] += glue(r, pi.β, pi.d) end return nothing end @@ -115,35 +141,82 @@ end end @inline @inbounds function force(inter::FinnisSinclair, s::Simulation, i::Integer, j::Integer) + + # parameters element_i = s.atoms[i].name element_j = s.atoms[j].name element_pair = string(sort([element_i, element_j])...) - pi = get_pair_params(element_i,element_i, inter) # inter.params[inter.element_map[element_i],:] - pj = get_pair_params(element_j, element_j, inter) # inter.params[inter.element_map[element_j],:] - pij = get_pair_params(element_i, element_j, inter) # inter.params[inter.element_map[element_pair],:] + pi = get_pair_params(element_i, element_i, inter) + pj = get_pair_params(element_j, element_j, inter) + pij = get_pair_params(element_i, element_j, inter) - r_vec = vector(s.coords[i], s.coords[j], s.box_size) + # distance i -> j + r_vec = vector(s.coords[i], s.coords[j], s.box_size) .* inter.f r2 = sum(abs2, r_vec) r = sqrt(r2) dr = r_vec / r # pair contribution - f_pair = pair_potential_derivative(r, pij.c, pij.c₀, pij.c₁, pij.c₂) -# dpairdR = f_pair * r_vec_norm -# dpairdR_j = - dpairdR_i + f_pair = ∂Upair_∂r(r, pij.c, pij.c₀, pij.c₁, pij.c₂) # glue contribution - dudρ_i = glue_energy_derivative(s.glue_densities[i], pi.A) - dudρ_j = glue_energy_derivative(s.glue_densities[j], pj.A) - dΦdr_i = glue_potential_derivative(r, pi.β, pi.d) - dΦdr_j = glue_potential_derivative(r, pj.β, pj.d) - + dudρ_i = ∂Uglue_∂ρ(s.glue_densities[i], pi.A) + dudρ_j = ∂Uglue_∂ρ(s.glue_densities[j], pj.A) + dΦdr_i = ∂glue_∂r(r, pi.β, pi.d) + dΦdr_j = ∂glue_∂r(r, pj.β, pj.d) f_glue = (dudρ_j * dΦdr_i + dudρ_i * dΦdr_j) -# println("f_glue ", f_glue, " f_pair ", f_pair) + + # total force contribution f = f_pair + f_glue return f * dr end +@inline @inbounds function potential_energy(inter::FinnisSinclair, s::Simulation, i) + # logger - general inters - computes the glue energy part only for a single atom + # note: assumes that the glue densities are up to date, only calcutes them when inter.glue_densities is empty + + # check if densities are zero, if so calculate current, otherwise assume they are current + no_glue = all(isapprox.(s.glue_densities, zeros(length(s.glue_densities)), atol=1e-4)) + if no_glue + update_glue_densities!(inter, s.coords, s, parallel=false) + # for (_i,_j) in s.neighbours + + # # parameters + # el_i = s.atoms[_i].name + # el_j = s.atoms[_j].name + # pi = get_pair_params(el_i,el_i,inter) + # pj = get_pair_params(el_j,el_j,inter) + + # # distance + # r_vec = vector(s.coords[_i], s.coords[_j], s.box_size) .* inter.f + # r2 = sum(abs2, r_vec) + # r = sqrt(r2) + + # # glue update + # s.glue_densities[_i] += glue(r, pj.β, pj.d) + # s.glue_densities[_j] += glue(r, pi.β, pi.d) + # end + end + + A = get_pair_params(s.atoms[i].name, s.atoms[i].name, inter).A + # u = Uglue(s.glue_densities[i], A) + # println("densities ", s.glue_densities) + # println("Uglue ", u) + # return u + return Uglue(s.glue_densities[i], A) +end + +@inline @inbounds function potential_energy(inter::FinnisSinclair, s::Simulation, i, j) + # logger - general inters - computes the pair energy part only for a single atom pair + r_vec = vector(s.coords[i], s.coords[j], s.box_size) .* inter.f + r2 = sum(abs2, r_vec) + r = sqrt(r2) + pij = get_pair_params(s.atoms[i].name, s.atoms[j].name, inter) + # u = Upair(r, pij.c, pij.c₀, pij.c₁, pij.c₂) + # println("Upair ", u) + # return u + return Upair(r, pij.c, pij.c₀, pij.c₁, pij.c₂) +end @inline @inbounds function force_old( inter::FinnisSinclair, @@ -197,15 +270,15 @@ end # r_vec_norm = r_vec_norms[[n],:] # pair contribution - f_pair = - pair_potential_derivative(r, pij.c, pij.c₀, pij.c₁, pij.c₂) + f_pair = - ∂Upair_∂r(r, pij.c, pij.c₀, pij.c₁, pij.c₂) # dpairdR = f_pair * r_vec_norm # dpairdR_j = - dpairdR_i # glue contribution - dudρ_i = glue_energy_derivative(ρs[i], pi.A) - dudρ_j = glue_energy_derivative(ρs[j], pj.A) - dΦdr_i = glue_potential_derivative(r, pi.β, pi.d) - dΦdr_j = glue_potential_derivative(r, pj.β, pj.d) + dudρ_i = ∂Uglue_∂ρ(ρs[i], pi.A) + dudρ_j = ∂Uglue_∂ρ(ρs[j], pj.A) + dΦdr_i = ∂glue_∂r(r, pi.β, pi.d) + dΦdr_j = ∂glue_∂r(r, pj.β, pj.d) f_glue = - (dudρ_j * dΦdr_i + dudρ_i * dΦdr_j) # println("\nr ", r, " f_pair ", f_pair, " f_glue ", f_glue, " f_i ", f_glue+f_pair) @@ -229,56 +302,6 @@ end return collect(1:n_atoms), fs end -""" - pair_potential(r, c, c₀, c₁, c₂) - -Energy contribution directly from atom pair distances. -""" -function pair_potential(r::T, c::T, c₀::T, c₁::T, c₂::T; f::Float64=10.)::T where T<:Real - return (r > c) ? 0 : ((r - c)*f)^2 * (c₀ + c₁*r + c₂*r^2) -end - -""" - glue_energy(ρ, A) - -Energy based on the glue density . -""" -function glue_energy(ρ::Float64, A::Float64)::Float64 - return -A * √ρ -end - -@inline @inbounds function potential_energy(inter::FinnisSinclair, s::Simulation, i) - # logger - general inters - computes the glue energy part only for a single atom - # note: assumes that the glue densities are up to date, only calcutes them when inter.glue_densities is empty - if length(s.glue_densities) == 0 - for (_i,j) in s.neighbours - element_i = s.atoms[_i].name - element_j = s.atoms[j].name - element_pair = string(sort([element_i, element_j])...) - pi = get_pair_params(element_i,element_i,inter) # inter.params[inter.element_map[element_i],:] - pj = get_pair_params(element_j,element_j,inter) - r_vec = vector(s.coords[_i], s.coords[j], s.box_size) - r2 = sum(abs2, r_vec) - r = sqrt(r2) - s.glue_densities[_i] += glue_potential(r, pj.β, pj.d) - s.glue_densities[j] += glue_potential(r, pi.β, pi.d) - end - end - - element = s.atoms[i].name - A = get_pair_params(element, element, inter).A - return glue_energy(s.glue_densities[i], A) -end - -@inline @inbounds function potential_energy(inter::FinnisSinclair, s::Simulation, i, j) - # logger - general inters - computes the pair energy part only for a single atom pair - r_vec = vector(s.coords[i], s.coords[j], s.box_size) - r2 = sum(abs2, r_vec) - r = sqrt(r2) - pij = get_pair_params(s.atoms[i].name, s.atoms[j].name, inter) - return pair_potential(r, pij.c, pij.c₀, pij.c₁, pij.c₂) -end - @inline @inbounds function potential_energy(inter::FinnisSinclair, s::Simulation) #logger - specific inters - computing the potential energy combining glue and pair components. @@ -298,16 +321,16 @@ end r2 = sum(abs2, r_vec) r = sqrt(r2) - e_pair += pair_potential(r, pij.c, pij.c₀, pij.c₁, pij.c₂) + e_pair += Upair(r, pij.c, pij.c₀, pij.c₁, pij.c₂) - ρs[i] += glue_potential(r, pj.β, pj.d) - ρs[j] += glue_potential(r, pi.β, pi.d) + ρs[i] += glue(r, pj.β, pj.d) + ρs[j] += glue(r, pi.β, pi.d) end es_glue = zeros(n_atoms) for (i, atom) in enumerate(s.atoms) A = get_pair_params(atom.name, atom.name, inter).A - es_glue[i] = glue_energy(ρs[i], A) + es_glue[i] = Uglue(ρs[i], A) end e_glue = sum(es_glue) return e_pair + e_glue From acf7ce856af8c6fd724f3ebca35c8063b761d35e Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Mon, 15 Mar 2021 21:06:22 +0100 Subject: [PATCH 16/49] updated algebraic signs --- src/interactions/glue_fs.jl | 30 ++++++++++++++++-------------- 1 file changed, 16 insertions(+), 14 deletions(-) diff --git a/src/interactions/glue_fs.jl b/src/interactions/glue_fs.jl index b4a5a402..4b82ea76 100644 --- a/src/interactions/glue_fs.jl +++ b/src/interactions/glue_fs.jl @@ -122,9 +122,8 @@ function update_glue_densities!( pij = get_pair_params(el_i,el_j,inter) # computing distance - r_vec = vector(s.coords[i], s.coords[j], s.box_size) .* inter.f - r2 = sum(abs2, r_vec) - r = sqrt(r2) + dr = vector(s.coords[i], s.coords[j], s.box_size) .* inter.f + r = norm(dr) # updating glue densities s.glue_densities[i] += glue(r, pj.β, pj.d) @@ -134,7 +133,7 @@ function update_glue_densities!( end @inline @inbounds function force!(forces, inter::FinnisSinclair, s::Simulation, i::Integer, j::Integer) - fdr = force(inter, s, i, j) + fdr = force(inter, s, i, j) .* inter.f forces[i] -= fdr forces[j] += fdr return nothing @@ -151,20 +150,19 @@ end pij = get_pair_params(element_i, element_j, inter) # distance i -> j - r_vec = vector(s.coords[i], s.coords[j], s.box_size) .* inter.f - r2 = sum(abs2, r_vec) - r = sqrt(r2) - dr = r_vec / r + dr = vector(s.coords[i], s.coords[j], s.box_size) .* inter.f + r = norm(dr) + dr = normalize(dr) # pair contribution - f_pair = ∂Upair_∂r(r, pij.c, pij.c₀, pij.c₁, pij.c₂) + f_pair = - ∂Upair_∂r(r, pij.c, pij.c₀, pij.c₁, pij.c₂) # glue contribution dudρ_i = ∂Uglue_∂ρ(s.glue_densities[i], pi.A) dudρ_j = ∂Uglue_∂ρ(s.glue_densities[j], pj.A) dΦdr_i = ∂glue_∂r(r, pi.β, pi.d) dΦdr_j = ∂glue_∂r(r, pj.β, pj.d) - f_glue = (dudρ_j * dΦdr_i + dudρ_i * dΦdr_j) + f_glue = - (dudρ_j * dΦdr_i + dudρ_i * dΦdr_j) # total force contribution f = f_pair + f_glue @@ -208,9 +206,10 @@ end @inline @inbounds function potential_energy(inter::FinnisSinclair, s::Simulation, i, j) # logger - general inters - computes the pair energy part only for a single atom pair - r_vec = vector(s.coords[i], s.coords[j], s.box_size) .* inter.f - r2 = sum(abs2, r_vec) - r = sqrt(r2) + dr = vector(s.coords[i], s.coords[j], s.box_size) .* inter.f + r = norm(dr) + # r2 = sum(abs2, r_vec) + # r = sqrt(r2) pij = get_pair_params(s.atoms[i].name, s.atoms[j].name, inter) # u = Upair(r, pij.c, pij.c₀, pij.c₁, pij.c₂) # println("Upair ", u) @@ -218,6 +217,7 @@ end return Upair(r, pij.c, pij.c₀, pij.c₁, pij.c₂) end +#= @inline @inbounds function force_old( inter::FinnisSinclair, coords, @@ -334,4 +334,6 @@ end end e_glue = sum(es_glue) return e_pair + e_glue -end \ No newline at end of file +end + +=# \ No newline at end of file From b99aca387ac53b3538b9fd56b828fe9bab0a9522 Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Tue, 16 Mar 2021 20:22:19 +0100 Subject: [PATCH 17/49] added convenience function to easily generate variables for the Finnis-Sinclair 1984 empirical model --- src/interactions/glue_fs.jl | 230 +++++++++++------------------------- 1 file changed, 66 insertions(+), 164 deletions(-) diff --git a/src/interactions/glue_fs.jl b/src/interactions/glue_fs.jl index 4b82ea76..b6986910 100644 --- a/src/interactions/glue_fs.jl +++ b/src/interactions/glue_fs.jl @@ -11,13 +11,72 @@ abstract type GlueInteraction <: SpecificInteraction end """ FinnisSinclairInteraction(nl_only,element_pair_map,params) -The Finnis-Sinclair interaction. +The Finnis-Sinclair interaction. This interaction expects units to be of +these https://lammps.sandia.gov/doc/units.html units (eV, Å, K, ps and so on). """ struct FinnisSinclair <: GlueInteraction nl_only::Bool element_pair_map::Dict params::DataFrame - f::Real +end + + +kb = 8.617333262145e-5 # eV/K + +""" + get_finnissinclair1984(nl_only) + +Finnis and Sinclair 1984 parameterization: https://doi.org/10.1080/01418618408244210 +""" +function get_finnissinclair1984(nl_only::Bool) + + elements = ["V", "Nb", "Ta", "Cr", "Mo", "W", "Fe"] + element_pairings = [string(el, el) for el in elements] + element_pair_map = Dict(pair => i for (i,pair) in enumerate(element_pairings)) + + df = DataFrame( + element_pair = element_pairings, + d = [3.692767, 3.915354, 4.076980, 3.915720, + 4.114825, 4.400224, 3.699579], + # (Å) + A = [2.010637, 3.013789, 2.591061, 1.453418, + 1.887117, 1.896373, 1.889846], + # (eV) + β = [0, 0, 0, 1.8, 0, 0, 1.8], + # (1) + c = [3.8, 4.2, 4.2, 2.9, 3.25, 3.25, 3.4], + # (Å) + c₀ = [-0.8816318, -1.5640104, 1.2157373, 29.1429813, + 43.4475218, 47.1346499, 1.2110601], + # (1) + c₁ = [1.4907756, 2.0055779, 0.0271471, -23.3975027, + -31.9332978, -33.7665655, -0.7510840], + # (1) + c₂ = [-0.3976370, -0.4663764, -0.1217350, 4.7578297, + 6.0804249, 6.2541999, 0.1380773], + # (1) + ) + + masses = Dict("V" => 50.9415, "Nb" => 92.9064, "Ta" => 180.9479, + "Cr" => 51.996, "Mo" => 95.94, "W" => 183.85, + "Fe" => 55.847) # g/mole + + # Å + bcc_lattice_constants = Dict( + "V" => 3.0399, "Nb" => 3.3008, "Ta" => 3.3058, + "Cr" => 2.8845, "Mo" => 3.1472, "W" => 3.1652, + "Fe" => 2.8665 + ) + + + fs84 = FinnisSinclair(nl_only, element_pair_map, df) + + reference_energies = DataFrame( + element_pair = element_pairings, + u = [5.31, 7.57, 8.1, 4.1, 6.82, 8.9, 4.28], + u_vac = [1.92, 2.64, 3.13, 1.97, 2.58, 3.71, 1.77] + ) + return fs84, elements, masses, bcc_lattice_constants, reference_energies end """ @@ -38,10 +97,6 @@ end Derivative of the glue density function. """ ∂glue_∂r(r, β, d) = ForwardDiff.derivative(r -> glue(r,β,d), r) -# function glue_potential_derivative(r::T, β::T, d::T)::T where T<:Real -# return r > d ? 0 : 2*(r-d) + 3*β*(r-d)^2/d -# end -# ∂glue_∂r(r,β,d) = ForwardDiff.derivative(r -> glue_potential(r,β,d), r) """ Uglue(ρ, A) @@ -59,9 +114,6 @@ end Energy derivative given glue density. """ ∂Uglue_∂ρ(ρ,A) = ForwardDiff.derivative(ρ -> Uglue(ρ,A), ρ) -# function glue_energy_derivative(ρ::Float64, A::Float64)::Float64 -# return - A/(2 * √ρ) -# end """ Upair(r, c, c₀, c₁, c₂; f=10.) @@ -79,9 +131,6 @@ Derivative of the pair potential. f is a fudge factor to help translate a Å model to a nm model. """ ∂Upair_∂r(r, c, c₀,c₁, c₂) = ForwardDiff.derivative(r -> Upair(r,c,c₀,c₁,c₂), r) -# function pair_potential_derivative(r::T, c::T, c₀::T, c₁::T, c₂::T; f::Float64=10.)::T where T<:Real -# return (r > c) ? 0 : 2 * ((r - c)*f) * (c₀ + c₁*r + c₂*r^2) + ((r - c)*f)^2 * (c₁ + 2*c₂*r) -# end """ get_pair_params(element1, element2, inter) @@ -122,7 +171,7 @@ function update_glue_densities!( pij = get_pair_params(el_i,el_j,inter) # computing distance - dr = vector(s.coords[i], s.coords[j], s.box_size) .* inter.f + dr = vector(s.coords[i], s.coords[j], s.box_size) r = norm(dr) # updating glue densities @@ -133,7 +182,7 @@ function update_glue_densities!( end @inline @inbounds function force!(forces, inter::FinnisSinclair, s::Simulation, i::Integer, j::Integer) - fdr = force(inter, s, i, j) .* inter.f + fdr = force(inter, s, i, j) forces[i] -= fdr forces[j] += fdr return nothing @@ -150,7 +199,7 @@ end pij = get_pair_params(element_i, element_j, inter) # distance i -> j - dr = vector(s.coords[i], s.coords[j], s.box_size) .* inter.f + dr = vector(s.coords[i], s.coords[j], s.box_size) r = norm(dr) dr = normalize(dr) @@ -171,169 +220,22 @@ end @inline @inbounds function potential_energy(inter::FinnisSinclair, s::Simulation, i) # logger - general inters - computes the glue energy part only for a single atom - # note: assumes that the glue densities are up to date, only calcutes them when inter.glue_densities is empty + # note: assumes that the glue densities are up to date, only calcutes them when inter.glue_densities are all 0 # check if densities are zero, if so calculate current, otherwise assume they are current no_glue = all(isapprox.(s.glue_densities, zeros(length(s.glue_densities)), atol=1e-4)) if no_glue update_glue_densities!(inter, s.coords, s, parallel=false) - # for (_i,_j) in s.neighbours - - # # parameters - # el_i = s.atoms[_i].name - # el_j = s.atoms[_j].name - # pi = get_pair_params(el_i,el_i,inter) - # pj = get_pair_params(el_j,el_j,inter) - - # # distance - # r_vec = vector(s.coords[_i], s.coords[_j], s.box_size) .* inter.f - # r2 = sum(abs2, r_vec) - # r = sqrt(r2) - - # # glue update - # s.glue_densities[_i] += glue(r, pj.β, pj.d) - # s.glue_densities[_j] += glue(r, pi.β, pi.d) - # end end A = get_pair_params(s.atoms[i].name, s.atoms[i].name, inter).A - # u = Uglue(s.glue_densities[i], A) - # println("densities ", s.glue_densities) - # println("Uglue ", u) - # return u return Uglue(s.glue_densities[i], A) end @inline @inbounds function potential_energy(inter::FinnisSinclair, s::Simulation, i, j) # logger - general inters - computes the pair energy part only for a single atom pair - dr = vector(s.coords[i], s.coords[j], s.box_size) .* inter.f + dr = vector(s.coords[i], s.coords[j], s.box_size) r = norm(dr) - # r2 = sum(abs2, r_vec) - # r = sqrt(r2) pij = get_pair_params(s.atoms[i].name, s.atoms[j].name, inter) - # u = Upair(r, pij.c, pij.c₀, pij.c₁, pij.c₂) - # println("Upair ", u) - # return u return Upair(r, pij.c, pij.c₀, pij.c₁, pij.c₂) end - -#= -@inline @inbounds function force_old( - inter::FinnisSinclair, - coords, - s::Simulation - ) - # computing the embedding densities - n_atoms = length(s.coords) -# ρs = zeros(n_atoms) -# rs = zeros(length(s.neighbours)) -# r_vec_norms = zeros(length(s.neighbours),3) - -# for (n,(i,j)) in enumerate(s.neighbours) -# element_i = s.atoms[i].name -# element_j = s.atoms[j].name -# element_pair = string(sort([element_i, element_j])...) -# pi = get_pair_params(element_i,element_i,inter) # inter.params[inter.element_map[element_i],:] -# pj = get_pair_params(element_j,element_j,inter) # inter.params[inter.element_map[element_j],:] -# pij = get_pair_params(element_i,element_j,inter) # inter.params[inter.element_map[element_pair],:] - -# r_vec = vector(s.coords[i], s.coords[j], s.box_size) -# r2 = sum(abs2, r_vec) -# r = sqrt(r2) -# # storing distance (vectors) so we don't need to recompute -# rs[n] = r -# r_vec_norms[[n],:] = r_vec / r -# # storing glue densities -# ρs[i] += glue_potential(r, pj.β, pj.d) -# ρs[j] += glue_potential(r, pi.β, pi.d) -# end - println("\n\nwuppety") - println("\nmax r ", maximum(rs), " min r ", minimum(rs), " avg r ", sum(rs)/length(rs), " box ", s.box_size) -# println("\ncoords \n", s.coords) - ρs = inter.glue_densities - println("\nrhos \n", ρs) - fs = [zeros(1,3) for _ in 1:n_atoms] - for (n,(i,j)) in enumerate(s.neighbours) - element_i = s.atoms[i].name - element_j = s.atoms[j].name - element_pair = string(sort([element_i, element_j])...) - pi = get_pair_params(element_i,element_i,inter) # inter.params[inter.element_map[element_i],:] - pj = get_pair_params(element_j,element_j,inter) # inter.params[inter.element_map[element_j],:] - pij = get_pair_params(element_i,element_j,inter) # inter.params[inter.element_map[element_pair],:] - - r_vec = vector(s.coords[i], s.coords[j], s.box_size) - r2 = sum(abs2, r_vec) - r = sqrt(r2) -# r = rs[n] - r2 = r^2 - r_vec_norm = r_vec / r -# r_vec_norm = r_vec_norms[[n],:] - - # pair contribution - f_pair = - ∂Upair_∂r(r, pij.c, pij.c₀, pij.c₁, pij.c₂) -# dpairdR = f_pair * r_vec_norm -# dpairdR_j = - dpairdR_i - - # glue contribution - dudρ_i = ∂Uglue_∂ρ(ρs[i], pi.A) - dudρ_j = ∂Uglue_∂ρ(ρs[j], pj.A) - dΦdr_i = ∂glue_∂r(r, pi.β, pi.d) - dΦdr_j = ∂glue_∂r(r, pj.β, pj.d) - - f_glue = - (dudρ_j * dΦdr_i + dudρ_i * dΦdr_j) -# println("\nr ", r, " f_pair ", f_pair, " f_glue ", f_glue, " f_i ", f_glue+f_pair) -# dgluedR = f_glue * r_vec_norm -# dgluedR_j = - dgluedR_i -# ## density change by moving the current atom -# dgluedR_i_curr = r_vec_norm * dudρ_i * dΦdr_j -# dgluedR_j_curr = r_vec_norm * dudρ_j * dΦdr_i -# ## density change by moving a neighbouring atom -# dgluedR_i_neigh = - r_vec_norm * dudρ_j * dΦdr_i -# dgluedR_j_neigh = - r_vec_norm * dudρ_i * dΦdr_j - - # updating the forces -# f_i = (dpairdR_i + dgluedR_i_curr + dgluedR_i_neigh) -# f_j = (dpairdR_j + dgluedR_j_curr + dgluedR_j_neigh) - f_i = (f_glue + f_pair) * r_vec_norm #dpairdR + dgluedR - fs[i] += f_i - fs[j] -= f_i - end -# println("\nfs \n", fs) - return collect(1:n_atoms), fs -end - -@inline @inbounds function potential_energy(inter::FinnisSinclair, s::Simulation) - #logger - specific inters - computing the potential energy combining glue and pair components. - - e_pair = 0. - e_glue = 0. - n_atoms = length(s.coords) - ρs = zeros(n_atoms) - for (i,j) in s.neighbours - element_i = s.atoms[i].name - element_j = s.atoms[j].name - element_pair = string(sort([element_i, element_j])...) - pi = get_pair_params(element_i,element_i,inter) # inter.params[inter.element_map[element_i],:] - pj = get_pair_params(element_j,element_j,inter) # inter.params[inter.element_map[element_j],:] - pij = get_pair_params(element_i,element_j,inter) # inter.params[inter.element_map[element_pair],:] - - r_vec = vector(s.coords[i], s.coords[j], s.box_size) - r2 = sum(abs2, r_vec) - r = sqrt(r2) - - e_pair += Upair(r, pij.c, pij.c₀, pij.c₁, pij.c₂) - - ρs[i] += glue(r, pj.β, pj.d) - ρs[j] += glue(r, pi.β, pi.d) - end - - es_glue = zeros(n_atoms) - for (i, atom) in enumerate(s.atoms) - A = get_pair_params(atom.name, atom.name, inter).A - es_glue[i] = Uglue(ρs[i], A) - end - e_glue = sum(es_glue) - return e_pair + e_glue -end - -=# \ No newline at end of file From fd7cb6f18f93093cd32fc478610185c4a75d9f87 Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Tue, 16 Mar 2021 21:47:42 +0100 Subject: [PATCH 18/49] refactored the ground state energy, forces and vacancy formation energy tests --- test/glue.jl | 272 ++++++++++----------------------------------------- 1 file changed, 54 insertions(+), 218 deletions(-) diff --git a/test/glue.jl b/test/glue.jl index 06c88ae1..f6e9d5f5 100644 --- a/test/glue.jl +++ b/test/glue.jl @@ -3,254 +3,90 @@ using Test using Crystal using DataFrames -elements = ["V", "Nb", "Ta", "Cr", "Mo", "W", "Fe"] -element_pairings = [string(el,el) for el in elements] -element_pair_map = Dict(pair => i for (i,pair) in enumerate(element_pairings)) +kb = Molly.kb +fs_inter, elements, masses, bcc_lattice_constants, reference_energies = Molly.get_finnissinclair1984(true) -df = DataFrame( - element_pair = element_pairings, - d = [3.692767, 3.915354, 4.076980, 3.915720, 4.114825, 4.400224, 3.699579], - A = [2.010637, 3.013789, 2.591061, 1.453418, 1.887117, 1.896373, 1.889846], - β = [0, 0, 0, 1.8, 0, 0, 1.8], - c = [3.8, 4.2, 4.2, 2.9, 3.25, 3.25, 3.4], - c₀ = [-0.8816318, -1.5640104, 1.2157373, 29.1429813, 43.4475218, 47.1346499, 1.2110601], - c₁ = [1.4907756, 2.0055779, 0.0271471, -23.3975027, -31.9332978, -33.7665655, -0.7510840], - c₂ = [-0.3976370, -0.4663764, -0.1217350, 4.7578297, 6.0804249, 6.2541999, 0.1380773], -) - -masses = Dict("V" => 50.9415, "Nb" => 92.9064, "Ta" => 180.9479, - "Cr" => 51.996, "Mo" => 95.94, "W" => 183.85, - "Fe" => 55.847) - -# Å -bcc_lattice_constants = Dict( - "V" => 3.0399, "Nb" => 3.3008, - "Ta" => 3.3058, "Cr" => 2.8845, "Mo" => 3.1472, - "W" => 3.1652, "Fe" => 2.8665 -) - -reference_energies = DataFrame( - element_pair = element_pairings, - u = [5.31, 7.57, 8.1, 4.1, 6.82, 8.9, 4.28], - u_vac = [1.92, 2.64, 3.13, 1.97, 2.58, 3.71, 1.77] -) - -fs_inter = FinnisSinclair(true, element_pair_map, df) - -function test_forces_zero(forces, n_atoms; dims=3) - zeros = [zero(rand(3)) for _ in 1:n_atoms] - return all(isapprox.(forces, zeros, atol=1e-4)) +function forces_are_zero(forces; dims=3) + return all([isapprox(f, zero(rand(3)), atol=1e-4) for f in forces]) end -function test_forces_for_element(element::String, fs_inter; nx::Integer=3, ny::Integer=3, nz::Integer=3, ) - a = bcc_lattice_constants[element] - atoms, coords, box, box_size, box_vectors = Crystal.make_bcc_unitcell(element, a=a) - sc_atoms, sc_coords, sc_box, sc_box_size = Crystal.make_supercell(atoms, coords, box, box_size, nx=nx, ny=ny, nz=nz) - n_atoms = length(sc_atoms) - - specific_inter_list = ((fs_inter,),) - velocities = [velocity(1., .01, dims=3) for i in 1:n_atoms] - sim = VelocityVerlet() - nb_matrix = trues(n_atoms,n_atoms) - n_steps = 1 - dist_cutoff = 2 * a - - nf = DistanceNeighbourFinder(nb_matrix, n_steps, dist_cutoff) - - loggers = Dict("temperature" => TemperatureLogger(1)) - - s = Simulation( - simulator=sim, - atoms=sc_atoms, - specific_inter_lists=specific_inter_list, - general_inters=(), - coords=[SVector{3}(v) for v in sc_coords], - velocities=velocities, - temperature=.01, - box_size=sc_box_size[1,1], - timestep=.002, - n_steps=1, - neighbour_finder=nf, - loggers=loggers, - ) - find_neighbours!(s, s.neighbour_finder, 0) - forces = accelerations(s, parallel=false) - return test_forces_zero(forces, n_atoms) +function glue_remains_reasonable(glue_init,glue_end) + return all(isapprox.(glue_init,glue_end,rtol=.1)) end -function test_bcc_vacancy_forces(fs_inter;element::String="Fe",nx::Int64=3,ny::Int64=3,nz::Int64=3,) - a = bcc_lattice_constants[element] - atoms, coords, box, box_size, box_vectors = Crystal.make_bcc_unitcell(element, a=a) - sc_atoms, sc_coords, sc_box, sc_box_size = Crystal.make_supercell(atoms, coords, box, box_size, nx=nx, ny=ny, nz=nz) - n_atoms = length(sc_atoms) - - # introducing a vacancy - sc_atoms_vac, sc_coords_vac = Crystal.add_vacancies(sc_atoms, sc_coords, ixs=[1]) - n_atoms_vac = length(sc_atoms_vac) - - specific_inter_list = ((fs_inter,),) - velocities = [velocity(1., .01, dims=3) for i in 1:n_atoms_vac] - sim = VelocityVerlet() - nb_matrix = trues(n_atoms_vac,n_atoms_vac) - dist_cutoff = 2 * a - n_steps = 1 - nf = DistanceNeighbourFinder(nb_matrix, n_steps, dist_cutoff) - - loggers = Dict("temperature" => TemperatureLogger(1)) +function groundstate_energy_as_expected(s, element, inter) + element_pair = string(element, element) + row = reference_energies[inter.element_pair_map[element_pair],:] + return isapprox(s.loggers["energy"].energies[1]/length(s.coords), -row.u, atol=1e-2) +end - s = Simulation( - simulator=sim, - atoms=sc_atoms_vac, - specific_inter_lists=specific_inter_list, - general_inters=(), - coords=[SVector{3}(v) for v in sc_coords_vac], - velocities=velocities, - temperature=.01, - box_size=sc_box_size[1,1], - timestep=.002, - n_steps=1, - neighbour_finder=nf, - loggers=loggers, - ) - find_neighbours!(s, s.neighbour_finder, 0) - forces = accelerations(s, parallel=false) - return !test_forces_zero(forces, n_atoms_vac) +function vacancy_formation_energy_as_expected(s, s_vac, element, inter) + element_pair = string(element, element) + row = reference_energies[inter.element_pair_map[element_pair],:] + u_gs = s.loggers["energy"].energies[1] + n_atoms = length(s.coords) + u_vac = s_vac.loggers["energy"].energies[1] + n_atoms_vac = length(s_vac.coords) + u_vac_form = u_vac - n_atoms_vac/n_atoms * u_gs + return isapprox(u_vac_form, row.u_vac, atol=7e-2) end -function test_energies_for_element(element::String, fs_inter, u; nx::Integer=3, ny::Integer=3, nz::Integer=3, ) +function run_bcc(element::String, fs_inter, T::Real=.01*kb; + nx::Integer=3, ny::Integer=3, nz::Integer=3, vac::Bool=false) + a = bcc_lattice_constants[element] atoms, coords, box, box_size, box_vectors = Crystal.make_bcc_unitcell(element, a=a) - sc_atoms, sc_coords, sc_box, sc_box_size = Crystal.make_supercell(atoms, coords, box, box_size, nx=nx, ny=ny, nz=nz) + sc_atoms, sc_coords, sc_box, sc_box_size = Crystal.make_supercell(atoms, coords, + box, box_size, nx=nx, ny=ny, nz=nz) + + if vac + # introducing a vacancy + sc_atoms, sc_coords = Crystal.add_vacancies(sc_atoms, sc_coords, ixs=[1]) + end n_atoms = length(sc_atoms) - specific_inter_list = ((fs_inter,),) - velocities = [velocity(1., .01, dims=3) for i in 1:n_atoms] - sim = VelocityVerlet() + velocities = [velocity(1., T, dims=3) for i in 1:n_atoms] nb_matrix = trues(n_atoms,n_atoms) - n_steps = 1 dist_cutoff = 2 * a - nf = DistanceNeighbourFinder(nb_matrix, n_steps, dist_cutoff); - - loggers = Dict("temperature" => TemperatureLogger(1)) - - s = Simulation( - simulator=sim, - atoms=sc_atoms, - specific_inter_lists=specific_inter_list, - general_inters=(), - coords=[SVector{3}(v) for v in sc_coords], - velocities=velocities, - temperature=.01, - box_size=sc_box_size[1,1], - timestep=.002, - n_steps=1, - neighbour_finder=nf, - loggers=loggers, + loggers = Dict( + "glue" => GlueDensityLogger(1), + "forces" => ForceLogger(5), + "energy" => EnergyLogger(1) ) - find_neighbours!(s, s.neighbour_finder, 0) - u_md = Molly.potential_energy(fs_inter, s)/n_atoms - return isapprox(u_md, u, atol=1e-2) -end - -function bcc_vacancy_formation_energy(fs_inter;element::String="Fe",nx::Int64=3,ny::Int64=3,nz::Int64=3,) - a = bcc_lattice_constants[element] - atoms, coords, box, box_size, box_vectors = Crystal.make_bcc_unitcell(element, a=a) - sc_atoms, sc_coords, sc_box, sc_box_size = Crystal.make_supercell(atoms, coords, box, box_size, nx=nx, ny=ny, nz=nz) - n_atoms = length(sc_atoms) - - # energy of the vacancy free system - specific_inter_list = ((fs_inter,),) - velocities = [velocity(1., .01, dims=3) for i in 1:n_atoms] - sim = VelocityVerlet() - nb_matrix = trues(n_atoms,n_atoms) - n_steps = 1 - dist_cutoff = 2 * a - - nf = DistanceNeighbourFinder(nb_matrix, n_steps, dist_cutoff); - - loggers = Dict("temperature" => TemperatureLogger(1)) + nf = DistanceNeighbourFinder(nb_matrix, 1, dist_cutoff) + s = Simulation( - simulator=sim, + simulator=VelocityVerlet(), atoms=sc_atoms, - specific_inter_lists=specific_inter_list, - general_inters=(), + general_inters=(fs_inter,), coords=[SVector{3}(v) for v in sc_coords], velocities=velocities, - temperature=.01, + temperature=T, box_size=sc_box_size[1,1], timestep=.002, - n_steps=1, + n_steps=100, neighbour_finder=nf, loggers=loggers, ) - find_neighbours!(s, s.neighbour_finder, 0) - u_md = Molly.potential_energy(fs_inter, s) - - # introducing a vacancy - sc_atoms_vac, sc_coords_vac = Crystal.add_vacancies(sc_atoms, sc_coords, ixs=[1]) - n_atoms_vac = length(sc_atoms_vac) - - velocities = [velocity(1., .01, dims=3) for i in 1:n_atoms_vac] - nb_matrix = trues(n_atoms_vac,n_atoms_vac) - - nf = DistanceNeighbourFinder(nb_matrix, n_steps, dist_cutoff); - - loggers = Dict("temperature" => TemperatureLogger(1)) - - s = Simulation( - simulator=sim, - atoms=sc_atoms_vac, - specific_inter_lists=specific_inter_list, - general_inters=(), - coords=[SVector{3}(v) for v in sc_coords_vac], - velocities=velocities, - temperature=.01, - box_size=sc_box_size[1,1], - timestep=.002, - n_steps=1, - neighbour_finder=nf, - loggers=loggers, - ) - find_neighbours!(s, s.neighbour_finder, 0) - u_md_vac = Molly.potential_energy(fs_inter, s) - return u_md_vac - n_atoms_vac/n_atoms * u_md + simulate!(s) + return s end - @testset "Finnis-Sinclair" begin - @testset "test forces" begin - for element in elements - @test test_forces_for_element(element, fs_inter) - end - end - - @testset "bcc + vacancy: test forces" begin - for element in elements - @test test_bcc_vacancy_forces(fs_inter, element=element) - end - end - - @testset "potential energies" begin - for element in elements - element_pair = string(element, element) - row = reference_energies[fs_inter.element_pair_map[element_pair],:] - @testset "$element" begin - @test test_energies_for_element(element, fs_inter, -row.u) - end + for element in elements + @testset "$element" begin + s = run_bcc(element, fs_inter, .01*kb) + @test forces_are_zero(s.loggers["forces"].forces[1]) + @test glue_remains_reasonable(s.loggers["glue"].glue_densities[1], + s.loggers["glue"].glue_densities[end]) + @test groundstate_energy_as_expected(s, element, fs_inter) + s_vac = run_bcc(element, fs_inter, .01*kb; vac=true) + @test !forces_are_zero(s_vac.loggers["forces"].forces[1]) + @test vacancy_formation_energy_as_expected(s, s_vac, element, fs_inter) end end - @testset "potential energies: bcc vacancy" begin - for element in elements - element_pair = string(element, element) - row = reference_energies[fs_inter.element_pair_map[element_pair],:] - u_vac = bcc_vacancy_formation_energy(fs_inter, element=element) - @testset "$element" begin - @test isapprox(u_vac, row.u_vac, atol=7e-2) - end - end - end -end \ No newline at end of file +end From f80004e85a44fb9ec5f91aaddd9451d0b942b0e9 Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Tue, 16 Mar 2021 22:41:59 +0100 Subject: [PATCH 19/49] refactored fs.ipynb: removed superfluous tests and reduced the math to the core components --- fs.ipynb | 877 +++++++++---------------------------------------------- 1 file changed, 132 insertions(+), 745 deletions(-) diff --git a/fs.ipynb b/fs.ipynb index 99857c7c..066add81 100644 --- a/fs.ipynb +++ b/fs.ipynb @@ -11,11 +11,10 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Content\n", - "1. Finnis-Sinclair forces for perfect bcc crystals\n", - "2. Finnis-Sinclair potential energies for perfect bcc crystals\n", - "3. Finnis-Sinclair forces for bcc crystals with 1 vacancy\n", - "4. Finnis-Sinclair potential energies for bcc crystals with 1 vacancy" + "Content:\n", + "1. the math of a glue potential in general & what's special in the Finnis-Sinclair version\n", + "2. show components of the model (pair energy, glue density, glue energy)\n", + "3. show how to run a simulation" ] }, { @@ -23,8 +22,8 @@ "metadata": {}, "source": [ "TODOs:\n", - "* check if MD trajectories around room temperature are expected (the crystals should stay relatively stable and not explode - currently it seems they evaporate)\n", - "* Find a more elegant solution for the force implementation that does not require altering `function accelerations`" + "* clean up this notebook: add some descriptions\n", + "* generate animation of vibrating lattice" ] }, { @@ -56,50 +55,16 @@ "# Pkg.add(url=\"https://github.com/eschmidt42/crystal\")" ] }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "# Pkg.add(\"Makie\")" - ] - }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "Pkg.status()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "using Molly" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ + "using Molly\n", "using DataFrames\n", - "# using Molly\n", "using Plots\n", "using Test\n", - "# using LaTeXStrings\n", - "# using LinearAlgebra\n", - "# using SparseArrays\n", "using Crystal" ] }, @@ -107,134 +72,66 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Parameterisation of V, Nb, Ta, Cr, Mo, W, Fe" + "## The math 💖" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Parameterisation by Finnis et al. 1984, _A simple empirical N-body potential for transition metals_" + "$$\n", + "U = U_{\\text{pair}} + U_{\\text{glue}}\n", + "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "| element | d | A | $\\beta$ | c | $c_0$ | $c_1$ | $c_2$ | \n", - "| --- | --- | --- | --- | --- | --- | --- | --- |\n", - "| V | 3.692767 | 2.010637 | 0 | 3.8 | -0.8816318 | 1.4907756 | -0.3976370 |\n", - "| Nb | 3.915354 | 3.013789 | 0 | 4.2 | -1.5640104 | 2.0055779 | -0.4663764 |\n", - "| Ta | 4.076980 | 2.591061 | 0 | 4.2 | 1.2157373 | 0.0271471 | -0.1217350 |\n", - "| Cr | 3.915720 | 1.453418 | 1.8 | 2.9 | 29.1429813 | -23.3975027 | 4.7578297 |\n", - "| Mo | 4.114825 | 1.887117 | 0 | 3.25 | 43.4475218 | -31.9332978 | 6.0804249 |\n", - "| W | 4.400224 | 1.896373 | 0 | 3.25 | 47.1346499 | -33.7665655 | 6.2541999 |\n", - "| Fe | 3.699579 | 1.889846 | 1.8 | 3.4 | 1.2110601 | -0.7510840 | 0.1380773 |" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "elements = [\"V\", \"Nb\", \"Ta\", \"Cr\", \"Mo\", \"W\", \"Fe\"]\n", - "element_pairings = [string(el,el) for el in elements]\n", - "element_pair_map = Dict(pair => i for (i,pair) in enumerate(element_pairings))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "df = DataFrame(\n", - " element_pair = element_pairings,\n", - " d = [3.692767, 3.915354, 4.076980, 3.915720, 4.114825, 4.400224, 3.699579],\n", - " A = [2.010637, 3.013789, 2.591061, 1.453418, 1.887117, 1.896373, 1.889846],\n", - " β = [0, 0, 0, 1.8, 0, 0, 1.8],\n", - " c = [3.8, 4.2, 4.2, 2.9, 3.25, 3.25, 3.4],\n", - " c₀ = [-0.8816318, -1.5640104, 1.2157373, 29.1429813, 43.4475218, 47.1346499, 1.2110601],\n", - " c₁ = [1.4907756, 2.0055779, 0.0271471, -23.3975027, -31.9332978, -33.7665655, -0.7510840],\n", - " c₂ = [-0.3976370, -0.4663764, -0.1217350, 4.7578297, 6.0804249, 6.2541999, 0.1380773],\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "masses = Dict(\"V\" => 50.9415, \"Nb\" => 92.9064, \"Ta\" => 180.9479,\n", - " \"Cr\" => 51.996, \"Mo\" => 95.94, \"W\" => 183.85,\n", - " \"Fe\" => 55.847)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Å\n", - "bcc_lattice_constants = Dict(\n", - " \"V\" => 3.0399, \"Nb\" => 3.3008, \n", - " \"Ta\" => 3.3058, \"Cr\" => 2.8845, \"Mo\" => 3.1472, \n", - " \"W\" => 3.1652, \"Fe\" => 2.8665\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "reference_energies = DataFrame(\n", - " element_pair = element_pairings,\n", - " u = [5.31, 7.57, 8.1, 4.1, 6.82, 8.9, 4.28],\n", - " u_vac = [1.92, 2.64, 3.13, 1.97, 2.58, 3.71, 1.77]\n", - ")" + "$$\n", + "U_{\\text{pair}} = \\sum_{i,j>i}V_{ij}(r_{ij})\n", + "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Interaction" + "$r_{ij} = \\left\\| \\vec{R}_j - \\vec{R}_i| \\right\\|_2$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Instantiating the interaction" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fs_inter = FinnisSinclair(true, element_pair_map, df)" + "$$\n", + "U_{\\text{glue}} = \\sum_i f(\\rho_i)\n", + "$$\n", + "\n", + "$\\rho_i = \\sum_j \\phi(r_{ij})$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Glue potential" + "$$\n", + "\\begin{split}\n", + "U &= U_{\\text{pair}} + U_{\\text{glue}} \\\\\n", + " &= \\sum_{i,j>i}V_{ij}(r_{ij}) + \\sum_i f(\\rho_i) \\\\\n", + "\\end{split}\n", + "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "The glue potential is the core component which makes the Finnis-Sinclair empirical potential and other similar approaches different to, for example, the Lennard-Jones potential. \n", - "TODO: some more explanation" + "$$\n", + "\\vec{F} = - \\partial_i U\n", + "$$\n", + "\n", + "$\\partial_i$ is the wiggling of $\\vec{R}_i$" ] }, { @@ -242,309 +139,61 @@ "metadata": {}, "source": [ "$$\n", - "\\phi(r) = (r-d)^2 + \\beta (r-d)^3/d\n", + "\\begin{split}\n", + "\\vec{F}_i &= - \\sum_j V^\\prime_{ij}(r_{ij}) \\left( -\\frac{\\vec{R}_{ij}}{r_{ij}} \\right) + f^\\prime_i(\\rho_i) \\sum_j \\phi_{j}^\\prime(r_{ij}) \\left( -\\frac{\\vec{R}_{ij}}{r_{ij}} \\right) + \\sum_j f^\\prime_j(\\rho_j) \\phi_i(r_{ji}) \\left( \\frac{\\vec{R}_{ji}}{r_{ji}} \\right) \\\\\n", + " &= \\sum_j \\left[ V^\\prime_{ij}(r_{ij}) + f^\\prime_i(\\rho_i) \\phi_{j}^\\prime(r_{ij}) + f^\\prime_j(\\rho_j) \\phi_i(r_{ji}) \\right] \\frac{\\vec{R}_{ij}}{r_{ij}}\n", + "\\end{split}\n", "$$" ] }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "element_pair, d, A, β, c, c₀, c₁, c₂ = df[1,:] # parameters for Vanadium" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "r = collect(range(0, stop=2*3.3058, length=1000));\n", - "ɸ = Molly.glue_potential.(r, β, d);" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "ɸs = [ɸ]\n", - "element_pairs = [element_pair]\n", - "for i in 2:nrow(df)\n", - " element_pair, d, A, β, c, c₀, c₁, c₂ = df[i,:]\n", - " ɸ = Molly.glue_potential.(r, β, d)\n", - " append!(ɸs,[ɸ])\n", - " element_pairs = hcat(element_pairs, string(element_pair))\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plot(r, ɸs, label=element_pairs)" - ] - }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Computing forces" + "## The model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "This is what we need for the simulation." + "We are going to look into $f$, $\\phi$ and $V$ " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "$\\partial_k = \\partial_{R_k} = $ change of atom position $k$, $r_{ij} = \\|R_{ij}\\|_2$, $R_{ij} = R_j - R_i \\in \\mathbb{R}^3$ \n", - "\n", - "$$\n", - "\\partial_k u = \\partial_k u_\\text{pair} + \\partial_k u_\\text{glue} \n", - "$$\n", - "\n", - "$$\n", - "\\partial_k u_\\text{pair} = \\sum_{i>j} V_{ij}^\\prime(r_{ij})\\partial_k r_{ij}\n", - "$$\n", - "\n", - "$$\n", - "\\partial_k u_\\text{glue} = \\sum_i f_i^\\prime(\\rho_i) \\cdot \\partial_k \\rho_i \n", - "$$\n", - "\n", - "$$\n", - "V_{ij}^\\prime(r_{ij}) = 2(r-c)(c_0 + c_1 r + c_2 r^2) + (r-c)^2(c_1 + 2c_2r)\n", - "$$\n", - "\n", - "$$\n", - "f_i^\\prime(\\rho_i) \\cdot \\partial_k \\rho_i =\n", - "\\begin{cases}\n", - " k = i, & f_k^\\prime(\\rho_k) \\sum_j\\phi_j^\\prime(r_{kj})\\partial_k r_{kj} \\\\\n", - " k \\ne i, & \\sum_{i\\ne k} f_i^\\prime(\\rho_i) \\partial_k \\phi_k^\\prime(r_{ik})\\partial_k r_{ik} \\\\\n", - "\\end{cases}\n", - "$$\n", - "\n", - "$$\n", - "f_i^\\prime = \\frac{1}{2}A_i\\rho_i^{-3/2}\n", - "$$\n", - "\n", - "$$\n", - "\\partial_k\\phi(r) = \\left[2(r-d) + 3\\beta (r-d)^2/d\\right] \\cdot\n", - "\\begin{cases}\n", - " k = i, &\\frac{R_{kj}}{r_{kj}} \\\\ \n", - " k = j, &\\frac{R_{ik}}{r_{ik}} \\\\\n", - "\\end{cases}\n", - "$$\n", - "\n", - "$$\n", - "\\partial_k r_{ij} =\n", - "\\begin{cases}\n", - " k = i, &\\frac{R_{kj}}{r_{kj}} \\\\ \n", - " k = j, &\\frac{R_{ik}}{r_{ik}} \\\\\n", - "\\end{cases}\n", - "$$" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "element = \"Fe\"\n", - "a = bcc_lattice_constants[element]\n", - "atoms, coords, box, box_size, box_vectors = Crystal.make_bcc_unitcell(element, a=a)\n", - "sc_atoms, sc_coords, sc_box, sc_box_size = Crystal.make_supercell(atoms, coords, box, box_size, nx=3, ny=3,nz=3)\n", - "n_atoms = length(sc_atoms);" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "specific_inter_list = ((fs_inter,),)\n", - "velocities = [velocity(1., .01, dims=3) for i in 1:n_atoms]\n", - "sim = VelocityVerlet()\n", - "nb_matrix = trues(n_atoms,n_atoms)\n", - "n_steps = 1\n", - "dist_cutoff = 2 * a\n", - "\n", - "nf = DistanceNeighbourFinder(nb_matrix, n_steps, dist_cutoff);\n", - "\n", - "loggers = Dict(\"temperature\" => TemperatureLogger(1))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "s = Simulation(\n", - " simulator=sim, \n", - " atoms=sc_atoms, \n", - " specific_inter_lists=specific_inter_list,\n", - " general_inters=(),\n", - " coords=[SVector{3}(v) for v in sc_coords], \n", - " velocities=velocities,\n", - " temperature=.01, \n", - " box_size=sc_box_size[1,1],\n", - " timestep=.002,\n", - " n_steps=10,\n", - " neighbour_finder=nf,\n", - " loggers=loggers,\n", - ")\n", - "find_neighbours!(s, s.neighbour_finder, 0)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "sparse_forces = force.((fs_inter,), (s.coords,), (s,))" + "units https://lammps.sandia.gov/doc/units.html" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Testing that all forces are about 0" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "function test_sparse_forces_zero(sparse_forces, n_atoms; dims=3)\n", - " zeros = [zero(rand(1,3)) for _ in 1:n_atoms]\n", - " forces = getindex.(sparse_forces,2)[1]\n", - " return all(isapprox.(forces, zeros, atol=1e-6))\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "@test test_sparse_forces_zero(sparse_forces, n_atoms)" + "### Parameterisations" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Modifying `Molly.accelerations` so the forces from the glue interaction are properly used to update the atom forces." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "function test_forces_zero(forces, n_atoms; dims=3)\n", - " zeros = [zero(rand(3)) for _ in 1:n_atoms]\n", - " return all(isapprox.(forces, zeros, atol=1e-4))\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "@test test_forces_zero(accelerations(s, parallel=false), n_atoms)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "simulate!(s, parallel=false)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "s.loggers[\"temperature\"].temperatures" + "Parameterisation by Finnis et al. 1984, _A simple empirical N-body potential for transition metals_" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Testing forces for all elements" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "function test_forces_for_element(element::String, fs_inter; nx::Integer=3, ny::Integer=3, nz::Integer=3, )\n", - " a = bcc_lattice_constants[element]\n", - " atoms, coords, box, box_size, box_vectors = Crystal.make_bcc_unitcell(element, a=a)\n", - " sc_atoms, sc_coords, sc_box, sc_box_size = Crystal.make_supercell(atoms, coords, box, box_size, nx=nx, ny=ny, nz=nz)\n", - " n_atoms = length(sc_atoms)\n", - " \n", - " specific_inter_list = ((fs_inter,),)\n", - " velocities = [velocity(1., .01, dims=3) for i in 1:n_atoms]\n", - " sim = VelocityVerlet()\n", - " nb_matrix = trues(n_atoms,n_atoms)\n", - " n_steps = 1\n", - " dist_cutoff = 2 * a\n", - "\n", - " nf = DistanceNeighbourFinder(nb_matrix, n_steps, dist_cutoff);\n", - "\n", - " loggers = Dict(\"temperature\" => TemperatureLogger(1))\n", - " \n", - " s = Simulation(\n", - " simulator=sim, \n", - " atoms=sc_atoms, \n", - " specific_inter_lists=specific_inter_list,\n", - " general_inters=(),\n", - " coords=[SVector{3}(v) for v in sc_coords], \n", - " velocities=velocities,\n", - " temperature=.01, \n", - " box_size=sc_box_size[1,1],\n", - " timestep=.002,\n", - " n_steps=1,\n", - " neighbour_finder=nf,\n", - " loggers=loggers,\n", - " )\n", - " find_neighbours!(s, s.neighbour_finder, 0)\n", - " forces = accelerations(s, parallel=false)\n", - " return test_forces_zero(forces, n_atoms)\n", - "end" + "| element | d | A | $\\beta$ | c | $c_0$ | $c_1$ | $c_2$ | \n", + "| --- | --- | --- | --- | --- | --- | --- | --- |\n", + "| V | 3.692767 | 2.010637 | 0 | 3.8 | -0.8816318 | 1.4907756 | -0.3976370 |\n", + "| Nb | 3.915354 | 3.013789 | 0 | 4.2 | -1.5640104 | 2.0055779 | -0.4663764 |\n", + "| Ta | 4.076980 | 2.591061 | 0 | 4.2 | 1.2157373 | 0.0271471 | -0.1217350 |\n", + "| Cr | 3.915720 | 1.453418 | 1.8 | 2.9 | 29.1429813 | -23.3975027 | 4.7578297 |\n", + "| Mo | 4.114825 | 1.887117 | 0 | 3.25 | 43.4475218 | -31.9332978 | 6.0804249 |\n", + "| W | 4.400224 | 1.896373 | 0 | 3.25 | 47.1346499 | -33.7665655 | 6.2541999 |\n", + "| Fe | 3.699579 | 1.889846 | 1.8 | 3.4 | 1.2110601 | -0.7510840 | 0.1380773 |" ] }, { @@ -553,32 +202,23 @@ "metadata": {}, "outputs": [], "source": [ - "@testset \"test forces\" begin\n", - " for element in elements\n", - " @test test_forces_for_element(element, fs_inter)\n", - " end\n", - "end" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Computing energies" + "kb = Molly.kb\n", + "fs84, elements, masses, bcc_lattice_constants, reference_energies = Molly.get_finnissinclair1984(true)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "This is only really interesting for logging / development of potentials." + "### Glue contribution - $\\phi(r)$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Pair energy" + "The glue potential is the core component which makes the Finnis-Sinclair empirical potential and other similar approaches different to, for example, the Lennard-Jones potential. \n", + "TODO: some more explanation" ] }, { @@ -586,12 +226,8 @@ "metadata": {}, "source": [ "$$\n", - "V_{ij}(r_{ij}) = \n", - "\\begin{cases} \n", - "r \\le c, & (r-c)^2 \\left( c_0 + c_1 r + c_2 r^2 \\right) \\\\\n", - "r > c, & 0 \\\\\n", - "\\end{cases}\n", - "$$\n" + "\\phi(r) = (r-d)^2 + \\beta (r-d)^3/d\n", + "$$" ] }, { @@ -600,7 +236,7 @@ "metadata": {}, "outputs": [], "source": [ - "element_pair, d, A, β, c, c₀, c₁, c₂ = df[1,:] # parameters for Vanadium" + "element_pair, d, A, β, c, c₀, c₁, c₂ = fs84.params[1,:] # parameters for Vanadium" ] }, { @@ -609,7 +245,8 @@ "metadata": {}, "outputs": [], "source": [ - "V = Molly.pair_potential.(r, c, c₀, c₁, c₂);" + "r = collect(range(0, stop=6, length=1000));\n", + "ɸ = Molly.glue.(r, β, d);" ] }, { @@ -618,12 +255,12 @@ "metadata": {}, "outputs": [], "source": [ - "Vs = [V]\n", + "ɸs = [ɸ]\n", "element_pairs = [element_pair]\n", - "for i in 2:nrow(df)\n", - " element_pair, d, A, β, c, c₀, c₁, c₂ = df[i,:]\n", - " V = Molly.pair_potential.(r, c, c₀, c₁, c₂)\n", - " append!(Vs,[V])\n", + "for i in 2:nrow(fs84.params)\n", + " element_pair, d, A, β, c, c₀, c₁, c₂ = fs84.params[i,:]\n", + " ɸ = Molly.glue.(r, β, d)\n", + " append!(ɸs,[ɸ])\n", " element_pairs = hcat(element_pairs, string(element_pair))\n", "end" ] @@ -634,14 +271,14 @@ "metadata": {}, "outputs": [], "source": [ - "plot(r, Vs, label=element_pairs)" + "plot(r, ɸs, label=element_pairs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Glue energy" + "### Glue energy - $f(\\rho)$" ] }, { @@ -671,7 +308,7 @@ "outputs": [], "source": [ "ρ = 4. # density that you get summing phi-contributions from neighbours\n", - "Molly.glue_energy(ρ, 1.)" + "Molly.Uglue(ρ, .15)" ] }, { @@ -682,9 +319,9 @@ "source": [ "ρ = collect(range(0, stop=50, length=100));\n", "\n", - "A = df.A[1] # Va\n", - "uₙ = Molly.glue_energy.(ρ, A)\n", - "element_pair = df.element_pair[1]" + "A = fs84.params.A[1] # Va\n", + "uₙ = Molly.Uglue.(ρ, A)\n", + "element_pair = fs84.params.element_pair[1]" ] }, { @@ -695,9 +332,9 @@ "source": [ "uₙs = [uₙ]\n", "element_pairs = [element_pair]\n", - "for i in 2:nrow(df)\n", - " element_pair, d, A, β, c, c₀, c₁, c₂ = df[i,:]\n", - " uₙ = Molly.glue_energy.(ρ, A)\n", + "for i in 2:nrow(fs84.params)\n", + " element_pair, d, A, β, c, c₀, c₁, c₂ = fs84.params[i,:]\n", + " uₙ = Molly.Uglue.(ρ, A)\n", " append!(uₙs,[uₙ])\n", " element_pairs = hcat(element_pairs, string(element_pair))\n", "end" @@ -716,18 +353,20 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "### Pair + glue energy = magic" + "### Pair energy - $V(r)$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "$$ u_\\text{tot} = u_N + u_P $$\n", - "\n", - "$$ u_P = \\frac{1}{2}\\sum_{i=1,j=1}^{n_\\text{atoms},n_\\text{atoms}} V(r_{ij}) $$\n", - "\n", - "$$ u_N = \\sum_{i=1}^{n_\\text{atoms}} u_\\text{glue}(\\rho_i) $$" + "$$\n", + "V_{ij}(r_{ij}) = \n", + "\\begin{cases} \n", + "r \\le c, & (r-c)^2 \\left( c_0 + c_1 r + c_2 r^2 \\right) \\\\\n", + "r > c, & 0 \\\\\n", + "\\end{cases}\n", + "$$\n" ] }, { @@ -736,7 +375,7 @@ "metadata": {}, "outputs": [], "source": [ - "u = Molly.potential_energy(fs_inter, s) / n_atoms" + "element_pair, d, A, β, c, c₀, c₁, c₂ = fs84.params[1,:] # parameters for Vanadium" ] }, { @@ -745,14 +384,7 @@ "metadata": {}, "outputs": [], "source": [ - "@test isapprox(u, -4.28, atol=1e-2)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Testing potential energies for all elements" + "V = Molly.Upair.(r, c, c₀, c₁, c₂);" ] }, { @@ -761,40 +393,13 @@ "metadata": {}, "outputs": [], "source": [ - "function test_energies_for_element(element::String, fs_inter, u; nx::Integer=3, ny::Integer=3, nz::Integer=3, )\n", - " a = bcc_lattice_constants[element]\n", - " atoms, coords, box, box_size, box_vectors = Crystal.make_bcc_unitcell(element, a=a)\n", - " sc_atoms, sc_coords, sc_box, sc_box_size = Crystal.make_supercell(atoms, coords, box, box_size, nx=nx, ny=ny, nz=nz)\n", - " n_atoms = length(sc_atoms)\n", - " \n", - " specific_inter_list = ((fs_inter,),)\n", - " velocities = [velocity(1., .01, dims=3) for i in 1:n_atoms]\n", - " sim = VelocityVerlet()\n", - " nb_matrix = trues(n_atoms,n_atoms)\n", - " n_steps = 1\n", - " dist_cutoff = 2 * a\n", - "\n", - " nf = DistanceNeighbourFinder(nb_matrix, n_steps, dist_cutoff);\n", - "\n", - " loggers = Dict(\"temperature\" => TemperatureLogger(1))\n", - " \n", - " s = Simulation(\n", - " simulator=sim, \n", - " atoms=sc_atoms, \n", - " specific_inter_lists=specific_inter_list,\n", - " general_inters=(),\n", - " coords=[SVector{3}(v) for v in sc_coords], \n", - " velocities=velocities,\n", - " temperature=.01, \n", - " box_size=sc_box_size[1,1],\n", - " timestep=.002,\n", - " n_steps=1,\n", - " neighbour_finder=nf,\n", - " loggers=loggers,\n", - " )\n", - " find_neighbours!(s, s.neighbour_finder, 0)\n", - " u_md = Molly.potential_energy(fs_inter, s)/n_atoms\n", - " return isapprox(u_md, u, atol=1e-2)\n", + "Vs = [V]\n", + "element_pairs = [element_pair]\n", + "for i in 2:nrow(fs84.params)\n", + " element_pair, d, A, β, c, c₀, c₁, c₂ = fs84.params[i,:]\n", + " V = Molly.Upair.(r, c, c₀, c₁, c₂)\n", + " append!(Vs,[V])\n", + " element_pairs = hcat(element_pairs, string(element_pair))\n", "end" ] }, @@ -804,245 +409,14 @@ "metadata": {}, "outputs": [], "source": [ - "@testset \"potential energies\" begin \n", - " for element in elements\n", - " element_pair = string(element, element)\n", - " row = reference_energies[fs_inter.element_pair_map[element_pair],:]\n", - " @testset \"$element\" begin \n", - " @test test_energies_for_element(element, fs_inter, -row.u) \n", - " end\n", - " end\n", - "end" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Running simulation with potential energy logger" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "loggers = Dict(\n", - " \"temperature\" => TemperatureLogger(1),\n", - " \"energy\" => EnergyLogger(1),\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "s = Simulation(\n", - " simulator=sim, \n", - " atoms=sc_atoms, \n", - " specific_inter_lists=specific_inter_list,\n", - " general_inters=(),\n", - " coords=[SVector{3}(v) for v in sc_coords], \n", - " velocities=velocities,\n", - " temperature=.01, \n", - " box_size=sc_box_size[1,1],\n", - " timestep=.002,\n", - " n_steps=10,\n", - " neighbour_finder=nf,\n", - " loggers=loggers,\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "simulate!(s, parallel=false)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "s.loggers[\"energy\"].energies" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Vacancies" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Ensuring that the vacancy formation energies are within .07 eV of the reference." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "function bcc_vacancy_formation_energy(fs_inter;element::String=\"Fe\",nx::Int64=3,ny::Int64=3,nz::Int64=3,)\n", - " a = bcc_lattice_constants[element]\n", - " atoms, coords, box, box_size, box_vectors = Crystal.make_bcc_unitcell(element, a=a)\n", - " sc_atoms, sc_coords, sc_box, sc_box_size = Crystal.make_supercell(atoms, coords, box, box_size, nx=nx, ny=ny, nz=nz)\n", - " n_atoms = length(sc_atoms)\n", - " \n", - " # energy of the vacancy free system\n", - " specific_inter_list = ((fs_inter,),)\n", - " velocities = [velocity(1., .01, dims=3) for i in 1:n_atoms]\n", - " sim = VelocityVerlet()\n", - " nb_matrix = trues(n_atoms,n_atoms)\n", - " n_steps = 1\n", - " dist_cutoff = 2 * a\n", - "\n", - " nf = DistanceNeighbourFinder(nb_matrix, n_steps, dist_cutoff);\n", - "\n", - " loggers = Dict(\"temperature\" => TemperatureLogger(1))\n", - "\n", - " s = Simulation(\n", - " simulator=sim, \n", - " atoms=sc_atoms, \n", - " specific_inter_lists=specific_inter_list,\n", - " general_inters=(),\n", - " coords=[SVector{3}(v) for v in sc_coords], \n", - " velocities=velocities,\n", - " temperature=.01, \n", - " box_size=sc_box_size[1,1],\n", - " timestep=.002,\n", - " n_steps=1,\n", - " neighbour_finder=nf,\n", - " loggers=loggers,\n", - " )\n", - " find_neighbours!(s, s.neighbour_finder, 0)\n", - " u_md = Molly.potential_energy(fs_inter, s)\n", - " \n", - " # introducing a vacancy\n", - " sc_atoms_vac, sc_coords_vac = Crystal.add_vacancies(sc_atoms, sc_coords, ixs=[1])\n", - " n_atoms_vac = length(sc_atoms_vac)\n", - " \n", - " velocities = [velocity(1., .01, dims=3) for i in 1:n_atoms_vac]\n", - " nb_matrix = trues(n_atoms_vac,n_atoms_vac)\n", - "\n", - " nf = DistanceNeighbourFinder(nb_matrix, n_steps, dist_cutoff);\n", - "\n", - " loggers = Dict(\"temperature\" => TemperatureLogger(1))\n", - "\n", - " s = Simulation(\n", - " simulator=sim, \n", - " atoms=sc_atoms_vac, \n", - " specific_inter_lists=specific_inter_list,\n", - " general_inters=(),\n", - " coords=[SVector{3}(v) for v in sc_coords_vac], \n", - " velocities=velocities,\n", - " temperature=.01, \n", - " box_size=sc_box_size[1,1],\n", - " timestep=.002,\n", - " n_steps=1,\n", - " neighbour_finder=nf,\n", - " loggers=loggers,\n", - " )\n", - " find_neighbours!(s, s.neighbour_finder, 0)\n", - " u_md_vac = Molly.potential_energy(fs_inter, s)\n", - " return u_md_vac - n_atoms_vac/n_atoms * u_md\n", - "end\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "@testset \"potential energies: bcc vacancy\" begin \n", - " for element in elements\n", - " element_pair = string(element, element)\n", - " row = reference_energies[fs_inter.element_pair_map[element_pair],:]\n", - " u_vac = bcc_vacancy_formation_energy(fs_inter, element=element)\n", - " @testset \"$element\" begin \n", - " @test isapprox(u_vac, row.u_vac, atol=7e-2) \n", - " end\n", - " end\n", - "end" + "plot(r, Vs, label=element_pairs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Testing that the forces in an unrelaxed crystal with a vacancy are not all zero." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "function test_bcc_vacancy_forces(fs_inter;element::String=\"Fe\",nx::Int64=3,ny::Int64=3,nz::Int64=3,)\n", - " a = bcc_lattice_constants[element]\n", - " atoms, coords, box, box_size, box_vectors = Crystal.make_bcc_unitcell(element, a=a)\n", - " sc_atoms, sc_coords, sc_box, sc_box_size = Crystal.make_supercell(atoms, coords, box, box_size, nx=nx, ny=ny, nz=nz)\n", - " n_atoms = length(sc_atoms)\n", - " \n", - " # introducing a vacancy\n", - " sc_atoms_vac, sc_coords_vac = Crystal.add_vacancies(sc_atoms, sc_coords, ixs=[1])\n", - " n_atoms_vac = length(sc_atoms_vac)\n", - " \n", - " velocities = [velocity(1., .01, dims=3) for i in 1:n_atoms_vac]\n", - " nb_matrix = trues(n_atoms_vac,n_atoms_vac)\n", - "\n", - " nf = DistanceNeighbourFinder(nb_matrix, n_steps, dist_cutoff);\n", - "\n", - " loggers = Dict(\"temperature\" => TemperatureLogger(1))\n", - "\n", - " s = Simulation(\n", - " simulator=sim, \n", - " atoms=sc_atoms_vac, \n", - " specific_inter_lists=specific_inter_list,\n", - " general_inters=(),\n", - " coords=[SVector{3}(v) for v in sc_coords_vac], \n", - " velocities=velocities,\n", - " temperature=.01, \n", - " box_size=sc_box_size[1,1],\n", - " timestep=.002,\n", - " n_steps=1,\n", - " neighbour_finder=nf,\n", - " loggers=loggers,\n", - " )\n", - " find_neighbours!(s, s.neighbour_finder, 0)\n", - " forces = accelerations(s, parallel=false)\n", - " return !test_forces_zero(forces, n_atoms_vac)\n", - "end\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "@testset \"bcc + vacancy: test forces\" begin\n", - " for element in elements\n", - " @test test_bcc_vacancy_forces(fs_inter, element=element)\n", - " end\n", - "end" + "## Simulating" ] }, { @@ -1054,40 +428,47 @@ "nx = 3\n", "ny = 3\n", "nz = 3\n", - "element = \"Fe\"\n", - "T = 300\n", - "n_steps = 1_000\n", - "dt = .002\n", + "element = \"W\"\n", + "T = 300. # Kelvin\n", + "T *= kb\n", + "n_steps = 5000\n", + "dt = .002 # ps; ns = 1e-9, ps = 1e-12, fs = 1e-15\n", "\n", "a = bcc_lattice_constants[element]\n", "atoms, coords, box, box_size, box_vectors = Crystal.make_bcc_unitcell(element, a=a)\n", "sc_atoms, sc_coords, sc_box, sc_box_size = Crystal.make_supercell(atoms, coords, box, box_size, nx=nx, ny=ny, nz=nz)\n", "n_atoms = length(sc_atoms)\n", "\n", - "# introducing a vacancy\n", - "sc_atoms_vac, sc_coords_vac = Crystal.add_vacancies(sc_atoms, sc_coords, ixs=[1])\n", - "n_atoms_vac = length(sc_atoms_vac)\n", - "\n", - "specific_inter_list = ((fs_inter,),)\n", - "velocities = [velocity(T, sc_atoms_vac[i].mass, dims=3) for i in 1:n_atoms_vac]\n", - "sim = VelocityVerlet()\n", - "nb_matrix = trues(n_atoms_vac,n_atoms_vac)\n", + "general_inters = (fs_inter,)\n", + "velocities = [velocity(sc_atoms[i].mass, T, dims=3) for i in 1:n_atoms]\n", + "nb_matrix = trues(n_atoms,n_atoms)\n", "dist_cutoff = 2 * a\n", - "nf = DistanceNeighbourFinder(nb_matrix, n_steps, dist_cutoff)\n", - "\n", + "nf = DistanceNeighbourFinder(nb_matrix, 1, dist_cutoff)\n", + "# thermostat = NoThermostat()\n", + "thermostat = AndersenThermostat(1.)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ "loggers = Dict(\n", " \"temperature\" => TemperatureLogger(1),\n", " \"pot\" => EnergyLogger(1),\n", - " \"coords\" => CoordinateLogger(5),\n", - " \"writer\" => StructureWriter(5,\"traj_fs_bcc_vac_fe.pdb\")\n", + " \"coords\" => CoordinateLogger(100),\n", + " \"velos\" => VelocityLogger(1),\n", + " \"forces\" => ForceLogger(100),\n", + " \"glue\" => GlueDensityLogger(1),\n", + "# \"writer\" => StructureWriter(5,\"traj_fs_bcc_vac_fe.pdb\")\n", ")\n", "\n", "s = Simulation(\n", - " simulator=sim, \n", - " atoms=sc_atoms_vac, \n", - " specific_inter_lists=specific_inter_list,\n", - " general_inters=(),\n", - " coords=[SVector{3}(v) for v in sc_coords_vac], \n", + " simulator=VelocityVerlet(), \n", + " atoms=sc_atoms, \n", + " general_inters=general_inters,\n", + " coords=[SVector{3}(v) for v in sc_coords], \n", " velocities=velocities,\n", " temperature=T, \n", " box_size=sc_box_size[1,1],\n", @@ -1101,9 +482,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "scrolled": false - }, + "metadata": {}, "outputs": [], "source": [ "simulate!(s)" @@ -1115,7 +494,9 @@ "metadata": {}, "outputs": [], "source": [ - "s.loggers[\"temperature\"].temperatures" + "x = collect(1:length(s.loggers[\"glue\"].glue_densities))\n", + "y = [[arr[i] for arr in s.loggers[\"glue\"].glue_densities] for i in 1:n_atoms]\n", + "plot(x, y, title=\"glue densities\", legend=false)" ] }, { @@ -1124,7 +505,9 @@ "metadata": {}, "outputs": [], "source": [ - "s.loggers[\"pot\"].energies" + "x = collect(1:length(s.loggers[\"temperature\"].temperatures))\n", + "y = s.loggers[\"temperature\"].temperatures/kb\n", + "plot(x,y,title=\"Temperature\",legend=false)" ] }, { @@ -1132,7 +515,11 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": [] + "source": [ + "x = collect(1:length(s.loggers[\"pot\"].energies))\n", + "y = s.loggers[\"pot\"].energies / length(s.coords)\n", + "plot(x,y,title=\"Potential energy per atom\",legend=false)" + ] } ], "metadata": { From 275c29658d9c8814f8fc4658c0dd8d8d5cbdea2e Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Wed, 17 Mar 2021 19:01:43 +0100 Subject: [PATCH 20/49] removed superfluous enumerate --- src/interactions/glue_fs.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/interactions/glue_fs.jl b/src/interactions/glue_fs.jl index b6986910..01e1f4ee 100644 --- a/src/interactions/glue_fs.jl +++ b/src/interactions/glue_fs.jl @@ -161,7 +161,7 @@ function update_glue_densities!( end # updating glue densities - for (n,(i,j)) in enumerate(s.neighbours) + for (i,j) in s.neighbours # collecting parameters el_i = s.atoms[i].name From b136903116ca06e2b08e9917a7560283848456f6 Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Wed, 17 Mar 2021 19:03:04 +0100 Subject: [PATCH 21/49] added docstrings --- fs.ipynb | 237 +++++++++++++++++++++++++++++++------------------------ 1 file changed, 135 insertions(+), 102 deletions(-) diff --git a/fs.ipynb b/fs.ipynb index 066add81..ed8800c0 100644 --- a/fs.ipynb +++ b/fs.ipynb @@ -7,6 +7,15 @@ "# The Finnis-Sinclair potential" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Reading:\n", + "- Finnis and Sinclair 1984, [A simple empirical N-body potential for transition metals](https://doi.org/10.1080/01418618408244210)\n", + "- Daw and Baskes 1984, [Embedded-atom method: Derivation and application to impurities, surfaces, and other defects in metals](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.29.6443)" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -22,7 +31,6 @@ "metadata": {}, "source": [ "TODOs:\n", - "* clean up this notebook: add some descriptions\n", "* generate animation of vibrating lattice" ] }, @@ -32,35 +40,12 @@ "metadata": {}, "outputs": [], "source": [ - "import Pkg" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "Pkg.activate(\".\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Pkg.add(url=\"https://github.com/eschmidt42/crystal\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ + "import Pkg\n", + "\n", + "Pkg.activate(\".\")\n", + "\n", + "# Pkg.add(url=\"https://github.com/eschmidt42/crystal\")\n", + "\n", "using Molly\n", "using DataFrames\n", "using Plots\n", @@ -79,6 +64,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ + "The potential energy is composed out of a pair and a \"glue\" contribution\n", "$$\n", "U = U_{\\text{pair}} + U_{\\text{glue}}\n", "$$" @@ -88,6 +74,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ + "The pair contribution\n", "$$\n", "U_{\\text{pair}} = \\sum_{i,j>i}V_{ij}(r_{ij})\n", "$$" @@ -97,24 +84,26 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "$r_{ij} = \\left\\| \\vec{R}_j - \\vec{R}_i| \\right\\|_2$" + "where $r_{ij} = \\left\\| \\vec{R}_j - \\vec{R}_i| \\right\\|_2$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ + "The energy contribution of the glue is taken for each atom\n", "$$\n", "U_{\\text{glue}} = \\sum_i f(\\rho_i)\n", "$$\n", "\n", - "$\\rho_i = \\sum_j \\phi(r_{ij})$" + "where $\\rho_i = \\sum_j \\phi(r_{ij})$ is the glue at atom $i$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ + "Combining both\n", "$$\n", "\\begin{split}\n", "U &= U_{\\text{pair}} + U_{\\text{glue}} \\\\\n", @@ -127,23 +116,27 @@ "cell_type": "markdown", "metadata": {}, "source": [ + "To compute the force on atom i with this potential we apply\n", "$$\n", - "\\vec{F} = - \\partial_i U\n", + "\\vec{F}_i = - \\partial_i U\n", "$$\n", "\n", - "$\\partial_i$ is the wiggling of $\\vec{R}_i$" + "where $\\partial_i$ is the wiggling of $\\vec{R}_i$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ + "Applying the chain rule we obtain\n", "$$\n", "\\begin{split}\n", "\\vec{F}_i &= - \\sum_j V^\\prime_{ij}(r_{ij}) \\left( -\\frac{\\vec{R}_{ij}}{r_{ij}} \\right) + f^\\prime_i(\\rho_i) \\sum_j \\phi_{j}^\\prime(r_{ij}) \\left( -\\frac{\\vec{R}_{ij}}{r_{ij}} \\right) + \\sum_j f^\\prime_j(\\rho_j) \\phi_i(r_{ji}) \\left( \\frac{\\vec{R}_{ji}}{r_{ji}} \\right) \\\\\n", " &= \\sum_j \\left[ V^\\prime_{ij}(r_{ij}) + f^\\prime_i(\\rho_i) \\phi_{j}^\\prime(r_{ij}) + f^\\prime_j(\\rho_j) \\phi_i(r_{ji}) \\right] \\frac{\\vec{R}_{ij}}{r_{ij}}\n", "\\end{split}\n", - "$$" + "$$\n", + "\n", + "The interesting part is the second equation, because this we can directly use to implement the force updates, looping the neighbouring pairs in parallel ($j$ is a neighbouring atom of atom $i$). " ] }, { @@ -157,14 +150,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We are going to look into $f$, $\\phi$ and $V$ " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "units https://lammps.sandia.gov/doc/units.html" + "Let's look into $f$, $\\phi$ and $V$. " ] }, { @@ -197,37 +183,47 @@ ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "kb = Molly.kb\n", - "fs84, elements, masses, bcc_lattice_constants, reference_energies = Molly.get_finnissinclair1984(true)" + "Note that this potential assumes *metal* units, e.g. eV, K and so on: https://lammps.sandia.gov/doc/units.html" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Glue contribution - $\\phi(r)$" + "Collecting the parameters" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "kb = Molly.kb \n", + "fs84, elements, masses, bcc_lattice_constants, reference_energies = Molly.get_finnissinclair1984(true)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "The glue potential is the core component which makes the Finnis-Sinclair empirical potential and other similar approaches different to, for example, the Lennard-Jones potential. \n", - "TODO: some more explanation" + "### Glue contribution - $\\phi(r)$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ + "To compute the glue density $\\rho$ for an atom, we need to compute the contributions from its neighbours (self-contribution is ignored). The neighbour contributions are obtained using $\\phi(r)$:\n", + "\n", "$$\n", "\\phi(r) = (r-d)^2 + \\beta (r-d)^3/d\n", - "$$" + "$$\n", + "\n", + "$d$ and $\\beta$ are model parameters." ] }, { @@ -271,7 +267,23 @@ "metadata": {}, "outputs": [], "source": [ - "plot(r, ɸs, label=element_pairs)" + "?plot" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plot(r, ɸs, label=element_pairs, xlabel=\"r\", ylabel=\"phi(r)\", title=\"Glue contributions\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Interestingly we find that some $\\phi$ do indeed take on negative values. This is interesting because of the choice for $f$ by Finnis and Sinclair." ] }, { @@ -302,13 +314,12 @@ ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "ρ = 4. # density that you get summing phi-contributions from neighbours\n", - "Molly.Uglue(ρ, .15)" + "From a physical point of view it makes sense to expect $\\rho > 0$, but in some situations we still might obtain $\\rho < 0$, and have a problem, unless we allow for complex energy values.\n", + "\n", + "Let's plot $u_{\\text{glue}}$ for $\\rho \\in \\mathbb{R}^+$." ] }, { @@ -317,19 +328,11 @@ "metadata": {}, "outputs": [], "source": [ - "ρ = collect(range(0, stop=50, length=100));\n", - "\n", + "ρ = collect(range(0, stop=50, length=100))\n", "A = fs84.params.A[1] # Va\n", "uₙ = Molly.Uglue.(ρ, A)\n", - "element_pair = fs84.params.element_pair[1]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ + "element_pair = fs84.params.element_pair[1];\n", + "\n", "uₙs = [uₙ]\n", "element_pairs = [element_pair]\n", "for i in 2:nrow(fs84.params)\n", @@ -346,7 +349,7 @@ "metadata": {}, "outputs": [], "source": [ - "plot(ρ, uₙs, label=element_pairs)" + "plot(ρ, uₙs, label=element_pairs, xlabel=\"Glue density\", ylabel=\"Glue energy\", title=\"Glue energy varying with glue density\")" ] }, { @@ -360,6 +363,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ + "The pair energy is defined as\n", "$$\n", "V_{ij}(r_{ij}) = \n", "\\begin{cases} \n", @@ -370,21 +374,10 @@ ] }, { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "element_pair, d, A, β, c, c₀, c₁, c₂ = fs84.params[1,:] # parameters for Vanadium" - ] - }, - { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "V = Molly.Upair.(r, c, c₀, c₁, c₂);" + "Let's have a look" ] }, { @@ -393,6 +386,10 @@ "metadata": {}, "outputs": [], "source": [ + "element_pair, d, A, β, c, c₀, c₁, c₂ = fs84.params[1,:] # parameters for Vanadium\n", + "\n", + "V = Molly.Upair.(r, c, c₀, c₁, c₂)\n", + "\n", "Vs = [V]\n", "element_pairs = [element_pair]\n", "for i in 2:nrow(fs84.params)\n", @@ -400,23 +397,30 @@ " V = Molly.Upair.(r, c, c₀, c₁, c₂)\n", " append!(Vs,[V])\n", " element_pairs = hcat(element_pairs, string(element_pair))\n", - "end" + "end\n", + "\n", + "plot(r, Vs, label=element_pairs, xlabel=\"r\", ylabel=\"Pair energy contribution\", title=\"Pair energy vs separation\")" ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "plot(r, Vs, label=element_pairs)" + "## Simulating" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Simulating" + "Now that we have looked at the individual components of the model, let's simulate Tungsten near room temperature for a $3\\times3\\times3$ supercell of a body centered cubic crystal." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Creating the crystal" ] }, { @@ -429,23 +433,18 @@ "ny = 3\n", "nz = 3\n", "element = \"W\"\n", - "T = 300. # Kelvin\n", - "T *= kb\n", - "n_steps = 5000\n", - "dt = .002 # ps; ns = 1e-9, ps = 1e-12, fs = 1e-15\n", "\n", "a = bcc_lattice_constants[element]\n", "atoms, coords, box, box_size, box_vectors = Crystal.make_bcc_unitcell(element, a=a)\n", "sc_atoms, sc_coords, sc_box, sc_box_size = Crystal.make_supercell(atoms, coords, box, box_size, nx=nx, ny=ny, nz=nz)\n", - "n_atoms = length(sc_atoms)\n", - "\n", - "general_inters = (fs_inter,)\n", - "velocities = [velocity(sc_atoms[i].mass, T, dims=3) for i in 1:n_atoms]\n", - "nb_matrix = trues(n_atoms,n_atoms)\n", - "dist_cutoff = 2 * a\n", - "nf = DistanceNeighbourFinder(nb_matrix, 1, dist_cutoff)\n", - "# thermostat = NoThermostat()\n", - "thermostat = AndersenThermostat(1.)" + "n_atoms = length(sc_atoms)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Setting simulation specs" ] }, { @@ -454,6 +453,19 @@ "metadata": {}, "outputs": [], "source": [ + "T = 300. # Kelvin\n", + "T *= kb\n", + "n_steps = 5000\n", + "dt = .002 # ps; ns = 1e-9, ps = 1e-12, fs = 1e-15\n", + "\n", + "general_inters = (fs84,)\n", + "velocities = [velocity(sc_atoms[i].mass, T, dims=3) for i in 1:n_atoms]\n", + "nb_matrix = trues(n_atoms,n_atoms)\n", + "dist_cutoff = 2 * a\n", + "nf = DistanceNeighbourFinder(nb_matrix, 1, dist_cutoff)\n", + "# thermostat = NoThermostat()\n", + "thermostat = AndersenThermostat(1.)\n", + "\n", "loggers = Dict(\n", " \"temperature\" => TemperatureLogger(1),\n", " \"pot\" => EnergyLogger(1),\n", @@ -479,6 +491,13 @@ ")" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Running the simulation (this should take about 30s)" + ] + }, { "cell_type": "code", "execution_count": null, @@ -488,6 +507,13 @@ "simulate!(s)" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Visualizing content of the loggers" + ] + }, { "cell_type": "code", "execution_count": null, @@ -520,6 +546,13 @@ "y = s.loggers[\"pot\"].energies / length(s.coords)\n", "plot(x,y,title=\"Potential energy per atom\",legend=false)" ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { From 931f21454c792e826a182f134d75706ba0fb1122 Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Wed, 17 Mar 2021 19:09:51 +0100 Subject: [PATCH 22/49] cleaned --- fs.ipynb | 8 -------- 1 file changed, 8 deletions(-) diff --git a/fs.ipynb b/fs.ipynb index ed8800c0..349f1938 100644 --- a/fs.ipynb +++ b/fs.ipynb @@ -26,14 +26,6 @@ "3. show how to run a simulation" ] }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "TODOs:\n", - "* generate animation of vibrating lattice" - ] - }, { "cell_type": "code", "execution_count": null, From 4e2c73b3c1f4af88606c7964fe73106df977b15f Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Wed, 17 Mar 2021 19:16:48 +0100 Subject: [PATCH 23/49] removed cell --- fs.ipynb | 9 --------- 1 file changed, 9 deletions(-) diff --git a/fs.ipynb b/fs.ipynb index 349f1938..e19a90cb 100644 --- a/fs.ipynb +++ b/fs.ipynb @@ -253,15 +253,6 @@ "end" ] }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "?plot" - ] - }, { "cell_type": "code", "execution_count": null, From f88c5b2dd4e25cd3437e8fb2961d16998ce7fb32 Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Sat, 20 Mar 2021 04:18:50 +0100 Subject: [PATCH 24/49] removing comments --- src/forces.jl | 6 ------ 1 file changed, 6 deletions(-) diff --git a/src/forces.jl b/src/forces.jl index e7f80349..45d70f49 100644 --- a/src/forces.jl +++ b/src/forces.jl @@ -103,10 +103,6 @@ function accelerations(s::Simulation; parallel::Bool=true) sparse_forces = force.(inter_list, (s.coords,), (s,)) ge1, ge2 = getindex.(sparse_forces, 1), getindex.(sparse_forces, 2) sparse_vec = SparseVector(n_atoms, reduce(vcat, ge1), reduce(vcat, ge2)) -# if typeof(inter_list[1]) <: GlueInteraction -# sparse_vec = [SVector{3}(v) for v in sparse_vec] -# end -# println("\nsparse_vec \n", sparse_vec) forces += Array(sparse_vec) end @@ -115,8 +111,6 @@ function accelerations(s::Simulation; parallel::Bool=true) s.forces[i] = forces[i] end -# println("\nforces \n", forces) - return forces end From 5a8100226bfe5c71ee00dcb0d9f0038b9b24d112 Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Sat, 27 Mar 2021 19:34:33 +0100 Subject: [PATCH 25/49] replaced ForwardDiff.derivate with analytical derivatives and added kb as an attribute of FinnisSinclair --- src/interactions/glue_fs.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/interactions/glue_fs.jl b/src/interactions/glue_fs.jl index 01e1f4ee..edf0efd6 100644 --- a/src/interactions/glue_fs.jl +++ b/src/interactions/glue_fs.jl @@ -18,11 +18,10 @@ struct FinnisSinclair <: GlueInteraction nl_only::Bool element_pair_map::Dict params::DataFrame + kb::Real end -kb = 8.617333262145e-5 # eV/K - """ get_finnissinclair1984(nl_only) @@ -68,8 +67,8 @@ function get_finnissinclair1984(nl_only::Bool) "Fe" => 2.8665 ) - - fs84 = FinnisSinclair(nl_only, element_pair_map, df) + kb = 8.617333262145e-5 # eV/K + fs84 = FinnisSinclair(nl_only, element_pair_map, df, kb) reference_energies = DataFrame( element_pair = element_pairings, @@ -96,7 +95,7 @@ end Derivative of the glue density function. """ -∂glue_∂r(r, β, d) = ForwardDiff.derivative(r -> glue(r,β,d), r) +∂glue_∂r(r, β, d) = 2*(r-d) + 3*β*(r-d)^2/d """ Uglue(ρ, A) @@ -113,7 +112,8 @@ end Energy derivative given glue density. """ -∂Uglue_∂ρ(ρ,A) = ForwardDiff.derivative(ρ -> Uglue(ρ,A), ρ) +∂Uglue_∂ρ(ρ,A) = - A / (2 * √ρ) + """ Upair(r, c, c₀, c₁, c₂; f=10.) From c0d37fa46938c679772208687964e2bfe4aeafc7 Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Sat, 27 Mar 2021 19:43:03 +0100 Subject: [PATCH 26/49] removed f comment in docstrings --- src/interactions/glue_fs.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/interactions/glue_fs.jl b/src/interactions/glue_fs.jl index edf0efd6..609b82e3 100644 --- a/src/interactions/glue_fs.jl +++ b/src/interactions/glue_fs.jl @@ -84,7 +84,6 @@ end The core component of the Finnis-Sinclair type GlueInteraction. Used to calculate contribution to the glue value of an atom, based on its neighbouring atoms. -f is a fudge factor to help translate a Å model to a nm model. """ function glue(r, β, d) return r > d ? 0 : (r-d)^2 + β*(r-d)^3/d @@ -128,7 +127,6 @@ end ∂Upair_∂r(r, c, c₀, c₁, c₂, f=10.) Derivative of the pair potential. -f is a fudge factor to help translate a Å model to a nm model. """ ∂Upair_∂r(r, c, c₀,c₁, c₂) = ForwardDiff.derivative(r -> Upair(r,c,c₀,c₁,c₂), r) From a71dc512a0a222ba903ec1d9abf140a68ac38727 Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Sat, 27 Mar 2021 20:50:46 +0100 Subject: [PATCH 27/49] replaced dataframes with dicts --- src/interactions/glue_fs.jl | 115 ++++++++++++++++++------------------ 1 file changed, 57 insertions(+), 58 deletions(-) diff --git a/src/interactions/glue_fs.jl b/src/interactions/glue_fs.jl index 609b82e3..93ee551c 100644 --- a/src/interactions/glue_fs.jl +++ b/src/interactions/glue_fs.jl @@ -8,16 +8,31 @@ SpecificInteraction. """ abstract type GlueInteraction <: SpecificInteraction end +struct FinnisSinclairPair + # container for two element interaction -> c c0 c1 c2 + c::Real + c₀::Real + c₁::Real + c₂::Real +end + +struct FinnisSinclairSingle + # container for single element stuff: glue density, glue energy -> A + A::Real + β::Real + d::Real +end + """ - FinnisSinclairInteraction(nl_only,element_pair_map,params) + FinnisSinclairInteraction(nl_only,pairs,singles,kb) The Finnis-Sinclair interaction. This interaction expects units to be of these https://lammps.sandia.gov/doc/units.html units (eV, Å, K, ps and so on). """ struct FinnisSinclair <: GlueInteraction nl_only::Bool - element_pair_map::Dict - params::DataFrame + pairs::Dict{String,FinnisSinclairPair} + singles::Dict{String,FinnisSinclairSingle} kb::Real end @@ -30,31 +45,26 @@ Finnis and Sinclair 1984 parameterization: https://doi.org/10.1080/0141861840824 function get_finnissinclair1984(nl_only::Bool) elements = ["V", "Nb", "Ta", "Cr", "Mo", "W", "Fe"] - element_pairings = [string(el, el) for el in elements] - element_pair_map = Dict(pair => i for (i,pair) in enumerate(element_pairings)) - df = DataFrame( - element_pair = element_pairings, - d = [3.692767, 3.915354, 4.076980, 3.915720, - 4.114825, 4.400224, 3.699579], - # (Å) - A = [2.010637, 3.013789, 2.591061, 1.453418, - 1.887117, 1.896373, 1.889846], - # (eV) - β = [0, 0, 0, 1.8, 0, 0, 1.8], - # (1) - c = [3.8, 4.2, 4.2, 2.9, 3.25, 3.25, 3.4], - # (Å) - c₀ = [-0.8816318, -1.5640104, 1.2157373, 29.1429813, - 43.4475218, 47.1346499, 1.2110601], - # (1) - c₁ = [1.4907756, 2.0055779, 0.0271471, -23.3975027, - -31.9332978, -33.7665655, -0.7510840], - # (1) - c₂ = [-0.3976370, -0.4663764, -0.1217350, 4.7578297, - 6.0804249, 6.2541999, 0.1380773], - # (1) - ) + d = [3.692767, 3.915354, 4.076980, 3.915720, + 4.114825, 4.400224, 3.699579] # (Å) + A = [2.010637, 3.013789, 2.591061, 1.453418, + 1.887117, 1.896373, 1.889846] # (eV) + β = [0, 0, 0, 1.8, 0, 0, 1.8] # (1) + c = [3.8, 4.2, 4.2, 2.9, 3.25, 3.25, 3.4] # (Å) + c₀ = [-0.8816318, -1.5640104, 1.2157373, 29.1429813, + 43.4475218, 47.1346499, 1.2110601] # (1) + c₁ = [1.4907756, 2.0055779, 0.0271471, -23.3975027, + -31.9332978, -33.7665655, -0.7510840] # (1) + c₂ = [-0.3976370, -0.4663764, -0.1217350, 4.7578297, + 6.0804249, 6.2541999, 0.1380773] # (1) + + fs_singles = Dict() + fs_pairs = Dict() + for (i,el) in enumerate(elements) + fs_singles[el] = FinnisSinclairSingle(A[i],β[i],d[i]) + fs_pairs[string(el,el)] = FinnisSinclairPair(c[i],c₀[i],c₁[i],c₂[i]) + end masses = Dict("V" => 50.9415, "Nb" => 92.9064, "Ta" => 180.9479, "Cr" => 51.996, "Mo" => 95.94, "W" => 183.85, @@ -68,13 +78,13 @@ function get_finnissinclair1984(nl_only::Bool) ) kb = 8.617333262145e-5 # eV/K - fs84 = FinnisSinclair(nl_only, element_pair_map, df, kb) + fs84 = FinnisSinclair(nl_only, fs_pairs, fs_singles, kb) - reference_energies = DataFrame( - element_pair = element_pairings, - u = [5.31, 7.57, 8.1, 4.1, 6.82, 8.9, 4.28], - u_vac = [1.92, 2.64, 3.13, 1.97, 2.58, 3.71, 1.77] - ) + u = [5.31, 7.57, 8.1, 4.1, 6.82, 8.9, 4.28] # eV/atom + u_vac = [1.92, 2.64, 3.13, 1.97, 2.58, 3.71, 1.77] # eV/atom + + reference_energies = Dict(el=>Dict("u"=>u[i], "u_vac"=>u_vac[i]) + for (i,el) in enumerate(elements)) return fs84, elements, masses, bcc_lattice_constants, reference_energies end @@ -94,7 +104,7 @@ end Derivative of the glue density function. """ -∂glue_∂r(r, β, d) = 2*(r-d) + 3*β*(r-d)^2/d +∂glue_∂r(r, β, d) = r > d ? 0 : 2*(r-d) + 3*β*(r-d)^2/d """ Uglue(ρ, A) @@ -130,16 +140,6 @@ Derivative of the pair potential. """ ∂Upair_∂r(r, c, c₀,c₁, c₂) = ForwardDiff.derivative(r -> Upair(r,c,c₀,c₁,c₂), r) -""" - get_pair_params(element1, element2, inter) - -Convenience function to generate element pair and return relevant model parameters. -""" -function get_pair_params(element1::String, element2::String, inter::FinnisSinclair) - pair = string(sort([element1, element2])...) - return inter.params[inter.element_pair_map[pair],:] -end - """ update_glue_densities!(inter, coords, s, parallel) @@ -164,9 +164,9 @@ function update_glue_densities!( # collecting parameters el_i = s.atoms[i].name el_j = s.atoms[j].name - pi = get_pair_params(el_i,el_i,inter) - pj = get_pair_params(el_j,el_j,inter) - pij = get_pair_params(el_i,el_j,inter) + pi = inter.singles[el_i] + pj = inter.singles[el_j] + pij = inter.pairs[string(el_i,el_j)] # computing distance dr = vector(s.coords[i], s.coords[j], s.box_size) @@ -189,12 +189,12 @@ end @inline @inbounds function force(inter::FinnisSinclair, s::Simulation, i::Integer, j::Integer) # parameters - element_i = s.atoms[i].name - element_j = s.atoms[j].name - element_pair = string(sort([element_i, element_j])...) - pi = get_pair_params(element_i, element_i, inter) - pj = get_pair_params(element_j, element_j, inter) - pij = get_pair_params(element_i, element_j, inter) + el_i = s.atoms[i].name + el_j = s.atoms[j].name + el_ij = string(el_i,el_j) + pi = inter.singles[el_i] + pj = inter.singles[el_j] + pij = inter.pairs[el_ij] # distance i -> j dr = vector(s.coords[i], s.coords[j], s.box_size) @@ -225,15 +225,14 @@ end if no_glue update_glue_densities!(inter, s.coords, s, parallel=false) end - - A = get_pair_params(s.atoms[i].name, s.atoms[i].name, inter).A - return Uglue(s.glue_densities[i], A) + + return Uglue(s.glue_densities[i], inter.singles[s.atoms[i].name].A) end @inline @inbounds function potential_energy(inter::FinnisSinclair, s::Simulation, i, j) # logger - general inters - computes the pair energy part only for a single atom pair dr = vector(s.coords[i], s.coords[j], s.box_size) r = norm(dr) - pij = get_pair_params(s.atoms[i].name, s.atoms[j].name, inter) + pij = inter.pairs[string(s.atoms[i].name, s.atoms[j].name)] return Upair(r, pij.c, pij.c₀, pij.c₁, pij.c₂) -end +end \ No newline at end of file From 4b55ae425ceb9a3f0ff59b1d6026c82dd4ed5257 Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Sat, 27 Mar 2021 21:01:14 +0100 Subject: [PATCH 28/49] replaced another derivative forgotten earlier --- src/interactions/glue_fs.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/interactions/glue_fs.jl b/src/interactions/glue_fs.jl index 93ee551c..416aeca5 100644 --- a/src/interactions/glue_fs.jl +++ b/src/interactions/glue_fs.jl @@ -138,7 +138,7 @@ end Derivative of the pair potential. """ -∂Upair_∂r(r, c, c₀,c₁, c₂) = ForwardDiff.derivative(r -> Upair(r,c,c₀,c₁,c₂), r) +∂Upair_∂r(r, c, c₀,c₁, c₂) = (r > c) ? 0 : 2*(r - c) * (c₀ + c₁*r + c₂*r^2) + (r - c)^2 * (c₁ + 2*c₂*r) """ update_glue_densities!(inter, coords, s, parallel) From b332f60b0e8733b06371e60d0bd7523f2a87f6cf Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Sat, 27 Mar 2021 21:01:59 +0100 Subject: [PATCH 29/49] updated tests --- test/glue.jl | 33 ++++++++++++++++++--------------- test/runtests.jl | 2 ++ 2 files changed, 20 insertions(+), 15 deletions(-) diff --git a/test/glue.jl b/test/glue.jl index f6e9d5f5..780259c0 100644 --- a/test/glue.jl +++ b/test/glue.jl @@ -3,7 +3,6 @@ using Test using Crystal using DataFrames -kb = Molly.kb fs_inter, elements, masses, bcc_lattice_constants, reference_energies = Molly.get_finnissinclair1984(true) function forces_are_zero(forces; dims=3) @@ -15,27 +14,30 @@ function glue_remains_reasonable(glue_init,glue_end) end function groundstate_energy_as_expected(s, element, inter) - element_pair = string(element, element) - row = reference_energies[inter.element_pair_map[element_pair],:] - return isapprox(s.loggers["energy"].energies[1]/length(s.coords), -row.u, atol=1e-2) + ref = reference_energies[element] + return isapprox(s.loggers["energy"].energies[1]/length(s.coords), -ref["u"], atol=1e-2) end function vacancy_formation_energy_as_expected(s, s_vac, element, inter) - element_pair = string(element, element) - row = reference_energies[inter.element_pair_map[element_pair],:] + ref = reference_energies[element] u_gs = s.loggers["energy"].energies[1] n_atoms = length(s.coords) u_vac = s_vac.loggers["energy"].energies[1] n_atoms_vac = length(s_vac.coords) u_vac_form = u_vac - n_atoms_vac/n_atoms * u_gs - return isapprox(u_vac_form, row.u_vac, atol=7e-2) + return isapprox(u_vac_form, ref["u_vac"], atol=7e-2) end -function run_bcc(element::String, fs_inter, T::Real=.01*kb; +function run_bcc(element::String, fs_inter, T::Real=.01; nx::Integer=3, ny::Integer=3, nz::Integer=3, vac::Bool=false) + masses = Dict("V" => 50.9415, "Nb" => 92.9064, "Ta" => 180.9479, + "Cr" => 51.996, "Mo" => 95.94, "W" => 183.85, + "Fe" => 55.847) a = bcc_lattice_constants[element] - atoms, coords, box, box_size, box_vectors = Crystal.make_bcc_unitcell(element, a=a) + elements = [element for _ in 1:2] + el2atom_map = Dict(element => Atom(name=element, mass=masses[element])) + atoms, coords, box, box_size, box_vectors = Crystal.make_bcc_unitcell(elements, a=a, el2atom_map=el2atom_map) sc_atoms, sc_coords, sc_box, sc_box_size = Crystal.make_supercell(atoms, coords, box, box_size, nx=nx, ny=ny, nz=nz) @@ -45,7 +47,7 @@ function run_bcc(element::String, fs_inter, T::Real=.01*kb; end n_atoms = length(sc_atoms) - velocities = [velocity(1., T, dims=3) for i in 1:n_atoms] + velocities = [velocity(sc_atoms[i].mass, T*fs_inter.kb, dims=3) for i in 1:n_atoms] nb_matrix = trues(n_atoms,n_atoms) dist_cutoff = 2 * a @@ -63,7 +65,7 @@ function run_bcc(element::String, fs_inter, T::Real=.01*kb; general_inters=(fs_inter,), coords=[SVector{3}(v) for v in sc_coords], velocities=velocities, - temperature=T, + temperature=T*fs_inter.kb, box_size=sc_box_size[1,1], timestep=.002, n_steps=100, @@ -74,16 +76,17 @@ function run_bcc(element::String, fs_inter, T::Real=.01*kb; return s end -@testset "Finnis-Sinclair" begin - +@testset "Glue: Finnis-Sinclair" begin + T = .01 # K for element in elements @testset "$element" begin - s = run_bcc(element, fs_inter, .01*kb) + s = run_bcc(element, fs_inter, T) @test forces_are_zero(s.loggers["forces"].forces[1]) @test glue_remains_reasonable(s.loggers["glue"].glue_densities[1], s.loggers["glue"].glue_densities[end]) @test groundstate_energy_as_expected(s, element, fs_inter) - s_vac = run_bcc(element, fs_inter, .01*kb; vac=true) + + s_vac = run_bcc(element, fs_inter, T; vac=true) @test !forces_are_zero(s_vac.loggers["forces"].forces[1]) @test vacancy_formation_energy_as_expected(s, s_vac, element, fs_inter) end diff --git a/test/runtests.jl b/test/runtests.jl index b3fabeba..2526d63c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -22,6 +22,8 @@ end CUDA.allowscalar(false) # Check that we never do scalar indexing on the GPU +include("glue.jl") + @testset "Spatial" begin @test vector1D(4.0, 6.0, 10.0) == 2.0 @test vector1D(1.0, 9.0, 10.0) == -2.0 From 825f1ad994c75346044481b23aea154154f50ef9 Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Sat, 27 Mar 2021 21:03:07 +0100 Subject: [PATCH 30/49] minor change to improve readability of ForceLogger --- src/loggers.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/loggers.jl b/src/loggers.jl index 17ce7599..3990da06 100644 --- a/src/loggers.jl +++ b/src/loggers.jl @@ -85,7 +85,7 @@ end function ForceLogger(T, n_steps::Integer; dims::Integer=3) return ForceLogger(n_steps, - Array{SArray{Tuple{dims}, T, 1, dims}, 1}[]) + Vector{SVector{dims}}[]) end function ForceLogger(n_steps::Integer; dims::Integer=3) From d646f9bc0a66eb6bf823338527b43b4536f3f645 Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Sat, 27 Mar 2021 21:15:37 +0100 Subject: [PATCH 31/49] removed DataFrames reference --- test/glue.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/test/glue.jl b/test/glue.jl index 780259c0..5edca2f0 100644 --- a/test/glue.jl +++ b/test/glue.jl @@ -1,7 +1,6 @@ using Molly using Test using Crystal -using DataFrames fs_inter, elements, masses, bcc_lattice_constants, reference_energies = Molly.get_finnissinclair1984(true) From 945ddb02fb829d6e130516d24ff5113dd0062be9 Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Sat, 27 Mar 2021 21:16:05 +0100 Subject: [PATCH 32/49] removed DataFrames and Plots reference --- Project.toml | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index b0143bd1..dec14060 100644 --- a/Project.toml +++ b/Project.toml @@ -7,15 +7,13 @@ version = "0.2.2" 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" -DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" +Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" KernelDensity = "5ab0869b-81aa-558d-bb23-cbf5423bbe9b" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce" -Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" Requires = "ae029012-a4dd-5104-9daa-d747884805df" From 115adf2c0b6a87d69437a0046a0e0d1f175c67a9 Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Sat, 27 Mar 2021 21:37:29 +0100 Subject: [PATCH 33/49] relocated packages from Molly.jl to docs --- docs/Project.toml | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/docs/Project.toml b/docs/Project.toml index f3233bec..8a325b7e 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,6 +1,10 @@ [deps] +Crystal = "3c6eccdf-2a89-4c24-a1d4-ff210daa8476" +DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" Molly = "aa0f7f06-fcc0-5ec4-a7f3-a573f33f9c4c" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] Documenter = "0.25" From 30a54b265169056aa3d27866b9df9713cfe389de Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Sat, 27 Mar 2021 21:38:54 +0100 Subject: [PATCH 34/49] made notebook compatible with changes in in src/interactions/glue_fs.jl --- docs/fs.ipynb | 635 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 635 insertions(+) create mode 100644 docs/fs.ipynb diff --git a/docs/fs.ipynb b/docs/fs.ipynb new file mode 100644 index 00000000..976a8fdf --- /dev/null +++ b/docs/fs.ipynb @@ -0,0 +1,635 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# The Finnis-Sinclair potential" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Reading:\n", + "- Finnis and Sinclair 1984, [A simple empirical N-body potential for transition metals](https://doi.org/10.1080/01418618408244210)\n", + "- Daw and Baskes 1984, [Embedded-atom method: Derivation and application to impurities, surfaces, and other defects in metals](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.29.6443)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Content:\n", + "1. the math of a glue potential in general & what's special in the Finnis-Sinclair version\n", + "2. show components of the model (pair energy, glue density, glue energy)\n", + "3. show how to run a simulation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "import Pkg\n", + "\n", + "Pkg.activate(\".\")\n", + "Pkg.instantiate()\n", + "\n", + "# Pkg.add(url=\"https://github.com/eschmidt42/crystal\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "using Molly\n", + "using Plots\n", + "using Test\n", + "using Crystal" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## The math 💖" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The potential energy is composed out of a pair and a \"glue\" contribution\n", + "$$\n", + "U = U_{\\text{pair}} + U_{\\text{glue}}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The pair contribution\n", + "$$\n", + "U_{\\text{pair}} = \\sum_{i,j>i}V_{ij}(r_{ij})\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "where $r_{ij} = \\left\\| \\vec{R}_j - \\vec{R}_i| \\right\\|_2$." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The energy contribution of the glue is taken for each atom\n", + "$$\n", + "U_{\\text{glue}} = \\sum_i f(\\rho_i)\n", + "$$\n", + "\n", + "where $\\rho_i = \\sum_j \\phi(r_{ij})$ is the glue at atom $i$." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Combining both\n", + "$$\n", + "\\begin{split}\n", + "U &= U_{\\text{pair}} + U_{\\text{glue}} \\\\\n", + " &= \\sum_{i,j>i}V_{ij}(r_{ij}) + \\sum_i f(\\rho_i) \\\\\n", + "\\end{split}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To compute the force on atom i with this potential we apply\n", + "$$\n", + "\\vec{F}_i = - \\partial_i U\n", + "$$\n", + "\n", + "where $\\partial_i$ is the wiggling of $\\vec{R}_i$." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Applying the chain rule we obtain\n", + "$$\n", + "\\begin{split}\n", + "\\vec{F}_i &= - \\sum_j V^\\prime_{ij}(r_{ij}) \\left( -\\frac{\\vec{R}_{ij}}{r_{ij}} \\right) + f^\\prime_i(\\rho_i) \\sum_j \\phi_{j}^\\prime(r_{ij}) \\left( -\\frac{\\vec{R}_{ij}}{r_{ij}} \\right) + \\sum_j f^\\prime_j(\\rho_j) \\phi_i(r_{ji}) \\left( \\frac{\\vec{R}_{ji}}{r_{ji}} \\right) \\\\\n", + " &= \\sum_j \\left[ V^\\prime_{ij}(r_{ij}) + f^\\prime_i(\\rho_i) \\phi_{j}^\\prime(r_{ij}) + f^\\prime_j(\\rho_j) \\phi_i(r_{ji}) \\right] \\frac{\\vec{R}_{ij}}{r_{ij}}\n", + "\\end{split}\n", + "$$\n", + "\n", + "The interesting part is the second equation, because this we can directly use to implement the force updates, looping the neighbouring pairs in parallel ($j$ is a neighbouring atom of atom $i$). " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## The model" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's look into $f$, $\\phi$ and $V$. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Parameterisations" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Parameterisation by Finnis et al. 1984, _A simple empirical N-body potential for transition metals_" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "| element | d | A | $\\beta$ | c | $c_0$ | $c_1$ | $c_2$ | \n", + "| --- | --- | --- | --- | --- | --- | --- | --- |\n", + "| V | 3.692767 | 2.010637 | 0 | 3.8 | -0.8816318 | 1.4907756 | -0.3976370 |\n", + "| Nb | 3.915354 | 3.013789 | 0 | 4.2 | -1.5640104 | 2.0055779 | -0.4663764 |\n", + "| Ta | 4.076980 | 2.591061 | 0 | 4.2 | 1.2157373 | 0.0271471 | -0.1217350 |\n", + "| Cr | 3.915720 | 1.453418 | 1.8 | 2.9 | 29.1429813 | -23.3975027 | 4.7578297 |\n", + "| Mo | 4.114825 | 1.887117 | 0 | 3.25 | 43.4475218 | -31.9332978 | 6.0804249 |\n", + "| W | 4.400224 | 1.896373 | 0 | 3.25 | 47.1346499 | -33.7665655 | 6.2541999 |\n", + "| Fe | 3.699579 | 1.889846 | 1.8 | 3.4 | 1.2110601 | -0.7510840 | 0.1380773 |" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note that this potential assumes *metal* units, e.g. eV, K and so on: https://lammps.sandia.gov/doc/units.html" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Collecting the parameters" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# kb = Molly.kb \n", + "fs84, elements, masses, bcc_lattice_constants, reference_energies = Molly.get_finnissinclair1984(true)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Glue contribution - $\\phi(r)$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To compute the glue density $\\rho$ for an atom, we need to compute the contributions from its neighbours (self-contribution is ignored). The neighbour contributions are obtained using $\\phi(r)$:\n", + "\n", + "$$\n", + "\\phi(r) = (r-d)^2 + \\beta (r-d)^3/d\n", + "$$\n", + "\n", + "$d$ and $\\beta$ are model parameters." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "el = \"V\"\n", + "el_pair = string(el,el)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "d, β = fs84.singles[el].d, fs84.singles[el].β" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "d, A, β = fs84.singles[el].d, fs84.singles[el].A, fs84.singles[el].β\n", + "c, c₀, c₁, c₂ = fs84.pairs[el_pair].c, fs84.pairs[el_pair].c₀, fs84.pairs[el_pair].c₁, fs84.pairs[el_pair].c₂" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "r = collect(range(0, stop=6, length=1000));\n", + "ɸ = Molly.glue.(r, β, d);" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ɸs = [ɸ]\n", + "element_pairs = [el_pair]\n", + "for el in elements[2:end]\n", + " el_pair = string(el,el)\n", + " d, β = fs84.singles[el].d, fs84.singles[el].β\n", + " ɸ = Molly.glue.(r, β, d)\n", + " append!(ɸs,[ɸ])\n", + " element_pairs = hcat(element_pairs, el_pair)\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "element_pairs" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plot(r, ɸs, label=element_pairs, xlabel=\"r\", ylabel=\"phi(r)\", title=\"Glue contributions\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Interestingly we find that some $\\phi$ do indeed take on negative values. This is interesting because of the choice for $f$ by Finnis and Sinclair." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Glue energy - $f(\\rho)$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Computing an energy based on local glue values" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "u_\\text{glue} = -A \\cdot \\sqrt{\\rho}\n", + "$$\n", + "\n", + "$$\n", + "\\rho = \\sum_{j \\in \\text{neighborhood}(i)} \\phi(r_{ij})\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "From a physical point of view it makes sense to expect $\\rho > 0$, but in some situations we still might obtain $\\rho < 0$, and have a problem, unless we allow for complex energy values.\n", + "\n", + "Let's plot $u_{\\text{glue}}$ for $\\rho \\in \\mathbb{R}^+$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "el = \"V\"\n", + "ρ = collect(range(0, stop=50, length=100)) \n", + "uₙ = Molly.Uglue.(ρ, fs84.singles[el].A)\n", + "\n", + "uₙs = [uₙ]\n", + "els = [el]\n", + "for el in elements[2:end]\n", + " uₙ = Molly.Uglue.(ρ, fs84.singles[el].A)\n", + " append!(uₙs,[uₙ])\n", + " els = hcat(els,el)\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plot(ρ, uₙs, label=els, xlabel=\"Glue density\", ylabel=\"Glue energy\", title=\"Glue energy varying with glue density\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Pair energy - $V(r)$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The pair energy is defined as\n", + "$$\n", + "V_{ij}(r_{ij}) = \n", + "\\begin{cases} \n", + "r \\le c, & (r-c)^2 \\left( c_0 + c_1 r + c_2 r^2 \\right) \\\\\n", + "r > c, & 0 \\\\\n", + "\\end{cases}\n", + "$$\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's have a look" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "d, A, β = fs84.singles[el].d, fs84.singles[el].A, fs84.singles[el].β\n", + "c, c₀, c₁, c₂ = fs84.pairs[el_pair].c, fs84.pairs[el_pair].c₀, fs84.pairs[el_pair].c₁, fs84.pairs[el_pair].c₂" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "el = \"V\"\n", + "el_pair = string(el,el)\n", + "c, c₀, c₁, c₂ = fs84.pairs[el_pair].c, fs84.pairs[el_pair].c₀, fs84.pairs[el_pair].c₁, fs84.pairs[el_pair].c₂\n", + "\n", + "V = Molly.Upair.(r, c, c₀, c₁, c₂)\n", + "\n", + "Vs = [V]\n", + "element_pairs = [el_pair]\n", + "for el in elements[2:end]\n", + " el_pair = string(el,el)\n", + " c, c₀, c₁, c₂ = fs84.pairs[el_pair].c, fs84.pairs[el_pair].c₀, fs84.pairs[el_pair].c₁, fs84.pairs[el_pair].c₂\n", + " V = Molly.Upair.(r, c, c₀, c₁, c₂)\n", + " append!(Vs,[V])\n", + " element_pairs = hcat(element_pairs, string(el_pair))\n", + "end\n", + "\n", + "plot(r, Vs, label=element_pairs, xlabel=\"r\", ylabel=\"Pair energy contribution\", title=\"Pair energy vs separation\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Simulating" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now that we have looked at the individual components of the model, let's simulate Tungsten near room temperature for a $3\\times3\\times3$ supercell of a body centered cubic crystal." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Creating the crystal" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "?Crystal.make_bcc_unitcell" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "nx = 3\n", + "ny = 3\n", + "nz = 3\n", + "el = \"W\"\n", + "\n", + "a = bcc_lattice_constants[el]\n", + "el2atom_map = Dict(el=>Atom(name=el, mass=masses[el]))\n", + "elements = [el for _ in 1:2]\n", + "atoms, coords, box, box_size, box_vectors = Crystal.make_bcc_unitcell(elements, a=a, el2atom_map=el2atom_map)\n", + "sc_atoms, sc_coords, sc_box, sc_box_size = Crystal.make_supercell(atoms, coords, box, box_size, nx=nx, ny=ny, nz=nz)\n", + "n_atoms = length(sc_atoms)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Setting simulation specs" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "T = 300. # Kelvin\n", + "T *= fs84.kb\n", + "n_steps = 5000\n", + "dt = .002 # ps; ns = 1e-9, ps = 1e-12, fs = 1e-15\n", + "\n", + "general_inters = (fs84,)\n", + "velocities = [velocity(atom.mass, T, dims=3) for atom in sc_atoms]\n", + "nb_matrix = trues(n_atoms,n_atoms)\n", + "dist_cutoff = 2 * a\n", + "nf = DistanceNeighbourFinder(nb_matrix, 1, dist_cutoff)\n", + "# thermostat = NoThermostat()\n", + "thermostat = AndersenThermostat(1.)\n", + "\n", + "loggers = Dict(\n", + " \"temperature\" => TemperatureLogger(1),\n", + " \"pot\" => EnergyLogger(1),\n", + " \"coords\" => CoordinateLogger(100),\n", + " \"velos\" => VelocityLogger(1),\n", + " \"forces\" => ForceLogger(100),\n", + " \"glue\" => GlueDensityLogger(1),\n", + "# \"writer\" => StructureWriter(5,\"traj_fs_bcc_vac_fe.pdb\")\n", + ")\n", + "\n", + "s = Simulation(\n", + " simulator=VelocityVerlet(), \n", + " atoms=sc_atoms, \n", + " general_inters=general_inters,\n", + " coords=[SVector{3}(v) for v in sc_coords], \n", + " velocities=velocities,\n", + " temperature=T, \n", + " box_size=sc_box_size[1,1],\n", + " timestep=dt,\n", + " n_steps=n_steps,\n", + " neighbour_finder=nf,\n", + " loggers=loggers,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Running the simulation (this should take about 30s)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "simulate!(s)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Visualizing content of the loggers" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "x = collect(1:length(s.loggers[\"glue\"].glue_densities))\n", + "y = [[arr[i] for arr in s.loggers[\"glue\"].glue_densities] for i in 1:n_atoms]\n", + "plot(x, y, title=\"glue densities\", legend=false)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "x = collect(1:length(s.loggers[\"temperature\"].temperatures))\n", + "y = s.loggers[\"temperature\"].temperatures/fs84.kb\n", + "plot(x,y,title=\"Temperature\",legend=false)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "x = collect(1:length(s.loggers[\"pot\"].energies))\n", + "y = s.loggers[\"pot\"].energies / length(s.coords)\n", + "plot(x,y,title=\"Potential energy per atom\",legend=false)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 1.5.3", + "language": "julia", + "name": "julia-1.5" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.5.3" + }, + "toc": { + "base_numbering": 1, + "nav_menu": { + "height": "167px", + "width": "309px" + }, + "number_sections": true, + "sideBar": true, + "skip_h1_title": false, + "title_cell": "Table of Contents", + "title_sidebar": "Contents", + "toc_cell": false, + "toc_position": {}, + "toc_section_display": true, + "toc_window_display": false + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} From 331cc70404d97cca736121b31afd4467df9dc7b5 Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Sat, 27 Mar 2021 21:40:31 +0100 Subject: [PATCH 35/49] removed because relocated to ./docs --- fs.ipynb | 572 ------------------------------------------------------- 1 file changed, 572 deletions(-) delete mode 100644 fs.ipynb diff --git a/fs.ipynb b/fs.ipynb deleted file mode 100644 index e19a90cb..00000000 --- a/fs.ipynb +++ /dev/null @@ -1,572 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# The Finnis-Sinclair potential" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Reading:\n", - "- Finnis and Sinclair 1984, [A simple empirical N-body potential for transition metals](https://doi.org/10.1080/01418618408244210)\n", - "- Daw and Baskes 1984, [Embedded-atom method: Derivation and application to impurities, surfaces, and other defects in metals](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.29.6443)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Content:\n", - "1. the math of a glue potential in general & what's special in the Finnis-Sinclair version\n", - "2. show components of the model (pair energy, glue density, glue energy)\n", - "3. show how to run a simulation" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import Pkg\n", - "\n", - "Pkg.activate(\".\")\n", - "\n", - "# Pkg.add(url=\"https://github.com/eschmidt42/crystal\")\n", - "\n", - "using Molly\n", - "using DataFrames\n", - "using Plots\n", - "using Test\n", - "using Crystal" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## The math 💖" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The potential energy is composed out of a pair and a \"glue\" contribution\n", - "$$\n", - "U = U_{\\text{pair}} + U_{\\text{glue}}\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The pair contribution\n", - "$$\n", - "U_{\\text{pair}} = \\sum_{i,j>i}V_{ij}(r_{ij})\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "where $r_{ij} = \\left\\| \\vec{R}_j - \\vec{R}_i| \\right\\|_2$." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The energy contribution of the glue is taken for each atom\n", - "$$\n", - "U_{\\text{glue}} = \\sum_i f(\\rho_i)\n", - "$$\n", - "\n", - "where $\\rho_i = \\sum_j \\phi(r_{ij})$ is the glue at atom $i$." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Combining both\n", - "$$\n", - "\\begin{split}\n", - "U &= U_{\\text{pair}} + U_{\\text{glue}} \\\\\n", - " &= \\sum_{i,j>i}V_{ij}(r_{ij}) + \\sum_i f(\\rho_i) \\\\\n", - "\\end{split}\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "To compute the force on atom i with this potential we apply\n", - "$$\n", - "\\vec{F}_i = - \\partial_i U\n", - "$$\n", - "\n", - "where $\\partial_i$ is the wiggling of $\\vec{R}_i$." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Applying the chain rule we obtain\n", - "$$\n", - "\\begin{split}\n", - "\\vec{F}_i &= - \\sum_j V^\\prime_{ij}(r_{ij}) \\left( -\\frac{\\vec{R}_{ij}}{r_{ij}} \\right) + f^\\prime_i(\\rho_i) \\sum_j \\phi_{j}^\\prime(r_{ij}) \\left( -\\frac{\\vec{R}_{ij}}{r_{ij}} \\right) + \\sum_j f^\\prime_j(\\rho_j) \\phi_i(r_{ji}) \\left( \\frac{\\vec{R}_{ji}}{r_{ji}} \\right) \\\\\n", - " &= \\sum_j \\left[ V^\\prime_{ij}(r_{ij}) + f^\\prime_i(\\rho_i) \\phi_{j}^\\prime(r_{ij}) + f^\\prime_j(\\rho_j) \\phi_i(r_{ji}) \\right] \\frac{\\vec{R}_{ij}}{r_{ij}}\n", - "\\end{split}\n", - "$$\n", - "\n", - "The interesting part is the second equation, because this we can directly use to implement the force updates, looping the neighbouring pairs in parallel ($j$ is a neighbouring atom of atom $i$). " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## The model" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Let's look into $f$, $\\phi$ and $V$. " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Parameterisations" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Parameterisation by Finnis et al. 1984, _A simple empirical N-body potential for transition metals_" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "| element | d | A | $\\beta$ | c | $c_0$ | $c_1$ | $c_2$ | \n", - "| --- | --- | --- | --- | --- | --- | --- | --- |\n", - "| V | 3.692767 | 2.010637 | 0 | 3.8 | -0.8816318 | 1.4907756 | -0.3976370 |\n", - "| Nb | 3.915354 | 3.013789 | 0 | 4.2 | -1.5640104 | 2.0055779 | -0.4663764 |\n", - "| Ta | 4.076980 | 2.591061 | 0 | 4.2 | 1.2157373 | 0.0271471 | -0.1217350 |\n", - "| Cr | 3.915720 | 1.453418 | 1.8 | 2.9 | 29.1429813 | -23.3975027 | 4.7578297 |\n", - "| Mo | 4.114825 | 1.887117 | 0 | 3.25 | 43.4475218 | -31.9332978 | 6.0804249 |\n", - "| W | 4.400224 | 1.896373 | 0 | 3.25 | 47.1346499 | -33.7665655 | 6.2541999 |\n", - "| Fe | 3.699579 | 1.889846 | 1.8 | 3.4 | 1.2110601 | -0.7510840 | 0.1380773 |" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Note that this potential assumes *metal* units, e.g. eV, K and so on: https://lammps.sandia.gov/doc/units.html" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Collecting the parameters" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "kb = Molly.kb \n", - "fs84, elements, masses, bcc_lattice_constants, reference_energies = Molly.get_finnissinclair1984(true)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Glue contribution - $\\phi(r)$" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "To compute the glue density $\\rho$ for an atom, we need to compute the contributions from its neighbours (self-contribution is ignored). The neighbour contributions are obtained using $\\phi(r)$:\n", - "\n", - "$$\n", - "\\phi(r) = (r-d)^2 + \\beta (r-d)^3/d\n", - "$$\n", - "\n", - "$d$ and $\\beta$ are model parameters." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "element_pair, d, A, β, c, c₀, c₁, c₂ = fs84.params[1,:] # parameters for Vanadium" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "r = collect(range(0, stop=6, length=1000));\n", - "ɸ = Molly.glue.(r, β, d);" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "ɸs = [ɸ]\n", - "element_pairs = [element_pair]\n", - "for i in 2:nrow(fs84.params)\n", - " element_pair, d, A, β, c, c₀, c₁, c₂ = fs84.params[i,:]\n", - " ɸ = Molly.glue.(r, β, d)\n", - " append!(ɸs,[ɸ])\n", - " element_pairs = hcat(element_pairs, string(element_pair))\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plot(r, ɸs, label=element_pairs, xlabel=\"r\", ylabel=\"phi(r)\", title=\"Glue contributions\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Interestingly we find that some $\\phi$ do indeed take on negative values. This is interesting because of the choice for $f$ by Finnis and Sinclair." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Glue energy - $f(\\rho)$" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Computing an energy based on local glue values" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "$$\n", - "u_\\text{glue} = -A \\cdot \\sqrt{\\rho}\n", - "$$\n", - "\n", - "$$\n", - "\\rho = \\sum_{j \\in \\text{neighborhood}(i)} \\phi(r_{ij})\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "From a physical point of view it makes sense to expect $\\rho > 0$, but in some situations we still might obtain $\\rho < 0$, and have a problem, unless we allow for complex energy values.\n", - "\n", - "Let's plot $u_{\\text{glue}}$ for $\\rho \\in \\mathbb{R}^+$." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "ρ = collect(range(0, stop=50, length=100))\n", - "A = fs84.params.A[1] # Va\n", - "uₙ = Molly.Uglue.(ρ, A)\n", - "element_pair = fs84.params.element_pair[1];\n", - "\n", - "uₙs = [uₙ]\n", - "element_pairs = [element_pair]\n", - "for i in 2:nrow(fs84.params)\n", - " element_pair, d, A, β, c, c₀, c₁, c₂ = fs84.params[i,:]\n", - " uₙ = Molly.Uglue.(ρ, A)\n", - " append!(uₙs,[uₙ])\n", - " element_pairs = hcat(element_pairs, string(element_pair))\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plot(ρ, uₙs, label=element_pairs, xlabel=\"Glue density\", ylabel=\"Glue energy\", title=\"Glue energy varying with glue density\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Pair energy - $V(r)$" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The pair energy is defined as\n", - "$$\n", - "V_{ij}(r_{ij}) = \n", - "\\begin{cases} \n", - "r \\le c, & (r-c)^2 \\left( c_0 + c_1 r + c_2 r^2 \\right) \\\\\n", - "r > c, & 0 \\\\\n", - "\\end{cases}\n", - "$$\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Let's have a look" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "element_pair, d, A, β, c, c₀, c₁, c₂ = fs84.params[1,:] # parameters for Vanadium\n", - "\n", - "V = Molly.Upair.(r, c, c₀, c₁, c₂)\n", - "\n", - "Vs = [V]\n", - "element_pairs = [element_pair]\n", - "for i in 2:nrow(fs84.params)\n", - " element_pair, d, A, β, c, c₀, c₁, c₂ = fs84.params[i,:]\n", - " V = Molly.Upair.(r, c, c₀, c₁, c₂)\n", - " append!(Vs,[V])\n", - " element_pairs = hcat(element_pairs, string(element_pair))\n", - "end\n", - "\n", - "plot(r, Vs, label=element_pairs, xlabel=\"r\", ylabel=\"Pair energy contribution\", title=\"Pair energy vs separation\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Simulating" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Now that we have looked at the individual components of the model, let's simulate Tungsten near room temperature for a $3\\times3\\times3$ supercell of a body centered cubic crystal." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Creating the crystal" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "nx = 3\n", - "ny = 3\n", - "nz = 3\n", - "element = \"W\"\n", - "\n", - "a = bcc_lattice_constants[element]\n", - "atoms, coords, box, box_size, box_vectors = Crystal.make_bcc_unitcell(element, a=a)\n", - "sc_atoms, sc_coords, sc_box, sc_box_size = Crystal.make_supercell(atoms, coords, box, box_size, nx=nx, ny=ny, nz=nz)\n", - "n_atoms = length(sc_atoms)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Setting simulation specs" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "T = 300. # Kelvin\n", - "T *= kb\n", - "n_steps = 5000\n", - "dt = .002 # ps; ns = 1e-9, ps = 1e-12, fs = 1e-15\n", - "\n", - "general_inters = (fs84,)\n", - "velocities = [velocity(sc_atoms[i].mass, T, dims=3) for i in 1:n_atoms]\n", - "nb_matrix = trues(n_atoms,n_atoms)\n", - "dist_cutoff = 2 * a\n", - "nf = DistanceNeighbourFinder(nb_matrix, 1, dist_cutoff)\n", - "# thermostat = NoThermostat()\n", - "thermostat = AndersenThermostat(1.)\n", - "\n", - "loggers = Dict(\n", - " \"temperature\" => TemperatureLogger(1),\n", - " \"pot\" => EnergyLogger(1),\n", - " \"coords\" => CoordinateLogger(100),\n", - " \"velos\" => VelocityLogger(1),\n", - " \"forces\" => ForceLogger(100),\n", - " \"glue\" => GlueDensityLogger(1),\n", - "# \"writer\" => StructureWriter(5,\"traj_fs_bcc_vac_fe.pdb\")\n", - ")\n", - "\n", - "s = Simulation(\n", - " simulator=VelocityVerlet(), \n", - " atoms=sc_atoms, \n", - " general_inters=general_inters,\n", - " coords=[SVector{3}(v) for v in sc_coords], \n", - " velocities=velocities,\n", - " temperature=T, \n", - " box_size=sc_box_size[1,1],\n", - " timestep=dt,\n", - " n_steps=n_steps,\n", - " neighbour_finder=nf,\n", - " loggers=loggers,\n", - ")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Running the simulation (this should take about 30s)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "simulate!(s)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Visualizing content of the loggers" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "x = collect(1:length(s.loggers[\"glue\"].glue_densities))\n", - "y = [[arr[i] for arr in s.loggers[\"glue\"].glue_densities] for i in 1:n_atoms]\n", - "plot(x, y, title=\"glue densities\", legend=false)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "x = collect(1:length(s.loggers[\"temperature\"].temperatures))\n", - "y = s.loggers[\"temperature\"].temperatures/kb\n", - "plot(x,y,title=\"Temperature\",legend=false)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "x = collect(1:length(s.loggers[\"pot\"].energies))\n", - "y = s.loggers[\"pot\"].energies / length(s.coords)\n", - "plot(x,y,title=\"Potential energy per atom\",legend=false)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Julia 1.5.3", - "language": "julia", - "name": "julia-1.5" - }, - "language_info": { - "file_extension": ".jl", - "mimetype": "application/julia", - "name": "julia", - "version": "1.5.3" - }, - "toc": { - "base_numbering": 1, - "nav_menu": { - "height": "167px", - "width": "309px" - }, - "number_sections": true, - "sideBar": true, - "skip_h1_title": false, - "title_cell": "Table of Contents", - "title_sidebar": "Contents", - "toc_cell": false, - "toc_position": {}, - "toc_section_display": true, - "toc_window_display": false - } - }, - "nbformat": 4, - "nbformat_minor": 4 -} From d87b4c6939ac885d339cdc8c6881e5041afabeb8 Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Sat, 27 Mar 2021 21:43:23 +0100 Subject: [PATCH 36/49] removed superfluous cells --- docs/fs.ipynb | 41 ----------------------------------------- 1 file changed, 41 deletions(-) diff --git a/docs/fs.ipynb b/docs/fs.ipynb index 976a8fdf..1899772a 100644 --- a/docs/fs.ipynb +++ b/docs/fs.ipynb @@ -203,7 +203,6 @@ "metadata": {}, "outputs": [], "source": [ - "# kb = Molly.kb \n", "fs84, elements, masses, bcc_lattice_constants, reference_energies = Molly.get_finnissinclair1984(true)" ] }, @@ -246,16 +245,6 @@ "d, β = fs84.singles[el].d, fs84.singles[el].β" ] }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "d, A, β = fs84.singles[el].d, fs84.singles[el].A, fs84.singles[el].β\n", - "c, c₀, c₁, c₂ = fs84.pairs[el_pair].c, fs84.pairs[el_pair].c₀, fs84.pairs[el_pair].c₁, fs84.pairs[el_pair].c₂" - ] - }, { "cell_type": "code", "execution_count": null, @@ -283,15 +272,6 @@ "end" ] }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "element_pairs" - ] - }, { "cell_type": "code", "execution_count": null, @@ -400,18 +380,6 @@ "Let's have a look" ] }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "d, A, β = fs84.singles[el].d, fs84.singles[el].A, fs84.singles[el].β\n", - "c, c₀, c₁, c₂ = fs84.pairs[el_pair].c, fs84.pairs[el_pair].c₀, fs84.pairs[el_pair].c₁, fs84.pairs[el_pair].c₂" - ] - }, { "cell_type": "code", "execution_count": null, @@ -458,15 +426,6 @@ "Creating the crystal" ] }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "?Crystal.make_bcc_unitcell" - ] - }, { "cell_type": "code", "execution_count": null, From bbc47f10d3fe29eda7173ff091b942166bba3913 Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Fri, 2 Apr 2021 13:25:55 +0200 Subject: [PATCH 37/49] adjusted glue tests to changes in Crystal.jl and introduced more named testsets for legitibility of the test summary in case of failures --- test/glue.jl | 64 ++++++++++++++++++++++++++++++---------------------- 1 file changed, 37 insertions(+), 27 deletions(-) diff --git a/test/glue.jl b/test/glue.jl index 5edca2f0..d3eff066 100644 --- a/test/glue.jl +++ b/test/glue.jl @@ -2,8 +2,6 @@ using Molly using Test using Crystal -fs_inter, elements, masses, bcc_lattice_constants, reference_energies = Molly.get_finnissinclair1984(true) - function forces_are_zero(forces; dims=3) return all([isapprox(f, zero(rand(3)), atol=1e-4) for f in forces]) end @@ -12,13 +10,11 @@ function glue_remains_reasonable(glue_init,glue_end) return all(isapprox.(glue_init,glue_end,rtol=.1)) end -function groundstate_energy_as_expected(s, element, inter) - ref = reference_energies[element] +function groundstate_energy_as_expected(s, ref) return isapprox(s.loggers["energy"].energies[1]/length(s.coords), -ref["u"], atol=1e-2) end -function vacancy_formation_energy_as_expected(s, s_vac, element, inter) - ref = reference_energies[element] +function vacancy_formation_energy_as_expected(s, s_vac, ref) u_gs = s.loggers["energy"].energies[1] n_atoms = length(s.coords) u_vac = s_vac.loggers["energy"].energies[1] @@ -27,28 +23,27 @@ function vacancy_formation_energy_as_expected(s, s_vac, element, inter) return isapprox(u_vac_form, ref["u_vac"], atol=7e-2) end -function run_bcc(element::String, fs_inter, T::Real=.01; +function run_bcc(element::String, fs_inter, a::Real; T::Real=.01, nx::Integer=3, ny::Integer=3, nz::Integer=3, vac::Bool=false) masses = Dict("V" => 50.9415, "Nb" => 92.9064, "Ta" => 180.9479, "Cr" => 51.996, "Mo" => 95.94, "W" => 183.85, "Fe" => 55.847) - a = bcc_lattice_constants[element] + elements = [element for _ in 1:2] el2atom_map = Dict(element => Atom(name=element, mass=masses[element])) - atoms, coords, box, box_size, box_vectors = Crystal.make_bcc_unitcell(elements, a=a, el2atom_map=el2atom_map) - sc_atoms, sc_coords, sc_box, sc_box_size = Crystal.make_supercell(atoms, coords, - box, box_size, nx=nx, ny=ny, nz=nz) + unitcell = Crystal.make_bcc_unitcell(elements, a, el2atom_map) + supercell = Crystal.make_supercell(unitcell, nx=nx, ny=ny, nz=nz) if vac # introducing a vacancy - sc_atoms, sc_coords = Crystal.add_vacancies(sc_atoms, sc_coords, ixs=[1]) + supercell = Crystal.add_vacancies(supercell, ixs=[1]) end - n_atoms = length(sc_atoms) + n_atoms = length(supercell.coords) - velocities = [velocity(sc_atoms[i].mass, T*fs_inter.kb, dims=3) for i in 1:n_atoms] + velocities = [velocity(atom.mass, T, dims=3) for atom in supercell.atoms] nb_matrix = trues(n_atoms,n_atoms) - dist_cutoff = 2 * a + dist_cutoff = 2. * a loggers = Dict( "glue" => GlueDensityLogger(1), @@ -60,12 +55,12 @@ function run_bcc(element::String, fs_inter, T::Real=.01; s = Simulation( simulator=VelocityVerlet(), - atoms=sc_atoms, + atoms=supercell.atoms, general_inters=(fs_inter,), - coords=[SVector{3}(v) for v in sc_coords], + coords=[SVector{3}(v) for v in supercell.coords], velocities=velocities, - temperature=T*fs_inter.kb, - box_size=sc_box_size[1,1], + temperature=T, + box_size=supercell.edge_lengths[1], timestep=.002, n_steps=100, neighbour_finder=nf, @@ -77,17 +72,32 @@ end @testset "Glue: Finnis-Sinclair" begin T = .01 # K + fs_inter, elements, masses, bcc_lattice_constants, reference_energies = Molly.get_finnissinclair1984(true) + + T *= fs_inter.kb for element in elements @testset "$element" begin - s = run_bcc(element, fs_inter, T) - @test forces_are_zero(s.loggers["forces"].forces[1]) - @test glue_remains_reasonable(s.loggers["glue"].glue_densities[1], - s.loggers["glue"].glue_densities[end]) - @test groundstate_energy_as_expected(s, element, fs_inter) + a = bcc_lattice_constants[element] + + s = run_bcc(element, fs_inter, a, T=T, vac=false) + @testset "forces zero" begin + @test forces_are_zero(s.loggers["forces"].forces[1]) + end + @testset "reasonable glue" begin + @test glue_remains_reasonable(s.loggers["glue"].glue_densities[1], + s.loggers["glue"].glue_densities[end]) + end + @testset "checking groundstate" begin + @test groundstate_energy_as_expected(s, reference_energies[element]) + end - s_vac = run_bcc(element, fs_inter, T; vac=true) - @test !forces_are_zero(s_vac.loggers["forces"].forces[1]) - @test vacancy_formation_energy_as_expected(s, s_vac, element, fs_inter) + s_vac = run_bcc(element, fs_inter, a, T=T, vac=true) + @testset "vacancy forces nonzero" begin + @test !forces_are_zero(s_vac.loggers["forces"].forces[1]) + end + @testset "vacancy formation energy" begin + @test vacancy_formation_energy_as_expected(s, s_vac, reference_energies[element]) + end end end From d073ef8527d830cdc0fb8bcbb57fc9eb92d7b53d Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Mon, 10 May 2021 18:43:21 +0200 Subject: [PATCH 38/49] updated test: switched from Crystal to SingleCrystal package --- test/glue.jl | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/test/glue.jl b/test/glue.jl index d3eff066..68c4d252 100644 --- a/test/glue.jl +++ b/test/glue.jl @@ -1,6 +1,6 @@ using Molly using Test -using Crystal +using SingleCrystal function forces_are_zero(forces; dims=3) return all([isapprox(f, zero(rand(3)), atol=1e-4) for f in forces]) @@ -23,23 +23,23 @@ function vacancy_formation_energy_as_expected(s, s_vac, ref) return isapprox(u_vac_form, ref["u_vac"], atol=7e-2) end +make_atom(name,mass) = Atom(name=name,mass=mass) + function run_bcc(element::String, fs_inter, a::Real; T::Real=.01, nx::Integer=3, ny::Integer=3, nz::Integer=3, vac::Bool=false) masses = Dict("V" => 50.9415, "Nb" => 92.9064, "Ta" => 180.9479, "Cr" => 51.996, "Mo" => 95.94, "W" => 183.85, "Fe" => 55.847) - - elements = [element for _ in 1:2] - el2atom_map = Dict(element => Atom(name=element, mass=masses[element])) - unitcell = Crystal.make_bcc_unitcell(elements, a, el2atom_map) - supercell = Crystal.make_supercell(unitcell, nx=nx, ny=ny, nz=nz) + + unitcell = SingleCrystal.make_bcc_unitcell(element, a, make_atom) + supercell = SingleCrystal.make_supercell(unitcell, nx=nx, ny=ny, nz=nz) if vac # introducing a vacancy - supercell = Crystal.add_vacancies(supercell, ixs=[1]) + supercell = SingleCrystal.add_vacancies(supercell, ixs=[1]) end - n_atoms = length(supercell.coords) + n_atoms = length(supercell.positions) velocities = [velocity(atom.mass, T, dims=3) for atom in supercell.atoms] nb_matrix = trues(n_atoms,n_atoms) @@ -57,7 +57,7 @@ function run_bcc(element::String, fs_inter, a::Real; T::Real=.01, simulator=VelocityVerlet(), atoms=supercell.atoms, general_inters=(fs_inter,), - coords=[SVector{3}(v) for v in supercell.coords], + coords=[SVector{3}(v) for v in supercell.positions], velocities=velocities, temperature=T, box_size=supercell.edge_lengths[1], @@ -86,7 +86,7 @@ end @testset "reasonable glue" begin @test glue_remains_reasonable(s.loggers["glue"].glue_densities[1], s.loggers["glue"].glue_densities[end]) - end + end @testset "checking groundstate" begin @test groundstate_energy_as_expected(s, reference_energies[element]) end From fbf262928108aa7842934dd15a86c3f36e4e9f21 Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Mon, 10 May 2021 18:43:41 +0200 Subject: [PATCH 39/49] switched toml from Crystal to SingleCrystal --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index dec14060..881409cc 100644 --- a/Project.toml +++ b/Project.toml @@ -7,7 +7,6 @@ version = "0.2.2" BioStructures = "de9282ab-8554-53be-b2d6-f6c222edabfc" CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" -Crystal = "3c6eccdf-2a89-4c24-a1d4-ff210daa8476" Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" @@ -17,6 +16,7 @@ NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" Requires = "ae029012-a4dd-5104-9daa-d747884805df" +SingleCrystal = "3c6eccdf-2a89-4c24-a1d4-ff210daa8476" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" From d95cfbe260bd848ac818b585a4e767c654839f61 Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Mon, 10 May 2021 19:06:02 +0200 Subject: [PATCH 40/49] switched to toml SingleCrystal.jl --- docs/Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/Project.toml b/docs/Project.toml index 8a325b7e..6befc670 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,9 +1,9 @@ [deps] -Crystal = "3c6eccdf-2a89-4c24-a1d4-ff210daa8476" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" Molly = "aa0f7f06-fcc0-5ec4-a7f3-a573f33f9c4c" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +SingleCrystal = "3c6eccdf-2a89-4c24-a1d4-ff210daa8476" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] From b276156041a34eb6c273795bc4b8ab2e2972baab Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Mon, 10 May 2021 19:07:14 +0200 Subject: [PATCH 41/49] updated notebook to work with the new data structure of the potential and changes in the package constructing the bcc crystal --- docs/fs.ipynb | 132 ++++++++++++++++---------------------------------- 1 file changed, 43 insertions(+), 89 deletions(-) diff --git a/docs/fs.ipynb b/docs/fs.ipynb index 1899772a..4ac5ded5 100644 --- a/docs/fs.ipynb +++ b/docs/fs.ipynb @@ -22,8 +22,8 @@ "source": [ "Content:\n", "1. the math of a glue potential in general & what's special in the Finnis-Sinclair version\n", - "2. show components of the model (pair energy, glue density, glue energy)\n", - "3. show how to run a simulation" + "2. components of the Finnis-Sinclair model (pair energy, glue density, glue energy)\n", + "3. running a simulation and visualizing properties" ] }, { @@ -37,9 +37,7 @@ "import Pkg\n", "\n", "Pkg.activate(\".\")\n", - "Pkg.instantiate()\n", - "\n", - "# Pkg.add(url=\"https://github.com/eschmidt42/crystal\")" + "Pkg.instantiate()" ] }, { @@ -51,7 +49,17 @@ "using Molly\n", "using Plots\n", "using Test\n", - "using Crystal" + "using SingleCrystal" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "make_atom(name,mass) = Atom(name=name,mass=mass)\n", + "SingleCrystal.make_bcc_unitcell(\"W\", 2.3, make_atom)" ] }, { @@ -203,7 +211,7 @@ "metadata": {}, "outputs": [], "source": [ - "fs84, elements, masses, bcc_lattice_constants, reference_energies = Molly.get_finnissinclair1984(true)" + "fs84, elements, masses, bcc_lattice_constants, reference_energies = Molly.get_finnissinclair1984(true);" ] }, { @@ -226,25 +234,6 @@ "$d$ and $\\beta$ are model parameters." ] }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "el = \"V\"\n", - "el_pair = string(el,el)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "d, β = fs84.singles[el].d, fs84.singles[el].β" - ] - }, { "cell_type": "code", "execution_count": null, @@ -252,33 +241,15 @@ "outputs": [], "source": [ "r = collect(range(0, stop=6, length=1000));\n", - "ɸ = Molly.glue.(r, β, d);" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "ɸs = [ɸ]\n", - "element_pairs = [el_pair]\n", - "for el in elements[2:end]\n", - " el_pair = string(el,el)\n", + "\n", + "ɸs = []\n", + "for el in elements\n", " d, β = fs84.singles[el].d, fs84.singles[el].β\n", " ɸ = Molly.glue.(r, β, d)\n", " append!(ɸs,[ɸ])\n", - " element_pairs = hcat(element_pairs, el_pair)\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plot(r, ɸs, label=element_pairs, xlabel=\"r\", ylabel=\"phi(r)\", title=\"Glue contributions\")" + "end\n", + "\n", + "plot(r, ɸs, label=hcat(elements...), xlabel=\"r\", ylabel=\"phi(r)\", title=\"Glue contributions\")" ] }, { @@ -330,26 +301,15 @@ "metadata": {}, "outputs": [], "source": [ - "el = \"V\"\n", "ρ = collect(range(0, stop=50, length=100)) \n", - "uₙ = Molly.Uglue.(ρ, fs84.singles[el].A)\n", "\n", - "uₙs = [uₙ]\n", - "els = [el]\n", - "for el in elements[2:end]\n", + "uₙs = []\n", + "for el in elements\n", " uₙ = Molly.Uglue.(ρ, fs84.singles[el].A)\n", " append!(uₙs,[uₙ])\n", - " els = hcat(els,el)\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plot(ρ, uₙs, label=els, xlabel=\"Glue density\", ylabel=\"Glue energy\", title=\"Glue energy varying with glue density\")" + "end\n", + "\n", + "plot(ρ, uₙs, label=hcat(elements...), xlabel=\"Glue density\", ylabel=\"Glue energy\", title=\"Glue energy varying with glue density\")" ] }, { @@ -386,23 +346,17 @@ "metadata": {}, "outputs": [], "source": [ - "el = \"V\"\n", - "el_pair = string(el,el)\n", - "c, c₀, c₁, c₂ = fs84.pairs[el_pair].c, fs84.pairs[el_pair].c₀, fs84.pairs[el_pair].c₁, fs84.pairs[el_pair].c₂\n", - "\n", - "V = Molly.Upair.(r, c, c₀, c₁, c₂)\n", - "\n", - "Vs = [V]\n", - "element_pairs = [el_pair]\n", - "for el in elements[2:end]\n", + "Vs = []\n", + "element_pairs = []\n", + "for el in elements\n", " el_pair = string(el,el)\n", " c, c₀, c₁, c₂ = fs84.pairs[el_pair].c, fs84.pairs[el_pair].c₀, fs84.pairs[el_pair].c₁, fs84.pairs[el_pair].c₂\n", " V = Molly.Upair.(r, c, c₀, c₁, c₂)\n", " append!(Vs,[V])\n", - " element_pairs = hcat(element_pairs, string(el_pair))\n", + " append!(element_pairs, [el_pair])\n", "end\n", "\n", - "plot(r, Vs, label=element_pairs, xlabel=\"r\", ylabel=\"Pair energy contribution\", title=\"Pair energy vs separation\")" + "plot(r, Vs, label=hcat(element_pairs...), xlabel=\"r\", ylabel=\"Pair energy contribution\", title=\"Pair energy vs separation\")" ] }, { @@ -438,11 +392,9 @@ "el = \"W\"\n", "\n", "a = bcc_lattice_constants[el]\n", - "el2atom_map = Dict(el=>Atom(name=el, mass=masses[el]))\n", - "elements = [el for _ in 1:2]\n", - "atoms, coords, box, box_size, box_vectors = Crystal.make_bcc_unitcell(elements, a=a, el2atom_map=el2atom_map)\n", - "sc_atoms, sc_coords, sc_box, sc_box_size = Crystal.make_supercell(atoms, coords, box, box_size, nx=nx, ny=ny, nz=nz)\n", - "n_atoms = length(sc_atoms)" + "unitcell = SingleCrystal.make_bcc_unitcell(el, a, make_atom)\n", + "supercell = SingleCrystal.make_supercell(unitcell, nx=nx, ny=ny, nz=nz)\n", + "n_atoms = length(supercell.atoms)" ] }, { @@ -464,9 +416,9 @@ "dt = .002 # ps; ns = 1e-9, ps = 1e-12, fs = 1e-15\n", "\n", "general_inters = (fs84,)\n", - "velocities = [velocity(atom.mass, T, dims=3) for atom in sc_atoms]\n", + "velocities = [velocity(atom.mass, T, dims=3) for atom in supercell.atoms]\n", "nb_matrix = trues(n_atoms,n_atoms)\n", - "dist_cutoff = 2 * a\n", + "dist_cutoff = 2. * a\n", "nf = DistanceNeighbourFinder(nb_matrix, 1, dist_cutoff)\n", "# thermostat = NoThermostat()\n", "thermostat = AndersenThermostat(1.)\n", @@ -478,17 +430,17 @@ " \"velos\" => VelocityLogger(1),\n", " \"forces\" => ForceLogger(100),\n", " \"glue\" => GlueDensityLogger(1),\n", - "# \"writer\" => StructureWriter(5,\"traj_fs_bcc_vac_fe.pdb\")\n", + "# \"writer\" => StructureWriter(5,\"traj_bcc.pdb\")\n", ")\n", "\n", "s = Simulation(\n", " simulator=VelocityVerlet(), \n", - " atoms=sc_atoms, \n", + " atoms=supercell.atoms, \n", " general_inters=general_inters,\n", - " coords=[SVector{3}(v) for v in sc_coords], \n", + " coords=[SVector{3}(v) for v in supercell.positions], \n", " velocities=velocities,\n", " temperature=T, \n", - " box_size=sc_box_size[1,1],\n", + " box_size=supercell.edge_lengths[1],\n", " timestep=dt,\n", " n_steps=n_steps,\n", " neighbour_finder=nf,\n", @@ -506,7 +458,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "scrolled": true + }, "outputs": [], "source": [ "simulate!(s)" From 71522fec08b9cd8b63b318e9e15d43af72aaf4a6 Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Mon, 10 May 2021 19:23:52 +0200 Subject: [PATCH 42/49] removed ForwardDiff --- Project.toml | 1 - src/Molly.jl | 1 - 2 files changed, 2 deletions(-) diff --git a/Project.toml b/Project.toml index 881409cc..0458dda1 100644 --- a/Project.toml +++ b/Project.toml @@ -9,7 +9,6 @@ CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" -ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" KernelDensity = "5ab0869b-81aa-558d-bb23-cbf5423bbe9b" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce" diff --git a/src/Molly.jl b/src/Molly.jl index 2f9058f0..864fd051 100644 --- a/src/Molly.jl +++ b/src/Molly.jl @@ -16,7 +16,6 @@ using Requires using Base.Threads using LinearAlgebra using SparseArrays -using ForwardDiff include("types.jl") include("cutoffs.jl") From b9ece8d502e3fdd70f302c32e302bbed9d3d726f Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Mon, 10 May 2021 19:31:07 +0200 Subject: [PATCH 43/49] removed empty spaces and updated doc strings --- src/forces.jl | 1 - src/interactions/glue_fs.jl | 10 +++++----- 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/src/forces.jl b/src/forces.jl index 45d70f49..b46a9e9e 100644 --- a/src/forces.jl +++ b/src/forces.jl @@ -110,7 +110,6 @@ function accelerations(s::Simulation; parallel::Bool=true) forces[i] /= s.atoms[i].mass s.forces[i] = forces[i] end - return forces end diff --git a/src/interactions/glue_fs.jl b/src/interactions/glue_fs.jl index 416aeca5..64bb8745 100644 --- a/src/interactions/glue_fs.jl +++ b/src/interactions/glue_fs.jl @@ -89,7 +89,7 @@ function get_finnissinclair1984(nl_only::Bool) end """ - glue(r,β,d; f=10.) + glue(r,β,d) The core component of the Finnis-Sinclair type GlueInteraction. Used to calculate contribution to the glue value of an atom, based on @@ -100,7 +100,7 @@ function glue(r, β, d) end """ - ∂glue_∂r(r, β, d; f=10.) + ∂glue_∂r(r, β, d) Derivative of the glue density function. """ @@ -125,16 +125,16 @@ Energy derivative given glue density. """ - Upair(r, c, c₀, c₁, c₂; f=10.) + Upair(r, c, c₀, c₁, c₂) -Energy contribution directly from atom pair distances. f is a fudge factor to help translate a Å model to a nm model. +Energy contribution directly from atom pair distances. """ function Upair(r, c, c₀, c₁, c₂) return (r > c) ? 0 : (r - c)^2 * (c₀ + c₁*r + c₂*r^2) end """ - ∂Upair_∂r(r, c, c₀, c₁, c₂, f=10.) + ∂Upair_∂r(r, c, c₀, c₁, c₂) Derivative of the pair potential. """ From d1f113afa77efd6e8f0072ef58fb21a9bfda7971 Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Fri, 14 May 2021 19:48:34 +0200 Subject: [PATCH 44/49] rm SingleCrystal package as direct dependency --- Project.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 0458dda1..58a43245 100644 --- a/Project.toml +++ b/Project.toml @@ -15,7 +15,6 @@ NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" Requires = "ae029012-a4dd-5104-9daa-d747884805df" -SingleCrystal = "3c6eccdf-2a89-4c24-a1d4-ff210daa8476" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" @@ -35,7 +34,8 @@ julia = "1.5" [extras] Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +SingleCrystal = "3c6eccdf-2a89-4c24-a1d4-ff210daa8476" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Statistics", "Test"] +test = ["Statistics", "SingleCrystal", "Test"] From 05d1990f39e9206a8b90991fd9d24b72f8c4fd4a Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Fri, 14 May 2021 20:14:07 +0200 Subject: [PATCH 45/49] replaced abstract types in structs --- src/interactions/glue_fs.jl | 31 ++++++++++++++++++------------- 1 file changed, 18 insertions(+), 13 deletions(-) diff --git a/src/interactions/glue_fs.jl b/src/interactions/glue_fs.jl index 64bb8745..a1191ff5 100644 --- a/src/interactions/glue_fs.jl +++ b/src/interactions/glue_fs.jl @@ -8,34 +8,39 @@ SpecificInteraction. """ abstract type GlueInteraction <: SpecificInteraction end -struct FinnisSinclairPair +struct FinnisSinclairPair{T<:Real} # container for two element interaction -> c c0 c1 c2 - c::Real - c₀::Real - c₁::Real - c₂::Real + c::T + c₀::T + c₁::T + c₂::T end -struct FinnisSinclairSingle +FinnisSinclairPair(c,c₀,c₁,c₂) = FinnisSinclairPair{typeof(c₀)}(c,c₀,c₁,c₂) + +struct FinnisSinclairSingle{T<:Real} # container for single element stuff: glue density, glue energy -> A - A::Real - β::Real - d::Real + A::T + β::T + d::T end +FinnisSinclairSingle(A,β,d) = FinnisSinclairPair{typeof(A)}(A,β,d) + """ FinnisSinclairInteraction(nl_only,pairs,singles,kb) The Finnis-Sinclair interaction. This interaction expects units to be of these https://lammps.sandia.gov/doc/units.html units (eV, Å, K, ps and so on). """ -struct FinnisSinclair <: GlueInteraction +struct FinnisSinclair{T<:Real} <: GlueInteraction nl_only::Bool - pairs::Dict{String,FinnisSinclairPair} - singles::Dict{String,FinnisSinclairSingle} - kb::Real + pairs::Dict{String,FinnisSinclairPair{T}} + singles::Dict{String,FinnisSinclairSingle{T}} + kb::T end +FinnisSinclair(nl_only, pairs, singles, kb) = FinnisSinclair{typeof(kb)}(nl_only, pairs, singles, kb) """ get_finnissinclair1984(nl_only) From 12c5a31da15ff75a9d933c0e85f09abbd57ea41a Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Fri, 14 May 2021 20:21:26 +0200 Subject: [PATCH 46/49] wrapped glue tests in a testset --- test/runtests.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 2526d63c..17ef3538 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -22,7 +22,9 @@ end CUDA.allowscalar(false) # Check that we never do scalar indexing on the GPU -include("glue.jl") +@testset "glue.jl" begin + include("glue.jl") +end @testset "Spatial" begin @test vector1D(4.0, 6.0, 10.0) == 2.0 From 3b9ce5197017d6167ce430f446ee1734d592d3e9 Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Fri, 14 May 2021 21:04:33 +0200 Subject: [PATCH 47/49] welp... --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 45a6f7c2..9c85bea5 100644 --- a/Project.toml +++ b/Project.toml @@ -33,8 +33,8 @@ StaticArrays = "0.12, 1" julia = "1.6" [extras] -Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" SingleCrystal = "3c6eccdf-2a89-4c24-a1d4-ff210daa8476" +Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] From 076ad3dc75fc290330528fd48926302e10edf639 Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Fri, 14 May 2021 21:05:15 +0200 Subject: [PATCH 48/49] adjusted to renaming of Simnulation struct --- src/interactions/glue_fs.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/interactions/glue_fs.jl b/src/interactions/glue_fs.jl index a1191ff5..8e91e934 100644 --- a/src/interactions/glue_fs.jl +++ b/src/interactions/glue_fs.jl @@ -164,7 +164,7 @@ function update_glue_densities!( end # updating glue densities - for (i,j) in s.neighbours + for (i,j) in s.neighbors # collecting parameters el_i = s.atoms[i].name From d2b3375bd231bfbdf2546e7d7beb5503183d76df Mon Sep 17 00:00:00 2001 From: eschmidt42 <11818904+eschmidt42@users.noreply.github.com> Date: Fri, 14 May 2021 21:05:38 +0200 Subject: [PATCH 49/49] updated test to new Simulation --- test/glue.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/glue.jl b/test/glue.jl index 68c4d252..baf919be 100644 --- a/test/glue.jl +++ b/test/glue.jl @@ -51,7 +51,7 @@ function run_bcc(element::String, fs_inter, a::Real; T::Real=.01, "energy" => EnergyLogger(1) ) - nf = DistanceNeighbourFinder(nb_matrix, 1, dist_cutoff) + nf = DistanceNeighborFinder(nb_matrix, 1, dist_cutoff) s = Simulation( simulator=VelocityVerlet(), @@ -63,7 +63,7 @@ function run_bcc(element::String, fs_inter, a::Real; T::Real=.01, box_size=supercell.edge_lengths[1], timestep=.002, n_steps=100, - neighbour_finder=nf, + neighbor_finder=nf, loggers=loggers, ) simulate!(s)