Generalized Green's function molecular dynamics for canonical ensemble simulations

The need of small integration time steps ($\sim 1$ fs) in conventional molecular dynamics simulations is an important issue that inhibits the study of physical, chemical, and biological systems in real timescales. Additionally, to simulate those systems in contact with a thermal bath, thermostating techniques are usually applied. In this work, we generalize the Green's function molecular dynamics technique to allow simulations within the canonical ensemble. By applying this technique to one-dimensional systems, we were able to correctly describe important thermodynamic properties such as the temperature fluctuations, the temperature distribution, and the velocity autocorrelation function. We show that the proposed technique also allows the use of time steps one order of magnitude larger than those typically used in conventional molecular dynamics simulations. We expect that this technique can be used in long-timescale molecular dynamics simulations.

Figure 1

Evolution of the normalized displacement of the central atom $[u_0\left(t\right)/d_0]$ obtained from the $\mathrm{e}\text{−}\mathrm{GFMD}$ and velocity-Verlet methods within the microcanonical ensemble. The exact solution for a linear chain with 23 atoms is also shown for comparison [4]. The inset shows the total energy variation, normalized by the initial total energy, for a time step of 1 fs. Data from $\mathrm{e}\text{−}\mathrm{GFMD}$ are scale up 100 times to a better visualization.

Figure 2

Average kinetic temperature (symbols) and temperature standard deviation (error bars) as a function of the damping time for different time steps ($N=203$) for $\mathrm{e}\text{−}\mathrm{GFMD}$. The equilibration time was 2000 steps and the average as well as the standard deviation were calculated for the following 1000 steps. The blue horizontal line indicates the target temperature. The inset shows the behavior for the BBK (red) and G-JF (blue) methods. The time steps are indicated with the same symbols used for $\mathrm{e}\text{−}\mathrm{GFMD}$.

Figure 3

Relative kinetic temperature fluctuation $\mathrm{\Delta }{T}_k/T_k$ as a function of the size of the linear chain for different target temperatures $T$. The straight line corresponds to the expected behavior for the canonical ensemble ($\mathrm{\Delta }{T}_k/T_k=\sqrt{2/N}$). For each chain, the system was equilibrated for 1 ns and the average was calculated for the following 1 ps ($\tau =100\mathrm{fs}$, $\mathrm{\Delta }t=1\mathrm{fs}$). The inset shows the kinetic temperature distribution for the final 1 ps for a chain of 223 atoms compared with the Boltzmann distribution (black curve).

Figure 4

Distributions of the velocities for a linear chain ($N=203$, $T=300\mathrm{K}$, $\tau =100\mathrm{fs}$) for time steps of 1 and 10 fs after ${10}^5$ time steps. The curves represent the Gaussian fittings.

Figure 5

Normalized velocity autocorrelation function obtained from the $\mathrm{e}\text{−}\mathrm{GFMD}$, BBK, and G-JF methods for $\mathrm{\Delta }t=10\mathrm{fs}$ ($N=203$, $\tau =100\mathrm{fs}$, $T=30\mathrm{K}$, and total simulation time of 10 ns).

Figure 6

Convergence of the kinetic temperature ($T_k$) and the potential energy ($E_p$) (and their respective standard deviations) with decreasing time step for $\mathrm{e}\text{−}\mathrm{GFMD}$, BBK, and G-JF methods. The properties were averaged for 25 000 time steps after an equilibration of $5\times {10}^5$ time steps for $\mathrm{e}\text{−}\mathrm{GFMD}$ and ${10}^6$ time steps for BBK, and G-JF methods ($N=203$, $\tau =100\mathrm{fs}$, and $T=300\mathrm{K}$).

Figure 7

Computation time of the $\mathrm{e}\text{−}\mathrm{GFMD}$ method for 100 time steps ($\mathrm{\Delta }t=1\mathrm{fs}$, $\tau =100\mathrm{fs}$, and $T=300\mathrm{K}$). The simulations were carried out on a 3.3 GHz-IBM Power 755 processor. The curves for costs proportional to $N^3$ and $N^2$ are also indicated for comparison.