All-atom, non-empirical, and tailor-made force field for α-RDX from first principles

Abstract

A general strategy is presented to develop an all-atom, non-empirical, and tailor-made force field (NETMFF) for high explosives (HEs). The central part of the strategy is a self-consistent force field (SCFF) optimization technique. The consistence of the force field is ensured by iterating the parameterization procedure. The generation of reliable ab initio reference data for optimizing the NETMFF parameters and the SCFF technique are discussed in detail. Starting with the crystal structure obtained from either experiment or crystal structure prediction (CSP), NETMFF can be developed only by first principles, including conventional DFT, the periodic DFT-D model, and the SAPT(DFT) plus Williams–Stone–Misquitta (WSM) methods. The development strategy of NEMTFF has been applied to α phase hexahydro-1,3,5-trinitro-1,3,5-triazine (α-RDX), an important high explosive. The force field obtained for α-RDX can yield the correct geometries for all RDX conformers found by the DFT calculations in a good way, and also can identify the transition states obtained from the DFT calculations by force field-based vibration analysis. More importantly, it can accurately predict densities of high explosives under the environmental conditions to which they are often subjected, a long-standing issue in the field of energetic materials. The parameterization strategy described in this paper can be easily generalized toward other known HEs or the new HEs whose crystal structures can be obtained by CSP.

---

1. Introduction

Of particular importance is the development of new explosives of high performance and low sensitivity in the field of energetic materials. Their main performance characteristics such as detonation velocity and pressure usually are dependent on their densities. The velocity and pressure at the Chapman–Jouguet point increase with the initial density linearly and quadratically, respectively.¹ Density is the primary physical parameter in detonation performance.² Engineering and molecular designs for high explosives (HEs) usually consider firstly their densities at ambient temperatures or high pressures. Therefore, accurately predicting densities of HEs under such conditions is an important theoretical task.

The early density prediction method is by the so-called “group or volume additivity”.^3–6 Since the volume additivity cannot readily account for molecular conformation, crystal packing efficiency or positional isomerism,⁷ it has been gradually replaced by crystal structure prediction (CSP) since 1990s. Excellent reviews in CSP are available from five blind test studies hosted in 1999, 2001, 2004, 2007, and 2010 by the Cambridge Crystallographic Data Center (CCDC).^8–12 The low success rate in the previous three blind tests (CSP1999,⁸ CSP2001,⁹ and CSP2004 (ref. 10)) reveals just how challenging CSP can be. However, CCDC concluded from the CSP2007 (ref. 11) and CSP2010 (ref. 12) blind tests that a state-of-the-art method developed by Neumann et al.¹³ is able to reliably predict crystal structures of small rigid and slightly flexible molecules. The central parts of the method are the dispersion-corrected density functional theory (DFT-D) and a tailor-made force fields (TMFF) for every molecule under consideration. TMFF, derived from the reference data generated by DFT-D, is used directly for crystal structure generation, and DFT-D is used to predict the final energy ranking which is essentially limited to 10−100 equilibrium crystal structures obtained from TMFF.¹⁴ As Podeszwa pointed out, one of the key issues in CSP is the accuracy of force field.¹⁵

In designing and formulating energetic materials, one always aims at maximizing performance, and also attempts to minimize sensitivity. The energetic materials of low sensitivity can reduce the risk of accidental stimuli such as shock or impact, and thus improve the safety. Unlike density, the sensitivity is not a well-defined physical or chemical terminology, which is associated with a number of factors, including molecular composition and structure, crystal structure, microstructure, and even the way of external stimulus. So-called microstructure, a catchall term that refers to the three-dimensional structure with various imperfections and grain boundaries, is dependent on particle size and preparation technique. Therefore, strictly speaking, molecular structure, crystal structure, and microstructures as well as a stimulus way should be taken into account in estimating sensitivity. Recently, Wu and Huang¹⁶ have developed a thermo-mechanical model in which the single crystal accommodates plastic deformation through the dislocation motion on operative slip systems. It can account for the role contributed by bulk thermal heating to hot spot mechanism. Unfortunately, the model can't automatically decide which slip systems can be operative under shock loading. Therefore, searching the operative slip systems from many potential slip systems of single crystal in response to shock become a prior consideration in developing our sensitivity estimation model. Our research group is endeavoring to reliably search the operative slip systems under shock loading with the high quality of force field.

In addition, since the reaction rate in the detonation growth process or in detonation reaction region depends upon temperature, it is necessary to know how the liberated thermal energy is partitioned within an explosive to increase its temperature. Thus, to better estimate its performance, it is important to understand how its specific heat capacity changes with respect to both temperature and pressure. The accuracy of force field is also critical for the quantitative calculations of the temperature- and pressure-dependent relations of heat capacity.

Although it has been demonstrated that the off-the-shelf force fields developed during the past decades can be used to predict many properties of energetic solids,^17–22 they are empirical, i.e., they have been fitted to experimental data or their parameters have been simply transferred from other force fields. The empirical is responsible for the fact that their accuracy is not only rather limited in CSP for HEs, but also extremely insufficient for predicting their densities within an accuracy of 1% at low to room temperature or at ambient to high pressure. It can be expected that the empirical force fields will encounter numerous obstacles in developing sensitivity evaluation and detonation models. Addressing this issue necessitates referring to the accurate force field obtained with first principles.

It has been well-known that organic HEs tend to exhibit polymorphism dominated by intermolecular forces and the quality of such force fields is paramount to their CSP. For this reason, one has first focused on developing the non-empirical intermolecular force fields for HEs. In 2007, Podeszwa²³ has developed a wave-function based intermolecular force field for hexahydro-1,3,5-trinitro-1,3,5-triazine (RDX) dimers from the combined symmetry-adapted perturbation theory and density functional theory [SAPT(DFT)].^24–30 In 2011, using the SAPT(DFT) plus the Williams–Stone–Misquitta (WSM) method,^31–35 two of authors have derived the tailor-made intermolecular force fields for RDX and 1,3,5,7-tetranitro-1,3,5,7-tetraazacyclooctane (HMX), respectively.^36,37 Their accuracy has a good progress, and therefore leading to a success in CSP. This is a major advancement in the development of accurate potentials for HEs. However, since they are derived from the potential energy surfaces for RDX or HMX dimers, they are unable to describe well the many-body nonadditive interaction that plays an important role in condensed phase; and fail to describe the intramolecular interactions and thus not being applied into the real-world MD simulation where the molecular flexibility is required to be taken into consideration. The two deficiencies weaken their capability in precise density prediction and slip system search.

Obviously, it is pressing to develop all-atom, non-empirical, and tailor-made force fields (NETMFFs) for HEs, which can be competent for crystal and density predictions and also can be reliably applied into the sensitivity and detonation estimations. Because α-RDX has many theoretical and experimental data, it becomes a case study of developing NETMFFs from first principles.

2. Mathematical model

The mathematical function used in NETMFF is shown in eqn (1)

The potential energy function includes valence terms consisting of diagonal and off-diagonal cross-coupling terms, and nonbond interaction terms made up of a damped Buckingham potential for treating the van der Waals (vdW) interaction and a point-charge Coulombic function for calculating the electrostatic interaction. The valence terms describe internal coordinates of bond (b), angle (θ), dihedral (torsion angle, ϕ), out-of-plane bending angle (χ), and the cross-coupling terms that reflect coupling between the internal coordinates, such as the bond–bond and bond–angle couplings. The cross terms have been found to be important for predicting vibrational spectra and structural variations associated with conformational changes.³⁸ The nonbond interactions are used for describing the interactions between the intermolecular atom–atom pairs and between all intramolecular atom–atom pairs but the valence-bonded atom–atom pairs (1–2 interactions) and the atom–atom pairs (1–3 interactions) separated by two covalent bonds. In most force fields, the vdW repulsion interaction is still often described by r⁻ⁿ (n = 9 or 12) atom–atom terms as in the Leonard-Jones potential. However, it has been found³⁹ that r⁻ⁿ give a rather poor account of the repulsive region, and the exponent form [A exp(−Br)] is much better. To avoid the C₆/r⁶ divergent character that leads possibly to unphysical attractive forces at short intermolecular distances, a damping function (f_damp) expressed in eqn (2) is introduced into the Buckingham potential of eqn (1).

where, R is atomic vdW radii. The function is close to 1 for r_ij > 0.5(R_i + R_j), but quickly and continuously goes to zero when r_ij decreases. The Buckingham potential uses the following combination rules to obtain the heteroatom parameters from homoatom ones:

The atomic partial charges for the electrostatic interaction can be estimated from the bond-increments (δ_ij) acting as the charge parameters of NETMFF. It represents the charge separation between two valence-bonded atoms i and j. For atom i, its partial charge (q_i) is the sum of all charge bond increments:

j stands for all atoms that are valence-bonded to atom i.

Five atom types were defined for the RDX molecule: 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. According to our previous experience in developing the intermolecular potential for RDX,³⁶ the n_3o and n_3r types share a set of common vdW interaction parameters in the present force field.

3. Parameterization strategy

In contrast to the empirical force fields derived by indirectly probing the energy surfaces with experiment data, the NETMFF parameters are derived by directly measuring the energy surfaces from first principles. However, it is not a trivial task to fit the complex energy surfaces due to the large number of parameters involved. The large dimensionality of parameter space usually leads to the multiple solutions (i.e., combination of parameters) that can reproduce given energy surfaces and other selected properties in a similar way. The problem is referred to as “parameter correlation”.⁴⁰ An arbitrary ‘best-fit’ is not necessarily the right answer because of the parameter correlation effect that may make some functional terms physically unreasonable.

In addition to potential functions, some factors in parameterization strategy can influence the quality of NETMFF. Among these factors are the parameterization procedure, reference data, initial conditions, parameter optimization method, and figure-of-merit function. A good parameterization strategy can substantially reduce the parameter correlation effect and thus effectively prevent the force field from unphysical parameters.

Presented in Fig. 1 is a flow diagram illustrating the parameterization procedure. The whole process can be divided into two parts: the first principles calculations aiming at providing necessary data for the force field optimizations and self-consistent force field (SCFF) optimization procedure aiming at yielding the desirable force field parameters.

Fig. 1 Flow diagram of the NEMTFF parameterization.

Crystal structure is a starting point in the parameterization strategy. For a new explosive whose crystal structure is unknown, its crystal structure must be derived by the CSP method before carrying out the strategy. A detailed description of CSP for new HEs, one of the key issues in HE design technique, goes beyond the scope of this paper and will be discussed elsewhere.

3.1 DFT calculations of conformers

Based on the molecular structure in the RDX crystal, the conformer and dimer training sets were constructed firstly. The former was designed to yield the reference data, acting as the bases for the charge and valence parameter optimizations, and the latter, discussed in next section, was designed for yielding the pair-additive vdW parameters that act as initial values for the vdW parameter optimization.

To set up the conformer training set, at least one optimized structure must be first provided, which was derived by optimizing the isolated molecule in the crystal conformation at the B3LYP/6-31G* level of theory. In most cases, there are different low-energy conformations in the conformational space of an organic HE. Some of them had better be added into the training set in order to make the set better reflect the molecular potential energy surface near the optimized structure. Therefore, a conformational space of the RDX molecule was searched with the Simulated Annealing algorithm⁴¹ based on a hybrid method (the semiempirical PM3 plus B3LYP/6-31G*) to locate the low-energy conformations. It was found that the AAE isomer has the lowest energy at the level of B3LYP/6-31G*, followed by the AAA, a pair of twist enantiomers (twist and twist_e), EEA, and EEE conformers which are higher in energy by 1.2, 2.7, 5.5, and 22.7 kJ mol⁻¹, respectively. They are minima on their respective potential energy surfaces with no imaginary frequencies. Some important transition states (TSs) with one imaginary frequency, AAE-AAA, a pair of AAE-EEA enantiomers (AAE-EEA and AAE-EEA_e), EEA-EEE, a pair of EEA-twist enantiomers (EEA-twist and EEA-twist_e), and boat, were also found. The result has a little difference from the previous publication⁴² where the boat conformer was reported to be a stable conformer without imaginary frequencies. Fig. 2 illustrates the energies of the various conformers and TSs relative to AAE and the barriers to interconversion. The energies don't include zero-point corrections. The atom coordinates for these conformers and TSs are given in the ESI.†

Fig. 2 B3LYP/6-31G* potential energy surface for RDX interconversion.

The optimized molecular structure from the RDX crystal is exactly the same as the lowest energy conformer (AAE) in the conformational space. Fig. 2 shows that the interconversion between AAE and AAA can occur in the highest probability as compared with those between AAE and other conformers. Therefore, AAA became the member of the conformer training set. Total energies, gradients and the Hessian matrixes (first and second derivatives of energies) of AAE and AAA were calculated at the B3LYP/6-31G* level of theory.

In order to sample the potential energy surfaces outside the minimum energy states of AAE and AAA, their different conformations were also added into the conformer training set. They were generated in Sobol pseudo-random sequence⁴³ by rotating three –NO₂ groups, or by rotating the N–C–N–NO₂ dihedrals, or by bending three out-of-planes made up of the triazine's nitrogen and three atoms of the –NO₂ group bonded to the nitrogen, and or by bending another three out-of-planes made up of the two triazine's carbon atoms, the triazine's nitrogen atom bonded to the two carbon atoms, and the nitrogen atom of –NO₂ bonded to the triazine's nitrogen. The conformations that are not in minimum energy states but can reasonably describe the potential energy surfaces outside the minimum energy states at which the AAE and AAA conformers lies, just energies and gradients for these conformations were calculated also using B3LYP/6-31G*, the Hessian matrixes were not required.

The reference data including energies, gradients, and Hessians were provided for the basis of the valence parameter optimization. In addition, the atomic ESP fitted charges for each conformers in the training set, a basis for the charge parameter optimization, were derived from the B3LYP/6-31G* wave-functions using the CHelpG scheme.⁴⁴ All DFT calculations above-mentioned were performed with use of the Dalton program.⁴⁵

3.2 Parameters of the pair-additive vdW potentials

A long-standing problem in the development of intermolecular potentials comes from the lack of reliable and accurate vdW potential parameters. However, accurately characterizing the vdW interactions using quantum mechanics method itself is not a trivial task. Especially for the organic HEs, such calculations are very challenging for a few of reasons. Firstly, a high level of theory and an extensive basis set must be used to calculate vdW interactions accurately. However, high level of theories like CCSD(T) and MP4, which can perform accurate calculation of the overall intermolecular energies between two interacting molecules, are computational demanding even for single point calculation. Secondly, In fact, it is almost impossible to perform the calculations of very high quality correlated wave functions for the dimers at thousands of relative orientations under the modern computational conditions. Thirdly, the CCSD(T) and MP4 methods are also not ideally suited for derivation of the vdW parameters due to their inabilities to decompose the overall interaction energy into physically significant contributions. This makes it hard to parameterize the analytic vdW potential functions composed of the dispersion and repulsion. Lastly, even if the ab initio calculations can be carried out accurately, the obtained energetic data only represents molecular interactions in the gas dimers, which can be significantly different than those found in the condensed phases.

In recent years, it has been found^36,37,46,47 that the hybrid method consisting of the SAPT(DFT) and the WSM method, as a good alternative to the CCSD(T) and MP4 methods, provides a feasible framework for the development of the pair-additive vdW potentials because it provides the interaction energy as a sum of physically significant contributions. This allows the separate parameterization of attractive and repulsive terms in the vdW potential function. It has been proved by our recent work^36,37 that the hybrid method can effectively address the first three issues mentioned above. Hence, solving the last issue is the key to the present work, which is discussed in the Section 3.3. A detailed description of how to derive the parameters of pair-additive vdW potentials, acting as initial values for the SCFF optimization, is presented here.

A set of previously reported vdW parameters for RDX from the hybrid method was coupled with the distributed multipole model,³⁶ and they can't be used here due to adopting the point-charge model in eqn (1). But the basic strategy for parameterizing vdW parameters used in the previous report³⁶ was applied in the present work.

The vdW interaction energy (E_vdW) for a HE dimer can be accurately calculated as follows:³⁷

E_vdW = E⁽¹⁾_elest(KS) − E_coul-class + E⁽¹⁾_exch(KS) + E⁽²⁾_ind,d-class + E⁽²⁾_disp,d(C₁₂) (5)

where, E⁽¹⁾_elest(KS) and E⁽¹⁾_exch(KS) are the first-order SAPT(DFT) electrostatic and exchange-repulsion energies, respectively, and E_coul-class is the coulombic energy calculated using the classical point-charge model. The first three terms constitute short-range repulsion energy (E_sr) that decays exponentially with increasing separation. The last two energy terms are damped induction and dispersion that were computed using the WSM distribution method. The two terms constitute the long-range attraction energy (E_lr) that can be described by a simple dispersion model (C₆/r⁶). The rank of the anisotropic dispersion coefficients included in eqn (5) is up to C₁₂.

Before constructing dimer training sets, if possible, the molecule had better be symmetrized as to be its potential point group of highest symmetry. This can cause greatly speeding up the WSM calculations. The molecular structure of RDX, taken from the experimental crystal structure,⁴⁸ belongs to the C₁ point group, but was found to be an approximated symmetry of C_s. Thus, we symmetrized the molecule as the C_s point group. The root-mean-square displacement of atomic coordinates resulted from the symmetrization was found to be less than 0.073. Therefore, the difference between the two structures is very little.

Considering a higher computational demanding of the first SAPT(DFT) energies than that of the WSM attraction energies, two training sets were set up based on the symmetrized molecule: one set that comprises 1000 dimers is for calculating ab initio repulsive energies (E_sr), and the other, comprising 2000 dimers, is for calculating ab initio attractive energies (E_lr). These dimers in the sets were generated using the same random algorithm as reported in previous publication.³⁶ Their intermolecular separations for each orientation were restricted to [R₀ − 1.5, R₀ + 1.2] a.u., R₀ is a molecular vdW contact radius that can be estimated from the Bondi vdW radii.⁴⁹

Before calculating the energies for RDX dimers, an asymptotically corrected (AC) PBE0 (ref. 30)/aug-Sadlej³⁶ wave function of the symmetrized RDX molecule was derived, which is necessary for the WSM and first-order SAPT(DFT) calculations. The distributed WSM static and frequency-dependent polarizabilities of rank 3 for the symmetrized molecule, which are essential for calculating the damped inductions and the distributed anisotropic dispersion coefficients that are used to calculate the damped dispersions, were obtained with the use of the AC PBE0/aug-Sadlej wave function. The damped induction and dispersion energies were calculated by applying the Tang–Toennies damping functions⁵⁰ with the damping parameter β = 21.777 a.u.³⁶ The atomic partial charges, used to compute the Coulombic energies (E_coul-class), were estimated from the charge parameters. These 1000 ab initio short-range energies and 2000 ab initio long-rang energies were obtained with the CamCASP⁵¹ and Orient,⁵² respectively.

Directly fitting ab initio long-rang or short-range energies is unwise and thus not recommended, since it leads easily to unphysical terms. Instead, we used two different multistage procedures^46,47 to obtain the attractive and repulsive parameters, respectively, in order to effectively keep off unphysical parameters.

The derivation of the C_6,add parameters is by a three-stage procedure:

Firstly, compute the distributed WSM frequency-dependent polarizabilities of rank 1 for all the atoms in the symmetrized RDX molecule.

Secondly, applying the combination rule of eqn (6), calculate the isotropic asymptotic dispersion parameters (C_6,asy) for four homoatom interaction pairs (C–C, H–H, N–N, and O–O) with the calculated frequency-dependent polarizabilities. This makes each C_6,asy parameters (Table 1) physical. However, the asymptotic dispersion model can't recover well the ab initio attractive energies.

Table 1 Parameters of the pair-additive vdW potentials

Atom pair	Aii,add (kcal mol^−1)	Bii,add (Å^−1)	C^ii6,add^a (kcal mol^−1)	C^ii6,asy (kcal mol^−1)
a S6 = 1.4271.
c_4–c_4	372 239.5015	4.6598	865.7856	606.6748
h_1–h_1	920.0656	3.7147	24.6642	17.2827
n_3r–n_3r/n_3o–n_3o	294 861.8691	4.7461	118.4293	82.9860
o_1–o_1	137 352.0936	4.0568	529.9247	371.3298

---

Finally, introduce a scaling factor by which the C_6,asy model dispersion energies can be scaled to recover the 2000 ab initio attractive energies. In order to find the scaling factor (S₆), a function of the following form is minimized:

where, n labels the dimer configurations comprising molecule A and B, E_lr(n) is the ab initio attractive energy for configuration n, atom i and j belong to molecules A and B, respectively, r_ij is the separation of the i–j pair. S₆ is determined by a least-squares fit. Thus, the final parameters (C^ij_6,add) of the pair-additive dispersions can be simply estimated by

C^ij_6,add = S₆C^ij_6,asy (7)

In order to derive the reasonable parameters (A_ij,add and B_ij,add) of the pair-additive repulsions, we introduced the anisotropic energies (E_sr,aniso) and performed the multistage procedure:

Firstly, compute the electron densities (ρ^A_e andρ^B_e) of molecules A and B with the AC PBE0/aug-Sadlej wave functions, and then partition them into their respective atomic contributions by the density-fitting method, that is, .

Secondly, derive the general parameter (K₀) of the overlap model by the least-squares fit:

where, S_ρ(n) is the density overlap between A and B in the dimer configuration n, which can be estimated by S_ρ = ∫ρ^A_eρ^B_edV. w_n is a weight that depends on the energy E_sr.

Thirdly, derive the distributed overlap model that can be described by the expression:

Where, i and j are the atoms belong to molecules A and B, respectively. S_ρ(ij) is the density overlap between atom i and j, which can be estimated from the derived ρⁱ_e and ρ^j_e, i.e., S_ρ(ij) = ∫ρⁱ_eρ^j_edV. To avoid the non-physical values of K_ij, a constraint minimization of the following function were carried out:

Fourthly, derive the anisotropic atom–atom repulsion potential model that can be represented by the Born–Mayer functions⁵³

G is a constant energy unit taken to be 10⁻³ a.u. α_ij describes the hardness of the interaction of atom pair i–j. ρ_ij describes the shape of the interaction between atoms i and j as a function of their relative orientation (Ω_ij), and it can be expressed as S̄-function expansion⁵³where the coefficients
depend on the atom–atom separation, and
are the non-normalized set of orthogonal orientation dependent functions developed by Stone.⁵³ Here, the S̄-functions are limited to those for which l_a + l_b = j, and the expansion converges rapidly, so that only a few terms are normally needed.

We fitted the contributions [E_sr(ij) = K_ijS_ρ(ij)] of individual atom–atom pairs to a single Born–Mayer term in eqn (11) via the obtained distributed overlap model. In this way, each fit involves only a small number of adjustable parameters and avoids the difficulties of fitting a large number of strongly coupled parameters. This fit procedure is relatively easy and ensures to derive the physically sensible parameters for each atom pairs. The Born–Mayer repulsion parameters are given in the ESI.†

ρ_ij(Ω_ij) can correctly represent the atom–atom interactions where the electron distribution around the atoms in molecule don't own spherically symmetry any more. Thus, the anisotropic Born–Mayer potential can more accurately recover E_sr than the isotropic potential. On the other hand, although the Born–Mayer potential is anisotropic, it closely approaches the isotropic exponent function in an explicitly mathematical way. For the two reasons, the Born–Mayer potential is an ideal bridge over the gap between E_sr and simple exponent potential, and thus provides a good insurance against the non-physical parameters of the exponent function.

Finally, compute the 2000 anisotropic repulsion energies using the obtained atom–atom Born–Mayer potential, and fit them to the exponent potential functions by applying the combination rules of eqn (3). The obtained values for α_ij and the isotropic part of ρ_ij(Ω_ij) act as the initial inputs of the fitting. Fitting results are given in Table 1.

Before employed in the next procedure, the parameters of the pair-additive vdW potential have to be validated by a CSP procedure where 230 space groups were searched through based on the molecule structure taken from the experimental crystal structure.⁴⁸ It was found that the global minimum structure is exactly the same as the experimental structure. The energy ranks of low energy structures for these space groups are given in the ESI.†

3.3 Periodic calculations of crystal with a tailor-made DFT-D model

Although in principle one can derive the vdW parameters by fitting to ab initio data calculated for molecular dimers, this approach does not work well for modeling molecules in condensed phases, because the pair-additive potential constructed from the dimers can't accurately represent molecular interactions in condensed phases. In parameterizing vdW terms of empirical force fields, a solution to this problem is fitting simulated thermophysical properties of molecules in condensed phases with experimental data. Obviously, the solution is empirical, and thus fails to cater for the aim of the present work to develop a non-empirical force field. Instead of this, a non-empirical solution is adopted by the procedure:

(1) Optimize the ab initio RDX crystal structure using the periodic DFT-D method;

(2) Compute the density (ρ₀) and intermolecular cohesive energy (ΔE^coh_inter) for the optimized crystal with the periodic DFT-D method. The latter are computed as:

ΔE^coh_inter = E(bulk)/Z − E(MIC) (13)

where E(bulk) is the total energy of the unit cell and Z is the number of molecules in the unit cell, and E(MIC) is the total energy of the molecule in the crystal.

(3) ρ₀ and ΔE^coh_inter act as the basis of parametrizing the vdW potential that can accurately represent the interactions of condensed phases, because the two physical quantities implicitly contain the many-body nonadditivity.

In ten years, there are a few DFT-D methods^54–57 developed for correcting the dispersion defect of density functionals that has been found in handling molecular crystals or clusters. The most well-known scheme among them is Grimme's approach⁵⁴ that has been widely inserted into the many programs and also successfully applied in energetic molecular crystals. However, these DFT-D methods can't achieve an enough high accuracy requested by the present work. Civalleri⁵⁸ has found that modifying the original Grimme model settings, S₆ and atomic vdW radii, can lead to better balanced results for both lattice constants and cohesive energies as well as for hydrogen bonded and dispersion bonded molecular crystals. Table 2 shows that for the C–H–N–O systems, the scaled Grimme's vdW radii by the factor of 1.05 used by Civalleri for heavy atoms are found to be close to the Bondi's ones.⁴⁹ The Bondi's radii have been arrived at by selecting from the most reliable X-ray diffraction data those which could be reconciled with crystal density at 0 K.⁴⁹ In fact, the original Grimme's radii⁵⁴ are designed for gas clusters, not for crystals.

Table 2 The three different atomic vdW radii

	Grimme's radii	Scaled Grimme's radii by 1.05	Bondi's radii
C	1.452	1.525	1.53
H	1.001	1.051	1.06
N	1.397	1.467	1.46
O	1.342	1.409	1.42

---

Therefore, using the Bondi's radii is more reasonable in the periodic Grimme DFT-D calculations. It was further found by us^36,37 that a tailor-made Grimme BLYP-D model can yield very precise densities and lattice energies for RDX and HMX crystal if the C^ij_6,add parameters and the S₆ coefficient as well as the Bondi's radii are used. This is because these parameters have greatly improvement over those in original Grimme's model in describing dispersion interactions. The RDX experimental crystalline structures were fully optimized at the BLYP/6-31G(d, p) level using the tailor-made BLYP-D, and then basis-set-superposition-error (BSSE) corrected intermolecular cohesive energy (ΔE^coh_inter) and intramolecular deformation energy [ΔE^def_intra = E(MIC) − E(mol)] were also calculated at the same level of theory, where E(mol) stands for the total energy of the free molecule. All periodic calculations were performed using Crystal 09 program.⁵⁹ Thus, we estimated the sublimation enthalpy (ΔH_sub) at 298 K according to the following equation.

where the term of –2RT represents an approximate correction for the difference between the gas phase enthalpy approximated by PV + 3RT (PV = RT), and the vibrational contribution to enthalpy for solid molecular crystal, which can be estimated with PV + 6RT (PV ≈ 0). The last term is the energy of zero-point vibrations in the Debye approximation, and the Debye temperature (Θ=198.3 K) in it can be obtained by the accurate thermal expansion MD simulations described below.

Table 3 collects these energies and density of the optimized RDX crystal. The estimated sublimation enthalpy is almost equal to the experimental value⁶⁰ at 298 K. This provides a good reason for believing that ρ₀ and ΔE^coh_inter can implicitly contain the true nonadditivity in the RDX crystal.

Table 3 The density at 0 K and the sublimation enthalpy at 298 K

ρ0 (g cm^−3)	ΔE^cohinter (kJ mol^−1)	ΔE^defintra (kJ mol^−1)	ΔHsub (kJ mol^−1)
			This work	Exp.
1.895	−155.4	13.8	134.7	134.3

---

3.4 Parameter optimization

As shown in Fig. 1, the NETMFF parameter optimization comprises charge parameter optimization, vdW parameter optimization, and valence parameter optimization. In principle, the calculated ab initio energies, first and second order derivatives of energies, and electrostatic potential (ESP) derived charges for the conformers in the conformer training set, as well as ab initio ΔE^coh_inter and ρ₀ for the optimized crystal structure are used to derive charge, vdW, and valence parameters by minimizing the least square methods:where, y^QM_i denotes ab initio reference data, y^FF_i(X̄_i;{k₁, k₂,…,k_M}) is the corresponding data computed with force field, X̄_i represents atom coordinates, and {k₁, k₂,…,k_M} are force field parameters. The “weighting factors” w_i normalize the quantities. Normally for valence and charge parameter optimizations, the number of data points (N) is much greater than the number of variables (M). The minimization algorithm is essentially a combination of Gauss–Newton and Levenberg–Marquardt methods.

It has been known⁵³ that the A_ij parameters are associated closely with the shapes of atoms i and j. The nonadditive effect of condensed phase always affects the atom shapes and thus the A_ij parameters. By contrast, the B_ij parameters, referred to as the repulsion hardness, have been found to fall in a quite narrow range,⁵³ and thus are almost independent of the nonadditive effect. On account of this, the hardness parameters were fixed to the values of B_ij,add during the vdW optimization process. Additionally, the pair-additive vdW potential has been validated by CSP. This means that the proportions of A_ij,add and C^ij_6,add between different atom–atom pairs are in reasonable balance. In this way, only two scale factors (k_A and k_C) of the A_ij,add and C^ij_6,add parameters need adjusting independently in the optimization.

The NETMFF parameterizations were carried out using the SCFF optimization scheme as shown in Fig. 1.

Firstly, the charge parameter optimization is relatively straightforward compared to the valence and vdW parameterizations. They are derived using a constrained ab initio ESP charge fit. Once deriving the bond increments, the partial atomic charges are exclusively determined by the atom types and defined via them. The bond increments are always fixed during the whole parameterization process.

Secondly, the vdW parameter optimization procedure is performed using initially guessed values of A_ij,add and C^ij_6,add. In order to minimize formula (15), the first order derivatives of the target values forρ₀ and ΔE^coh_inter, with respect to each of the adjustable parameters (k_A and k_C) need to be evaluated. According to the derivatives, the program automatically adjusts k_A and k_C. Based on previously MM-optimized target values, it estimates “on-the-fly” criterions forρ₀ and ΔE^coh_inter at each vdW parameter optimization stage (note: the criterions at the first optimization stage are those obtained from the periodic DFT-D optimization, see Table 3). Once the target values, obtained from the new vdW parameters, coincide with the corresponding criterions within the allowed errors, this procedure gives a set of new vdW parameters for the following valence parameter optimization.

Thirdly, given the charge and vdW parameters, the valence parameter procedure is based on the reproduction of the ab initio data including energies as well as first and second order derivatives of energies obtained from the conformer training set. Via an iterative approach, the bond, angle, and torsion parameters were primarily optimized with the initial values estimated by the UFF method.⁶¹ The parameters decide the molecular structures, and therefore have a strong impact on other properties. Then, the out-of-plane and cross coupling parameters are optimized in sequence to better reproduce the ab inito reference data. The procedure involves the iterative process where, upon changing one class of parameters, a set of previously optimized parameters are refined again if necessary. This approach yields a valence parameter set that accurately reproduces the ab initio reference data.

Lastly, the derived force field is used to fully optimize crystal structure and compute the target values of ΔE^coh_inter and ρ₀ for the optimized structure. The differences between the target values and the reference values are estimated. According to the differences, the SCFF program decides either to complete the parameter optimization or to do the next SCFF loop with previously optimized force field parameters.

All parameter optimizations can be performed in a suite of SCFF codes. The final NETMFF parameters for RDX molecule are given in Table 4.

Table 4 The NETMFF parameters for RDX

Bond		b0	kb1	kb2	kb3
n_3r	c_4	1.4991	263.1539	526.3078	614.0258
n_3r	n_3o	1.3701	350.6097	701.2194	818.0893
c_4	h_1	1.0908	391.7927	783.5854	914.1830
n_3o	o_1	1.1999	870.2045	1740.4090	2030.4772

---

Angle			θ0	kθ1	kθ2	kθ3
c_4	n_3r	c_4	109.5528	106.0244	20.5205	22.5493
c_4	n_3r	n_3o	108.3641	86.5467	15.6410	18.2999
n_3r	c_4	n_3	106.7516	76.2406	12.4776	16.0040
n_3r	c_4	h_1	100.9465	55.2234	5.7869	11.3586
h_1	c_4	h_1	103.8408	43.0038	5.7496	8.9255
n_3r	n_3o	o_1	150.5536	35.0371	36.5600	12.9724
o_1	n_3o	o_1	155.3308	43.2967	56.3925	18.5463

---

Dihedral				kϕ1	kϕ2	kϕ3
c_4	n_3r	c_4	n_3r	21.5972	13.2134	4.6374
c_4	n_3r	c_4	h_1	2.4911	0.1266	0.0041
n_3o	n_3r	c_4	n_3r	0.1931	2.4156	0.5234
n_3o	n_3r	c_4	H_1	3.7386	0.1228	0.1798
c_4	n_3r	n_3o	o_1	1.8245	2.5493	0.4440

---

Out-of-plane				kχ
c_4	n_3r	c_4	n_3o	0.9973
n_3r	n_3o	o_1	o_1	3.5010

---

Bond–bond			kbb′
c_4	n_3r	c_4	10.8266
c_4	n_3r	n_3o	75.3645
n_3r	c_4	n_3r	131.2915
n_3r	c_4	h_1	32.1571
h_1	c_4	h_1	26.6591
n_3r	n_3o	o_1	172.0110
o_1	n_3o	o_1	218.2496

---

Bond–angle			kbθ1	kbθ2
c_4	n_3r	c_4	29.6257	29.6257
c_4	n_3r	n_3o	34.7088	52.2388
n_3r	c_4	n_3r	54.4618	54.4618
n_3r	c_4	h_1	40.4474	13.5139
h_1	c_4	h_1	1.4994	1.4994
n_3r	n_3o	o_1	32.5971	42.7624
o_1	n_3o	o_1	56.0449	56.0449

---

Bond increment		δ
h_1	c_4	0.0829
n_3r	c_4	−0.1999
n_3o	n_3r	0.0097
o_1	n_3o	−0.3698

---

vdW		B	A	C6
c_4	c_4	4.6598	410318.0391	931.9609
h_1	h_1	3.7147	1014.1844	26.5494
n_3r	n_3r	4.7461	325024.9999	127.4813
o_1	o_1	4.0568	151402.6359	570.4289

---

4. Results

4.1 Gas phase molecular geometries and properties

The six stable conformers were optimized with NETMFF, and then their vibration analyses were also performed using the same molecular mechanics (MM) force field method. The corresponding data at the level of B3LYP/6-31G*, referred to as the QM reference data, were available for the validation of NETMFF. Fig. 3(a)–(d) illustrate charts of correlations between the MM-calculated and QM reference data for bond lengths, angles, dihedrals, and vibration frequencies, respectively.

Fig. 3 Comparison between force-field calculated results and the QM reference values.

There is a total of 125 data points for bond lengths, 216 for bond angles, and 288 for dihedrals. The data set of bond lengths ranges from ca. 1.08 to 1.47 Å. The maximum absolute percentage deviation is 2.8%, the average percentage deviation is 0.1%, and the root mean squares (rms) deviation is 0.8%. For the angles, their data set ranges from ca. 107° to 127°. The maximum absolute percentage deviation, the average percentage deviation, and rms are 6.5%, −0.3%, and 1.7%, respectively. As illustrated in Fig. 3(a) and (b), excellent agreement between the calculated and reference data for bond lengths and angles is obtained. As can be seen from Fig. 3(c), deviations obtained for the dihedrals are obviously higher than those obtained for the bond lengths and angles. The average percentage deviation and rms are 2.8% and 13.5%, respectively. The maximum absolute percentage deviation is up to ca. 83%, which is ascribed to too small torsion angles of 5.4° and 9.9° for the QM-optimized and MM-optimized twist geometries, respectively. Therefore, the force filed also can give reasonable values for torsions. Fig. 3(d) shows that the agreement of molecular vibration frequencies between the QM-optimized and MM-optimized geometries of the six conformers is reasonably good. There are 342 data points, the average deviation is −8.36 cm⁻¹, the rms deviation is 42.6 cm⁻¹, and the maximum absolute deviation is 133.1 cm⁻¹. These deviations are considered to be acceptable for the force-filed based calculations. We also performed the vibration frequency calculations for the QM-optimized transition states with NETMFF, showing that all of them are only one imaginary frequency. Therefore, the force field can be reliably used to vibration frequency calculation.

In addition, the O₂N–N–C–C out-of-plane bending angle (χ_{N–N–C–C}), as shown in Fig. 4, was reviewed, which describes the tilt of the N–N bond out of the plane formed by the C–N–C ring atoms. The most significant difference in the molecular structures in the crystalline and gaseous states is the orientation of the NO₂ groups relative to the plane of the ring, which can be described by the structural parameter χ_{N–N–C–C}. There are two types of NO₂ groups in the conformers: one belongs to A type whose NO₂ group is in a pseudo-axial (A) arrangement relative to the ring, the other belongs to E type whose NO₂ group is in a pseudo-equatorial (E) arrangement relative to the ring. Table 5 lists the structural parameters for the QM-optimized and MM-optimized geometries of the six conformers.

Fig. 4 Illustration of the out-of-plane bending angle (χ_{N–N–C–C}).

Table 5 The χ_{N–N–C–C} angles for the QM-optimized and MM-optimized geometries of the six conformers^a

AAA		AAE		EEA		EEE		Twist/twist_e
a Italic is for E type of NO2 group, bold is for A type of NO2 group.
MM	QM	MM	QM	MM	QM	MM	QM	MM	QM
31.3	30.1	34.7	34.0	38.5	40.6	46.4	34.0	38.6	33.5
31.2	30.2	34.6	34.0	46.8	35.2	46.4	34.0	41.4	33.0
31.3	30.2	47.4	37.9	46.8	35.2	46.5	34.0	37.9	36.7

---

Table 5 lists the structural parameters for the QM-optimized and MM-optimized geometries of the six conformers. The A type absolute deviation (χ_MM–χ_QM) for EEA, AAE, AAA, twist/tiwst_e are −2.1°, 0.7°, 1.2°, 5.1°, respectively, and the E type absolute deviations for twist/tiwst_e, AAE, EEA, EEE are 8.4°, 9.5°, 11.6°, 12.5°, respectively. Obviously, The ability of NETMFF to describe the A type and E type χ_{N–N–C–C} angles is quite different, and the former is much better than the latter. It was also noted that the number (N_E) for the E type χ_{N–N–C–C} angle increases with the sequence: AAA (0), twist/twist_e (1), AAE (1), EEA (2), EEE (3).

In fact, the different ability of NETMFF in describing the angle and torsion with these NO₂ groups was also found. For example, the maximum absolute percentage deviation (6.5%) in Fig. 3(b) belongs to the three N–C–N angles on the chair ring, they are bonded to the E type NO₂. A more detailed analysis about Fig. 3(b) was performed. The three N–C–N angles of AAA are bonded to the A type NO₂ and their absolute percentage deviation (APD) are 0.6%; For AAE, APD of the N–C–N angle bonded to the A type NO₂ is 0.6%, APDs of two other angles bonded to the A and E type NO₂ are 3%, the second maximum absolute percentage deviation in Fig. 3(b); For EEA, there are one N–C–N angle bonded to the E type NO₂ and its PAD is 6.2%, and the others bonded to the A and E type NO₂ and their PADs are 3%. The E type NO₂ group is responsible for large deviation in Fig. 3(b). Obviously, the higher N_E, the larger deviation. A similar conclusion for torsion can also be attained. Since the two types of NO₂ was designed to share common force field parameters, the force field can't reflect well the different influence of different spatial orientation of NO₂ on the relevant angle bending, torsion, and out-of-plane bending.

We further analyzed the maximum distance (MAXD) and the root-mean-square distance (RMSD) between two equivalent atoms belonging to the QM-optimized and MM-optimized geometries, respectively, and the degrees of overlay between them. Fig. 5 illustrates the analyzing results for the degrees of overlay, N_E, MAXD, and RMSD. With N_E increasing, the MAXD and RMSD increase but the degree of overlay decreases. In order to analyze the influence of N_E on the energies, we calculated the two classes of energies with NETMFF for comparison. One class is the energies of the MM-optimized conformers relative to the MM-optimized AAA, the other is those of the QM-optimized conformers relative to the QM-optimized AAA. We converted the QM energies listed in Fig. 2 into the relative energies for comparison. These relative energies are collected in Table 6. Fig. 5 and Table 6 show that an increase in N_E leads to the increase in difference between the MM-optimized and the QM-optimized geometries, and thus the increase in difference of energy. Among these MM-optimized geometries, the EEE one has a maximum deviation from its QM-optimized geometry, but it still takes on a reasonable structure that is similar to the QM-optimized. In other words, NETMFF can yield the correct geometry for every RDX conformer within a reasonable deviation.

Fig. 5 Illustration of degrees of overlay, as well as NE, MAXD, and RMSD for the conformers (red for the QM-optimized and green for the MM-optimized).

Table 6 Relative energies for the QM-optimized and MM-optimized conformers (kJ mol⁻¹)

		AAA^a	AAE	EEA	EEE	Twist/twist_e
a Since the force-field based energies for QM-optimized and MM-optimized AAA geometries are 76.4 and 76.2 kJ mol^−1, respectively, the relative energies for AAA are set to zero.
QM-optimized	QM calculation	0.0	−1.2	4.3	21.5	1.5
	MM calculation	0.0	−1.2	5.9	19.5	5.3
MM-optimized	MM calculation	0.0	−7.7	−11.4	−10.1	−5.0

---

Based on the validated intramolecular properties, we come to the conclusion that NETMFF can reasonably describe the intramolecular potential surface of the gas RDX molecule.

4.2 Molecular Crystal's properties

To validate the capability of the present force field in describing crystal properties, the fully NETMFF-based optimization of crystal structure of RDX was performed, and then ρ₀, ΔE^coh_inter, and the Γ-point frequencies for the MM-optimized crystal structure were also calculated. As shown in Fig. 6 and 7 as well as Table 7, they are compared with those obtained from the tailor-made DFT-D Model at the 6-31G**/BLYP level of theory. These calculations are used to evaluate if the depths and positions of potential wells are in a good state, which can be reflected by ΔE^coh_inter and ρ₀, respectively. The rms deviation of the geometrical parameters including cell parameters and atomic positions between the two optimized structure is estimated to be 0.2534. The average deviation for 504 Γ-point frequencies is 35.7 cm⁻¹, the rms deviation is 42.2 cm⁻¹, and the maximum absolute deviation is 189.5 cm⁻¹. These frequency deviations are similar to those estimated for molecular vibration frequencies. The results suggests that the structural parameters (cell parameters and atomic positions) and properties (density, energy, and frequencies) obtained from NETMFF are in good agreement with those obtained from the tailor-made DFT-D model.

Fig. 6 The overlay between the QM-optimized (red) and MM-optimized (green) structures.

Fig. 7 A comparison of the Γ-point frequencies between the QM-optimized and MM-optimized structures.

Table 7 A comparison of cell parameters, ρ₀ and ΔE^coh_inter between the QM- and MM-optimized structures

	a (Å)	B (Å)	C (Å)	α (°)	β (°)	γ (°)	ρ0 (g cm^−3)	ΔE^cohinter (kJ mol^−1)
MM-optimized	13.2120	11.2409	10.4814	90.0	90.0	90.0	1.895	−155.2
QM-optimized	13.6585	10.9818	10.3789	90.0	90.0	90.0	1.895	−155.4

---

A series of NPT MD simulations of the isothermal compressions at 298 K was carried out at the same pressures as the experimental measurements.⁶² It has been known that the α-RDX polymorph is stable below 4 GPa.⁶² The purpose for these simulations is to evaluate if the present force field is able to accurately describe the repulsion interactions. These simulations started with a 3 × 3 × 3 supercell built from the MM-optimized structure. The simulated results are tabulated in Table 8 along with experimental data for comparison. The maximum density discrepancy occurs at 0.73 GPa, but the magnitude of the error is only 0.66%. The errors at other pressures are less than 0.3%. The average percentage deviation of cell lengths is almost null, the rms percentage deviation is 3.4%, and the maximum percentage deviation is 5.2%. Based on the compression simulations, the Murnaghan equation of state⁶³ gives the reasonable values for the bulk modulus (B₀ = 13.36 GPa) and the pressure derivative of bulk modulus (B^′₀ = 6.60). They excellently agree with the corresponding experimental values (B₀ = 13.0 GPa andB^′₀ = 6.6).⁶²

Table 8 A comparison of 298 K isothermal compressions between the simulations and experiments

Pressure
Pressure (GPa)	a (Å)		b (Å)		c (Å)		Density (g cm^−3)
	Simu.	Expt.	Simu.	Expt.	Simu.	Expt.	Simu.	Expt.
0	13.890	13.200	11.168	11.600	10.555	10.720	1.802	1.798
0.73	13.655	13.020	10.979	11.420	10.376	10.530	1.897	1.885
1.75	13.446	12.890	10.811	11.220	10.217	10.270	1.987	1.987
2.75	13.295	12.800	10.690	11.070	10.103	10.160	2.055	2.050
3.36	13.221	12.740	10.630	11.000	10.046	10.080	2.090	2.089
3.95	13.157	12.660	10.578	10.930	9.998	10.030	2.121	2.125

---

These results indicate that agreement with the experiment values on hydrostatic compression is quite good for the force field and validate that the fore field can accurately describe the intermolecular repulsion interactions in the RDX condensed phase.

We also carried out another series of NPT MD simulations based on the 3 × 3 × 3 supercell to derive thermal expansion data at 1 atm pressure, which is used to evaluate if the present force field is able to accurately describe the attraction interaction in condensed phase. The simulated results are summarized in Table 9 along with the X-ray powder diffraction experimental data (on RDX with purity of 99.9%)⁶⁴ for comparison. The maximum density discrepancy occurs at 90 K, but the relative error is 0.71%. The errors at ambient temperatures or higher are less than 0.5%. The average percentage deviation of cell lengths is very small (∼0.02%), the rms percentage deviation is 4.0%, and the maximum percentage deviation is less than 6.0%.

Table 9 A comparison of 1 atm pressure thermal expansion data between the NPT simulations and experiments

Temperature (K)	a (Å)		b (Å)		c (Å)		Density (g cm^−3)
	Simu.	Expt.	Simu.	Expt.	Simu.	Expt.	Simu.	Expt.	Debye
a The experimental data at 90 K are from ref. 65.
90.15^a	13.718	13.140	11.030	11.419	10.424	10.586	1.871	1.858	1.879
100.15	13.727	—	11.037	—	10.431	—	1.867	—	1.876
150.15	13.765	—	11.068	—	10.460	—	1.852	—	1.857
200.15	13.808	—	11.102	—	10.492	—	1.835	—	1.837
250.15	13.847	—	11.133	—	10.522	—	1.819	—	1.817
273.15	13.867	—	11.149	—	10.537	—	1.811	—	1.808
298.15	13.888	—	11.166	—	10.553	—	1.803	—	1.798
303.15	13.893	13.200	11.170	11.609	10.557	10.724	1.801	1.796	1.796
323.15	13.914	13.209	11.187	11.629	10.573	10.743	1.793	1.788	1.788
343.15	13.939	13.216	11.207	11.648	10.592	10.760	1.783	1.781	1.780
363.15	13.953	13.223	11.218	11.667	10.602	10.778	1.778	1.775	1.772
383.15	13.978	13.231	11.239	11.687	10.622	10.798	1.768	1.767	1.764
403.15	14.001	13.238	11.257	11.704	10.639	10.817	1.760	1.761	1.756
423.15	14.022	13.250	11.274	11.726	10.655	10.841	1.752	1.752	1.748
443.15	14.046	13.257	11.293	11.743	10.673	10.862	1.743	1.745	1.740

---

Within the temperature ranging from 330–443 K, the volume coefficients of thermal expansion (CTE) obtained from the simulations is 2.37 × 10⁻⁴ K⁻¹, which is very close to the experimental value of 2.07 × 10⁻⁴ K⁻¹.⁶⁴ It has been well known that the asymmetric shape of the intermolecular potential causes thermal expansion. These results show that the present force field can accurately recover the asymmetric shape decided by intermolecular attractions and repulsions. However, as compared with the anisotropic experimental linear CTEs along the a-, b- and c-axes (3.07 × 10⁻⁵, 8.28 × 10⁻⁵, and 9.19 × 10⁻⁵, respectively),⁶⁴ the corresponding predicted values (7.85 × 10⁻⁵, 7.82 × 10⁻⁵, and 7.82 × 10⁻⁵, respectively) tend toward isotropy. Our previous simulations using the intermolecular SPOT potential,³⁶ where the molecular geometry was kept rigid and same as the experimental one,⁴⁸ give the values of 3.07 × 10⁻⁵, 8.28 × 10⁻⁵, and 9.19 × 10⁻⁵ K⁻¹, respectively, which take on the similar anisotropy in experiment. The isotropic linear CTEs illustrate the valence parts of the present force field can achieve reasonable but not high precision in describing intramolecular interactions.

Despite this, the NETMFF for RDX, consisting of the reasonable valence parameters and the accurate intermolecular parameters, has a very good ability to predict densities as shown in Tables 8 and 9. This suggests that the NETMFF approach become a new technical means of great promise for predicting accurately densities of high explosives at the temperatures and pressures that they are often subjected to, which is a long-standing issue in the field of energetic materials.

According to the expansion simulations, we obtained the Debye expansion model for RDX:

where, D(x) is the Debye function, Θ is Debye temperature of 198.3 K, α_Θ = 2.28 × 10⁻⁴ K⁻¹. In the Debye approximation, the limiting temperature dependency of the volume CTE (α) will be α ∼ T³ (T → 0), and α ∼ const (α_Θ) at high temperature (T ≫ Θ), α_Θ is very close to the above-mentioned value of 2.37 × 10⁻⁴ K⁻¹ over the temperatures of 330–443 K. The densities estimated by the Debye expansion model are also listed in Table 9.

5 Conclusions

Presented is the all-atom, non-empirical, and tailor-made force field (NETMFF) for molecular modeling and simulation studies of α-RDX in the condensed phases. Its parameters are divided into three categories: valence, charge, and vdW terms. The vdW term adopts the damped Buckingham function that describes very well the vdW interactions whose physical nature is short-range forces and long-range dispersion forces. After definition of atom types was introduced, the charge parameters were first derived using a constrained ab initio ESP charge fit. Then with a set of initial vdW parameters derived from the SAPT(DFT) and WSM calculations, the vdW and valence parameters were derived based on ab initio data using a SCFF optimization technique. The final NETMFF parameters are required to accurately recover the density and intermolecular cohesive energy obtained from the periodic tailor-made DFT-D optimization, and thus implicitly containing the many-body nonadditive interactions in condensed phase.

The central part of this work is SCFF optimization procedure. The consistence of various parts of the force field except charge parameters is ensured by iterating the parameterization procedure. The procedure aims to avoid the “parameter correlation” trap as effectively as possible. The “parameter correlation” effect is due to that the force field functions are usually inadequate to represent some very complicated energy surfaces, and in some regions these functions are often redundantly used. The effect, coupled with incomplete sampling of the potential energy surfaces and unreasonable initial parameters, may result in an arbitrary ‘best-fit’ that is valid only in a very small region where the final parameters with unphysical values manage to cancel each other. When these parameters are used in real simulations, once coordinates fall even slightly outside the valid region, these simulations could completely fail. Therefore, another three critical issues of the work had to be addressed before performing the SCFF procedure:

(1) Using ab initio method, construct a desirable training set that acts as the basis of the charge and valence parameter optimization. It not only correctly reflect ab initio minimum energy states of interest, but also reasonably cover the ab initio potential energy surfaces outside the minimum energy states.

(2) Using wavefunction-based theory of intermolecular forces like the SAPT(DFT) and WSM method, derive accurate the pair-additive vdW parameters from the dimer's potential surface. The parameters act as initial input into the SCFF procedure, and however, their accuracy must be validated by the CSP method.

(3) Using the accurate periodic ab initio method, optimize the crystal structure and compute its density and intermolecular cohesive energy. The two properties act as the basis of the SCFF optimization and ensure that the final vdW parameters can correctly recover the non-additive interactions in condensed phase.

Based on the above-mentioned strategy, the obtained NETMFF for α-RDX can yield the correct geometries for each RDX conformers found by the DFT calculations in a good way, and also can identify the transition states obtained from the DFT calculations by force-field-based vibration analysis. However, a main defect of the force field was also found, that is, its ability in describing the out-of-plane bending angle of E type NO₂ group is obviously inferior to the ability in describing other geometrical types including the out-of-plane bending angle of A type NO₂ group. This causes the present force field describes intramolecular interactions only in a reasonable way rather than in a precise way. The most successful thing of the force field lies in that without any assumptions it can accurately predict densities of high explosives at the environmental conditions to which they are often subjected, a long-standing issue in the field of energetic materials. This ascribes to the accurate intermolecular parameters and reasonable valence parameters in the force field.

It also should be emphasized that the force field was developed only by means of first principles and starting with the crystal structure. Therefore, the parameterization strategy described in this paper can be easily generalized toward other known HEs or the new HEs whose crystal structures can be obtained by CSP.

Acknowledgements

The project was supported by the National Natural Science Foundation of China (11372053 and 11172044) and by the Opening Project of the State Key Laboratory of Explosion Science and Technology (KFJJ14-06M).

References

1. M. J. Kamlet and S. J. Jacobs, Chemistry of detonations. I. Simple method for calculating detonation properties of C–H–N–O explosives, J. Chem. Phys., 1968, 48, 23–35.
2. C. L. Mader, Organic Energetic Compounds, Nova Science Publishers, New York, 1996, p. 193.
3. C. M. Tarver, Density estimations for explosives and related compounds using the group additivity approach, J. Chem. Eng. Data, 1979, 24, 136–145.
4. J. R. Stine, Prediction of crystal densities of organic explosives by group additivity, Los Alamos National Laboratory's Report, New Mexico, 1981.
5. H. L. Ammon and S. Mitchell, A new atom/functional group volume additivity data base for the calculation of the crystal densities of C, H, N, O and F-containing compounds, Propellants, Explos., Pyrotech., 1998, 23, 260–265.
6. H. L. Ammon, New atom/functional group volume additivity data bases for the calculation of the crystal densities of C-, H-, N-, O-, F-, S-, P-, Cl-, and Br-containing compounds, Struct. Chem., 2001, 12, 205–212.
7. J. R. Holden, Z. Du and H. L. Ammon, Structure and density predictions for energetic materials, in Energetic Materials Part 1. Decomposition, Crystal and Molecular Properties, Elsevier, 2003, p. 188.
8. J. P. M. Lommerse, W. D. S. Motherwell and H. L. Ammon, et al., A test of crystal structure prediction of small organic molecules, Acta Crystallogr., Sect. B: Struct. Sci., 2000, 56, 697–714.
9. W. D. S. Motherwell, H. L. Ammon and J. D. Dunitz, et al., Crystal structure prediction of small organic molecules: a second blind test, Acta Crystallogr., Sect. B: Struct. Sci., 2002, 58, 647–661.
10. G. M. Day, W. D. S. Motherwell and H. L. Ammon, et al., A third blind test of crystal structure prediction, Acta Crystallogr., Sect. B: Struct. Sci., 2005, 61, 511–527.
11. G. M. Day, T. G. Cooper and A. J. Cruz-Cabeza, et al., Significant progress in predicting the crystal structures of small organic molecules – a report on the fourth blind test, Acta Crystallogr., Sect. B: Struct. Sci., 2009, 65, 107–125.
12. D. A. Bardwell, C. S. Adjiman and Y. A. Arnautova, et al., Towards crystal structure prediction of complex organic compounds – a report on the fifth blind test, Acta Crystallogr., Sect. B: Struct. Sci., 2011, 67, 535–551.
13. M. A. Neumann, F. J. J. Leusen and J. Kendrick, A Major Advance in Crystal Structure Prediction, Angew. Chem., Int. Ed., 2008, 47, 2427–2430.
14. M. A. Neumann, Tailor-Made Force Fields for Crystal-Structure Prediction, J. Phys. Chem. B, 2008, 112, 9810–9829.
15. R. Podeszwa, B. M. Rice and K. Szalewicz, Predicting Structure of Molecular Crystals from First Principles, Phys. Rev. Lett., 2008, 101, 115503.
16. Y.-Q. Wu and F.-L. Huang, Thermal mechanical anisotropic constitutive model and numerical simulations for shocked b-HMX single crystals, Eur. J. Mech. Solid., 2012, 36, 66–82.
17. D. C. Sorescu, B. M. Rice and D. L. Thompson, Intermolecular potential for the hexahydro-1,3,5-trinitro-1,3,5-s-triazine crystal (RDX): A crystal packing, Monte Carlo, and molecular dynamics study, J. Phys. Chem. B, 1997, 101, 798–808.
18. G. D. Smith and R. K. Bharadwaj, Quantum Chemistry Based Force Field for Simulations of HMX, J. Phys. Chem. B, 1999, 103, 3570–3575.
19. 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.
20. E. P. Wallis and D. L. Thompson, Molecular dynamics simulations of ring inversion in RDX, J. Chem. Phys., 1993, 99, 2661–2673.
21. S. Boyd, M. Gravelle and P. Politzer, Nonreactive molecular dynamics force field for crystalline hexahydro-1,3,5-trinitro-1,3,5 triazine, J. Chem. Phys., 2006, 124, 104508.
22. S. W. Bunte and H. Sun, Molecular Modeling of Energetic Materials: The Parameterization and Validation of Nitrate Esters in the COMPASS Force Field, J. Phys. Chem. B, 2000, 104, 2477–2489.
23. R. Podeszwa, R. Bukowski, B. M. Rice and K. Szalewicz, Potential energy surface for cyclotrimethylene trinitramine dimer from symmetry-adapted perturbation theory, Phys. Chem. Chem. Phys., 2007, 9, 5561–5569.
24. H. L. Williams and C. F. Chabalowski, Using Kohn–Sham orbitals in symmetry-adapted perturbation theory to investigate intermolecular interactions, J. Phys. Chem. A, 2001, 105, 646–659.
25. A. J. Misquitta and K. Szalewicz, Intermolecular forces from asymptotically corrected density functional description of monomers, Chem. Phys. Lett., 2002, 357, 301–306.
26. A. J. Misquitta, B. Jeziorski and K. Szalewicz, Dispersion energy from density-functional theory description of monomers, Phys. Rev. Lett., 2003, 91, 033201.
27. A. J. Misquitta and K. Szalewicz, Symmetry-adapted perturbation-theory calculations of intermolecular forces employing density-functional description of monomers, J. Chem. Phys., 2005, 122, 214109.
28. A. J. Misquitta, R. Podeszwa, B. Jeziorski and K. Szalewicz, Intermolecular potentials based on symmetry-adapted perturbation theory with dispersion energies from time-dependent density-functional calculations, J. Chem. Phys., 2005, 123, 214103.
29. A. Hesselmann and G. Jansen, First-order intermolecular interaction energies from Kohn–Sham orbitals, Chem. Phys. Lett., 2002, 357, 464–470.
30. H. Song, H. Xiao and H. Dong, Correlated intermolecular interaction components from asymptotically corrected Kohn–Sham orbitals, Sci. China, Ser. B: Chem., 2004, 47, 466–479.
31. J. G. Williams and A. J. Stone, Distributed dispersion: A new approach, J. Chem. Phys., 2003, 119, 4620–4628.
32. A. J. Misquitta and A. J. Stone, Distributed polarizabilities obtained using a constrained density-fitting algorithm, J. Chem. Phys., 2006, 124, 024111.
33. A. J. Misquitta and A. J. Stone, Accurate Induction Energies for Small Organic Molecules: 1. Theory, J. Chem. Theory Comput., 2008, 4, 7–18.
34. A. J. Misquitta, A. J. Stone and S. L. Price, Accurate Induction Energies for Small Organic Molecules. 2. Development and Testing of Distributed Polarizability Models against SAPT(DFT) Energies, J. Chem. Theory Comput., 2008, 4, 19–32.
35. A. J. Misquitta and A. J. Stone, Dispersion energies for small organic molecules: first row atoms, Mol. Phys., 2008, 106, 1631–1643.
36. H.-J. Song and F.-L. Huang, Accurately Predicting the Density and Hydrostatic Compression of Hexahydro-1,3,5-Trinitro-1,3,5-Triazine from First Principles, Chin. Phys. Lett., 2011, 28, 096103.
37. H.-J. Song and F. Huang, Accurately predicting the structure, density, and hydrostatic compression of crystalline β-1,3,5,7-tetranitro-1,3,5,7-tetraazacyclooctane based on its wave-function–based potential, EPL, 2011, 95, 53001.
38. A. R. Leach, Molecular Modeling Principles and Applications, Addison Wesley Longman, London, 1996, p. 145.
39. A. J. Stone, Intermolecular Potentials, Science, 2008, 321, 787–789.
40. N. Foloppe and A. D. Mackerell, All-Atom Empirical Force Field for Nucleic Acids: I. Parameter Optimization Based on Small Molecule and Condensed Phase Macromolecular Target Data, J. Comput. Chem., 2000, 21, 86–104.
41. S. Kirkpatrick, C. D. Gelatt and M. P. Vecchi, Optimization by Simulated Annealing, Science, 1983, 220, 671–680.
42. T. Vladimiroff and B. M. Rice, Reinvestigation of the Gas-Phase Structure of RDX Using Density Functional Theory Predictions of Electron-Scattering Intensities, J. Phys. Chem. A, 2002, 106, 10437–10443.
43. I. M. Sobol, On the distribution of Points in a Cube and The Approximate Evaluation of Integrals, USSR Comput. Math. Math. Phys., 1967, 7, 86–112.
44. C. M. Breneman and K. B. Wiberg, Determining atom-centered monopoles from molecular electrostatic potentials – the need for high sampling density in formamide conformational-analysis, J. Comput. Chem., 1990, 11, 361–373.
45. Dalton, a molecular electronic structure program, Release 2.0, 2005, see http://www.kjemi.uio.no/software/dalton/dalton.html.
46. A. J. Misquitta, G. W. A. Welch, A. J. Stone and S. L. Price, A first principles prediction of the crystal structure of C₆Br₂ClFH₂, Chem. Phys. Lett., 2008, 456, 105–109.
47. T. S. Totton, A. J. Misquitta and M. Kraft, A First Principles Development of a General Anisotropic Potential for Polycyclic Aromatic Hydrocarbons, J. Chem. Theory Comput., 2010, 6, 683–695.
48. C. S. Choi and E. Prince, The Crystal Structure of Cyelotrimethylene-trinitramine, Acta Crystallogr., Sect. B: Struct. Crystallogr. Cryst. Chem., 1972, 28, 2857–2862.
49. A. Bondi, van der Waals Volumes and Radii, J. Phys. Chem., 1964, 68, 441–451.
50. K. T. Tang and J. P. Toennies, J. Chem. Phys., 1984, 80, 3726–3741.
51. A. J. Misquitta and A. J. Stone, CamCASP: a program for studying intermolecular interactions and for the calculation of molecular properties in distributed form, University of Cambridge, 2009.
52. A. J. Stone, A. Dullweber, M. P. Hodges, P. L. A. Popelier and D. J. Wales, ORIENT 3.2, University of Cambridge, 1996.
53. A. J. Stone, The Theory of Intermolecular Forces, Clarendon Press, Oxford, 2000.
54. S. Grimme, Semiempirical GGA-Type Density Functional Constructed with a Long-Range Dispersion Correction, J. Comput. Chem., 2006, 27, 1787–1799.
55. P. Jurecka, J. Cerny, P. Hobza and P. D. R. Salahub, Density Functional Theory Augmented with an Empirical Dispersion Term. Interaction Energies and Geometries of 80 Noncovalent Complexes Compared with Ab Initio Quantum Mechanics Calculations, J. Comput. Chem., 2007, 28, 555–569.
56. F. Ortmann, F. Bechstedt and W. G. Schmidt, Semiempirical van der Waals correction to the density functional description of solids and molecular structures, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 205101.
57. A. Tkatchenko and M. Scheffler, Accurate Molecular van der Waals Interactions from Ground-State Electron Density and Free-Atom Reference Data, Phys. Rev. Lett., 2009, 102, 073005.
58. B. Civalleri, C. M. Zicovich-Wilson, L. Valenzano and P. Ugliengo, B3LYP augmented with an empirical dispersion term (B3LYP-D*) as applied to molecular crystals, CrystEngComm, 2008, 10, 405–410.
59. R. Dovesi, V. R. Saunders, C. Roetti, R. Orlando, C. M. Zicovich-Wilson, F. Pascale, B. Civalleri, K. Doll, N. M. Harrison, I. J. Bush, P. D'Arco and M. Llunell, CRYSTAL09 User's Manual, University of Torino, Torino, 2009.
60. R. B. Cundall, T. F. Palmer and C. E. C. Wood, Vapour Pressure Measurements on Some Organic High Explosives, J. Chem. Soc., Faraday Trans. 1, 1978, 74, 1339–1345.
61. A. K. Rappe, C. J. Casewit, K. S. Colwell, W. A. Goddard III and W. M. Skiff, UFF, a full periodic table force field for molecular mechanics and molecular dynamics simulations, J. Am. Chem. Soc., 1992, 114, 10024–10035.
62. B. Olinger, B. Roof and H. H. Cady, Symposium (Int'l) on High Dynamic Pressures, C. E. A., Paris, France, 1978, pp. 3–8.
63. F. D. Murnaghan, Finite deformations of an elastic solid, Am. J. Math., 1937, 49, 235–260.
64. J. Sun, X. Shu, Y. Liu, H. Zhang, X. Liu, Y. Jiang, B. Kang, C. Xue and G. Song, Investigation on the Thermal Expansion and Theoretical Density of 1,3,5-Trinitro-1,3,5-Triazacyclohexane, Propellants, Explos., Pyrotech., 2010, 35, 1–6.
65. P. Hakey, W. Ouellette, J. Zubieta and T. Korter, Redetermination of cyclotrimethylenetrinitramine, Acta Crystallogr., Sect. E: Struct. Rep. Online, 2008, 64, o1428.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c4ra07195f

---