Skip to content

Commit

Permalink
added overdamped langevin and sketched potentialfunction
Browse files Browse the repository at this point in the history
  • Loading branch information
axsk committed Sep 28, 2023
1 parent 011bcc3 commit 7445ca5
Show file tree
Hide file tree
Showing 3 changed files with 74 additions and 0 deletions.
1 change: 1 addition & 0 deletions src/Molly.jl
Original file line number Diff line number Diff line change
Expand Up @@ -59,6 +59,7 @@ include("interactions/periodic_torsion.jl")
include("interactions/rb_torsion.jl")
include("interactions/implicit_solvent.jl")
include("interactions/muller_brown.jl")
include("interactions/potential.jl")
include("energy.jl")
include("constraints.jl")
include("simulators.jl")
Expand Down
22 changes: 22 additions & 0 deletions src/interactions/potential.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,22 @@
struct PotentialFunction{F}
f::F
# fwdiff gradient caches
end

function Molly.potential_energy(V::PotentialFunction,
sys,
neighbors=nothing;
n_threads=Threads.nthreads())

return V.f(sys.coords)
end

function Molly.forces(V::PotentialFunction,
sys,
neighbors=nothing;
n_threads=Threads.nthreads())

return ForwardDiff.gradient(V.f, sys.coords)
end


51 changes: 51 additions & 0 deletions src/simulators.jl
Original file line number Diff line number Diff line change
Expand Up @@ -930,3 +930,54 @@ function random_unit_vector(float_type, dims)
vec = randn(float_type, dims)
return vec / norm(vec)
end


"""
EulerMaruyama(; <keyword arguments>)
A simulator for a stochastic differential equation (SDE) using the Euler-Maruyama method.
# Arguments
- `dt::S`: the time step of the simulation.
- `temperature::K`: the equilibrium temperature of the simulation.
- `friction::F`: the friction coefficient of the simulation.
- `remove_CM_motion::Bool=true`: whether to remove the center of mass motion
every time step.
"""
@kwdef struct OverdampedLangevin{T,K,F}
dt::S
temperature::K
friction::F
remove_CM_motion::Bool
end

function simulate!(sys,
sim::OverdampedLangevin,
n_steps::Integer;
n_threads::Integer=Threads.nthreads(),
rng=Random.GLOBAL_RNG)

sys.coords = wrap_coords.(sys.coords, (sys.boundary,))

sim.remove_CM_motion && remove_CM_motion!(sys)
neighbors = find_neighbors(sys, sys.neighbor_finder; n_threads=n_threads)
run_loggers!(sys, neighbors, 0; n_threads=n_threads)

for step_n in 1:n_steps
a = accelerations(sys, neighbors; n_threads=n_threads)
# noise = kB * T / M
noise = random_velocities(sys, sim.temperature; rng=rng)
sys.coords += a ./ sim.friction .* sim.dt .+ sqrt(2 / sim.friction * sim.dt) .* noise

sys.coords = wrap_coords.(sys.coords, (sys.boundary,))
sim.remove_CM_motion && remove_CM_motion!(sys)

run_loggers!(sys, neighbors, step_n; n_threads=n_threads)

if step_n != n_steps
neighbors = find_neighbors(sys, sys.neighbor_finder, neighbors, step_n;
n_threads=n_threads)
end
end
return sys
end

0 comments on commit 7445ca5

Please sign in to comment.