Molecular Dynamics: Theory

Microcanonical (NVE) ensemble

The Hamiltonian for the microcanonical ensemble is,

\[\mathcal{H} = \sum_{i=1}^N \frac{\mathbf{p}_i^2}{2m_i} + U(\mathbf{r}_i)\]

where \(\mathbf{p}_i\) and \(\mathbf{r}_i\) are the position and momentum of particle \(i\) and \(U\) is the DFT total (potential) energy. Hamilton’s equations can be solved to give the following equations of motion:

\[\begin{split}\mathbf{\dot{r}}_i &= \frac{\mathbf{p}_i}{m_i} \\ \mathbf{\dot{p}}_i &= \frac{\partial U(\mathbf{r}_i)}{\partial\mathbf{r}_i} = \mathbf{F_i}\end{split}\]

In order to construct a time-reversible algorithm from these equations, the Liouvillian formulation is employed [Ta1] (trivially, in this case). The Liouville operator \(L\) can be defined in terms of position and momentum components:

\[iL = \mathbf{\dot{r}}\frac{\partial}{\partial\mathbf{r}} + \mathbf{\dot{p}}\frac{\partial}{\partial\mathbf{p}} = i(L_r + L_p).\]

The Liouvillian can be used to construct the classical propagator, which relates the state \(f\) of the system at time 0 to its state at time \(t\):

\[f[\mathbf{p}^N(t),\mathbf{r}^N(t)] = e^{iLt}f[\mathbf{p}^N(0),\mathbf{r}^N(0)]\]

Taking the individual position and momentum parts of the Liouvillian \(L_r\) and \(L_p\), it can be shown that applying it to the state \(f\) result in a simple linear shift in coordinates and a simple linear shift in momentum respectively:

\[\begin{split}iL_rf(t) &= f[\mathbf{p}^N(0),\mathbf{r}^N(0) + \mathbf{\dot{r}}^N(0)t] \\ iL_pf(t) &= f[\mathbf{p}^N(0) + \mathbf{F}^N(0)t,\mathbf{r}^N(0)]\end{split}\]

However, we cannot simply replace \(e^{iLt}\) with \(e^{iL_rt}\) because \(iL_r\) and \(iL_p\) are non-commuting operators, so we must employ the Trotter-Suzuki identity:

\[e^{A+B} = \lim_{P\rightarrow\infty}\left(e^{A/2P}e^{B/P}e^{A/2P}\right)^P\]

Thus for a small enough time step \(\Delta t = t/P\) and to first order, a discrete time step corresponds to the application of the discrete time propagator \(G\),

\[G(\Delta t) = U_1\left(\frac{\Delta t}{2}\right)U_2\left(\Delta t\right)U_1\left(\frac{\Delta t}{2}\right) = e^{iL_1\frac{\Delta t}{2}}e^{iL_2\Delta t}e^{iL_1\frac{\Delta t}{2}},\]

which can be shown to be unitary and therefore time-reversible. Applying the operators \(U\) in the sequence determined by the Trotter decomposition generates the velocity Verlet algorithm, which is used to integrate microcanonical molecular dynamics in CONQUEST. For a detailed derivation of the algorithm, refer to Frenkel & Smit [Ta1].

Go to top.

Extended Lagrangian Born-Oppenheimer MD (XL-BOMD)

If the electronic density from the previous ionic step is used as an initila guess for the next SCF cycle, a problem arises because this process breaks the time-reversibility of the dynamics. This is manifested as a gradual drift in the total energy in the case of a NVE simulation, or the conserved quantity in the case of non-Hamiltonian dynamics. The solution proposed by Niklasson [Ta2][Ta3] is to introduce auxilliary electronic degrees of freedom into the Lagrangian, which can be propagated via time-reversible integrators.

The extended Lagrangian used in CONQUEST is [Ta4],

\[\mathcal{L}^\mathrm{XBO}\left(\mathbf{X}, \mathbf{\dot{X}}, \mathbf{R}, \mathbf{\dot{R}}\right) = \mathcal{L}^\mathrm{BO}\left(\mathbf{R}, \mathbf{\dot{R}}\right) + \frac{1}{2}\mu\mathrm{Tr}\left[\mathbf{\dot{X}}^2\right] - \frac{1}{2}\mu\omega^2\mathrm{Tr}\left[(\mathbf{LS} - \mathbf{X})^2\right],\]

where \(\mathbf{X}\) is a sparse matrix with the same range as \(\mathbf{LS}\), \(\mu\) is the fictitious electron mass and \(\omega\) is the curvature of the auxiliary harmonic potential. The Euler-Lagrange equations of motion are then,

\[\begin{split}m_i\mathbf{\ddot{r}_i} &= -\frac{\partial U[{{\mathbf{R;LS}}}]}{\partial\mathbf{r}_i} = \mathbf{F_i} \\ \mathbf{\ddot{X}} &= \omega^2(\mathbf{LS} - \mathbf{X}),\end{split}\]

The first of these is simply Newton’s second law, and the velocity update equation of motion in the microcanonical ensemble. The second can be integrated using a time-reversible algorithm, the velocity Verlet scheme in the case of CONQUEST [Ta4]:

\[\begin{split}\mathbf{X}(t+\delta t) &= 2\mathbf{X}(t) -\mathbf{X}(t-\delta t) + \delta t^2\omega^2\left[\mathbf{L}(t)\mathbf{S}(t)-\mathbf{X}(t)\right] \\ &+ a\sum_{m=0}^M c_m\mathbf{X}(t-m\delta t)\end{split}\]

i.e. the trajectory of \(\mathbf{X}(t)\) is time-reversible, and evolves in a harmonic potential centred on the ground state density \(\mathbf{L}(t)\mathbf{S}(t)\). The matrix \(\mathbf{XS}^{-1}\) is a good guess for the \(\mathbf{L}\) matrix in the Order(N) scheme.

Despite the time-reversitibility, the \(\mathbf{X}\) matrix tends in practice to gradually drift from the harmonic centre over time, increasing the number of SCF iterations required to reach the minimum over the course of the simulation. To remove such numerical errors, the final dissipative term is included, and is found to have a minimal effect on the time-reversibility. We note that since the auxiliary variable \(X\) is used to generate an intial guess for the SCF process, it does not appear in the conserved (pseudo-Hamiltonian) quantity for the dynamics.

Go to top.

Non-Hamiltonian dynamics

Extended system method

Hamiltonian dynamics generally describes systems that are isolated from their surroundings, but in the canonical and isobaric-isothermal ensembles, we need to couple the system to an external heat bath and/or stress. It is possible to model such systems by positing a set of equations of non-Hamiltonian equations of motion, and proving that they generate the correct statistical ensemble [Ta5]. This is the extended system approach: we modify the Hamiltonian to include the thermostat and/or barostat degrees of freedom, derive the (pseudo-) Hamiltonian equations of motion, and demostrate that the correct phase space distribution for the ensemble is recovered.

Go to top.

Canonical (NVT) ensemble

In the Nose-Hoover formulation [Ta6][Ta7], the Hamiltonian for a system in the canonical ensemble can be written,

\[\mathcal{H} = \sum_i \frac{1}{2}m_i s^2\mathbf{\dot{r}}_i^2 + U(\mathbf{r}_i) + \frac{1}{2}Q\dot{s}^2 - (n_f + 1)k_B T \ln s,\]

where \(\mathbf{r}_i\) and \(\mathbf{\dot{r}_i}\) are respectively the position and velocity of particle \(i\), \(U\) is the potential energy, in this case the DFT total energy, \(s\) is a dimensionless quantity that can be interpreted post-hoc as a time step scaling factor, \(Q\) is the fictitious mass of the heat bath and \(n_f\) is the number of ionic degrees of freedom. Hamilton’s equations can be solved to generate the Nose-Hoover equations of motion. However Martyna et al. demonstrate that this method does not generate an ergodic trajectory, and proposed an alternative formulation [Ta8] in which the temperature is controlled by a chain of \(M\) coupled thermostats of mass \(Q_k\), notional position \(\eta_k\) and conjugate momentum \(p_{\eta_k}\):

\[\begin{split}\mathbf{\dot{r}_i} &= \frac{\mathbf{p}_i}{m_i} \\ \mathbf{\dot{p}_i} &= -\frac{\partial U(\mathbf{r})}{\partial \mathbf{r}_i} - \frac{p_{\eta_1}}{Q_1}\mathbf{p}_i \\ \dot{\eta}_k &= \frac{p_{\eta_k}}{Q_k} \\ \dot{p}_{\eta_1} &= \left(\sum_{i=1}^N\frac{\mathbf{p}_i}{m_i} - n_fk_BT\right) - \frac{p_{\eta_{2}}}{Q_{\eta_{2}}}p_{\eta_1} \\ \dot{p}_{\eta_k} &= \left(\frac{p^2_{\eta_{k-1}}}{Q_{k-1}} - k_BT\right) - \frac{p_{\eta_{k+1}}}{Q_{k+1}}p_{\eta_k} \\ \dot{p}_{\eta_M} &= \left(\frac{p^2_{\eta_{M-1}}}{Q_{M-1}} - k_BT\right)\end{split}\]

The Liouvillian for these equations of motion can be non-uniquely decomposed into components of ionic position (\(iL_r\)) and momentum (\(iL_p\)) as in the microcanonical case, the extended Lagrangian (\(iL_\mathrm{XL}\), and a Nose-Hoover chain component (\(iL_\mathrm{NHC}\))

\[iL = iL_\mathrm{NHC} + iL_p + iL_{\mathrm{XL}} + iL_r,\]

which is directly translated into an algorithm with the Trotter-Suzuki expansion,

\[\begin{split}\exp(iL\Delta t) = &\exp\left(iL_\mathrm{NHC}\frac{\Delta t}{2}\right)\exp\left(iL_p\frac{\Delta t}{2}\right) \times \\ &\exp\left(iL_\mathrm{XL}\frac{\Delta t}{2}\right)\exp\left(iL_r\Delta t\right)\exp\left(iL_\mathrm{XL}\frac{\Delta t}{2}\right) \times \\ &\exp\left(iL_p\frac{\Delta t}{2}\right)\exp\left(iL_\mathrm{NHC}\frac{\Delta t}{2}\right)\end{split}\]

This is recognisable as the velocity Verlet algorithm with extended Lagrangian integration which can be reduced to a single step, as described in Extended Lagrangian Born-Oppenheimer MD (XL-BOMD), with a half time step integration of the Nose-Hoover chain equations of motion before and after. For full details of the integration scheme, see Hirakawa et al. [Ta9].

Go to top.

Isobaric-Isothermal (NPT) ensemble

The Parinello-Rahman equations of motion [Ta10] extend the fixed cell equations of motion to include the cell degrees of freedom in the extended system approach. We use the Martyna-Tobias-Tuckerman-Klein modification [Ta11], which couples the variable cell equations of motion to a Nose-Hoover chain the themrostat the system, recovering the isobaric-isothermal (NPT) ensemble. For an unconstrained cell (i.e. the lattice vectors can change freely), the equations of motion are,

\[\begin{split}\mathbf{\dot{r}}_i &= \frac{\mathbf{p}_i}{m_i} + \frac{\mathbf{p}_g}{W_g}\mathbf{r}_i \\ \mathbf{\dot{p}}_i &= \mathbf{F}_i - \frac{\mathbf{p}_g}{W_g}\mathbf{p}_i - \left(\frac{1}{N_f}\right)\frac{\mathrm{Tr}[\mathbf{p}_g]}{W_g}\mathbf{p}_i - \frac{p_\xi}{Q}\mathbf{p}_i \\ \mathbf{\dot{h}} &= \frac{\mathbf{p}_g\mathbf{h}}{W_g} \\ \mathbf{\dot{p}_g} &= V(\mathbf{P}_\mathrm{int}-\mathbf{I}P_\mathrm{ext}) + \left[\frac{1}{N_f}\sum_{i=1}^N\frac{\mathbf{p}_i^2}{m_i}\right]\mathbf{I} - \frac{p_\xi}{Q}\mathbf{p}_g \\ \dot{\xi} &= \frac{p_\xi}{Q} \\ \mathbf{\dot{p}}_g &= \sum_{i=1}^N\frac{\mathbf{p}_i^2}{m_i} + \frac{1}{W_g}\mathrm{Tr}[\mathbf{p}_g^T\mathbf{p}_g] - (N_f + d^2)kT\end{split}\]

Here, \(\mathbf{r}_i\), \(\mathbf{p}_i\) and \(m_i\) are the position, momentum and mass of particle \(i\) respectively, \(\xi\), \(p_\xi\) and \(Q\) are the position, momentum and mass of the thermostat, and \(\mathbf{h}\), \(\mathbf{p}_g\) and \(W_g\) are the matrix of lattice vectors, matrix of cell velocities and cell mass respectively. Note that these equations only include one Nose-Hoover thermostat for simplicity. Conquest uses the Shinoda-Shiga-Mikami splitting of the Liouvillian [Ta12] to propagate the system. The Liouvillian is decomposed as,

\[iL = iL_r + iL_h + iL_v + iL_\mathrm{bath},\]

which can be further split,

\[\begin{split}iL_\mathrm{bath} &= iL_\mathrm{box} + iL_\mathrm{particles} \\ iL_\mathrm{box} &= iL_\mathrm{vbox} + iL_\xi + iL_{v_{\xi_1}} + iL_{v_{\xi_k}} + iL_{v_{\xi_M}} \\ iL_\mathrm{particles} &= iL_\mathrm{vpart} + iL_\xi + iL_{v_{\xi_1}} + iL_{v_{\xi_k}} + iL_{v_{\xi_M}}\end{split}\]

Using Liouville’s theorem, we have,

\[\begin{split}iL_r &= \sum_{i=1}^N[\mathbf{v}_i + \mathbf{v}_g\mathbf{r}_i]\cdot\nabla_{\mathbf{r}_i} \\ iL_h &= \sum_{\alpha,\beta}\mathbf{v}_{g,\alpha\beta}\mathbf{h}_{\alpha\beta}\frac{\partial}{\partial\mathbf{h}_{\alpha\beta}} \\ iL_v &= \sum_{i=1}^N\left(\frac{\mathbf{F}_i}{m_i}\right)\cdot\nabla_{\mathbf{v}_i} \\ iL_\mathrm{bath} &= iL_\mathrm{vpart} + iL_\mathrm{vbox} + iL_\xi + iL_{v_{\xi_1}} + iL_{v_{\xi_k}} + iL_{v_{\xi_M}} \\ &= \sum_{i=1}^N\left[-\left\{\mathbf{v}_g + \frac{1}{N_f}\mathrm{Tr}(\mathbf{v}_g) + v_{\xi_1}\right\}\mathbf{v}_i\right]\nabla_{\mathbf{v}_i} \\ &+ \sum_{\alpha,\beta}\left[\frac{F_\mathrm{box}}{W} - v_{\xi_1}\mathbf{v}_{g,\alpha\beta}\right]\frac{\partial}{\partial\mathbf{v}_{g,\alpha\beta}} \\ &+ \sum_{k=1}^M v_{\xi_k}\frac{\partial}{\partial\xi_k} \\ &+ \left[\frac{F_{\mathrm{NHC}_1}}{Q_1} - v_{\xi_1}v_{\xi_2}\right]\frac{\partial}{\partial v_{\xi_1}} \\ &+ \sum_{k=2}^M\left[\frac{1}{Q_k}(Q_{k-1}v_{\xi_{k-1}}^2 - kT_\mathrm{ext}) - v_{\xi_k}v_{\xi_{k+1}}\right]\frac{\partial}{\partial v_{\xi_k}} \\ &+ \left[\frac{1}{Q_M}(Q_{M-1}v_{\xi_{M-1}}^2 - kT_\mathrm{ext})\right]\frac{\partial}{\partial v_{\xi_M}}\end{split}\]

Here we use \(M\) heat baths in a Nose-Hoover chain. The Trotter-Suzuki expansion is,

\[e^{iL\Delta t} = e^{iL_\mathrm{bath}\frac{\Delta t}{2}}e^{iL_v\frac{\Delta t}{2}}e^{iL_h\frac{\Delta t}{2}}e^{iL_r\Delta t}e^{iL_h\frac{\Delta t}{2}}e^{iL_v\frac{\Delta t}{2}}e^{iL_\mathrm{bath}\frac{\Delta t}{2}}.\]

The Liouvillian for the heat baths can be further expanded:

\[e^{iL_\mathrm{particles}\frac{\Delta t}{2}} = e^{\left(iL_{v_{\xi_1}} + iL_{v_{\xi_k}} + iL_{v_{\xi_M}}\right)\frac{\Delta t}{4}}e^{\left(iL_\xi + iL_\mathrm{vpart}\right)\frac{\Delta t}{2}}e^{\left(iL_\xi + iL_{v_{\xi_1}} + iL_{v_{\xi_k}} + iL_{v_{\xi_M}}\right)\frac{\Delta t}{4}}\]

Finally, expanding the first propagator in the previous expression, we have,

\[\begin{split}e^{\left(iL_{v_{\xi_1}} + iL_{v_{\xi_k}} + iL_{v_{\xi_M}}\right)\frac{\Delta t}{4}} &= e^{-i\left(-v_{\xi_1}v_{\xi_2}\frac{\partial}{\partial \xi_1} - \sum_{k=2}^Mv_{\xi_k}v_{\xi_{k+1}}\frac{\partial}{\partial \xi_k} - v_{\xi_{M-1}}v_{\xi_M}\frac{\partial}{\partial \xi_M}\right)\frac{\Delta t}{8}} \\ &\times e^{i\left(F_{\mathrm{NHC}_1}\frac{\partial}{\partial v_{\xi_1}} + F_{\mathrm{NHC}_k}\frac{\partial}{\partial v_{\xi_k}} + F_{\mathrm{NHC}_M}\frac{\partial}{\partial v_{\xi_M}}\right)\frac{\Delta t}{4}} \\ &\times e^{-i\left(-v_{\xi_1}v_{\xi_2}\frac{\partial}{\partial \xi_1} - \sum_{k=2}^Mv_{\xi_k}v_{\xi_{k+1}}\frac{\partial}{\partial \xi_k} - v_{\xi_{M-1}}v_{\xi_M}\frac{\partial}{\partial \xi_M}\right)\frac{\Delta t}{8}}\end{split}\]

These expressions are directly translated into the integration algorithm.

Go to top.

Weak coupling thermostat/barostat

Instead of modifying the Hamiltonian, the Berendsen-type weak coupling method [Ta13] involves coupling the ionic degrees of freedom to a an external temperature and/or pressure bath via “the principle of least local perturbation consistent with the required global coupling.” Thermostatting is acheived via a Langevin-type equation of motion, in which the system is globally coupled to a heat bath and subjected to random noise:

\[m_i\ddot{\mathbf{r}}_i = \mathbf{F}_i + m_i \gamma\left(\frac{T_0}{T}-1\right)\dot{\mathbf{r}}_i,\]

where \(\gamma\) is a global friction constant chosen to be the same for all particles. This can be acheived in practice by rescaling the velocities \(\mathbf{v}_i \rightarrow \lambda\mathbf{v}_i\), where \(\lambda\) is,

\[\lambda = \left[ 1 + \frac{\Delta t}{\tau_T}\left(\frac{T_0}{T}-1\right)\right]^{\frac{1}{2}}\]

A similar argument can be applied for weak coupling to an external pressure bath. In the isobaric-isoenthalpic ensemble, the velocity of the particles can be expressed,

\[\dot{\mathbf{r}} = \mathbf{v} - \frac{\beta(P_0 - P)}{3\tau_P}\mathbf{r},\]

i.e. the fractional coordinates are scaled by a factor determined by the difference between the internal and external pressures, the isothermal compressibility \(\beta\) and a pressure coupling time constant $tau_P$. In the isotropic case, the cell scaling factor \(\mu\) can be expressed,

\[\mu = \left[ 1 - \frac{\Delta t}{\tau_P}(P_0 - P)\right]^{\frac{1}{3}},\]

where the compressibility is absorbed into the time time constant \(\tau_P\). Allowing for fluctuations of all cell degrees of freedom, the scaling factor becomes,

\[\mathbf{\mu} = \mathbf{I} - \frac{\beta\Delta t}{3\tau_P}(\mathbf{P}_0 - \mathbf{P})\]

While trivial to implement and in general stable, the weak-coupling method does not recover the correct phase space distribution for the canonical or isobaric-isothermal ensembles, for which the extended system method is required. It is no longer supported in CONQUEST, but the concepts are useful.

Go to top.

Stochastic velocity rescaling

Stochastic velocity rescaling (SVR) [Ta14] is a modification of the weak coupling method, in which a correctly constructed random force is added to enforce the correct NVT (or NPT) phase space distribution. The kinetic energy is rescaled such that the change in kinetic energy between thermostatting steps is,

\[dK = (\bar{K} - K)\frac{dt}{\tau} + 2\sqrt{\frac{K\bar{K}}{N_f}}\frac{dW}{\sqrt{\tau}}\]

where \(\bar{K}\) is the target kinetic energy (external temperature), \(dt\) is the time step, \(\tau\) is the time scale of the thermostat, \(N_f\) is the number of degrees of freedom and \(dW\) is a Wiener process. Practically, the particle velocities are rescaled by a factor of \(\alpha\), defined via,

\[\alpha^2 = e^{-\Delta t/\tau} + \frac{\bar{K}}{N_fK}\left(1-e^{-\Delta t/\tau}\right)\left(R_1^2 + \sum_{i=2}^{N_f}R_i^2\right) + 2e^{-\Delta t/2\tau}\sqrt{\frac{\bar{K}}{N_fK}\left(1-e^{-\Delta t/\tau}\right)R_1}\]

Where \(R_i\) is a set of \(N_f\) normally distributed random numbers with unitary variance. This method can be applied to thermostat the NPT ensemble by barostatting the system with the Parinello-Rahman method, and using the above expressions, but with additional \(R_i\)’s for the cell degrees of freedom, and thermostatting the cell velocities as well as the particle velocities [Ta15].


D. Frenkel and B. Smit. Understanding molecular simulation: from algorithms to application. Academic Press, 2002.


A. M. N. Niklasson. Extended Born-Oppenheimer Molecular Dynamics. Phys. Rev. Lett., 100:123004, 2008. doi:10.1103/PhysRevLett.100.123004.


A. M. N. Niklasson and M. J. Cawkwell. Generalized extended Lagrangian Born-Oppenheimer molecular dynamics. J. Chem. Phys., 141:164123, 2014. doi:10.1063/1.4898803.


M. Arita, D. R. Bowler, and T. Miyazaki. Stable and Efficient Linear Scaling First-Principles Molecular Dynamics for 10000+ Atoms. J. Chem. Theor. Comput., 10:5419, 2014. doi:10.1021/ct500847y.


M. E. Tuckerman. Statistical mechanics: theory and molecular simulations. Oxford Graduate Texts, 2010.


S. Nosé. A unified formulation of the constant temperature molecular dynamics methods. J. Chem. Phys., 81:511, 1984. doi:10.1063/1.447334.


W. G. Hoover. Canonical dynamics: Equilibrium phase-space distributions. Phys. Rev. A, 31:1695, 1985. doi:10.1103/PhysRevA.31.1695.


G. J. Martyna, M. L. Klein, and M. Tuckerman. Nosé–hoover chains: the canonical ensemble via continuous dynamics. J. Chem. Phys., 97:2635, 1992. doi:10.1063/1.463940.


T. Hirakawa, T. Suzuki, D. R. Bowler, and T. Miyazaki. Canonical-ensemble extended lagrangian born-oppenheimer molecular dynamics for the linear scaling density functional theory. J. Phys.: Condens. Matter, 29:405901, 2017. doi:10.1088/1361-648X/aa810d.


M. Parrinello and A. Rahman. Polymorphic transitions in single crystals: A new molecular dynamics method. J. Appl. Phys., 52:7182–7190, December 1981. doi:10.1063/1.328693.


G. J. Martyna, M. E. Tuckerman, D. J. Tobias, and M. L. Klein. Explicit reversible integrators for extended systems dynamics. Mol. Phys., 87:1117, 1996. doi:10.1080/002689799163235.


W. Shinoda, M. Shiga, and M. Mikami. Rapid estimation of elastic constants by molecular dynamics simulation under constant stress. Phys. Rev. B, 69:4396, 2004. doi:10.1103/PhysRevB.69.134103.


H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. Dinola, and J. R. Haak. Molecular dynamics with coupling to an external bath. J. Chem. Phys., 81:3684, 1984. doi:10.1063/1.448118.


G. Bussi, D. Donadio, and M. Parrinello. Canonical sampling through velocity rescaling. J. Chem. Phys., 126:014101, 2007. doi:10.1063/1.2408420.


G. Bussi, T. Zykova-Timan, and M. Parrinello. Isothermal-isobaric molecular dynamics using stochastic velocity rescaling. J. Chem. Phys., 130:074101, 2009. doi:10.1063/1.3073889.

Go to top.