Efficient Exploration of Reactive Potential Energy Surfaces Using Car-Parrinello Molecular Dynamics

The possibility of observing chemical reactions in ab initio molecular dynamics runs is severely hindered by the short simulation time accessible. We propose a new method for accelerating the reaction process, based on the ideas of the extended Lagrangian and coarse-grained non-Markovian metadynamics. We demonstrate that by this method it is possible to simulate reactions involving complex atomic rearrangements and very large energy barriers in runs of a few picoseconds.

Figure 1

C1-C4 distance during the MD run. The masses are set to 50 amu and the $k_{\alpha }$’s to 0.3. $V(t,\mathbf{s})$ is updated every 0.012 ps, $\Delta s^{\perp }$ is 0.15, $\Delta s^{\parallel }$ fluctuates in the interval $[0.07:0.1]$, and $W$ between 1.0 and $8.0\mathrm{k}\mathrm{c}\mathrm{a}\mathrm{l}/\mathrm{m}\mathrm{o}\mathrm{l}$. The arrows indicate which configuration the system is in: (a) is the cycle, (b) the trans, and (c) the cis configuration. With larger Gaussians [$\Delta s^{\perp }=0.2$, and consequently larger $\Delta s^{\parallel }$ and $W$ ($\simeq 6.5\mathrm{k}\mathrm{c}\mathrm{a}\mathrm{l}/\mathrm{m}\mathrm{o}\mathrm{l}$)] and smaller $M_{\alpha }$ (1 amu), the transition is observed within 1 ps, at a cost of a more coarse evaluation of the energy profile.

Figure 2

Examples of minimum energy structures in the ${\mathrm{S}\mathrm{i}}_6{\mathrm{H}}_8$ stoichiometry. (a) is the structure reported in Ref. 25, while (b) and (c) are encountered in the MD run performed with $\Delta s^{\perp }=0.3$ and $M_{\alpha }=1\mathrm{a}\mathrm{m}\mathrm{u}$, where $V(t,\mathbf{s})$ is updated with a new Gaussian every 0.003 ps. (b) and (c) are, respectively, 12.2 and $8.5\mathrm{k}\mathrm{c}\mathrm{a}\mathrm{l}/\mathrm{m}\mathrm{o}\mathrm{l}$ lower in energy than (a).

Figure 3

Potential energy during the MD run at 300 K. The $M_{\alpha }$ are equal to 50 amu and $V(t,\mathbf{s})$ is incremented by a new Gaussian every 0.006 ps. $\Delta s^{\perp }$ is equal to 0.1, $\Delta s^{\parallel }$ fluctuates between 0.02 and 0.055, and $W$ between 2 and $10\mathrm{k}\mathrm{c}\mathrm{a}\mathrm{l}/\mathrm{m}\mathrm{o}\mathrm{l}$. Some minimum energy configurations, encountered along the trajectory, are reported. (a) and (b) are ${\mathrm{S}\mathrm{i}}_6{\mathrm{H}}_8$, (c) is ${\mathrm{S}\mathrm{i}}_6{\mathrm{H}}_6$, and (d) and (e) are ${\mathrm{S}\mathrm{i}}_6{\mathrm{H}}_4$