In order to investigate the precise structure and properties of energetic molecular crystals, classic molecular FFs were refitted for EMs (Figure 2 and Table 1). In the 1990s, Sorescu et al. [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 [26]. Smith et al. [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 [28,29,30,31].

In the 2000s, more classic FFs were applied. For RDX crystal, Agrawal et al. [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. [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. [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. [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 [37].

In the 2010s, the symmetry-adapted perturbation theory (SAPT)-based potential was also applied to describe thermal properties for FOX-7 [38], and thermal expansion and pressure response were investigated thereby. Neyertz et al. [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. [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 [41,42].

In 2021, Wang et al. developed a tailor-made OPLS all-atom FF (OPLS-AA) [43] with polymorphism of CL-20 being predicted successfully [44].

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 [20]. The energy expression of SAPT potential is similar to the SRT pair potential, and it was also applied on FOX-7 [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 [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 [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 [34]. In comparison, NETMFF is more complex, with more valence terms included, as listed in Table 1 [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 E_el and a Lennard–Jones LJ 12-6 term [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 [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 sp² ring and NO₂ structures maintain planarity by using harmonic function [39]. Moreover, some general FFs were also refitted for energetic crystals, such as Dreiding for TATB [41] and OPLS-AA for CL-20 [43].

To date, the effects of different molecular structures—especially functional groups—were considered to make precise predictions. Initially, SRT [20] and SAPT FFs [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 [32] accurately predicted the chair and inverted chair conformation, bond lengths, and bond angles of the RDX molecule; however, the parameters for the N-NO₂ improper dihedral angles are not available in AMBER, and the parameters for O-NO₂ were used in SRT-AMBER, and there are some inaccuracies in the calculated orientations of the NO₂ groups; thus, modifications in the torsional parameters are needed. SB potential [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 [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-NO₂ rotational barrier. Moreover, in NETMFF for RDX [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 NO₂ 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 NO₂ group are well described, and the lattice parameters and thermal expansion are well predicted. In GRBF for TATB [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 [39], as well as in Bedrov’s [36] and Dreiding potentials for TATB [41]. The differences in the functional forms, including descriptions of functional groups for the FFs, bring differences in prediction precision and areas of application.