diff --git a/Project.toml b/Project.toml
index ea2db14ef..1de9df9bd 100644
--- a/Project.toml
+++ b/Project.toml
@@ -21,6 +21,7 @@ GPUArrays = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LAMMPS = "ee2e13b9-eee9-4449-aafa-cfa6a2dbe14d"
 NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
 PeriodicTable = "7b2266bf-644c-5ea3-82d8-af4bbd25a884"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
@@ -41,6 +42,7 @@ Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9"
 GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
 KernelDensity = "5ab0869b-81aa-558d-bb23-cbf5423bbe9b"
 PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d"
+LAMMPS = "ee2e13b9-eee9-4449-aafa-cfa6a2dbe14d"
 
 [extensions]
 MollyCUDAExt = "CUDA"
@@ -48,6 +50,7 @@ MollyEnzymeExt = "Enzyme"
 MollyGLMakieExt = ["GLMakie", "Colors"]
 MollyKernelDensityExt = "KernelDensity"
 MollyPythonCallExt = "PythonCall"
+MollyLAMMPSExt = "LAMMPS"
 
 [compat]
 AcceleratedKernels = "0.4.3"
@@ -71,6 +74,7 @@ GPUArrays = "11"
 Graphs = "1.8"
 KernelAbstractions = "0.9"
 KernelDensity = "0.5, 0.6"
+LAMMPS = "0.9"
 LinearAlgebra = "1.9"
 NearestNeighbors = "0.4"
 PeriodicTable = "1.1"
diff --git a/docs/src/examples.md b/docs/src/examples.md
index 220cb0f22..ed904f7ea 100644
--- a/docs/src/examples.md
+++ b/docs/src/examples.md
@@ -533,6 +533,84 @@ sys = System(
 
 potential_energy(sys) # -2080.2391023908813 eV
 ```
+## LAMMPS Calculator
+
+[`LAMMPSCalculator`](@ref) is made available by using [LAMMPS.jl](https://github.com/cesmix-mit/LAMMPS.jl) with Molly. This allows force/energy calculation using a potential from LAMMPS. This calculator is limited to fully periodic, 3D systems, with rectangular or cubic cells (i.e., no triclinic boundaries). These limitations are not fundamental and can be supported through a pull request. Furthermore, this calculator will always promote floats to Float64 as LAMMPS always uses double precision, and the data must reside on the CPU. 
+
+This calculator will always use a single thread as Molly does not yet depend on MPI. Despite this, LAMMPS.jl expects MPI to be initialized on systems that support MPI. You must call `LAMMPS.MPI.Init()` at the start of your script for [`LAMMPSCalculator`](@ref) to initialize properly. Do NOT try to call Molly with `mpirun` (until the date that is supported)
+
+To use the [`LAMMPSCalculator`](@ref) requires a Molly [`System`](@ref) whose energy and force units are commensurate with the unit system chosen for LAMMPS. For example, if `metal` units are chosen the system should have energy of eV and force units of eV / angstrom. The mass unit is automatically converted. Next you must define your potential. This is simply the string used to define the potential in a LAMMPS input file. For example, for Lennard-Jones you could define:
+
+```julia
+lj_cmds = ["pair_style lj/cut 8.5", "pair_coeff * * 0.0104 3.4", "pair_modify shift yes"]
+```
+
+The optional argument, `label_type_map` allows you to define the atom types that Molly will pass to LAMMPS. This makes it possible to define multi-atomic potentials which require the types to be specified in the definition of the potential. By default, types are assigned in the order they are discovered in the system. For monoatomic systems (this never needs to be defined) the `label_type_map` could be defined as `Dict(:Al => 1)`. 
+
+The `calculate_potential` tells LAMMPS to calculate potential energy. This will be turned on automatically if Molly detects that your system has a `PotentialEnergyLogger` or `TotalEnergyLogger`. For other, less explicit uses of potential energy (e.g. Monte-Carlo) this flag must be set manually. 
+
+
+```julia
+
+using LAMMPS
+using Molly
+using SimpleCrystals
+
+LAMMPS.MPI.Init()
+
+dt = 1.0u"fs"
+damping = 0.5u"ps^-1"
+
+ar_crystal = FCC(5.2468u"Å", :Ar, SVector(4,4,4))
+diamond_crystal = Diamond(5.43u"Å", :Si, SVector(3, 3, 3))
+al_crystal = FCC(4.041u"Å", :Al, SVector(4,4,4))
+
+pot_basepath = abspath(dirname(LAMMPS.locate()), "..", "share", "lammps", "potentials")
+
+eam_pot = joinpath(pot_basepath, "Al_zhou.eam.alloy")
+sw_pot = joinpath(pot_basepath, "Si.sw")
+
+# these are LJ argon params, but will use with silicon just to see if energy the same.
+lj_cmds = ["pair_style lj/cut 8.5", "pair_coeff * * 0.0104 3.4", "pair_modify shift yes"]
+sw_cmds = ["pair_style sw", "pair_coeff * * \"$(sw_pot)\" Si"]
+eam_cmds = ["pair_style eam/alloy", "pair_coeff * * \"$(eam_pot)\" Al"]
+
+pots = (
+    (lj_cmds, ar_crystal, "LJ", -19.85644u"eV"),
+    (sw_cmds, diamond_crystal, "SW", -936.70522u"eV"),
+    (eam_cmds, al_crystal, "EAM", -915.27403u"eV")
+)
+
+for (pot_cmd, crys, pot_type, E_pot) in pots
+
+    sys = System(
+        crys,
+        energy_units = u"eV",
+        force_units = u"eV / angstrom",
+        loggers = (PotentialEnergyLogger(typeof(1.0 * u"eV"), 10),)
+    )
+
+    inter = LAMMPSCalculator(
+        sys,
+        "metal", 
+        pot_cmd;
+        logfile_path = joinpath(@__DIR__, "log.lammps"),
+        calculate_potential = true
+    )
+
+    random_velocities!(sys, 100u"K")
+
+    sys = System(sys; general_inters = (inter, ))
+
+    sim = Langevin(dt = dt, temperature = 100u"K", friction = damping)
+
+    simulate!(sys, sim, 1)
+    PE_LAMMPS = values(sys.loggers[1])[1]
+    println("LAMMPS Calculator through Molly got: $(PE_LAMMPS)")
+    println("Potential energy calculated via C++ LAMMPS: $(E_pot)")
+end
+```
+
 
 ## Density functional theory
 
diff --git a/ext/MollyLAMMPSExt.jl b/ext/MollyLAMMPSExt.jl
new file mode 100644
index 000000000..d91f3afe8
--- /dev/null
+++ b/ext/MollyLAMMPSExt.jl
@@ -0,0 +1,211 @@
+module MollyLAMMPSExt
+
+using LAMMPS
+using Molly
+using Unitful
+import AtomsCalculators
+import AtomsBase
+
+const lammps_mass_units_map = Dict("metal" => u"g/mol", "real" => u"g/mol", "si" => u"kg")
+
+function convert_mass(m, lammps_units)
+    sys_mass_unit = unit(m)
+    lammps_mass_units = get(lammps_mass_units_map, lammps_units, sys_mass_unit)
+    lammps_mass_is_molar = (dimension(lammps_mass_units) == u"𝐌* 𝐍^-1")
+    sys_mass_is_molar = (dimension(sys_mass_unit) == u"𝐌* 𝐍^-1")
+
+    lammps_molar_sys_not = (lammps_mass_is_molar && (!sys_mass_is_molar))
+    sys_molar_lammps_not = (sys_mass_is_molar && (!lammps_mass_is_molar))
+
+    if lammps_molar_sys_not
+        return ustrip.(lammps_mass_units, Unitful.Na * m)
+    elseif sys_molar_lammps_not
+        return ustrip.(lammps_mass_units, m / Unitful.Na) 
+    else # both molar or both non-molar
+        return ustrip.(lammps_mass_units, m)
+    end
+end
+
+function check_lammps_units(energy_unit, force_unit, lammps_units)
+
+    length_unit = inv(force_unit / energy_unit)
+
+    err = (U, LE, LA, EE, EA) -> ArgumentError("You picked $(U) units. Expected length units $(LE) and got $(LA). Expected energy units $(EE) and got $(EA)")
+
+    if lammps_units == "metal"
+        if length_unit != u"angstrom" || energy_unit != u"eV"
+            error(err("metal", u"angstrom", length_unit, u"eV", energy_unit))
+        end
+    elseif lammps_units == "lj"
+        if length_unit != NoUnits || energy_unit != NoUnits
+            error(err("lj", "NoUnits", length_unit, "NoUnits", energy_unit))
+        end
+    elseif lammps_units == "real"
+        if length_unit != u"angstrom" || energy_unit != u"kcal/mol"
+            error(err("real", u"angstrom", length_unit, u"kcal/mol", energy_unit))
+        end
+    elseif lammps_units == "si"
+        if length_unit != u"m" || energy_unit != u"J"
+            error(err("si", u"m", length_unit, u"J", energy_unit))
+        end
+    else
+        error(ArgumentError("Unsupported LAMMPS unit system, $(lammps_units). Expected one of (metal, real, lj, si)."))
+    end
+
+end
+
+function has_logger_with_pe(sys)
+    for logger in sys.loggers
+        if logger.observable == Molly.potential_energy_wrapper || logger.observable == Molly.total_energy_wrapper
+            return true
+        end
+    end
+    return false
+end
+
+function Molly.LAMMPSCalculator(
+        sys::System{3, AT, T},
+        lammps_unit_system::String,
+        potential_definition::Union{String, Array{String}};
+        label_type_map::Dict{Symbol, Int} = Dict{Symbol, Int}(),
+        extra_lammps_commands::Union{String, Array{String}} = "",
+        logfile_path::String = "none",
+        calculate_potential::Bool = false
+    ) where {AT, T}
+
+    # check that we're on CPU
+    if AT != Array
+        error(ArgumentError("LAMMPSCalculator only supports CPU execution."))
+    end
+
+    if sys.boundary == TriclinicBoundary
+        error(ArgumentError("LAMMPSCalculator does not support triclinic systems yet. PRs welcome :)"))
+    end
+
+    if Molly.has_infinite_boundary(sys.boundary)
+        error(ArgumentError("LAMMPSCalculator does not support systems with infinite boundaries. Must be fully periodic. PRs welcome :)"))
+    end
+
+    if length(potential_definition) == 0
+        error(ArgumentError("Cannot pass emptry string as potential definition for LAMMPSCalculator"))
+    end
+
+    if T != Float64
+        @warn "LAMMPS uses Float64, you are using $T. You might incur a penalty from type promotion."
+    end
+
+    check_lammps_units(sys.energy_units, sys.force_units, lammps_unit_system)
+
+    # OpenMP doesnt seem to make a difference...
+    # lmp = LMP(["-screen","none", "-sf", "omp", "-pk", "omp", "$(n_threads)"], LAMMPS.MPI.COMM_WORLD)
+    lmp = LMP(["-screen","none"], LAMMPS.MPI.COMM_WORLD)
+
+    all_syms = Molly.atomic_symbol(sys)
+    unique_syms = unique(all_syms)
+    unique_sym_idxs = Dict(sym => findfirst(x -> x == sym, all_syms) for sym in unique_syms)
+
+    if any(unique_syms .== :unknown)
+        error(ArgumentError("All atoms must have atomic symbols to use LAMMPSCalculator"))
+    end
+    
+    ids = collect(Int32, 1:length(sys))
+    xhi, yhi, zhi = ustrip.(sys.boundary)
+
+    if length(label_type_map) == 0
+        label_type_map = Dict(sym => Int32(i) for (i, sym) in enumerate(unique_syms))
+        types = [label_type_map[sym] for sym in all_syms]
+    else 
+        unique_sym_user = keys(label_type_map)
+        if Set(unique_sym_user) != Set(unique_syms)
+            error(ArgumentError("You provided a label_type_map with $(unique_sym_user) symbols, but" *
+                " the system has $(unique_syms). They must match exactly if you pass label_type_map."))
+        end
+        types = [Int32(label_type_map[sym]) for sym in all_syms]
+    end
+
+    m_lmp = Dict(label_type_map[sym] => convert_mass(sys.masses[i], lammps_unit_system) for (sym, i) in unique_sym_idxs)
+
+    label_map_cmd = "labelmap atom " * join(["$(i) $(sym)" for (sym,i) in label_type_map], " ") 
+
+    setup_cmd = """
+            log $(logfile_path)
+            units $(lammps_unit_system)
+            atom_style atomic
+            atom_modify map array sort 0 0
+        """
+    
+    cell_cmd = """
+            boundary p p p
+            region cell block 0 $(xhi) 0 $(yhi) 0 $(zhi) units box
+            create_box $(length(unique_syms)) cell
+            $(label_map_cmd)
+        """
+    
+    mass_cmd = join(["mass $(type) $(m)" for (type,m) in  m_lmp], "\n")
+
+    if length(extra_lammps_commands) > 0
+        command(lmp, extra_lammps_commands)
+    end
+
+    calculate_potential |= has_logger_with_pe(sys)
+    if calculate_potential
+        command(lmp, "compute pot_e all pe")
+    end
+
+    command(lmp, setup_cmd)
+    command(lmp, cell_cmd)
+    command(lmp, mass_cmd)
+
+    LAMMPS.create_atoms(
+        lmp,
+        reinterpret(reshape, Float64, sys.coords),
+        ids,
+        types
+    )   
+
+    try
+        command(lmp, potential_definition)
+    catch e
+        if startswith(e.msg, "Number of element to type mappings does")
+            @info "Ensure path to potential definition is wrapped in quotes if there are spaces in path."
+        end
+        rethrow(e)
+    end
+
+    # This allows LAMMPS to register the computes/fixes
+    # and build the neighbor list. 
+    command(lmp, "run 0 post no")
+
+    length_unit = inv(sys.force_units / sys.energy_units)
+    return LAMMPSCalculator{typeof(lmp), typeof(sys.energy_units), typeof(length_unit)}(
+                            lmp, -1, sys.energy_units, length_unit)
+end
+
+function maybe_run_lammps_calc!(lammps_calc, r::AbstractVector{T}, step_n) where T
+    if lammps_calc.last_updated != step_n
+        # Send current coordinates to LAMMPS
+        scatter!(lammps_calc.lmp, "x", reinterpret(reshape, Float64, r))
+        # Run a step, will execute the registered computes/fixes
+        command(lammps_calc.lmp, "run 0 pre no post no")
+        lammps_calc.last_updated = step_n
+    end
+end
+
+AtomsCalculators.energy_unit(calc::LAMMPSCalculator) = calc.energy_unit
+AtomsCalculators.length_unit(calc::LAMMPSCalculator) = calc.length_unit
+
+function AtomsCalculators.forces!(fs::AbstractVector{T}, sys, inter::LAMMPSCalculator; step_n=0, kwargs...) where T
+    maybe_run_lammps_calc!(inter, sys.coords, step_n)
+    gather!(inter.lmp, "f", reinterpret(reshape, Float64, fs)) 
+    return fs
+end
+
+function AtomsCalculators.potential_energy(sys, inter::LAMMPSCalculator; step_n=0, kwargs...)
+    maybe_run_lammps_calc!(inter, sys.coords, step_n)
+    return extract_compute(inter.lmp, "pot_e", STYLE_GLOBAL, TYPE_SCALAR)[1] * sys.energy_units
+end
+
+
+end # module MollyLAMMPSExt
+
+
diff --git a/src/types.jl b/src/types.jl
index 2fbde1597..5d1d59945 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -21,7 +21,8 @@ export
     masses,
     charges,
     MollyCalculator,
-    ASECalculator
+    ASECalculator,
+    LAMMPSCalculator
 
 const DefaultFloat = Float64
 
@@ -1484,6 +1485,56 @@ end
 
 function update_ase_calc! end
 
+"""
+    LAMMPSCalculator(
+        sys::System{3, AT, T},
+        lammps_unit_system::String,
+        potential_definition::Union{String, Array{String}};
+        extra_lammps_commands::Union{String, Array{String}} = "",
+        logfile_path::String = "none",
+        calculate_potential::Bool = false
+    )
+
+Defines a general interaction that will call LAMMPS to calculate forces and energies. Forces
+and energies are calculated on a single thread. You must call LAMMPS.MPI.Init() for LAMMPS.jl
+to load the LAMMPS executable on systems where MPI is available.
+
+The LAMMPS potential files can be found at:
+`abspath(dirname(LAMMPS.locate()), "..", "share", "lammps", "potentials")`
+
+Restrictions:
+-------------
+- CPU only
+- Floats promote to Float64
+- No TriclinicBoundary
+- 3D systems only
+- Fully periodic systems only
+
+Arguments:
+----------
+- `sys::System{3}`: The system object this interaction will be applied to. You still have to add this
+    interaction to the system object yourself after constructing the LAMMPSCalculator.
+- `lammps_unit_system::String` : String representing unit system passed to lammps unit command (e.g., metal)
+- `potential_definition::Union{String, Array{String}}` : Commands passed to lammps which define your interaction.
+    For example, to define LJ you pass:
+    `lj_cmds = ["pair_style lj/cut 8.5", "pair_coeff * * 0.0104 3.4", "pair_modify shift yes"]`
+- `label_type_map::Dict{Symbol, Int} = Dict{Symbol, Int}()` : By default atom types are assigned in the 
+    order they appear in the system. This can make defining the potential for multi-atomic systems
+    difficult. By providing this dictionary you can overide the type label assigned to each unique species. 
+- `extra_lammps_commands::Union{String, Array{String}}` : Any other commands you want to pass to LAMMPS. 
+    This feature is untested, but allows flexability to modify the LAMMPS system. 
+- `logfile_path::String = "none"` : Path where LAMMPS logfile is written. Defaults to no log file.
+- `calculate_potential::Bool = false`: Flag to set if LAMMPS should calculate potential energy. Will
+    automatically be turned on if loggers requiring potential energy are detected. Other use cases
+    like Monte-Carlo (and its associated barostat) require setting this flag to true. 
+"""
+mutable struct LAMMPSCalculator{T, E, L}
+    lmp::T # T will be LMP but that is not available here
+    last_updated::Int
+    energy_unit::E
+    length_unit::L
+end
+
 iszero_value(x) = iszero(x)
 
 # Only use threading if a condition is true
diff --git a/test/Project.toml b/test/Project.toml
index 2a8079aae..a6a90c408 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -14,6 +14,7 @@ FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
 GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
 GPUArrays = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7"
 KernelDensity = "5ab0869b-81aa-558d-bb23-cbf5423bbe9b"
+LAMMPS = "ee2e13b9-eee9-4449-aafa-cfa6a2dbe14d"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7"
 Metal = "dde4c033-4e86-420c-a63e-0dd931031962"
diff --git a/test/basic.jl b/test/basic.jl
index 0079b23c6..5169e909c 100644
--- a/test/basic.jl
+++ b/test/basic.jl
@@ -526,3 +526,115 @@ end
         test_forces(ab_sys, calc)
     end
 end
+
+
+@testset "LAMMPSCalculator" begin
+    
+    LAMMPS.MPI.Init()
+
+    dt = 1.0u"fs"
+    damping = 0.5u"ps^-1"
+
+    ar_crystal = FCC(5.2468u"Å", :Ar, SVector(4,4,4))
+    ar_crystal_real = FCC(5.2468u"Å", 39.95u"g/mol", SVector(4,4,4))
+    diamond_crystal = Diamond(5.43u"Å", :Si, SVector(3, 3, 3))
+    al_crystal = FCC(4.041u"Å", :Al, SVector(4,4,4))
+
+    pot_basepath = abspath(dirname(LAMMPS.locate()), "..", "share", "lammps", "potentials")
+
+    eam_pot = joinpath(pot_basepath, "Al_zhou.eam.alloy")
+    sw_pot = joinpath(pot_basepath, "Si.sw")
+
+    # these are LJ argon params, but will use with silicon just to see if energy the same.
+    lj_cmds = ["pair_style lj/cut 8.5", "pair_coeff * * 0.0104 3.4", "pair_modify shift yes"]
+    lj_cmds_real = ["pair_style lj/cut 8.5", "pair_coeff * * 0.24037 3.4", "pair_modify shift yes"]
+    sw_cmds = ["pair_style sw", "pair_coeff * * \"$(sw_pot)\" Si"]
+    eam_cmds = ["pair_style eam/alloy", "pair_coeff * * \"$(eam_pot)\" Al"]
+
+    pots = (
+        (lj_cmds, ar_crystal, "LJ", -19.85644u"eV"),
+        (lj_cmds_real, ar_crystal_real, "LJ-real", -458.93197u"kcal/mol"),
+        (sw_cmds, diamond_crystal, "SW", -936.70522u"eV"),
+        (eam_cmds, al_crystal, "EAM", -915.27403u"eV")
+    )
+
+    for (pot_cmd, crys, pot_type, E_pot) in pots
+
+        unit_sys = pot_type == "LJ-real" ? "real" : "metal"
+
+        sys = System(
+            crys,
+            energy_units = unit(E_pot),
+            force_units = unit(E_pot) / u"angstrom",
+            loggers = (PotentialEnergyLogger(typeof(1.0 * unit(E_pot)), 10),)
+        )
+
+        if pot_type == "LJ_real"
+            new_atoms_data = []
+            for i in eachindex(sys)
+                push!(new_atoms_data, Molly.AtomData(element = "Ar"))
+            end
+
+            sys = System(sys; atoms_data=[new_atoms_data...])
+        end
+
+        inter = LAMMPSCalculator(
+            sys,
+            unit_sys, 
+            pot_cmd;
+            logfile_path = joinpath(@__DIR__, "log.lammps"),
+            calculate_potential = true
+        )
+
+        random_velocities!(sys, 100u"K")
+
+        sys = System(sys; general_inters = (inter, ))
+
+        sim = Langevin(dt = dt, temperature = 100u"K", friction = damping)
+
+        simulate!(sys, sim, 1)
+        PE_LAMMPS = values(sys.loggers[1])[1]
+
+        @test PE_LAMMPS ≈ E_pot
+
+        LAMMPS.close!(inter.lmp)
+    end
+
+    sys = System(
+            ar_crystal,
+            energy_units = u"eV",
+            force_units = u"eV / angstrom",
+            loggers = (PotentialEnergyLogger(typeof(1.0u"eV"), 10),)
+    )
+
+    @test_throws ArgumentError LAMMPSCalculator(
+        sys,
+        "metal", 
+        eam_cmds;
+        label_type_map = Dict(:Al => 1, :Si => 2, :Ar => 3),
+        logfile_path = joinpath(@__DIR__, "log.lammps"),
+        calculate_potential = true
+    )
+
+    @test_throws ArgumentError LAMMPSCalculator(
+        sys,
+        "real", 
+        eam_cmds;
+        label_type_map = Dict(:Al => 1),
+        logfile_path = joinpath(@__DIR__, "log.lammps"),
+        calculate_potential = true
+    )
+
+    @test_throws ArgumentError LAMMPSCalculator(
+        sys,
+        "fake-unit-system", 
+        eam_cmds;
+        label_type_map = Dict(:Al => 1),
+        logfile_path = joinpath(@__DIR__, "log.lammps"),
+        calculate_potential = true
+    )
+
+    LAMMPS.MPI.Finalize()
+
+
+end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 373328e71..19536dba3 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,3 +1,4 @@
+using LAMMPS
 using Molly
 using Molly: from_device, to_device
 using AMDGPU