Classic Forcefield Refitted for Energetic Materials: Comparison
Please note this is a comparison between Version 1 by Wen Qian and Version 2 by Beatrix Zheng.

Energetic materials (EM) act as the key chemical activity component in weapon systems, and energetic crystals are the main and functional part therein. Until now, the most widely used energetic crystals have been molecular crystals, including TNT, RDX, HMX, and TATB. The thermal properties, mechanical properties, and chemical reactivity of energetic molecular crystals are of great importance for a whole EM. The classic forcefields (FF) are still useful for typical energetic crystals—i.e., even though they generally have simple potential functions, the accuracy can be ensured by refitting for specific compound. 

  • molecular forcefield
  • energetic materials
  • molecular dynamics simulation

1. Development of Refitted FFs

In order to investigate the precise structure and properties of energetic molecular crystals, classic molecular classic forcefields (FF)s were refitted for Energetic materials (EMs)Ms (Figure 12 and Table 1). In the 1990s, Sorescu et al. [1][2][3][4][5][6][20,21,22,23,24,25] developed a classic FF (SRT) with pair potential refitted for predicting lattice parameters, density, bulk modulus, and shear modulus of aliphatic energetic compounds—including nitramine heterocycles RDX, HMX, and CL-20, and linear conjugated energetic compound FOX-7. It was also applied on linear nitrate PETN which has overloaded substituents with lower accuracy [7][26]. Smith et al. [8][27] developed another kind of multi-particle potential, the Smith–Bharadwaj (SB) potential, which is widely used in the calculations of shock compression, shear bands, elastic constants, and modulus for HMX, CL-20, and RDX [9][10][11][12][28,29,30,31].
Figure 12.
Development of classic FFs refitted for EMs.
In the 2000s, more classic FFs were applied. For RDX crystal, Agrawal et al. [13][14][32,33] combined the advantages of traditional AMBER FF with SRT FF, and applied the refitted FF (SRT-AMBER) which can effectively describe the flexibility of molecules on the calculations of cell parameters, melting point, and mechanics for TNAZ and RDX. Boyd et al. [15][34] also developed a multi-particle potential for RDX; thereby, the vibration spectra, thermodynamics, thermal expansion, and mechanics were calculated. For aromatic nitro compound TATB, a multi-particle FF (GRBF) with Lennard–Jones portion is developed by Gee et al. [16][35] using an ab initio approach for large molecular system TATB crystal, with cell parameters, density, thermal expansion, and pressure–volume isotherm prediction. Bedrov et al. [17][36] developed a potential to calculate the vibration spectra, elastic stiffness coefficients, and isotropic moduli; moreover, thermal conductivity can also be obtained using non-equilibrium molecular dynamics (NEMD) under Bedrov’s potential [18][37].
In the 2010s, the symmetry-adapted perturbation theory (SAPT)-based potential was also applied to describe thermal properties for FOX-7 [19][38], and thermal expansion and pressure response were investigated thereby. Neyertz et al. [20][39] developed a classic potential for aromatic nitro compounds TNT and DNT, which possess conjugated structure, and calculated the cell parameters, density, tensile modulus, bulk modulus, and shear modulus. Song et al. [21][40] developed an all-atom, non-empirical, tailor-made FF (NETMFF) for RDX, and cell parameters, density, and thermal expansion were predicted thereby; the Dreiding FF was also applied on TATB, with the anisotropic thermal expansion simulated successfully [22][23][41,42].
Table 1.
Energy expressions and application scopes for the classic FFs.
FFs Valence Terms van der Waals

Interaction Term
Electrostatic Interaction Term Applications
SRT [1]SRT [20] \ Buckingham 6-exp form Coulomb function RDX HMX CL-20 FOX-7 PETN

(lattice parameter, density, mechanics)
SAPT [19]SAPT [38] \ Buckingham 6-exp form Simplified Coulomb function FOX-7 (thermal properties, pressure responses, isothermo)
SRT-AMBER [13]SRT-AMBER [32] Harmonic bond stretching, harmonic angle bending, cosine torsion term Buckingham 6-exp form Coulomb function RDX (lattice parameter, density, melting point, mechanics)

TNAZ (lattice parameter, density, melting point)
SB [8]SB [27] Harmonic bond term, angle term, dihedral term, anharmonic torsion term Buckingham form Coulomb function RDX HMX CL-20

(shock compression, shear bands, elastic constants and modulus)
Boyd’s [15]Boyd’s [34] Bond stretching described by Morse function, angle bending described by harmonic function Buckingham LJ 6-12 form Coulomb function RDX (lattice parameter, density, thermodynamics, vibration spectra, thermal expansion, mechanics)
NETMFF [21]NETMFF [40] Bond term, angle term, dihedral (torsion angle) term, out-of-plane bending angle term, cross-coupling terms of bond–bond, bond–angle couplings Damped Buckingham form Coulomb function RDX (lattice parameter, density, thermal expansion)
GRBF [16]GRBF [35] Harmonic bond stretch term, bond-angle bend term, dihedral angle torsion term Lennard–Jones 12-6 form Coulomb function TATB (lattice parameters, density, thermal expansion, isotherm)
Bedrov’s [17]Bedrov’s [36] Harmonic functions of covalent bonds, three-center bends, and improper dihedrals Buckingham 6-exp form Coulomb function TATB (lattice parameter, thermal expansion, mechanics, vibration spectra, thermal conductivity)
Neyertz’s [20]Neyertz’s [39] Angle-bending deformations described by harmonic function, torsional motions around the dihedral angles τ, sp2 ring and NO2 structures kept planar described by harmonic function Lenard-Jones 12-6 form Coulomb function TNT DNT (lattice parameter, density, tensile, bulk and shear modulus)
Dreiding [22]Dreiding [41] Bond stretching interaction term, angle bending interaction term, dihedral angle interaction term, inversion interaction term Lennard–Jones 12-6 form Coulomb function TATB (geometries, crystal packing, thermal expansion)
OPLS-AA [24]OPLS-AA [43] Bond term, angle term, dihedral term Lennard–Jones 12-6 form Coulomb function CL-20 (lattice parameter, density, polymorph prediction)
In 2021, Wang et al. developed a tailor-made OPLS all-atom FF (OPLS-AA) [24][43] with polymorphism of CL-20 being predicted successfully [25][44].

2. Functional Forms of the Refitted FFs

The energy expressions of the classic FFs have something in common, i.e., most of FFs have valence terms, van der Waals interaction terms, and electrostatic interaction terms, and they differ only in functional forms (Table 1). For aliphatic energetic compounds including heterocycles RDX, HMX, CL-20, and linear PETN, FOX-7, SRT potential was applied. The energy expression of SRT only has a pairwise Buckingham 6-exp form with repulsion and dispersion interactions, and a Coulombic potential of electrostatic interactions [1][20]. The energy expression of SAPT potential is similar to the SRT pair potential, and it was also applied on FOX-7 [19][38].
For the multi-particle potentials of nitramine heterocycles, the energy terms of SB potential cover harmonic bonds, angles, dihedrals terms, anharmonic torsions, and non-bonded Buckingham and Coulomb interactions [8][27]. SRT-AMBER combines the forms of SRT and AMBER FFs, with harmonic bond stretching and angle bending terms, cosine torsion term, and nonbond term replaced by Buckingham (exp-6) potential and Coulombic potential in SRT [13][32]. Boyd’s multi-particle potential contains bond stretching, angle bending, vdW term, and electrostatic term; described by Morse function, harmonic function, Buckingham potential LJ 6-12, and Coulomb function, respectively [15][34]. In comparison, NETMFF is more complex, with more valence terms included, as listed in Table 1 [21][40].
For the aromatic nitro compounds, GRBF and Bedrov’s potential were applied for TATB, and Neyertz’s FF was applied for TNT and DNT. GRBF FF includes valence terms containing harmonic bond stretch, bond-angle bend, dihedral angle torsion, and intermolecular pair potentials, including a Coulombic electrostatic term Eel and a Lennard–Jones LJ 12-6 term [16][35]. Bedrov’s potential includes the harmonic functions of covalent bonds, three-center bends and improper dihedrals, and nonbonded interactions of Buckingham (exp-6) potential and Coulomb interactions [17][36]. In Neyertz’s potential, angle-bending deformations are described by harmonic function, and torsional motions around the dihedral angles τ are represented by a polynomial in cos τ, while sp2 ring and NO2 structures maintain planarity by using harmonic function [20][39]. Moreover, some general FFs were also refitted for energetic crystals, such as Dreiding for TATB [22][41] and OPLS-AA for CL-20 [24][43].
To date, the effects of different molecular structures—especially functional groups—were considered to make precise predictions. Initially, SRT [1][20] and SAPT FFs [19][38] are simple pair potentials which only describe the intermolecular interactions, with assumption of rigid molecules, thus they cannot describe the effect of different molecular structures, resulting in lower accuracy. Combing the rigid-molecule SRT potential with intramolecular interactions from AMBER, the SRT-AMBER potential [13][32] accurately predicted the chair and inverted chair conformation, bond lengths, and bond angles of the RDX molecule; however, the parameters for the N-NO2 improper dihedral angles are not available in AMBER, and the parameters for O-NO2 were used in SRT-AMBER, and there are some inaccuracies in the calculated orientations of the NO2 groups; thus, modifications in the torsional parameters are needed. SB potential [8][27] with anharmonic torsions terms that display extrema at the torsion angles that correspond to stationary points on the conformational energy surface can effectively predict lattice parameters and elastic tensors for RDX and HMX. Furthermore, Boyd’s potential [15][34] managed to stabilize the RDX crystal lattice with flexible molecules in the correct conformation, with Morse bond stretching, harmonic angle bending, cosine torsions terms, and in which the torsion parameter C-N-N-O values were modified to adjust the N-NO2 rotational barrier. Moreover, in NETMFF for RDX [21][40], many more parameters regarding functional groups were considered, including parameters c_4 for the carbon atom, n_3r for the nitrogen atom of the triazine ring, n_3o for the nitrogen atom of the NO2 group, o_1 for the oxygen atom, and h_1 for the hydrogen atom; as well as torsion parameters for six RDX conformers. Consequently, the angles and torsions for the NO2 group are well described, and the lattice parameters and thermal expansion are well predicted. In GRBF for TATB [16][35], the amino C-C-N-H and nitro C-C-N-O group rotational barriers were modified, and the nonbonded interactions for amino groups and nitro groups were also considered; similar parameters were considered in Neyertz’s potential for TNT and DNT [20][39], as well as in Bedrov’s [17][36] and Dreiding potentials for TATB [22][41]. The differences in the functional forms, including descriptions of functional groups for the FFs, bring differences in prediction precision and areas of application.

3. Application of Prediction

3.1. Cell Parameters and Density

The refitted FFs can effectively predict cell parameters, and further density. For example, the cell parameter and density of RDX have been separately predicted using SRT, SRT-AMBER, Boyd’s FF, SB, and NETMFF potentials [1][8][13][15][21][20,27,32,34,40]. In comparison, the density prediction by NETMFF is much closer to the experimental results [26][45], while those by SRT and SRT-AMBER are not so accurate. It is attributed to SRT’s lack of valence terms, while NETMFF has the greatest abundance of valence terms. It shows that the valence terms are important in describing the structural properties of cell parameters and density. Furthermore, the cell parameters and density of TATB can be described well by GRBF, Bedrov’s, and Dreiding FFs [15][16][23][34,35,42]; and those of TNT, DNT, and CL-20 can be described well by Neyertz’s FF [20][39], and SB [8][27] and OPLS-AA FFs [25][44], respectively.

3.2. Polymorphism

Classic FFs have been refitted to predict various crystalline phases and their properties. For example, the Neyertz’s FF [20][39] was applied to monoclinic and orthorhombic TNT, with density, tensile modulus, bulk modulus, and shear modulus predicted, in agreement with experimental results. Wang et al. [25][44] also applied the OPLS-AA FF for the polymorph prediction of CL-20, in which the FF parameters were refitted based on DFT calculations with corrected density functional M06-2X, with most polymorphs being reproduced successfully.

3.3. Vibration Spectra

Kroonblawd et al. [18][37] used the modified Bedrov’s potential to investigate the vibration spectra of TATB. The modified FF overcame the gap between the prediction by the original one and experimental observations. The calculated spectrum can effectively assign the vibration modes, including amine antisymmetric stretch, amine symmetric stretch, amine scissoring plus C-NH2 stretch, nitro antisymmetric stretch plus amine scissor/rock, and ring stretch. The Boyd’s potential [15][34] has also been effectively applied to investigating the vibration spectra for RDX, in which special attention has been paid to the vibrational states between 200 and 700 cm−1 described as “doorway modes” for the transfer of energy from lattice phonons to the molecular vibrations involved in bond breaking [27][28][46,47].

3.4. Thermal Property

Among the thermal properties, thermal expansion is an important feature related to performance. The refitted FFs effectively described the linear and volume thermal expansion for series of energetic crystals: the SRT FF has been refitted for FOX-7, and the linear and volume thermal expansion coefficients (CTE) of FOX-7 crystal were determined from the averages of lattice dimensions extracted from trajectories of MD calculations; the results of linear CTE values show high anisotropy because of the layer wave-like stacking style inside the FOX-7 crystal, in agreement with experiment [7][26]. The MD simulations for CTE of FOX-7 were also carried out with the fitted SAPT potential, showing a significant anisotropy along different crystalline directions [19][38]. Compared with experiment [29][48], SRT is better than SAPT, since SRT and SAPT have the same vdW term, while SAPT has a simplified Coulomb term. It demonstrates the importance of electrostatic interaction in the thermal expansion of FOX-7 crystal. Furthermore, SRT, SRT-AMBER, Boyd’s, and NETMFF FFs were refitted for thermal expansion of RDX crystal. Compared with the experiments [29][48], it was found that SRT and Boyd’s FFs describe the anisotropy better, while NETMFF tends to isotropy [1][8][15][21][20,27,34,40]. The GRBF FF [16][35] was used to predict the CTE for TATB, agreeing well with the experiments [30][49]. Dreiding FF was also applied to the TATB crystal, with the planar stacking maintained and the anisotropic thermal expansion reproduced successfully [22][23][41,42]. Thermal conductivity of TATB was calculated using NEMD with Bedrov’s potential [17][36]. Therein, the thermal flux, temperature gradient, as well as resulting thermal conductivity, were obtained from direct MD simulations of the slabs, the results from different directions show high anisotropy [18][37].

3.5. Mechanical Property

The SRT pair potential was refitted for RDX using quantum chemistry calculations at the theoretical level of MP2, and NPT-MD with this refitted potential were performed to predict mechanical properties—including bulk modulus and volume compressibility which match with experiments [1][20]. The SRT-AMBER potential has also been used to investigate the mechanics for crystalline RDX; however, the prediction accuracy of SRT-AMBER has no advantages compared with other FFs, thus the FF parameters need further modification [13][32]. Additionally, the SB potential has widely been used to calculate the mechanical properties for energetic crystals of HMX, CL-20, and RDX, with elastic constants and bulk modulus being obtained [8][27]; the Neyertz’s potential was applied to predict the tensile modulus, bulk modulus, and shear modulus for TNT and DNT [20][39]; and the elastic stiffness coefficient, bulk modulus, and shear modulus of TATB were accurately predicted using Bedrov’s potential [17][18][36,37].

3.6. Shock Responses

For shock responses, the SB potential can effectively describe the shock-induced shear bands in RDX crystal, in which the intermolecular and intramolecular temperatures of a crystalline region and adjacent shear band can be calculated [9][28]. Moreover, it can describe the shock compression in the RDX crystal, including stress along compression direction, temperature evolution under shock compression, compression ratio, and temperature as a function of shock pressure [10][29].
Video Production Service