Theoretical Study of Shocked Formic Acid: Born–Oppenheimer MD

Dec 10, 2015 - Although the shock Hugoniot curve calculated using nonreactive classical MD for formic acid is in reasonable agreement with one set of ...
0 downloads 5 Views 2MB Size
Download PDF
Article pubs.acs.org/JPCB

Theoretical Study of Shocked Formic Acid: Born−Oppenheimer MD Calculations of the Shock Hugoniot and Early-Stage Chemistry Betsy M. Rice* and Edward F. C. Byrd U.S. Army Research Laboratory (ARL), Aberdeen Proving Ground, Aberdeen 21005, Maryland, United States S Supporting Information *

ABSTRACT: Quantum and classical molecular dynamics simulations are used to explore whether chemical reactivity of shocked formic acid occurs at pressures greater than 15 GPa, a question arising from results of different shock compression experiments. The classical molecular dynamics simulations were performed using a quantum-based nonreactive pair additive interaction potential whereas the full resolution quantum mechanical molecular dynamics simulations allow chemical reactions. Although the shock Hugoniot curve calculated using nonreactive classical MD for formic acid is in reasonable agreement with one set of experimental results, shock Hugoniot points calculated using Born− Oppenheimer MD at 30 GPa are in agreement with the set of experimental data that suggests chemical reactivity at these elevated temperatures and pressures. Examination of atomic positions throughout the Born−Oppenheimer MD trajectories clearly indicates extensive and complex chemical reaction, chiefly involving hydrogen-atom transfer and intermolecular complexation.

sponding to a shock state.7,8 The earlier of the two studies used an electronic-structure-based MD method and an O(N) tightbinding method to understand a “bend” in the shock Hugoniot curve of shocked benzene.7 Their calculations demonstrate polymerization of a system of 576 benzene molecules within 1 ps upon shock compression. Longer (7 ps) MD simulations of systems of 54 molecules at thermodynamic states corresponding to those of the NEMD shocked states produced similar chemical moieties as found in the NEMD simulations. A similar and more recent study8 used SCC-DFTB and a computationally efficient extended Lagrangian Born−Oppenheimer MD (denoted XL-BOMD) method7,8 to simulate the shock response of a system of 32 phenylacetylene molecules for 50 ps; these calculations also showed shocked-induced polymerization. In this work, we provide a comparison of the experimental data with predicted shock Hugoniot values from atomistic simulations using both classical9,10 and Born−Oppenheimer molecular dynamics (BOMD) simulations using dispersion corrected density-functional theory (DFT). While this approach is computationally more expensive than those using the SCC-DFTB method, we desired to minimize any empiricism in our method. Therefore, our simulation times are more limited; however, our goal is to simply determine whether chemical reactivity occurs at conditions consistent with those of a shocked state similar to that of Trunin et al.1 The classical mechanical molecular dynamics simulations9,10 used a

1. INTRODUCTION The question of whether formic acid undergoes chemical reaction at shock conditions above 15 GPa arose when results of shock compression experiments by two groups produced conflicting results. The shock compression of formic acid by the reflection method (Trunin et al.1) produced a shock Hugoniot that exhibits a discontinuity at pressures corresponding to ∼15 GPa. This behavior has been interpreted as due to chemical reaction, such as compositional change or polymerization, and is not unexpected. Formic acid is known to polymerize under high pressure at room temperature,2 as well as decompose at high temperatures and pressures.3 A second study in which formic acid was shock compressed in gas gun-driven plate impact experiments (Manner et al.)4 produced a shock Hugoniot curve that did not exhibit the discontinuity at pressures above 15 GPa, indicating that no reaction occurred at the pressures observed by Trunin et al.1 It is unlikely that any atomistic simulation could directly mimic the experimental conditions that might reveal the origin of the discrepancy between the two experimental studies; however, atomistic simulation using quantum mechanical methods can confirm whether reactivity occurs at the elevated temperatures and pressures corresponding to the Trunin et al.1 studies. Numerous examples exist in the literature in which electronic-structure-based MD methods are being used to study shock-initiated chemistry,5,6 with many utilizing the semiempirical self-consistent-charge density-functional tightbinding (SCC-DFTB) approach coupled with the multiresolution multiscale shock technique (MSST). Of particular relevance to this study are two papers that demonstrate shockinduced polymerization; both use direct nonequilibrium MD (NEMD) simulations and shrinking boundary conditions to mimic compression to study materials at conditions correThis article not subject to U.S. Copyright. Published XXXX by the American Chemical Society

Special Issue: Bruce C. Garrett Festschrift Received: September 10, 2015 Revised: December 10, 2015

A

DOI: 10.1021/acs.jpcb.5b08845 J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

composed of an equilibration simulation of 25 ps, followed by a long production run from which averages were obtained. The length of the production portion of each trajectory ranged from 100 to 180 ps, sufficient to ensure the averaged thermodynamic properties remained constant with increased time. We also calculated the universal Hugoniot curve (UHC) for liquids15 for this model

nonreactive interaction potential for formic acid to determine the unreacted shock Hugoniot of formic acid for comparison with experiment. The BOMD simulations were used to generate the shock Hugoniot state at 30 GPa for comparison with experiment, as well as to explore chemical reactivity at conditions near this point. Section 2 will describe the computational methods used in the study, while section 3 will discuss the results. We will provide a summary of the results and conclusions in section 4.

U = 1.37c0 − 0.37c0 exp( −2u/c0) + 1.62u

(1)

where U is the shock velocity, u is the particle velocity, and c0 is the sound speed at 1 bar. The UHC is an empirical description that relates the sound speed of a material in its initial state with its unreacted shock Hugoniot. Deviations from this curve, such as those corresponding to the Trunin et al. data,1 are assumed to be indicative of chemical reaction or phase transition.16 The sound speed for this model (1.47 km/s) was estimated from results of NPT-MD simulations (see procedure described in ref 17), and is faster than those used by Manner et al.4 (1.28 km/s) and Trunin et al.1 (1.34 km/s). 2.2. Quantum MD (Born−Oppenheimer MD) Simulations. Born−Oppenheimer MD simulations (hereafter denoted as BOMD) within the isothermal−isobaric (NPT) ensemble using the CP2K program suite18 were used to determine shock Hugoniot points for comparison with experiment, and to explore reactivity at a pressure higher than that of the discontinuity in the Trunin shock Hugoniot curve;1 we chose 30 GPa. Several challenges exist for reactive molecular dynamics; molecular simulations are to mimic the experimental condition, in this case, a system at a specific thermodynamic state (i.e., temperature and pressure). For such a case, NPT-BOMD simulations would be appropriate. However, such simulations involve coupling the system to a barostat and a thermostat, through which the system is driven to the desired thermodynamic state in the simulation. Any energy release or absorption from the reactions might be masked by the imposed temperature and pressure coupling, potentially affecting any chemical reaction occurring in the system. Further, molecular simulations usually require an equilibration period during which the system relaxes from the initial state. In any reactive MD simulation, reactive events during the equilibration period can occur, and cannot be assumed to be representative of the chemical events at the desired thermodynamic state. Furthermore, the chemistry occurring within the nonequilibrated system could influence subsequent chemistry at later times in the simulations when the target thermodynamic state is reached. An additional challenge is faced when using BOMD as simulation times are limited due to computational constraints; thus, we can only say that we have reached an equilibrium on the time scales of the simulation (a few picoseconds at the current state of the art). Recognizing these challenges, we proceeded with our study by first performing an NPT-BOMD simulation at 0 GPa, 298 K (denoted RTP), for 2.1 ps to provide the thermodynamic reference state (P0, V0) needed to calculate a shock Hugoniot point. NPT-BOMD simulations were then performed at 30 GPa for two temperatures, 1550 and 2050 K for 1.3 and 3.2 ps, respectively; we determined that these temperatures bounded the Hugoniot temperature. The system was composed of 157 molecules, and the initial configurations of the system at 1550 and 2050 K, 30 GPa, were generated as follows: first, classical isothermal−isobaric simulations using the Roszek et al. nonreactive potential energy function12 were performed at 1550 K, 30 GPa, to generate a starting structure for the BOMD

2. COMPUTATIONAL METHODS 2.1. Classical MD Simulations. The shock Hugoniot curve of formic acid was calculated using classical molecular dynamics simulations and the procedure developed by J. J. Erpenbeck.11 In this method, the conditions for which the Hugoniot relation is satisfied [i.e., HG = 0 where HG = E − E0 − 1/2(P + P0)(V0 − V)] are extracted from the formic acid equation of state (EOS) evaluated over a range of temperature (T), volume (V), and pressure (P) values using isothermal−isochoric molecular dynamics (NVT-MD). Each point on the shock Hugoniot curve is obtained by first generating a set of EOS points, evaluating the Hugoniot function HG for each point in the set, and then fitting this set of HG to a polynomial in temperature, whose root is the Hugoniot temperature TH. The Hugoniot pressure, PH, is obtained by fitting the system pressures within each set to a polynomial in temperature, and evaluating at TH. The thermodynamic reference state (P0, V0) used in the Hugoniot relation was predicted for liquid formic acid at T = 298 K, 1 atm, using isothermal−isobaric (NPT-MD) simulations. The intra- and intermolecular potentials used for simulation of formic acid were developed by Roszak et al.12 and are parametrized using information from quantum mechanical calculations of small formic acid clusters. The intramolecular terms consist of harmonic functions to describe bond stretching and angle bending and cosine-type functions to describe torsional motion. The intermolecular potential consists of the superposition of pairwise Lennard-Jones 9−6 and Coulombic potential terms. The Lennard-Jones 9−6 parameters and mixing rules for various types of atomic pairs as previously published12 have been used in the present study without change. Calculations by Roszak et al.12 of the melting temperature and of the 250 K pressure−volume isotherm of the crystal were in good agreement with experiment;12 we have reproduced the isotherm to ensure that our understanding and implementation of the force field were correct. All classical MD simulations were carried out using the DL_POLY_2.0 set of codes,13 and all simulations used a cubic cell containing 648 molecules. Periodic boundary conditions were imposed in all dimensions, and a cutoff distance of 10 Å was used for the calculation of interatomic forces. Standard long-range corrections for the energy and virial contributions were included, and Ewald’s method was used to treat the Coulombic long-range interactions.14 A time step of 1 fs was used in all the simulations except for those at higher densities; the time steps were set to 0.5 and 0.3 fs for simulations at 0.3772 and 0.3601 cc/g, respectively. Leap-frog Verlet algorithms as implemented in DL_POLY_2.0 were used to integrate the Nosé−Hoover equations of motion for both the NPT-MD simulations used in generating the reference state and the NVT-MD simulations used to generate the Hugoniot curve. The Nosé−Hoover thermostat and barostat relaxation constants are 0.1 and 1.0 ps, respectively. All trajectories were B

DOI: 10.1021/acs.jpcb.5b08845 J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

is close to the derived shock Hugoniot point, and is not actually on the shock Hugoniot curve. All fully periodic BOMD simulations were done using the Quickstep method19 implemented in CP2K with the Perdew− Burke−Ernzerhof (PBE)20,21 generalized gradient approximation (GGA) density-functional theory (DFT). The DZVPGTH basis set was used in conjunction with the GTH-DCACP pseudopotential22 and a 2000 Ry cutoff for the auxiliary plane wave basis set. The van der Waals correcting DCACP pseudopotentials were used as they have been shown to accurately predict crystal structural features of a variety of energetic molecular crystals at ambient conditions.23,24 The electronic energy was converged to an accuracy of 1 × 10−6 hartree using the orbital transformation (OT) direct energy minimization scheme.25,26 The OT method is guaranteed convergence if the energy is reduced each step and scales very well for large systems. However, this method precludes the use of orbital mixing; we also did not employ any electronic smearing. A time step of 0.25 fs for the Nosé−Hoover equations of motion was used in all simulations (time constants of 50 and 100 fs for the barostat and thermostat, respectively).

simulations. Initial atomic velocities for the corresponding temperatures were randomly assigned by CP2K. Using the Erpenbeck procedure11 described earlier, the HG values for the T = 1550 and 2050 K results were determined, and a linear fit to these two points produced the Hugoniot temperature TH. The Hugoniot volume, VH, was obtained by fitting the system volumes for the two temperatures to a linear equation in temperature, and evaluating this equation at TH. As will be described hereafter, we observed chemical reactions during the equilibration periods of the NPT-BOMD simulations used to calculate the Hugoniot point, with reactions continuing for the duration of the trajectory. It was not known whether the barostat or thermostat significantly influenced the observed chemistry, or whether the chemistry was the result of a highly nonequilibrium initial state. To explore the first possibility, we performed a microcanonical (NVE-BOMD) simulation using a system configuration from the NPT-BOMD simulation at 1550 K to determine (a) whether chemistry would continue and (b) whether there would be any changes in energy, temperature, and/or pressure accompanying any subsequent chemical reactivity in the absence of the barostat and thermostat. This NVE-BOMD simulation was initiated using the configuration at 1.0 ps in the 1550 K 30 GPa NPTBOMD simulation, and integrated for 1.3 ps. We also performed a second NPT-BOMD simulation at 1550 K, 30 GPa, using different initial conditions than the first, to determine if the initial conditions affected subsequent chemistry and the final thermodynamic state. The initial configuration of the system used in the second NPT-BOMD simulation was similar to the first, except each edge of the simulation cell was larger by 20%; the molecular mass centers were located at the same fractional positions as in the first simulation, and the body-fixed coordinates remained the same. We noted the initial state of the first NPT-BOMD simulation led to an extremely high temperature rise at the beginning of the trajectory, with chemical reaction observed as the system relaxed to the target temperature (i.e., 1550 K). As we were concerned this early time “hot” portion of the simulation would initiate reactions that would normally not be seen at 1550 K (and thus influence the late time trajectory), we enforced a “cold” start for the second simulation by stopping the trajectory every time the instantaneous temperature exceeded the target temperature of 1550 K. We restarted it with rescaled velocities using the configuration corresponding to the last time step at which the temperature was under 1000 K; the atomic velocities corresponding to this configuration were scaled to 300 K. Once the temporal behavior of the system temperature did not exceed 1550 K over a period of 250 fs after rescaling (this took five rescaling runs to accomplish), we then proceeded with the NPT-BOMD simulation having a target state of 1550 K, 30 GPa, and no further rescaling of atomic velocities. The total integration time for the second NPT-BOMD simulation is 3.3 ps, not counting the five restarted scaling runs (these totaled ∼6.5 fs). A second shock Hugoniot point was calculated using the results of this simulation. All trajectories were integrated until the thermodynamic properties appeared to reach an equilibrium. Averages were obtained from the last 200 fs of each BOMD trajectory and used in the calculation of the shock Hugoniot point. For convenience, we will refer to the first and second NPT-BOMD simulations at 1550 K, 30 GPa, as “hot” and “cold” (referring to the high initial temperature in the first and the rescaled lower temperature of the second). We emphasize that these two simulations target an EOS point that

3. RESULTS The various EOS points used to generate the shock Hugoniot curve reported in refs 9 and 10 were calculated using classical NVT-MD simulations for densities ranging from 0.5974 to 0.3601 cc/g and temperatures ranging from 325 up to 2900 K, and are given in Supporting Information Table S1. For each density, the HG and PH were fitted to a linear equation in temperature, and the Hugoniot temperature and corresponding Hugoniot pressure were obtained as described in section 2. These values are given in Table 1, along with shock velocity, particle velocity, and relative volume change of the solution using the relations given in Manner et al.4 Table 1. Shock Hugoniot States for Formic Acid V (cc/g)

P (GPa)

0.5974 0.5548 0.5122 0.4747 0.4355 0.4070 0.3920 0.3772 0.3601

1.099 1.985 4.050 7.399 14.130 22.860 29.570 38.320 51.940

0.4046b 0.4042c

30.000 30.000

T (K)

Ust (mm/μs)

Classical NVT-MD 332 2.266 389 2.572 456 3.238 577 4.002 834 5.109 1219 6.178 1552 6.855 2010 7.624 2794 8.652 DFT NPT-BOMD 1865 6.975 1874 6.978

up (mm/μs)

ΔV/V0a (up/Ust)

0.341 0.543 0.879 1.300 1.945 2.602 3.033 3.534 4.221

0.150 0.211 0.272 0.325 0.381 0.421 0.443 0.464 0.488

3.698 3.681

0.527 0.528

a

V0 = 0.7032 cc/g for the classical model, V0 = 0.8574 cc/g for the DFT system. bDFT value obtained using the “hot” 1550 K, 30 GPa data point. cDFT value obtained using the “cold” 1550 K, 30 GPa data point.

BOMD results for the two shock Hugoniot points at 30 GPa are also given in Table 1; BOMD data used to generate the shock Hugoniot points are given in Table S2. Table S2 also reports results of the NVE-BOMD simulations that were performed to determine if the thermodynamic state would significantly change in the absence of a barostat and thermostat C

DOI: 10.1021/acs.jpcb.5b08845 J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

Figure 1. (a) Shock velocity versus particle velocity resulting from shock Hugoniot calculations using classical (denoted MD)9,10 and quantum (denoted DFT) MD simulations. Experimental data from Manner et al.4 (denoted Manner) and Trunin et al.1 (denoted Trunin) are also included. Black and red lines show the UHC for the classical MD results and Manner et al.4 data, respectively. (b) Pressure versus relative volume change for all data shown in part a. Note there are two DFT results on top of one another.

MD results, clearly supporting the case for nonreactivity at shock conditions up to ∼23 GPa in these experiments. The better agreement of the BOMD data with the Trunin et al.1 data at shock conditions well above 15 GPa, the regime where the Trunin et al.1 data clearly diverge with the UHCs and the nonreactive MD results, is suggestive of a phase transition or chemical reaction in these simulations. Therefore, we searched for this possibility in the BOMD-simulated material from the “hot” and “cold” simulations. We monitored all interatomic distances of the original covalently bound atom pairs in each formic acid molecule over the duration of the trajectories. As discussed below, for almost all molecules in the trajectories, a hydrogen atom left the parent molecule. Therefore, we also monitored the distance between any nearest neighbor hydrogen and each of the heavy atoms in the molecule. For convenience, we will denote the two oxygen atoms of each parent by their state at the beginning of the trajectories, i.e., the carbonyl oxygen and the hydroxyl oxygen. Additionally, to test for intermolecular chemistry, we monitored the distances between the carbon atom of the parent and its nearest neighbor carbon atom or the nearest-neighbor carboxyl and hydroxyl oxygen atoms from adjacent molecules. Finally, we monitored the distances between the hydroxyl and carbonyl oxygen atoms with nearest-neighbor hydroxyl or carbonyl oxygen atoms on adjacent molecules. It is challenging to define chemical species in a system that is in a hot dense state using interatomic distances as the criteria to define bonding, as highly excited molecules demonstrate large vibrational amplitudes of the bonds.6 In some cases, it was not clear if the bond distance trace reflected dissociation followed by a swift recombination (perhaps as a result of being trapped) or whether the trace reflected an extremely hot vibration. In other cases, dissociation was clear, as the departing atom would move progressively farther from its original bonding partner. Of the 314 molecules monitored between the “hot” and “cold” trajectories, all but three reacted. Of the 311 molecules that reacted, 306 molecules experienced cleavage of at least one of the original covalent bonds of the parent. There were 305 of the total 306 dissociating parent molecules that experienced cleavage of the O−H bond; 26 had other bonds (i.e., C−O

at 1550 K, 30 GPa. Although the NVE-BOMD simulations were simulated at a slightly larger constant volume (and corresponding lower pressure) than the average value from NPT-BOMD simulation from which the initial configuration was taken, the average energy difference is ∼0.02 eV/molecule. Also, chemical reactivity of the system as described below continued in this trajectory. The theoretical results are illustrated in Figure 1 for comparison with experimental information. Note that the two quantum data points are almost superimposed. Both frames also show the UHC for liquids for the classical MD results and the Manner et al. data4 to provide a convenient visual comparison of the classical MD and experimental data. The UHC for the Trunin data1 is not shown since it differs only slightly from the Manner et al.4 curve. As reported in ref 4, the Manner et al.4 data follow the UHC, while the Trunin et al.1 data deviate from this curve at pressures greater than 15 GPa. The classical and quantum MD densities of the liquid under room conditions are 1.42 and 1.17 g/cm3, respectively, whereas the experimental value is 1.22 g/cm3. In Figure 1a, the classical MD simulations are in fair agreement with both the Manner et al. results4 and the lower pressure data of Trunin et al.1 (up ∼ 2.5 km/s, corresponding to P < 15 GPa). Above that pressure, however, the deviation of the classical MD results from the Trunin et al.1 results increases with increasing pressure. The differences in the classical MD results from experiment are more pronounced in Figure 1b, showing that the simulated classical system is less compressible than the experimental systems under shock. The classical MD results are in closer agreement with the Manner et al.4 data, and are in good agreement with the UHC calculated using the model’s sound speed and ambient state density. Linear fits to the Manner et al.4 and MD UHC curves in Figure 1a produce slopes of 1.66 and 1.67, respectively, with corresponding intercepts of 1.57 and 1.79 km/s, respectively, indicating that the UHCs are only shifted by a constant. Linear fits to the Manner et al.4 and MD data points in Figure 1a are 1.73 and 1.66, respectively, with corresponding intercepts of 1.45 and 1.77 km/s, respectively. These fits highlight the approximate equidistance of the Manner et al.4 data in Figure 1a with the explicitly nonreactive D

DOI: 10.1021/acs.jpcb.5b08845 J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

Figure 2. Interatomic distances versus time during a trajectory in which extensive hydrogen migration and reaction occurs. (a) Interatomic distances between atoms pairs comprising the covalent bonds in the parent formic acid molecule. (b) Interatomic distances between each atom in the parent molecule with its nearest neighbor. The atom on the parent is the left-most partner of each pair given in the legend. (c) Interatomic distance of the parent hydroxyl oxygen with its nearest neighbor hydrogen. The identity of the parent of the nearest neighbor hydrogen is given by a color and index in the legend. (d) Same as part c except it shows the interatomic distance of the parent carbonyl oxygen with its nearest neighbor hydrogen.

near neighbor hydrogens. Taken together, it is clear that, as the original hydrogen departs from the parent at ∼600 fs, another hydrogen approaches the carboxyl oxygen to form a new OH bond, and the molecule remains “formic acid” (Figure 2b). At ∼1750 fs, a hydrogen exchange reaction occurs in which a new hydrogen replaces the hydrogen bonded to the carbonyl atom (Figure 2d). At 1950 fs, the original hydrogen returns to the hydroxyl oxygen, and protonated formic acid is formed and persists for ∼500 fs. The hydrogen on the carboxyl atom departs, and a hydrogen-atom exchange reaction occurs for the hydrogen bonded to the hydroxyl oxygen. Finally, at 3000 fs, the protonated formic acid is reformed with yet a different hydrogen atom adding to the carbonyl oxygen. We frequently observed all of these behaviors in the reactions involving hydrogen chemistry. The excessive hydrogen-atom transfer between oxygen atoms is not surprising, as this is a weak acid (at RTP) in a hot, dense environment. However, this is not the totality of chemical events in the system, as we see formation of “long-lived” protonated formic acid, formate moieties, and exotic complexes formed after H−C and O−C bond scission. The five reactive events that did not result in cleavage of any parent bond all resulted in the formation of protonated formic acid through hydrogen addition to the carbonyl oxygen atom. To illustrate, Figure 3 shows the time history of bonds between the parent carbon with its nearest neighbor oxygen and hydrogen atoms, and each of the parent oxygen atoms with its nearest hydrogen neighbor. Because there was no hydrogen exchange reaction, the hydroxyl O−H bond is the original O−H of the parent. At ∼750 fs, a second hydrogen adds to the carbonyl oxygen, forming protonated formic acid, which remains until the end of

and/or C−H) broken during the trajectories. The O−H dissociation was often preceded or accompanied by the addition of a hydrogen atom to the parent molecule. We observed that where hydrogen addition and dissociation occur simultaneously, the incoming hydrogen either binds with the carbonyl atom as the O−H bond cleaves, or a hydrogen exchange reaction occurs with the original hydroxyl oxygen. For those cases in which dissociation is preceded by hydrogen addition, a protonated formic acid is formed via addition to the carbonyl oxygen. The protonated formic acid can now (a) reform formic acid via hydrogen loss from either the carbonyl or the hydroxyl oxygen, (b) form a formate moiety via hydrogen loss from both oxygen atoms, or (c) undergo additional hydrogen exchange reactions with approaching hydrogen atoms. In these simulations, the hydrogen atoms are extremely mobile and reactive. This is in direct contrast to the RTP simulation, in which no reactive events occurred for the duration of the trajectory (2.1 ps). A typical trajectory showing significant hydrogen exchange is shown in Figure 2. Figure 2a shows the time history of the original bonds in the parent molecule; clearly, the O−H bond dissociates at ∼600 fs, recombines at ∼1950 fs, and then dissociates again around 2500 fs. Figure 2b is similar to part a, except the bonds that are shown correspond to the parent carbon with its nearest neighbor oxygen and hydrogen atoms (both defined at the start of the trajectory as a hydroxyl or carbonyl oxygen) and each of the parent oxygen atoms with its corresponding nearest hydrogen neighbor. Figure 2c,d shows the same time histories for the near neighbor hydrogens to the parent hydroxyl and carbonyl oxygen atoms, respectively, that are shown in Figure 2b, except that the colors of the points denote the parent molecules of the E

DOI: 10.1021/acs.jpcb.5b08845 J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

vibrational periods. The system described in Figure 4a,b is used to illustrate this. In Figure 4b, the C−C interaction is stable for the duration of the trajectory and has an average value of 1.538 Å. Simple transient close encounters are assumed to be minimally interacting. We found that Figure 4a,b is not sufficient to characterize the molecular species due to multiple neighbor interactions not captured in these plots, and therefore, we isolated the final complex for visual analysis (Figure 5), using a cutoff distance of 1.6 Å for determining the terminal ends of the complex. A similarly complex heavy-atom reactive event was also found after examining the time traces for a different molecule (Figure 6). In Figure 6a, the bond between the C and the hydroxyl group is highly excited at the beginning of the trajectory, and is the first to break, followed shortly thereafter by cleavage of the C−H bond, leaving CO which remained as an isolated molecule until 800 fs after which a stray proton from a different parent was added. Figure 6b shows the interatomic distance between this proton and the carbon of the parent atom in Figure 6a. At the end of the simulation, this HCO moiety complexed with another formic acid via its carbonyl oxygen to form the structure shown in Figure 7. Meanwhile, the departing OH from the original molecule undergoes a hydrogen exchange reaction shortly before 200 fs and eventually forms an extended complex (shown in Figure 8) involving atoms from eight parent molecules. The hydrogen from the original C−H cleavage eventually becomes part of a recombined formic acid molecule (not shown). These complexation reactions were not isolated examples observed in our examination of the atomic behavior of various moieties throughout the trajectories. Over both the “hot” and “cold” trajectories, we observed multiple pair interactions that are suggestive of complexation. We observed two instances of complexation formation via C−C pair interactions (“hot” trajectory only), 62 instances via C−hydroxyl O interactions (43 for the “cold” and 19 for the “hot”), and 78 instances via C−carbonyl O interactions (49 for the “cold” and 29 for the “hot”). These interactions are similar to those seen in experimental2,27 and theoretical studies2 of formic acid under high pressure static compression, where both papers observed hydrogen bond symmetrization at pressures of ∼20 GPa and polymerization occurring at ∼40 GPa, although those studies were performed at room temperature.

Figure 3. Interatomic distances versus time during a trajectory between each atom in the parent molecule with its nearest neighbor as in Figure 2b.

the simulation. No additional hydrogen exchange reactions occur after addition of the second hydrogen. A few heavy-atom chemical reactions were observed, and the chemistry is quite complex, as evident in the time traces shown in Figure 4a,b. Considering only Figure 4a alone, one would conclude formation of CO2 after cleavage of the C−H bond early in the trajectory followed by elimination of the remaining hydrogen atom at ∼600 fs. However, examination of the distance of the parent carbon atom with its nearest neighbor carbon (Figure 4b) over the course of the trajectory shows that the CO2 is interacting with another molecule via a C−C bond shortly after scission of the C−H bond. Following the time trace of the second molecule, we observed the elimination of an OH at ∼200 fs, tens of femtoseconds after complexation with the first fragment (not shown). After that, extensive additional complexation and dissociation reactions between the first CO2-like fragment with other fragments and free hydrogens occur for the remainder of the trajectory, with the complex at the end of the trajectory shown in Figure 5. This complex is composed of atoms from five parent molecules. We have defined a complexation reaction as one in which an interatomic distance between atoms on molecular fragments is no larger than 1.6 Å on average, and is stable over three

Figure 4. Interatomic distances versus time for (a) atoms pairs comprising the covalent bonds in the parent formic acid molecule and (b) between the C atom in the parent molecule and its near neighbor C or O atom. The left-most atom label in the legend is from the parent molecule, while the right-most indicates the nearest neighbor atom. F

DOI: 10.1021/acs.jpcb.5b08845 J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

Figure 5. Representation of complex at the end of the trajectory described in Figure 4.

Figure 6. Interatomic distances versus time for (a) atoms pairs comprising the covalent bonds in the parent formic acid molecule and (b) those between the C atom in the parent molecule and its nearest H neighbor (not including the original C−H hydrogen).

As stated earlier, only 26 molecules displayed CO or CH scission reactions (for both the “cold” and “hot” trajectories), while there were 142 instances in which complexes were formed. This indicates that the chemistry at these conditions is not direct formation of simple products (e.g., CO, CO2, H2O, etc.) but rather formation of exotic polyatomic species or potentially early stages of polymerization. Some of these species are transient in nature over the course of the simulation and highly reactive, and the lengths of the simulations preclude us from drawing any conclusions about steady state chemical equilibrium. Further, simple time traces such as shown in Figures 4 and 6 are not sufficient to fully capture all of the chemistry of this system due to interactions with multiple neighbors. These examples clearly demonstrate the problematic nature of identifying species in a hot dense reacting environment using interatomic distance criteria.6

Figure 7. Representation of the first complex at the end of the trajectory described in Figure 6.

Figure 8. Representation of the second complex at the end of the trajectory described in Figure 6. G

DOI: 10.1021/acs.jpcb.5b08845 J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

4. SUMMARY AND CONCLUSIONS The original purpose of this study was to address a discrepancy between experimental studies of shocked formic acid. The Trunin et al. study1 using the reflection method produced a shock Hugoniot with a feature that suggested chemical reactivity above 15 GPa; a second study by Manner et al.4 using gas gun-driven plate impact experiments produced a shock Hugoniot that indicated no reaction for shock conditions above 15 GPa. We performed BOMD simulations to generate shock Hugoniot points at 30 GPa (well above the pressure corresponding to the discontinuity in the shock Hugoniot of Trunin et al.1) and compared these to results from experiment1,4 and classical nonreactive MD simulations.9,10 The behavior of the Manner et al.4 data was in better agreement with the nonreactive classical MD shock Hugoniot curve. The reactive QM Hugoniot points at 30 GPa were in closer agreement with the Trunin et al.1 data, which showed a distinct discontinuity in the shock Hugoniot curve at ∼15 GPa. We can only speculate that differences in the experiments produce the reactive versus nonreactive shock response. Interrogation of BOMD simulations at 1550 K, 30 GPa, conditions close to the shock Hugoniot points, indicated extensive chemical reaction. The vast majority of reactions involve hydrogen-atom migration. There were fewer, but still significant, instances of heavy-atom chemistry. Most species were transient; however, there existed several examples of longlived protonated and molecular fragment species (e.g., HCOO) and intermolecular complexation involving heavy atoms from multiple parents. However, the majority of the species at any one moment in the trajectories remained either formic acid or protonated formic acid, although they were not the parent molecule, but rather recombinant from hydrogen migration events. Determination of the true identity of species based on bond distance analysis is difficult due to the hot, dense state of the material. Additionally, the small sample size (157 molecules) limits any statistical analysis of species and reactions. We are confident in stating that chemistry is occurring, that nonhydrogen chemistry is occurring, but anything more conclusive than that requires larger and longer simulations, as well as more formal electron-density-based means of conclusively determining chemical bonds versus stabilization of hot vibrational modes due to interactions with neighbor molecules. While we clearly observed chemistry, we do not know whether or not we have reached chemical equilibrium in such a limited simulation time. We can only state that we have reached a thermodynamic equilibrium on the time scale of our simulation. The question arises as to how the thermodynamic equilibrium will change as chemical equilibrium is achieved, thus potentially influencing the final shock Hugoniot value. If chemical equilibrium has not been reached, then the shock Hugoniot value for the system at chemical equilibrium could be different. Thus, it is possible that the reasonable agreement of the BOMD predictions of the shock Hugoniot with the Trunin et al.1 data at 30 GPa is fortuitous. Regardless, this study has shown initial chemical reactions at conditions near that of shocked formic acid at 30 GPa. This should spur further experimental explorations targeting smaller temporal and spatial regimes in order to better identify early-stage chemistry. It should also motivate development of novel quantum mechanical methods to model reactive events on suitably long time scales, such as the XL-BOMD method,7,8 for use in

determining the EOS for reactive events, and methods to unambiguously identify transient species and chemical bonding in systems at extreme conditions.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.5b08845. Thermodynamic properties from classical and quantum MD simulations (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: (410) 306-1904. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was partially supported by the DoD High Performance Computing Modernization Program Software Application Institute for Multi-Scale Reactive Modeling of Insensitive Munitions (2008-2014).

■

REFERENCES

(1) Trunin, R. F.; Zhernokletov, M. V.; Kuznetsov, N. F.; Radchenko, O. A.; Sychevskaya, V. V.; Shutov, V. V. Shock Wave Compression Of Liquid Organic Compounds. Sov. J. Chem. Phys. 1992, 11, 606−619. (2) Goncharov, A. F.; Manaa, M. R.; Zaug, J. M.; Gee, R. H.; Fried, L. E.; Montgomery, W. B. Polymerization Of Formic Acid Under High Pressure. Phys. Rev. Lett. 2005, 94, 065505−065505-4. (3) Montgomery, W.; Zaug, J. M.; Howard, W. M.; Goncharov, A. F.; Crowhurst, J. C.; Jeanloz, R. Melting Curve And High-Pressure Chemistry Of Formic Acid To 8 GPa And 600 K. J. Phys. Chem. B 2005, 109, 19443−19447. (4) Manner, V. W.; Sheffield, S. A.; Dattelbaum, D. M.; Stahl, D. B. Shock Compression of Formic Acid. In Shock Compression of Condensed Matter-2011; Elert, M. L., Butler, W. T., Borg, J. P., Jordan, J. L., Vogler, T. J., Eds.; AIP Conference Proceedings; American Institute of Physics: Melville, NY, 2012; Vol. 1426, pp 201− 204. (5) Taylor, D. E.; Rice, B. M. Quantum-Informed Multiscale M&S For Energetic Materials. Adv. Quantum Chem. 2014, 69, 171−219. (6) Ge, N.-N.; Wei, Y.-K.; Zhao, F.; Chen, X.-R.; Ji, G.-F. PressureInduced Metallization Of Condensed Phase B-HMX Under Shock Loadings Via Molecular Dynamics Simulations In Conjunction With Multi-Scale Shock Technique. J. Mol. Model. 2014, 20, 2350. (7) Niklasson, A. M. N.; Tymczak, C. J.; Challacombe, M. TimeReversible Born-Oppenheimer Molecular Dynamics. Phys. Rev. Lett. 2006, 97, 123001−123001-4. (8) Niklasson, A. M. N.; Cawkwell, M. J. Fast Method For Quantum Mechanical Molecular Dynamics. Phys. Rev. B: Condens. Matter Mater. Phys. 2012, 86, 174308−174308-10. (9) Baker, E. L.; Stiel, L. I.; Capellos, C.; Rice, B. M.; Bunte, S. W.; Byrd, E. F. C. Formic Acid Investigation For The Prediction Of High Explosive Detonation Properties And Performance. Int. J. Energ. Mater. Chem. Propul. 2010, 9, 377−384. (10) Baker, E. L.; Stiel, L. I.; Capellos, C.; Rice, B. M.; Bunte, S. W.; Byrd, E. F. C. Formic Acid Effects On The Prediction Of High Explosive Detonation Properties. Proc. of the 14th Int. Detonation Symp. 2010, 1079−1086. (11) Erpenbeck, J. J. Molecular Dynamics Of Detonation. I. Equation Of State And Hugoniot Curve For A Simple Reactive Fluid. Phys. Rev. A: At., Mol., Opt. Phys. 1992, 46, 6406−6416. (12) Roszak, S.; Gee, R. H.; Balasubramanian, K.; Fried, L. E. New Theoretical Insight Into The Interactions And Properties Of Formic H

DOI: 10.1021/acs.jpcb.5b08845 J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

Acid: Development Of A Quantum-Based Pair Potential For Formic Acid. J. Chem. Phys. 2005, 123, 144702−144702−10. (13) Smith, W.; Forester, T. R. J. DL_POLY_2.0: A General-Purpose Parallel Molecular Dynamics Simulation Package, v.2.20. J. Mol. Graphics 1996, 14, 136−14, accessed 2011. (14) Allen, M. P.; Tindesley, D. J. Computer Simulation of Liquids; Oxford University Press: New York, 1989. (15) Woolfolk, R. W.; Cowperthwaite, M.; Shaw, R. A ‘Universal’ Hugoniot For Liquids. Thermochim. Acta 1973, 5, 409−414. (16) Sheffield, S. A.; Alcon, R. R. Shock-Induced Reactions in Several Liquids. In Shock Compression of Condensed Matter-1989; Schmidt, S. C., Johnson, J. N., Davison, L. W., Eds.; Elsevier Science Publishers B.V.: New York, 1990; pp 683−686. (17) Hooper, J. B.; Bedrov, D.; Smith, G. D.; Hanson, B.; Borodin, O.; Dattelbaum, D. M.; Kober, E. M. A Molecular Dynamics Simulation Study Of The Pressure-Volume-Temperature Behavior Of Polymers Under High Pressure. J. Chem. Phys. 2009, 130, 144904− 144904−11. (18) http://www.cp2k.org/ for CP2K Open Source Molecular Dynamics, v 2.0.0, accessed Dec 5, 2012. (19) VandeVondele, J.; Krack, J.; Mohamed, F.; Parrinello, M.; Chassaing, T.; Hutter, J. Quickstep: Fast And Accurate Density Functional Calculations Using A Mixed Gaussian And Plane Waves Approach. Comput. Phys. Commun. 2005, 167, 103−128. (20) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868; Phys. Rev. Lett. 1997, 78, 1396 (erratum). (21) Perdew, J. P.; Ruzsinszky, A.; Csonka, G. I.; Vydrov, O. A.; Scuseria, G. E.; Constantin, L. A.; Zhou, X.; Burke, K. Restoring The Density-Gradient Expansion For Exchange In Solids And Surfaces. Phys. Rev. Lett. 2008, 100, 136406−136406-4; Phys. Rev. Lett. 2009, 102, 039902 (erratum). (22) Lin, I.-C.; Coutinho-Neto, M. D.; Felsenheimer, C.; von Lilienfeld, O. A.; Tavernelli, I.; Rothlisberger, U. Library Of Dispersion-Corrected Atom-Centered Potentials For Generalized Gradient Approximation Functionals: Elements H, C, N, O, He, Ne, Ar, And Kr. Phys. Rev. B: Condens. Matter Mater. Phys. 2007, 75, 205131−205131-5. (23) Balu, R.; Byrd, E. F. C.; Rice, B. M. Assessment of Dispersion Corrected Atom Centered Pseudopotentials: Application to Energetic Molecular Crystals. J. Phys. Chem. B 2011, 115, 803−810. (24) Sorescu, D. C.; Byrd, E. F. C.; Rice, B. M.; Jordan, K. D. Assessing the Performances of Dispersion-Corrected Density Functional Methods for Predicting the Crystallographic Properties of High Nitrogen Energetic Salts. J. Chem. Theory Comput. 2014, 10, 4982− 4994. (25) VandeVondele, J.; Hutter, J. An Efficient Orbital Transformation Method For Electronic Structure Calculations. J. Chem. Phys. 2003, 118, 4365−4369. (26) Weber, V.; VandeVondele, J.; Hutter, J.; Niklasson, A. M. N. Direct Energy Functional Minimization Under Orthogonality Constraints. J. Chem. Phys. 2008, 128, 084113−084113−9. (27) Manner, V. W.; Chellappa, R. S.; Sheffield, S. A.; Liu, Z.; Dattelbaum, D. M. High-Pressure Far-Infrared Spectroscopic Studies Of Hydrogen Bonding In Formic Acid. Appl. Spectrosc. 2013, 67, 1080−1086.

I

DOI: 10.1021/acs.jpcb.5b08845 J. Phys. Chem. B XXXX, XXX, XXX−XXX