ISOKANN.jl

OpenMMSimulation(; pdb, steps, ...)
OpenMMSimulation(; py, steps)

Constructs an OpenMM simulation object. Either use OpenMMSimulation(;py, steps) where pyis the location of a .py python script creating a OpenMM simulation object or supply a .pdb file viapdb` and the following parameters (see also defaultsystem):

Arguments

• pdb::String: Path to the PDB file.
• ligand::String: Path to ligand file.
• forcefields::Vector{String}: List of force field XML files.
• temp::Float64: Temperature in Kelvin.
• friction::Float64: Friction coefficient in 1/picosecond.
• step::Float64: Integration step size in picoseconds.
• steps::Int: Number of simulation steps.
• features: Which features to use for learning the chi function. - A vector of Int denotes the indices of all atoms to compute the pairwise distances from. - A vector of CartesianIndex{2} computes the specific distances between the atom pairs. - A number denotes the radius below which all pairs of atoms will be used (computed only on the starting configuration) - If nothing all pairwise distances are used.
• minimize::Bool: Whether to perform energy minimization on first state.
• nthreads: The number of threads to use for parallelization of multiple simulations.
• mmthreads: The number of threads to use for each OpenMM simulation. Set to "gpu" to use the GPU platform.

Returns

• OpenMMSimulation: An OpenMMSimulation object.

struct SimulationData{S,D,C,F}

A struct combining a simulation with the simulated coordinates and corresponding ISOKANN trainingsdata

Fields

• sim::S: The simulation object.
• data::D: The ISOKANN trainings data.
• coords::C: The orginal coordinates of the simulations.
• featurizer::F: A function mapping coordinates to ISOKANN features.

Iso(data; opt=NesterovRegularized(), model=defaultmodel(data), gpu=false, kwargs...)

Iso(sim::IsoSimulation; nx=100, nk=10, nd=1, kwargs...)

Convenience constructor which generates the SimulationData from the simulation sim and constructs the Iso object. See also Iso(data; kwargs...)

Arguments

• sim::IsoSimulation: The IsoSimulation object.
• nx::Int: The number of starting points.
• nk::Int: The number of koopman samples.
• nout::Int: Dimension of the χ function.

propagate(sim::OpenMMSimulation, x0::AbstractMatrix, nk)

Propagates nk replicas of the OpenMMSimulation sim from the inintial states x0.

Arguments

• sim: An OpenMMSimulation object.
• x0: Matrix containing the initial states as columns
• nk: The number of replicas to create.

run!(iso::Iso, n=1, epochs=1)

Run the training process for the Iso model.

Arguments

a iso::Iso: The Iso model to train.

• n::Int: The number of (outer) Koopman iterations.
• epochs::Int: The number of (inner) epochs to train the model for each Koopman evaluation.

Fully connected neural network with layers layers from n to nout dimensions. activation determines the activation function for each but the last layer lastactivation can be used to modify the last layers activation function

densenet(; layers::Vector{Int}, activation=Flux.sigmoid, lastactivation=identity, layernorm=true) -> Flux.Chain

Construct a fully connected neural network (Flux.Chain) with customizable layer sizes, activations, and optional input layer normalization.

Arguments

• layers::Vector{Int}: List of layer dimensions. For example, [10, 32, 16, 1] creates a network with input size 10, two hidden layers of size 32 and 16, and an output layer of size 1.
• activation: Activation function applied to all layers except the last. Defaults to Flux.sigmoid.
• lastactivation: Activation function for the final layer. Defaults to identity.
• layernorm::Bool: Whether to prepend a Flux.LayerNorm layer to normalize the input. Defaults to true.

Returns

A Flux.Chain composed of dense layers (and optionally a leading LayerNorm).

AdamRegularized(adam=1e-3, reg=1e-4)

Constructs an optimizer that combines weight decay regularization with ADAM. Uses reg for the weight decay parameter and lr as the learning rate for ADAM. Note that this is different from AdamW (Adam+WeightDecay) (c.f. Decay vs L2 Reg.).

NesterovRegularized(; lr=1e-3, reg=1e-4)

Constructs an optimizer that combines weight decay regularization with Nesterov momentum. Uses reg for the weight decay parameter and lr as the learning rate for Nesterov acceleration. This worked well as alternative where ADAM had problems.

Iso(data; opt=NesterovRegularized(), model=defaultmodel(data), gpu=false, kwargs...)

Iso(sim::IsoSimulation; nx=100, nk=10, nd=1, kwargs...)

Convenience constructor which generates the SimulationData from the simulation sim and constructs the Iso object. See also Iso(data; kwargs...)

Arguments

• sim::IsoSimulation: The IsoSimulation object.
• nx::Int: The number of starting points.
• nk::Int: The number of koopman samples.
• nout::Int: Dimension of the χ function.

SimulationData(sim::IsoSimulation, nx::Int, nk::Int; ...)
SimulationData(sim::IsoSimulation, xs::AbstractMatrix, nk::Int; ...)
SimulationData(sim::IsoSimulation, (xs,ys); ...)
SimulationData(xs, ys; pdb="", ...)  # for external simulation data

Generates SimulationData from a simulation with either

• nx initial points and nk Koopman samples
• xs as initial points and nk Koopman sample
• xs as inintial points and ys as Koopman samples
• xs and ys from external simulations
• xs a trajectory of an external simulation, implicitly computing ys via data_from_trajectory of succesive samples

AdamRegularized(adam=1e-3, reg=1e-4)

Constructs an optimizer that combines weight decay regularization with ADAM. Uses reg for the weight decay parameter and lr as the learning rate for ADAM. Note that this is different from AdamW (Adam+WeightDecay) (c.f. Decay vs L2 Reg.).

NesterovRegularized(; lr=1e-3, reg=1e-4)

Constructs an optimizer that combines weight decay regularization with Nesterov momentum. Uses reg for the weight decay parameter and lr as the learning rate for Nesterov acceleration. This worked well as alternative where ADAM had problems.

ca_rmsd(cainds, pdb="data/villin nowater.pdb", pdbref="data/villin/1yrf.pdb")

Returns a ReactionCoordsRMSD object which is used to calculate the Root Mean Square Deviation (RMSD) of the provided C-alpha atoms.

Inputs: - cainds: Indices of the C-alpha atoms to consider for the RMSD - target: PDB File containing the target structure to which the RMSD is computed - source: Alternative PDB File for the source coordinates in the case that the indices differ (i.e. when matching different topologies)

Example: rsmd = ca_rmsd(3:10, "data/villin/1yrf.pdb", "data/villin nowater.pdb") rmsd(rand(300,10))

data_from_trajectory(xs::AbstractMatrix; reverse=true, stride=1, lag=1)

Generate the (x,y) data pairs for ISOKANN from the trajectory xs.

stride controls the stride of the starting positions x and lag the lag for the end positions y in terms of trajectory frames. If reverse is true, construct also the time-reversed pairs (recomended for stable ISOKANN training).

laggedtrajectory(data::SimulationData, n) = laggedtrajectory(data.sim, n, x0=coords(data)[:, end])

Simulate a trajectory comprising of n simulations from the last point in data

load_trajectory(filename; top=nothing, kwargs...)

wrapper around Python's mdtraj.load(). Returns a (3 * natom, nframes) shaped array.

localpdistinds(coords::AbstractMatrix, radius)

Given coords of shape ( 3n x frames ) return the pairs of indices whose minimal distance along all frames is at least once lower then radius

mergedata(d1::SimulationData, d2::SimulationData)

Merge the data and features of d1 and d2, keeping the simulation and features of d1. Note that there is no check if simulation features agree.

Fully connected neural network with layers layers from n to nout dimensions. activation determines the activation function for each but the last layer lastactivation can be used to modify the last layers activation function

pdists(coords::AbstractArray, inds::Vector{<:Tuple})

Compute the pairwise distances between the particles specified by the tuples inds over all frames in traj. Assumes a column contains all 3n coordinates.

picking(X, n; dists = pairwise_one_to_many)

The picking algorithm, i.e. greedy farthest point sampling, for n points on the columns of X. A custom distance function (::Vector, ::Matrix)->(::Vector) may be passed through dists.

Returns X[:,qs], i.e. the picked samples, their former indices qs and their distances d to all other points.

rates(iso::Iso)

Return the coarse grained rate matrix Q satisfying Kχ = exp(τQ)χ

In the 1D ISOKANN case return the rates for χ and 1-χ.

reactionpath_minimum(iso::Iso, x0; steps=100)

Compute the reaction path by integrating ∇χ with orthogonal energy minimization.

Arguments

• iso::Iso: The isomer for which the reaction path minimum is to be computed.
• x0: The starting point for the reaction path computation.
• steps=100: The number of steps to take along the reaction path.

reactionpath_ode(iso, x0; steps=101, extrapolate=0, orth=0.01, solver=OrdinaryDiffEq.Tsit5(), dt=1e-3, kwargs...)

Compute the reaction path by integrating ∇χ as well as orth * F orthogonal to ∇χ where F is the original force field.

Arguments

• iso::Iso: The isomer for which the reaction path minimum is to be computed.
• x0: The starting point for the reaction path computation.
• steps=100: The number of steps to take along the reaction path.
• minimize=false: Whether to minimize the orthogonal to ∇χ before integration.
• extrapolate=0: How fast to extrapolate beyond χ 0 and 1.
• orth=0.01: The weight of the orthogonal force field.
• solver=OrdinaryDiffEq.Tsit5(): The solver to use for the ODE integration.
• dt=1e-3: The initial time step for the ODE integration.

reactive_path(xi::AbstractVector, coords::AbstractMatrix; sigma, minjump=0, maxjump=1, method=QuantilePath(0.05), normalize=false, sortincreasing=true)

Find the maximum likelihood path (under the model of brownion motion with noise sigma) through coords with times xi. Supports either CPU or GPU arrays.

Arguments

• coords: (ndim x npoints) matrix of coordinates.
• xi: time coordinate of the npoints points
• sigma: spatial noise strength of the model.
• minjump, maxjump: lower and upper bound to the jump in time xi along the path. Tighter bounds reduce the computational cost.
• method: either FromToPath, QuantilePath, FullPath or MaxPath, specifying the end points of the path
• normalize: whether to normalize all coords first
• sortincreasing: return the path from lower to higher xi values

resample_kde(xs, ys, n; kwargs...)

Return n indices of ys such that the corresponding points "fill the gaps" in the KDE of xs. For possible kwargs see kde_needles.

resample_kde(data::SimulationData, model, n; bandwith, unique)

add new samples to data by running new simulations starting at some ys (i.e. the propagated points) of data where these points are iteratively selected to be closest to the minimum of a KDE of the current chi values from xs. If unique is true, start simulations from point only where there were no simulations before. bandwith controls the bandwidth of the KDE.

restricted_localpdistinds(coords, radius, atoms)

Like localdists, but consider only the atoms with index in atoms

run!(iso::Iso, n=1, epochs=1)

Run the training process for the Iso model.

Arguments

a iso::Iso: The Iso model to train.

• n::Int: The number of (outer) Koopman iterations.
• epochs::Int: The number of (inner) epochs to train the model for each Koopman evaluation.

runadaptive!(iso; generations=1, nx=10, iter=100, cutoff=Inf)

Train iso with adaptive sampling. Sample kde new data points followed by iter isokann iterations and repeat this generations times. cutoff specifies the maximal data size, after which new data overwrites the oldest data. unique enforces resampling from yet unsampled ys only.

save_reactive_path(iso::Iso,
    coords::AbstractMatrix=coords(iso.data) |> cpu;
    sigma=1,
    maxjump=1,
    out="out/reactive_path.pdb",
    source=pdbfile(iso.data),
    kwargs...)

Extract and save the reactive path of a given iso.

Computes the maximum likelihood path with parameter sigma along the given data points, aligns it and saves it to the out path.

See also reactive_path.

Arguments

• iso::Iso: The Iso for which the reactive path is computed.

• coords: The coordinates amongst which the shortest path is calculated

• fullcoords: The full coordinates which are used to save the computed path frames (defaults to coords).

• weights: A vector of weights applied to each atom of the coordinates, determining the relative weight for the distance calculations.

• out="out/reactive_path.pdb": The output file path for saving the reactive path.

• source: The source .pdb file providing the topology

• kwargs...: additional parameters passed to reactive_path.

Returns

• ids: The IDs of the reactive path.

save_trajectory(filename, coords::AbstractMatrix; top::String)

save the trajectory given in coords to filename with the topology provided by the file top using mdtraj.

savecoords(path::String, iso::Iso, coords::AbstractMatrix=coords(iso.data); sorted=true, aligned=true)

Save the coordinates of the specified matrix of coordinates to a file, using the molecule in iso as a template. If sorted the sort the coordinates by their increasing χ value. If align then align each frame to the previous one.

trajectory(l::AbstractLangevin; T=lagtime(l), x0=randx0(l), save_start=false, saveat=lagtime(l), dt=dt(l)) generate a trajectory of length T, starting at x0 with stepsize dt, saving the output every saveat time.

trajectorydata_bursts(sim::IsoSimulation, steps, nk; kwargs...)

Simulate a single long trajectory of steps times the lagtime and start nk burst trajectories at each step for the Koopman samples.

x0–-x––x–- / | / | y y y y

trajectorydata_linear(sim::IsoSimulation, steps; reverse=false, kwargs...)

Simulate a single long trajectory of steps times the lagtime and use this "chain" to generate the corresponding ISOKANN data. If reverse is true, also add the time-reversed transitions

x (<)–> x (<)–> x

writechemfile(filename, data::Array{<:Any,2}; source)

Save the coordinates in data to filename with source as template using the Chemfiles library

FeaturesPairs(pairs::Vector{Tuple{Int,Int}})
FeaturesPairs(system; selector="all", maxdist=Inf, maxfeatures=Inf)

Creates a FeaturesPairs object from either:

• a list of index pairs (Vector{Tuple{Int,Int}}) passed directly.
• an OpenMMSimulation or PDB file path (String), selecting atom pairs using MDTraj selector syntax (selector), optionally filtered by maxdist (in nm) and limited to maxfeatures randomly sampled pairs.

propagate(sim::OpenMMSimulation, x0::AbstractMatrix, nk)

Propagates nk replicas of the OpenMMSimulation sim from the inintial states x0.

Arguments

• sim: An OpenMMSimulation object.
• x0: Matrix containing the initial states as columns
• nk: The number of replicas to create.

savecoords(path, sim::OpenMMSimulation, coords::AbstractArray{T})

Save the given coordinates in a .pdb file using OpenMM

trajectory(sim::OpenMMSimulation{Nothing}, steps=steps(sim); saveevery=1, x0=coords(sim), resample_velocities=false, throw=false, showprogress=true, reclaim=true)

Simulates the trajectory of an OpenMM simulation.

Arguments

• sim::OpenMMSimulation{Nothing}: The OpenMM simulation object.
• steps: The number of steps to simulate. Defaults to the number of steps defined in the simulation object.
• saveevery: Interval at which to save the trajectory. Defaults to 1.
• x0: Initial coordinates for the simulation. Defaults to the current coordinates of the simulation object.
• sample_velocities: Whether to sample velocities at the start of the simulation.
• resample_velocities: Whether to resample velocities after each saveevery steps. Defaults to false.
• throw: Whether to throw an error if the simulation fails. If false it returns the trajectory computed so far. Defaults to false.
• showprogress: Whether to display a progress bar during the simulation. Defaults to true.
• reclaim: Whether to reclaim CUDA memory before the simulation. Defaults to true.

Returns

• The trajectory of the simulation as a matrix of coordinates.

DataTuple = Tuple{Matrix{T},Array{T,3}} where {T<:Number}

We represent data as a tuple of xs and ys.

xs is a matrix of size (d, n) where d is the dimension of the system and n the number of samples. ys is a tensor of size (d, k, n) where k is the number of koopman samples.

abstract type IsoSimulation

Abstract type representing an IsoSimulation. Should implement the methods coords, propagate, dim

LazyTrajectory(path::String)

Represents the trajectory path as matrix whose columns are lazily loaded from disk.

struct ReactionCoordsRMSD

Instances of this object allow to compute the Root Mean Square Deviation (RMSD) to a part of a reference molecule. See also ca_rmsd.

TransformStabilize(target, last=nothing)

Wraps another target and permutes its target to match the previous target

Currently we also have the stablilization (wrt to the model though) inside most Transforms. TODO: Decide which to keep

TransformGramSchmidt()

Compute the target through a Gram-Schmidt orthonormalisation.

TransformISA(permute)

Compute the target via the inner simplex algorithm (without feasiblization routine). permute specifies whether to apply the stabilizing permutation

TransformPseudoInv(normalize, direct, eigenvecs, permute)

Compute the target by approximately inverting the action of K with the Moore-Penrose pseudoinverse.

If direct==true solve chi * pinv(K(chi)), otherwise inv(K(chi) * pinv(chi))). eigenvecs specifies whether to use the eigenvectors of the schur matrix. normalize specifies whether to renormalize the resulting target vectors. permute specifies whether to permute the target for stability.

TransformShiftscale()

Classical 1D shift-scale (ISOKANN 1)

scales with n=length(hs)

bootstrap(sim, nx, ny)

compute initial data by propagating the molecules initial state to obtain the xs and propagating them further for the ys

centercoords any given states by shifting their individual 3d mean to the origin

compute the chi exit rate as per Ernst, Weber (2017), chap. 3.3

constrained_free_energy(iso, xs; sim, steps)

Compute the free energy using Thermodynamic Integration. Starting from the levelset samples xs orhtogonal simulations estimate the mean force along χ, which is integrated to yield the PMF.

Arguments

iso the Iso object. xs the starting points (which should be well distributed in state space). sim the simulation used for the orthongal sampling. steps the number of steps in each orthogonal simulation.

Returns

F the free energy energy surface of χ in kJ/mol up to an additive constant.

delta_G(PMF,chi_vals)

Convenience function to compute free energy differences in a double well free energy surface.

densenet(; layers::Vector{Int}, activation=Flux.sigmoid, lastactivation=identity, layernorm=true) -> Flux.Chain

Construct a fully connected neural network (Flux.Chain) with customizable layer sizes, activations, and optional input layer normalization.

Arguments

• layers::Vector{Int}: List of layer dimensions. For example, [10, 32, 16, 1] creates a network with input size 10, two hidden layers of size 32 and 16, and an output layer of size 1.
• activation: Activation function applied to all layers except the last. Defaults to Flux.sigmoid.
• lastactivation: Activation function for the final layer. Defaults to identity.
• layernorm::Bool: Whether to prepend a Flux.LayerNorm layer to normalize the input. Defaults to true.

Returns

A Flux.Chain composed of dense layers (and optionally a leading LayerNorm).

energyminimization_chilevel(iso, x0; f_reltol=1e-3, alphaguess=1e-5, iterations=20, show_trace=false, skipwater=false, algorithm=Optim.GradientDescent, xtol=nothing)

Local energy minimization on the current levelset of the chi function

expectation(f, xs)

Computes the expectation value of f over xs. Supports WeightedSamples through extra dispatch

exportdata(data::AbstractArray, model, sys, path="out/data.pdb")

Export data to a PDB file.

This function takes an AbstractArray data, sorts it according to the model evaluation, removes duplicates, transforms it to standard form and saves it as a PDB file to path.

fixperm(new, old)

Permutes the rows of new such as to minimize L1 distance to old.

Arguments

• new: The data to match to the reference data.
• old: The reference data.

flatpairdists(x)

Assumes each col of x to be a flattened representation of multiple 3d coords. Returns the flattened pairwise distances as columns.

collapse the first and second dimension of the array A into the first dimension

flattenlast(x)

Concatenate all but the first dimension of x. Usefull to convert a tensor of samples into a matrix

Given a model and return a copy with its last layer replaced with given output dimension n

obtain the input dimension of a Flux model

integratechi(f, chivals)

Cumulative integral of the mean force with respect to χ using the trapezoid rule.

Arguments

f The mean force. chi_vals The levelset χ values.

Returns

F the (rigid) free energy surface of χ.

load(path::String, iso::Iso)

Load the Iso object from a JLD2 file Note that it will be loaded to the CPU, even if it was saved on the GPU. An OpenMMSimulation will be reconstructed anew from the saved pdb file.

localmeanforce(iso, xs; sim, steps)

Compute the free energy using Thermodynamic Integration. Bins the samples into levelsets, computes the mean force along χ locally in every levelset. (Extremely extensive sampling necessary.)

Arguments

iso the Iso object. xs the starting points (which should be well distributed in state space). nbins The number of bins/levelsets.

Returns

F the free energy surface of χ in kJ/mol up to an additive constant.

 marginal_free_energy(iso::Iso;nbins)

Compute the free energy from the density of chi values.

Arguments

iso the Iso object. nbins the number of bins of the histogram used for estimation.

Returns

F the free energy energy surface of χ in kJ/mol up to an additive constant.

convenience wrapper returning the provided model with the default AdamW optimiser

Obtain the output dimension of a Flux model

pairdistfeatures(inds::AbstractVector)

Returns a featurizer function which computes the pairwise distances between the particles specified by inds

pickclosest(haystack, needles)

Return the indices into haystack which lie closest to needles without duplicates by removing haystack candidates after a match. Note that this is not invariant under pertubations of needles

scales with n log(n) m where n=length(haystack), m=length(needles)

picking_aligned(x::AbstractMatrix, m::Integer)

The picking algorithm using pairwise aligned distances, e.g. for molecular coordinates. Assumes the columnes of x to be vectors of size 3xN holding the cartesian coordinates.

see also picking()

scatter plot of all first "O" atoms of the starting points xs as well as the "O" atoms from the koopman samples to the first point from ys

combined plot of loss graph and current atoms

reactionforce(iso, sim, x, direction, orth=1)

Compute the vector f with colinear component to dchi/dx such that dchi/dx * f = 1 and orth*forcefield in the orthogonal space

computation of the real invariant subspace from complex eigenvectors

resample_kde_ash(xs, ys, n=10; m=20, target=Distributions.Uniform())

Pick n indices iys of the samples ys such that [xs; ys[iys]] approximates the target distribution. Approximates the density through an AverageShiftedHistogram on the range rng with window width m.

resample_strat(d::SimulationData, model, n)

χ-stratified subsampling. Select n samples amongst the provided ys/koopman points of d such that their χ-value according to model is approximately uniformly distributed and propagate them. Returns a new SimulationData which has the new data appended.

function sample_coords(iso,n_points;xs)

Convenience function to uniformly sample npoints out of the χ distribution of xs coordinates.

save(path::String, iso::Iso)

Save the complete Iso object to a JLD2 file

saveextrema(path::String, iso::Iso)

Save the two extermal configurations (metastabilities) to the file path.

fixed point plot, i.e. x vs model(x)

given an array of arbitrary shape, select the rows inds in the first dimension

simulationtime(iso::Iso)

print and return the total simulation time contained in the data of iso in nanoseconds.

sqpairdist(x::AbstractArray)

Compute the squared pairwise distances between the columns of x. If x has 3 dimensions, the computation is batched along the 3rd dimension.

subsample(model, data::Array, n) :: Matrix
subsample(model, data::Tuple, n) :: Tuple

Subsample n points of data uniformly in model. If model returns multiple values per sample, subsample along each dimension.

subsample(data, nx)

return a random subsample of nx points from data

subsample_inds(model, xs, n) :: Vector{Int}

Returns n indices of xs such that model(xs[inds]) is approximately uniformly distributed.

subsbample_uniformgrid(ys, n) -> is

given a list of values ys, return nindicesissuch thatys[is]` are approximately uniform by picking the closest points to a randomly perturbed grid in [0,1].

updatehistory(L, l)

Insert the newest observations l of size (n,d) into columns 2:d+1 of the history matrix L of size (n,h). Automatically grow the n dimension of L if l is bigger.

Pairwise distances between all provided atoms

featurepairs(d::ISOKANN.SimulationData)

returns pairs of atom indices corresponding to the pairwise distance features

integrate_girsanov(sim::OpenMMSimulation; x0=coords(sim), steps=steps(sim), bias, reclaim=true)

Integrate overdamped Langevin dynamics while computing Girsanov weights for a given bias:

Xt = F dt + σ u dt + sigma dBt
where F = force(sim, Xt), σ = sqrt(2 kB T / (γ M)), u = bias(Xt; t, sigma, F)

Arguments

• sim::OpenMMSimulation : The simulation object.
• x0 : Initial coordinates (default: coords(sim)).
• steps : Number of integration steps (default: steps(sim)).
• bias : Function x -> u that returns the bias u for the Girsanov reweighting (note that the actual force perturbation is sigma * u).
• reclaim : Whether to preallocate memory for performance.

Returns

• x : Final coordinates.
• g : Cumulative Girsanov weight (logarithmic).
• z : Full trajectory of positions.

Notes

The method assumes overdamped Langevin dynamics and accumulates the corresponding Girsanov weight g.

integrate_langevin(sim::OpenMMSimulation, x0=coords(sim); steps=steps(sim), perturbation!::Union{Function,Nothing}=nothing, saveevery::Union{Int,Nothing}=nothing, reclaim=true)

Perform standard Langevin dynamics integration for a given system.

Arguments

• sim::OpenMMSimulation : The simulation object with system parameters.
• x0 : Initial coordinates (default: coords(sim)).
• steps : Number of integration steps (default: steps(sim)).
• perturbation! : Optional function (F, x) -> nothing to perturb the force F in-place.
• saveevery : Interval to save coordinates (default: saves only final state).
• reclaim : Whether to preallocate/reclaim memory for performance.

Returns

• Array of positions: either the final positions or a trajectory if saveevery is specified.

Notes

Uses a simple Euler-Maruyama step for (underdamped) Langevin dynamics with optional bias.

langevin_girsanov!(sim::OpenMMSimulation, steps=steps(sim); bias=sim.bias, saveevery=1, x0=coords(sim), resample_velocities=false, showprogress=true, throw=true, reclaim=true)

Perform underdamped Langevin dynamics using the ABOBA integrator with Girsanov reweighting.

Arguments

• sim::OpenMMSimulation : Simulation object containing masses, friction, forces, and temperature.
• steps : Number of integration steps (default: steps(sim)).
• bias : Function q; t, sigma, F -> B returning perturbation force for reweighting.
• saveevery : Interval to save coordinates and weights (default: 1).
• x0 : Initial coordinates (default: coords(sim)).
• resample_velocities : If true, resample momenta when saving (invalidates Girsanov weights).
• showprogress : Whether to show a progress bar.
• throw : Whether to throw errors on invalid states.
• reclaim : Whether to preallocate memory for performance.
• sigmascaled : Whether the bias is scaled by the noise amplitude, i.e. whether the resulting force is σ * bias (true) or only bias (false) (default: true).

Returns

• WeightedSamples containing:
  • qs : Saved coordinates at intervals of saveevery.
  • weights : Corresponding Girsanov weights exp(-g).

Notes

• Uses the ABOBA splitting scheme for stable underdamped Langevin integration.
• Computes Girsanov weights to account for the applied bias.

laggedtrajectory(sim::OpenMMSimulation, lags; steps=steps(sim), resample_velocities=true, kwargs...)

Generate a lagged trajectory for a given OpenMMSimulation. E.g. x0–x–x–x for lags=3 and steps=2

Arguments

• sim::OpenMMSimulation: The simulation object.
• lags: The number of steps.
• steps: The lagtime, i.e. number of steps to take in the simulation.
• resample_velocities: Whether to resample velocities according to Maxwell-Boltzman for each lag.
• kwargs...: Additional keyword arguments to pass to the trajectory function.

Returns

• A matrix of lags samples which each have steps simulation-steps inbetween them.

generate n random inintial points for the simulation mm