I would like to approximate the acceleration of a molecular dynamics system.

I'm following an online tutorial to solve a set of equations for molecular dynamics. I can use $F=ma$ to calculate the acceleration, but I don't understand how in that tutorial the tutor uses the two equations stated below with initial geometry and coordinates to get the acceleration.

We've been told to use the equation below. I only have the initial mass, speed and coordinates of the system and that $V$ has the Lennard Jones parameters. I do not know how to solve this problem as there seems to be too many unknowns. Any help will be appreciated.

\begin{align} a_{x_i}&=-\frac{1}{m_i}\frac{\partial V}{\partial x_i}\\ \frac{\partial V}{\partial x_i}&\sim\frac{V(x_i+\Delta x_i)-V(x_i-\Delta x_i)}{2\Delta x_i} \end{align}


closed as unclear what you're asking by Aaron Stevens, John Rennie, Jon Custer, alephzero, AccidentalFourierTransform Nov 6 '18 at 17:43

Please clarify your specific problem or add additional details to highlight exactly what you need. As it's currently written, it’s hard to tell exactly what you're asking. See the How to Ask page for help clarifying this question. If this question can be reworded to fit the rules in the help center, please edit the question.

  • 3
    $\begingroup$ This site isn't for homework problems but conceptual ones. You have a formula, what is your conceptual difficulty applying it? Also, I think some of the text in your post is missing and it makes it hard to tell what you mean when you say "I do not Any help will be grateful. " $\endgroup$ – jacob1729 Nov 5 '18 at 22:27
  • $\begingroup$ Do you have a function for $V$? $\endgroup$ – PM 2Ring Nov 5 '18 at 22:34
  • 1
    $\begingroup$ Note that $V$ seems to be a potential energy, not a velocity $\endgroup$ – Aaron Stevens Nov 5 '18 at 22:42
  • $\begingroup$ Sorry I fixed the missing part of the post. I am not looking to solve homework, I'm following an online tutorial to solve a set of equations for molecular dynamics. I can use F=ma to calculate the acceleration, but I don't understand how in that tutorial the tutor uses these two equations with initial geometry and coordinates to get the acceleration. $\endgroup$ – Raymond Ghaffarian Shirazi Nov 5 '18 at 22:44
  • $\begingroup$ If V is potential I have the Lennard Jones parameters. $\endgroup$ – Raymond Ghaffarian Shirazi Nov 6 '18 at 12:18

From your question and your comments, I have inferred that you would like to understand how to construct a numerical method for deterministic MD. You are given initial velocities and positions of the atoms and that they interact through pairwise Lennard-Jones potential.

Part 1 - Gradient of the Potential Energy

Assume that

  • $n_a \in \mathbb{N}_{+}$ is the number of atoms
  • $d \in \{ 1, 2, 3\}$ is the number of physical dimension of the system

In what follows, by abuse of notation, if $q \in \mathbb{R}^{d\times n_a}$, then $q_{\cdot, i} = (q_{0,i},\ldots,q_{d-1,i})^T \in \mathbb{R}^d$ for all $q \in \mathbb{R}^{d\times n_a}$ and for all admissible $i$. Furthermore, $\mathcal{D}^j_i$ denotes a partial derivative of order $j$ with respect to the variable with the index $i$, and $\mathcal{D}_{\cdot, i}$ represents the gradient with respect to the variable with the index $(\cdot, i)$ using the index notation as described above.

Furthermore, assume

  • $C=\{ q \in \mathbb{R}^{d\times n_a} : (\exists i < n_a) (\exists j < n_a) (i \neq j \wedge q_{\cdot,i} = q_{\cdot,j}) \}$
  • $q_0 \in \mathbb{R}^{d\times n_a}\setminus C$ is the vector of the initial positions
  • $v_0 \in \mathbb{R}^{d\times n_a}$ is the vector of the initial velocities
  • $d_r:\mathbb{R}^d \times \mathbb{R}^d \longrightarrow \mathbb{R}$ is the distance function that is given by $d_r(x_0, x_1)=\lVert x_0 - x_1 \rVert_2$ for all $x_0, x_1 \in \mathbb{R}^d$.
  • $r:\mathbb{R}^d \times \mathbb{R}^d \longrightarrow \mathbb{R}^d$ is the position difference function given by $r(x_0, x_1)=x_0 - x_1$ for all $x_0, x_1 \in \mathbb{R}^d$.

Given that the system interacts through Lennard-Jones potential, the potential energy function is defined below (note: all indices start from $0$):

  • $V:\mathbb{R}^{d\times n_a} \longrightarrow \mathbb{R}$ is the potential energy function given by $V(q)=\sum_{i=0}^{i=n_a-2}\sum_{j=i+1}^{j=n_a-1}\phi (d_r(q_{\cdot, i}, q_{\cdot, j} ))$ for any vector of positions $q \in \mathbb{R}^{d\times n_a}$
  • $\phi:\mathbb{R}_{+}\longrightarrow \mathbb{R}$ is the pairwise Lennard-Jones potential given by $\phi (r) = 4 \epsilon ( (\frac{\sigma}{r})^{12} - (\frac{\sigma}{r})^{6} ) $ for any distance $r\in\mathbb{R}_{+}$.

Suppose that $q:[0,t_e]\longrightarrow \mathbb{R}^{d\times n_a}$ represents the positions of the atoms of the system and $M \in \mathbb{R}^{dn_a \times dn_a}$ is the diagonal mass matrix. Then it is of interest to solve the IVP (Newton's equation of motion)

  • $(\mathcal{D}_1^2 q)(t) = -M^{-1} (\nabla V) (q(t))$ on $t \in (0,t_e)$ subject to $q(0) = q_0$ and $(\mathcal{D}_1 q)(0) = v_0$.

I believe that the equation above is identical to the equation in your question. Indeed, restating it in the index notation results in the system of equations

  • $(\mathcal{D}_1^2 q_{\cdot, i})(t) = -\frac{1}{M_{id,id}} ((\mathcal{D}_{\cdot, i} V) (q (t)))$,

where $i \in \{0,\ldots, n_a-1\}$, i.e. each $i$ represents a particular atom.

However, it is possible to find a closed-form analytical expression for $\mathcal{D}_{\cdot, i} V: \mathbb{R}^{d\times n_a}\setminus C \longrightarrow \mathbb{R}^d$ for all $i \in \{0,\ldots,n_a-1\}$ using the chain rule:

$(\mathcal{D}_{\cdot, i} V) (q) = \sum_{j=0,j\neq i}^{j=n_a - 1} ((\mathcal{D} \phi) (\mathcal{D}_{q_{\cdot, i}} d_r)) (q)$ for all $q \in \mathbb{R}^{d \times n_a} \setminus C$.

Thus, if $\mathbf{r}_{ij}= r (q_{\cdot, i}, q_{\cdot, j})$ and $r_{ij}= d_r(q_{\cdot, i}, q_{\cdot, j})$, then

$(\mathcal{D}_{\cdot, i} V) (q) = \sum_{j=0,j\neq i}^{j=n_a - 1} (4\epsilon ( -12 r_{ij}^{-13}\sigma^{12} + 6 r_{ij}^{-7}\sigma^6)) (r_{ij}^{-1} (q_{\cdot, i} - q_{\cdot, j}))$

or, equivalently,

$(\mathcal{D}_{\cdot, i} V) (q) = -\sum_{j=0,j\neq i}^{j=n_a - 1} 48\epsilon ( r_{ij}^{-14}\sigma^{12} - \frac{1}{2} r_{ij}^{-8}\sigma^6 ) \mathbf{r}_{ij}$.

Therefore, you can restate the Newton equation (in components) as

$(\mathcal{D}_1^2 q_{\cdot, i})(t) = \frac{1}{M_{id,id}} \sum_{j=0,j\neq i}^{j=n_a - 1} 48\epsilon ( r_{ij}^{-14}\sigma^{12} - \frac{1}{2} r_{ij}^{-8}\sigma^6 ) \mathbf{r}_{ij}$.

It should be noted that in MD usually a truncated potential energy function is used. All of them share the common property $V(r) = 0$ for all $r \geq r_c$, where $r_c$ is called cutoff radius. There are many different analytical expressions that achieve this. Ideally, it is of interest that $V$ is as smooth as possible without being too complex from the computational perspective. Most of these functions have closed form analytical expressions for their partial derivatives and, of course, the procedure for the calculation of these derivatives is identical to the one presented above.

In conclusion, for deterministic MD with the Lennard-Jones potential it is not necessary to use the numerical approximation to evaluate the gradient of the potential energy function, because an exact analytical expression can be used instead.

Part 2 - Numerical Discretisation of the Newton Equation

One of the most common numerical methods for the discretisation of the equation from Part 1 is called Velocity Verlet. There are several commonly used alternative formulations of the method. Some of these may have certain advantages over the others.

The main reason why Velocity Verlet is popular for the numerical simulation of the Newton equation is that it preserves certain important geometric characteristics of the flow of the (non-discretised) equation. Of course, there exist methods other than the Velocity Verlet that preserve the same characteristics, but they are usually more complex from the perspective of analysis. For more information see the book "Molecular Dynamics With Deterministic and Stochastic Numerical Methods" by Ben Leimkuhler and Charles Matthews and related publications. There is also a nice introductory article on Wikipedia: "Verlet integration". Unfortunately, it does not seem to have the information about the potential function that is presented above.

In what follows $f_i:\mathbb{R}^{d \times n_a}\setminus C \longrightarrow \mathbb{R}^d$ is used to denote $f_i = -(\mathcal{D}_{\cdot,i} V)$, $m_i = M_{id,id}$, $h \in \mathbb{R}_{+}$ is the size of the time step of the method and $k \in \mathbb{N}$ is the index of the time step. One formulation of the Velocity Verlet method is

$ v_{\cdot, i}(t_{k + \frac{1}{2}}) = v_{\cdot, i}(t_k) + \frac{h}{2} m_i^{-1} f_i (q(t_k))\\ q_{\cdot, i}(t_{k + 1}) = q_{\cdot, i}(t_k) + hv_{\cdot, i}(t_{k + \frac{1}{2}})\\ v_{\cdot, i}(t_{k + 1}) = v_{\cdot, i}(t_{k + \frac{1}{2}}) + \frac{h}{2}m_i^{-1}f_i(q(t_{k + 1})) $

Of course, I can provide more details if you are interested - please let me know.


Not the answer you're looking for? Browse other questions tagged or ask your own question.