Article pubs.acs.org/JPCC
Dissociation of Sarin on a Cement Analogue Surface: Effects of Humidity and Confined Geometry Christopher J. O’Brien,*,† Jeffery A. Greathouse,‡ and Craig M. Tenney‡ †
Department of Computational Materials and Data Sciences and ‡Department of Geochemistry, Sandia National Laboratory, Albuquerque, New Mexico 87185, United States ABSTRACT: First-principles molecular dynamics simulations were used to investigate the dissociation of sarin (GB) on the calcium silicate hydrate (CSH) mineral tobermorite (TBM), a surrogate for cement. CSH minerals (including TBM) and amorphous materials of similar composition are the major components of Portland cement, the binding agent of concrete. Metadynamics simulations were used to investigate the effect of the TBM surface and confinement in a microscale pore on the mechanism and free energy of dissociation of GB. Our results indicate that both the adsorption site and the humidity of the local environment significantly affect the sarin dissociation energy. In particular, sarin dissociation in a low-water environment occurs via a dealkylation mechanism, which is consistent with previous experimental studies.
1. INTRODUCTION Restoration of critical civilian or military infrastructure after a chemical warfare agent (CWA) attack will almost certainly involve decontamination of concrete. Developing a more detailed understanding of how a CWA binds and dissociates within concrete will allow the development of more efficient decontamination strategies. This approach extends prior modeling efforts from highly idealized systems to more realistic surfaces and materials1−4 while still keeping complexity sufficiently low to investigate fundamental properties of adsorption and dissociation. As a fundamental understanding of CWA behavior on model surfaces is gained, that knowledge will be used to guide the development of more complex pore models, which can be used to further probe CWA fate and transport in varying environmental scenarios. A principle motivation of this work is to determine the effect of cement analogue surfaces on the binding and dissociation of isopropyl methylphosphonofluoridate, commonly known as sarin or GB. The structure of GB is provided as the first reactant in Scheme 1. Composed of aggregate and cement, concrete is chemically and physically heterogeneous across multiple length scales. While individual aggregate particles may be relatively nonporous and uniform, cement is composed of several mineral species and has pores that vary in size from nanometers to micrometers. Prior experiments suggest the cement phase, which also contains a significant amount of chemically and physically bound water, primarily determines the transport and fate of chemical species within concrete.5−7 The most commonly employed cement is Portland cement, which itself has a complex mineralogy.8 It includes phases such as calcium-silicate-hydrates (CSH), calcium hydroxides, calcium carbonates, and alumina-ferric-oxide sulfates. Exposed pore surface chemistry will vary with initial cement composition, age,9 and environment. Realistic molecular models of cement, © 2016 American Chemical Society
calcium-silicate-hydrate, oxide, and hydroxide mineral phases have been developed10−12 and used to study adsorption of water10,13 and docking of radionuclides.14 A common feature of cement models (see Chapter 5.7 of ref 8) is the presence of the calcium silicate hydrate mineral tobermorite (TBM), Ca2.25[Si3O7.5(OH)1.5]·H2O.12,15 More realistic models include deviation from the ideal stoichiometry,11 the presence of broken Si−O chains,16 or details of mesoscale pore structure.17 In this study, a microporous model of TBM is used as a surrogate for cement. 1.1. Adsorption of Sarin to Surfaces. Prior research on the adsorption or degradation of GB (or simulant) has mostly involved empirical measurement and models of interactions with bulk materials. Because of the chemical and physical complexity of materials such as concrete, the primary mechanisms controlling the rate of decomposition remain unknown. Molecular and mesoscale simulations are ideal for predicting fundamental interactions between chemical warfare agent (CWA) molecules, decontamination materials, and the interior surfaces of porous materials. Classical molecular dynamics (MD) and first-principles molecular dynamics (FPMD)18 models have been used to examine the behavior of CWAs in water and near surfaces. The mechanisms of adsorption are best understood on inorganic oxides (e.g., MgO,4,19 γ-Al2O3,20−23 aluminosilicates,2 and silica1,24,25). Metal oxides are of particular interest as protective coatings because of their ability to destroy GB upon adsorption. GB binds most strongly to metals or hydroxide groups through the phosphoryl oxygen (PGBOGB). Adsorption preferentially Received: October 4, 2016 Revised: November 22, 2016 Published: November 22, 2016 28100
DOI: 10.1021/acs.jpcc.6b10046 J. Phys. Chem. C 2016, 120, 28100−28109
Article
The Journal of Physical Chemistry C
Scheme 1. SN2 Reaction Scheme Illustrating the Reaction between GB and H2O, Forming an Activated State, That Decomposes via the Departure of FGB
Scheme 2. SN1 Reaction Scheme Illustrating a Two-Step Mechanism for the Decomposition of GB in the Presence of Watera
a
The key difference with that of the SN2 scheme is that the binding of water does not occur before FGB leaves.
GB, is dissociation via the SN1 mechanism (Scheme 2). The SN1 mechanism proceeds via replacement of the FGB by a water molecule, which then gives up a proton. The hydrolysis reaction occurs in neutral water but is greatly accelerated in basic solutions.28 This is particularly relevant to decomposition of GB in cement as it is a highly alkaline environment. Experiments have also demonstrated the existence of a dealkylation mechanism (removal of the propyl group) as illustrated in Scheme 3. Dealkylation was found to be very slow on alumina but could be accelerated by removal of hydroxyl groups, which results in the presence of acidic adsorption sites.21
occurs to basic surface sites, when available, such as nonhydroxylated Al in alumina or Mg in MgO.21 Of particular interest to this work are interactions of GB with minerals similar to those present in cement. The adsorption of GB and soman (GD) in the absence of water to a model aluminosilicate surface was examined by Michalkova et al.,2 who demonstrated the formation of multiple hydrogen bonds, particularly between the propyl or methyl groups and the surface O atoms. These simulations also showed an interaction between F GB and surface OH groups, although such configurations were less strongly adsorbed than those in which the phosphoryl oxygen was involved. Related calculations of adsorption of GB and GD to [MO(OH)n]n−5 (for n = 3 or 4 and M = Si or Al) fragments found multiple hydrogen bonds forming between the propyl and methyl groups and the cluster oxygen atoms, even if P was chemically bound to an unterminated O in the fragment.3 On nonplanar surfaces, these studies demonstrate that hydrogen bonding plays a major role in adsorption beyond the principle hydrogen bond between the phosphoryl oxygen and surface OH groups. A detailed ab initio study by Alam et al.26 suggests that watermediated hydrogen bonding interactions between GB and inorganic surfaces are actually more energetically favorable than direct bonding of GB to inorganic surfaces. Such a prediction is difficult to verify as most experimental GB surface interaction studies (including those discussed herein) are conducted under high vacuum. In such a case, surfaces are typically hydroxylterminated or denuded of water. Also important is the potential effect of simulants versus actual agents. An experimental study of GB adsorbing to γAl2O3 found that GB forms strong hydrogen bonds with the surface through the phosphoryl oxygen and readily dissociates through hydrolysis,21,22 whereas the simulant dimethyl methylphosphonate (DMMP) only weakly physioadsorbs to the surface and dissociates very slowly,4 or not at all.20 1.2. Stability of Adsorbed Sarin. Simulations and experiments using idealized model systems have shown that CWA degradation pathways and rates can vary significantly depending upon the adsorption or solvation state,27 pH,28 and surface properties (hydrophilic vs hydrophobic),1,7 including water content.7,21 These findings underscore the importance of surface effects in the transport and fate of CWA in porous materials. In water, GB is known to dissociate via an SN2 mechanism that involves water forming an activated state by binding to PGB, as illustrated in Scheme 1.29 As a result, FGB leaves GB and receives a proton from water, forming aqueous HF. An alternative reaction, which also results in FGB leaving
Scheme 3. Dealkylation Reaction Scheme Illustrating the Decomposition of Sarin in the Absence of Water
Calculations and experiments show that GB may undergo destructive adsorption where GB will hydrolyze after bonding to certain surfaces such as MgO4,19 and γ-Al2O3.3,21,22 It is notable that GB does not readily hydrolyze under vacuum conditions on a hydroxylated silica surface.24 The role of water in these reactions is complex. Experimental work with γ-Al2O3 has shown that although the phosphoryl oxygen binds to nonhydroxylated Al sites, hydrolysis of GB would then likely proceed by replacing F with the oxygen of the neighboring OH group.21 The authors partially based this conclusion on the observation that hydroxylated surfaces are not necessary for adsorption of GB but are necessary for decomposition of adsorbed GB. Similarly, an experimental study of GB adsorbed onto cement paste also showed rapid decomposition and acceleration of the reaction with an increase in water content.7 Although bonding to a surface greatly affects the favorability of decomposition of GB, recent work by Kuo et al.1 indicates that even being near a surface can modify the activation energy. Kuo et al. calculated that the activation energy of decomposition is reduced by ∼15 kcal mol−1 when GB is in the presence of, but not bound to, a hydrophilic surface (amorphous SiO2 glass). The authors also determined that if GB is constrained to be near a model hydrophobic surface (a layer of butane molecules) the reaction mechanism changed to SN1 with a resulting increase in reaction energy of ∼10 kcal mol−1 with respect to that in bulk water. 28101
DOI: 10.1021/acs.jpcc.6b10046 J. Phys. Chem. C 2016, 120, 28100−28109
Article
The Journal of Physical Chemistry C
narrow (∼11 Å) slit-pore is generated and filled with water to obtain the appropriate room-temperature density. MD simulations were conducted with the cell size fixed to match the lattice constants calculated with first-principles density functional theory at 300 K to prevent unreasonable strain in the surface when performing FPMD calculations. 2.2. Classical and First-Principles Molecular Dynamics. The TBM pore, water molecules, and GB were modeled with interatomic potentials taken from ClayFF,31 flexible SPC water,32 and GAFF.33 The PACKMOL34 code was used to generate initial configurations of GB with varying amounts of water in the TBM slit-pore. The PACKMOL34 code and ClayFF interatomic potential31 were used to generate initial configurations of GB with varying amounts of water in the TBM slit-pore. Water was treated as an explicit solvent. Classical MD calculations were performed with the LAMMPS package,35 while visualization and figure generation were performed utilizing VMD.36 The CP2K37 code was used to examine the dissociation energy and mechanism in water and when attached to various binding sites at the surface. The mixed Gaussian and planewave method38 was employed using the Becke39 exchange and Lee−Yang−Parr40 correlation functionals. The cutoff for plane waves representing the density was set to 300 Ry, with a CP2Kspecific relative cutoff of 60 Ry. Plane waves were calculated only at the Γ point in reciprocal space. The wave functions were calculated in the unrestricted Kohn−Sham scheme, and self-consistency was defined with a threshold of 1 × 10−6 Hartree using the Orbital Transform41 scheme to determine the ground state. MOLOPT basis sets of triple-ζ quality with two polarization functions42 were employed, whereas ionic core states were modeled using the Goedecker−Teter−Hutter pseudopotentials.43 The QuickStep algorithm44 for performing FPMD, within the CP2K code, is employed to determine the dissociation energy of GB in water and attached to various sites on the surface. MD was performed using a time step of 0.48 fs. All runs were performed at 300 K under NVT conditions employing a Langevin thermostat. To better account for dispersion interactions, critically important to systems containing organic molecules, the DFT-D2 dispersion correction was employed as developed by Grimme et al.45 Determining the conformation with which GB binds to TBM is crucial to understanding not only the mechanism but also the energetics of dissociation. The mechanism and energetics of the dissociation reactions were calculated using metadynamics.46,47 Metadynamics is a powerful technique for efficiently exploring a free energy surface (FES) while performing a FPMD run. Metadynamics involves the use of a predefined parameter, called a CV (collective variable), which is manipulated to force a reaction to proceed in the desired manner and monitored to measure how a reaction progresses. Examples of CVs include coordination, interatomic distance, torsion, and bond angle deviation. One or more CVs may be employed in a model. This involves an intelligent choice of the best ways to characterize the reaction of interest but does not require specification of the exact path of the reaction a priori. Literature suggests that phosphororganics, including GB, degrade via the leaving of an O or F group. In this work, the choice of CVs was made in a manner consistent with the work of Kuo et al.1 Three CVs were chosen: (c1) FGB−PGB and (c2) PGB−OGB interatomic distances and (c3) the coordination of water about PGB. The coordination of water about PGB is defined in ref 1 as
This study takes GB both bound to and in the proximity of the surface of hydrophilic TBM and determines if the nature of the surface interaction affects the mechanism or energetics of decomposition.
2. METHODOLOGY Molecular dynamics was used to survey the preferred sites of binding of GB to TBM and obtain a qualitative sense of the relative stability of the sites. Once an understanding of the binding sites was obtained, FPMD (also known as Ab Initio MD18) was utilized to calculate the dissociation energy of solvated and adsorbed GB and, with the aid of metadyamics,30 the reaction path. FPMD calculations were undertaken, beginning with configurations generated by classical MD, so that a better treatment of the energetics of bond breaking would be possible. Both FPMD and classical MD techniques describe the motion of atomic nuclei using Newtonian mechanics but differ in the level of theory used in calculating the interatomic forces. In FPMD, the forces originate from explicit calculation of the electronic structure of the system, typically employing density functional or Hartree−Fock theory, whereas classical MD utilizes parametrized descriptions of the interatomic interactions. 2.1. Generation of Tobermorite Surfaces. The orthorhombic representation of 11 Å TBM by Murray et al.,12 illustrated in Figure 1a, was employed rather than the originally
Figure 1. Unit cell of tobermorite (TBM) from (a) ref 12 at 0 K and (b) the slit-pore used as a model surface at 300 K. Dark lines indicate the boundaries of the periodic unit cell that were chosen for generation of a pore approximately 11 Å in width. Surfaces were chosen to contain Ca ions with interstitial water molecules near the surface. Ca ions are colored green, O atoms red, H atoms white, and Si atoms gold.
characterized monoclinic form.15 Because of the large size of the cell and the high cost of FPMD calculations, it is desirable to find a way to reduce the cell size that preserves the unique surfaces that could be exhibited by the TBM crystal. First, the roto-inversion symmetry about the z axis is exploited to cut the sample in half at the low-bond density x−y plane between layers of Ca ions. A likely environmentally stable surface, illustrated in Figure 1b, was constructed by cutting at a {001} plane separating layers of interstitial waters and Ca. MD simulations were conducted in a periodic simulation cell with the dimensions of a single uncut TBM unit cell. The Ca/water surfaces wrapped around the periodic boundary so that a 28102
DOI: 10.1021/acs.jpcc.6b10046 J. Phys. Chem. C 2016, 120, 28100−28109
Article
The Journal of Physical Chemistry C
c3
rij 6
() = ∑ 1−( ) 1−
NH2O
3.1. Binding of Sarin to Tobermorite. Classical MD was used to sample a large number of configurations for the binding of GB to TBM surfaces. Because FPMD calculations are very costly, they can be run for only very short times and will not sample the configurations that MD models can by modeling a much longer period of time. For the systems studied, FPMD calculations would cover a few picoseconds using similar levels of resources it took MD models to progress tens of nanoseconds. Multiple starting configurations and long simulation times were used to qualitatively determine favorable binding sites. In this work, the RP chiral form is employed but does not change the conclusions presented because of the lack of chiral binding sites on the TBM surface. It is observed that GB moves between the more stable binding sites, even at room temperature, as GB undergoes conformational changes. It must be noted that GAFF fixes partial charges with O set at −0.843 e− and F set at −0.341 e−. The Ca ion at the surface of TBM is treated as an aqueous ion, with a partial charge of +2.0 e−. In contrast, DFT calculations yield a Mulliken charge of −0.63 on OGB, −0.42 on FGB, and +1.2 on the surface Ca ions. These differences have a substantial impact on the stability of binding of GB to TBM surfaces. Radial distribution functions provide at least a relative measure of the GB−TBM binding strength that can be inferred by the presence of a sharp peak in the distribution. An increased degree of correlation at a particular distance indicates a stronger bond due to reduced displacement of the atom about that point. Figure 2, generated from MD
r0
rij 12 r0
(1)
where r0 = 1.59 Å and rij is the distance between PGB and an O in a water molecule. Certain cases that required that c2 be changed to monitor the OGB−CGB bond distance arose. Our results indicate that the PGB−OGB bond length does not vary significantly during reactions studied herein, and the possibility of dealkylation (via leaving of the propanol group) is more likely in many cases. The free energy surface was sampled by adding local repulsive potentials (hills) at the present coordinates of CVs in phase space to force the system to sample other regions of phase space. Hills of Gaussian shape with a height of 3 × 10−3 H and a width of 0.10 were added every 15 MD time steps. Metadynamics has been previously employed to study the decomposition of VX {the nerve agent O-ethyl S-[2(diisopropylamino)ethyl] methylphosphonothioate} with the hypochlorite ion by Gee et al.27 and the hydrolysis of GB in water on hydrophilic and hydrophobic surfaces by Kuo et al.1 A limitation of metadynamics is that it is known to underestimate the reaction barrier of SN2 reactions, such as the hydrolysis of GB in water.48 Various sources of error for this technique have been determined, including the observation that GGA techniques are known to underestimate the activation energy,48,49 and the limited duration of FPMD calculations may prevent products from becoming properly coordinated with solvent.48 A FPMD equilibration was conducted under NVT conditions for a minimum of 2 ps before beginning metadynamics calculations to allow bond angles and distances to adjust to those resulting from the FPMD calculations. The free energy surface (FES) generated by metadynamics was plotted as a contour surface for visualization. For analysis, the FES is represented on a grid with divisions of 0.0106 Å. Energy minima in the FES were used to assign product and reactant states that determined the end points of the reaction path. The FES and minimum energy path (MEP) plotted have undergone thermodynamic integration to determine the free energy change between points rather than showing the minimal potential energy values. Because of the desire to sample a larger number of dissociation reactions rather than determine quantitative values for reaction energy, the reaction barrier was simply estimated as the highest-energy point on the MEP between the reactant and product local minima. More quantitative estimates would require climbing image nudged elastic-band calculations to determine the free energy of the reaction barrier. However, for this initial work, it was desirable to survey a larger number of possible binding states and estimate how the binding and the amount of water affect the dissociation of GB.
Figure 2. Radial distribution functions for FGB−Ca (solid blue line) and OGB−Ca (dashed green line) in a TBM pore with 50 H2O molecules calculated using the MD model that illustrate a strong preference for binding of the phosphoryl oxygen to Ca atoms in the surface.
models, illustrates the stability of the bond between the phosphoryl oxygen and surface Ca, and the comparatively weaker interaction of the FGB−Ca bond in the presence of 50 H2O molecules within a TBM slit-pore. However, as reported in Table 1, the FPMD model did not predict that the OGB−Ca bond would be stable at 300 K, even for an interval of
