Theoretical Study of the Microhydration the Chemical Warfare Agent Sulfur Mustard.

A microhydration study of sulfur mustard (SM) was carried out using M06-2X, B3LYP, B3LYP-D3, and MP2 levels of theory with the 6-311++G(2d,2p) basis set. The changes in energetics, structural parameters and vibrational wavenumbers following the addition of up to three discrete water molecules to SM were analyzed. We observed slight changes in the geometry of SM upon microhydration. The stability of hydrated clusters is due to weak C-H···O-H hydrogen bonds. The free energy change for the formation of the clusters is positive at room temperature and becomes exergonic when the temperature decreases. The infrared stretchings of C-Cl of SM and O-H of water are redshifted upon the addition of water molecules. The findings from this work add to the literature of hydrated SM and can be useful in its detection and subsequent destruction.

Introduction

Chemical warfare agents (CWAs) have plagued the world since their initial use during World War I. These dangerous compounds have been classified based on their mode of action on living organisms.[1,2] One such category is known as blister agents or vesicants and regroups chemicals which, upon contact with the skin, invariably causes the formation of blisters.[1,2] The first blister agent used is bis(2-chloroethyl) sulfide and is commonly known as sulfur mustard (SM),[1,2] the structure of which is shown in Figure . SM is known for its delayed but long-lasting harmful effects.[3] It can easily be synthesized and can penetrate various types of fabrics.

Figure 1

Structure of SM.

The toxicity of SM results from the alkylation reaction between the sulfonium ion and essential biomolecules.[1,2,4−6] The most damaging effect of SM is its rapid and irreversible SN2 reaction with the DNA nucleobases,[1,2,4−6] a reaction for which no definite cure has been found till now.[1,2,7,8]

Different methods have been used to neutralize SM in the environment, and one of them is hydrolysis.[9] The hydrolysis of SM in water is preceded by a relatively slow hydration step but results in the formation of the nontoxic thiodiglycol (TDG) and hydrochloric acid under pseudo-first-order conditions.[10,11] Conversely, SM gives rise to relatively toxic polymeric sulfonium products if a limited amount of water is used for the hydrolysis.[11] As well as influencing the toxicity and reactivity[12] of SM, water strongly affects the rate of the SM hydrolysis reaction.[11,13−18] For instance, on activated carbon, the rate of hydrolysis of SM vapor is enhanced by the presence of increasing amounts of water in the micropores of the textile material.[13,14] Similar results have been observed on nanosized metal oxides[15−17] and peroxides[18] used in the decontamination of SM and other CWAs via hydrolysis. In some cases, the absence of water on the metal oxide surfaces resulted in the formation of an alkoxide derivative which remained strongly bound to the surface[16,17] and thus reducing their ability to neutralize the CWA. On the other hand, the volatility and evaporation rate of SM is barely affected by humidity,[19,20] and adsorptivity onto SiO2 nanoparticles.[21]

Interactions between SM and water can occur in the atmosphere,[22] water bodies,[23] and soil[24,25] and on solid surfaces such as paint,[26] concrete,[24,27] or even in the human body.[4−6] SM forms a highly unstable sulfonium ion intermediate via anchimeric assistance, which is a step accelerated in the presence of water.[4−6,12] In addition, during the detection of SM in air samples using mass spectrometry, the presence of water was found to influence detector responses,[22] making the detection of its oxidized derivatives easier than the detection of SM.

Hydrated microclusters formed by the nerve agent sarin have been probed to understand the impact of adding discrete water molecules to its binding energies[28] and the vibrational wavenumber of the phosphate-water[29] hydrogen bond (H-bond). Different nerve agents and their simulants have been microhydrated to determine the extent to which the simulants mimic the properties of the original CWAs in solution.[30] The elucidation of the reaction mechanism through which discrete water molecules assist in reversing the lethal effects of the blister agent lewisite by dithiols have also been studied.[31,32]

Over the past decades, theoretical studies on the DNA alkylation reaction by SM and its nitrogen mustard analogue, mechlorethamine (mec), have also been carried out to understand the mechanism involved in the toxicity of SM.[33−36] A comparison of the reaction in the gas and aqueous phases revealed that the activation barrier for the SM-DNA adduct formed increases when the reaction was carried out in an implicit water solvation model. Ab initio molecular dynamics studies by Mann[35] on mec showed that the specific binding and positioning of explicit water molecules in the first hydration sphere of mec strongly influences the formation of the cyclic aziridinium ion, analogous to the sulfonium ion intermediate. A recent study on the alkylation of adenine and guanine by mec by Larrañaga and co-workers[36] showed that the use of the polarizable continuum model (PCM)[37] overestimates the barrier energy and free energy for the formation of aziridinium. A comparison of the reactivity of SM and mec toward DNA purine bases[33] concluded that the sulfonium ion is more reactive than the aziridinium ion, especially in the aqueous media.[36]

Although the addition of explicit water molecules to other blister agents such as lewisite,[31,32] mec,[35] and the aziridinium ion[36] has been studied, the microhydration of SM is not reported. This has driven us toward the aim of this study, which was to investigate the interaction of up to three discrete water molecules with SM. Within this aim, we were interested in the binding energy, structural parameters, and vibrational wavenumbers of the hydrated SM clusters using theoretical methods.

Computational Details

All computations were performed using Gaussian 09W[38,39] at 298 K and 1 atm running on SEAGrid.[40−43] GaussView 05[44] and Chemcraft 1.7[45] software were used to visualize the results, and all figures were generated using the CYLview[46] program.

Lach and co-workers[47] reported a conformational analysis on SM using the BLYP with Grimme’s D2 dispersion correction in static and molecular dynamics simulations using the Car–Parrinello method with Troullier–Martins pseudopotentials for the core electrons of the molecule. These reported conformers of SM were revisited by using the B3LYP[48,49] functional together with the 6-311++G(2d,2p) basis set for all atoms. This particular method was chosen to compare our results to the microhydration of sarin.[28] The most stable conformer of SM, shown in Figure , was identified as A. The other minimum structures are labeled alphabetically based on their relative stability and are listed in the Supporting Information in Table S1.

Figure 2

Three lowest energy conformers of SM and their ZPE-corrected relative stability in kcal/mol computed at B3LYP/6-311++G(2d,2p) given in parentheses.

Hydrated SM clusters originating from A were optimized using the B3LYP, M06-2X,[50] B3LYP-D3,[51] and MP2[52−57] levels of theory with the 6-311++G(2d,2p) basis set. For the monohydrated complexes, a single water molecule was added to A at different positions and orientations. For the di- and trihydrated clusters, an ensemble of structures was considered whereby the water dimer and trimer were added to A.

All complexes incorporating n water molecules (for n = 1–3) are regrouped under the label G (Generation n), a term introduced by Michaux and co-workers.[58] All complexes are labeled as An-m, where A refers to the A conformer of SM, n represents the number of water molecules, and m is an integer of increasing value representing the relative stability of the complex within the generation considered. The effect of water as bulk solvation on SM was also investigated using the PCM[37] solvation model.

Frequency computations were carried out to check the nature of the stationary points. Selected IR wavenumbers influenced by the formation of H-bonds in the most stable SM·nH2O clusters of each G were identified using VEDA4[59,60] and GaussView 05[44] programs. The computed harmonic wavenumbers are generally overestimated compared to their experimental counterparts, and thus scale factors of 0.947[61] for M06-2X/6-311++G(2d,2p) and 0.963[61] for B3LYP/6-311++G(2d,2p) methods were used for better agreement between the computed harmonic wavenumbers and their experimental counterparts. For the wavenumbers computed using B3LYP-D3 and MP2 with the 6-311++G(2d,2p) basis set, the unscaled values are reported. We report the results obtained from M06-2X/6-311++G(2d,2p) computations as this method performs well for systems involving organic molecules and noncovalent interactions.[62] The results obtained from the other three methods are provided as part of the Supporting Information.

The binding energies (BEs) were calculated using the equations

ΔEBEZPE + BSSE is the binding electronic energy of the complex corrected for zero-point energy (ZPE) and basis set superposition error (BSSE). The BSSE correctional value was obtained using the counterpoise method.[63,64] The ΔGBE is the binding free energy of the hydrated complexes, and E(SM·nH2O) and G(SM·nH2O) are the electronic and free energies of the hydrated SM, respectively. The E(SM) and G(SM) and the E(H2O) and G(H2O) are the electronic and free energies of isolated SM and isolated n water molecules, respectively.

The enthalpy and entropy of the hydrated complexes are also calculated using equations similar to eq . The ΔHBE and ΔSBE are used to represent the change in enthalpy and entropy upon formation of the hydrated complexes, respectively. The H(SM·nH2O), H(SM), and H(H2O) are the enthalpy of the hydrated SM, the isolated SM, and the isolated n water molecules, respectively. The S(SM·nH2O), S(SM), and S(H2O) are used to represent the entropy of the hydrated SM, the isolated SM, and the isolated n water molecules, respectively.

Results and Discussion

Conformational Search

Twenty-two conformers of SM were located using the B3LYP/6-311++G(2d,2p) method (Figure S1). Conformer A is the most stable one, and this is in agreement with previous studies.[47,65,66] However, the relative stability of the other conformers differ from the results of Lach and co-workers.[47]A was used as the free SM molecule to which discrete water molecules were subsequently added. The atom labelings and torsional angles of SM are defined in Scheme S1. The structure of all the conformers of SM along with their relative stability obtained using B3LYP/6-311++G(2d,2p) are shown in Figure S1, and their torsional angles are listed in Table S1.

Effect of Discrete Solvent Molecules and Bulk Solvation on the Properties of SM

Microhydration: Effect on Energetics

The addition of discrete water molecules to A gave rise to three monohydrated, seven dihydrated, and nine trihydrated complexes when the M06-2X/6-311++G(2d,2p) method was used. The SM·nH2O (for n = 1–3) complexes obtained are shown in Figure , and their binding energy, enthalpy, entropy, and free energy values are summarized in Table . The hydrated complexes obtained using the other three levels of theory are displayed in Figures S2–S4, and their energy and thermodynamic parameters are listed in Tables S3–S5.

Figure 3

Complexes of SM·nH2O (for n = 1–3) optimized using the M06-2X/6-311++G(2d,2p) method, and these are arranged in order of decreasing relative stability within each generation. Atom labelings are also included.

Table 1

Energetics of All the Complexes Obtained Using the M06-2X/6-311++G(2d,2p) Methoda

generation	cluster	ΔEBEZPE + BSSE (kcal/mol)	ΔEBEZPE (kcal/mol)	ΔHBE (kcal/mol)	ΔSBE (cal/mol/K)	ΔGBE (kcal/mol)
G0	A	0	0	0	0	0
G1	A1-1	–6.75	–7.44	–7.67	–30.30	1.37
	A1-2	–3.19	–3.70	–3.92	–27.80	4.37
	A1-3	–3.01	–3.48	–3.67	–27.05	4.40
G2	A2-1	–12.14	–13.44	–14.38	–60.54	3.67
	A2-2	–11.98	–13.27	–14.28	–61.27	3.98
	A2-3	–9.52	–10.76	–11.44	–56.74	5.47
	A2-4	–9.43	–10.61	–11.41	–58.70	6.09
	A2-5	–8.95	–10.05	–10.84	–58.52	6.60
	A2-6	–8.86	–9.97	–10.65	–56.25	6.12
	A2-7	–8.64	–9.74	–10.59	–58.11	6.73
G3	A3-1	–19.63	–21.88	–23.64	–93.37	4.20
	A3-2	–19.41	–21.70	–23.68	–97.76	5.47
	A3-3	–17.25	–19.36	–20.82	–94.08	7.23
	A3-4	–17.22	–19.19	–20.67	–89.98	6.16
	A3-5	–17.19	–19.40	–21.04	–93.29	6.78
	A3-6	–16.90	–19.06	–20.56	–89.24	6.05
	A3-7	–16.15	–18.18	–19.98	–93.52	7.91
	A3-8	–14.51	–16.11	–17.65	–87.07	8.31
	A3-9	–14.47	–16.36	–17.52	–87.63	8.61

The , , , , and were all calculated at 298.15 K and 1 atm.

The monohydrated structures obtained at the MP2 level of theory are comparable to the M06-2X results, while those obtained at B3LYP and B3LYP-D3 show slight variations in the number of H-bonds formed (Figures S3–S5), which affects the relative stability and binding electronic, enthalpy, entropy, and free energies of the hydrated complexes. It is worth pointing out that one of the starting structures of A with the dimer of water leads to a C dihydrated complex (Figure S6). All the complexes obtained have positive ΔGBE values, which is due to the decrease in entropy accompanying the combination of n water molecules with A to form a single entity.

The A1-1 structure is the most stable of the three monohydrated complexes at M06-2X, B3LYP, B3LYP-D3, and MP2 levels, with a relatively large ΔEBEZPE + BSSEvalues of −6.75, −2.34, −5.83, and −4.45 kcal/mol, respectively. The greater stability of the A1-1 complex compared to the other two monohydrated structures at all levels of theory is due to the larger number of H-bonding formed and, in particular, the involvement of four C–H···O–H interactions at M06-2X, B3LYP-D3, and MP2 levels and three C–H···O–H for the B3LYP level of theory.

Increasing the number of explicit water molecules in the complexes increases the magnitude of the ΔEBEZPE + BSSE, ΔEBEZPE, and ΔGBE values. The increase in binding energy may be explained on the basis of the effect of cooperativity between the water dimer and A, which forms larger cyclic networks of H-bonds than the water monomer does with A. As a result of the larger ring sizes resulting from H-bonding, a reinforcement of the H-bonds occurs, leading to the lower ΔEBEZPE + BSSEvalues of the dihydrated complexes compared to the monohydrated ones. Similar observations are made when the hydrated clusters are compared at B3LYP, B3LYP-D3, and MP2 levels of theory.

For the trihydrated complexes, all four ab initio methods give similar structures with very few clusters showing slight variations in the distances and orientation of the water molecules. The three solvent molecules occupy the first hydration shell in all clusters except for A3-4 and A3-8, where two and one water molecules, respectively, lie in the first shell, while the remaining water molecule(s) are in the second hydration shell. These complexes are the third and fourth most unstable structures of the G group with relatively large and positive ΔGBE values. The relative stability among the SM·3H2O aggregates differ for each level of theory. The M06-2X, B3LYP-D3, and MP2 levels of theory predict the A3-1 complex as the most stable structure when three water molecules are added to A, while B3LYP predicts the A3-4 structure of M06-2X as the most stable trihydrated complex.

The trends in ΔEBEZPE with and without BSSE and ΔGBE for the most stable SM·nH2O for each G are shown in Figure a,b. Increasing the number of water molecules present around SM affects the ΔEBEZPE + BSSE, ΔEBEZPE, and the ΔGBE in a nonlinear way. Along with the increase in the number of molecules, a rise in the BSSE correctional values can also be noted. Applying the correction for BSSE decreases the magnitudes of the ΔEBEZPEof the clusters as shown in Table and Figure a.

Figure 4

(a) Variation of the binding energy with the number of water molecules present around SM for the most stable hydrated complex of each generation computed using the M06-2X/6-311++G(2d,2p) method. (b) Variation of the free binding energy with the number of water molecules present around SM for the most stable hydrated complex of each generation computed using the M06-2X/6–311++G(2d,2p) method.

Microhydration: Effect on Structure

In all the cases considered, the addition of discrete water molecules has a negligible effect on the isolated structure of A. Similar observations have been made for the aziridinium ion of mec when it was microhydrated with one to four water molecules.[36] For the most stable mono-, di-, and trihydrated complexes of SM, the C–Cl bonds involved in the H-bonding are extended by 0.01, 0.02, and 0.01 Å, respectively, at all four levels of theory. This bond is also lengthened in all other complexes where the C–Cl fragment is bound to water molecules by similar values. Likewise, the bond and torsional angles of A are barely affected when water molecules are added, changing by less than 9° for all the methods used. In the A3–7 trihydrated complex, the two gauche torsional angles are both changed from 82° to 51° and 61° when M06-2X was used. A similar large decrease in torsional angle is observed at the MP2 level, but computations using B3LYP and B3LYP-D3 for the same trihydrated complex results in less drastic changes in the torsional angles of A. The structural parameters of the most stable complexes are reported in Table for M06-2X and Tables S6–S8 for the other three methods used.

Table 2

Structural Parameters of the Most Stable SM·nH2O Clusters Using the M06-2X/6-311++G(2d,2p) Method

	isolated A	A1-1	A2-1	A3-1
Bond length (Å)
C2–Cl	1.80	1.81	1.82	1.81
C2′–Cl′	1.80	1.81	1.81	1.81
C1–C2	1.51	1.51	1.51	1.51
C1′–C2′	1.51	1.51	1.52	1.51
S–C1	1.82	1.82	1.82	1.82
S–C1′	1.82	1.82	1.82	1.82
C1–H1	1.09	1.09	1.09	1.09
C1′–H1′	1.09	1.09	1.09	1.09
C1–H2	1.09	1.09	1.09	1.09
C1′–H2′	1.09	1.09	1.09	1.09
C2–H3	1.09	1.09	1.08	1.09
C2′–H3′	1.09	1.09	1.08	1.08
C2–H4	1.08	1.08	1.08	1.09
C2′–H4′	1.08	1.08	1.08	1.08
Bond angle (°)
C1–C2–Cl	110.2	110.0	110.1	110.2
C1′–C2′–Cl′	110.2	110.0	109.8	110.4
S–C1–C2	111.1	111.2	110.7	110.5
S–C1′–C2′	111.1	111.2	110.9	110.9
C1′–S–C1	100.2	100.0	100.9	100.8
Torsional angle (°)
Cl–C2–C1–S (ϕ1)	–179.5	–174.2	–176.6	–174.9
C2′–C1–S–C1 (ϕ2)	82.4	76.3	77.2	79.2
C2′–C1′–S–C1 (ϕ3)	82.4	76.3	77.7	78.0
Cl′–C2′–C1′–S (ϕ4)	–179.5	–174.2	174.2	170.5

In all the complexes formed, the H-bonding motifs are cyclic in nature and consist of five- to eight-membered rings systems between SM and water. Cyclic H-bonding networks are usually favored over other possible arrangements due to the occurrence of cooperative effects between the single donors and acceptors of the ring, which strengthen the interaction and further stabilize the complex relative to the isolated starting molecules. Moreover, anticooperativity, which occur due to the formation of H-bonds between a single donor (or acceptor) and multiple acceptors (or donors), is also known to contribute to the stability of H-bonded complexes.[67] The H-bonding distances between the H atom of the donor (X–H) and acceptor atoms (Y) and the H-bonding angle (X–H···Y) are reported in Table for the M06-2X level of theory. The same parameters for B3LYP, B3LYP-D3, and MP2 are reported in Tables S9–S11 and are comparable to the values obtained at the M06-2X level for the same complexes. Based on the distances and angles of the noncovalent interactions observed, most of the H-bonds are weak. In particular, the C–H···O–H interaction has H-bond lengths of greater than 2.50 Å and bond angles of less than 130° in the majority of the clusters, emphasizing its weak nature. Nevertheless, the involvement of multiple C–H···O–H bonds at the same time directly accounts for the greater relative stability of certain complexes within the same generation, as demonstrated clearly by the A1-1, A2-1, A2-2, and A3-1 complexes.

Table 3

H-Bond Lengths and Angles for the Most Stable SM·nH2O Complexes in Each Generation Using the M06-2X/6-311++G(2d,2p) Method

H-bonds	H-bond length (Å)	H-bond angle (°)
A1-1
C1–H2···O1–H5	2.49	117.7
C2–H4···O1–H5	2.53	108.3
C1′–H2′···O1–H5	2.49	120.4
C2′–H4′···O1–H5	2.53	119.2
O1–H5···Cl	2.82	110.9
O1–H6···Cl′	2.82	110.9
A2-1
O1–H6···O2–H7	1.88	161.2
C1–H2···O1–H5	2.58	117.3
C2–H4···O1–H5	2.46	120.6
C1′–H2′···O1–H5	2.40	120.4
C2′–H4′···O1–H5	2.55	106.7
O1–H5···Cl–C2	2.62	127.4
O2–H7···Cl′–C2′	2.37	161.1
A3-1
O1–H6···O2–H7	1.82	154.2
O2–H7···O3–H9	2.02	140.6
O3–H9···O1–H6	1.89	152.3
C2–H4···O2–H7	2.55	175.8
C1′–H2′···O2–H7	2.33	159.2
C1–H2···O1–H5	2.91	114.9
C2–H4···O1–H5	2.57	128.6
C1′–H2′···O1–H5	2.50	128.4
C2′–H4′···O1–H5	2.93	100.4
O1–H5···Cl–C2	2.48	145.6
O3–H10···Cl′–C2′	2.91	108.6

Microhydration: Effect on Vibrational Modes and Wavenumbers

The formation of H-bonds between water and SM can be confirmed via IR spectroscopy. For the nonhydrated A structure and the isolated n water clusters, experimental and theoretical vibrational wavenumbers are reported in Table .

Table 4

Selected Vibrational Wavenumbers of the Most Stable Complexes Computed Using M06-2X/6-311++G(2d,2p) Method and Scaled by 0.947[61]a

	IR wavenumbers (cm–1)
bonds and modes	A	(H2O)n	A1-1	A2-1	A3-1
Water
ν(as)O–H		3783	3753 (−30)
ν(s)O–H		3687	3651 (−36)
δ(s)H–O–H		1551	1518 (−33)
Water dimer
ν(as)O–H		3773		3749 (−24)
ν(as)O′–H′		3766		3699 (−67)
ν(s)O–H		3678		3609 (−69)
ν(s)O′–H′		3606		3520 (−86)
δ(s)H–O–H		1571		1575 (4)
δ(s)H′–O′–H′		1554		1560 (6)
Water trimer
ν(as)O–H		3757			3736 (−21)
ν(as)O′–H′		3754			3725 (−29)
ν(as)O″–H″		3752			3686 (−66)
ν(s)O–H		3526			3586 (61)
ν(s)O′–H′		3520			3488 (−32)
ν(s)O″–H″		3463			3360 (−103)
δ(s)H–O–H		1587			1594 (7)
δ(s)H′–O′–H′		1566			1575 (9)
δ(s)H″–O″–H″		1563			1558 (−5)
SM·nH2O
ν(s)Cl–C	704 (650, 702)b		685 (−19), 676 (−27)	685 (−19), 673 (−30)	684 (−20), 674 (−29)
ν(s)C–H (internal C)	2931 (2915), 2932		2949 (18), 2951 (20)	2937 (6), 2943 (11)	2920 (−11), 2945 (13)
ν(s)C–H (external C)	2957 (2933), 2957		2970 (12), 2971 (13)	2966 (9), 2971 (14)	2967 (9), 2972 (14)

Values given in parentheses are the changes in vibrational wavenumbers relative to the isolated A and n water cluster, respectively.

Experimental vibrational wavenumbers of isolated SM.[68]

In general, harmonic wavenumbers calculated theoretically are overestimated compared to experimental IR wavenumbers. The values obtained for gaseous A using M06-2X and MP2 are found to differ by at least 100 cm–1 from those obtained experimentally. The use of a scale factor becomes necessary for better agreement with experimental or highly accurate theoretically predicted spectroscopic data. Table shows the scaled wavenumbers obtained using M06-2X/6-311++G(2d,2p) for the water clusters, the nonhydrated A, and the most stable hydrated complexes of SM in each generation. The scaled IR vibrational wavenumbers obtained using B3LYP and those unscaled obtained using B3LYP-D3 and MP2 for the compounds are provided in Tables S11–S13.

The IR spectra of A and the most stable hydrated SM·nH2O complexes along with the spectra of isolated n water clusters are shown in Figure .

Figure 5

IR spectra of the isolated A and water clusters and the most stable SM·nH2O complexes.

Microhydrating SM with one to three water molecules results in the appearance of new peaks in the spectra of the CWA in the 0–500 cm–1 region, which corresponds to the intermolecular H-bonds formed in all the hydrated complexes. The C–Cl stretching modes of SM are split and redshifted by 19 and 27 cm–1 in A1-1, forming two very close peaks that remain at almost the same wavenumbers even as the number of water molecules in the complex increases. On the other hand, the C–H bonds next to S (internal C–H) and Cl (external C–H) atoms are blueshifted by relatively small values. The increase in wavenumber for stretching C–H is a phenomenon often encountered when sp3-hybridized alkyl moieties act as a H-bond donor in complexes due to small reductions in the C–H bond length.[69]

The most notable changes in IR wavenumbers occur for the stretching O–H modes of water in the hydrated clusters of SM when compared to the isolated water clusters in the gas phase. For the n = 1 solvent molecule, the symmetric, asymmetric stretching and bending modes of the bonded O–H of water in A1-1 are redshifted by 36, 30, and 33 cm–1, respectively, compared to a single isolated water molecule.

In the A2-1 complex, the stretching and bending modes of O–H of water are split and shifted to a different extent depending on the binding positions of the water molecule and whether the O–H is free or bonded to an acceptor or donor in SM. The peak corresponding to the O–H H-bonded to the second water molecule is redshifted by 86 cm–1 upon binding to A. The single free O–H bonding in A2-1 vibrates at a higher wavenumber than the remaining three bonded O–H stretch peaks with a shift of −24 cm–1 compared to the isolated water dimer. The more pronounced change in wavenumber is indicative of the stronger H-bonding between the two polar water molecules than between the hydrophobic SM and water.

For the trihydrated A3-1 complex, the stretching O–H modes are redshifted to a greater extent compared to the A2-1 complex. The free and bonded O–H stretch modes of three water molecules in the isolated water trimer formed two separate peaks, which are split in the spectrum of A3-1 due to their participation or nonparticipation in H-bonding. The two O–H bonds that remain free in A3-1 are redshifted by 21 and 29 cm–1 and form two extremely close peaks, while the O–H bound to Cl is redshifted by 66 cm–1 and forms a separate peak of greater intensity.

The bonded O–H modes, which were previously bound to other water molecules in the isolated trimer, give rise to three distinct peaks in the spectrum of A3-1. Of these, two stretching O–H modes are redshifted, while one is blueshifted. The extent of the redshift is influenced by the number of H-donors and acceptors of A bound to the water molecule. The water molecule forming weak acceptor H-bonds with two C–H and strong donor and acceptor bonds with the other two water molecules forms the most intense peak in the trihydrated complex spectrum and shows the largest redshift of 103 cm–1.

In contrast, the water molecule forming a single strong donor H-bond with another water molecule and weaker bonds with five different H-donors forms the second most intense peak, which is shifted by −33 cm–1. The blueshifted O–H peak arises from the water molecule acting as a double donor, forming weak H-bonds with Cl and another water molecule and a strong bond with the third water molecule in the complex. The increase in the IR wavenumber of O–H is also observed in A3-2 and A3-7 both of which also have lower dipole moment than A. This unusual increase in wavenumber is observed at all the levels of theory considered.

Bulk Solvation of SM Using PCM

The structural parameters of A in water change negligibly compared to the gas phase. The largest change in the geometry is an increase of 0.01 Å in the two C–Cl bonds at all levels of theory. The electronic solvation and Gibbs free solvation energies of A with the M06-2X functional are both −4.66 kcal/mol. The corrected solvation free energy[70] using the same functional is −8.44 kcal/mol. The solvation free energy is comparable to those obtained using the three other methods (Table S14). These values indicate the stabilization of the A molecule in the presence of a bulk solvent. When discrete water molecules were added, the free binding energies are all positive values due to the lower entropy of the microsolvated complex with respect to the separate SM and water molecules.

Comparison with Other Chemical Warfare Agents

Compared to the microhydrated sarin conformers at B3LYP/6-311++G(2d,2p) studied by Alam and co-workers,[28] the addition of discrete water molecules to SM give rise to more varied H-bonding networks and binding positions for n = 1–3 water molecules. The structures of both CWAs are negligibly affected by the presence of discrete water molecules in their hydration spheres. The ΔEBEZPE + BSSE values of the most stable SM·nH2O clusters have slightly lower magnitudes than those of sarin·nH2O, while the ΔGBE values of the blister agent are higher. The reason for the difference in energies is directly related to the chemical structure of the CWAs and the binding positions of water in the hydrated complexes. For sarin, water binds to the polar P=O bond, while in SM, the water molecules bind to less polar bonds and atoms. The energy difference is further illustrated when the solvation of the two CWAs is considered: sarin dissolves much more readily in aqueous solutions than SM at 25 °C.[71]

Nanoparticles can hydrolyze SM into TDG, and the presence of water is known to accelerate the hydrolysis reaction. The general model proposed[17,18,20] is that SM is first adsorbed onto the surface of the nanoparticle and reacts with free weakly bound water molecules present on the surface of the nanoparticle. The water molecules assist in the breaking of C–Cl bond and results in the formation of the sulfonium ion. In the molecular dynamics simulations of Mann,[35] the cleavage of C–Cl for the mec molecule was observed to occur after the closest solvent molecules formed H-bonds with the Cl atom after which it lengthened and broke. The reactive sulfonium intermediate formed on the nanoparticle then reacted with hydroxyl or water nucleophiles and was liberated as TDG.

The binding energies of SM with DNA bases were calculated using the B3LYP/6-31+G(d) method and found to be −36 and −22 kcal/mol in the gas and aqueous phases (PCM) at the N7 of the guanine DNA base, respectively.[33] The large difference in energy between SM-DNA adducts and SM·nH2O complies with in vivo experimental observations: SM tends to attack purine DNA bases in the aqueous cellular environment rather than undergo hydrolysis despite its instability in aqueous media after hydration.[12] Larrañaga and co-workers[36] modeled the alkylation of DNA by aziridinium microhydrated with two water molecules and concluded that the presence of water in the hydration shells lowers the activation barrier for the aziridinium formation, which precedes the DNA alkylation reaction step. Further, the water present in the first hydration sphere of aziridinium then form H-bonds with the approaching DNA molecule, aiding in the SN2 reaction.

Conclusions

The microhydration of the harmful blister agent sulfur mustard with n = 1 to 3 explicit water molecules was studied theoretically. From the symmetrical global minimum A, three monohydrated, seven dihydrated, and nine trihydrated complexes were obtained at all levels of theory used. The n water molecules added occupied the first hydration shell of the CWA in all the clusters optimized except for two trihydrated structures where the second and third water molecules entered the second solvation shell.

The water molecules form weak hydrogen bonds with the different donors and acceptors of hydrogen in A, resulting in negligible effects on its geometry. The only notable exception observed is the lengthening of the C–Cl bond by 0.01–0.02 Å in the different complexes where water binds to the Cl atom. This particular bond is responsible for the reactivity of SM, and its cleavage is known to be enhanced in aqueous media. A case of conformational change was also observed in which A adopted the structure of the third most stable conformer of SM. Nevertheless, this complex has a relatively higher positive free binding energy compared to the remaining dihydrated complexes. The binding energies of the clusters become increasingly negative as the number of water molecules added is raised. This highlights the role of water molecules in stabilizing the hydrated complexes compared to the isolated SM molecule. On the other hand, the free binding energy becomes more positive as the change in entropy becomes more negative. The occurrence of the noncovalent interactions is reflected in the appearance of new peaks in the 0–500 cm–1 region of the IR spectra of the hydrated SM and the noticeable shifts in the O–H stretching wavenumbers of water. Increasing the number of water molecules enhances theses effects and can aid in the identification of SM by IR spectroscopy in situations where the CWA is suspected to have been released in the environment.

Overall, the fundamental research on the microhydration of SM complexes can be particularly useful in creating hydrolysis reaction models, especially on nanoparticles and metal–organic frameworks, to support experimental research aimed toward the development of detection and elimination techniques for SM. It can also be used in computational studies on the toxicology of SM to provide a more complete description of the situation and aid in combatting the CWA.