An Empirically Optimized Classical Force-Field for Molecular

Article pubs.acs.org/JPCA

An Empirically Optimized Classical Force-Field for Molecular Simulations of 2,4,6-Trinitrotoluene (TNT) and 2,4-Dinitrotoluene (DNT) S. Neyertz,*,† D. Mathieu,‡ S. Khanniche,†,‡ and D. Brown† †

LEPMI(LMOPS), UMR 5279 CNRS-University of Savoie-Grenoble INP-University J. Fourier, Bât. IUT, Savoie Technolac, 73376 Le Bourget du Lac Cedex, France ‡ CEA, DAM, Le Ripault, 37260 Monts, France ABSTRACT: An empirical classical all-atom specific force-field for use in molecular dynamics simulations (MD) has been developed to reproduce the experimental densities and structures of trinitrotoluene (TNT) in its crystalline and liquid phases at six different temperatures, as well as its enthalpies of sublimation and fusion. The average structural parameters and partial charges were obtained from density functional theory optimizations of single molecules at the B3LYP/6-311+G** level. The other constants for the potential were adjusted in order to obtain a classical force-field, which is able to reproduce the aforementioned properties for TNT with a high degree of accuracy. This force-field was also found to predict closely the experimental densities and structures of 2,4-dinitrotoluene (2,4-DNT) in its crystalline and liquid phases as well as its enthalpy of sublimation. It was a bit less successful for its enthalpy of fusion, but it still remained reasonable, and the model mechanical properties were of the right order of magnitude. As such, this fairly simple force-field can be used for MD simulations of both TNT and 2,4-DNT nitroaromatic compounds. are generic.15,24−27 Other works include Monte Carlo simulations of solubilities in supercritical CO2 with TNT united-atom models (either a monosite or a seven-site model),28 or the transfer of an intermolecular potential developed for cyclotrimethylene trinitramine (RDX) to crystalline TNT while keeping the molecules rigid.29 In the present article, we develop a classical all-atom specific potential for TNT, following a commonly used empirical approach based on quantum chemical calculations, literature values and experimental data.30 The average structural parameters and partial charges obtained from ab initio optimizations of the NAC single-molecules and values issued from the literature are used as starting points to parametrize this classical force-field, which is able to reproduce the experimental densities of TNT both in the crystal and liquid phases, the experimental structures of crystalline TNT at various temperatures, as well as its enthalpies of sublimation and fusion. We show that this relatively simple force-field can be applied with only a slight change to 2,4-DNT as well. The starting configurations for the model test systems are described in section 2, the empirically optimized force-field is given in section 3, and the results of MD simulations of TNT and 2,4DNT in their different phases are compared to experimental evidence in section 4.

1. INTRODUCTION 2,4,6-Trinitrotoluene (TNT), C7N3O6H5, is well-known as an explosive material. For health, security, and environmental reasons, there is a need for gas sensors that are sensitive to very low concentrations of such nitroaromatics compounds (NAC).1 The detection device also has to be both selective and robust, and several strategies for the development of NAC or related molecules gas sensors have been proposed.2−4 Some of them are based on, e.g., functionalized mesoporous films coated on a quartz crystal microbalance.5−7 However, it is not straightforward to assess what really happens at the atomic level when the nitroaromatic molecules meet the sensor interface and how this can be related to specific sensor chemical structures and morphologies. Such information is important for designing more selective and more efficient sensors. In addition, the NAC−sensor interactions have to be weak enough for the adsorptions to be reversible, which allows for the sensor to be reused again.4 Computer simulation techniques such as fully atomistic classical molecular dynamics (MD)8 can, in principle, be used to complement experimental evidence. The prerequisites of such simulations are realistic starting configurations as well as appropriate force-fields for both the sensors and the NACs. For the latter, the literature is rather scarce on the subject of specific nitroaromatic classical force-fields, which are able to accurately reproduce their physical properties. TNT has been extensively studied using ab initio calculations with a large variety of basis sets.9−20 This has also been the case, albeit to a lesser extent, for 2,4-DNT,21−23 which is a degradation product of TNT.12 However, most classical force-fields applied to nitroaromatics © 2012 American Chemical Society

Received: June 1, 2012 Revised: July 23, 2012 Published: July 23, 2012 8374

dx.doi.org/10.1021/jp305362n | J. Phys. Chem. A 2012, 116, 8374−8381

The Journal of Physical Chemistry A

Article

Table 1. Optimized Force-Field Parameters Related to eqs 1−6; They Apply to Both TNT and 2,4-DNT except for the Car− Car−Northo−O Torsion Explicitely Defined As (T) for TNT and (D) for DNT bond types

b0 (Å)

Car−Car Car−Cmet Car−N Car−H Cmet−Hmet N−O Dihedral type

1.393 Cmet−Car−Car 1.506 Car−Car−Car 1.483 Car−Car−N 1.081 Car−Car−H 1.090 Car−N−O 1.222 O−N−O a0 (kJ mol−1) a1 (kJ mol−1)

*−Car−Car−* Hmet−Cmet−Car−Car Car−Car−Northo−O (T) Car−Car−Northo−O (D) Car−Car−Npara−O

bond-angle type

16.74 0.50 8.28 11.97 13.39

0 −1.51 0 0 0

θ0 (deg)

kθ (kJ mol−1)

out-of-plane atom

koop (kg s−2)

122.47 120.00 118.60 120.54 117.30 125.39 a2 (kJ mol−1)

879.1 879.1 2270.9 879.1 879.1 732.5 a3 (kJ mol−1)

Car N

667 417

−16.74 0 −8.28 −11.97 −13.39

0 2.01 0 0 0

2. TEST SYSTEMS All ab initio calculations were performed with the Gaussian09 code,31 while the classical simulations were carried out using the molecular dynamics codes of the gmq package.32 The TNT and 2,4-DNT molecules were first geometryoptimized using density functional theory (DFT) at the B3LYP/6-311+G** level. The average equilibrium bond lengths and bond angles were taken from the optimized single-molecule structures. The partial charges qi/e were calculated using both the CHELPG33 and the MK34 methods. They were extracted by an ESP-fitting procedure35 and symmetrized. Although attempts to solve the crystal structures of TNT have been made for quite a long time,36 they were hindered by the polymorphic behavior of this molecule. Several crystal structures have now been proposed, and there is a general consensus that TNT can crystallize in both a monoclinic and an orthorhombic form.37−39 According to Vrcelj et al.,39 the most stable form is the monoclinic crystal with cell parameters of a = 14.9113 Å, b = 6.0340 Å, c = 20.8815 Å, β = 110.365°, and the P21/a space group at the temperature T = 100 K. There is also a less stable orthorhombic form at T = 123 K, with cell parameters of a = 14.910 Å, b = 6.034 Å, c = 19.680 Å, and the Pca21 space group. The phase transformation from the orthorhombic to the monoclinic form occurs a few degrees below the melting point situated at ∼352−354 K,40,41 as shown by endotherms obtained by differential scanning calorimetry.39 The corresponding densities are 1713 kg m−3 for the monoclinic and 1704 kg m−3 for the orthorhombic form. Although to our knowledge, no unit-cell parameters have been reported at room temperature, the density of TNT is given as being ∼1654 kg m−3 at 298 K.42 Crystalline simulation boxes were built from the experimental data of Vrcelj et al.39 by scaling the unit-cell by 3 × 7 × 2 (7056 atoms and 336 molecules) for the monoclinic structure at 100 K and by 2 × 3 × 2 (2016 atoms and 96 molecules) for the orthorhombic structure at 123 K. The monoclinic TNT simulations were carried out at 100, 298, and 352 K, while the orthorhombic TNT simulations were carried out at both 123 and 298 K. 2,4-DNT has also been reported as crystallizing in a monoclinic form. According to Hanson et al.,43 the monoclinic crystal has cell parameters of a = 8.0057 Å, b = 15.1273 Å, c = 12.8853 Å, β = 95.877°, a density of 1552 kg m−3, and the P21/ n (No. 14) space group at T = 173 K. There is another monoclinic 2,4-DNT crystal with the same space group but at T

like-atom pair i···i

σii (Å)

εii /kB (K)

Car···Car Cmet···Cmet N···N O···O H···H Hmet···Hmet

3.029 3.029 2.762 2.950 2.673 2.673

53.8 53.8 47.8 70.0 21.1 40.0

= 298 K, which has cell parameters of a = 8.126 Å, b = 15.217 Å, c = 12.998 Å, β = 94.77°, and a density of 1511 kg m−3.44 Crystalline simulation boxes were built using the aforementioned experimental data at both 173 and 298 K by scaling the respective unit-cells by 3 × 2 × 2 (1824 atoms and 96 molecules). In addition, with the melting point of 2,4-DNT being at ∼341−343 K,40,41 an additional crystalline simulation was carried out at 343 K. Liquid simulation boxes were constructed by placing 64 TNT molecules at a density close to the experimental value of 1443 kg m−3 at 373 K.45 From the same source, the liquid density at the melting point 352 K and at 423 K can be estimated as being 1464 and 1391 kg m−3, respectively. Liquid TNT simulations were thus carried out at 352, 373, and 423 K. In the case of 2,4-DNT, 64 molecules were placed at the experimentally estimated density of 1327 kg m−3, corresponding to the melting temperature of 343 K.45 In all cases, the simulation boxes were subsequently allowed to relax toward their equilibrium densities.

3. EMPIRICALLY-OPTIMIZED MOLECULAR MODEL The force-field used in this study describes the potential energy of the system as a superposition of simple analytical functions (eq 1): Upot =

∑ Ubend(θ) + ∑ Utors(τ) + ∑ Uoop(i) θ

+

i−sp2

τ

∑ (i , j)nb

Uvdw(r ) +

∑ (i , j)nb

Ucoul(r ) (1)

In eq 1, the first three terms arise from near-neighbor connections in the structure and can be referred to as the bonded potentials. The last two terms, referred to as the nonbonded potentials, depend on the distance r between two interacting sites, which are situated either on the same molecule (but are separated by more than two bonds) or on two different molecules. All parameters related to eq 1 are given in Table 1, with the exception of the partial charges, which are shown in Figure 1. Each specific term in eq 1 will be separately explained in the following description. The acronyms used in Table 1 for the different atom types in TNT and 2,4-DNT are Car for the aromatic carbon, Cmet for the methyl carbon, O for the oxygen, N for the nitrogens that can be considered separately depending on their position in the 8375



Article

usual approach found in generic force-fields, i.e., the description of the torsional energy via a combination of torsional and nonbonded potentials. The sp2 ring and NO2 structures are kept planar by using a harmonic function in the perpendicular distance d from the isp2 atom to the plane defined by its three attached atoms: Uoop(i−sp2 ) = Figure 1. Optimized partial charges qi/e for (A) TNT and (B) 2,4DNT.

(4)

(5)

where ε is the well-depth of the potential and σ is the distance at which the potential is zero. The parameters were first taken from TRIPOS 5.2 and empirically adjusted to maintain the crystalline structures of TNT. They are given in Table 1 for like-atom pairs, while the cross-terms for unlike-atom pairs can be obtained from the standard combination rules:8 σij =

σii + σjj 2

εij =

(εii × εjj)

(6)

The last term of eq 1 accounts for the Coulombic interactions with qi and qj being the charges on atoms i and j, respectively, and ε0 being the vacuum permittivity qiqj Ucoul(r ) = 4πε0r (7) The electrostatic potential is calculated using the Ewald summation method because of its long-range character.49,50 Trials were initially carried out with both the CHELPG and the MK partial charges. It was found that the CHELPG partial charges scaled by a factor of 0.95 for both TNT and 2,4-DNT best reproduced the experimental densities and enthalpies. Both sets of scaled charges are given in Figure 1. For each system under study, the MD simulations were carried out using the starting structures described in section 2 and the optimized force-field (eq 1, Table 1, and Figure 1) for ∼10 ps under constant volume NVT conditions (controlled number of atoms N, controlled volume V, and controlled temperature T). The liquid models were then switched to NpT conditions, in which the isotropic pressure p was maintained at 1 bar by loose-coupling with a constant equal to 10 ps,51 and the simulation box was kept cubic. The crystalline models were run under NPT conditions in which the on-diagonal and offdiagonal components of the pressure tensor P were maintained at 1 and 0 bar, respectively, by loose-coupling, but the simulation box was allowed to relax toward its equilibrium shape and density. The crystalline models at 298, 343, and 352 K for TNT and 343 K for 2,4-DNT were obtained by heating up the experimental crystal structures at a rate of 1 K ps−1, and the temperature was maintained at the desired value by loosecoupling to a heat bath52 with a constant equal to 0.1 ps. The truncation radii used for both the real-space electrostatic and van der Waals contributions were set to 10 Å, and an optimal convergence of the Ewald sums53 was obtained for (α, Kmax) parameters equal to (0.32, 14) for the monoclinic TNT crystal,

3 n=0

d2

⎛⎛ σ ⎞12 ⎛ σ ⎞6 ⎞ Uvdw(r ) = ULJ(r ) = 4ε⎜⎜ ⎟ − ⎜ ⎟ ⎟ ⎝r⎠ ⎠ ⎝⎝ r ⎠

k Ubend(θ ) = θ (cos θ − cos θ0)2 (2) 2 where kθ is a constant determining the flexibility of the angle, and θ0 is the equilibrium bond angle. The latter were also extracted from the average geometries optimized using Gaussian09 at the B3LYP/6-311+G** level. The kθ were taken from the freely available TRIPOS 5.2 force-field48 and converted to the analytical form of eq 2 for all angles, except for those related to the methyl group. Indeed, the CH3 angle constraints fixed Hmet−Cmet−Hmet to its equilibrium θ0 value of 108.06° and Hmet−Cmet−Car to θ0 = 110.85°. The other θ0 and their associated kθ can be found in Table 1. The torsional motions around the dihedral angles τ are represented by a polynomial in cos τ:

∑ an cosn τ

2

The koop constants (Table 1) were converted from TRIPOS 5.2 and also adjusted when necessary to reproduce the Gaussian09optimized structures. The nonbonded van der Waals interactions are described by the Lennard-Jones form:

molecule as Northo and Npara, H for the aromatic hydrogens, and Hmet for the methyl hydrogens. High-frequency intramolecular bond stretching modes are removed by rigidly constraining all bond lengths, while the methyl groups are kept in a tetrahedral geometry according to the SHAKE routine.46,47 It allows for the use of a reasonable time step in the MD integration algorithm, i.e., Δt = 10−15 s, and avoids problems of equipartition of kinetic energy. The relative tolerance for all geometrical constraints is set to 10−7. The bond lengths b0 (Table 1) have been obtained from the average ab initio optimized TNT and 2,4-DNT geometries. In eq 1, the angle-bending deformations are described by a harmonic function in the cosine of the bond angles θ:

Utors(τ ) =

koop

(3)

with the dihedral angle τ varying from −180° to +180° and τ = 0° being the trans conformation. The an coefficients were again derived from TRIPOS 5.2, with the exception of the ortho Car− Car−N−O coefficients, which were adjusted in order to reproduce the actual Gaussian09-optimized torsions in the gas phase. The an torsional coefficients used in eq 3 are given in Table 1. The wild card denotes any atom-type, except when the coefficients for a specific dihedral are defined elsewhere in the table. It should be noted that, since two sites that are separated by more than two bonds also interact via nonbonded potentials, the total torsional energy for a 1−2−3−4 atom sequence is actually described by the torsional potential around bond 2−3 plus the 1···4 nonbonded interactions. Another parametrization approach would have been to describe the total torsional energy by eq 3 on its own by using directly ab initio calculated torsional profiles. Although this was attempted, there were such strong 1···6 ortho O···Hmet interactions that the ab initiocalculated torsion potential was not enough to keep the gasphase optimized torsions on its own. As such, we adopted the 8376



Article

Figure 2. Schematic representations of the TNT simulation cells at the end of the MD production runs. (A) Monoclinic crystal at 100 K, (B) orthorhombic crystal at 123 K, and (C) liquid phase at 373 K. In panels A and B, the different atom-types are displayed with the following color code: red = O; cyan = C; blue = N; white = H. In panel C, each TNT molecule is shown in a different color in order to visualize the disorder characteristic of the liquid phase.

the one used here.11 In the DFT-optimized 2,4-DNT structure, the ortho nitro-group is also twisted with respect to the aromatic cycle, but only by ∼30°, in agreement with a calculation carried out with the same basis set and an earlier version of Gaussian.22 As noted for TNT, lower basis sets tend to underestimate the twisting angle by ∼10°.21 In order for the classical force-field to reproduce correctly both the DFToptimized structures and the specificities of the ortho groups, specific torsional parameters for the Car−Car−Northo−O dihedral had to be defined depending on whether the molecule is TNT or 2,4-DNT (see Table 1). Along with the partial charges, this is the only difference between both types of molecular models since the para-nitro group was always found in the same plane as the aromatic ring. The minimum intramolecular distance between an orthonitro O and a methyl H in both TNT and 2,4-DNT is ∼2.55 Å, which qualifies it for an H-bond interaction.55 This has already been pointed out in other ab initio characterizations,22 and it is in agreement with the proposed mechanisms of decomposition involving the rupture of the C−NO2 bond weakened by hydrogen bonding with one methyl H.13,56 The average experimental ⟨ρexp⟩ and model ⟨ρmod⟩ densities for the TNT and 2,4-DNT systems under study are given in Table 2, along with the average model intermolecular energies ⟨Uinter pot ⟩. It is clear from Table 2 that the empirically optimized potential reproduces rather well the crystalline and liquid

(0.34, 12) for the orthorhombic TNT crystal, the monoclinic 2,4-DNT crystal, and liquid 2,4-DNT, as well as (0.32, 8) for liquid TNT. Long-range corrections to the van der Waals contributions to the energy and the pressure were added using the approximation that the radial distribution functions are equal to unity beyond the cutoff.8 Equilibration runs were carried out for 500 ps in order to allow for the densities to stabilize after which the production run continued up to 1500 ps for the crystalline simulations and up to 2000 ps for the liquid simulations. Configurations were stored every 10 ps and thermodynamic properties every 1 ps for postanalyses. Schematic representations of TNT in two crystalline and the liquid phases are displayed in Figure 2 using the VMD 1.8.2 visualization software.54

4. PREDICTED PROPERTIES The optimized TNT structure calculated with density functional theory at the B3LYP/6-311+G** level shows a planar aromatic ring symmetry with the planes of the two equivalent ortho-nitro groups being rotated with respect to the aromatic ring by ∼40°. It is well-known that the level of theory used in the optimization affects the molecular geometries. Low basis sets such as 6-31G* tend to give values around 30°,10 but larger basis sets such as 6-311+(d,p) predict values around 40°.10,11 A study of TNT with several ab initio methods and basis sets actually showed that the most stable conformation is indeed obtained with the same DFT/B3LYP/6-311+G** approach as 8377



Article

Table 2. Average Experimental ⟨ρexp⟩ and Model ⟨ρmod⟩ Densities in kg m−3 for TNT and 2,4-DNT at Different Temperatures T in Ka phase

T

⟨ρexp⟩

⟨ρmod⟩

⟨Uinter pot ⟩

monoclinic TNT monoclinic TNT monoclinic TNT orthorhombic TNT orthorhombic TNT liquid TNT liquid TNT liquid TNT monoclinic DNT monoclinic DNT monoclinic DNT liquid DNT

100 298 352 123 298 352 373 423 173 298 343 343

1713 1654

1720.7 ± 0.1 1654.1 ± 0.2 1632.6 ± 0.2 1715.1 ± 0.1 1654.7 ± 0.3 1452 ± 2 1439 ± 2 1386 ± 1 1542.7 ± 0.2 1497.3 ± 0.3 1477.9 ± 0.3 1334 ± 2

−131.9 ± 0.1 −121.9 ± 0.1 −118.7 ± 0.1 −130.8 ± 0.1 −121.8 ± 0.1 −95.1 ± 0.2 −94.0 ± 0.2 −87.5 ± 0.2 −104.8 ± 0.1 −99.0 ± 0.1 −96.5 ± 0.1 −81.2 ± 0.2

1704 1654 1464 1443 1391 1552 1511 1327

DFT-optimized angles because of the crystal packing interactions.10,15 The corresponding crystalline MD angle probability density distributions are much smoother because of the thermal vibrations, but Figure 3 shows that the model results agree very well with the static peaks obtained from the refined X-ray data.39 Please note that the angles referring to the ortho-nitro groups are displayed with full symbols, while those referring to the para-nitro groups are in open symbols. Similar to the angles that are smeared out by the effect of motion in the MD runs, the crystalline TNT simulation boxes are found to relax slightly from the static crystalline cells.39 While most of these deformations amount to less than 0.6%, the MD-averaged b and c axes tend to vary a bit more (∼3− 4%), albeit in opposite directions, while the β angle of the monoclinic structure increases by ∼1.8%. Such variations are not unexpected since they are dependent on the complex interplay of potential functions, and indeed much larger differences between model and experimental crystals (5−6%) have been considered as satisfactory in the literature.61 In addition, these relaxations of the crystalline boxes only lead to very subtle changes in the actual structures as shown by the pseudo X-ray powder diffractograms (Debye formula), which were generated from the coordinates of the static39 and the MD-averaged TNT crystalline structures using the DISCUS program62 with a wavelength of λ = 1.5406 Å. The corresponding intensity vs 2θ plots are given in Figure 4 and confirm that the MD-averaged structures remain in excellent agreement with experimental data. As shown in Table 2, the empirically optimized force field is also found to describe very well the liquid phase of TNT (Figure 2C). The differences between ⟨ρmod⟩ and ⟨ρexp⟩ are 0.8% at 352 K, 0.3% at 373 K, and 0.4% at 423 K. While we are not aware of any experimental data about the vaporisation enthalpy of TNT, its enthalpy of fusion ΔHfus at the melting point temperatures of 352−354 K has been reported as being 21.2 kJ mol−140 and 23.4 kJ mol−1.41 The difference in model enthalpies at 352 K between the liquid and the monoclinic crystal is ΔHfus = 20.1 ± 0.4 kJ mol−1, which is clearly very close to the former value. For crystalline 2,4-DNT, the differences in densities between model and experimental structures are ∼0.6% and ∼0.9% for the monoclinic phases at 173 and 298 K, respectively. The experimental enthalpy of sublimation for 2,4-DNT of 99.6 ± 1.3 kJ mol−1 at 298 K63 is also in excellent agreement with our

The model intermolecular potential energies ⟨Uinter pot ⟩ are given in kJ mol−1. References to the experimental densities are provided in section 2. All model results are displayed with their associated standard errors. a

densities of both TNT and 2,4-DNT. For crystalline TNT, the differences in densities between model and experiment are only ∼0.4% for the monoclinic phase at 100 K (Figure 2A) and ∼0.7% for the orthorhombic phase at 123 K (Figure 2B). This agreement is even better at 298 K, where ⟨ρmod⟩ is within 0.1% from ⟨ρexp⟩. As can be seen in the average intermolecular potential energies ⟨Uinter pot ⟩, the experimental observation of the monoclinic form being more stable than the orthorhombic form, but with only a small energetic difference between both crystalline forms,39 is also seen in the models. In the literature, the experimental enthalpy of sublimation for TNT is reported as being in the range ∼100−125 kJ mol−1 depending on the sources.57−60 Our model enthalpy of sublimation can be approximated from the solid intermolecular potential energy by −1 ΔHsub ≈ RT − ⟨Uinter at 352 pot ⟩, and our ΔHsub of 121.6 kJ mol K is in agreement with the value of 118.4 ± 4.2 kJ mol−1 given by Edwards for temperatures of 323 to 353 K.57 In the experimental X-ray diffraction studies, several authors37−39 report two types of molecules in the asymmetric unit that differ very slightly (less than 10°) in the angles between the planes of the nitro groups and the aromatic rings. As such, the experimentally refined static structures display a spread of Car−Car−N−O angles, which are different from the

Figure 3. Probability densities for the TNT Car−Car−N−O dihedral angles in both (A) the monoclinic crystal at 100 K and (B) the orthorhombic crystal at 123 K. The experimental X-ray refined static angles39 are shown with symbols and lines, while the MD-averaged angles are displayed with symbols only. They have been further resolved into ortho-NO2 (full symbols) and para-NO2 (open symbols). 8378



Article

Figure 4. Model pseudo X-ray diffractograms for both (A) the TNT monoclinic crystal at 100 K and (B) the TNT orthorhombic crystal at 123 K. The experimental data are displayed with lines and the MD-averaged data with symbols.

Figure 5. (A) Experimental (symbols and lines) and MD-averaged probability densities (symbols only) for the 2,4-DNT Car−Car−N−O dihedral angles in the monoclinic crystal at 298 K. (B) Pseudo X-ray diffractograms for the 2,4-DNT monoclinic crystal at 173 K. The experimental data are displayed with lines and the MD-averaged data with symbols.

TNT properties. As such, it was decided to keep the relatively simple force-field as similar as possible for both NACs under study. Mechanical properties were not used directly to develop the force-field, as there are, to our knowledge, no related experimental data available on single-crystals of either TNT or 2,4-DNT. However, the model crystalline tensile (Young’s) moduli E, bulk moduli K, and shear moduli G at room temperature were obtained by ramping up or down the appropriate required pressure tensor P components at a typical (absolute) rate of 10 bar ps−1 and measuring the associated changes in strain at low deformations. Although the induced strain rates are much higher than in experiment, due to the limited time scale available to MD simulations, checks were made to ensure that there was no significant rate effect by fixing pressure tensor components at their final ramped-up or -down values and subsequently monitoring the strain. No creep behavior was detected in these relatively fast relaxing small molecule systems. The results for the moduli (Table 3) can thus be considered as being equivalent to those obtained at limiting low deformation rates. The moduli are typical of those for pure organic crystals and indeed the values compare well with the mechanical properties of the pure crystalline phase of another explosive, bicyclo-HMX (cis-1,3,4,6-tetranitrooctahydroimidazo-[4,5-d]imidazole), which were similarly obtained from MD simulations.64 The Poisson’s ratio are also typically of the order of 0.3−0.4.64 For the tensile and shear moduli, there are variations due to the anisotropy of the systems, but in general, the values are of the same order of magnitude for both monoclinic and orthorhombic TNT. However, the mechanical properties of 2,4-

model approximated sublimation enthalpy at standard conditions, ΔHsub ≈ 101.5 ± 0.1 kJ mol−1. As for TNT (Figures 3 and 4), the MD angle probability density distributions of the nitro Car−Car−N−O angles are smoothed by thermal motion but remain close to those of the static crystal structures. This is displayed for T = 298 K by Figure 5A. The box length deformations are less than ∼3%, although it should be noted that the β angle of both monoclinic structures increases from ∼95° to ∼100°. The empirical force-field is thus a bit less accurate for 2,4-DNT than for TNT. However, as shown at 173 K by Figure 5B, the X-ray pseudodiffractograms intensity vs 2θ plots confirm that the MD-averaged structures also remain crystalline in the case of 2,4-DNT. As far as the liquid phase of 2,4-DNT is concerned, the model liquid density is ∼0.5% from the experimentally extrapolated value,45 which is very close. At the melting point temperature of 343 K, the enthalpy of fusion ΔHfus for 2,4DNT has been reported as being 20.1 kJ mol−1.40 The difference in enthalpy between our model at 343 K for the liquid and for the monoclinic crystal is ΔHfus = 15.4 ± 0.2 kJ mol−1, which differs more from the experimental value than the other properties. However, it should be noted that there can be large uncertainties in the experimental measurements of such thermodynamic data in nitroaromatics.63 In addition, considering the fact that the classical force-field was primarily optimized with respect to the TNT experimental data, its agreement for 2,4-DNT in terms of densities and enthalpy of sublimation is still very good. A few attempts were carried out in order to increase the model approximated enthalpy of fusion for 2,4-DNT, but they involved changing other parameters in Table 1, which in turn led to poorer representations of the 8379



Article

Table 3. Tensile E, Bulk K, and Shear G Moduli in GPa for the Model Crystalline Structures of TNT and 2,4-DNT at 298 K; Standard Errors on the Values Are of the Order of 0.1 GPa phase P tensor component changed monoclinic TNT orthorhombic TNT monoclinic DNT

tensile modulus, E (GPa)

bulk modulus, K (GPa)

other related explosives, although interactions that are not present in either TNT or 2,4-DNT would have to be added and specifically optimized.

■

AUTHOR INFORMATION

Corresponding Author

shear modulus, G (GPa)

*Tel: +33 4 79758697. Fax: +33 4 79758614. E-mail: sylvie. [email protected].

Pxx

Pyy

Pzz

Pxx, Pyy, Pzz

Pxy

Pxz

Pyz

Notes

10.8

8.4

6.9

10.4

7.5

5.6

7.9

9.6

7.9

6.5

10.3

7.5

4.8

7.3

3.5

6.7

3.8

9.0

1.6

4.3

3.5

■

The authors declare no competing financial interest.

ACKNOWLEDGMENTS This work was granted access to the HPC resources of CCRT/ CINES/IDRIS under the allocations 2011- and 2012-095053 made by GENCI (Grand Equipement National de Calcul Intensif), France. The MUST cluster at the University of Savoie (France) is also acknowledged for the provision of computer time.

DNT are lower. In the literature, E and G have been reported as being ∼6.5 and ∼2.5 GPa, respectively, at room temperature for cast TNT using the disturbance resonance technique. Stress−strain experiments on the same cast samples gave E ≈ 5.9 GPa. It was, however, noted by the authors that the results were dependent upon the state and preparation conditions of the samples, and their TNT bar densities of ∼1625 kg m−3 were lower than the crystalline ones.65 An earlier work reported E to be ∼5.3 GPa with a TNT density of ∼1570 kg m−3.66 Other authors obtain values of ∼3.1 GPa for the TNT tensile modulus but point out that, during the casting and cooling procedure, defects such as dislocations, cracks, porosity, and larger cavities are often introduced.67 Within these conditions, it is not surprising that the experimentally reported moduli values are lower than the model single-crystal ones. However, the force-field gives values that are very reasonable. It also reproduces the experimental observation that G < E < K for TNT.68

■

REFERENCES

(1) Montméat, P.; Thery-Merland, F.; Hairault, L. Techniques de l’Ingénieur 2003, IN 14, 1−8. (2) Pasquinet, E.; Bouvier, C.; Thery-Merland, F.; Hairault, L.; Lebret, B.; Méthivier, C.; Pradier, C. M. J. Colloid Interface Sci. 2004, 272, 21−27. (3) Adhikari, B.; Majumdar, S. Prog. Polym. Sci. 2004, 29, 699−766. (4) Klingenfus, J.; Palmas, P. Phys. Chem. Chem. Phys. 2011, 13, 10661−10669. (5) Vallé, K.; Belleville, P.; Pereira, F.; Sanchez, C. Nat. Mater. 2006, 5, 107−111. (6) Palaniappan, A.; Su, X.; Tay, F. E. H. J. Electroceram. 2006, 16, 503−505. (7) Palaniappan, A.; Li, X.; Tay, F. E. H.; Li, J.; Su, X. Sens. Actuators, B 2006, 119, 220−226. (8) Allen, M. P.; Tildesley, D. J., Computer Simulation of Liquids; Clarendon Press: Oxford, U.K., 1987. (9) Cox, J. R.; Hillier, I. H. Chem. Phys. 1988, 124, 39−46. (10) Clarkson, J.; Ewen Smith, W.; Batchelder, D. N.; Alastair Smith, D.; Coats, A. M. J. Mol. Struct. 2003, 648, 203−214. (11) Alzate, L. F.; Ramos, C. M.; Hernández, N. M.; Hernández, S. P.; Mina, N. Vibrational Spectroscopy 2006, 42, 357−368. (12) Qasim, M.; Kholod, Y.; Gorb, L.; Magers, D.; Honea, P.; Leszczynski, J. Chemosphere 2007, 69, 1144−1150. (13) Cohen, R.; Zeiri, Y.; Wurzberg, E.; Kosloff, R. J. Phys. Chem. A 2007, 111, 11074−11083. (14) Politzer, P.; Murray, J. S.; Koppes, W. M.; Concha, M. C.; Lane, P. Cent. Eur. J. Energ. Mater. 2009, 6, 167−182. (15) Qasim, M.; Gorb, L.; Magers, D.; Honea, P.; Leszczynski, J.; Moore, B.; Taylor, L.; Middleton, M. J. Hazard. Mater. 2009, 167, 154−163. (16) Chen, X.-F.; Hou, C.-Y.; Han, K.-L. J. Phys. Chem. A 2010, 114, 1169−1177. (17) Sorescu, D. C.; Rice, B. M. J. Phys. Chem. C 2010, 114, 6734− 6748. (18) Saloni, J.; Lipkowski, P.; Dasary, S.; Anjaneyulu, Y.; Yu, H.; Hill, G., Jr. Polymer 2011, 52, 1206−1216. (19) Zhou, S.-Q.; Ju, X.-H.; Gu, X.; Zhao, F.-Q.; Yi, J.-H. Struct. Chem. 2012, 23, 921−930. (20) Chua, C. K.; Pumera, M.; Rulisek, L. J. Phys. Chem. C 2012, 116, 4243−4251. (21) Chen, P. C.; Lo, W.; Hu, K. H. J. Mol. Struct. 1997, 389, 91−96. (22) Chen, Y.; Liu, H.; Deng, Y.; Schauki, D.; Fitch, M. J.; Osiander, R.; Dodson, C.; Spicer, J. B.; Shur, M.; Zhang, X.-C. Chem. Phys. Lett. 2004, 400, 357−361. (23) Fayet, G.; Rotureau, P.; Joubert, L.; Adamo, C. J. Hazard. Mater. 2009, 171, 845−850.

5. CONCLUSIONS In this work, an empirical classical all-atom specific force-field has been developed primarily for molecular dynamics simulations of nitroaromatic TNT and applied later to 2,4DNT. The average structural parameters and partial charges were obtained from density functional theory optimizations of the single molecules at the B3LYP/6-311+G** level. Other constants for the classical force-field were initially taken from the literature. They were then adjusted in order to obtain a relatively simple force-field, which is able to reproduce at six different temperatures and within an error of

An Empirically Optimized Classical Force-Field for Molecular

Recommend Documents