Theoretical Study of the Adsorption of Organophosphorous

Subscriber access provided by Washington & Lee University Library and VIVA (Virtual Library of Virginia)

Article

Theoretical Study of the Adsorption of Organophosphorous Compounds to Models of a Silica Surface Diego Troya, Angela C Edwards, and John R. Morris J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/jp404065n • Publication Date (Web): 20 Jun 2013 Downloaded from http://pubs.acs.org on June 22, 2013

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry C is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Theoretical Study of the Adsorption of Organophosphorous Compounds to Models of a Silica Surface Diego Troya,* Angela C. Edwards, and John R. Morris Department of Chemistry, Virginia Tech, 107 Davidson Hall, 24061 Blacksburg, VA *

Corresponding author: [email protected]

Abstract We present a study of the interaction of seven organophosphorous (OP) compounds with models of a silica surface using electronic structure calculations. Two of the OP species are chemical-warfare nerve agents (sarin and soman), and the rest are organophosphate (methyl dichlorophosphate, dimethyl chlorophosphate, trimethyl phosphate)

or

organophosphonate

(dimethyl

methylphosphonate,

diisopropyl

methylphosphonate) mimics utilized in recent experiments. Calculations that approximate an amorphous silica surface by a simple gas-phase silanol molecule qualitatively reproduce the trend in binding energies determined in recent ultrahigh vacuum experiments. These calculations also show that of the various atoms in the OP species that can accept hydrogen bonds from silanol groups, the oxygen atom that forms a double bond with the central phosphorous atom provides the largest binding energy, followed first by the remaining oxygen atoms, and then by the halogen atoms. More quantitative comparisons with experiments are established by using calculations that more rigorously appoximate an amorphous silica surface with silica clusters that contain one or two silanol groups. MP2 calculations of the binding energies extrapolated to the completebasis-set limit in a cluster model exhibiting a single isolated silanol group differ from the recent ultrahigh vacuum experimental values by less than 1 kcal/mol for many of the OP

ACS Paragon Plus Environment

1


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 38

compounds examined. Calculations on a silica cluster with two silanol groups reveal the possibility of formation of two hydrogen bonds with the OP compounds and produce binding energies that are also consistent with other recent measurements. Symmetryadapted perturbation theory is used to provide deeper insight into the origin of the exceptionally strong hydrogen bonds established between silica surfaces and OP compounds.

Keywords: adsorption energy, chemical-warfare agents, silica clusters


2

Page 3 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Introduction The adsorption of gas-phase molecules on surfaces continues to be of interest due to its broad implication in fields ranging from catalysis to sensing. This adsorption is primarily governed by the strength of the intermolecular interaction between the functional groups of the gas-phase species and the surface. Silica is among the most prevalent solids in the environment, and molecular adsorption on its surface has long been of interest.1 The silica surface is characterized by silanol groups that can establish relatively strong hydrogen-bonding interactions with a variety of gas-phase species, including biomolecules.2 The ability of silica surfaces to strongly adsorb molecular species via hydrogen bonds implies that the residence time of the adsorbate on the surface can be very long even at ambient temperatures, and this can represent an environmental hazard depending on the toxicity of the adsorbed species. Such is the case of gas-phase chemical-warfare agents (CWAs), whose adsorption and survivability on common surfaces present a challenge beyond the battlefield. In order to better understand the adsorption of CWAs to silica and design more effective sensing methods3 and mitigation strategies,4 vigorous research has been carried out in recent time at both the experimental and theoretical levels. A well-known class of CWAs is that formed by organophosphonate compounds, which include the nerve agents sarin, soman, tabun, and VX.5 The difficulties that investigating these substances in a laboratory setting entail have drawn attention to other organophosphorous (OP) compounds that mimic the structure and surface chemistry of the nerve agents, but do not possess their toxicity. The mimics contain many of the functional groups of the live agents, including a phosphoryl or phosphonyl group in which a central phosphorous atom


3


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 38

is linked through a covalent double bond to an oxygen atom with formally sp2 hybridization, with additional alkyl, alkoxy, and halogen substituents around that central atom. Five of such mimics for which the binding energies to silica have been measured are methyl dichlorophosphate (MDCP), dimethyl chlorophosphate (DMCP), trimethyl phosphate

(TMP),

dimethyl

methylphosphonate

(DMMP),

and

diisopropyl

methylphosphonate (DIMP). In the work presented below, we describe a computational study of the adsorption of all five compounds, in addition to sarin (propan-2-yl methylphosphonofluoridate) and soman (O-pinacolyl methylphosphonofluoridate) to silica models. All seven molecules have at least two hydrogen-bond acceptors with which they can establish strong interactions with the surface: the sp2 oxygen atom, and the sp3 oxygen atom of the alkoxy substituents. Early detailed studies determined that, unlike common metal-oxide surfaces, silica allowed for the molecular adsorption of DMMP.6 Following that work, measurements of the adsorption of OP compounds on silica surfaces were performed by Kanan and Tripp.7 The binding of MDCP, DMMP, TMP, and trichlorophosphate (TCP) to an amorphous silica thin film pretreated at 673 K was tracked by infrared (IR) spectroscopy. Spectral signatures indicated that the adsorption of the OP compounds was mostly molecular (i.e., dissociation does not take place). Regarding the mechanism of adsorption, the authors concluded that MDCP formed two hydrogen bonds (with both the sp2 and sp3 oxygen atoms acting as acceptors) with two silanol groups on the surface, DMMP also formed two hydrogen bonds with the surface, but both involving the two sp3 oxygen atoms in the adsorbate, and all three sp3 oxygen atoms of TMP participated in hydrogen bonds with the surface. The strength of the interaction, determined from high-


4

Page 5 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


temperature evacuation times and chemical displacers,8 followed the MDCP < DMMP < TMP trend. Subsequent computational studies by Bermudez, however, convincingly showed that for cluster models containing up to 72 atoms of a silica surface with a silanol density of 4 OH nm-2, the largest binding energy corresponded to the formation of two hydrogen bonds between two vicinal silanol groups and the sp2 oxygen atom of the OP compound; binding energies involving the sp3 oxygen atoms appeared to be lower.9 For that computational work, Bermudez used density functional theory (B3LYP) in combination with mostly double-zeta split-valence basis sets with various levels of polarization. The challenge of computing the energy for a molecule as large as DMMP interacting with a silica cluster composed of up to 72 atoms was elegantly solved by applying hybrid ONIOM methods, in which ‘inactive’ regions of the surface that provide only structural support to the adsorption site were treated with Hartree-Fock theory. Further insight into the binding mechanism of DMMP to amorphous silica surfaces was provided by Quenneville et al. using molecular dynamics.10 A particularly valuable aspect of those simulations is that they examined the effect of the silanol surface density on the adsorption of DMMP. This is of special importance because the surface chemistry of amorphous silica is extremely sensitive to the preparation and treatment of the sample: while a fully hydroxylated silica surface at ambient conditions reaches a silanol surface density of nearly 5 OH nm-2,11 thermal treatment of the surface can significantly decrease this value such that annealing at 700 K generates a density of about 2 OH nm-2, and annealing beyond 1000 K causes the density to fall below 1 OH nm-2.12 Furthermore, the surface reactions that take place during thermal treatment also affect the


5


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 38

type of silanol groups at the surface.13,14 Thus, while under ambient conditions vicinal silanols that can establish hydrogen bonds with each other dominate the surface, followed by isolated silanols, and geminal silanols (in which two surface OH groups are coordinated to one Si atom), at 700 K the silica surface is dominated by isolated silanol groups. The major conclusions of the molecular dynamics simulations relevant to the work presented in this paper were that at high OH densities, about 4 times as many molecule-surface hydrogen bonds involved the sp2 oxygen atom as the sp3 oxygens. Only 4% of the interactions involved hydrogen bonding to both sp2 and sp3 oxygen atoms, which is in sharp contrast with the conclusions of the experiments of Kanan and Tripp. At a smaller surface density, more comparable to those experiments, double binding involving both sp2 and sp3 oxygen atoms increased, but only to 6% of the total events. In addition, the lowest-energy binding geometry discovered by Bermudez, in which the sp2 oxygen atom engages two silanol groups simultaneously, was only observed in 8% of the events. In a very recent paper, Taylor et al. reported a combined theory-experiment study of the adsorption of sarin, DMMP, DIMP, and diisopropyl fluorophosphate (DFP) to amorphous silica.15 The adsorption enthalpies were measured in the experiment using inverse gas chromatography. The experimental conditions of these experiments involved relatively high temperatures (above 400 K), and this required the ab initio calculations (performed using MP2 theory and a double-zeta split-valence basis set augmented with polarization and diffusion functions, 6-31++G**) to supplement the standard calculations at 0 K with thermal corrections. Another challenge encountered in the comparison of the experimental and calculated enthalpies is that the surface conditioning enables for the


6

Page 7 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


participation of isolated, geminal, and vicinal silanol groups in the binding. To overcome this challenge, the MP2/6-31++G** calculations considered binding to all three types of silanols using various cluster models ranging from 4 to 6 Si atoms, and the binding energy compared to experiments was obtained from an ensemble average assuming equal populations of the three silanol types. For all OP compounds examined, the calculations showed the strongest binding to involve exclusively the sp2 oxygen atom, which can accept two hydrogen bonds from vicinal silanol groups. Overall, the difference in the adsorption enthalpies between theory and experiment was small for DIMP and sarin (2-3 kcal/mol), but became larger for DMMP (6 kcal/mol) and DFP (7 kcal/mol). Intriguingly, while the calculations with the mimics underestimate the experimental adsorption enthalpy, they overestimated the measured energy for sarin. The somewhat erratic differences between theory and experiment led the authors to call for additional experiments and computations.15 In addition to the potential contribution of various surface sites to the adsorption of OP compounds to silica, further complicating the comparisons between theory and experiment is the fact that relative humidity affects the binding, as shown in experiments by Bermudez.16 In this context, we turn our attention to very recent experiments by Wilmsmeyer et al.,17,18 which avoided many of the experimental issues noted above and thus present an excellent benchmark for accurate theoretical work. In these experiments, a surface made of particulate amorphous silica was conditioned at 700 K in ultrahigh vacuum (UHV). The surface was then exposed to MDCP, DMCP, TMP, DMMP, and DIMP also in UHV, and the binding was monitored via IR spectroscopy. The activation energies for desorption of these OP compounds from silica were subsequently determined


7


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 38

via temperature-programmed desorption. The thermal treatment in the experiment affords a surface that is dominated by isolated silanol groups, and the UHV environment avoids the presence of water, which can affect the surface morphology after treatment and significantly influence the binding of the OP compounds. MP2 calculations using very large basis sets (cc-pVTZ and cc-pVQZ), in which isolated silanol groups were described using a simple gas-phase model (H3SiOH) or a cluster, were also performed to accompany the experiments. The calculations of the binding energies faithfully reproduced the experimental trend in the activation energies of desorption (MDCP < DMCP < TMP < DMMP < DIMP), and showed a direct correlation between the binding energies and the atomic charge of the sp2 oxygen atom in that series of compounds. In the current work, we utilize the experiments by Wilmsmeyer et al. as an accurate benchmark that serves to calibrate our computational approach and further build on the preliminary calculations reported in that contribution.17 Once the electronic structure calculation methods are validated, we utilize them to further quantify the various binding environments discussed in the literature with unprecedented levels of accuracy and to explore the fundamental nature of the forces that control the binding of OP compounds to silica. Calculation methods All of the electronic structure calculations were performed with the Gaussian09 package of programs,19 except for the symmetry-adapted perturbation theory results, which were obtained using PSI4.20 The electronic Schrodinger equation was solved using density functional (B3LYP) and MP2 theory, in combination with split-valence and correlation-consistent basis sets. A subset of calculations were performed using coupled-


8

Page 9 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


cluster theory (CCSD(T)). The binding energies were computed using a dual-level approach, in which the geometry is obtained at an affordable level of theory (MP2 or B3LYP with a 6-31G* basis set), and the energy is refined using those geometries at the MP2/cc-pVTZ and MP2/cc-pVQZ levels. The binding energy was computed as the difference in electronic energy between the OP-silica complex and separated reagents. No zero-point corrections were included, but an earlier study showed this correction to be small.9 The basis set superposition error was removed in all energy calculations with counterpoise method.21 For the systems where both MP2/cc-pVTZ and MP2/cc-pVQZ calculations were possible, an extrapolation to the complete-basis-set limit (MP2/CBS) was performed utilizing the procedure of Halkier et al.22 Method Validation A fundamental aspect of the calculations that we show in this work is the use of a dual-level approach. In a dual-level calculation, the geometry of the system is obtained at an affordable level of theory and subsequently used directly to compute the energy with a higher level of theory. The advantage of this approach is that computationally expensive first and second derivatives required in the geometry optimization and subsequent frequency calculations that verify the location of a minimum are avoided at the higher level. The main assumption of the approach is that the differences in the geometries predicted by the low and high levels are not too large, or equally affect all reagent and product species. To quantify the error of our dual-level approach, we have computed the MP2/cc-pVTZ energy of the dimer formed by MDCP and a gas-phase silanol molecule (SiH3OH) utilizing B3LYP/6-31G*, MP2/6-31G*, and MP2/cc-pVTZ geometries. This dimer corresponds to a geometry in which the sp2 oxygen atom of MDCP acts as the


9


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 38

hydrogen-bond acceptor, H3SiO-H---O=P(OCH3)Cl2. The MP2/cc-pVTZ binding energies calculated using the B3LYP/6-31G*, MP2/6-31G*, and MP2/cc-pVTZ geometries are 5.94, 5.96, and 5.97 kcal/mol, respectively. Clearly, the small differences between the three energies are negligible compared with the total interaction energy, which serves to validate the dual-level approach used in this paper. We note that in the case of intermolecular interactions such as the ones considered in this work, the success of a dual-level approach is rooted in the relative isotropy of the intermolecular potential energy surface around the long-range attractive well. This can be illustrated by the small differences in the MP2/cc-pVTZ energies with B3LYP/6-31G*, MP2/6-31G*, and MP2/cc-pVTZ geometries noted above, even though the SiOH---OR hydrogen bond distances calculated with the low-level methods are slightly different (1.89, 1.92, and 1.87 Å, respectively). The use of dual-level strategies is especially advantageous in the calculations involving cluster models of a silica surface, for which performing multiple geometry optimizations and frequency calculations at the MP2/ccpVQZ, or even MP2/cc-pVTZ levels is currently prohibitive with standard computing hardware. Our choice of MP2 theory to capture the interaction energy between OP compounds and silica is informed by its general success in computing intermolecular interactions.23,24 Nonetheless, given the relatively large interaction energies reported in prior computational work on some OP-silica systems,9,15 we have also benchmarked the accuracy of MP2 calculations with respect to the gold-standard method for intermolecular interactions, coupled-cluster theory with single, double, and perturbative triple excitations, CCSD(T).25 CCSD(T)/cc-pVTZ and MP2/cc-pVTZ energies (with MP2/6-


10

Page 11 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


31G* geometries) for the MDCP-silanol dimer mentioned above differ only by 0.04 kcal/mol (6.00 vs. 5.96 kcal/mol, respectively), which corroborates the superb accuracy of MP2 calculations for the intermolecular interactions investigated in this work. Treatment of conformers A significant challenge in the calculations reported in this paper is the vast conformational space of the investigated OP compounds. Figure 1 shows an example of this space for DMMP by displaying the dependence of the potential energy of DMMP as a function of one of the O=P–O–C dihedral angles. The calculations correspond to a relaxed scan at the MP2/6-31G* level in which all coordinates except for the dihedral angle are allowed to optimize, with energies subsequently refined at the MP2/cc-pVTZ level. The figure shows two stable conformers. In the higher-energy one, the O=P–O–C dihedral angle is ~180˚, corresponding to an ‘anti’ alignment of the P=O and O–CH3 bonds (see inset for a corresponding structure). The lower-energy conformer corresponds to a ‘syn’-type arrangement of those groups. Note that in these calculations, the O=P–O– C dihedral of the other O–CH3 moiety in DMMP was always ‘syn’ to the P=O group, and this conformation also possessed to the lowest energy. These results are consistent with the prior work by Cuisset et al.,26 and similar calculations in all other OP compounds showed that the lowest-energy conformers always correspond to geometries in which the O=P–O–CH3 moiety is not anti, but tends to be roughly syn. Interestingly, the calculations of DMMP adsorption to silica clusters of Taylor et al.15 seem to have used the high-energy conformer of Figure 1, which raises the question of whether the conformer employed in the calculation plays a significant role in the adsorption energy. To illuminate this issue, we have calculated the binding energy of


11


DMMP to a gas-phase silanol molecule using both conformers shown in Figure 1. At the MP2/cc-pVTZ level, the binding energies are 9.53 and 8.97 kcal/mol for the low- and high-energy conformer, respectively. While this difference of slightly over 0.5 kcal/mol between the binding energies of the two conformers is not very large compared to the total binding energies, it becomes significant if one considers that the differences between theoretical and experimental binding energies that we pursue in this work reach values of less than 1 kcal/mol (see below). 6 5

Energy / kcal/mol

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 38

4 3 2 1 0

0

60

120

180

240

O=P-O-C dihedral / deg

300

360

Figure 1. MP2/cc-pVTZ relaxed potential energy scans of one of the O=P–O–C dihedral angles in DMMP. Insets show stable conformers. Color code: P: yellow, C: brown, O: red, H: white. Further complications in the identification of the lowest-energy conformer emerge with the presence of the large, fluxional isopropyl and pinacolyl moieties in DIMP, sarin, and soman. In the live agents, prior work showed that some of the conformers are very close in energy.27,28 In this work, we attempted to calculate binding energies utilizing the


12

Page 13 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


lowest-energy conformers of the OP compounds. In isolated cases, some of the highenergy conformers resulted in larger binding energies, and in those cases, the asymptotic energy of the OP compound utilized in the calculation of the binding energy corresponded to the same conformer adsorbed to the model of silica. Sarin and soman contain a chiral center, but in the absence of a chiral environment the enantiomers are isoenergetic,28 and thus were not separated in the calculations. After benchmarking the computational approach of this work, we now proceed to report the results of our calculations of the interaction energy between MDCP, DMCP, TMP, DMMP, DIMP, sarin, and soman, and models of an amorphous silica surface. Results Gas-phase models As mentioned before, there are various groups on the OP compounds studied in this work that can act as acceptors of a hydrogen bond donated by silanol groups on a silica surface, including the sp2 and sp3 oxygen atoms, and the halogen atoms in MDCP, DMCP, sarin, and soman. An unequivocal ranking for the interaction energy when the OP compounds establish a hydrogen bond with a silanol group via either one of the three acceptors would be provided by standard geometry optimizations, but these failed to locate a minimum in the interaction involving the chlorine atoms in MDCP and DMCP. Instead, the optimizations using a gas-phase silanol molecule started with geometries favorable for a SiO–H---Cl–P hydrogen bond finished with the silanol group interacting with the sp2 or sp3 oxygen atoms in the OP molecules. Therefore, while we can rank which of the sp2 or sp3 oxygen atoms accepting hydrogen bonds in the OP compounds offers a larger binding energy with silanol directly via geometry optimizations, an


13


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 38

estimation of the strength of the SiO–H---Cl–P hydrogen bond cannot be accomplished this way. Nevertheless, a comparison of the relative strength of the intermolecular attraction when Cl atoms act as hydrogen-bond acceptors can be established by computing the intermolecular potential energy surface along a reasonable approach coordinate of the OP compounds to silanol. Figure 2 shows intermolecular potential energy surfaces for the approach of MDCP to a silanol molecule at the MP2/cc-pVTZ level. For the approaches involving hydrogen bonds to the sp2 and sp3 oxygen atoms of MDCP, the potential energy surface corresponds to a frozen scan of the SiOH---OR coordinate in which the rest of variables are held fixed to their values at the respective silanol-MDCP MP2/6-31G* minima. For the approach to the halogen atom, the potential surface corresponds to a frozen scan of the SiOH---ClR coordinate in which the relative orientation of silanol to MDCP is taken from the energy minimum of a silanol/PCl5 geometry optimization. Figure 2 clearly shows that the Cl atom in MDCP is a relatively weak hydrogen-bond acceptor. In addition, hydrogen bonding to the sp2 oxygen atom is more stabilizing than to the sp3 oxygen atom, and this is shown by both the lower energy of the minimum in the potential energy surface and the fact that the SiOH---OR distance at the minimum is longer for the sp3 oxygen atom. For sarin, geometry optimizations aimed at locating a complex involving the halogen atom as an acceptor found a stable minimum. However, the energy of such an interaction is smaller than when either of the two dissimilar oxygen atoms in sarin acts as a hydrogen-bond acceptor. At the MP2/cc-pVTZ level, binding energies involving the O(sp2), O(sp3), and F atom are 8.90, 5.03, and 3.53 kcal/mol, respectively. Therefore, since hydrogen-bonding geometries in which the halogen atom serves as an acceptor do


14

Page 15 of 38

not seem to be competitive, in the following we focus on comparing the hydrogen-bond strengths for the sp3 and sp2 O atoms acting as acceptors.

4

Energy / kcal/mol

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


2

0

-2

-4

-6

1.5

2

2.5

r(H-O,Cl) / Angstrom

3

3.5

Figure 2. MP2/cc-pVTZ frozen intermolecular potential energy scans of the approach of MDCP to silanol. The insets show the geometry of the minima in the corresponding potential curves. Color code: Si and P: yellow, C: brown, O: red, H: white, Cl: green. Figure 3 shows the optimum geometries for the five mimics considered in this work and sarin binding to silanol via the sp2 or sp3 oxygen atoms, and Table 1 lists their relative energies. The calculations show that the strongest interaction between OP compounds and silanol occurs when the sp2 oxygen atom is the hydrogen-bond acceptor in all molecules studied in this work. The binding energies are ~2.5-4.0 kcal/mol smaller when the sp3 oxygen atom is the hydrogen-bond acceptor, and the H---O hydrogen bond distances are on average 0.1 Angstrom longer. Taking the binding energies of DMMP as an example and assuming thermal desorption from the surface, the residence time of a


15


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 38

DMMP molecule bound to silanol through the sp2 O-atom at 300 K is more than 100 times longer than through the sp3 O-atom. This estimation suggests that if access of isolated surface silanol groups to the O atoms of the OP compounds is equally favorable in the experiments, the measured binding energies should be dominated by direct interaction with the O(sp2) of the OP compounds. This finding agrees with the experiments of Wilmsmeyer et al.,17 which reported IR spectral signatures consistent with this binding through the O(sp2) atom (the P=O stretching mode of the OP compounds was noticeably red-shifted upon adsorption to a silica surface compared to its location in the gas-phase spectrum of the pure OP compounds).

&'()%

"!#

%$!"#

%$)'+(%

)',(%

)'+-%

)',3%

)'+&%

)',-%

)'--%

*!

%$"!!

%$)',)% )'-4%

)'+(%

)',-% ".!

%$/01.2%

Figure 3. Optimum geometries of six organophosphorous compounds hydrogen-bonding to silanol via the sp2 (left figures) and sp3 (right) oxygen atoms. Numbers are O-H distances (Å). Colors: Si and P: yellow, C: brown, O: red, H: white, Cl: green, F: blue.


16

Page 17 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 1. Binding energies (kcal/mol) between organophosphorous compounds and silanol calculated at the MP2 level. H3SiO-H---O(sp2) H3SiO-H---O(sp3)

a

cc-pVTZ

CBSa

cc-pVTZ

CBSa

MDCP

5.93

6.90

3.65

4.41

DMCP

7.93

9.06

4.84

5.75

TMP

9.03

10.19

5.63

6.76

DMMP

9.53

10.64

6.00

7.21

DIMP

10.54

11.97

6.79

8.17

Sarin

8.90

10.41

5.03

5.95

MP2 energies extrapolated to the complete-basis-set limit. According to the data in Table 1, the ranking of binding energies for the OP

compounds examined in this work is MDCP < DMCP < TMP < DMMP < DIMP, with the live agent, sarin, being very close to DMCP and TMP. The trend in the five nerveagent mimics reproduces exactly the measured activation energies of desorption by Wilmsmeyer et al.,17,18 which serves to validate the calculations. As discussed in that work,17 the ranking can be rationalized solely based on the polarization of the O(sp2) atoms in each molecule, as the trend in the binding energies follows precisely the trend in the atomic charge in those hydrogen-bond-accepting atoms. Cluster models: Isolated silanol. To delve deeper into the binding energies of OP compounds to silanol groups and explore additional binding geometries proposed in the literature, we have examined a more realistic model of a silica surface. This model considers a grain of silica that contains a single silanol group and is allowed to fully relax its geometry to the lowest


17


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 38

energy conformation. The dangling bonds are saturated with hydrogen atoms. The cluster in this work (H8O9Si6) is slightly larger than the one used for isolated-silanol calculations by Taylor et al.,15 (H1O4Si4, saturated with 9 effective-core-potential F atoms), and similar to the smaller clusters used in analogous calculations by Bermudez.9 All clusters are likely a faithful representation of an isolated silanol group in amorphous silica because the electronic structure of the O-H group bound to a surface silicon atom should be affected only to a small extent by slight changes in the locations of other cluster atoms removed from the group. Figure 4 shows the optimum geometries of the organophosphorous compounds bound to the silica cluster with an isolated silanol group considered in this work, with the O(sp2) atom acting as the hydrogen-bond acceptor. The corresponding MP2/cc-pVTZ binding energies (from geometries obtained at the B3LYP/6-31G* level) are presented in Table 2. Even though the gas-phase calculations shown in Table 1 indicate that binding with an sp3 oxygen atom is substantially weaker than through its sp2 counterpart, we have also considered this possibility in the cluster calculations, and the resulting binding energies are shown also in Table 2. These binding energies through O(sp3) are between 2 and 4 kcal/mol lower than the SiO-H---O(sp2) ones, which ratifies the hypothesis that on a silica surface with isolated silanol groups, the residence time of an O(sp2)-bonded OP compound should be far longer than when an O(sp3) atom is the hydrogen-bond acceptor.


18

Page 19 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


,-/0% ,-.,% !"#

%$,-/1% "!#

%$,-02%

"!!

%$&!

%$,-0.% "'!

%$,-/3% ()*'+%

Figure 4. Optimum geometries of six organophosphorous compounds hydrogen-bonding to an isolated silanol group in a silica cluster via the sp2 oxygen atom. Numbers are O-H distances in Angstrom. Same color code as in Fig. 3. The trend in binding energies for the five OP compounds that mimic nerve agents is the same as in the gas-phase models, and both exactly reproduce the experiments.17,18 In addition, quantitative comparison of the binding energies between the gas-phase and cluster models in Tables 1 and 2 clearly shows that at the same level of theory (MP2/ccpVTZ), the cluster calculations predict binding energies 1-2 kcal/mol larger than the gasphase model. Just as the difference in the binding energies among various OP compounds can be explained via the electronic structure of the O(sp2) atoms, this difference between cluster and gas-phase models can be rationalized on the basis of the electronic structure of the silanol group. In effect, MP2/cc-pVQZ Mulliken atomic charges for the O and H


19


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 38

atoms of the silanol group are -0.624 and +0.268 in the cluster model, and -0.607 and +0.248 in the gas-phase model, respectively. Therefore, while the presence of an extended surface might allow for stabilizing dispersion interactions with regions of the OP compound removed from the hydrogen-bond location (e.g., the isopropyl groups in DIMP), much of the 1-2 kcal/mol difference between the gas-phase and cluster models can be attributed to the difference in the electronic structure of the hydrogen-bonddonating silanol hydroxyl groups. That the hydrogen bond is stronger when the cluster model is involved can be further corroborated via comparison of the O---H distances at the minima. At the same level of theory, these distances are about 0.05-0.10 Å shorter in the cluster model than in the gas phase for each OP compound.

Table 2. Binding energies (kcal/mol) between organophosphorous compounds and a silica cluster calculated at the MP2 level. SiO-H---O(sp2)

SiO-H---O(sp3)

2(SiO-H)---O(sp2)

2(SiO-H)---2O(sp2,sp3)

Exp.d

cc-pVTZ

CBSa

cc-pVTZ

cc-pVTZ

cc-pVTZ

MDCP

7.96

9.07

4.14

12.12 (14.21c)

11.10

10.3±0.2

DMCP

8.97

10.09

7.09

14.60

13.07

11.5±0.2

TMP

10.57

11.63b

8.10

16.50

15.10

12.5±0.1

DMMP

10.95

12.21

8.17

18.04

15.71

13.0±0.1

DIMP

12.39

13.42b

8.66

19.71

17.18

13.8±0.2

Sarin

9.92

10.81b

7.26

15.20

14.23

a

MP2 energies extrapolated to the complete-basis-set limit. MP2/cc-pVQZ energies required for these calculation employed a cc-pVTZ basis set for atoms of the cluster removed from the silanol group. See text. c Extrapolation to the CBS limit using cc-pVTZ and cc-pVQZ energies. d Activation energy for desorption (see the text for details). Ref. 17. b


20

Page 21 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Cluster models: Vicinal silanols. The results presented so far consider that the interaction with OP compounds involves only one isolated silanol group on the silica surface. This approximation is likely sensible for silica surfaces with a high-temperature treatment that have been exposed to OP compounds in the absence of water (UHV conditions). Under such conditions, the density of surface silanol groups can fall below 1.5 nm-2, and the average distance between silanol groups can reach 1.0 nm.12 However, UHV conditions are not accessible in all experiments, and larger silanol surface densities, up to full saturation (4.5 nm-2), are common.7,15 Therefore, it is of interest to investigate binding of OP compounds to silica surfaces in which a high density of silanol groups might permit to establish two simultaneous hydrogen bonds to the same adsorbate. Computational precedent for such bindings is provided by the work of Bermudez, who evaluated two silanol groups binding to the same O(sp2) atom, or to one sp2 and one sp3 O atom, at the B3LYP level with a split-valence double-zeta basis set for TCP, DMMP, and sarin,9 and by Taylor et al.,15 who only considered the first binding possibility, but at the MP2 level (and a slightly different split-valence double-zeta basis set) for DMMP, DIMP, DFP, and sarin. An interesting difference between these two works, aside from the calculation level, is that while Taylor et al. considered binding with two silanol groups for which the Si atoms are nearest neighbors (i.e., the two Si atoms are linked by just one O atom in the cluster), Bermudez’s calculations were performed with silanol groups for which the corresponding Si atoms are next nearest neighbors (i.e., the Si atoms are separated by an O-Si-O moiety). In this work, we considered both possibilities and found that next-nearestneighbor silanol groups produced significantly larger binding energies with our cluster


21


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 38

model. Therefore, our discussion is focused on that arrangement. It is likely that the lower energy obtained with the next-nearest-neighbor geometry is a result of the larger geometric flexibility of that arrangement compared to the tighter nearest-neighbor configuration. Figure 5 shows optimum geometries in complexes between various OP compounds and a silica cluster, where next-nearest-neighbor silanol groups donate hydrogen bonds to the unique O(sp2) atom of the compounds. MP2/cc-pVTZ binding energies are listed in Table 2, and are clearly larger than the ones obtained at the same level of theory with the cluster exhibiting only one isolated silanol (Figure 4). The differences between both calculations range from 4 to 7 kcal/mol and, for all molecules, the binding energies with two silanol groups are slightly smaller than twice the binding energy of those molecules to the silica cluster with just one isolated silanol. This small virtual loss of binding energy from the theoretical maximum is probably due to the geometric requirements for the double binding to occur, which might not overlap perfectly with the next-nearest-neighbor silanols in our cluster. In an actual silica surface, it is likely that there is a distribution of distances between neighboring silanol groups on the surface, which will yield a distribution of binding energies. We have not considered the possibility of that distribution in this work, and present results only for the chosen silica cluster, which we hope to be a reasonable representation of the distribution average. A second reason that the calculated binding energies when two silanols donate hydrogen bonds to a single O(sp2) is smaller than twice the binding energy involving a single silanol group is that in the latter situation the available electronic density around the O(sp2) atom has to be shared between two hydrogen bonds, and this might weaken each


22

Page 23 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


individual bond. Evidence for this hypothesis is provided via analysis of the hydrogenbond lengths in Figures 4 and 5. Using TMP as a reference, we see that when there is a single hydrogen bond, its length is 1.70 Å. On the other hand, that same hydrogen bond becomes 1.84 Å when there is a second hydrogen bond drawing electron density from the same O(sp2) atom.

,-4.%

,-./%

!"#

%$,-.4%

"!!

%$,-.5%

"!#

%$,-.2% ,-/6%

,-.0%

,-.5%

"'!

%$,-.1%

&!

%$,-.3%

,-/0%

,-45%

()*'+%

Figure 5. Optimum geometries of six organophosphorous compounds hydrogen-bonding to two silanol groups in a silica cluster via the sp2 oxygen atom. Numbers are O-H distances in Angstrom. Same color code as in Fig. 3. Regardless of the absolute values of the binding energies, the MDCP < DMCP < TMP < DMMP < DIMP trend is maintained. This is not surprising because, as mentioned before, the strength of the hydrogen bond through the O(sp2) atom correlates with the electronic structure of that atom, and changes to that atomic property due to substituents


23


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 38

from compound to compound should be unaffected by whether one or two surface hydrogen-bond donors are involved. Such a notion can be substantiated by comparison of MP2/cc-pVTZ Mulliken charges for the O(sp2) atoms accepting one hydrogen bond (0.587, -0.656, and -0.668 for DMCP, TMP, and DIMP, respectively) and two hydrogen bonds (-0.637, -0.709, -0.730). This comparison shows how, while donation of two hydrogen bonds exacerbates the polarization of the accepting atom, the trend among OP compounds is not affected by the number of donors. We have also examined the possibility that the OP compounds act as bidentate hydrogen-bond acceptors via utilization of both the O(sp2) and an O(sp3) atom. Geometries for this type of bidentate binding are shown in Figure 6, and energies are in Table 2. Overall, the binding energies of the complexes shown in Figure 6 are only slightly smaller than those in Figure 5, where two hydrogen bonds engage both lone pairs of a single O(sp2) atom. These results imply that even if the distribution of adsorbates bound to a high-silanol-density silica surface is dominated by double binding to the O(sp2) atom at room temperature, a non-negligible fraction of adsorbates will exhibit bidentate binding as shown in Fig. 6. The binding energy when the OP compounds act as bidentate acceptors in Table 2 is slightly larger than the sum of individual binding energies to the sp2 and sp3 oxygen atoms, also presented in that table. As discussed above, the inability of these calculations to yield the theoretical maximum binding energy represented by the sum of the energies of the two individual hydrogen bonds is likely a result of the geometric constraints required for bidentate binding. Even though it was suggested from early experimental work,7 we have not calculated the binding energies of OP compounds acting as bidentate


24

Page 25 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


hydrogen-bond acceptors utilizing two O(sp3) atoms because the foregoing discussion strongly suggests that O(sp3) atoms are not as good hydrogen-bond acceptors as O(sp2) atoms.

4-31%

,-10%

,-./%

,-.,%

,-0/% ,-14%

!"#

%$,-.1%

"!#

%$&!

%$,-.3%

,-0/%

,-./%

"!!

%$"'!

%$,-12% ,-02%

()*'+%

Figure 6. Optimum geometries of six organophosphorous compounds hydrogen-bonding to two silanol groups in a silica cluster via the sp2 and one of the sp3 oxygen atoms. Numbers are O-H distances in Angstrom. Same color code as in Fig. 3. For sarin, the possibility exists for bidentate binding through an O(sp2) and an F atom, as noted by Bermudez.9 However, our MP2/cc-pVTZ calculations show that the binding energy of this geometry (11.68 kcal/mol) is clearly smaller than double binding to the O(sp2) atom (15.20 kcal/mol) or bidentate binding through O(sp2) and O(sp3) oxygen atoms (14.23 kcal/mol).


25


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 38

Comparison with experiment After computing the most representative interactions of OP compounds with models of silica, we now turn our attention to comparing the results of the calculations with experiment as a benchmark of the computational approach. Table 2 shows the experimental results of Wilmsmeyer et al.17,18 Other experimental results are available,7,15 but we prefer to focus first on the ones of Wilmsmeyer et al. because they were performed in UHV, and thus water and other small contaminants do not interfere with the adsorption process and the silica surface is more likely to retain its structure after thermal treatment. Of the four sets of MP2/cc-pVTZ binding energies in Table 2, the ones involving two silanol groups on the cluster overestimate experiment, and those involving an isolated silanol group underestimate it. The use of correlation-consistent basis sets (ccpVXZ) with counterpoise correction for intermolecular calculations is known to produce larger binding energies with a larger cardinal number X=D, T, Q, etc.24 Therefore, the MP2/cc-pVTZ calculations are in all cases an underestimation of the complete-basis-set limit. With this in mind, it can be noted that while calculations with larger basis sets will bring the computed results using an isolated-silanol group closer to experiment, those involving two silanol groups will be further from it. These results are consistent with the fact that the amorphous silica surface employed in the experiments of Wilmsmeyer et al. is dominated by isolated silanols as a result of the thermal history of the surface in those measurements, which should provide a density of ~2 OH nm-2 and an average separation between silanol groups of about 8 Å.12 The adsorption geometries explored in this work for which two silanol groups simultaneously hydrogen-bond to an adsorbate have a distance between silanol groups of about 4.5 Å. Therefore, it is unlikely that double


26

Page 27 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


binding to an O(sp2) atom or bidentate binding are present to a significant extent in the experiments of Wilmsmeyer et al. A remarkable aspect of the MP2/cc-pVTZ results for O(sp2) interaction with the cluster containing one isolated silanol is that not only that they reproduce the experimental trend (MDCP < DMCP < TMP < DMMP < DIMP), but also that they underestimate the experiment by a rather constant amount of about 2 kcal/mol. This somewhat uniform difference between theory and experiment points to a systematic effect, instead of one that depends on the system, and makes a significant difference with earlier comparisons between calculations and measurements where the deviations ranged from -6 to +3 kcal/mol.15 Clearly, the use of basis sets beyond cc-pVTZ in the calculations should bring them closer to experiments, but many of these calculations are at the limit of our current computational power. Nevertheless, cc-pVQZ calculations were possible for MDCP, DMCP, and DMMP, and this enabled us to obtain complete-basis-set estimates for those compounds. In order to be able to obtain this CBS estimate for the larger molecules (TMP, DIMP, and sarin), we drew inspiration from the ONIOM calculations of Bermudez in which the regions deemed to be important for the interaction were treated with a better level of theory than the rest.9 After extensive testing, we determined that the use of a cc-pVQZ basis set for all of the atoms of the adsorbate and the O3Si-OH moiety in the cluster to which it is bound, and a cc-pVTZ basis set for the rest of the ‘inactive’ atoms of the cluster more removed from the adsorbate represented a good compromise for the calculation and entailed little error. For instance, the full counterpoise-corrected MP2/cc-pVQZ energy for the MDCP interaction with the isolated-silanol cluster is 8.60 kcal/mol, and 8.54 kcal/mol when a cc-pVTZ basis set is


27


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 38

used for the ‘inactive’ atoms as defined above. The use of smaller (double-zeta) basis sets for the inactive atoms resulted in larger error. We note that despite adoption of this timesaving strategy, the first of the five calculations involved in computing the counterpoisecorrected binding energy for DIMP entailed 1846 contracted basis functions from 3033 primitive gaussians and 232 electrons, and this carried a significant computational expenditure. The MP2/CBS binding energies obtained from the MP2/cc-pVTZ and MP2/ccpVQZ calculations for the isolated-silanol cluster models are shown in Table 2. As expected, extrapolation to the complete basis set brings the calculated energies closer to the measurements. In fact, the agreement with experiment reaches 1 kcal/mol in some cases, and this represents an unprecedented level of accuracy for this type of systems. Due to the expense of these calculations, MP2/CBS calculations for binding possibilities involving two silanol groups on the cluster were performed for only MDCP accepting two hydrogen bonds via the O(sp2) atom, and as expected, the MP2/CBS energy is further away from experiment than the MP2/cc-pVTZ value (see value between parentheses in Table 2). A note of caution in the comparison between theory and experiment is that while the calculations refer to the difference in electronic energy between the bound species and the reagents’ asymptote, the temperature-programmed desorption experiments provide an activation energy of desorption. Computing the activation energy of desorption from ab initio calculations requires various corrections. First, the zero-point energy needs to be added to the electronic energy to generate the internal energy at 0 K. Then, vibrational, rotational, and translational thermal corrections are needed to bring the


28

Page 29 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0 K energy to the temperature of desorption in the experiment. Subsequently, the internal energy needs to be transformed to enthalpy, and finally, the adsorption enthalpy needs to be transformed to the activation energy of desorption.29 We have performed such sequence of calculations for DMCP, and find that the calculated activation energy of desorption for the interaction with the isolated-silanol silica cluster through the O(sp2) atom is 0.65 kcal/mol larger than the binding energy reported in Table 2. This correction improves agreement with experiment, and sample calculations on other compounds provide similar results. We note that some of the calculations involved in the transformation from binding energies to activation energies of desorption involve significant approximations and can lead to sizeable errors. Especially significant is the relatively low level of theory used in the calculation of vibrational frequencies after geometry optimization (B3LYP/6-31G*), which will affect the accuracy of both the zeropoint energies and vibrational thermal corrections, and the fact that a large number of the vibrations in the extended systems considered in this work are likely highly anharmonic, which will affect the vibrational corrections to a significant extent. A second set of available experiments is that provided by Taylor et al.15 Even though, as mentioned before, the experiments were performed with a different surface of likely higher silanol density, the ranking of adsorption energies for common molecules in that work matches closely our calculations. Indeed, the experiments find a sarin < DMMP < DIMP trend in the adsorption enthalpies, which agrees with our calculations (Table 2) on all binding possibilities explored in this work, including isolated silanol, double binding to the O(sp2) atom, and bidentate binding. The agreement of our calculations


29


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 38

with this second set of experiments further lends confidence to the computational approach followed in this work. Soman Soman is a nerve gas with a chemical structure that resembles quite closely that of sarin. In effect, the only difference between these two agents is that one of the substituents of the central phosphorous atom is a pinacoloxy group (-OCH(CH3)C(CH3)3) in soman, and an isopropoxy group (-OCH(CH3)2) in sarin. The remaining substituents (methyl, O(sp2), and fluoride) are identical. Electronic structure calculations using the gas-phase silanol model reveal that the differences between the binding energies of these two agents are not very large. For instance, at the MP2/cc-pVTZ level, the H3SiO-H--O(sp2) energies for both sarin and soman (structures shown in Figure 7) are within 0.2 kcal/mol of each other. Similar small differences in binding energies are obtained for the H3SiO-H---O(sp3) binding possibility. These small variations in the binding energies from sarin to soman suggest that the electronic structures of the hydrogen-bond acceptors in these agents should be similar, and this can be corroborated by a simple analysis of Mulliken charges. Also at the MP2/cc-pVTZ level, these charges are -0.594(-0.595), -0.491(-0.468), and -0.348(0.350) for the O(sp2), O(sp3), and F atoms in soman(sarin). Given the remarkably close energies obtained for sarin and soman using the gas-phase models, we expect that the differences will not be much larger in the cluster models and have consequently not carried out the pertinent calculations. In the experiment, the more extended nature of the pinacolyl moiety of soman compared to the isopropyl group of sarin might give rise to a slightly stronger attraction with a silica surface via dispersion interactions, or via dipole-


30

Page 31 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


dipole induced interactions if a second surface silanol group is involved. This potential additional stabilization should nonetheless be rather small compared with the strong hydrogen bonds that dominate the adsorbate-surface interaction.

!"#$%&

!'("%&

Figure 7. Optimum geometries of sarin and soman hydrogen-bonding to a gas-phase silanol molecule via the sp2 oxygen atom. Same color code as in Fig. 3. Symmetry-adapted perturbation theory calculations To provide deeper insight into the fundamental aspects that control the strength of the interaction between OP compounds and silica, we have carried out an analysis of the contributions to the binding energy using symmetry-adapted perturbation theory (SAPT).30 SAPT methods are particularly useful in that they provide an accurate estimate of the interaction energy of a complex directly, i.e., without needing to calculate the energy of the asymptote, and that they enable the decomposition of the total interaction energy into electrostatic, exchange, induction, and dispersion contributions. On the other hand, SAPT calculations are far more computationally demanding than regular MP2 calculations. Therefore, here we present only the results of the interaction of the nerveagent mimics with a gas-phase silanol molecule using SAPT2 theory31 with the aug-ccpVDZ basis set.


31


-1

20

Energy / kcal mol

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 32 of 38

10 Electrostatic Induction Dispersion Exchange Total

0

-10

-20

!"#

%$"!#

%$&!

%$"!!

%$"'!

%$Figure 8. SAPT2/aug-cc-pVDZ calculations of the interaction energy of OP compounds with silanol. Figure 8 shows the total SAPT2 interaction energies, segregated into electrostatic, induction, dispersion, and exchange contributions. The total energies follow the trend found in the experiment and the regular MP2 calculations, further validating the SAPT approach. In addition, we note that all separate contributions evolve smoothly along the same trend. Electrostatic interactions account for the majority of the energy, as expected for a hydrogen bond between atoms having relatively large partial charges. The increasing electrostatic energy along the MDCP < DMCP < TMP < DMMP < DIMP trend adds a foundation to the discussion of the experimental trends in terms of atomic charges by Wilmsmeyer et al.17 Beyond electrostatics, the SAPT calculations reveal that


32

Page 33 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


induction and dispersion are comparable to each other and both represent a nonnegligible contribution to the overall binding energy. This result is of importance, as these contributions are often neglected when discussing the factors that make up a hydrogen bond. Electronic repulsion, as characterized by the exchange component, counters the attractive terms, and also grows along the same trend as seen in experiment. Concluding remarks We have employed high-level ab initio calculations to gain insight into the adsorption of seven organophosphorous compounds (MDCP, DMCP, TMP, DMMP, DIMP, sarin, and soman) to amorphous silica. The adsorption is governed by strong hydrogen bonds established between surface silanol groups, which act as hydrogen-bond donors, and hydrogen-bond acceptors in the adsorbates. Of the various hydrogen-bondaccepting atoms in the OP compounds, the oxygen atom linked to the central phosphorous atoms via a double bond establishes the strongest hydrogen bond with silanol groups. The hydrogen-bond strengths when the rest of oxygen atoms in the OP compounds are involved are significantly lower, and binding through halogen atoms is even weaker. Modeling surface isolated silanol groups crudely with a gas-phase SiH3OH molecule captures the essential trend in binding energies among the various OP compounds measured in several experiments, but lacks the accuracy required for quantitative comparison with the measurements. Cluster models are a more faithful representation of an extended surface, and enable study of binding possibilities that involve more than one silanol group on the surface, which likely occur in amorphous silica surfaces of high silanol density. The formation of two hydrogen bonds between the


33


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 38

cluster and the adsorbates unsurprisingly yields higher binding energies than those provided by a single isolated silanol group, and in all OP compounds, double binding to the O(sp2) is stronger than bidentate binding through the O(sp2) and an O(sp3) atom of the adsorbate. Direct comparison between the calculations of this work and recent UHV experiments suggests that the surface used in the measurements is dominated by isolated silanol groups, which is consistent with its thermal history. The differences between the calculated binding energies and the measured desorption energies are rather uniform for the five CWA mimics for which they are available, and reach 1 kcal/mol in many cases. This represents an unprecedented level of accuracy for calculations on OP compounds adsorbing to silica models, and emerges from the use of high levels of ab initio theory (MP2, calibrated with sample CCSD(T) energies) in combination with very large basis sets (cc-pVTZ and cc-pVQZ, with extrapolation to the complete-basis-set limit). Sarin is shown to possess binding energies intermediate between the DMCP and TMP mimics in most cases, and soman gives binding energies similar to those of sarin in isolated-silanol calculations. The high accuracy shown by the calculations of this work will be invaluable in calibrating force fields that can be used in future molecular dynamics simulations. Acknowledgments The support of the Army Research Office, W911NF-09-1-0150, the Defense Threat Reduction Agency, W911NF-06-1-0111, and the National Science Foundation CHE-0547543 is gratefully acknowledged. The authors also acknowledge Advanced Research Computing at Virginia Tech for providing computational resources and technical support that have contributed to the results reported within this paper. URL: http://www.arc.vt.edu


34

Page 35 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


REFERENCES (1) Papirer, E. Adsorption on Silica Surfaces; Marcel Dekker: New York, 2000; Vol. 90. (2) Rimola, A.; Costa, D.; Sodupe, M.; Lambert, J.-‐F.; Ugliengo, P. Silica Surface Features and Their Role in the Adsorption of Biomolecules: Computational Modeling and Experiments. Chem. Rev. 2013, 113, 4216-‐4313 (3) Kim, K.; Tsay, O. G.; Atwood, D. A.; Churchill, D. G. Destruction and Detection of Chemical Warfare Agents. Chem. Rev. 2011, 111, 5345-‐5403 (4) Yang, Y. C.; Baker, J. A.; Ward, J. A. Decontamination of Chemical Warfare Agents. Chem. Rev. 1992, 92, 1729-‐1743 (5) Delfino, R. T.; Ribeiro, T. S.; Figueroa-‐Villar, J. D. Organophosphorous Compounds as Chemical Warfare Agents. J. Braz. Chem. Soc. 2009, 20, 407-‐428 (6) Henderson, M. A.; Jin, T.; White, J. M. A TPD/AES Study of the Interaction of Dimethyl Methylphosphonate with Iron Oxide (.alpha.-Fe2O3) and Silicon Dioxide. J. Phys. Chem. 1986, 90, 4607-‐4611 (7) Kanan, S. M.; Tripp, C. P. An Infrared Study of Adsorbed Organophosphonates on Silica: A Prefiltering Strategy for the Detection of Nerve Agents on Metal Oxide Sensors. Langmuir 2001, 17, 2213-‐2218 (8) Kanan, S. M.; Tripp, C. P. Prefiltering Strategies for Metal Oxide Based Sensors: The Use of Chemical Displacers to Selectively Dislodge Adsorbed Organophosphonates from Silica Surfaces. Langmuir 2002, 18, 722-‐ 728 (9) Bermudez, V. M. Computational Study of the Adsorption of Trichlorophosphate, Dimethyl Methylphosphonate, and Sarin on Amorphous SiO2. J. Phys. Chem. C 2007, 111, 9314-‐9323 (10) Quenneville, J.; Taylor, R. S.; van Duin, A. C. T. Reactive Molecular Dynamics Studies of DMMP Adsorption and Reactivity on Amorphous Silica Surfaces. J. Phys. Chem. C 2010, 114, 18894-‐18902 (11) Kinney, D. R.; Chuang, I. S.; Maciel, G. E. Interior Hydroxyls of the Silica Gel System as Studied by Silicon-29 CP-MAS NMR Spectroscopy. J. Am. Chem. Soc. 1993, 115, 8695-‐8705 (12) Zhuravlev, L. T. The Surface Chemistry of Amorphous Silica. Zhuravlev Model. Colloids Surf. A 2000, 173, 1-‐38 (13) Bolis, V.; Cavenago, A.; Fubini, B. Surface Heterogeneity on Hydrophilic and Hydrophobic Silicas: Water and Alcohols as Probes for HBonding and Dispersion Forces. Langmuir 1997, 13, 895-‐902 (14) Bolis, V.; Fubini, B.; Marchese, L.; Martra, G.; Costa, D. Hydrophilic and Hydrophobic Sites on Dehydrated Crystalline and Amorphous Silicas. J. Chem. Soc., Faraday Trans. 1991, 87, 497-‐505 (15) Taylor, D. E.; Runge, K.; Cory, M. G.; Burns, D. S.; Vasey, J. L.; Hearn, J. D.; Griffith, K.; Henley, M. V. Surface Binding of Organophosphates on Silica: Comparing Experiment and Theory. J. Phys. Chem. C 2013, 117, 2699-‐2708 (16) Bermudez, V. M. Effect of Humidity on the Interaction of Dimethyl Methylphosphonate (DMMP) Vapor with SiO2 and Al2O3 Surfaces, Studied


35


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 36 of 38

Using Infrared Attenuated Total Reflection Spectroscopy. Langmuir 2010, 26, 18144-‐18154 (17) Wilmsmeyer, A. R.; Gordon, W.; Davis, E.; Troya, D.; Mantooth, B.; Lalain, T.; Morris, J. R. Infrared Spectra and Binding Energies of Chemical Warfare Nerve Agent Simulants on the Surface of Amorphous Silica. Submitted. (18) Wilmsmeyer, A. R.; Uzarski, J.; Barrie, P. J.; Morris, J. R. Interactions and Binding Energies of Dimethyl Methylphosphonate and Dimethyl Chlorophosphate with Amorphous Silica. Langmuir 2012, 28, 10962-‐10967 (19) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennuci, B.; Petersson, G. A.; et. al Gaussian 09, Revision A.1 2009. (20) Turney, J. M.; Simmonett, A. C.; Parrish, R. M.; Hohenstein, E. G.; Evangelista, F. A.; Fermann, J. T.; Mintz, B. J.; Burns, L. A.; Wilke, J. J.; Abrams, M. L.; Russ, N. J.; Leininger, M. L.; Janssen, C. L.; Seidl, E. T.; Allen, W. D.; Schaefer, H. F.; King, R. A.; Valeev, E. F.; Sherrill, C. D.; Crawford, T. D. Psi4: an Open-Source Ab Initio Electronic Structure Program. Wiley Interdisciplinary Reviews: Computational Molecular Science 2012, 2, 556-‐565 (21) Boys, S. F.; Bernardi, F. Calculation of Small Molecular Interactions by Differences of Separate Total Energies - Some Procedures with Reduced Errors. Mol. Phys. 1970, 19, 553-‐566 (22) Halkier, A.; Helgaker, T.; Jorgensen, P.; Klopper, W.; Koch, H.; Olsen, J.; WIlson, A. K. Basis-Set Convergence in Correlated Calculations on Ne, N2, and H2O. Chem. Phys. Lett. 1998, 286, 243-‐252 (23) Sinnokrot, M. S.; Valeev, E. F.; Sherrill, C. D. Estimates of the Ab Initio Limit for pi-pi Interactions: the Benzene Dimer. J. Am. Chem. Soc. 2002, 124, 10887-‐10893 (24) Alexander, W. A.; Troya, D. Theoretical Study of the Ar-, Kr-, and XeCH4, -CF4 Intermolecular Potential-Energy Surfaces. J. Phys. Chem. A 2006, 110, 10834-‐10843 (25) Crawford, T. D.; Schaefer, H. F. An Introduction to Coupled Cluster Theory for Computational Chemists. In Reviews in Computational Chemistry, Vol. 14; Lipkowitz, K. B., Boyd, D. B., Eds. 2000; Vol. 14, p 33-‐136 (26) Cuisset, A.; Mouret, G.; Pirali, O.; Roy, P.; Cazier, F.; Nouali, H.; Demaison, J. Gas-Phase Vibrational Spectroscopy and Ab Initio Study of Organophosphorous Compounds: Discriminaton between Species and Conformers. J. Phys. Chem. B 2008, 112, 12516-‐12525 (27) Walker, A. R. H.; Suenram, R. D.; Samuels, A.; Jensen, J.; Ellzy, M. W.; Lochner, J. M.; Zeroka, D. Rotational Spectrum of Sarin. J. Mol. Spectrosc. 2001, 207, 77-‐82 (28) Kaczmarek, A.; Gorb, L.; Sadlej, A. J.; Leszczynski, J. Sarin and Soman: Structure and properties, Struct. Chem. 2004, 15, 517-‐525 (29) Tosoni, S.; Sauer, J. Accurate Quantum Chemical Energies for the Interaction of Hydrocarbons with Oxide Surfaces: CH4/MgO(001)w. Phys. Chem. Chem. Phys. 2010, 12, 14330-‐14340


36

Page 37 of 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(30) Jeziorski, B.; Moszynski, R.; Szalewicz, K. Perturbation Theory Approach to Intermolecular Potential Energy Surfaces of van der Waals Complexes. Chem. Rev. 1994, 94, 1887-‐1930 (31) Hohenstein, E. G.; Sherrill, C. D. Density Fitting of Intramonomer Correlation Effects in Symmetry-Adapted Perturbation Theory. J. Chem. Phys. 2010, 133, 014101-‐12


37


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 38 of 38

Graphic for TOC entry


38

Theoretical Study of the Adsorption of Organophosphorous

Recommend Documents