diff --git a/.github/dependabot.yml b/.github/dependabot.yml
new file mode 100644
index 00000000..dc98455d
--- /dev/null
+++ b/.github/dependabot.yml
@@ -0,0 +1,7 @@
+# https://docs.github.com/code-security/dependabot/dependabot-version-updates/configuration-options-for-the-dependabot.yml-file
+version: 2
+updates:
+  - package-ecosystem: "github-actions"
+    directory: "/" # Location of package manifests
+    schedule:
+      interval: "weekly"
diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 472f6886..3dbe474c 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -28,37 +28,29 @@ jobs:
           - NotZygote
           - Zygote
     steps:
-      - uses: actions/checkout@v2
-      - uses: julia-actions/setup-julia@v1
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
-      - uses: actions/cache@v1
-        env:
-          cache-name: cache-artifacts
-        with:
-          path: ~/.julia/artifacts
-          key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }}
-          restore-keys: |
-            ${{ runner.os }}-test-${{ env.cache-name }}-
-            ${{ runner.os }}-test-
-            ${{ runner.os }}-
+      - uses: julia-actions/cache@v1
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
         env:
           GROUP: ${{ matrix.test }}
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v1
+      - uses: codecov/codecov-action@v4
         with:
-          file: lcov.info
+          token: ${{ secrets.CODECOV_TOKEN }}
   docs:
     name: Documentation
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v2
-      - uses: julia-actions/setup-julia@v1
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v2
         with:
           version: '1'
+      - uses: julia-actions/cache@v1
       - run: |
           julia --project=docs -e '
             using Pkg
diff --git a/Project.toml b/Project.toml
index 07d58d0c..cfc28237 100644
--- a/Project.toml
+++ b/Project.toml
@@ -55,7 +55,7 @@ Combinatorics = "1"
 DataStructures = "0.18"
 Distances = "0.10"
 Distributions = "0.23, 0.24, 0.25"
-Enzyme = "0.11.15"
+Enzyme = "0.11.15, 0.12"
 EzXML = "1"
 FLoops = "0.2"
 ForwardDiff = "0.10.35"
diff --git a/docs/src/documentation.md b/docs/src/documentation.md
index 2b45bd18..c558e6fd 100644
--- a/docs/src/documentation.md
+++ b/docs/src/documentation.md
@@ -889,6 +889,7 @@ function Molly.simulate!(sys,
         remove_CM_motion!(sys)
 
         # Apply the loggers like this
+        # Computed quantities can also be given as keyword arguments to run_loggers!
         run_loggers!(sys, neighbors, step_n, run_loggers; n_threads=n_threads)
 
         # Find new neighbors like this
@@ -1064,14 +1065,32 @@ Running loggers can be disabled entirely with `run_loggers=false`, which is the
 Loggers are currently ignored for the purposes of taking gradients, so if a logger is used in the gradient calculation the gradients will appear to be nothing.
 
 Many times, a logger will just record an observation to an `Array` containing a record of past observations.
-For this purpose, you can use the [`GeneralObservableLogger`](@ref) without defining a custom logging function. Simply define your observation function as
+For this purpose, you can use the [`GeneralObservableLogger`](@ref) without defining a custom logging function.
+Define your observation function as
 ```julia
-function my_observable(sys::System, neighbors; n_threads::Integer)
+function my_observable(sys::System, neighbors; n_threads::Integer, kwargs...)
     # Probe the system for some desired property
     return observation
 end
 ```
-A logger which records this property every `n_steps` can be constructed through 
+Keyword arguments `current_forces` and `current_potential_energy` can also be used here to avoid recomputing values that are passed from the simulator:
+```julia
+function my_pe_observable(sys::System, neighbors; n_threads::Integer,
+                          current_potential_energy=nothing, kwargs...)
+    if isnothing(current_potential_energy)
+        # Compute potential energy
+        return potential_energy(sys, neighbors; n_threads=n_threads)
+    else
+        # Potential energy was passed from simulator, reuse
+        return current_potential_energy
+    end
+end
+```
+These keyword arguments are also available in [`log_property!`](@ref).
+Which values are passed depends on the simulator being used, for example [`SteepestDescentMinimizer`](@ref) passes `current_potential_energy` because it uses it for minimization.
+Note that loggers are called after [`apply_coupling!`](@ref), so the coordinates may have changed since the potential energy or forces were computed.
+
+A logger which records the property every `n_steps` can be constructed through 
 ```julia
 my_logger = GeneralObservableLogger(my_observable, T, n_steps)
 ```
diff --git a/src/loggers.jl b/src/loggers.jl
index eb26480d..283a80fd 100644
--- a/src/loggers.jl
+++ b/src/loggers.jl
@@ -79,9 +79,9 @@ Custom loggers should implement this function.
 Additional keyword arguments can be passed to the logger if required.
 """
 function log_property!(logger::GeneralObservableLogger, s::System, neighbors=nothing,
-                        step_n::Integer=0; n_threads::Integer=Threads.nthreads(), kwargs...)
+                        step_n::Integer=0; kwargs...)
     if (step_n % logger.n_steps) == 0
-        obs = logger.observable(s, neighbors; n_threads=n_threads)
+        obs = logger.observable(s, neighbors; kwargs...)
         push!(logger.history, obs)
     end
 end
@@ -92,7 +92,7 @@ function Base.show(io::IO, gol::GeneralObservableLogger)
             gol.observable)
 end
 
-temperature_wrapper(s, neighbors=nothing; n_threads::Integer=Threads.nthreads()) = temperature(s)
+temperature_wrapper(sys, neighbors; kwargs...) = temperature(sys)
 
 """
     TemperatureLogger(n_steps)
@@ -100,7 +100,10 @@ temperature_wrapper(s, neighbors=nothing; n_threads::Integer=Threads.nthreads())
 
 Log the [`temperature`](@ref) throughout a simulation.
 """
-TemperatureLogger(T::DataType, n_steps::Integer) = GeneralObservableLogger(temperature_wrapper, T, n_steps)
+function TemperatureLogger(T::DataType, n_steps::Integer)
+    return GeneralObservableLogger(temperature_wrapper, T, n_steps)
+end
+
 TemperatureLogger(n_steps::Integer) = TemperatureLogger(typeof(one(DefaultFloat)u"K"), n_steps)
 
 function Base.show(io::IO, tl::GeneralObservableLogger{T, typeof(temperature_wrapper)}) where T
@@ -108,7 +111,7 @@ function Base.show(io::IO, tl::GeneralObservableLogger{T, typeof(temperature_wra
             tl.n_steps, ", ", length(values(tl)), " temperatures recorded")
 end
 
-coordinates_wrapper(sys, neighbors=nothing; n_threads::Integer=Threads.nthreads()) = sys.coords
+coordinates_wrapper(sys, neighbors; kwargs...) = sys.coords
 
 """
     CoordinateLogger(n_steps; dims=3)
@@ -116,7 +119,14 @@ coordinates_wrapper(sys, neighbors=nothing; n_threads::Integer=Threads.nthreads(
 
 Log the coordinates throughout a simulation.
 """
-CoordinateLogger(T, n_steps::Integer; dims::Integer=3) = GeneralObservableLogger(coordinates_wrapper, Array{SArray{Tuple{dims}, T, 1, dims}, 1}, n_steps)
+function CoordinateLogger(T, n_steps::Integer; dims::Integer=3)
+    return GeneralObservableLogger(
+        coordinates_wrapper,
+        Array{SArray{Tuple{dims}, T, 1, dims}, 1},
+        n_steps,
+    )
+end
+
 CoordinateLogger(n_steps::Integer; dims::Integer=3) = CoordinateLogger(typeof(one(DefaultFloat)u"nm"), n_steps; dims=dims)
 
 function Base.show(io::IO, cl::GeneralObservableLogger{T, typeof(coordinates_wrapper)}) where T
@@ -125,7 +135,7 @@ function Base.show(io::IO, cl::GeneralObservableLogger{T, typeof(coordinates_wra
             length(values(cl)) > 0 ? length(first(values(cl))) : "?", " atoms")
 end
 
-velocities_wrapper(sys::System, neighbors=nothing; n_threads::Integer=Threads.nthreads()) = sys.velocities
+velocities_wrapper(sys, neighbors; kwargs...) = sys.velocities
 
 """
     VelocityLogger(n_steps; dims=3)
@@ -133,7 +143,14 @@ velocities_wrapper(sys::System, neighbors=nothing; n_threads::Integer=Threads.nt
 
 Log the velocities throughout a simulation.
 """
-VelocityLogger(T, n_steps::Integer; dims::Integer=3) = GeneralObservableLogger(velocities_wrapper, Array{SArray{Tuple{dims}, T, 1, dims}, 1}, n_steps)
+function VelocityLogger(T, n_steps::Integer; dims::Integer=3)
+    return GeneralObservableLogger(
+        velocities_wrapper,
+        Array{SArray{Tuple{dims}, T, 1, dims}, 1},
+        n_steps,
+    )
+end
+
 VelocityLogger(n_steps::Integer; dims::Integer=3) = VelocityLogger(typeof(one(DefaultFloat)u"nm * ps^-1"), n_steps; dims=dims)
 
 function Base.show(io::IO, vl::GeneralObservableLogger{T, typeof(velocities_wrapper)}) where T
@@ -142,21 +159,7 @@ function Base.show(io::IO, vl::GeneralObservableLogger{T, typeof(velocities_wrap
             length(values(vl)) > 0 ? length(first(values(vl))) : "?", " atoms")
 end
 
-"""
-    TotalEnergyLogger(n_steps)
-    TotalEnergyLogger(T, n_steps)
-
-Log the [`total_energy`](@ref) of a system throughout a simulation.
-"""
-TotalEnergyLogger(T::DataType, n_steps) = GeneralObservableLogger(total_energy, T, n_steps)
-TotalEnergyLogger(n_steps) = TotalEnergyLogger(typeof(one(DefaultFloat)u"kJ * mol^-1"), n_steps)
-
-function Base.show(io::IO, el::GeneralObservableLogger{T, typeof(total_energy)}) where T
-    print(io, "TotalEnergyLogger{", eltype(values(el)), "} with n_steps ",
-            el.n_steps, ", ", length(values(el)), " energies recorded")
-end
-
-kinetic_energy_wrapper(s::System, neighbors=nothing; n_threads::Integer=Threads.nthreads()) = kinetic_energy(s)
+kinetic_energy_wrapper(sys, neighbors; kwargs...) = kinetic_energy(sys)
 
 """
     KineticEnergyLogger(n_steps)
@@ -164,7 +167,10 @@ kinetic_energy_wrapper(s::System, neighbors=nothing; n_threads::Integer=Threads.
 
 Log the [`kinetic_energy`](@ref) of a system throughout a simulation.
 """
-KineticEnergyLogger(T::Type, n_steps::Integer) = GeneralObservableLogger(kinetic_energy_wrapper, T, n_steps)
+function KineticEnergyLogger(T::Type, n_steps::Integer)
+    return GeneralObservableLogger(kinetic_energy_wrapper, T, n_steps)
+end
+
 KineticEnergyLogger(n_steps::Integer) = KineticEnergyLogger(typeof(one(DefaultFloat)u"kJ * mol^-1"), n_steps)
 
 function Base.show(io::IO, el::GeneralObservableLogger{T, typeof(kinetic_energy_wrapper)}) where T
@@ -172,35 +178,82 @@ function Base.show(io::IO, el::GeneralObservableLogger{T, typeof(kinetic_energy_
             el.n_steps, ", ", length(values(el)), " energies recorded")
 end
 
+function potential_energy_wrapper(sys, neighbors; n_threads::Integer,
+                                  current_potential_energy=nothing, kwargs...)
+    if isnothing(current_potential_energy)
+        return potential_energy(sys, neighbors; n_threads=n_threads)
+    else
+        return current_potential_energy
+    end
+end
+
 """
     PotentialEnergyLogger(n_steps)
     PotentialEnergyLogger(T, n_steps)
 
 Log the [`potential_energy`](@ref) of a system throughout a simulation.
 """
-PotentialEnergyLogger(T::Type, n_steps::Integer) = GeneralObservableLogger(potential_energy, T, n_steps)
+function PotentialEnergyLogger(T::Type, n_steps::Integer)
+    return GeneralObservableLogger(potential_energy_wrapper, T, n_steps)
+end
+
 PotentialEnergyLogger(n_steps::Integer) = PotentialEnergyLogger(typeof(one(DefaultFloat)u"kJ * mol^-1"), n_steps)
 
-function Base.show(io::IO, el::GeneralObservableLogger{T, typeof(potential_energy)}) where T
+function Base.show(io::IO, el::GeneralObservableLogger{T, typeof(potential_energy_wrapper)}) where T
     print(io, "PotentialEnergyLogger{", eltype(values(el)), "} with n_steps ",
             el.n_steps, ", ", length(values(el)), " energies recorded")
 end
 
+function total_energy_wrapper(sys, neighbors; kwargs...)
+    return kinetic_energy(sys) + potential_energy_wrapper(sys, neighbors; kwargs...)
+end
+
+"""
+    TotalEnergyLogger(n_steps)
+    TotalEnergyLogger(T, n_steps)
+
+Log the [`total_energy`](@ref) of a system throughout a simulation.
+"""
+TotalEnergyLogger(T::DataType, n_steps) = GeneralObservableLogger(total_energy_wrapper, T, n_steps)
+TotalEnergyLogger(n_steps) = TotalEnergyLogger(typeof(one(DefaultFloat)u"kJ * mol^-1"), n_steps)
+
+function Base.show(io::IO, el::GeneralObservableLogger{T, typeof(total_energy_wrapper)}) where T
+    print(io, "TotalEnergyLogger{", eltype(values(el)), "} with n_steps ",
+            el.n_steps, ", ", length(values(el)), " energies recorded")
+end
+
+function forces_wrapper(sys, neighbors; n_threads::Integer, current_forces=nothing, kwargs...)
+    if isnothing(current_forces)
+        return forces(sys, neighbors; n_threads=n_threads)
+    else
+        return current_forces
+    end
+end
+
 """
     ForceLogger(n_steps; dims=3)
     ForceLogger(T, n_steps; dims=3)
 
 Log the [`forces`](@ref) throughout a simulation.
 """
-ForceLogger(T, n_steps::Integer; dims::Integer=3) = GeneralObservableLogger(forces, Array{SArray{Tuple{dims}, T, 1, dims}, 1}, n_steps)
+function ForceLogger(T, n_steps::Integer; dims::Integer=3)
+    return GeneralObservableLogger(
+        forces_wrapper,
+        Array{SArray{Tuple{dims}, T, 1, dims}, 1},
+        n_steps,
+    )
+end
+
 ForceLogger(n_steps::Integer; dims::Integer=3) = ForceLogger(typeof(one(DefaultFloat)u"kJ * mol^-1 * nm^-1"), n_steps; dims=dims)
 
-function Base.show(io::IO, fl::GeneralObservableLogger{T, typeof(forces)}) where T
+function Base.show(io::IO, fl::GeneralObservableLogger{T, typeof(forces_wrapper)}) where T
     print(io, "ForceLogger{", eltype(eltype(values(fl))), "} with n_steps ",
             fl.n_steps, ", ", length(values(fl)), " frames recorded for ",
             length(values(fl)) > 0 ? length(first(values(fl))) : "?", " atoms")
 end
 
+virial_wrapper(sys, neighbors; n_threads, kwargs...) = virial(sys, neighbors; n_threads=n_threads)
+
 """
     VirialLogger(n_steps)
     VirialLogger(T, n_steps)
@@ -211,14 +264,16 @@ This should only be used on systems containing just pairwise interactions, or
 where the specific interactions, general interactions and constraints do not
 contribute to the virial.
 """
-VirialLogger(T::Type, n_steps::Integer) = GeneralObservableLogger(virial, T, n_steps)
+VirialLogger(T::Type, n_steps::Integer) = GeneralObservableLogger(virial_wrapper, T, n_steps)
 VirialLogger(n_steps::Integer) = VirialLogger(typeof(one(DefaultFloat)u"kJ * mol^-1"), n_steps)
 
-function Base.show(io::IO, vl::GeneralObservableLogger{T, typeof(virial)}) where T
+function Base.show(io::IO, vl::GeneralObservableLogger{T, typeof(virial_wrapper)}) where T
     print(io, "VirialLogger{", eltype(values(vl)), "} with n_steps ",
             vl.n_steps, ", ", length(values(vl)), " virials recorded")
 end
 
+pressure_wrapper(sys, neighbors; n_threads, kwargs...) = pressure(sys, neighbors; n_threads=n_threads)
+
 """
     PressureLogger(n_steps)
     PressureLogger(T, n_steps)
@@ -229,10 +284,10 @@ This should only be used on systems containing just pairwise interactions, or
 where the specific interactions, general interactions and constraints do not
 contribute to the pressure.
 """
-PressureLogger(T::Type, n_steps::Integer) = GeneralObservableLogger(pressure, T, n_steps)
+PressureLogger(T::Type, n_steps::Integer) = GeneralObservableLogger(pressure_wrapper, T, n_steps)
 PressureLogger(n_steps::Integer) = PressureLogger(typeof(one(DefaultFloat)u"bar"), n_steps)
 
-function Base.show(io::IO, pl::GeneralObservableLogger{T, typeof(pressure)}) where T
+function Base.show(io::IO, pl::GeneralObservableLogger{T, typeof(pressure_wrapper)}) where T
     print(io, "PressureLogger{", eltype(values(pl)), "} with n_steps ",
             pl.n_steps, ", ", length(values(pl)), " pressures recorded")
 end
@@ -404,9 +459,9 @@ end
 
 function log_property!(logger::TimeCorrelationLogger, s::System, neighbors=nothing,
                         step_n::Integer=0; n_threads::Integer=Threads.nthreads(), kwargs...)
-    A = logger.observableA(s, neighbors; n_threads=n_threads)
+    A = logger.observableA(s, neighbors; n_threads=n_threads, kwargs...)
     if logger.observableA != logger.observableB
-        B = logger.observableB(s, neighbors; n_threads=n_threads)
+        B = logger.observableB(s, neighbors; n_threads=n_threads, kwargs...)
     else
         B = A
     end
@@ -504,9 +559,9 @@ function Base.values(aol::AverageObservableLogger; std::Bool=true)
 end
 
 function log_property!(aol::AverageObservableLogger{T}, s::System, neighbors=nothing,
-                        step_n::Integer=0; n_threads::Integer=Threads.nthreads(), kwargs...) where T
+                        step_n::Integer=0; kwargs...) where T
     if (step_n % aol.n_steps) == 0
-        obs = aol.observable(s, neighbors; n_threads=n_threads)
+        obs = aol.observable(s, neighbors; kwargs...)
         push!(aol.current_block, obs)
 
         if length(aol.current_block) == aol.current_block_size
@@ -602,7 +657,6 @@ function log_property!(mcl::MonteCarloLogger{T},
                         step_n::Integer=0;
                         success::Bool,
                         energy_rate::T,
-                        n_threads::Integer=Threads.nthreads(),
                         kwargs...) where T
     mcl.n_trials += 1
     if success
diff --git a/src/simulators.jl b/src/simulators.jl
index 56d8057a..2b74e2b7 100644
--- a/src/simulators.jl
+++ b/src/simulators.jl
@@ -65,9 +65,9 @@ function simulate!(sys,
                     run_loggers=false)
     sys.coords = wrap_coords.(sys.coords, (sys.boundary,))
     neighbors = find_neighbors(sys, sys.neighbor_finder; n_threads=n_threads)
-    run_loggers!(sys, neighbors, 0, run_loggers; n_threads=n_threads)
-    using_constraints = length(sys.constraints) > 0
     E = potential_energy(sys, neighbors; n_threads=n_threads)
+    run_loggers!(sys, neighbors, 0, run_loggers; n_threads=n_threads, current_potential_energy=E)
+    using_constraints = length(sys.constraints) > 0
     println(sim.log_stream, "Step 0 - potential energy ",
             E, " - max force N/A - N/A")
     hn = sim.step_size
@@ -98,7 +98,8 @@ function simulate!(sys,
                     E_trial, " - max force ", max_force, " - rejected")
         end
 
-        run_loggers!(sys, neighbors, step_n, run_loggers; n_threads=n_threads)
+        run_loggers!(sys, neighbors, step_n, run_loggers; n_threads=n_threads,
+                     current_potential_energy=E)
 
         if max_force < sim.tol
             break
@@ -136,9 +137,11 @@ function simulate!(sys,
     sys.coords = wrap_coords.(sys.coords, (sys.boundary,))
     !iszero(sim.remove_CM_motion) && remove_CM_motion!(sys)
     neighbors = find_neighbors(sys, sys.neighbor_finder; n_threads=n_threads)
-    run_loggers!(sys, neighbors, 0, run_loggers; n_threads=n_threads)
-    accels_t = accelerations(sys, neighbors; n_threads=n_threads)
+    forces_t = forces(sys, neighbors; n_threads=n_threads)
+    forces_t_dt = zero(forces_t)
+    accels_t = forces_t ./ masses(sys)
     accels_t_dt = zero(accels_t)
+    run_loggers!(sys, neighbors, 0, run_loggers; n_threads=n_threads, current_forces=forces_t)
     using_constraints = length(sys.constraints) > 0
     if using_constraints
         cons_coord_storage = similar(sys.coords)
@@ -155,7 +158,8 @@ function simulate!(sys,
                                                          sim.dt; n_threads=n_threads)
         sys.coords = wrap_coords.(sys.coords, (sys.boundary,))
 
-        accels_t_dt = accelerations(sys, neighbors; n_threads=n_threads)
+        forces_t_dt = forces(sys, neighbors; n_threads=n_threads)
+        accels_t_dt = forces_t_dt ./ masses(sys)
 
         sys.velocities += ((accels_t .+ accels_t_dt) .* sim.dt / 2)
         using_constraints && apply_velocity_constraints!(sys; n_threads=n_threads)
@@ -169,12 +173,15 @@ function simulate!(sys,
         neighbors = find_neighbors(sys, sys.neighbor_finder, neighbors, step_n, recompute_forces;
                                    n_threads=n_threads)
         if recompute_forces
-            accels_t = accelerations(sys, neighbors; n_threads=n_threads)
+            forces_t = forces(sys, neighbors; n_threads=n_threads)
+            accels_t = forces_t ./ masses(sys)
         else
+            forces_t = forces_t_dt
             accels_t = accels_t_dt
         end
 
-        run_loggers!(sys, neighbors, step_n, run_loggers; n_threads=n_threads)
+        run_loggers!(sys, neighbors, step_n, run_loggers; n_threads=n_threads,
+                     current_forces=forces_t)
     end
     return sys
 end
@@ -218,7 +225,8 @@ function simulate!(sys,
     end
 
     for step_n in 1:n_steps
-        accels_t = accelerations(sys, neighbors; n_threads=n_threads)
+        forces_t = forces(sys, neighbors; n_threads=n_threads)
+        accels_t = forces_t ./ masses(sys)
 
         sys.velocities += accels_t .* sim.dt
 
@@ -243,7 +251,8 @@ function simulate!(sys,
         neighbors = find_neighbors(sys, sys.neighbor_finder, neighbors, step_n, recompute_forces;
                                    n_threads=n_threads)
 
-        run_loggers!(sys, neighbors, step_n, run_loggers; n_threads=n_threads)
+        run_loggers!(sys, neighbors, step_n, run_loggers; n_threads=n_threads,
+                     current_forces=forces_t)
     end
     return sys
 end
@@ -284,7 +293,8 @@ function simulate!(sys,
     end
 
     for step_n in 1:n_steps
-        accels_t = accelerations(sys, neighbors; n_threads=n_threads)
+        forces_t = forces(sys, neighbors; n_threads=n_threads)
+        accels_t = forces_t ./ masses(sys)
 
         coords_copy = sys.coords
         if using_constraints
@@ -312,7 +322,8 @@ function simulate!(sys,
                                    n_threads=n_threads)
         coords_last = coords_copy
 
-        run_loggers!(sys, neighbors, step_n, run_loggers; n_threads=n_threads)
+        run_loggers!(sys, neighbors, step_n, run_loggers; n_threads=n_threads,
+                     current_forces=forces_t)
     end
     return sys
 end
@@ -367,7 +378,8 @@ function simulate!(sys,
     end
 
     for step_n in 1:n_steps
-        accels_t = accelerations(sys, neighbors; n_threads=n_threads)
+        forces_t = forces(sys, neighbors; n_threads=n_threads)
+        accels_t = forces_t ./ masses(sys)
 
         sys.velocities += accels_t .* sim.dt
         apply_velocity_constraints!(sys; n_threads=n_threads)
@@ -395,7 +407,8 @@ function simulate!(sys,
         neighbors = find_neighbors(sys, sys.neighbor_finder, neighbors, step_n, recompute_forces;
                                    n_threads=n_threads)
 
-        run_loggers!(sys, neighbors, step_n, run_loggers; n_threads=n_threads)
+        run_loggers!(sys, neighbors, step_n, run_loggers; n_threads=n_threads,
+                     current_forces=forces_t)
     end
     return sys
 end
@@ -573,7 +586,8 @@ function simulate!(sys,
     run_loggers!(sys, neighbors, 0, run_loggers; n_threads=n_threads)
 
     for step_n in 1:n_steps
-        accels_t = accelerations(sys, neighbors; n_threads=n_threads)
+        forces_t = forces(sys, neighbors; n_threads=n_threads)
+        accels_t = forces_t ./ masses(sys)
 
         noise = random_velocities(sys, sim.temperature; rng=rng)
         sys.coords += (accels_t ./ sim.friction) .* sim.dt .+ sqrt((2 / sim.friction) * sim.dt) .* noise
@@ -586,7 +600,8 @@ function simulate!(sys,
         neighbors = find_neighbors(sys, sys.neighbor_finder, neighbors, step_n;
                                    n_threads=n_threads)
 
-        run_loggers!(sys, neighbors, step_n, run_loggers; n_threads=n_threads)
+        run_loggers!(sys, neighbors, step_n, run_loggers; n_threads=n_threads,
+                     current_forces=forces_t)
     end
     return sys
 end
@@ -632,9 +647,11 @@ function simulate!(sys, sim::NoseHoover, n_steps::Integer;
     sys.coords = wrap_coords.(sys.coords, (sys.boundary,))
     !iszero(sim.remove_CM_motion) && remove_CM_motion!(sys)
     neighbors = find_neighbors(sys, sys.neighbor_finder; n_threads=n_threads)
-    run_loggers!(sys, neighbors, 0, run_loggers; n_threads=n_threads)
-    accels_t = accelerations(sys, neighbors; n_threads=n_threads)
+    forces_t = forces(sys, neighbors; n_threads=n_threads)
+    forces_t_dt = zero(forces_t)
+    accels_t = forces_t ./ masses(sys)
     accels_t_dt = zero(accels_t)
+    run_loggers!(sys, neighbors, 0, run_loggers; n_threads=n_threads, current_forces=forces_t)
     v_half = zero(sys.velocities)
     zeta = zero(inv(sim.dt))
 
@@ -649,7 +666,8 @@ function simulate!(sys, sim::NoseHoover, n_steps::Integer;
         T_half = uconvert(unit(sim.temperature), 2 * KE_half / (sys.df * sys.k))
         zeta = zeta_half + (sim.dt / (2 * (sim.damping^2))) * ((T_half / sim.temperature) - 1)
 
-        accels_t_dt = accelerations(sys, neighbors; n_threads=n_threads)
+        forces_t_dt = forces(sys, neighbors; n_threads=n_threads)
+        accels_t_dt = forces_t_dt ./ masses(sys)
 
         sys.velocities = (v_half .+ accels_t_dt .* (sim.dt / 2)) ./
                          (1 + (zeta * sim.dt / 2))
@@ -663,12 +681,15 @@ function simulate!(sys, sim::NoseHoover, n_steps::Integer;
         neighbors = find_neighbors(sys, sys.neighbor_finder, neighbors, step_n, recompute_forces;
                                     n_threads=n_threads)
         if recompute_forces
-            accels_t = accelerations(sys, neighbors; n_threads=n_threads)
+            forces_t = forces(sys, neighbors; n_threads=n_threads)
+            accels_t = forces_t ./ masses(sys)
         else
+            forces_t = forces_t_dt
             accels_t = accels_t_dt
         end
 
-        run_loggers!(sys, neighbors, step_n, run_loggers; n_threads=n_threads)
+        run_loggers!(sys, neighbors, step_n, run_loggers; n_threads=n_threads,
+                     current_forces=forces_t)
     end
     return sys
 end
@@ -985,14 +1006,15 @@ function simulate!(sys::System{D, G, T},
 
         ΔE = E_new - E_old
         δ = ΔE / (sys.k * sim.temperature)
-        if δ < 0 || rand() < exp(-δ)
-
-            run_loggers!(sys, neighbors, i, run_loggers; n_threads=n_threads, success=true,
+        if δ < 0 || (rand() < exp(-δ))
+            run_loggers!(sys, neighbors, i, run_loggers; n_threads=n_threads,
+                         current_potential_energy=E_new, success=true,
                          energy_rate=(E_new / (sys.k * sim.temperature)))
             E_old = E_new
         else
             sys.coords = coords_old
-            run_loggers!(sys, neighbors, i, run_loggers; n_threads=n_threads, success=false,
+            run_loggers!(sys, neighbors, i, run_loggers; n_threads=n_threads,
+                         current_potential_energy=E_old, success=false,
                          energy_rate=(E_old / (sys.k * sim.temperature)))
         end
     end
diff --git a/test/simulation.jl b/test/simulation.jl
index 0472ca92..06fc09ae 100644
--- a/test/simulation.jl
+++ b/test/simulation.jl
@@ -4,7 +4,7 @@
     temp = 298.0u"K"
     boundary = RectangularBoundary(2.0u"nm")
     simulator = VelocityVerlet(dt=0.002u"ps", coupling=AndersenThermostat(temp, 10.0u"ps"))
-    gen_temp_wrapper(s, neighbors=nothing; n_threads::Integer=Threads.nthreads()) = temperature(s)
+    gen_temp_wrapper(s, neighbors; kwargs...) = temperature(s)
 
     s = System(
         atoms=[Atom(charge=0.0, mass=10.0u"g/mol", σ=0.3u"nm", ϵ=0.2u"kJ * mol^-1") for i in 1:n_atoms],
@@ -800,7 +800,7 @@ end
         loggers=(
             coords=CoordinateLogger(10),
             mcl=MonteCarloLogger(),
-            avgpe=AverageObservableLogger(potential_energy, typeof(atoms[1].ϵ), 10),
+            avgpe=AverageObservableLogger(Molly.potential_energy_wrapper, typeof(atoms[1].ϵ), 10),
         ),
     )