A flexible-molecule force field to model and study hexanitrohexaazaisowurtzitane (CL-20) – polymorphism under extreme conditions

Abstract

The quantum-chemistry based force field (FF) developed for HMX by Smith and Bharadwaj (SB) [G. D. Smith and R. K. Bharadwaj, J. Phys. Chem. B, 1999, 103(18), 3570–3575] is transferred to another nitramine of different stoichiometry: hexanitrohexaazaisowurtzitane (CL-20 or HNIW). The modification of a single parameter alongside a very small number of add-ons related to carbon–carbon bonds, angles and dihedrals lead to two SB FF variants denoted SB-CL20 and SB-CL20 + CCNN. These flexible-molecule FFs should inherit the predictive capabilities of SB FF. For this purpose, we perform Molecular Dynamics simulations at ambient temperature and selected pressures. The modeled structures of the various CL-20 polymorphs are consistent with experimental data. Focusing on the ε-polymorph, we determine an equation of state which consolidates the general trend underpinned by most published results, and we confirm the increasing stiffness of the crystal under pressures up to 90 GPa. Moreover, we link some subtle pressure-induced changes of the elastic and structural properties to the flexibility and mobility of well-identified nitro groups. Finally, the simulations of the γ ↔ ζ phase transition suggest different multiple-step direct and reverse thermodynamic paths.

---

1. Introduction

hexanitrohexaazaisowurtzitane (denoted HNIW or CL-20) is a polycyclic nitramine with six nitro groups developed as a highly energetic, compact and efficient explosive for industrial and military applications. With the formula C₆H₆O₁₂N₁₂, it has five polymorphs: α, β, γ, ε¹ and high-pressure ζ.² Fig. 1 shows the phase transformation diagram using γ-polymorph as starting material (data from ref. 3), and the corresponding molecular conformations and crystal symmetries. Due to the large size of the molecule, the phase transformations can be challenging, with still unknown mechanisms and rates ranging from “slow” to “sluggish” or even “one-directional”. The quickest one is the reversible γ ↔ ζ transformation occurring at 0.7 GPa. Due to thermal decomposition upon heating (above 260 °C, with a melting point depending on pressure), there is no known liquid phase. At ambient conditions, γ-phase can be considered the stable one since ε-phase transforms to γ upon heating at 120–140 °C, and ε cannot be subsequently recovered when returned to ambient conditions.⁴ α-phase is stabilized by its ability to trap small molecules like H₂O,¹ CO or CO₂.³ Among α, β and ε metastable phases, ε is considered the most stable and β the least one.^5,6

Fig. 1 Phase transformation diagram (left) using γ-polymorph as starting material (data from ref. 3), corresponding molecular conformations showing the characteristic orientations of the nitro groups, and (right) crystal symmetries.

Experimental studies to understand the properties of ε-CL20 under pressure have been published but they show high discrepancies. Ciezak et al.⁷ have interpreted some Raman and Far-Infrared spectra changes as an ε → γ phase transition, which is not supported by the more recent XRD works from Millar⁸ and Gump and Peiris.⁹ The equations of state determined by XRD experiments by Millar,⁸ Gump and Peiris⁹ and Pinkerton¹⁰ are very different, especially the bulk modulus derivative which defines the high-pressure behavior. To reduce this lack of agreement, one can rely on theoretical studies of CL-20, using ab initio methods like DFT,^6,11–13 DFT-D,^10,14 DFTB,¹⁵ ab initio Molecular Dynamics (AIMD),¹⁶ and Molecular Dynamics (MD) using various force fields (FF) such as Sorescu–Rice–Thompson SRT (rigid molecules),¹⁷ COMPASS (flexible molecules, but not free-to-use)¹⁸ or ReaxFF (a reactive force field with no assumption on molecule topologies).^19,20 Resulting data are often model-sensitive, but the point is that both experimental and theoretical findings bolster a general trend, a consensus to be consolidated.

In order to study phase transformations, elastic and plastic deformations and shock-induced modifications of CL-20 in inert regimes, a non-reactive flexible-molecule FF is required. Due to the transferability of SRT FF over many nitramine crystals,²¹ its intermolecular parameters have been merged to the generic intramolecular topology of the general Amber FF (GAFF), yielding Amber-SRT FF.²² However the crystal densities of such modeled nitramines are systematically underestimated. Agrawal et al.²³ have reported at least for RDX that the orientation of the nitro groups relative to RDX ring significantly deviates from experimental data, altering in turn the crystal symmetry, and eventually requiring an improvement of the torsion and dihedral parameters. Fortunately, these parameters are correctly trained in the quantum-chemistry based force field of Smith & Bharadwaj (SB) developed for HMX,^24,25 and moreover remarkably transferable to RDX.^26,27

So the purpose of this paper is to present two FFs denoted SB-CL20 and SB-CL20 + CCNN, as two extensions of the SB FF to be used with CL-20. To achieve this transfer, only one parameter is modified, based on stoichiometry considerations, and a very few other missing ones are added. These FFs are then assessed using MD simulations at ambient temperature. We first compare the modeled and experimental structures of the various CL-20 polymorphs. Then, focusing on ε-phase, we determine its equation of state, study the evolution of the elastic properties with pressure and the induced modifications of the molecular conformation. Finally, the γ ↔ ζ phase transition being the quickest among the reported ones for CL-20, we check the ability of both FFs to reproduce it.

2. Computational methods

2.1. SB-CL20 force field

Smith & Bharadwaj²⁴ have developed a quantum-chemistry based force field (SB FF) to model HMX (C₄H₈O₈N₈, Fig. 2-left) with flexible molecules, where atoms have fixed partial atomic charges. It includes energy terms for bonds, angles, torsions, dihedrals and non-bonded Buckingham and Coulomb interactions. For CL-20, we use hereafter a cutoff radius of 8 Å for the Buckingham part and 10 Å the Coulomb part. For efficiency, the latter is evaluated using the damped shifted force method,²⁸ with a damping parameter of 0.25 Å⁻¹.

Fig. 2 Molecules of HMX (left) and CL-20 (middle), and CCNN dihedral interaction for SB-CL20 FF (right).

Due to the same stoichiometry, SB FF is straightforwardly transferable to RDX (C₃H₆O₆N₆), but not to CL-20 (C₆H₆O₁₂N₁₂, Fig. 2-middle). Indeed, each carbon atom in RDX/HMX binds to two hydrogen atoms, whereas it binds to only one in CL-20 molecule. In order to restore electroneutrality, we simply subtract one opposite hydrogen charge to the carbon one, yielding q^CL20_C = −0.2160e. The charges used for the other atoms are the ones prescribed by Smith & Bharadwaj in their original paper.²⁴ The Buckingham parameters remain unchanged.

To make SB FF work with CL-20 topology, we add the C–C bond parameters k = 536 kcal mol⁻¹ and r₀ = 1.53 Å, corresponding to harmonic bond energy U = k(r − r₀)²/2. The angular energy term has the harmonic form U = k(θ − θ₀)²/2. We choose to set C–C–N angle parameters equal to C–N–C in SB FF and we add the H–C–C parameters k = 75 kcal mol⁻¹ and θ₀ = 110.70°. These bond and angle parameters are simply taken from OPLS-AA FF.²⁹

We also propose an alternative version of SB-CL20 FF including CCNN dihedral interaction. Starting from the functional form proposed in OPLS-AA FF through eqn (1), we refine the three k parameters using GULP code (General Utility Lattice Program).³⁰ The optimal fit of the k parameters simultaneously minimizes the structural errors in β-, γ- and ε-CL20. The training set includes the experimental lattice parameters and atomic positions of the three selected polymorphs, using single cells and the respective crystal symmetries. Using SB-CL20 FF, this minimization (BFGS algorithm) yields the refined parameters presented in Table 1. They result in a dihedral interaction less constrained than with OPLS-AA FF, as shown in Fig. 2-right with a larger well in the range [90°, 180°] U [−180°, −90°], but steeper slopes within [−90°, 90°].

Table 1 Parameters for CCNN dihedral interaction

n	k (kcal mol^−1), refined for CL-20 FF	k (kcal mol^−1), OPLS-AA FF	d
1	8.936	11.034	1
2	−3.872	−0.968	−1
3	1.736	0.270	1

---

2.2. Molecular dynamics

To assess SB-CL20 FF, we perform MD simulations with LAMMPS code (large-scale atomic/molecular massively parallel simulator).³¹ We use a velocity-verlet integrator with a time step of 0.2 fs, periodic boundary conditions (PBC) and the NPT Parrinello–Rahman algorithm (constant number of atoms, controlled temperature and hydrostatic pressure with damping parameters of 10 fs and 100 fs, respectively). The unit cells of the CL-20 polymorphs have been replicated 13 × 9 × 5, 13 × 9 × 9, 9 × 15 × 9, 15 × 9 × 9 and 9 × 15 × 9 times, respectively for α, β, γ, ε and ζ, so that the simulation box contain 4680 CL-20 molecules for α, 4212 for β and 4860 for γ, ε and ζ. Note that we use the anhydrous form of α polymorph. Further specific details concerning the determination of the elastic properties and the equation of state of ε-CL20 are given in the following sections. The FF parameter files and some configuration files for LAMMPS are provided in ESI.†

3. Results and discussions

3.1. Lattice parameters

Table 2 presents the lattice parameters of the five CL-20 polymorphs modeled at 300 K and 1 atm, or 3.3 GPa for ζ. While the energy and structure convergence for α, γ, ε and ζ phases is achieved after 10 ps, a longer time of 160 ps is needed for β-phase. This does not seem anomalous since amidst the metastable polymorphs, β is the least stable one.^5,6 After this convergence is achieved, every lattice parameter is averaged for an additional 10 ps for α, γ, ε and ζ phases, and 40 ps for β-phase (sampling every 0.2 ps). Note that using the charges prescribed by Bedrov et al.²⁵ leads to β-phase instability (disordered molecular orientations). This is the reason of using original SB charges. All the converged structures are in good agreement with experimental data. The main improvement with CCNN dihedral add-on is on the “b” parameter of γ-phase, with a deviation from the experimental value lowered from 7.4% to 5.6%. The implication to the γ ↔ ζ phase transition is addressed in Section 3.5.

Table 2 Lattice parameters of CL-20 polymorphs modeled with SB-CL20 force field and with CCNN dihedral add-on, at 300 K and 1 atm (90° angles below 0.1% deviation are not shown)

CL20		a (Å)	b (Å)	c (Å)	α (°)	β (°)	γ (°)	ρ (g cm^−3)
α^anhydrous	exp.^1	9.485	13.225	23.673				1.961
	SB-CL20	9.313	12.969	24.237				1.989
	%dev	−1.8	−1.9	+2.4				+1.4
	+CCNN	9.426	13.001	23.875				1.990
	%dev	−0.6	−1.7	+0.9				+1.5
β	exp.^1	9.676	13.006	11.649				1.985
	SB-CL20	9.834	12.691	12.760				1.983
	%dev	+1.6	−2.4	+1.0				−0.1
	+CCNN	9.683	12.572	12.101				1.976
	%dev	+0.1	−3.3	+3.9				−0.5
γ	exp.^1	13.231	8.170	14.876		109.17		1.916
	SB-CL20	13.752	7.563	15.179		108.71		1.947
	%dev	+3.9	−7.4	+2.0		−0.4		+1.6
	+CCNN	13.611	7.714	15.163		109.03		1.934
	%dev	+2.9	−5.6	+1.9		+0.1		+0.9
ε	exp.^1	8.852	12.556	13.386		106.82		2.044
	SB-CL20	8.851	12.268	13.484		104.93		2.057
	%dev	−0.0	−2.3	+0.7		−1.8		+0.6
	+CCNN	8.951	12.284	13.437		104.95		2.039
	%dev	+1.1	−2.2	+0.4		−1.8		−0.2
ζ at 3.3 GPa	exp.^2	12.579	7.7219	14.126		111.218		2.275
	SB-CL20	12.433	7.583	14.435		112.25		2.311
	%dev	−1.2	−1.8	+2.2		+0.9		+1.6
	+CCNN	12.328	7.630	14.472		112.28		2.313
	%dev	−2.0	−1.2	+2.4		+1.0		+1.7

---

3.2. Elastic constants of ε-CL20

To compute the components C_ij of the elastic matrix C of ε-CL20 as modeled in Table 2, the simulation box is deformed with a small strain ε_j of 10⁻³. After energy minimization, the induced stresses σ_i are computed and allow the determination of the elastic constants, using the Hooke's law σ_i = C_ijε_j (i and j from 1 to 6, according to Voigt notation). Bulk and shear moduli (B and G, respectively) are computed according to Reuss, Voigt and Hill definitions.³² The results are summarized in Table 3, along with some theoretical and experimental data.

Table 3 Elastic constants of ε-CL20

ε-CL20	MD 300 K, SB-CL20	MD 300 K, SB-CL20 + CCNN	DFT-D, Zhong (2019)^14^a	MD 298 K, COMPASS, Tan (2011)^18	Exp. 296 K, Brillouin scatt., Haycraft (2009)^32
a DFT-D values transformed from P21/a to more common P21/n.
C11	26.39	22.48	30.78	18.89	7.70
C22	27.70	26.54	33.26	18.95	28.29
C33	26.72	28.34	34.26	32.16	28.05
C12	12.90	10.46	0.44	7.33	5.69
C13	6.71	6.13	7.29	4.80	9.21
C23	4.23	2.85	5.95	1.59	−1.22
C44	3.47	2.98	8.61	8.63	12.64
C55	3.93	4.13	12.64	4.24	3.86
C66	9.87	8.40	10.63	3.67	4.73
C15	3.37	3.08	0.51	1.29	1.23
C25	2.87	1.40	3.31	1.78	1.01
C35	−0.89	1.62	−2.04	0.07	3.07
C46	−0.50	−1.27	0.62	−0.66	0.74
B (Reuss)	13.44	11.87	13.68	10.23	7.15
B (Voigt)	14.28	12.92	13.96	10.83	10.16
B (Hill)	13.86	12.40	13.82	10.53	8.66
G (Reuss)	5.33	5.04	11.29	5.57	3.76
G (Voigt)	7.25	6.96	12.01	7.06	7.60
G (Hill)	6.29	6.00	11.65	6.32	5.68
ρ (g cm^−3)	2.055	2.039	1.972	2.067	2.044

---

Both models using SB-CL20 and SB-CL20 + CCNN have quite similar elastic properties. The main difference concerns C₁₁, which is significantly lower with CCNN add-on, but still far from the experimental value. Among the other constants, only C₁₃, C₁₅, C₂₅ and C₅₅ (DFT-D set aside) are found comparable to the other methods. C₃₅ has a different sign without or with CCNN add-on, the sign agreement being with CCNN add-on. C₄₆ is found negative with MD methods, but every value is rather close to zero, which may be the reason for the sign discrepancy. The overall comparison of the C_ij values yields no real consensus. Only C₁₁ seems consistently lower than C₂₂ and C₃₃. The densities of our models are the closest to the experimental one, so a straightforward comparison seems relevant. Nevertheless it should be noted that the elastic constants of β-HMX and α-RDX computed with SB FF are in better agreement with ultrasonic or impulse stimulated light/thermal scattering methods than with Brillouin scattering ones.^33–36 Unfortunately, and to the best of our knowledge, only Brillouin scattering experiments are currently available for ε-CL20. Meanwhile, our bulk moduli and DFT-D ones are in good agreement, but not with COMPASS MD and Haycraft experiment which are too soft. All shear moduli but DFT-D one are in good agreement. Moreover, our bulk moduli are in good agreement with the XRD experimental values from Gump and Peiris⁹ and Pinkerton,¹⁰ which are 13.6 GPa and 13.1 GPa, respectively.

3.3. Equation of state of ε-CL20

Still focusing on ε-CL20, we start again from the modeled structures reported in Table 2 and we perform MD simulations using SB-CL20 FF without or with CCNN dihedral add-on to determine the equation of state (EoS) at 300 K. For 900 ps, we apply a hydrostatic-pressure ramp from 1 atm up to 9 GPa. The resulting pressure–volume curves are shown in Fig. 3-top (left for SB-CL20 and right for SB-CL20 + CCNN). The two light-grey dashed curves are the fits with a 3^rd order Birch–Murnaghan (BM) EoS:where B₀ is the bulk modulus at zero pressure, its derivative and V₀ the volume at zero pressure. Both fits perfectly match the MD results.

Fig. 3 Equation of state (EoS) of ε-CL20. Top: MD at 300 K with SB-CL20 (left) or CCNN add-on (right) and 3^rd order BM EoS fits (dashed grey lines). Bottom: comparison (left) of the 3^rd order BM EoS from theoretical and experimental studies (parameters and references in Table 2), and high-pressure extrapolations (right). The solid or dashed lines stand for the pressure and volume ranges as specified in the respective studies, and the dotted lines stand for the high-pressure extrapolations.

In Fig. 3-bottom-left, we compare our results to other theoretical or experimental 3^rd order BM EoS, of which the parameters are reported in Table 4. Every curve is shown in the pressure and volume ranges as specified in every respective study (except for Zhong DFT-D study which ranges up to 75 GPa). Fig. 3-bottom-right displays high-pressure extrapolations (except for Zhong DFT-D study). Both curves using SB-CL20 FF without (black) or with (red) CCNN add-on are really close to each other. Both agree very well with Zhong DFT-D and Gump XRD experiment, with consistent values of B₀ and . Whatever the given pressure or given volume (including within the extrapolated range), our results deviate by no more than 7% from these two curves. Specifically, this implies a good behavior of our model under pressure. Sorescu DFT-D study confirms this good trend. However, the SRT rigid model seems too stiff at low pressure and the bulk derivatives obtained from COMPASS model, Pinkerton's and Millar's XRD experiments seem too low or too high, and inconsistent with the general trend underpinned by the other estimates.

Table 4 3^rd order BM EoS parameters for ε-CL20: bulk modulus B₀ (GPa), bulk modulus derivative and volume V₀ (ų) at zero pressure

	B0 (GPa)		V0^a (Å^3)
a V0 is taken from zero-pressure simulations or experiments in the corresponding studies (therefore not a result of the fitting procedure but a constrained parameter).b B0 and presently determined for a 3^rd order BM EoS, using data published in the supplemental material of the corresponding studies.c B0 and presently determined for a 3^rd order BM EoS, using data up to 45 GPa of the corresponding study.
MD SB-CL20	13.98	11.05	1414.698
MD SB-CL20 + CCNN	12.16	12.29	1427.593
DFT-D (Sorescu 2010)^10^b	11.61	10.85	1434.675
DFT-D (Zhong 2019)^14^c	11.83	9.81	1476.13
MD COMPASS (Tan 2011)^18	12.8	6.5	1407.80
MD Rigid SRT (Sorescu 1999)^17^b	16.34	8.64	1464.75
Exp. XRD (Gump 2008)^9	13.6	11.7	1423.127
Exp. XRD (Pinkerton 1999)^10^b	13.1	6.9	1430.322
Exp. XRD (Millar 2012)^8	9.5	27	1431.77

---

3.4. Pressure-induced modifications of ε-CL20

3.4.1. Evolution of elastic properties and lattice parameters. We perform the same simulations than in previous section, 900 ps pressure ramps but now up to 90 GPa. Every 20 ps (every increase of 2 GPa), a configuration is extracted to compute the elastic matrix and the bulk and shear moduli (Hill), using the same the process as in Section 3.2 but at the corresponding pressure. The evolution of the elastic properties and the pressure–volume curves extracted from the MD simulations are displayed in Fig. 4. We also display in Fig. 5 the evolution of the lattice parameters during this hydrostatic compression. Note that despite the faster pressure increase in these simulations, the pressure–volume behavior up to 9 GPa perfectly matches the slower compression of Section 3.3.

Fig. 4 Evolution of the elastic properties of ε-CL20 (top-left, bottom-left and bottom-right), and pressure–volume curves extracted from the MD simulations (top-right). The stiffness increases rapidly between 0 and 6 GPa. The arrows and ovals indicate sudden changes of slope (same reference numbers in Fig. 4 and 5).

Fig. 5 Evolution of the lattice parameters of ε-CL20 with pressure. The arrows indicate sudden changes of slope (same reference numbers in Fig. 4 and 5).

The increase of the bulk modulus with pressure (Fig. 4-top-left) confirms that ε-CL20 stiffens under pressure, as stated by Zhong et al. in their DFT-D study.¹⁴ But the elastic properties of our models (Fig. 4-top-left, bottom-left and bottom-right) show some slope changes (indicated by arrows or ovals and referenced by circled numbers) which are not present in the DFT-D study. These changes (see ② and ③ in Fig. 4) happen in the pressure ranges 14–18 GPa and 36–40 GPa with SB-CL20, and at lower pressures with SB-CL20 + CCNN, 12–16 GPa and 30–36 GPa. The pressure–volume curves (Fig. 4-top-right) show no noticeable changes in these respective ranges. However some are visible on the lattice parameter–pressure curves (Fig. 5). From 0 to 6 GPa, every lattice parameter undergoes a marked decrease, with a magnitude found in the same order than in the experiments of Pinkerton,¹⁰ Millar⁸ and Gump and Peiris⁹ (see Fig. S1 in ESI†). In this pressure range, ε-CL20 is found quite soft, but with a rapidly-increasing stiffness. Around 6 GPa (see ① in Fig. 5), there is a first significant variation of the trend, with a practically linear slope and a slope-sign change for “a” and “β”, respectively. This new increase of β is explained by “a” and “c” decreasing differently as from 6 GPa: linearly and still with a varying slope, respectively (β is by definition the angle between the lattice vectors a and c). Due to their limited pressure range, the above experiments cannot confirm this particular trend change. Then the β-slope variation (see ② in Fig. 5) around 16 GPa (red) or 18 GPa (black) is due to “a” and “c” now both varying linearly. The last change (see ③ in Fig. 5) occurring in the 30–40 GPa range is mainly due to another slope variation of “a”.

The molecule of CL-20 forms a quite rigid cage with flexible nitro groups. When the molecular crystal is compressed, the first effects are likely to occur on the orientation of these nitro groups and we address this structural origin in the next section.

3.4.2. Orientation of the nitro groups of the molecule (wag angles). The six N–NO₂ bonds of the CL-20 molecule are quite flexible, hence the different conformers. To study the pressure-induced modifications of the orientation of the nitro groups, we use the wag angle δ, as defined by Munday et al. to study high-pressure RDX phases.³⁶ The wag angle δ is the angle formed by an N–N bond and the corresponding C–N–C triangle, as represented in Fig. 6-left. With six wag angles per CL-20 molecule and four molecules per unit cell (see Fig. 1), a sign convention is required. For the 1^st molecule we start from the first carbon atom listed in the .cif file (ref. 117779 on CCDC). We follow the atom order and assign colors to each N–NO₂ groups: red, green, blue, cyan, pink and yellow (resp. the 1^st, 2^nd, 3^rd, 4^th, 5^th and 6^th groups), as shown in Fig. 6-left. From this representation, an N–NO₂ group above its C–N–C triangle is set positive, and negative otherwise. Then the symmetry-generated 2^nd, 3^rd and 4^th molecules (same atom order as for the 1^st molecule, so same color order of the N–NO₂ groups) are symmetry-corrected to get the same sign convention for the same N–NO₂ group.

Fig. 6 Left: wag angle δ and sign convention. Right: normalized distributions of the wag angles of the six nitro groups (averaged over the four molecular symmetries) in ε-CL20 modeled at 300 K and 1 atm (Δδ = 1°).

As a result, we display in Fig. 6-right the normalized distributions of the wag angles of the six nitro groups (averaged over the four molecular symmetries) in ε-CL20 at 300 K and 1 atm, using SB-CL20 or SB-CL20 + CCNN (see structures in Table 2). The results for the other polymorphs are presented in ESI (Fig. S2–S5).† For ε-CL20, the distributions globally peak at expected angles. However, the mean values reported in Table 5 indicate that the model with CCNN dihedral add-on has the best agreement with the experimental molecular conformation. Consequently, this is the preferred model to investigate the origin of the pressure-induced modifications of ε-CL20. The results presented hereafter are for SB-CL20 + CCNN model of ε-CL20, and we have checked that they are qualitatively the same for SB-CL20 model, but shifted to slightly higher pressures in agreement with Section 3.4.1.

Table 5 Mean wag angles in ε-CL20, from Fig. 6-right

Nitro group	1 – red	2 – green	3 – blue	4 – cyan	5 – pink	6 – yellow
δ from exp. .cif file	−39.6	−1.4	−30.3	−32.7	+23.7	−40.0
δ – MD SB-CL20	−36.1	−4.9	−38.1	−36.3	+18.2	−33.8
δ – MD SB-CL20 + CCNN	−35.9	−2.1	−30.8	−28.3	+25.7	−36.9

---

The wag angle distributions at selected pressures are displayed in Fig. 7, where the four molecular symmetries are systematically shown (and not averaged like in Fig. 6-right). For convenience, a given group has the same color for the four symmetries. First, one can notice that the green and pink distributions shift simultaneously when the pressure increases up to 12 GPa. This corresponds to a concerted rotation of the 2^nd and 5^th nitro groups, as schematized on the molecule at the bottom-left corner. When the pressure comes around 16 GPa, these two distributions begin to split, each of both showing two preferred angles. Since all green and all pink distributions remain superimposed, this degeneracy concerns the four molecular symmetries evenly. Focusing on molecular symmetry #1 (the animation from 0 to 90 GPa is provided in ESI†), we present in Fig. 8 the nature of this degeneracy at 30 GPa, observing the crystal along vector [2a, b, 0] (values in Fig. 5). From this perspective, there are actually two distinct conformations: conformation A (positive “green” peak and negative “pink” peak) and conformation B (slightly negative “green” peak and positive “pink” peak). Moreover, a clear AB pattern appears and, when observed along the adequate vector, we have found the same alternation for the other molecular symmetries.

Fig. 7 Evolution of the wag angles of ε-CL20 with pressure (SB-CL20 + CCNN). For convenience, a given group has the same color for the four molecular symmetries. All same-color curves are superimposed, so the molecular symmetry #1, #2, #3 and #4 undergo the same changes. From 0 to 12 GPa, the 2^nd (green) and 5^th (pink) nitro groups smoothly rotate together (toward the right, as represented on the molecule at the bottom left corner). At 16 GPa, the distributions of these same groups start to split, both showing two preferential wag angles. Degeneracy occurs once more at 40 GPa, involving also the 4^th (cyan) and 6^th (yellow) nitro groups, yielding a partial loss of the ε-conformer identity.

Fig. 8 Degeneracy of the molecular symmetry #1 in ε-CL20 at 30 GPa (SB-CL20 + CCNN), observed along vector [2a, b, 0]. The dotted lines indicate the orientation of the respective C–N–C triangles (see Fig. 6-left).

In Fig. 9, the crystal is observed along vector [0, b, 0] and we superimpose the symmetries #2, #3 and #4 (slightly faded for clarity) to the symmetry #1 (the animation from 0 to 90 GPa is provided in ESI†). Though some defects, a supercell can be evidenced, with the parameters: a′ = 2a = 16.164 Å, b′ = b = 10.343 Å, c′ = 2c = 23.764 Å and β = 104.21°. In every subcell there is one hetero-conformation, yielding an AABB layering system where the layers spread along the supercell's longest diagonal and lattice parameter “b”. Note that using an even replication 14 × 8 × 8 for this simulation (instead of odd 15 × 9 × 9), we obtain at this pressure a supercell-framed crystal with no defect (see Fig. S6 in ESI†). Also, the same simulation performed on a single cell yields four non-degenerated molecules (and with an erroneous triclinic symmetry in the case of SB-CL20 + CCNN FF). This suggests that proper ab initio or lattice dynamics studies at these pressures and above should be conducted at least on 2 × 1 × 2 replications.

Fig. 9 Structure of ε-CL20 supercell, induced by ε-conformer degeneracy at 30 GPa (SB-CL20 + CCNN). The molecular symmetries #2, #3 and #4 (slightly faded for clarity) are superimposed to the molecular symmetry #1, evidencing an AABB layering system.

Up to this pressure, the red, blue, cyan and yellow distributions remain quite unchanged, meaning that the associated nitro groups keep the identity of ε-conformer. Nevertheless, part of this identity is lost over 40 GPa since the cyan and yellow distributions degenerate, along with a second split of the green and pink ones associated to a degeneracy of conformation B. Focusing once more on molecular symmetry #1, we present in Fig. 10 the structure of the degeneracy of the 4^th and 6^th nitro groups at 90 GPa. Both are correlated since there are also two distinct conformations: conformation E (negative “yellow” peak and negative “cyan” peak) and conformation F (slightly negative “yellow” peak and positive “cyan” peak). However, no clear pattern appears, this degeneracy seems spatially random and involves less than 20% of the molecules. Consequently, the study above 40 GPa of the structure, behavior and properties of ε-CL20 with ab initio methods, generally involving small-sized samples, should be considered with caution.

Fig. 10 Degeneracy of the molecular symmetry #1 in ε-CL20 at 90 GPa (SB-CL20 + CCNN).

3.4.3. Discussion. Experimentally, Ciezak et al.⁷ have shown pressure-induced variations in Raman and Far-Infrared spectra of ε-CL20, interpreted as an ε → γ phase transition between 4.1 and 6.4 GPa and another unknown transition between 14.8 and 18.7 GPa. However the ε → γ phase transition is not supported by Millar's more recent experiment,⁸ and not by our study either. Anyway, we support the pressure-induced variations, but with a different interpretation: a smooth wag angle variation of the 2^nd and 5^th nitro groups of ε-conformer between 0 and 12 GPa, their degeneracy as from 16 GPa, and their implication to the sudden slope-variations of the lattice parameters and elastic properties. Besides, the pressure values at which these changes happen (around 6 GPa and 16 GPa for the changes referenced as ① and ② in Section 3.4.1, respectively) agree well with the experimental pressures reported above.

The pressure-induced rotation and degeneracy of some flexible nitro groups are correlated to the evolution of some of the elastic and structural properties. The concerned nitro groups are well-identified and may trigger shock-induced explosive events. While it is currently unclear whether this may come from the sharp degeneracy of the flexible groups (2^nd and 5^th groups) or from a neat bond breaking of the rigid others (1^st, 3^rd, 4^th and 6^th groups), Xue et al.¹⁵ give us some hints in their DFTB-based MD study. Looking at their Fig. 5–7, representing ε-polymorph shocked between 8 and 11 km s⁻¹, it appears that NO₂ fission (resp. ring opening) first occurs at (resp. close to) the position of the rigid groups.

3.5. Implication of CCNN dihedral add-on to γ ↔ ζ phase transitions

To study the γ ↔ ζ phase transition we start from the structures obtained in Table 2. On one hand, the γ-phase undergoes a hydrostatic pressure ramp from 1 atm up to 3.3 GPa over 660 ps (+0.5 GPa/100 ps). We also investigate negative pressure and we apply a descending pressure ramp from 1 atm at the same rate (−0.5 GPa/100 ps). On the other hand, the ζ-phase undergoes a descending pressure ramp from 3.3 GPa to negative pressure at the same rate (−0.5 GPa/100 ps). The pressure–volume and the lattice parameter–pressure curves are displayed in Fig. 11, using SB-CL20 without or with CCNN dihedral add-on.

Fig. 11 Pressure–volume (top) and lattice parameter–pressure (bottom – a, b, c and β) curves during the γ ↔ ζ phase transition simulations. Left: SB-CL20 FF; the ζ → γ transition is observed at −0.5 GPa, whereas the γ → ζ transition fails around 1.4 GPa (not the lattice parameters of ζ). Right: SB-CL20 + CCNN; the complete γ → ζ transition is observed around 1.1 GPa, whereas the ζ → γ transition fails (both ζ and γ samples disintegrate when the pressure becomes lower than −0.7 GPa – see the animation in ESI†).

Using SB-CL20 FF, we are not able to observe the γ → ζ transition (Fig. 11-left, black curves). A first step occurs around 1.4 GPa but the transformation is not complete and the lattice parameters do not match with ζ-phase. In particular, it would have been expected to see “a” decreasing and “b” increasing both more, and the “c” jump is not in the expected direction. The configuration at 1.4 GPa has been extracted and resumed at this constant pressure for additional 1.72 ns without significant improvement. Nevertheless, a complete ζ → γ transition occurs (Fig. 11-left, red curves) but at −0.5 GPa, which is a negative pressure.

Conversely, using SB-CL20 + CCNN FF, which improves the “b” parameter of γ-polymorph (Table 2) and yields the correct sign for the mean wag angle of the 2^nd nitro group (see Fig. S2 in ESI†), we are able to observe a complete γ → ζ transition (Fig. 11-right, black curves) around 1.1 GPa (experimentally 0.7 GPa). But we are not able to observe the ζ → γ transition (Fig. 11-right, red curves) since our ζ sample (and also γ) disintegrates when the pressure becomes lower than −0.7 GPa and the volume larger than 1680 ų per cell (see the animation, the density and pressure curves in ESI†). Longer simulations at selected pressures around −0.7 GPa confirm this observation: either the pressure is slightly above −0.7 GPa and there is no transformation, or it is slightly below −0.7 GPa and the density goes abruptly down to zero and the whole crystal disintegrates.

Experimentally, the γ ↔ ζ phase transition occurs at 0.7 GPa. However, Russell et al.³⁷ have observed that the direct γ → ζ transition preserves the monocrystal whereas the reverse ζ → γ generates crackings, resulting in polycrystals. This behavior may suggest a better agreement with the one of SB-CL20 + CCNN FF, for which the direct transition is straightforward and the reverse one seems more difficult to achieve. This also suggests the possibility of different thermodynamic paths for the direct and reverse transformations, including multiple steps. Anyhow, our results emphasize the complexity of the transitions involving such a large molecule, even in the case of the most kinetically favorable one.

4. Conclusion

We have successfully transferred the quantum-chemistry based force field (FF) developed for HMX by Smith and Bharadwaj (SB)²⁴ to a nitramine of different stoichiometry: hexanitrohexaazaisowurtzitane (CL-20 or HNIW). To assess this flexible-molecule FF, denoted SB-CL20, or SB-CL20 + CCNN for the variant with CCNN dihedral add-on, we have performed Molecular Dynamics simulations at 300 K and under pressures up to 90 GPa. The structures at ambient conditions of the various CL-20 polymorphs are consistent with experimental data. Applying pressure to ε-polymorph, we have determined a relevant equation of state when compared to the trend underpinned by most of the published results. We have also confirmed the pressure-induced increase of the stiffness of ε-CL20.

We have shown that some elastic and structural changes induced at high pressure are related to the flexibility of some of the nitro groups of the molecule, leading to concerted rotations of these groups or ε-conformation degeneracy with a spatial periodicity. At higher pressure, a second degeneracy occurs, spatially random and yielding a partial loss of ε-conformer identity. The implication of the flexible/rigid groups to trigger explosive events is to be studied, along with the effects of co-crystallization on their stabilization/flexibilization.

The compression of ε-CL20 using SB-CL20 FF without or with CCNN dihedral add-on results in very similar qualitative behaviors. This is not the case concerning the reversible γ ↔ ζ phase transition. At best, SB-CL20 + CCNN FF reproduces the direct γ → ζ transition and SB-CL20 FF the reverse ζ → γ one. The simulations suggest the possibility of different multiple-step, thermodynamic paths for the direct and reverse transitions. A comparison with the experimental behavior may suggest a better agreement with SB-CL20 FF + CCNN, but the full mechanism remains to be elucidated.

At last, the transferability of SB-CL20 FFs over the five CL-20 polymorphs and of SB FF over HMX and RDX allows investigating the properties of co-crystals designed for modified power and sensitivity. For instance in the case of the 2 : 1 co-crystal of CL-20 : HMX,³⁸ the present flexible-molecule model should enable understanding of how HMX structurally and dynamically alters CL-20 to make it less sensitive than but as powerful as the single crystal.

Conflicts of interest

There are no conflicts of interest to declare.

Acknowledgements

This work was supported by the Air Force Office of Scientific Research (AFOSR) under grant number FA9550-18-1-0236, and used the HPC resources from the Advanced Cyberinfrastructure for Education and Research (ACER) at the University of Illinois at Chicago (UIC).

References

1. A. T. Nielsen, A. P. Chafin, S. L. Christian, D. W. Moore, M. P. Nadler, R. A. Nissan, D. J. Vanderah, R. D. Gilardi, C. F. George and J. L. Flippen-Anderson, Synthesis of polyazapolycyclic caged polynitramines, Tetrahedron, 1998, 54, 11793–11812.
2. D. I. A. Millar, H. E. Maynard-Casely, A. K. Kleppe, W. G. Marshall, C. R. Pulham and A. S. Cumming, Putting the squeeze on energetic materials—structural characterisation of a high-pressure phase of CL-20, CrystEngComm, 2010, 12, 2524–2527.
3. T. P. Russell, P. J. Miller, G. J. Piermarini and S. Block, Pressure/temperature phase diagram of hexanitrohexaazaisowurtzitane, J. Phys. Chem., 1993, 97, 1993–1997.
4. J. C. Gump, C. A. Stoltz and S. M. Peiris, Phase stability of epsilon and gamma hniw (cl-20) at high-pressure and temperature, AIP Conf. Proc., 2007, 955, 127–132.
5. M. F. Foltz, C. L. Coon, F. Garcia and A. L. Nichols, The thermal stability of the polymorphs of hexanitrohexaazaisowurtzitane, Part I, Propellants, Explos., Pyrotech., 1994, 19, 19–25.
6. Y. Kholod, S. Okovytyy, G. Kuramshina, M. Qasim, L. Gorb and J. Leszczynski, An analysis of stable forms of CL-20: A DFT study of conformational transitions, infrared and Raman spectra, J. Mol. Struct., 2007, 843, 14–25.
7. J. A. Ciezak, T. A. Jenkins and Z. Liu, Evidence for a High-Pressure Phase Transition of ε-2,4,6,8,10,12-Hexanitrohexaazaisowurtzitane (CL-20) Using Vibrational Spectroscopy, Propellants, Explos., Pyrotech., 2007, 32, 472–477.
8. D. I. A. Millar, in Energetic Materials at Extreme Conditions, ed. D. I. A. Millar, Springer Berlin Heidelberg, Berlin, Heidelberg, 2012, pp. 95–124.
9. J. C. Gump and S. M. Peiris, Phase transitions and isothermal equations of state of epsilon hexanitrohexaazaisowurtzitane (CL-20), J. Appl. Phys., 2008, 104, 083509.
10. D. C. Sorescu and B. M. Rice, Theoretical Predictions of Energetic Molecular Crystals at Ambient and Hydrostatic Compression Conditions Using Dispersion Corrections to Conventional Density Functionals (DFT-D), J. Phys. Chem. C, 2010, 114, 6734–6748 (Pinkerton did the experiments mentioned as a 1999 private communication, of which the data are however published in their ESI).
11. X.-J. Xu, W.-H. Zhu and H.-M. Xiao, DFT Studies on the Four Polymorphs of Crystalline CL-20 and the Influences of Hydrostatic Pressure on ε-CL-20 Crystal, J. Phys. Chem. B, 2007, 111, 2090–2097.
12. J. Cabalo and R. Sausa, Theoretical and Experimental Study of the C–H Stretching Overtones of 2,4,6,8,10,12-hexanitro-2,4,6,8,10,12 hexaazaisowurtzitane (CL20), J. Phys. Chem. A, 2013, 117, 9039–9046.
13. E. F. C. Byrd, G. E. Scuseria and C. F. Chabalowski, An ab Initio Study of Solid Nitromethane, HMX, RDX, and CL20: Successes and Failures of DFT, J. Phys. Chem. B, 2004, 108, 13100–13106.
14. M. Zhong, H. Qin, Q.-J. Liu, F.-S. Liu and B. Tang, Structures, Elasticity, and Sensitivity Characteristics of ε-CL-20 under High Pressure from First-Principles Calculations, Phys. Status Solidi B, 2019, 256(5), 1800440.
15. X. Xue, Y. Wen and C. Zhang, Early Decay Mechanism of Shocked ε-CL-20: A Molecular Dynamics Simulation Study, J. Phys. Chem. C, 2016, 120, 21169–21177.
16. O. Isayev, L. Gorb, M. Qasim and J. Leszczynski, Ab Initio Molecular Dynamics Study on the Initial Chemical Events in Nitramines: Thermal Decomposition of CL-20, J. Phys. Chem. B, 2008, 112, 11005–11013.
17. D. C. Sorescu, B. M. Rice and D. L. Thompson, Theoretical Studies of the Hydrostatic Compression of RDX, HMX, HNIW, and PETN Crystals, J. Phys. Chem. B, 1999, 103, 6783–6790.
18. J.-J. Tan, G.-F. Ji, X.-R. Chen and Z. Li, Structure, equation of state and elasticity of crystalline HNIW by molecular dynamics simulations, Phys. B, 2011, 406, 2925–2930.
19. F. Wang, L. Chen, D. Geng, J. Lu and J. Wu, Effect of density on the thermal decomposition mechanism of ε-CL-20: a ReaxFF reactive molecular dynamics simulation study, Phys. Chem. Chem. Phys., 2018, 20, 22600–22609.
20. C. Ren, X. Li and L. Guo, Reaction Mechanisms in the Thermal Decomposition of CL-20 Revealed by ReaxFF Molecular Dynamics Simulations, Acta Phys.-Chim. Sin., 2018, 34, 1151–1162.
21. D. C. Sorescu, B. M. Rice and D. L. Thompson, A Transferable Intermolecular Potential for Nitramine Crystals, J. Phys. Chem. A, 1998, 102, 8386–8392.
22. M. Bergh and C. Caleman, A Validation Study of the General Amber Force Field Applied to Energetic Molecular Crystals, J. Energ. Mater., 2016, 34, 62–75.
23. P. M. Agrawal, B. M. Rice, L. Zheng and D. L. Thompson, Molecular Dynamics Simulations of Hexahydro-1,3,5-trinitro-1,3,5-s-triazine (RDX) Using a Combined Sorescu−Rice−Thompson AMBER Force Field, J. Phys. Chem. B, 2006, 110, 26185–26188.
24. G. D. Smith and R. K. Bharadwaj, Quantum Chemistry Based Force Field for Simulations of HMX, J. Phys. Chem. B, 1999, 103, 3570–3575.
25. D. Bedrov, C. Ayyagari, G. D. Smith, T. D. Sewell, R. Menikoff and J. M. Zaug, Molecular dynamics simulations of HMX crystal polymorphs using a flexible molecule force field, J. Comput.-Aided Mater. Des., 2001, 8, 77–85.
26. M. J. Cawkwell, T. D. Sewell, L. Zheng and D. L. Thompson, Shock-induced shear bands in an energetic molecular crystal: Application of shock-front absorbing boundary conditions to molecular dynamics simulations, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 78, 014107.
27. D. Bedrov, J. B. Hooper, G. D. Smith and T. D. Sewell, Shock-induced transformations in crystalline RDX: A uniaxial constant-stress Hugoniostat molecular dynamics simulation study, J. Chem. Phys., 2009, 131, 034712.
28. C. J. Fennell and J. D. Gezelter, Is the Ewald summation still necessary? Pairwise alternatives to the accepted standard for long-range electrostatics, J. Chem. Phys., 2006, 124, 234104.
29. W. L. Jorgensen, D. S. Maxwell and J. Tirado-Rives, Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids, J. Am. Chem. Soc., 1996, 118, 11225–11236.
30. J. D. Gale and A. L. Rohl, The General Utility Lattice Program (GULP), Mol. Simul., 2003, 29, 291–341.
31. S. Plimpton, Fast Parallel Algorithms for Short-Range Molecular Dynamics, J. Comput. Phys., 1995, 117, 1–19.
32. J. J. Haycraft, The elastic constants and related properties of the epsilon polymorph of the energetic material CL-20 determined by Brillouin scattering, J. Chem. Phys., 2009, 131, 214501.
33. L. L. Stevens and C. J. Eckhardt, The elastic constants and related properties of β-HMX determined by Brillouin scattering, J. Chem. Phys., 2005, 122, 174701.
34. B. Sun, J. M. Winey, Y. M. Gupta and D. E. Hooks, Determination of second-order elastic constants of cyclotetramethylene tetranitramine (β-HMX) using impulsive stimulated thermal scattering, J. Appl. Phys., 2009, 106, 053505.
35. B. Sun, J. M. Winey, N. Hemmi, Z. A. Dreger, K. A. Zimmerman, Y. M. Gupta, D. H. Torchinsky and K. A. Nelson, Second-order elastic constants of pentaerythritol tetranitrate and cyclotrimethylene trinitramine using impulsive stimulated thermal scattering, J. Appl. Phys., 2008, 104, 073517.
36. L. B. Munday, P. W. Chung, B. M. Rice and S. D. Solares, Simulations of High-Pressure Phases in RDX, J. Phys. Chem. B, 2011, 115, 4378–4386.
37. T. P. Russell, P. J. Miller, G. J. Piermarini and S. Block, High-pressure phase transition in .gamma.-hexanitrohexaazaisowurtzitane, J. Phys. Chem., 1992, 96, 5509–5512.
38. O. Bolton, L. R. Simke, P. F. Pagoria and A. J. Matzger, High Power Explosive with Good Sensitivity: A 2:1 Cocrystal of CL-20:HMX, Cryst. Growth Des., 2012, 12, 4311–4314.

---

Footnote

† Electronic supplementary information (ESI) available: FF parameter and configuration files to get started with ε-CL20 simulations using LAMMPS and our SB-CL20 FFs. Model/experiment comparison of the pressure-induced variations of ε-CL20 lattice parameters. Wag angle distributions of α, β, γ and ζ polymorphs modeled with our SB-CL20 FFs at ambient conditions (or 3.3 GPa for ζ), and comparison with experimental data. Replication-induced defects in supercell-framed ε-CL20 at 30 GPa, using SB-CL20 + CCNN FF. Two animations of the compression of ε-CL20 modeled with SB-CL20 + CCNN FF from 0 to 90 GPa. Density and pressure during the decompression of the γ sample modeled with SB-CL20 + CCNN FF, up to its disintegration. One animation of the disintegration at negative pressure of the ζ sample modeled with SB-CL20 + CCNN FF. See DOI: 10.1039/c9ra07645j

---