Received 1 October 2013
Crystal structure prediction and hydrogen-bond symmetrization of solid hydrazine under high pressure: a first-principles study
Hua-Di Zhang,a Song-Kuan Zheng,a Xi-Lian Jin,a Shu-Qing Jiang,a Zhi He,a Bing-Bing Liua and Tian Cuia*
In this article, the crystal structure of solid hydrazine under pressure has been extensively investigated using ab initio evolutionary simulation methods. Calculations indicate that hydrazine remains both insulating and stable up to at least 300 GPa at low temperatures. A structure with P21 symmetry is found for the first time through theoretical prediction in the pressure range 0-99 GPa and it is consistent with previous experimental results. Two novel structures are also proposed, in the space groups Cc and C2/c, postulated to be stable in the range 99-235 GPa and above 235 GPa, respectively. Below 3.5 GPa, C2 symmetry is found originally, but it becomes unstable after adding the van der Waals interactions. The P21Cc transition is first order, with a volume discontinuity of 2.4%, while the CcC2/c transition is second order with a continuous volume change. Pressure-induced hydrogen-bond symmetrization occurs at 235 GPa during the CcC2/c transition. The underlying mechanism of hydrogen-bond symmetrization has also been investigated by analysis of electron localization functions and vibrational Raman/IR spectra.
Hydrazine is a binary inorganic chemical compound with the formula N2H4 that has been widely used in various rocket fuels and fuel cells and in the organic chemical industry. However, there is some doubt as to the crystal structure of solid hydrazine at atmospheric pressure. In 1951, Collin & Lipscomb (1951) indicated that the symmetry group is P21/m with two molecules in the unit cell. Busing et al. (1961) reported that the symmetry is P21 with two N2H4 molecules per unit cell. Baglin et al. (1967) analyzed the optical modes of the crystal lattice in the far-IR region and pointed out that the crystal symmetry should be the P21 space group. Then in 1969, Guay & Savoie (1969) analyzed the Raman spectra of solid hydrazine and concluded that the symmetry group was one of the C2h space groups with four N2H4 molecules in the unit cell. In 1975, Durig et al. (1975) analyzed the Raman spectra of solid hydrazine and hydrazine-d4 and found themselves in agreement with Guay and Savoie's suggestion. However, they pointed out that the number of molecules per primitive cell might be larger than four. Compared with the structure of solid hydrazine at atmospheric pressure, the high-pressure structure has been less well investigated and still remains unknown. In this paper, we conclude that the phase in space group P21 is the correct phase at atmospheric pressure, consistent with the findings of Busing and co-workers. Beyond that, we also predict the high-pressure structure of solid hydrazine using an ab initio evolutionary algorithm (Oganov & Glass, 2006; Glass et al., 2006) and first-principles density functional theory (DFT).
Solid hydrazine is a hydrogen-bonded molecular compound. Hydrogen bonds are typically weaker than covalent, ionic and metallic bonds, but stronger than van der Waals interactions. Hydrogen-bond symmetrization has been found in hydrogen-bonded molecular compounds such as H2O, NH3, HF, HCl and HBr (Pickard & Needs, 2008; Duan et al., 2010; Zhang et al., 2010). For example, Pickard & Needs (2008) reported that half of the hydrogen bonds between NH2- anions in solid ammonia become symmetric at 331 GPa. In this paper, we investigate high-pressure-induced hydrogen-bond symmetrization using first-principles DFT.
Searches for stable crystal structures (with the lowest Gibbs free energy) were performed using the evolutionary algorithm USPEX, as implemented in the USPEX program (Oganov & Glass, 2006; Glass et al., 2006). In this approach, no initial experimental information is required except for chemical composition, and thus it provides fully non-empirical crystal structure predictions. The first-generation structures were produced randomly. Each subsequent generation was produced from 50% of the lowest-enthalpy structures of the preceding generation and the lowest-enthalpy structure always survived into the next generation. The variation operators used for producing offspring included heredity (40% of structures), atomic permutation (20%), atom-position mutations (20%) and lattice mutations (20%).
Geometry optimizations were performed within the framework of the generalized gradient approximation (GGA) with Perdew-Burke-Ernzerhof (Paier et al., 2005) parameterization for the exchange-correlation functional for DFT using the projector-augmented wave (PAW) method (Kresse & Joubert, 1999), as implemented in the Vienna ab initio simulation package (VASP) (Kresse & Furthmüller, 1996). For the searches, we used a plane-wave cut-off energy of 520 eV, and a Monkhorst-Pack (Monkhorst & Pack, 1976) Brillouin zone sampling grid of spacing 2 × 0.06 Å-1. The structures obtained in the searches were re-relaxed at a higher level using a plane-wave cut-off energy of 800 eV and a Brillouin zone sampling grid of spacing 2 × 0.025 Å-1. The total energy was minimized until the energy convergence became less than 1 × 10-8 eV. To obtain the equilibrium structures of the unit cells at a given pressure, the internal atomic positions were optimized until the residual forces became less than 1 × 10-3 eV Å-1.
The phonon calculations were carried out using a supercell approach, implemented in the Phonopy program (Togo et al., 2008). A convergence check gave the use of the 4 × 3 × 3, 3 × 3 × 4 and 3 × 3 × 4 supercells for the P21, Cc and C2/c structures, respectively. A cut-off energy of 800 eV was used consistently throughout all phonon calculations. Brillouin zone integrations were carried out using only the point, after more complete k-point samplings had been tested showing no significant differences from the -point only results.
The first-principles calculations of vibrational Raman/IR spectra were performed using a second-order response (Lazzeri & Mauri, 2003) to DFT, as implemented in the QUANTUM-ESPRESSO package (Giannozzi et al., 2009). We used the local-density approximation for the exchange-correlation functional, norm-conserving pseudopotentials and a plane-wave expansion up to a 180 Ry cut-off. Brillouin zone sampling was performed on 13 × 7 × 9 Monkhorst-Pack meshes for the Cc and C2/c structures.
We explored the crystal structures of solid hydrazine using the evolutionary algorithm USPEX with one to four N2H4 formula units per cell at 5, 40, 100, 200 and 240 GPa. Analysis of the predicted structures gave us a shortlist of candidate structures in the space groups C2 (two molecules per cell), Aba2 (four molecules per cell), P21 (two molecules per cell), Cc (four molecules per cell), P21/m (four molecules per cell), Pm (four molecules per cell), P21/c (two molecules per cell) and C2/c (four molecules per cell). We then re-calculated the enthalpies of the most stable phases at a larger number of pressures ranging from 0 to 300 GPa, generating the data shown in Fig. 1. In the range from 0 to 300 GPa, C2, P21, Cc and C2/c are the most favourable structures, with lowest enthalpies below 3.7 GPa, in the range 3.7-100.5 GPa, in the range 100.5-235 GPa and above 235 GPa, respectively. We also calculated the decomposition (N2 + 2H2) enthalpies by adopting the , P63/m, C2/c and structures (Kohanoff et al., 1997; Johnson & Ashcroft, 2000; Pickard & Needs, 2007) for H2, and the P41212, I213, Pba2 and structures (Ma et al., 2009; Wang et al., 2012; Pickard & Needs, 2009) for N2. The decomposition enthalpies are at least 0.5 and 2.0 eV higher than that of the P21 structure in the pressure ranges 0-20 and 20-300 GPa, respectively, so the enthalpies of N2 + 2H2 are not shown in Fig. 1. This eliminates the possibility that N2H4 decomposes into N2 and H2 up to at least 300 GPa.
| || Figure 1 |
Calculated enthalpies per N2H4 unit of various structures as a function of pressure. Differences in enthalpy from the P21 phase are plotted in the pressure ranges (a) 0-20 GPa and (b) 20-300 GPa. Insets: enthalpies with zero-point motion corrections, relative to the P21 phase. (c) Differences in enthalpy after adding the van der Waals interactions to the P21 and C2 phases from 0 to 6 GPa.
The favoured structures of solid hydrazine obtained in this work are shown in Fig. 2, and the corresponding parameters at selected pressures are listed in Table 1. The C2 structure consists of two formula units (i.e. N4H8) in a monoclinic crystal lattice, as shown in Fig. 2(a). The single N-N bond length is 1.454 Å at 1 GPa. The P21 structure also consists of two formula units in a monoclinic crystal lattice, as shown in Fig. 2(b). The single N-N bond length is 1.440 Å at 8 GPa. The Cc structure consists of four formula units in a monoclinic crystal lattice, as shown in Fig. 2(c). The single N-N bond length is 1.362 Å at 120 GPa, which is longer than the value of 1.335 Å for the single N-N bond in the I213 structure of solid nitrogen at the same pressure. The C2/c structure consists of four formula units in a monoclinic crystal lattice, as shown in Fig. 2(d). The single N-N bond length is 1.302 Å at 300 GPa, which is longer than the value of 1.259 Å for the single N-N bond in the structure of solid nitrogen at the same pressure.
| || Figure 2 |
The structures of solid hydrazine under high pressure. (a) The C2 phase at 1 GPa, (b) the P21 phase at 8 GPa, (c) the Cc phase at 120 GPa and (d) the C2/c phase at 300 GPa. Large red atoms are N and small green atoms are H.
It is known that quantum effects related to H atoms can be very important. In particular, due to the low mass of the H atom, the zero-point (ZP) energy (Kim et al., 2008) is expected to be large and might significantly alter the structural stability. We have thus estimated the ZP vibrational energies for the C2, P21 and Cc structures, and the results are plotted in the inset of Fig. 1. The ZP motion corrections only lower the transition pressure from 3.7 to 3.5 GPa for the C2P21 phase transition and from 100.5 to 99.0 GPa for the P21Cc phase transition.
Van der Waals interactions between atoms and molecules play an important role in many chemical systems (Grimme, 2006). However, a general drawback of all common GGA functions is that they cannot describe van der Waals interactions correctly. The enthalpies of the C2 and P21 structures are very similar, so the van der Waals interactions become significant. We have estimated the van der Waals interactions of the C2 and P21 structures from 0 to 6 GPa (Fig. 1c). The enthalpy of the P21 phase is lower than that of the C2 phase, so we believe that C2 may be only a metastable phase at low pressure.
An essential requirement for a stable structure is lattice dynamic stability. We thus calculated the phonon dispersions for the P21, Cc and C2/c phases at selected pressures, as shown in Fig. 3. It is quite clear that these three structures are dynamically stable, due to the absence of any imaginary phonons in the whole Brillouin zone. The mechanical stability of a structure provides a useful insight into the stability of the crystals. To evaluate the mechanical stability of these three phases, elastic constants were calculated and these are listed in Table 2. As shown in Table 2, the elastic constants of these three structures satisfy the mechanical stability criteria of Kresse & Furthmüller (1996), indicating that these three structures are mechanically stable.
| || Figure 3 |
Calculated phonon dispersion curves for (a) the C2 phase at 1 GPa, (b) the P21 phase at 8 GPa, (c) the Cc phase at 120 GPa and (d) the C2/c phase at 300 GPa.
The calculated equations of state shown in Fig. 4 suggest that the P21Cc phase transition is of the first-order reconstructive type with a volume discontinuity of 2.4%, while the CcC2/c transition is characterized as the second-order displacive type without a volume discontinuity. We have also calculated the density of electronic states (DOS) for the P21, Cc and C2/c structures at 80, 200 and 300 GPa, respectively, as shown in Fig. 5. The band gaps of these three structures are 5.2, 3.8 and 3.4 eV, respectively. In fact, the calculated band gaps of these three structures in the pressure range 0-300 GPa are always larger than 3.4 eV, and the real gap is probably even larger, as the PBE density functional tends to underestimate band gaps (Pickard & Needs, 2009). Thus, these calculations indicate that, at low temperatures, solid hydrazine remains insulating up to at least 300 GPa.
| || Figure 4 |
Volume plotted as a function of pressure for the P21, Cc and C2/c structures. Insets: the volume discontinuity at the pressures of (a) the P21Cc and (b) the CcC2/c phase transitions.
| || Figure 5 |
Calculated DOS projected on N and H atoms for (a) the C2 phase at 3 GPa, (b) the P21 phase at 80 GPa, (c) the Cc phase at 200 GPa and (d) the C2/c phase at 300 GPa.
As shown in Fig. 2(d), we find that half of the hydrogen bonds become symmetric in the C2/c structure, which is stable above 235 GPa. In an attempt to present a more transparent picture of hydrogen-bond symmetrization during the CcC2/c transition, we have calculated the variation in covalent bond lengths and hydrogen-bond lengths for the two structures under pressure, depicted in Fig. 6. As for the Cc phase, the intramolecular covalent bonds (H2-N1) elongate gradually and the intermolecular hydrogen bonds (H2N4) shrink gradually with increasing pressure. This clearly indicates that, under compression, the covalent bonds are weakened while the hydrogen bonds are strengthened. When the covalent bonds and hydrogen bonds become equal, hydrogen-bond symmetrization occurs at about 235 GPa. In fact, the Cc structure naturally takes the form of the C2/c structure after hydrogen-bond symmetrization. It is worth noting that the critical pressures corresponding to hydrogen-bond symmetrization almost perfectly match the transition pressure derived from the enthalpy calculations (Fig. 1b), further implying that the CcC2/c transition originates from pressure-induced hydrogen-bond symmetrization. The continuous change in structural parameters (Fig. 6), the features of the enthalpy curves (Fig. 1b) and the continuous volume change at the transformation point (Fig. 4c) indicate that the CcC2/c phase transition is of a second-order nature.
| || Figure 6 |
Evolution of bond lengths with pressure for covalent bonds (H2-N1) and hydrogen bonds (H2N4) in the Cc and C2/c structures. Insets: the atomic positions of atoms H2, N1 and N4 in (a) the Cc and (b) the C2/c structures.
To gain further insight into the origin of pressure-induced hydrogen-bond symmetrization in solid hydrazine, we performed an additional analysis of electron localization functions (Becke & Edgecombe, 1990) (ELFs), as shown in Fig. 7. The topological analysis of ELFs is widely used to characterize the degree of electron localization and the nature of chemical bonds in molecules and solids. Specifically, ELF = 1 corresponds to an extreme localization, while ELF = 0.5 reflects the behaviour of a homogeneous electron gas. As noted above, the intermolecular hydrogen bonds (H2N4) are weaker than the intramolecular covalent bonds (H2-N1) in the Cc structure at 200 GPa. Upon compression, as atoms H2 move toward the mid-point between N1 and N4, electrons are partially transferred from the intramolecular to the intermolecular region, which results in weakened covalent bonds and strengthened hydrogen bonds. After the CcC2/c transition, the bonds between the H atoms and the two N atoms form symmetric bonds.
| || Figure 7 |
Electron localization function maps for (a) the C2 phase at 3 GPa, (b) the P21 phase at 80 GPa, (c) the Cc phase at 200 GPa and (d) the C2/c phase at 300 GPa.
We have also calculated the vibrational Raman/IR spectra of the Cc and C2/c structures at 200 and 300 GPa, respectively. The vibrational frequencies of Raman/IR spectra can be obtained by calculating optical phonons at the zone centre (the point), as shown in Table 3. As for the Cc phase, all the optical modes (17A'' + 16A') are both Raman and IR active, but for the C2/c phase, each of the optical modes (7Ag + 8Bg + 9Au + 9Bu) is either Raman or IR active, i.e. Bg and Ag are Raman active, while Bu and Au are IR active. According to IR and Raman selection rules (Duval, 1992), the C2/c structure is centrosymmetric because the Raman and IR active modes are mutually exclusive, while the Cc structure is noncentrosymmetric, further implying that the C2/c phase is the structure with hydrogen-bond symmetrization. The H atoms that form symmetric bonds are at exactly the symmetry centres of the C2/c structure.
In summary, we have extensively investigated the crystal structure of solid hydrazine under pressure using ab initio evolutionary simulation methods. Our calculations indicate that, at low temperatures, it remains both insulating and stable up to at least 300 GPa. The P21 phase is confirmed as a good candidate for the low-pressure phase, and the Cc and C2/c structures have been found to be stable in the pressure range 99-235 GPa and above 235 GPa, respectively. The P21Cc phase transition in hydrazine is of the first-order reconstructive type with a volume discontinuity of 2.4%, while the CcC2/c transition is of the second-order displacive type without a volume discontinuity. The C2/c phase is found to be a structure with symmetric hydrogen bonds and hydrogen-bond symmetrization occurs at 235 GPa. Analysis of the ELFs and vibrational Raman/IR spectra provides a deep insight into the underlying mechanism of hydrogen-bond symmetrization.
This work was supported by the National Basic Research Program of China (grant No. 2011CB808200), the Program for Changjiang Scholars and Innovative Research Team in University (grant No. IRT1132), the National Natural Science Foundation of China (grant Nos. 51032001, 11074090, 10979001, 51025206 and 11174102) and the National Foundation for Fostering Talents of Basic Science (grant No. J1103202). Some of the calculations were performed at the High Performance Computing Center (HPCC) of Jilin University.
Baglin, F. G., Bush, S. F. & Durig, J. R. (1967). J. Chem. Phys. 47, 2104-2108.
Becke, A. D. & Edgecombe, K. E. (1990). J. Chem. Phys. 92, 5397-5403.
Busing, W. R., Zocchi, M. & Levy, H. A. (1961). Program of the Annual Meeting of the American Crystallographic Association. Paper N-3, August.
Collin, R. L. & Lipscomb, W. N. (1951). Acta Cryst. 4, 10-14.
Duan, D., Tian, F., He, Z., Meng, X., Wang, L., Chen, C., Zhao, X., Liu, B. & Cui, T. (2010). J. Chem. Phys. 133, 074509.
Durig, J. R., Griffin, M. G. & Macnamee, R. W. (1975). J. Raman Spectrosc. 3, 133-141.
Duval, E. (1992). Phys. Rev. B, 46, 5795-5797.
Giannozzi, P., et al. (2009). J. Phys. Condens. Matter, 21, 95502.
Glass, C. W., Oganov, A. R. & Hansen, N. (2006). Comput. Phys. Commun. 175, 713-720.
Grimme, S. (2006). J. Comput. Chem. 27, 1787-1799.
Guay, M. & Savoie, R. (1969). Can. J. Chem. 47, 201-208.
Johnson, K. A. & Ashcroft, N. W. (2000). Nature, 403, 632-635.
Kim, D. Y., Scheicher, R. H., Lebègue, S., Prasongkit, J., Arnaud, B., Alouani, M. & Ahuja, R. (2008). Proc. Natl Acad. Sci. USA, 105, 16454-16459.
Kohanoff, J., Scandolo, S., Chiarotti, G. L. & Tosatti, E. (1997). Phys. Rev. Lett. 78, 2783-2786.
Kresse, G. & Furthmüller, J. (1996). Phys. Rev. B, 54, 11169-11186.
Kresse, G. & Joubert, D. (1999). Phys. Rev. B, 59, 1758-1775.
Lazzeri, M. & Mauri, F. (2003). Phys. Rev. Lett. 90, 036401.
Ma, Y., Oganov, A. R., Li, Z., Xie, Y. & Kotakoski, J. (2009). Phys. Rev. Lett. 102, 065501.
Monkhorst, H. J. & Pack, J. D. (1976). Phys. Rev. B, 13, 5188-5192.
Oganov, A. R. & Glass, C. W. (2006). J. Chem. Phys. 124, 244704.
Paier, J., Hirschl, R., Marsman, M. & Kresse, G. (2005). J. Chem. Phys. 122, 234102.
Pickard, C. J. & Needs, R. J. (2007). Nature Phys. 3, 473-476.
Pickard, C. J. & Needs, R. J. (2008). Nature Mater. 7, 775-779.
Pickard, C. J. & Needs, R. J. (2009). Phys. Rev. Lett. 102, 125702.
Togo, A., Oba, F. & Tanaka, I. (2008). Phys. Rev. B, 8, 134106.
Wang, X., Wang, Y., Miao, M., Zhong, X., Lv, J., Cui, T., Li, J., Chen, L., Pickard, C. J. & Ma, Y. (2012). Phys. Rev. Lett. 109, 175502.
Zhang, L., Wang, Y., Zhang, X. & Ma, Y. (2010). Phys. Rev. B, 82, 014108.