Diffusive Molecular Dynamics

Diffusive Molecular Dynamics (DMD) introduced by Li et al. in 2011 [1] offers a new perspective of atomistic simulations, by examinig phase-field theory on an atomistic resolution grid. Li et al. derived the free energy equation for EAM. Later on Simpson and colleagues [2] proved rigoursly that as a matter of fact DMD is an proximation of an ensemble they introduced and offered a general formulation of it. It is highly recommended that interested reader examine the abovementioned references. Here we present a very generic formulation of DMD.

Suppose that we have a complicated system with an incalculable partition function Z and potentail energy U. On the other hand we have a similar system with a calculable partition function, \widetilde{Z} and potential energy \widetilde{U}. Gibbs-Bogoliubov and coincidently relative entropy formulation of information theory prove the following inequality:

-\beta^{-1} \log Z \le   \langle U- \widetilde{U} \rangle_{\widetilde{U}} -\beta^{-1} \log \widetilde{Z}

Here \beta is 1/k_B T where k_B, and T denote Boltzmann factor and the temperature respectively. Now suppose that asside from degrees of freedom (\mathbf{x}) our simplified system has extra adjustable parameters (\mathbf{x} and \mathbf{\alpha}), above inequality becomes

-\beta^{-1} \log Z \le   \langle U- \widetilde{U}(\mathbf{x},\mathbf{\alpha}) \rangle_{\widetilde{U}(\mathbf{x},\mathbf{\alpha})} -\beta^{-1} \log \widetilde{Z}(\mathbf{x},\mathbf{\alpha})

Naturally, minimizing the right side of the inequality with repect to the adjustable parameters, will get us as close as possible to -\beta^{-1} \log Z or otherwise known as free energy (F). That is the main idea of DMD approximating a compilicated system with a simplified system with extra parameters.

Let us be more specific, in a general DMD simulation there are d + 2\times n degrees of freedom per atom (or more accurately per site), where d is the dimension of simulation, and n is number of elements present in system. n of degrees of freedom represent the propbablity of it being occupuied by each elements (\mathbf{c}), d of them denote the position of the site (\mathbf{x}), and last n stiffness of the site associated with each element (\mathbf{\alpha}). It is evident that for every site.

0 \leq \sum_{\gamma=0}^{n-1} c_i^\gamma \leq 1,\quad \quad 0 \leq i <N

where N is number of sites. Here \gamma is merely a label presented in form of a superscript for conveniance. It is worth noting that unlike MD, there is no velocity degrees of freedom.


\alpha ‘s definition in mapp4py and this manual is different from the one Li et al. provided:

\alpha=\frac{1}{\sqrt{\alpha_{\mathrm{Li\  et\  al.}}}}

this way \alpha like x is of unit of length.

As was mentioned before \mathbf{x} and \mathbf{\alpha} are adjustable parameters and they need to adjusted in order to minimize error of our simplification. However, \mathbf{c} is a different story and it should not be adjusted to minimize the free energy. Instead it is evolved using a master equation. The proposed algorithm to perform a dmd simulation is outlined in the figure below:


DMD algorithm outline


This is done by minimizing the relative entropy by means of varying \mathbf{x} and \mathbf{\alpha}. It is tempting to achieve this byminimizing so called approximate free energy using traditional CG//L-BFGS methods. Although these methods are very effective for VG they are insufficient when it comes to DMD, especifically when we have sites that are vacant or close to being vacant. Suppose that

\widetilde{F}=\langle U- \widetilde{U} \rangle_{\widetilde{U}} -\beta^{-1} \log \widetilde{Z}

working out the derivitives with respect to \mathbf{x} and \mathbf{\alpha} of a site, namely i it can be shown that these dervatives are directely proportional to c_i. Therefore when the concentration of site i is zero or close to zero these derivatives vanish, leading to intial values of \mathbf{x}_i and \mathbf{\alpha}_i being unchanged. To remedy the situation instead we solve the following set of nonlinear equations

\frac{1}{c_i}\frac{\partial\widetilde{F}}{\partial \mathbf{x}_i}=0, \quad \frac{1}{c_i}\frac{\partial\widetilde{F}}{\partial \alpha_i}=0, \quad i=0,\cdots N-1

effictively removing the proportionality with respect to c_i. These nonlinear equations in MAPP are solved using newton method which in turn employes generalized minimal residual method (GMRES) linear solver. Although this approach seems more costly and elaborate due to calculation of Hessian matrix and inversion of it compared to minimization methods, faster rate convergence and accuracy of the results do compensate for it.

Evolution of Concentration

The master equation by means of which the site concentraions evolve are described in details here. Due to the stiff nature of master equations regular explicit methods such as forward euler are not viable options. Instead in MAPP we use backward differentiation formula algorithm to evaolve the concentration. This numerical method like all other explicit integrators requires solving a set of nonlinear equations at evrey step. However, the fact that timestep is variable and in a sense is adaptive, enables us to perform simulations for a long period. For more details see mapp4py.dmd.bdf


BDF integrator



Ju Li, S. Sarkar, W.t. Cox, T.j. Lenosky, E. Bitzek, and Yunzhi Wang. Diffusive molecular dynamics and its application to nanoindentation and sintering. Physical Review B (Condensed Matter and Materials Physics), 84(5):054103, August 2011. doi:10.1103/PhysRevB.84.054103.


G. Simpson, M. Luskin, and D. Srolovitz. A Theoretical Examination of Diffusive Molecular Dynamics. SIAM Journal on Applied Mathematics, 76(6):2175–2195, January 2016. URL: http://epubs.siam.org/doi/abs/10.1137/15M1024858, doi:10.1137/15M1024858.