diff --git a/Project.toml b/Project.toml
index 430d0dcf..d8c00e4d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,9 +7,14 @@ 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"
+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"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
diff --git a/fs.ipynb b/fs.ipynb
new file mode 100644
index 00000000..349f1938
--- /dev/null
+++ b/fs.ipynb
@@ -0,0 +1,581 @@
+{
+ "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"
+   ]
+  },
+  {
+   "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
+}
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")
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)
diff --git a/src/forces.jl b/src/forces.jl
index bcf09e95..e7f80349 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)
@@ -52,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)
@@ -72,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)
@@ -92,12 +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
+#             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
@@ -134,3 +152,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..01e1f4ee
--- /dev/null
+++ b/src/interactions/glue_fs.jl
@@ -0,0 +1,241 @@
+"""
+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. 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
+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
+
+"""
+    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.
+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
+end
+
+"""
+    ∂glue_∂r(r, β, d; f=10.)
+
+Derivative of the glue density function.
+"""
+∂glue_∂r(r, β, d) = ForwardDiff.derivative(r -> glue(r,β,d), r)
+
+"""
+    Uglue(ρ, A)
+
+Energy based on the glue density .
+"""
+function Uglue(ρ, A)
+   return -A * √ρ 
+end
+
+
+"""
+    ∂Uglue_∂ρ(ρ, A)
+
+Energy derivative given glue density.
+"""
+∂Uglue_∂ρ(ρ,A) = ForwardDiff.derivative(ρ -> Uglue(ρ,A), ρ)
+
+"""
+    Upair(r, c, c₀, c₁, c₂; f=10.)
+
+Energy contribution directly from atom pair distances. f is a fudge factor to help translate a Å model to a nm model. 
+"""
+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)
+
+"""
+    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)
+
+Convenience function to update the densities before the forces are computed in serial/parallel.
+"""
+function update_glue_densities!(
+        inter::FinnisSinclair, 
+        coords, 
+        s::Simulation;
+        parallel::Bool=false
+    )
+    n_atoms = length(s.coords)
+    
+    # wiping old glue densities
+    for i in 1:n_atoms
+        s.glue_densities[i] = 0
+    end
+    
+    # updating glue densities
+    for (i,j) in s.neighbours
+
+        # 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) 
+        
+        # computing distance
+        dr = vector(s.coords[i], s.coords[j], s.box_size)
+        r = norm(dr)
+        
+        # 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
+
+@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)
+    
+    # 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)
+
+    # distance i -> j
+    dr = vector(s.coords[i], s.coords[j], s.box_size)
+    r = norm(dr)
+    dr = normalize(dr)
+
+    # pair contribution
+    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)
+
+    # 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 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)
+    end
+        
+    A = get_pair_params(s.atoms[i].name, s.atoms[i].name, inter).A
+    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)
+    r = norm(dr)
+    pij = get_pair_params(s.atoms[i].name, s.atoms[j].name, inter)
+    return Upair(r, pij.c, pij.c₀, pij.c₁, pij.c₂)
+end
diff --git a/src/loggers.jl b/src/loggers.jl
index f418f327..17ce7599 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::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")
+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)
 
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)
diff --git a/test/glue.jl b/test/glue.jl
new file mode 100644
index 00000000..f6e9d5f5
--- /dev/null
+++ b/test/glue.jl
@@ -0,0 +1,92 @@
+using Molly
+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)
+    return all([isapprox(f, zero(rand(3)), atol=1e-4) for f in forces])
+end
+
+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)
+    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
+
+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 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)
+    
+    if vac
+        # introducing a vacancy
+        sc_atoms, sc_coords = Crystal.add_vacancies(sc_atoms, sc_coords, ixs=[1])
+    end
+    n_atoms = length(sc_atoms)
+    
+    velocities = [velocity(1., T, dims=3) for i in 1:n_atoms]
+    nb_matrix = trues(n_atoms,n_atoms)
+    dist_cutoff = 2 * a
+
+    loggers = Dict(
+        "glue" => GlueDensityLogger(1),
+        "forces" => ForceLogger(5),
+        "energy" => EnergyLogger(1)
+    )
+
+    nf = DistanceNeighbourFinder(nb_matrix, 1, dist_cutoff)
+    
+    s = Simulation(
+        simulator=VelocityVerlet(), 
+        atoms=sc_atoms, 
+        general_inters=(fs_inter,),
+        coords=[SVector{3}(v) for v in sc_coords], 
+        velocities=velocities,
+        temperature=T, 
+        box_size=sc_box_size[1,1],
+        timestep=.002,
+        n_steps=100,
+        neighbour_finder=nf,
+        loggers=loggers,
+    )
+    simulate!(s)
+    return s
+end
+
+@testset "Finnis-Sinclair" begin
+    
+    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
+    
+end