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

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. This journal is © The Royal Society of Chemistry.

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 C6H6O12N12, it has five polymorphs: α, β, γ, ε[1] and high-pressure ζ.[2]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.[4] α-phase is stabilized by its ability to trap small molecules like H2O,[1] CO or CO2.[3] 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.[7] 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[8] and Gump and Peiris.[9] The equations of state determined by XRD experiments by Millar,[8] Gump and Peiris[9] and Pinkerton[10] 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,[15] ab initio Molecular Dynamics (AIMD),[16] and Molecular Dynamics (MD) using various force fields (FF) such as Sorescu–Rice–Thompson SRT (rigid molecules),[17] COMPASS (flexible molecules, but not free-to-use)[18] 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,[21] its intermolecular parameters have been merged to the generic intramolecular topology of the general Amber FF (GAFF), yielding Amber-SRT FF.[22] However the crystal densities of such modeled nitramines are systematically underestimated. Agrawal et al.[23] 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.

Computational methods

SB-CL20 force field

Smith & Bharadwaj[24] have developed a quantum-chemistry based force field (SB FF) to model HMX (C4H8O8N8, 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,[28] with a damping parameter of 0.25 Å−1.

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 (C3H6O6N6), but not to CL-20 (C6H6O12N12, 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 qCL20C = −0.2160e. The charges used for the other atoms are the ones prescribed by Smith & Bharadwaj in their original paper.[24] 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−1 and r0 = 1.53 Å, corresponding to harmonic bond energy U = k(r − r0)2/2. The angular energy term has the harmonic form U = k(θ − θ0)2/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−1 and θ0 = 110.70°. These bond and angle parameters are simply taken from OPLS-AA FF.[29]

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).[30] 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°].

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

Molecular dynamics

To assess SB-CL20 FF, we perform MD simulations with LAMMPS code (large-scale atomic/molecular massively parallel simulator).[31] 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.†

Results and discussions

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.[25] 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.

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

Elastic constants of ε-CL20

To compute the components C of the elastic matrix C of ε-CL20 as modeled in Table 2, the simulation box is deformed with a small strain ε of 10−3. After energy minimization, the induced stresses σ are computed and allow the determination of the elastic constants, using the Hooke's law σ = Cε (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.[32] The results are summarized in Table 3, along with some theoretical and experimental data.

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]
C 11	26.39	22.48	30.78	18.89	7.70
C 22	27.70	26.54	33.26	18.95	28.29
C 33	26.72	28.34	34.26	32.16	28.05
C 12	12.90	10.46	0.44	7.33	5.69
C 13	6.71	6.13	7.29	4.80	9.21
C 23	4.23	2.85	5.95	1.59	−1.22
C 44	3.47	2.98	8.61	8.63	12.64
C 55	3.93	4.13	12.64	4.24	3.86
C 66	9.87	8.40	10.63	3.67	4.73
C 15	3.37	3.08	0.51	1.29	1.23
C 25	2.87	1.40	3.31	1.78	1.01
C 35	−0.89	1.62	−2.04	0.07	3.07
C 46	−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

DFT-D values transformed from P21/a to more common P21/n.

Both models using SB-CL20 and SB-CL20 + CCNN have quite similar elastic properties. The main difference concerns C11, which is significantly lower with CCNN add-on, but still far from the experimental value. Among the other constants, only C13, C15, C25 and C55 (DFT-D set aside) are found comparable to the other methods. C35 has a different sign without or with CCNN add-on, the sign agreement being with CCNN add-on. C46 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 values yields no real consensus. Only C11 seems consistently lower than C22 and C33. 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[9] and Pinkerton,[10] which are 13.6 GPa and 13.1 GPa, respectively.

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 3rd order Birch–Murnaghan (BM) EoS:where B0 is the bulk modulus at zero pressure, its derivative and V0 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 3rd order BM EoS fits (dashed grey lines). Bottom: comparison (left) of the 3rd 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 3rd 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 B0 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.

3rd order BM EoS parameters for ε-CL20: bulk modulus B0 (GPa), bulk modulus derivative and volume V0 (Å3) at zero pressure

	B 0 (GPa)		V 0 a (Å3)
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

V 0 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 0 and presently determined for a 3rd order BM EoS, using data published in the supplemental material of the corresponding studies.

B 0 and presently determined for a 3rd order BM EoS, using data up to 45 GPa of the corresponding study.

Pressure-induced modifications of ε-CL20

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.[14] 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,[10] Millar[8] and Gump and Peiris[9] (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.

Orientation of the nitro groups of the molecule (wag angles)

The six N–NO2 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.[36] 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 1st 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–NO2 groups: red, green, blue, cyan, pink and yellow (resp. the 1st, 2nd, 3rd, 4th, 5th and 6th groups), as shown in Fig. 6-left. From this representation, an N–NO2 group above its C–N–C triangle is set positive, and negative otherwise. Then the symmetry-generated 2nd, 3rd and 4th molecules (same atom order as for the 1st molecule, so same color order of the N–NO2 groups) are symmetry-corrected to get the same sign convention for the same N–NO2 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.

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 2nd and 5th 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 2nd (green) and 5th (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 4th (cyan) and 6th (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 4th and 6th 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).

Discussion

Experimentally, Ciezak et al.[7] 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,[8] 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 2nd and 5th 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 (2nd and 5th groups) or from a neat bond breaking of the rigid others (1st, 3rd, 4th and 6th groups), Xue et al.[15] 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−1, it appears that NO2 fission (resp. ring opening) first occurs at (resp. close to) the position of the rigid groups.

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 2nd 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 Å3 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.[37] 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.

Conclusion

We have successfully transferred the quantum-chemistry based force field (FF) developed for HMX by Smith and Bharadwaj (SB)[24] 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,[38] 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.