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$. 804