# Simulation of Nuclear Quantum Effects in Condensed Matter Systems via Quantum Baths
## Appendix A Colored Noise Generation  
While in the QT formalism, by [PERSON] et al., the GLE arises from a Markovian dynamics in an extended phase space, the Quantum Thermal Bath (QTB) explicitly relies on a non-Markovian random force, and requires the numerical generation of a colored noise, that is a _time-correlated_ noise, to mimic the quantum frequency-dependent energy distribution. 1. A priori noise generation  
In the original paper [12] as in the former work by [PERSON] [64], the noise is generated in the frequency domain: the real and imaginary parts of the discrete Fourier transform components of \(\bar{R}(\omega)\) are drawn as Gaussian random numbers with an amplitude corresponding to Equation (14). Backward Fourier transform then provides the force in the time domain \(R(t)\) that is stored in a large file before starting the actual simulation. During the simulation, the numbers are then simply read as needed. This method is well suited for _ab initio_ molecular dynamics [1; 53], where electronic calculations are computationally expensive so that the number of atoms and the simulation length remain relatively small. However, when simulating long trajectories with quickly evaluated analytical interatomic forces fields, storage and access to very large random force fields can become problematic. Furthermore, this _a priori_ generation method for the colored random force is not suited to adQTB simulations where the memory kernel is adapted on-the-fly during the molecular dynamics in order to compensate for ZPE leakage. 2. _On-the-fly_ noise generation  
An alternative, _on-the-fly_ noise-generation method was proposed [39], based on the observation that the QTB colored noise can be obtained by filtering a white noise [87] with the filter defined by:  
\[\bar{H}(\omega)=\sqrt{2m\gamma\theta(\omega,T)} \tag{11}\]  
The random force \(R(t)\) can then be obtained _on-the-fly_ in real time by convolution of the deterministic function \(H(t)\) (the inverse Fourier transform of \(\bar{H}(\omega)\)) with a Gaussian normalized white noise \(r(t)\) (with a spectral density \(\bar{C}_{rr}(\omega)=1\)). Discretization yields:  
\[\bar{H}_{k}=\bar{H}(k\,\delta\omega),\;k\in[-n,n-1]\]  
This, of course, introduces a cutoff frequency \(\omega_{max}\) related to the discretization time step: \(\omega_{max}=n\delta\omega=\frac{\pi}{\delta t}\). Since \(H(t)\) must be real,  
\[H(\ell\,\delta t)=H_{\ell}=\frac{1}{2n}\sum_{k=-n}^{n-1}\bar{H}_{k}\cos\frac{ \pi k\ell}{n} \tag{12}\]  
and,  
\[R(\ell\,\delta t)=R_{\ell}=\sum_{m=-n}^{n-1}H_{m}r_{\ell-m} \tag{13}\]  
where \(r_{m}\) is a random number with Gaussian distribution, zero mean and variance \(\sqrt{\delta t}\). In practice, the time step used for the noise generation does not have to be equal to the dynamical time step \(\delta t\): [PERSON] and [PERSON] introduce the possibility to choose a larger time step \(h=M\delta t\), by keeping the random force constant for \(M\) consecutive dynamical steps. The convolution formula then becomes:  
\[H(\ell h)=H_{\ell}=\frac{1}{2n}\sum_{k=-n}^{n-1}\hat{H}_{k}\cos\frac{\pi k \ell}{n}\;\;\text{and}\;R(\ell h)=R_{\ell}=\sum_{m=-n}^{n-1}H_{m}r_{\ell-m}\]  
This has two advantages: firstly it allows updating the random force every \(M\) time steps only, and dividing the number of convolution points \(2n\) by \(M\) while keeping the same time interval \(T=2 nM\delta t\) for the convolution. This considerably reduces the computationload. Secondly, it naturally decreases the cutoff frequency to \(\omega_{max}=\pi/h=\pi/M\delta t\), which can be useful to avoid kinetic energy overestimation, as indicated in Section 4.1.