Ab initio molecular dynamic study of solid-state transitions of ammonium nitrate

High-pressure polymorphism and phase transitions have wide ranging consequences on the basic properties of ammonium nitrate. However, the phase diagram of ammonium nitrate at high pressure and high temperature is still under debate. This study systematically investigates the phase transitions and structural properties of ammonium nitrate at a pressure range of 5–60 GPa and temperature range of 250–400 K by ab initio molecular dynamics simulations. Two new phases are identified: one corresponds to the experimentally observed phase IV’ and the other is named AN-X. Simultaneously, the lattice strains play a significant role in the formation and stabilization of phase IV’, providing a reasonable explanation for experimental observation of phase IV-IV’ transition which only appears under nonhydrostatic pressure. In addition, 12 O atoms neighboring the NH (N atom in ammonium cation) atom are selected as reference system to clearly display the tanglesome rotation of ammonium cation.

that of phase IV with subtle difference in the hydrogen-bonding network. Subsequently, Dunuwille et al. 21 further affirmed the phase IV-IV' transition using pressure-induced Raman spectral changes. Sorescu and Thompson studied the isotropic compression of phase IV using ab initio total energy calculations and found that phase IV was stable up to 600 GPa. 22 More recently, phase IV is considered to be stable up to 45 GPa at room temperature using synchrotron X-ray diffraction (XRD) and Raman spectroscopy measurements. 23 There are still numerous uncertainties in phase IV-IV' transition and higher pressure phases of AN.
Theoretical studies can make a detailed understanding of AN at the atomic level, such as melting point, energy bands, density of states, the nature of hydrogen bonding, and the structural and electronic properties. 22,[24][25][26] In this work, we use the supercell (Fig. 1a) of phase IV as initial structure to explore the phase transition of AN from 5 to 60 GPa through ab initio molecular dynamics simulations. The results shows that the phase IV' is confirmed and a new phase, AN-X, is discovered. In addition, we select 12 O atoms neighboring the N H (N atom in [NH 4 ] + ) as reference system to show the rotation details of [NH 4 ] + in different phases.

Results and Discussion
Before obtaining crystal structure from MD trajectories, it is indispensable to investigate the atomic vibration of a system, which can be done by calculating the atomic mean square displacement (MSD) as a function of time using: where N is the total number of atoms in the system, r i (t) is the atomic position at a time of t and r i (0) is the initial position of atom i, and < … > represents the ensemble average.
The MSD derived from MD simulations at 5 GPa and 300 K are shown in Fig. 2. It can be seen that the MSD of N and O atoms are approaching to a constant, indicating that they only slightly vibrate near the equilibrium position. So, the equilibrium positions of N and O atoms can be obtained by averaging their MD trajectories after the systems reach equilibrium state (about 5 ps). However, the MSD of H atoms can be clearly seen to fluctuate wildly with time, indicating that the above method doesn't work for hydrogen atoms. Observing the particles' trajectories over time, it is found that this trouble is caused by the rotation of [NH 4 ] + . Therefore, another method is adopted to get the equilibrium positions of H atoms. Firstly, H atoms are shifted into the supercell according to periodicity condition at every MD step to obtain the "H atom cloud" (see Fig. 3). Then, the positions of H atoms can be obtained by making statistical average over the "H atom cloud".
For convenience in discussion, the supercell is divided into four layers which are labeled as A, B, C and D layer, as shown in Fig. 1a. The "H atom cloud" in C layer along b-axis are shown in Fig. 3. In Fig. 3a, it is broadest near the middle and tapering toward both ends of the circled "H atom cloud". Moreover, the H atoms are symmetric about the dotted line between two N atoms. As the pressure increased to 20 GPa (see Fig. 3b), H atoms in the middle of circled "H atom cloud" gradually move to both ends, indicating that the vector of N-H bond points to one O atom more probably than the middle of two O atoms. At 30 and 60 GPa, the circled "H atom cloud" gather closely together and stay away from the dotted line, corresponding to two new structures, as shown in Fig. 3c,d.
Much more detailed information about the "H atom cloud" can be explored by examining the relative orientations of four N H -H i (i = 1, 2…4) bonds in [NH 4 ] + at every MD steps. Without loss of generality, we trace the ammonia cation in the central area of Fig. 3a  Three different primitive cells are obtained by making statistical average over the "H atom cloud", as shown in Fig. 3d-f. The first structure (Fig. 3d) is in agreement with the experimentally proposed AN-IV. The second structure (Fig. 3e), which is considered to be the metastable phase IV' observed in experiments, has a slightly distorted monoclinic lattice with P2 1 /m symmetry. It is the first time to uncover the third structure (Fig. 3f), which is defined as AN-X and identified as orthorhombic structure with the space group Pnma. The lattice parameters of AN-IV' and AN-X at 30 and 60 GPa are listed in Table 1, respectively. At the same time, the β angles of AN-IV, AN-IV' and AN-X as a function of pressure at 300 K are shown in Fig. 5a. From this graph, one can see that the β angle in AN-IV' is slightly smaller than 90°, which is different from AN-IV and AN-X. We speculate that the slightly change in β angle is the inducement of phase IV-IV' transition.
To investigate the effect of β angle in the phase IV-IV' transition at 300 K, both the simulations in the NPT and NVT ensemble are performed at approximately 20 GPa. In the NPT ensemble, the structure of AN-IV' is chosen as the initial structure, and then MD simulation is performed at 20 GPa and 300 K. For most of time, the crystal structures are AN-IV. However, in some time period, when mean value of β angles is slightly smaller than 90°, crystal structures transform to AN-IV' . For the NVT ensemble, the initial structure is obtained by adjusting the lattice parameters of AN-IV to be the same as AN-IV' (a = 10.59 Å, b = 9.23 Å, c = 9.05 Å, and β = 88.95°) at 20 GPa and keeping the relative positions of atoms in AN-IV invariant. By observing the snapshots during this NVT simulation, we find that the crystal structure is the same as AN-IV' . These two simulation results prove that some deviation from 90° in β angle indeed plays a crucial role in the formation of AN-IV' , which can explain why phase IV-IV' transition only appears under conditions of nonhydrostatic pressure.
To investigate the effect of lattice strains in the stabilization of AN-IV' at 300 K, we perform simulations in the NVT ensemble at approximately 50 GPa. Firstly, we optimize the structures of AN-IV' and AN-X at 0 K and 50 GPa. Then, we keep the relative position of atoms in the structures generated from the first step and exchange the lattice parameters of two phases at 300 K and 50 GPa. By this means, we obtain two kinds of initial structures for MD simulations: one has the relative position of atoms in AN-IV' at 0 K and the lattice parameters of AN-X at 300 K; another has the consistent relative position of atoms with AN-X at 0 K and the lattice parameters of AN-IV' at 300 K. Then, we perform MD simulations with above structures in the NVT ensemble. We discover that the configuration change from AN-IV' to AN-X in the first simulation, oppositely, there is few variation in the second simulation. Compare the results between these two simulations: Note that some deviation from 90° in β angle plays a crucial role in the stabilization of AN-IV' .
The phase IV-IV' transition mechanism can be further explored by observing the total energy of AN-IV and AN-IV' as a function of the rotational angle at different pressures (see Fig. 5b). In the calculations, we keep the   orientations of ammonium cations in the same layer at high pressure. Then the crystal structure will transform from AN-IV to AN-IV' .
To demonstrate that the structure of AN-IV' proposed in this work is the metastable phase IV' found in experiment, we simulate the Raman spectra of AN-IV and AN-IV' at 21.7 GPa (Fig. 6a). Because the orientations of [NH 4 ] + is a major difference between AN-IV and AN-IV' in our simulations, the Raman peaks involved [NH 4 ] + are particularly significant. It is noticed that the ν ' 2 ([NH 4 ] + ) mode at around 1730 cm −1 , which is used to confirm the formation of the phase IV' in the experiment 11 , dose exist in our simulated Raman spectra of AN-IV' at 21.7 GPa. However, there is no ν ' 2 ([NH 4 ] + ) mode in our calculated Raman spectra of AN-IV. In addition, the new peak for R([NO 3 ] − ), which is observed in the low frequency region around 400 cm −1 by experiment, appears in our calculated Raman spectra of AN-IV' . Highly agreement between the theoretical and experimental results provides a powerful support for the existence of the AN-IV' .
At first glance, the monoclinic structure of AN-IV' (P2 1 /m) proposed in the present work is different from the experimentally suggested orthorhombic structure (Pmmn). It's worth noting that the XRD can't confirm the position of H atoms and the β angle of our proposed AN-IV' is extreme approaching 90° (see Fig. 5a). If the H atoms in our proposed AN-IV' are eliminated, it has a same symmetry (Pmmn) with the experimentally suggested structure. Therefore, it is convinced that the phase AN-IV' proposed in this study is the metastable phase IV' observed in experiments.
Geometry optimizations are performed with full structural relaxation including atomic positions and lattice constants at 0 K to obtain the enthalpy difference curves of AN-IV' and AN-X relative to AN-IV (Fig. 6b). AN-IV' and AN-X have lower energy than AN-IV, suggesting that AN-IV' and AN-X are more stable than AN-IV in the pressure range of this study. This result is different from the experimental results that phase IV transformed into phase IV' at about 20 GPa and room temperature. However, extending the phase boundaries of our MD simulations to low temperature and low pressure (see Supplementary Figure S2 online), one can see that AN-IV' is stable until 0 GPa and 0 K. This is consistent with the results of calculated enthalpies at 0 K. Therefore, the contradiction between calculated (enthalpies) and measured stable pressure region of phase IV can be attributed to temperature effect. The inset in Fig. 6b displays the enthalpy difference of AN-IV' and AN-X at 0 K, showing that AN-X is more stable than AN-IV' above 35 GPa.
The simulation protocol, together with majority of the simulation results for AN, is schematically shown on Fig. 7. At 300 K, phase IV-IV' transition pressure is approximately 20 GPa, which is almost identical to the experimentally observed 18 GPa. During MD simulation at 43 GPa and 400 K, the mutual transformation between AN-IV' and AN-X are observed, suggesting that 43 GPa are very close to the boundary of phase IV′ -X transition at 400 K. There are also yellow spots (AN-IV′ ) in the red area (AN-X), which can be explained by Fig. 5b. In this graph, it is obvious that the total energies of AN-IV' and AN-X increase sharply with the change of angle at 30 GPa, indicating that the spatial orientation of ammonium cation will be strongly restrained when the pressure is over 30 GPa. Thus, once the AN-IV' formed, it is difficult to transform into AN-X at low temperature.

Conclusion
We explore the high-pressure phase transition of ammonium nitrate through ab initio molecular dynamics simulations. The existence of phase IV' is confirmed and a new phase, AN-X, is discovered. In addition, through the NVT ensemble simulations, we reveal that lattice strains play an important role in formation and stabilization of phase IV' . The above fact can explain the proposed conclusion in experiment that phase IV-IV' transition is induced by shear stress under nonhydrostatic pressure. Finally, the phase diagram of AN is determined at pressure range of 5-60 GPa and temperature range of 250-400 K.

Methods
In this study, we use ab initio molecular dynamics simulations based on density functional theory implemented in the Vienna ab initio simulation package (VASP) code [27][28][29] to schematically investigate AN. The projector-augmented wave (PAW) method 30,31 is adopted and exchange-correlation functions are treated within the Perdew-Burke-Ernzerhof generalized gradient approximation (GGA) 32 . A plane-wave basis-set cutoff energy  of 520 eV and 2 × 2 × 2 Monkhorst− Pack k-point mesh 33 are employed. Simulations are implemented both in the NPT (N-constant number of particles, P-constant pressure, and T-constant Temperature) ensemble 34 and NVT (V-constant volume) ensemble 35-37 for a system containing 16 AN molecules.
Considering the computational efficiency and the feature of MD, the simulations are performed in two ways. First, we optimize the supercell of AN-IV (Fig. 1) at 0 K and a series of pressures (5, 20, 40 and 60 GPa) as initial structures. We heat the structures at a series of temperatures (250, 300, 350 and 400 K) using their corresponding pressures. All the simulation conditions and results using this method are shown by small circles without cross lines in Fig. 7. Second, we choose the results of finished simulations as initial structures and perform MD simulations at new conditions. For example, at temperature of 300 K, we use the result of 5 GPa as initial structure for the simulation at 10 GPa, then we gradually increase pressure in 10 GPa steps until 60 GPa. In a similar manner, we use the result of 60 GPa as initial structure and decrease pressure to 15 GPa in the interval pressure of 10 or 5 GPa at 300 K. The whole simulation process of increasing and decreasing pressure is shown on Fig. 7 and marked as blue arrows. In addition, we also show the annealing process by the orange arrows, which are signed by small circles with cross lines. In addition, we employ simulation times of 20-30 ps (with time step of 1 fs) at the vicinity of the phase boundary, and at least 10 ps at the other conditions.
We employ density functional perturbation theory calculations with CASTEP (Cambridge Serial Total Energy Package) 38 code to study the Raman spectra of AN-IV and AN-IV' at 21.7 GPa. GGA-PBE for the exchange-correlation functional, norm-conserving potentials, the cutoff energy of 750 eV and the Monkhorst-Pack k-points of 5 × 6 × 5 for corresponding Brillouin zone have been applied.