Merge branch 'JuliaMolSim:master' into master

.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"

.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

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"

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)
 ```

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,31 +92,41 @@ 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)
     TemperatureLogger(T, n_steps)
 
 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
     print(io, "TemperatureLogger{", eltype(values(tl)), "} with n_steps ",
             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)
     CoordinateLogger(T, n_steps; dims=3)
 
 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,15 +135,22 @@ 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)
     VelocityLogger(T, n_steps; dims=3)
 
 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,65 +159,101 @@ 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)
     KineticEnergyLogger(T, n_steps)
 
 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
     print(io, "KineticEnergyLogger{", eltype(values(el)), "} with n_steps ",
             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