diff --git a/Project.toml b/Project.toml
index 02a06440..9c85bea5 100644
--- a/Project.toml
+++ b/Project.toml
@@ -33,8 +33,9 @@ StaticArrays = "0.12, 1"
 julia = "1.6"
 
 [extras]
+SingleCrystal = "3c6eccdf-2a89-4c24-a1d4-ff210daa8476"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Statistics", "Test"]
+test = ["Statistics", "SingleCrystal", "Test"]
diff --git a/docs/Project.toml b/docs/Project.toml
index f3233bec..6befc670 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,6 +1,10 @@
 [deps]
+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]
 Documenter = "0.25"
diff --git a/docs/fs.ipynb b/docs/fs.ipynb
new file mode 100644
index 00000000..4ac5ded5
--- /dev/null
+++ b/docs/fs.ipynb
@@ -0,0 +1,548 @@
+{
+ "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. components of the Finnis-Sinclair model (pair energy, glue density, glue energy)\n",
+    "3. running a simulation and visualizing properties"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "import Pkg\n",
+    "\n",
+    "Pkg.activate(\".\")\n",
+    "Pkg.instantiate()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "using Molly\n",
+    "using Plots\n",
+    "using Test\n",
+    "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)"
+   ]
+  },
+  {
+   "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": [
+    "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": [
+    "r = collect(range(0, stop=6, length=1000));\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",
+    "end\n",
+    "\n",
+    "plot(r, ɸs, label=hcat(elements...), 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",
+    "\n",
+    "uₙs = []\n",
+    "for el in elements\n",
+    "    uₙ = Molly.Uglue.(ρ, fs84.singles[el].A)\n",
+    "    append!(uₙs,[uₙ])\n",
+    "end\n",
+    "\n",
+    "plot(ρ, uₙs, label=hcat(elements...), 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": [
+    "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",
+    "    append!(element_pairs, [el_pair])\n",
+    "end\n",
+    "\n",
+    "plot(r, Vs, label=hcat(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",
+    "el = \"W\"\n",
+    "\n",
+    "a = bcc_lattice_constants[el]\n",
+    "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)"
+   ]
+  },
+  {
+   "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 supercell.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_bcc.pdb\")\n",
+    ")\n",
+    "\n",
+    "s = Simulation(\n",
+    "    simulator=VelocityVerlet(), \n",
+    "    atoms=supercell.atoms, \n",
+    "    general_inters=general_inters,\n",
+    "    coords=[SVector{3}(v) for v in supercell.positions], \n",
+    "    velocities=velocities,\n",
+    "    temperature=T, \n",
+    "    box_size=supercell.edge_lengths[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": {
+    "scrolled": true
+   },
+   "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
+}
diff --git a/src/analysis.jl b/src/analysis.jl
index ec3bb5e8..c468040a 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
             neighbors = s.neighbors
             @inbounds for ni in 1:length(neighbors)
diff --git a/src/forces.jl b/src/forces.jl
index c8611665..97cd2c00 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
                 neighbors = s.neighbors
                 @threads for ni in 1:length(neighbors)
@@ -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
                 neighbors = s.neighbors
                 for ni in 1:length(neighbors)
@@ -97,8 +108,8 @@ function accelerations(s::Simulation; parallel::Bool=true)
 
     for i in 1:n_atoms
         forces[i] /= s.atoms[i].mass
+        s.forces[i] = forces[i]
     end
-
     return forces
 end
 
@@ -134,3 +145,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..8e91e934
--- /dev/null
+++ b/src/interactions/glue_fs.jl
@@ -0,0 +1,243 @@
+"""
+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
+
+struct FinnisSinclairPair{T<:Real}
+    # container for two element interaction -> c c0 c1 c2
+    c::T
+    c₀::T
+    c₁::T
+    c₂::T
+end
+
+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::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{T<:Real} <: GlueInteraction
+    nl_only::Bool
+    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)
+
+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"]
+    
+    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,
+                  "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
+    )
+
+    kb = 8.617333262145e-5 # eV/K
+    fs84 = FinnisSinclair(nl_only, fs_pairs, fs_singles, kb)
+    
+    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
+
+"""
+    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 
+its neighbouring atoms.
+"""
+function glue(r, β, d)
+    return r > d ? 0 : (r-d)^2 + β*(r-d)^3/d
+end
+
+"""
+    ∂glue_∂r(r, β, d)
+
+Derivative of the glue density function.
+"""
+∂glue_∂r(r, β, d) = r > d ? 0 : 2*(r-d) + 3*β*(r-d)^2/d
+
+"""
+    Uglue(ρ, A)
+
+Energy based on the glue density .
+"""
+function Uglue(ρ, A)
+   return -A * √ρ 
+end
+
+
+"""
+    ∂Uglue_∂ρ(ρ, A)
+
+Energy derivative given glue density.
+"""
+∂Uglue_∂ρ(ρ,A) = - A / (2 * √ρ)
+
+
+"""
+    Upair(r, c, c₀, c₁, c₂)
+
+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₂)
+
+Derivative of the pair potential.
+"""
+∂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)
+
+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.neighbors
+
+        # collecting parameters
+        el_i = s.atoms[i].name
+        el_j = s.atoms[j].name
+        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)
+        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
+    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)
+    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
+    
+    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 = inter.pairs[string(s.atoms[i].name, s.atoms[j].name)]
+    return Upair(r, pij.c, pij.c₀, pij.c₁, pij.c₂)
+end
\ No newline at end of file
diff --git a/src/loggers.jl b/src/loggers.jl
index f418f327..3990da06 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,
+                       Vector{SVector{dims}}[])
+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 0b07d3d7..5b3d89da 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, neighbors, neighbor_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..baf919be
--- /dev/null
+++ b/test/glue.jl
@@ -0,0 +1,104 @@
+using Molly
+using Test
+using SingleCrystal
+
+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, 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, ref)
+    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, 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)
+    
+    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 = SingleCrystal.add_vacancies(supercell, ixs=[1])
+    end
+    n_atoms = length(supercell.positions)
+    
+    velocities = [velocity(atom.mass, T, dims=3) for atom in supercell.atoms]
+    nb_matrix = trues(n_atoms,n_atoms)
+    dist_cutoff = 2. * a
+
+    loggers = Dict(
+        "glue" => GlueDensityLogger(1),
+        "forces" => ForceLogger(5),
+        "energy" => EnergyLogger(1)
+    )
+
+    nf = DistanceNeighborFinder(nb_matrix, 1, dist_cutoff)
+    
+    s = Simulation(
+        simulator=VelocityVerlet(), 
+        atoms=supercell.atoms,  
+        general_inters=(fs_inter,),
+        coords=[SVector{3}(v) for v in supercell.positions],  
+        velocities=velocities,
+        temperature=T, 
+        box_size=supercell.edge_lengths[1], 
+        timestep=.002,
+        n_steps=100,
+        neighbor_finder=nf,
+        loggers=loggers,
+    )
+    simulate!(s)
+    return s
+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
+            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, 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
+    
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index df13b627..14767d66 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -22,6 +22,10 @@ end
 
 CUDA.allowscalar(false) # Check that we never do scalar indexing on the GPU
 
+@testset "glue.jl" begin
+    include("glue.jl")
+end
+
 @testset "Spatial" begin
     @test vector1D(4.0, 6.0, 10.0) ==  2.0
     @test vector1D(1.0, 9.0, 10.0) == -2.0