Langevin Dynamics In Hamiltonian Phase Space

Langevin dynamics is usually written like

M \ddot{\vec{X}} = - \nabla_{X} U (X) - γ M \dot{\vec{X}} + \sqrt{2 γ M k_{B} T} d \vec{W}

which is really just two first order differential form relations

d {\vec{x}}_{i} = {\vec{v}}_{i} d t

M_{i} d {\vec{v}}_{i} = - \frac{\partial}{\partial x_{i}} U (\vec{x}) d t - γ M_{i} {\vec{v}}_{i} d t + \sqrt{2 γ M_{i} k_{B} T} d {\vec{W}}_{i}

but statistical mechanics is formulated in Hamiltonian phase space.

Reformulating in Hamiltonian Phase space

Instead if we can generalize mass then can redefine the process as

d {\vec{q}}_{i} = \frac{\partial H}{\partial p_{i}} d t

d {\vec{p}}_{i} = - \frac{\partial H}{\partial q_{i}} d t - γ_{i} {\vec{p}}_{i} d t + K_{i} \sqrt{2 γ_{i} k_{B} T} d {\vec{W}}_{i}

where we will be solving for $K_{i}$ to converge to the canonical ensemble equilibrium distribution.

Drift Diffusion Process

We can then combine both equations to a single drift diffusion process and generalize the matrix so the momentum can share noise

d [\begin{matrix} q_{1} \\ ⋮ \\ q_{n} \\ p_{1} \\ ⋮ \\ p_{n} \end{matrix}] = [\begin{matrix} \frac{\partial H}{\partial p_{1}} \\ ⋮ \\ \frac{\partial H}{\partial p_{n}} \\ - \frac{\partial H}{\partial q_{1}} - γ_{1} p_{1} \\ ⋮ \\ - \frac{\partial H}{\partial q_{n}} - γ_{n} p_{n} \end{matrix}] d t + [\begin{matrix} 0 & \dots & 0 & \dots & 0 \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ 0 & \dots & 0 & \dots & 0 \\ 0 & \dots & K_{1, 1} \sqrt{2 γ_{1} k_{B} T} & \dots & K_{1, n} \sqrt{2 γ_{1} k_{B} T} \\ ⋮ & ⋮ & ⋮ & ⋱ & ⋮ \\ 0 & \dots & K_{n, 1} \sqrt{2 γ_{n} k_{B} T} & \dots & K_{n, n} \sqrt{2 γ_{n} k_{B} T} \end{matrix}] d \vec{W}

Fokker Plank

From the Fokker Plank equation any drift diffusion process

d \vec{z} = \vec{A} (\vec{z}) d t + 𝐁 (\vec{z}) d \vec{W}

The probability density at each point will evolve as

\frac{\partial ρ}{\partial t} = - \sum_{i} \frac{\partial}{\partial z_{i}} (ρ {\vec{A}}_{i}) + \frac{1}{2} \sum_{i} \sum_{j} \frac{\partial^{2}}{\partial z_{i} \partial z_{j}} ((𝐁^{T} 𝐁)_{i, j} ρ)

Substitute in the Langevin equation.

\frac{\partial ρ}{\partial t} = - \sum_{i = 1}^{n} \frac{\partial}{\partial q_{i}} (ρ \frac{\partial H}{\partial p_{i}}) - \frac{\partial}{\partial p_{i}} (ρ (\frac{\partial H}{\partial q_{i}} + γ_{i} p_{i})) + \frac{1}{2} \sum_{i = 1}^{n} \sum_{j = 1}^{n} \frac{\partial^{2}}{\partial p_{i} \partial p_{j}} (ρ 2 k_{B} T \sum_{k} K_{i, k} K_{j, k} \sqrt{γ_{i} γ_{j}})

We can combine the two sums and simplify.

\sum_{i = 1}^{n} - \frac{\partial}{\partial q_{i}} (ρ \frac{\partial H}{\partial p_{i}}) + \frac{\partial}{\partial p_{i}} (ρ (\frac{\partial H}{\partial q_{i}} + γ_{i} p_{i})) + k_{B} T \sum_{j = 1}^{n} \frac{\partial^{2}}{\partial p_{i} \partial p_{j}} (ρ \sum_{k} K_{i, k} K_{j, k} \sqrt{γ_{i} γ_{j}})

Let $F_{i, j} = \sum_{k} K_{i, k} K_{j, k} \sqrt{γ_{i} γ_{j}}$ so that

\sum_{i = 1}^{n} - \frac{\partial}{\partial q_{i}} (ρ \frac{\partial H}{\partial p_{i}}) + \frac{\partial}{\partial p_{i}} (ρ (\frac{\partial H}{\partial q_{i}} + γ_{i} p_{i})) + k_{B} T \sum_{j = 1}^{n} \frac{\partial^{2}}{\partial p_{i} \partial p_{j}} (ρ F_{i, j})

and use the chain rule to expand the derivatives in the first part.

= \sum_{i = 1}^{n} \frac{\partial H}{\partial q_{i}} \frac{\partial ρ}{\partial p_{i}} - \frac{\partial H}{\partial p_{i}} \frac{\partial ρ}{\partial q_{i}} + \frac{\partial ρ}{\partial p_{i}} γ_{i} p_{i} + ρ γ_{i} + k_{B} T \sum_{j = 1}^{n} \frac{\partial^{2}}{\partial p_{i} \partial p_{j}} (ρ F_{i, j})

Then for the second

= \sum_{i = 1}^{n} \frac{\partial H}{\partial q_{i}} \frac{\partial ρ}{\partial p_{i}} - \frac{\partial H}{\partial p_{i}} \frac{\partial ρ}{\partial q_{i}} + \frac{\partial ρ}{\partial p_{i}} γ_{i} p_{i} + ρ γ_{i} + k_{B} T \sum_{j = 1}^{n} \frac{\partial}{\partial p_{i}} (\frac{\partial ρ}{\partial p_{i}} F_{i, j} + \frac{\partial F_{i, j}}{\partial p_{i}} ρ)

and again.

= \sum_{i = 1}^{n} \frac{\partial H}{\partial q_{i}} \frac{\partial ρ}{\partial p_{i}} - \frac{\partial H}{\partial p_{i}} \frac{\partial ρ}{\partial q_{i}} + \frac{\partial ρ}{\partial p_{i}} γ_{i} p_{i} + ρ γ_{i} + k_{B} T \sum_{j = 1}^{n} (\frac{\partial^{2} ρ}{\partial p_{i} \partial p_{j}} F_{i, j} + ρ \frac{\partial^{2} F_{i, j}}{\partial p_{i} \partial p_{j}} + \frac{\partial ρ}{\partial p_{i}} \frac{\partial F_{i, j}}{\partial p_{j}} + \frac{\partial ρ}{\partial p_{j}} \frac{\partial F_{i, j}}{\partial p_{i}})

Enforcing canonical distribution convergence

In the canonical ensemble it should converge to the below.

ρ_{0} (\vec{q}, \vec{p}) = \frac{1}{Z} e^{- \frac{H (\vec{q}, \vec{p})}{k_{B} T}}

So it should be a stable distribution by definition of equilibrium

\frac{\partial ρ_{0}}{\partial t} = 0

and putting it into the Fokker Plank equation we derived

= \frac{1}{Z} \sum_{i = 1}^{n} \frac{\partial H}{\partial q_{i}} \frac{\partial e^{- \frac{H (\vec{q}, \vec{p})}{k_{B} T}}}{\partial p_{i}} - \frac{\partial H}{\partial p_{i}} \frac{\partial e^{- \frac{H (\vec{q}, \vec{p})}{k_{B} T}}}{\partial q_{i}} + \frac{\partial e^{- \frac{H (\vec{q}, \vec{p})}{k_{B} T}}}{\partial p_{i}} γ_{i} p_{i} + e^{- \frac{H (\vec{q}, \vec{p})}{k_{B} T}} γ_{i}

+ k_{B} T \sum_{j = 1}^{n} (\frac{\partial^{2} e^{- \frac{H (\vec{q}, \vec{p})}{k_{B} T}}}{\partial p_{i} \partial p_{j}} F_{i, j} + e^{- \frac{H (\vec{q}, \vec{p})}{k_{B} T}} \frac{\partial^{2} F_{i, j}}{\partial p_{i} \partial p_{j}} + \frac{\partial e^{- \frac{H (\vec{q}, \vec{p})}{k_{B} T}}}{\partial p_{i}} \frac{\partial F_{i, j}}{\partial p_{j}} + \frac{\partial e^{- \frac{H (\vec{q}, \vec{p})}{k_{B} T}}}{\partial p_{j}} \frac{\partial F_{i, j}}{\partial p_{i}})

the Possion bracket part goes to zero.

= \frac{1}{Z} \sum_{i = 1}^{n} \frac{\partial e^{- \frac{H (\vec{q}, \vec{p})}{k_{B} T}}}{\partial p_{i}} γ_{i} p_{i} + e^{- \frac{H (\vec{q}, \vec{p})}{k_{B} T}} γ_{i}

+ k_{B} T \sum_{j = 1}^{n} (\frac{\partial^{2} e^{- \frac{H (\vec{q}, \vec{p})}{k_{B} T}}}{\partial p_{i} \partial p_{j}} F_{i, j} + e^{- \frac{H (\vec{q}, \vec{p})}{k_{B} T}} \frac{\partial^{2} F_{i, j}}{\partial p_{i} \partial p_{j}} + \frac{\partial e^{- \frac{H (\vec{q}, \vec{p})}{k_{B} T}}}{\partial p_{i}} \frac{\partial F_{i, j}}{\partial p_{j}} + \frac{\partial e^{- \frac{H (\vec{q}, \vec{p})}{k_{B} T}}}{\partial p_{j}} \frac{\partial F_{i, j}}{\partial p_{i}})

Then do the chain rule a bunch.

= \frac{1}{Z} \sum_{i = 1}^{n} - \frac{e^{- \frac{H}{k_{B} T}}}{k_{B} T} \frac{\partial H}{\partial p_{i}} γ_{i} p_{i} + e^{- \frac{H}{k_{B} T}} γ_{i}

+ k_{B} T \sum_{j = 1}^{n} ((\frac{e^{- \frac{H}{k_{B} T}}}{(k_{B} T)^{2}} \frac{\partial H}{\partial p_{i}} \frac{\partial H}{\partial p_{j}} - \frac{e^{- \frac{H}{k_{B} T}}}{k_{B} T} (\frac{\partial^{2} H}{\partial p_{i} \partial p_{j}})) F_{i, j}

+ e^{- \frac{H}{k_{B} T}} \frac{\partial^{2} F_{i, j}}{\partial p_{i} \partial p_{j}}

- \frac{e^{- \frac{H}{k_{B} T}}}{k_{B} T} (\frac{\partial H}{\partial p_{i}} \frac{\partial F_{i, j}}{\partial p_{j}} + \frac{\partial H}{\partial p_{j}} \frac{\partial F_{i, j}}{\partial p_{i}}))

Factor out the probability.

= p_{0} \sum_{i = 1}^{n} - \frac{1}{k_{B} T} \frac{\partial H}{\partial p_{i}} γ_{i} p_{i} + γ_{i}

+ k_{B} T \sum_{j = 1}^{n} ((\frac{1}{(k_{B} T)^{2}} \frac{\partial H}{\partial p_{i}} \frac{\partial H}{\partial p_{j}} - \frac{1}{k_{B} T} (\frac{\partial^{2} H}{\partial p_{i} \partial p_{j}})) F_{i, j}

+ \frac{\partial^{2} F_{i, j}}{\partial p_{i} \partial p_{j}}

- \frac{1}{k_{B} T} (\frac{\partial H}{\partial p_{i}} \frac{\partial F_{i, j}}{\partial p_{j}} + \frac{\partial H}{\partial p_{j}} \frac{\partial F_{i, j}}{\partial p_{i}}))

Multiply through the temperature.

= p_{0} \sum_{i = 1}^{n} - \frac{1}{k_{B} T} \frac{\partial H}{\partial p_{i}} γ_{i} p_{i} + γ_{i}

+ \sum_{j = 1}^{n} (\frac{1}{k_{B} T} \frac{\partial H}{\partial p_{i}} \frac{\partial H}{\partial p_{j}} - (\frac{\partial^{2} H}{\partial p_{i} \partial p_{j}})) F_{i, j}

+ k_{B} T \frac{\partial^{2} F_{i, j}}{\partial p_{i} \partial p_{j}}

- (\frac{\partial H}{\partial p_{i}} \frac{\partial F_{i, j}}{\partial p_{j}} + \frac{\partial H}{\partial p_{j}} \frac{\partial F_{i, j}}{\partial p_{i}})

Then group based on the temperature.

0 = \sum_{i = 1}^{n}

\frac{1}{k_{B} T} (- \frac{\partial H}{\partial p_{i}} γ_{i} p_{i} + \sum_{j = 1}^{n} \frac{\partial H}{\partial p_{i}} \frac{\partial H}{\partial p_{j}} F_{i, j})

+ k_{B} T \sum_{j = 1}^{n} \frac{\partial^{2} F_{i, j}}{\partial p_{i} \partial p_{j}}

+ γ_{i} + \sum_{j = 1}^{n} (- \frac{\partial H}{\partial p_{i}} \frac{\partial F_{i, j}}{\partial p_{j}} - \frac{\partial H}{\partial p_{j}} \frac{\partial F_{i, j}}{\partial p_{i}} - (\frac{\partial^{2} H}{\partial p_{i} \partial p_{j}}) F_{i, j})

Any Hamiltonian and matrix $𝐅 = 𝐊^{⊤} 𝐊$ that satisfy the above equation will correctly have equilbibrium as a stable distribution

Quadratic Momentum

Lets assume that each term of the sum should be zero. Since it should be valid for any temperature there are 3 equalities that must hold

\frac{\partial H}{\partial p_{i}} γ_{i} p_{i} = \sum_{j = 1}^{n} \frac{\partial H}{\partial p_{i}} \frac{\partial H}{\partial p_{j}} F_{i, j}

0 = \frac{\partial^{2} F_{i, j}}{\partial p_{i} \partial p_{j}}

γ_{i} = \sum_{j = 1}^{n} (\frac{\partial H}{\partial p_{i}} \frac{\partial F_{i, j}}{\partial p_{j}} + \frac{\partial H}{\partial p_{j}} \frac{\partial F_{i, j}}{\partial p_{i}} + (\frac{\partial^{2} H}{\partial p_{i} \partial p_{j}}) F_{i, j})

From here or here a somewhat general Hamiltonian form is

H (q, p) = \frac{1}{2} p^{⊤} 𝐌 (q) p + U (q)

where the "mass matrix" is only a function of position, then one can define a matrix

𝐒 = \frac{1}{2} (𝐌 + 𝐌^{⊤})

where

\frac{\partial^{2} H}{\partial p_{i} \partial p_{j}} = S_{i, j} = S_{j, i}

Then if we assume that the $F$ 's are independent of momentum then the third equality becomes

γ_{i} = \sum_{j = 1}^{n} (\frac{\partial^{2} H}{\partial p_{i} \partial p_{j}}) F_{i, j}

γ_{i} = \sum_{j = 1}^{n} F_{i, j} S_{j, i} = (F S)_{i, j}

if all the gammas are equal then

𝐅 = γ 𝐒^{- 1}

F_{i, j} = (γ S^{- 1})_{i, j}

\sum_{k} γ K_{i, k} K_{j, k} = (γ S^{- 1})_{i, j}

\sum_{k} K_{i, k} K_{j, k} = (S^{- 1})_{i, j}

Which requires

𝐊^{⊤} 𝐊 = 𝐒^{- 1}

since $S$ is symmetric so it has an orthonormal eigenbasis

𝐒 = 𝐐^{⊤} 𝐃 𝐐

and since we assumed it to be invertable then the inverse is

𝐒^{- 1} = 𝐐^{⊤} 𝐃^{- 1} 𝐐

and also assume it's positive definite then all the eigenvalues are positive so one can take the square root

𝐒^{-} 1 = 𝐐^{⊤} 𝐃^{- 1 / 2} 𝐃^{- 1 / 2} 𝐐

𝐒^{-} 1 = (𝐃^{- 1 / 2} 𝐐)^{⊤} 𝐃^{- 1 / 2} 𝐐

So one possible solution is

𝐊 = 𝐃^{- 1 / 2} 𝐐