Canonical ensemble by Generalized Langevin Dynamics

In order to simulate the canonical ensemble we use a Langevin dynamics and we assume that our system is coupled with a thermal bath due to Quantum Monte Carlo statistical noise plus an additional friction term:

$$\left\{ \begin{array}{l} \dot{v}_i(t)= -\gamma_{ij}(x) v_j(t) +\frac{f_i(x(t))}{m_i} + \Gamma_i(t)\\ \dot{x}_i(t)= v_i(i) \end{array} \right.$$ (5.4)

with

$$\langle \Gamma_i(t) \Gamma_j(t') \rangle = \alpha_{ij}(x) \delta(t-t')$$ (5.5)

$$\langle \Gamma_i(t)\rangle = 0.$$ (5.6)

where $\Gamma$ is a generic Gaussian noise that implicitly contains the one associated to the statistical evaluation of the forces by QMC. Notice that in this case an explicit dependence on $x$ of the noise has been taken in account. This is a realistic case that has not been considered so far not even in (89). In the following we determinate a form for the friction matrix $\gamma_{ij}(x)$ that allows to converge the usual Boltzmann distribution at a given temperature.

To this purpose we write down the corresponding Fokker-Planck equation. Following (90) we have to evaluate the drift and the diffusion coefficient:

$$D_i(x,t) = \lim_{\tau \to 0} \frac{1}{\tau} \langle x_i(t+\tau) -x_i(t) \rangle$$ (5.7)

$$D_{ij}(x,t) = \frac{1}{2}\lim_{\tau \to 0} \frac{1}{\tau} \langle [x_i(t+\tau) -x_i(t)][x_j(t+\tau) -x_j(t)] \rangle$$ (5.8)

A straightforward calculation shows that:

$$D_{x_i}(t) = v_i(t)$$ (5.9)

$$D_{v_i}(t) = -\gamma_{ij}(x)v_j(t)+\frac{f_i(t)}{m_i}$$ (5.10)

$$D_{x_i,x_j} = D_{x_i,v_j}=0$$ (5.11)

$$D_{v_i,v_j} = \frac{\alpha_{ij}(x)}{2}$$ (5.12)

And so the corresponding Fokker-Plank equation will be:

$$\frac{\partial W(x,v,t)}{\partial t} = \sum_i \left \{ \frac{\par... ...rac{\alpha_{ij}(x)}{2} \frac{\partial}{\partial v_j} \right ]\right \} W(x,v,t)$$ (5.13)

Then the friction matrix $\hat{\gamma}$ is chosen in a way that the stationary solution of the Fokker-Planck equation is the canonical distribution:

$$p_{eq}(v_1,...,v_n,x_1,...,x_n) \simeq e^{-\beta H}.$$ (5.14)

More precisely by substituting the Boltzmann distribution

$$W_{eq}(x,v) = e^{\frac{ -\sum_i\frac{m_i v_i^2}{2}-V(x)}{KT}}$$ (5.15)

in the equation 5.13 we obtain:

$$\gamma_{ij}(x) = \frac{\alpha_{ij}}{2}(x) \beta m_j$$ (5.16)

So for a given noise on the forces $\alpha_{ij}$ and the desired temperature we can set the friction tensor using eq. 5.16 in order to obtain the Boltzmann distribution.
Notice that the external random noise scale as $\sqrt \tau$ (see also Ref. (90)) whereas the QMC noise scales as $\tau$ . Therefore in the limit $\tau \rightarrow 0$ if we do not add any external noise and we set $\hat{\gamma}$ according to eq. 5.16, the system will converge to the Newton dynamics at zero temperature. Nevertheless adding an appropriate external noise it is possible equilibrate the system to the canonical ensemble at the desired temperature.
A peculiar feature of this approach is that in the limit of small $\tau$ the statistical error on the forces becomes irrelevant because $\tau \leq \sqrt \tau$ .
The stability of this approach compared with the first order Langevin Dynamic is now evident, in fact there is no need to calculate the inverse of the covariance matrix. Moreover in the second order Langevin Dynamics the temperature can be estimated at posteriori by equality:

$$\frac{3}{2}K_B T = \frac{1}{2} M \langle V^2 \rangle,$$ (5.17)

compensating the error in the integration of the GLE.