Deep Potential Molecular Dynamics: a scalable model with the accuracy of quantum mechanics

Linfeng Zhang, Jiequn Han, Han Wang, Roberto Car, Weinan E

References

Appendix A Supplementary Materials

The data used for training and/or testing are extracted from the AIMD simulations summarized in Tab. 4. All the simulations adopt a time step of 0.48 fs. The PI-AIMD simulations use the CPMD codes of Quantum Espresso http://www.quantum-espresso.org/ for the DFT part and are interfaced with the i-PI code Ceriotti, More, and Manolopoulos (2014) for the path-integral part. The generalized Langevin equation with color noise Ceriotti, Manolopoulos, and Parrinello (2011) in i-PI requires 8 beads for a converged representation of the Feynman paths. The classical AIMD simulations use the CPMD codes of Quantum Espresso and adopt the Nosé-Hoover thermostat Martyna, Klein, and Tuckerman (1992) for thermalization. The Parrinello-Rahman technique Parrinello and Rahman (1980) for variable cell dynamics is adopted in all cases. The training datasets include 95000 snapshots (from 105000 total snapshots) randomly selected along the liquid water trajectory, 19500 snapshots (from 24000 total snapshots) randomly selected along the ice (b) trajectory, 9500 snapshots (from 12000 total snapshots) randomly selected along the ice (c) trajectory, and 9500 snapshots (from 12000 total snapshots) randomly selected along the ice (d) trajectory. The remaining snapshots in the database are used for testing purposes.

A.1.2 molecules

The data and their complete description for the organic molecules (benzene, uracil, napthalene, aspirin, salicylic acid, malonaldehyde, ethanol, and toluene) can be found at http://quantum-machine.org/. For each molecule, 95000 snapshots, randomly selected from the database, are used to train the DeePMD model. The remaining snapshots in the database are used for testing purposes.

A.2 Implementation of the method

The TensorFlow r1.0 software library (http://tensorflow.org/) is interfaced with our C++ codes for data training and for calculating the energy, the forces, and the virial.

We consider a system consisting of $N$ atoms. The global coordinates of the atoms, in the laboratory frame, are $\left\{\bm{R}_{1},\bm{R}_{2},\dots,\bm{R}_{N}\right\}$ , where $\bm{R}_{i}=\left\{x_{i},y_{i},z_{i}\right\}$ for each $i$ . The neighbors of atom $i$ are denoted by $\mathcal{N}(i)=\left\{j:|\bm{R}_{ij}|<R_{c}\right\}$ , where $\bm{R}_{ij}=\bm{R}_{i}-\bm{R}_{j}$ , and $R_{c}$ is the cut-off radius. The neighbor list $\mathcal{N}(i)$ is sorted according to the scheme illustrated in Fig. 1. In extended systems, the number of neighbors at different snapshots inside $R_{c}$ fluctuates. Let $N_{c}$ be the largest fluctuating number of neighbors. The two atoms used to define the axes of the local frame of atom $i$ are called the axis-atoms and are denoted by $a(i)\in\mathcal{N}(i)$ and $b(i)\in\mathcal{N}(i)$ , respectively. In general we choose two closest atoms, independently of their species, together with the center atom, to define the local frame. Thus, in all the water cases, we choose the other two atoms belonging to the same water molecule. We apply the same rule to the organic molecules, but in this case we exclude the hydrogen atoms in the definition of the axis-atoms.

Next, we define the rotation matrix $\mathcal{R}(\bm{R}_{ia(i)},\bm{R}_{ib(i)})$ for the local frame of atom $i$ ,

where $\bm{e}[\bm{x}]\equiv\frac{\bm{x}}{||\bm{x}||}$ . In this local frame of reference, we obtain the new set of coordinates:

and we define $R_{ij}^{\prime}=||\bm{R}_{ij}^{\prime}||$ . Then the spacial information for $j\in\mathcal{N}(i)$ is

When $\alpha=0,1,2,3$ , full (radial plus angular) information is provided. When $\alpha=0$ , only radial information is used. Note that for $j\in\mathcal{N}(i)$ , $D_{ij}^{\alpha}$ is a function of the global coordinates of three or four atoms:

This formula is useful in the derivation of the formulae for the forces and the virial tensor given below.

The neural network uses a fixed input data size. Thus, when the size of $\mathcal{N}(i)$ is smaller than $N_{c}$ , we temporarily set to zero the input nodes not used for storing the $D_{ij}^{\alpha}$ . The nodes set to zero are still labeled by $D_{ij}^{\alpha}$ .

The $D_{ij}^{\alpha}$ are then standardized to be the input data for the neural networks. In this procedure, the $D_{ij}^{\alpha}$ are grouped according to the different atomic species. Within each group we calculate the mean and standard deviation of each $D_{ij}^{\alpha}$ by averaging over the snapshots of the training sample and over all the atoms in the group. Then we shift the $D_{ij}^{\alpha}$ by their corresponding means, and divide them by their corresponding standard deviations. Because of the weight $1/R$ in the $D_{ij}^{\alpha}$ and because the unoccupied nodes are set to zero, some standard deviations are very small or even zero. This causes an ill-posed training process. Therefore, after the shift operations, we divide by 0.01 Å-1 the $D_{ij}^{\alpha}$ with standard deviation smaller than 0.01 Å-1. For simplicity, we still use the same notation for the standardized $D_{ij}^{\alpha}$ .

A.2.2 deep neural network for the energy

For atom $i$ , the “atomic energy” is represented as

where $N_{\bm{w}(i)}$ is the network that computes the atomic contribution to the total energy, and $\bm{w}(i)$ are the weights used to parametrize the network, which depend on the chemical species of atom $i$ .

In this work, $N_{\bm{w}(i)}$ is constructed as a feedforward network in which data flows from the input layer as $\{D_{ij}^{\alpha}\}$ , through multiple fully connected hidden layers, to the output layer as the atomic energy $E_{i}$ . More specifically, a feedforward neural network with $N_{h}$ hidden layers is a mapping

where the symbol “ $\circ$ ” denotes function composition. Here $\mathcal{L}_{i}^{p}$ is the mapping from layer $p-1$ to $p$ , a composition of a linear transformation and a non-linear transformation

It should be stressed that, to guarantee the permutational symmetry, atoms of the same species share the same parameters $\bm{w}$ .

A.2.3 forces and virial tensor

The total potential energy is the sum of the $E_{i}$ . Thus the forces are

The virial tensor is defined as $\Xi_{\alpha\beta}=-\frac{1}{2}\sum_{i}R_{i\alpha}F_{i\beta}$ , where the indices $\alpha$ and $\beta$ indicate Cartesian components in the lab reference frame. Due to the periodic boundary conditions, one cannot directly use the absolute coordinates $R_{i\alpha}$ to compute the virial tensor. Rather, in the AIMD framework, the virial tensor is defined with an alternative but equivalent formula, i.e.,

where $h$ is the cell tensor. In our framework, due to the decomposition of the local energy $E_{i}$ , one computes the virial tensor by:

$x_{\alpha}^{(i,j)}$ is the $\alpha$ -th component of the vector oriented from the $i$ -th to the $j$ -th atom in the difference:

${f}_{\beta}^{(i,j)}$ is the $\beta$ -th component of the negative gradient of $E_{i}$ w.r.t. $x_{j}$ , i.e.,

Together with the energy representation described above, all the quantities needed for training and MD simulations, although complicated, have been analytically defined. In particular, it is noted that the derivatives of the total energy with respect to the atomic positions, appearing in both the forces and the viral tensor, are computed by the chain rule through the backpropagation algorithm, provided by TensorFlow. To make it work, we additionally implement the computation of $\nabla_{\bm{R}_{i}}D^{\alpha}_{jk}$ in C++ and interface it with TensorFlow.

A.2.4 Training Details

During the training process, one minimizes the family of loss functions defined in the paper:

The network weights are optimized with the Adam stochastic gradient descent method Kingma and Ba (2015). An initial learning rate $r_{l0}=0.001$ is used with the Adam parameters set to $\beta_{1}$ =0.9, $\beta_{2}$ =0.999, and $\epsilon$ =1.0 $\times{10}^{-8}$ , which are the default settings in TensorFlow. The learning rate $r_{l}$ decays exponentially with the global step:

where $d_{r}$ , $c_{s}$ , and $d_{s}$ are the decay rate, the global step, and the decay step, respectively. In this paper, the batch size is 4 in all the training processes. The decay rate is 0.95. For liquid water, the training process undergoes 4000000 steps in total, and the learning rate is updated every 20000 steps. For molecules, the training process undergoes 8000000 steps in total, and the learning rate is updated every 40000 steps.

We remark that, for the prefactors, a proper linear evolution with the learning rate speeds up dramatically the training process. We define this process by:

in which $p_{start}$ is the prefactor at the beginning of the training process, and $p_{limit}$ is approximately the prefactor at the end. We define $p_{start}$ for the energy, the forces, and the virial as $p_{estart}$ , $p_{fstart}$ , and $p_{vstart}$ , respectively. Similarly, we define $p_{limit}$ for the energy, the forces, and the virial as $p_{elimit}$ , $p_{flimit}$ , and $p_{vlimit}$ , respectively. In this paper, we use the following scheme:

The above scheme is based on the following considerations. Each snapshot of the AIMD trajectories provides 1 energy, $3N$ forces, and 6 independent virial tensor elements. The number of force components is much larger than the number of energy and virial tensor components. Therefore, matching the forces at the very beginning of the training process makes the training efficient. As the training proceeds, increasing the prefactors of the energy and the virial tensor allows us to achieve a well balanced training in which the energy, the forces, and the virial are mutually consistent.

In the original Deep Potential paper Han et al. (2017), only the energy was used to train the network, requiring in some cases the use of Batch Normalization techniques Ioffe and Szegedy (2015) to deal with issues of overfitting and training efficiency. Adding the forces and/or the virial tensor provides a strong regularization of the network and makes training significantly more efficient. Thus Batch Normalization techniques are not necessary within the DeePMD framework.

A.2.5 DeePMD details

In the path-integral/classical $NPT$ simulations of liquid water and ice, we integrate our codes with the i-PI software. The DeePMD simulations are performed at the same thermodynamic conditions, and use the same temperature and pressure controls, of the corresponding AIMD simulations. All DeePMD trajectories for water and ice are 300 ps long and use the same time step of the AIMD simulations.

We use our own code to perform the constant temperature MD simulations of the organic molecules. In each DeePMD simulation the temperature is the same of that of the corresponding AIMD simulation. The time step and time length of the trajectories in these simulations are the same of those in the corresponding AIMD trajectories.

A.3 Additional Results

The radial distribution functions (RDFs) of ice Ih (b), (c) ,and (d) are reported in Figs. 6, 7, and 8, respectively.

The probability distribution function of the O-O bond orientation order parameter $Q_{6}$ is reported in Fig. 9. The bond orientation order parameter for oxygen $i$ , as proposed in Ref. Lechner and Dellago (2008), is defined by

In this work we take $r_{min}=0.31$ nm and $r_{max}=0.36$ nm.