Molecular Dynamics Simulations of Liquid Nitromethane Howard E. Alper, Fakhr Abu-Awwad, and Peter Politzer* Department of Chemistry, UniVersity of New Orleans, New Orleans, Louisiana 70148 ReceiVed: June 24, 1999; In Final Form: September 8, 1999 A potential energy function with harmonic intramolecular and Lennard-Jones plus Coulombic intermolecular terms was tested in molecular dynamics simulations of liquid nitromethane. Parameter values were adjusted iteratively until satisfactory agreement with density functional pair calculations and experimental data was achieved. The properties computed using the NVT and NPT ensembles were the heat of vaporization, dielectric constant, self-diffusion coefficient, density, heat capacity at constant pressure, pair correlation functions, single molecule and collective dipole moment reorientation times, the vibrational spectrum, and the effect of increasing pressure upon the C-N stretching frequency. Overall, the results were in reasonable accord with experimental results, the greatest discrepancy being for the dielectric constant. It was concluded, on the basis of the reorientation times and the calculated molecular surface electrostatic potential, that the intermolecular interactions in liquid nitromethane at 1 atm are not highly directional and site-specific. Introduction Molecular dynamics simulations have come to be an impor- tant tool for studying the initiation and propagation of detonation in energetic compounds. (For a recent brief review, see ref 1.) The general objectives are to achieve a better understanding of the process and the factors that are involved, e.g., the localization and transfer of energy, vibrational excitation and bond-breaking, lattice disruption and the role of defects, etc. An eventual goal is to minimize vulnerability to unintended detonations caused by accidental external stimuli, such as impact and shock. Some of the simulations have dealt with idealized model systems; for example, the lattice might be monatomic 2-4 or two- dimensional. 2,3,5-8 In studies of nitromethane, it has sometimes been represented simply by the diatomic CN. 5,9-11 There have also been efforts to simulate actual energetic compounds. In principle, this requires appropriate potential functions for describing both the inter- and intramolecular interactions in the system. A possible first step, however, is to treat just an isolated molecule, thereby eliminating the need for intermolecular potentials. This has been done for several nitramines. 12-15 Another option is to build the crystal lattice and to view it as composed of rigid molecules (i.e., fixed structures). Then only intermolecular interactions need be considered. Such a procedure has now been applied to a large number of energetic com- pounds, 16,17 using a potential function that was shown to be transferable between different chemical classes. The next step, and our objective in the present work, is to develop an approach that takes account of both inter- and intramolecular interactions in order that the simulations can eventually include such features as vibrational excitation, molecular rearrangement, and bond breaking. (These are excluded by the rigid molecule assumption.) Our initial focus is on nitromethane. While it is a liquid at room temperature and pressure, a satisfactory treatment of it should be applicable as well to solids. (A study with a similar objective has very recently been carried out for liquid dimethylnitramine. 18 ) Methods A. Potential Function. All simulations were performed with the molecular dynamics code CHARMM, 19,20 version c25b1, using the potential energy function, The first two summations in eq 1 describe intramolecular interactions between bonded atoms; the quantities (R - R 0 ) and (θ - θ 0 ) are the displacements from the equilibrium bond lengths and angles. The third summation encompasses Lennard- Jones and Coulombic interactions and would involve (for nitromethane) each pair of atoms i and j on different molecules A and B separated by r ij and having charges q i and q j . (For larger molecules, some intramolecular Lennard-Jones and Cou- lombic interactions are also taken into account. 19,20 ) We did not include Urey-Bradley, torsional, or out-of-plane angular terms 20 in eq 1 because (a) the rotational barrier of nitromethane is nearly zero 21 and (b) the vibrational spectrum of nitromethane was found to be reproduced satisfactorily without their inclusion. The assignment of the parameters K R , K θ , ǫ i , ǫ j , σ i and σ j , as well as the atomic charges, will be discussed below. ǫ was set equal to 1. Since the bond-related terms in eq 1 are harmonic, dissocia- tion cannot be simulated directly. It can be inferred, however, since a bond that would in reality undergo dissociation would in the present simulations be observed to have considerably increased vibrational energy and large-amplitude oscillations. 15 B. Parametrization. The starting nitromethane molecular geometry was optimized by a B3P86/6-31+G* density func- V ) ∑ bonds K R (R - R 0 ) 2 + ∑ bond angles K θ (θ - θ 0 ) 2 + ∑ i on A ∑ j on B { (ǫ i ǫ j ) 1/2 [( 0.5(σ i + σ j ) r ij 29 12 - ( 0.5(σ i + σ j ) r ij 29 6 ] + q i q j ∈ r ij } (1) 9738 J. Phys. Chem. B 1999, 103, 9738-9742 10.1021/jp9921102 CCC: $18.00 © 1999 American Chemical Society Published on Web 10/14/1999