Molecular dynamics simulations of surface-initiated melting of nitromethane

The melting of nitromethane initiated at solid-vacuum interfaces has been investigated using molecular dynamics nvt simulations with a realistic force field [D. C. Sorescu et al., J. Phys. Chem. B 104, 8406 (2000)]. The calculated melting point (251+/-5 K) is in good agreement with experiment (244.73 K) and values obtained previously (approximately 255.5 and 266.5+/-8 K) using other simulation methods [P. M. Agrawal et al., J. Chem. Phys. 119, 9617 (2003)]. Analyses of the molecular orientations and diffusion during the simulations as functions of the distance from the exposed surfaces show that the melting is a direct crystal-to-liquid transition, in which the molecules first gain rotational freedom, then mobility. There is a slight dependence of the melting temperature on the exposed crystallographic face.     

Entropy and dynamics of water in hydration layers of a bilayer

We compute the entropy and transport properties of water in the hydration layer of dipalmitoylphosphatidylcholine bilayer by using a recently developed theoretical scheme [two-phase thermodynamic model, termed as 2PT method; S.-T. Lin et al., J. Chem. Phys. 119, 11792 (2003)] based on the translational and rotational velocity autocorrelation functions and their power spectra. The weights of translational and rotational power spectra shift from higher to lower frequency as one goes from the bilayer interface to the bulk. Water molecules near the bilayer head groups have substantially lower entropy (48.36 J/mol/K) than water molecules in the intermediate region (51.36 J/mol/K), which have again lower entropy than the molecules (60.52 J/mol/K) in bulk. Thus, the entropic contribution to the free energy change (TΔS) of transferring an interface water molecule to the bulk is 3.65 kJ/mol and of transferring intermediate water to the bulk is 2.75 kJ/mol at 300 K, which is to be compared with 6.03 kJ/mol for melting of ice at 273 K. The translational diffusion of water in the vicinity of the head groups is found to be in a subdiffusive regime and the rotational diffusion constant increases going away from the interface. This behavior is supported by the slower reorientational relaxation of the dipole vector and OH bond vector of interfacial water. The ratio of reorientational relaxation time for Legendre polynomials of order 1 and 2 is approximately 2 for interface, intermediate, and bulk water, indicating the presence of jump dynamics in these water molecules.     

Molecular dynamics simulations of hexahydro-1,3,5-trinitro-1,3,5-s-triazine (RDX) using a combined Sorescu-Rice-Thompson AMBER force field

We present the results of molecular dynamics simulations of crystalline hexahydro-1,3,5-trinitro-1,3,5-s-triazine (RDX) using the SRT-AMBER force field (P. M. Agrawal et al., J. Phys. Chem. B 2006, 110, 5721), which combines the rigid-molecule force field developed by Sorescu-Rice-Thompson (D. C. Sorescu, B. M. Rice, and D. L. Thompson, J. Phys. Chem. B 1997, 101, 798) with the intramolecular interactions obtained from the Generalized AMBER Force Field (Wang et al., J. Comput. Chem. 2004, 25, 1157). The calculated crystal density at room conditions is about 10% lower than the measured value, while the lattice parameters and thermodynamic melting point are within about 5% at ambient pressure. The chair and inverted chair conformation, bond lengths, and bond angles of the RDX molecule are accurately predicted; however, there are some inaccuracies in the calculated orientations of the NO2 groups. The SRT-AMBER force field predicts overall reasonable results, but modifications, probably in the torsional parameters, are needed for a more accurate force field.     

Polarizable and flexible model for ethanol

A polarizable, flexible model for ethanol is obtained based on an extensive series of B3LYP/6-311++G(d,p) calculations and molecular dynamics simulations. The ethanol model includes electric-field dependence in both the atomic charges and the intramolecular degrees of freedom. Field-dependent intramolecular potentials have been attempted only once previously, for OH and HH stretches in water [P. Cicu et al., J. Chem. Phys. 112, 8267 (2000)]. The torsional potential involving the hydrogen-bonding hydrogen in ethanol is found to be particularly field sensitive. The methodology for developing field-dependent potentials can be readily generalized to other molecules and is discussed in detail. Molecular dynamics simulations of bulk ethanol are performed and the results are assessed based on comparisons with the self-diffusion coefficient [N. Karger et al., J. Chem. Phys. 93, 3437 (1990)], dielectric constant [J. T. Kindt and C. A. Schmuttenmaer, J. Phys. Chem. 100, 10373 (1996)], enthalpy of vaporization [R. C. Wilhoit and B. J. Zwolinski, J. Phys. Chem. Ref. Data, Suppl. 2, 2 (1973)], and experimental interatomic distributions [C. J. Benmore and Y. L. Loh, J. Chem. Phys. 112, 5877 (2000)]. The simultaneous variation of the atomic charges and the intramolecular potentials requires modified equations of motion and a multiple time step algorithm has been implemented to solve these equations. The article concludes with a discussion of the bulk structure and properties with an emphasis on the hydrogen bonding network.     

Molecular dynamics simulations of shock waves in oriented nitromethane single crystals

The structural relaxation of crystalline nitromethane initially at T = 200 K subjected to moderate (~15 GPa) supported shocks on the (100), (010), and (001) crystal planes has been studied using microcanonical molecular dynamics with the nonreactive Sorescu-Rice-Thompson force field [D. C. Sorescu, B. M. Rice, and D. L. Thompson, J. Phys. Chem. B 104, 8406 (2000)]. The responses to the shocks were determined by monitoring the mass density, the intermolecular, intramolecular, and total temperatures (average kinetic energies), the partitioning of total kinetic energy among Cartesian directions, the radial distribution functions for directions perpendicular to those of shock propagation, the mean-square displacements in directions perpendicular to those of shock propagation, and the time dependence of molecular rotational relaxation as a function of time. The results show that the mechanical response of crystalline nitromethane strongly depends on the orientation of the shock wave. Shocks propagating along [100] and [001] result in translational disordering in some crystal planes but not in others, a phenomenon that we refer to as plane-specific disordering; whereas for [010] the shock-induced stresses are relieved by a complicated structural rearrangement that leads to a paracrystalline structure. The plane-specific translational disordering is more complete by the end of the simulations (~6 ps) for shock propagation along [001] than along [100]. Transient excitation of the intermolecular degrees of freedom occurs in the immediate vicinity of the shock front for all three orientations; the effect is most pronounced for the [010] shock. In all three cases excitation of molecular vibrations occurs more slowly than the intermolecular excitation. The intermolecular and intramolecular temperatures are nearly equal by the end of the simulations, with 400-500 K of net shock heating. Results for two-dimensional mean-square molecular center-of-mass displacements, calculated as a function of time since shock wave passage in planes perpendicular to the direction of shock propagation, show that the molecular translational mobility in the picoseconds following shock wave passage is greatest for [001] and least for the [010] case. In all cases the root-mean-square center-of-mass displacement is small compared to the molecular diameter of nitromethane on the time scale of the simulations. The calculated time scales for the approach to thermal equilibrium are generally consistent with the predictions of a recent theoretical analysis due to Hooper [J. Chem. Phys. 132, 014507 (2010)].     

Molecular nitrogen-N2 properties: the intermolecular potential and the equation of state

Quantum mechanical (QM) high precision calculations were used to determine N(2)-N(2) intermolecular interaction potential. Using QM numerical data the anisotropic potential energy surface was obtained for all orientations of the pair of the nitrogen molecules in the rotation invariant form. The new N(2)-N(2) potential is in reasonably good agreement with the scaled potential obtained by van der Avoird et al. using the results of Hartree-Fock calculations [J. Chem. Phys. 84, 1629 (1986)]. The molecular dynamics (MD) of the N(2) molecules has been used to determine nitrogen equation of state. The classical motion of N(2) molecules was integrated in rigid rotor approximation, i.e., it accounted only translational and rotational degrees of freedom. Fincham [Mol. Simul. 11, 79 (1993)] algorithm was shown to be superior in terms of precision and energy stability to other algorithms, including Singer [Mol. Phys. 33, 1757 (1977)], fifth order predictor-corrector, or Runge-Kutta, and was therefore used in the MD modeling of the nitrogen pressure [S. Krukowski and P. Strak, J. Chem. Phys. 124, 134501 (2006)]. Nitrogen equation of state at pressures up to 30 GPa (300 kbars) and temperatures from the room temperature to 2000 K was obtained using MD simulation results. Results of MD simulations are in very good agreement (the error below 1%) with the experimental data on nitrogen equation of state at pressures below 1 GPa (10 kbars) for temperatures below 1800 K [R. T. Jacobsen et al., J. Phys. Chem. Ref. Data 15, 735 (1986)]. For higher temperatures, the deviation is slightly larger, about 2.5% which still is a very good agreement. The slightly larger difference may be attributed to the vibrational motion not accounted explicitly by rigid rotor approximation, which may be especially important at high temperatures. These results allow to obtain reliable equation of state of nitrogen for pressures up to 30 GPa (300 kbars), i.e., close to molecular nitrogen stability limit, determined by Nellis et al. [Phys. Rev. Lett. 53, 1661 (1984)].     

Temperature of systems out of thermodynamic equilibrium

Two phenomenological approaches are currently used in the study of the vitreous state. One is based on the concept of fictive temperature introduced by Tool [J. Res. Natl. Bur. Stand. 34, 199 (1945)] and recently revisited by Nieuwenhuizen [Phys. Rev. Lett. 80, 5580 (1998)]. The other is based on the thermodynamics of irreversible processes initiated by De Donder at the beginning of the last century [L'Affinite (Gauthier-Villars, Paris, 1927)] and recently used by Moller et al. for a thorough study of the glass transition [J. Chem. Phys. 125, 094505 (2006)]. This latter approach leads to the possibility of describing the glass transition by means of the freezing-in of one or more order parameters connected to the internal structural degrees of freedom involved in the vitrification process. In this paper, the equivalence of the two preceding approaches is demonstrated, not only for glasses but in a very general way for any system undergoing an irreversible transformation. This equivalence allows the definition of an effective temperature for all systems departed from equilibrium generating a positive amount of entropy. In fact, the initial fictive temperature concept of Tool leads to the generalization of the notion of temperature for systems out of thermodynamic equilibrium, for which glasses are just particular cases.     

Communication: Thermodynamic scaling of the Debye process in primary alcohols

The molecular dynamics of hydrogen-bonded liquids usually does not satisfy the thermodynamic scaling. However, very recently, two opposite conclusions about validity of thermodynamical scaling in monohydroxy alcohol, 2-ethyl-1-hexanol, were presented by Reiser et al. [J. Chem. Phys. 132, 181101 (2010)] and Fragiadakis et al. [J. Chem. Phys. 132, 144505 (2010)]. In this communication we present new experimental results that can explain this ostensible contradiction.     

Role of defects in the nonmonotonic behavior of secondary relaxation of polypropylene glycols

A nonmonotonic relaxation kinetic model [Ya. Ryabov et al., J. Phys. Chem. B 105, 1845 (2001)] is successfully applied to describe an intriguing slow down in the dielectric secondary gamma relaxation of polypropylene glycols (PPGs) with increasing temperature near the glass transition. The anomalous behavior is interpreted as a result of two simultaneous events: A thermal activation and a defect formation in the hydrogen bonded network formed by molecules of PPGs. This new insight into the molecular mechanism, which is responsible for the suggested sensitivity of the secondary process in PPGs to the glass transition phenomenon, is compared to our previous results obtained in terms of the minimal model for secondary relaxations.     

On the accuracy of force fields for predicting the physical properties of dimethylnitramine

The accuracy of three force fields for predicting the physical properties of dimethylnitramine (DMNA) has been investigated by using molecular dynamics simulations. The Sorescu, Rice, and Thompson (SRT) (J. Phys. Chem. B 1997, 101, 798) rigid-molecule, flexible generalized AMBER (J. Comput. Chem. 2004, 25, 1157), and Smith et al. flexible force fields (J. Phys. Chem. B 1999, 103, 705) were tested. The density, lattice parameters, isotherm, and melting point of DMNA are calculated using classical molecular dynamics. Except for the melting point, the predictions of the three force fields are in reasonable agreement with experimental values. The calculated thermodynamic melting points (Tmp) for the SRT, AMBER, and Smith et al. force fields are 380, 360, and 260 K, respectively. The experimental value is 331 K. Modifications of the torsional barriers in the AMBER force field resulted in Tmp = 346 K, in good agreement with the experimental value of 331 K. The calculated lattice parameters and bulk modulus are also improved with the modifications of the AMBER potential. The results indicate that, although not sufficiently accurate without modifications, the general force fields such as AMBER provide the basis for developing force fields that correctly predict the physical properties of nitramines.     

A molecular-dynamics study of structural and physical properties of nitromethane nanoparticles

The structural and physical properties of nanoparticles of nitromethane are studied by using molecular dynamics methods with a previously developed force field. [Agrawal et al., J. Chem. Phys. 119, 9617 (2003).] This force field accurately predicts solid- and liquid-state properties as well as melting of bulk nitromethane. Molecular dynamics simulations of nanoparticles with 480, 240, 144, 96, 48, and 32 nitromethane molecules have been carried out at various temperatures. The carbon-carbon radial distribution function, dipole-dipole correlation function, core density, internal enthalpy, and atomic diffusion coefficients of the nanoparticles were calculated at each temperature. These properties were used to characterize the physical phases and thus determine the melting transitions of the nanoparticles. The melting temperatures predicted by the various properties are consistent with one another and show that the melting temperature increases with particle size, approaching the bulk limit for the largest particle. A size dependence of melting points has been observed in experimental and theoretical studies of atomic nanoparticles, and this is a further demonstration of the effect for large nanoparticles of complex molecular materials.     

Affinity and its derivatives in the glass transition process

The thermodynamic treatment of the glass transition remains an issue of intense debate. When associated with the formalism of non-equilibrium thermodynamics, the lattice-hole theory of liquids can provide new insight in this direction, as has been shown by Schmelzer and Gutzow [J. Chem. Phys. 125, 184511 (2006)], by Mo?ller et al. [J. Chem. Phys. 125, 094505 (2006)], and more recently by Tropin et al. [J. Non-Cryst. Solids 357, 1291 (2011); ibid. 357, 1303 (2011)]. Here, we employ a similar approach. We include pressure as an additional variable, in order to account for the freezing-in of structural degrees of freedom upon pressure increase. Second, we demonstrate that important terms concerning first order derivatives of the affinity-driving-force with respect to temperature and pressure have been previously neglected. We show that these are of crucial importance in the approach. Macroscopic non-equilibrium thermodynamics is used to enlighten these contributions in the derivation of C(p),κ(T), and α(p). The coefficients are calculated as a function of pressure and temperature following different theoretical protocols, revealing classical aspects of vitrification and structural recovery processes. Finally, we demonstrate that a simple minimalist model such as the lattice-hole theory of liquids, when being associated with rigorous use of macroscopic non-equilibrium thermodynamics, is able to account for the primary features of the glass transition phenomenology. Notwithstanding its simplicity and its limits, this approach can be used as a very pedagogical tool to provide a physical understanding on the underlying thermodynamics which governs the glass transition process.     

The photodissociation of the water dimer in the A band: a twelve-dimensional quasiclassical study

The quasiclassical absorption spectrum of the water dimer in the A band was calculated taking into account motion in all degrees of freedom of the system. The ab initio excited state potentials employed were interpolated by the modified Shepard interpolation method using QMRCI energies and state-averaged MCSCF gradients and Hessians. The ground state vibrational wavefunction was variationally calculated using an adiabatic separation between the high and low frequency normal modes of the system. The calculated spectrum of water dimer shows a clear blueshift with respect to the monomer, but also a small red tail, in agreement with the prediction by Harvey et al. [J. Chem. Phys. 109, 8747 (1998)]. Previous three-dimensional model studies of the photodissociation of the water dimer by Valenzano et al. [J. Chem. Phys. 123, 034303 (2005)] did not show this red tail. A thorough analysis of the dependence of the spectrum on the modes coupled explicitly in the calculation of the spectrum shows that the red tail is due to coupling between the intramolecular stretch vibrations on different monomers.     

Simulations of the solid, liquid, and melting of 1-n-butyl-4-amino-1,2,4-triazolium bromide

Molecular dynamics simulations are used to study the solid and liquid properties and to predict the melting point of 1-n-propyl-4-amino-1,2,4-triazolium bromide ([patr][Br]) using a force field based on the one developed by Canongia Lopes et al. (J. Phys. Chem. B 2004, 108, 2038) for dialkyl substituted imidazolium salts, which was modified by including terms from the general AMBER force field. Electrostatic charges for the intermolecular interactions were determined from gas-phase ab initio electron structure calculations of the triazolium cation. Simulations of the solid state at 100 K reproduced the experimental density to within 4%. Simulations from 100 K to the melting point and the liquid from 333 to 500 K were performed to determine the temperature dependence of the densities of the two phases. The structures of the solid and liquid phases are characterized with radial distribution functions, which show that there are strong spatial correlations among neighboring ion pairs in liquid [patr][Br]. The dynamic behavior of the ions in the liquid state is also studied by computing velocity autocorrelation functions and the mean-square displacements between the ions. The melting point is determined by simulating void-induced melting. Changes in the density, intermolecular energy, and Lindemann index are used as indicators of the melting transition. The computed melting point is 360 +/- 10 K, which is within 10% of the experimental value 333 K.     

Three-dimensional potential energy surface of the Ar-OH(2Pi i) complex

Pure rotational transitions in the ground state for Ar-OH and Ar-OD [Y. Ohshima et al., J. Chem. Phys. 95, 7001 (1991) and Y. Endo et al., Faraday Discuss. 97, 341 (1994)], those in the excited states of the OH vibration, nu(s)=1 and 2, observed by Fourier-transform microwave spectroscopy in the present study, rotation-vibration transitions observed by infrared-ultraviolet double-resonance spectroscopy [K. M. Beck et al., Chem. Phys. Lett. 162, 203 (1989) and R. T. Bonn et al., J. Chem. Phys. 112, 4942 (2000)], and the P-level structure observed by stimulated emission pumping spectroscopy [M. T. Berry et al., Chem. Phys. Lett. 178, 301 (1991)] have been simultaneously analyzed to determine the potential energy surface of Ar-OH in the ground state. A Schrodinger equation, considering all the freedom of motions for an atom-diatom system in the Jacobi coordinate, R, theta, and r, was numerically solved to obtain energies of the rovibrational energy levels using the discrete variable representation method. A three-dimensional potential energy surface is determined by a least-squares fitting. In the analysis the potential parameters, obtained by ab initio calculations at the RCCSD(T) level of theory with a set of basis functions of aug-cc-pVTZ and midbond functions, are used as initial values. The determined intermolecular potential energy surface and its dependence on the OH monomer bond length are compared with those of an isovalent radical complex, Ar-SH.     

Reactivity enhancement of ultracold O(3P)+H2 collisions by van der Waals interactions

The role of van der Waals forces in O((3)P)+H(2)(upsilon=1,j=0) collisions is investigated theoretically at low and ultralow temperatures. Quantum scattering calculations have been performed for zero total angular momentum using the lowest London-Eyring-Polanyi-Sato double-polynomial (3)A(") potential-energy surface reported by [Rogers et al., J. Phys. Chem. A 104, 2308 (2000)] and its recent BMS1 and BMS2 extensions developed by [Brandao et al., J. Chem. Phys. 121, 8861 (2004)] which provide a more accurate treatment of the van der Waals interaction. Our calculations show that van der Waals forces strongly influence chemical reactivity at ultracold translational energies. The presence of a zero-energy resonance for the BMS1 surface is found to enhance reactivity in the ultracold regime and shift the Wigner threshold to lower temperatures.     

Molecular dynamics simulations of melting of perfect crystalline hexahydro-1,3,5-trinitro-1,3,5-s-triazine

The melting mechanism of superheated perfect crystalline hexahydro-1,3,5-trinitro-1,3,5-s-triazine (alpha-RDX) has been investigated using molecular dynamics simulations with the fully flexible force field developed by Smith and Bharadwaj [J. Phys. Chem. B 103, 3570 (1999)]. Sequential 50 ps equilibration simulations of the constant stress-constant temperature ensemble were performed at 10 K intervals over the range of 300-650 K, corresponding to a heating rate of 2.0 x 10(11) Ks. A solid-solid phase transition is observed between 480 and 490 K, followed by melting, which occurs between 500 and 510 K. The solid-solid phase transition, both displacive and rotational, is characterized by an abrupt decrease in the lengths of the unit cell edges a and b and an increase of the length of edge c. The molecular conformation in the new phase is AAE, although the axial nitro groups have different changes: one shift is more axial and the other is more equatorial. Phases other than alpha-RDX have been observed experimentally, however, there are insufficient data for comparisons to ascertain that the new phase observed here corresponds to a real phase. At the high heating rate (2.0 x 10(11) Ks) used in the simulations, the melted RDX reaches full orientational disorder at about 540 K and translational freedom at around 580 K. If the simulation at the melting temperature (510 K) is run sufficiently long complete rotational freedom is achieved in a few hundreds of picoseconds, while complete translational freedom requires much longer. These results show that given a sufficiently high heating rate, the system can exist for significant periods of time in a near-liquid state in which the molecules are not as free to rotate and diffuse as in the true liquid state. The bond lengths and bond angles undergo little change upon melting, while there are significant changes in the dihedral angles. The molecular conformation of RDX changes from AAE to EEE upon melting. The ramification of this for formulating force fields that accurately describe melting is that it is important that the torsional motions are accurately described.     

Correlation between thermodynamic and kinetic fragilities in nonpolymeric glass-forming liquids

A phenomenological relationship between reduced excess heat capacity of supercooled liquid DeltaC(p)(exc)(T(g))DeltaS(m) at the glass transition temperature T(g), fragility index m, and reduced glass transition temperature T(rg)=T(g)T(m), where T(m) is the melting (liquidus) temperature, was derived for fragile nonpolymeric glass-forming liquids under the assumptions that the fragile behavior of these liquids is described by the Vogel-Fulcher-Tammann (VFT) equation; the excess heat capacity of liquid is inversely proportional to the absolute temperature and the VFT temperature T(0) is equal to the Kauzmann temperature T(K). It was found that DeltaC(p)(exc)(T(g))DeltaS(m) is a composite function of m and T(rg), which indicates that the empirical correlation DeltaC(p)(exc)(T(g))DeltaS(m)=0.025m recently identified by Wang et al. [J. Chem Phys. 125, 074505 (2006)] is probably valid only for liquids which have nearly the same values of T(rg).     

System-size dependence of the free energy of crystalline solids

We investigate the system-size dependence of the Helmholtz free energy of crystalline solids from computer simulation. We employ a standard thermodynamic integration technique along a reversible path that links the crystalline solid with a noninteracting Einstein crystal with the same structure. The key contribution to the free energy is computed by using the so-called expanded-ensemble technique and the results are compared with those obtained from conventional integration of the derivative of the free energy along the path using Gaussian-Legendre quadrature. We find that both methods yield fully consistent results. The free energy is found to exhibit a strong dependence with system size, in agreement with the behavior found by Polson et al. [J. Chem. Phys. 112, 5339 (2000)] but at variance with the dependence reported more recently by Chang and Sandler [J. Chem. Phys. 118, 8390 (2003)]. This has been tested for the face-centered cubic (fcc) and hexagonal close-packed phases of a crystal of hard spheres at a density close to the melting point. We also investigate any possible dependence of the free energy of the solid phase with the shape of the simulation box. We find that this contribution may not be as important as previous investigations suggest. The present results seem to indicate that there is a non-negligible contribution to the free energy arising from the orientation of the closed-packed crystal layers with respect to the simulation cell. This contribution is particularly noticeable for small system sizes and is believed to be an effect of the periodic boundary conditions used in the simulations. The results presented here corroborate the stability of the fcc phase of the hard-sphere solid close to melting.     

Long-range-corrected hybrids using a range-separated Perdew-Burke-Ernzerhof functional and random phase approximation correlation

We build on methods combining a short-range density functional approximation with a long-range random phase approximation [B. G. Janesko, T. M. Henderson, and G. E. Scuseria, J. Chem. Phys. 130, 081105 (2009)] or second-order screened exchange [J. Paier et al., J. Chem. Phys. 132, 094103 (2010)] by replacing the range-separated local density approximation functional with a range-separated generalized gradient approximation functional in the short range. We present benchmark results that show a marked improvement in the thermodynamic tests over the previous local density approximation-based methods while retaining those methods' excellent performance in van der Waals interactions.