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

Project.toml

@@ -8,17 +8,15 @@ 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"
 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"
+SingleCrystal = "3c6eccdf-2a89-4c24-a1d4-ff210daa8476"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 

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"

fs.ipynb → docs/fs.ipynb

@@ -29,17 +29,26 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "scrolled": true
+   },
    "outputs": [],
    "source": [
     "import Pkg\n",
     "\n",
     "Pkg.activate(\".\")\n",
+    "Pkg.instantiate()\n",
     "\n",
-    "# Pkg.add(url=\"https://github.com/eschmidt42/crystal\")\n",
-    "\n",
+    "# Pkg.add(url=\"https://github.com/eschmidt42/crystal\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
     "using Molly\n",
-    "using DataFrames\n",
     "using Plots\n",
     "using Test\n",
     "using Crystal"
@@ -194,7 +203,6 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "kb = Molly.kb \n",
     "fs84, elements, masses, bcc_lattice_constants, reference_energies = Molly.get_finnissinclair1984(true)"
    ]
   },
@@ -224,7 +232,17 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "element_pair, d, A, β, c, c₀, c₁, c₂ = fs84.params[1,:] # parameters for Vanadium"
+    "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].β"
    ]
   },
   {
@@ -244,12 +262,13 @@
    "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",
+    "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, string(element_pair))\n",
+    "    element_pairs = hcat(element_pairs, el_pair)\n",
     "end"
    ]
   },
@@ -311,18 +330,16 @@
    "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",
+    "el = \"V\"\n",
+    "ρ = collect(range(0, stop=50, length=100)) \n",
+    "uₙ = Molly.Uglue.(ρ, fs84.singles[el].A)\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",
+    "els = [el]\n",
+    "for el in elements[2:end]\n",
+    "    uₙ = Molly.Uglue.(ρ, fs84.singles[el].A)\n",
     "    append!(uₙs,[uₙ])\n",
-    "    element_pairs = hcat(element_pairs, string(element_pair))\n",
+    "    els = hcat(els,el)\n",
     "end"
    ]
   },
@@ -332,7 +349,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plot(ρ, uₙs, label=element_pairs, xlabel=\"Glue density\", ylabel=\"Glue energy\", title=\"Glue energy varying with glue density\")"
+    "plot(ρ, uₙs, label=els, xlabel=\"Glue density\", ylabel=\"Glue energy\", title=\"Glue energy varying with glue density\")"
    ]
   },
   {
@@ -369,17 +386,20 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "element_pair, d, A, β, c, c₀, c₁, c₂ = fs84.params[1,:] # parameters for Vanadium\n",
+    "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 = [element_pair]\n",
-    "for i in 2:nrow(fs84.params)\n",
-    "    element_pair, d, A, β, c, c₀, c₁, c₂ = fs84.params[i,:]\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(element_pair))\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\")"
@@ -415,10 +435,12 @@
     "nx = 3\n",
     "ny = 3\n",
     "nz = 3\n",
-    "element = \"W\"\n",
+    "el = \"W\"\n",
     "\n",
-    "a = bcc_lattice_constants[element]\n",
-    "atoms, coords, box, box_size, box_vectors = Crystal.make_bcc_unitcell(element, a=a)\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)"
    ]
@@ -437,12 +459,12 @@
    "outputs": [],
    "source": [
     "T = 300. # Kelvin\n",
-    "T *= kb\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(sc_atoms[i].mass, T, dims=3) for i in 1:n_atoms]\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",
@@ -515,7 +537,7 @@
    "outputs": [],
    "source": [
     "x = collect(1:length(s.loggers[\"temperature\"].temperatures))\n",
-    "y = s.loggers[\"temperature\"].temperatures/kb\n",
+    "y = s.loggers[\"temperature\"].temperatures/fs84.kb\n",
     "plot(x,y,title=\"Temperature\",legend=false)"
    ]
   },