Girsanov Integrator

[axsk]

For implementing the optimal control from our article.

[axsk]

We have a first draft in 593178f

I am not sure if the units/magnitudes are correct yet.

[axsk]

Whereas the integrator should work (still untested), we also need a way to expose it to propagate etc.
I currently think of an alternative OpenMMSimulationControlled that has an additional force::Function field.

@joramkuntze: keeping you in the loop :)

[axsk]

It should be working now for the overdamped case, if dt/gamma <= 1e-6

Unfortunately OD isn't directly comparable to UD, so maybe we should implement this now for UD Langevin.
In principle it works the same, just the update step is a little more complicated, cf
https://arxiv.org/abs/2303.14696 (eqs 23 and table I) or
https://pubs.acs.org/doi/full/10.1021/acs.jpcb.4c01702 (eq. 16,17 or table I)

[axsk]

function langevin_girsanov(sim, q, steps, u)
    # from https://pubs.acs.org/doi/full/10.1021/acs.jpcb.4c01702

    kB = 0.008314463
    dt = stepsize(sim)
    ξ = friction(sim)
    T = temp(sim)
    M = repeat(masses(sim), inner=3)

    # Maxwell-Boltzmann distribution
    p = randn(length(M)) .* sqrt(M .* kB .* T)

    t2 = dt / 2
    a = @. t2 / M # eq 18
    d = @. exp(-ξ * dt) # eq 17
    f = @. sqrt(kB * T * M * (1 - exp(-2 * γ * dt))) # eq 17

    b = similar(p)
    η = similar(p)
    Δη = similar(p)
    g = 0

    for i in 1:steps
        randn!(η)
        @. q += a * p # A

        # girsanov part
        F = u(q) # perturbation force ∇U_bias = -F
        @. Δη = (d + 1) / f + dt / 2 * F
        g += η' * Δη + Δη' * Δη / 2
        F .+= force(sim, q) # total force: -∇V - ∇U_bias

        @. b = t2 * F
        @. p += b  # B
        @. p = d * p + f * η # O
        @. p += b  # B
        @. q += a * p # A
    end

    return q, exp(-g)
end