Article pubs.acs.org/JPCA

Molecular Dynamics Simulations of Small Clusters and Liquid Hydrogen Sulfide at Different Thermodynamic Conditions M. Albertí,*,† A. Amat,‡ A. Aguilar,† and F. Pirani§ †

IQTCUB, Departament de Química Física, Universitat de Barcelona, 08028 Barcelona, Spain Computational Laboratory for Hybrid/Organic Photovoltaics (CLHYO), CNR-ISTM, Via Elce di Sotto 8, 06123 Perugia, Italy § Dipartimento di Chimica, Biologia e Biotecnologie, Università di Perugia, 06123 Perugia, Italy ‡

S Supporting Information *

ABSTRACT: A new force field for the intermolecular H2S−H2S interaction has been used to study the most relevant properties of the hydrogen sulfide system from gaseous to liquid phases by means of molecular dynamics (MD) simulations. In order to check the validity of the interaction formulation, ab initio CCSD(T)/aug-cc-pVTZ calculations, including the counterpoise correction on the H2S, (H2S)2, and (H2S)3 structures optimized at the MP2/aug-cc-pVDZ level, have been performed. The (H2S)2,3 systems have been characterized by performing NVE MD simulations at decreasing values of the temperature, while the liquid sulfide behavior has been investigated considering a NpT ensemble of 512 molecules at several thermodynamic states, defined by different pressure and temperature values. Additional calculations using an ensemble of 2197 molecules at two different temperatures have been performed to investigate the liquid/vapor interface of the system. The S−S, S−H, and H−H radial distribution functions and the coordination number, calculated at the same conditions used in X-ray and neutron diffraction experiments, and the evaluated thermodynamic and structural properties have been compared successfully with experimental data, thus confirming the reliability of the force field formulation and of the MD predictions.

1. INTRODUCTION Hydrogen sulfide (H2S) is a toxic component of natural gas and crude oil1 that, at ambient conditions, is colorless, very poisonous, corrosive, flammable, and explosive. Natural gas can contain high quantities of hydrogen sulfide, while only small amounts occur in crude petroleum. The most important source of the H2S emissions due to the human activity comes from the petroleum refineries, and even though the levels of H2S in crude oil are low, they present serious health and safety risks to workers.2 Despite the risks associated with the H2S emissions, its use as gaseous signaling molecule for several applications has been pointed out in the past.3−5 Moreover, H2S prodrugs and novel agents might be efficacious for acute myocardial infarction, stroke, diabetes, arthritis, peripheral artery disease, metabolic syndrome, organ transplantation, erectile dysfunction, diabetes, inflammatory bowel disease, and pulmonary hypertension,3,6,7 though in such applications the highly toxicity of H2S has to be taken into account. From an environmental point of view, the phase behavior in water−H2S systems is crucial for the design and operation mode of pipelines and production/processing facilities.8 Although understanding the H2S interaction with other molecules is crucial for both environment and development of H2S-based therapeutic agents, due to its high toxicity, experiments are very scarce, in particular under high pressure and temperature conditions as those of petroleum reservoirs. For all these reasons, theoretical studies, addressed to the characterization of © XXXX American Chemical Society

the hydrogen sulfide properties under various conditions, become particularly important and interesting. Such studies represent also a priority to properly extend the investigation to more complex systems involving H2S and other molecules. From a theoretical point of view, molecular dynamics (MD) can be a useful tool to investigate hydrogen sulfide at different p−T thermodynamic conditions. Clearly, the reliability of the simulations depends on that of the force field used, and while a variety of force fields have been constructed for water, including the well-known SPC,9 SPC/E,10 TIP3P,11 TIP4P,12 TIP5P,13 and AMOEBA14 models, less information exists for hydrogen sulfide. One of the first attempts to describe hydrogen sulfide interactions was undertaken Jorgensen 3 decades ago,15 followed by the four-site model of Forester et al.,16 lately reparametrized by Kristóf and Liszi17 in order to reproduce the vapor/liquid coexistence curve. Potoff and co-workers,18 considered four sets of different partial charges, from which the relevant parameters of the Lennard-Jones (LJ) were selected. More recently, Drude polarizable force fields have been also proposed. 19,20 In order to overcome some Special Issue: Piergiorgio Casavecchia and Antonio Lagana Festschrift Received: December 3, 2015 Revised: February 2, 2016

A

DOI: 10.1021/acs.jpca.5b11843 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

that include those close to the critical point, are performed. Moreover, in order to analyze the description of the liquid/ vapor interface, a system containing a total of 2197 molecules of H2S has been simulated at pressure and temperature close to those of the boiling point. The paper is organized as follows: in section 2 the methodology used is described. In section 3, ab initio calculations on the hydrogen sulfide dimer and trimer are presented and compared with the results of MD simulations performed on the small clusters. In section 4, simulations on liquid sulfide at selected conditions of pressure and temperature are discussed. Concluding remarks are given in section 5.

deficiencies of the previous potentials, one of the last force fields adapted a transferable potential model for phase equilibria (TraPPE).21 Most of the previous potentials proposed for H2S have been parametrized to describe accurately only some macroscopic thermodynamic and structural properties of the system. In order to extensively exploit the predictive capabilities of a given force field, going from simple to more complex aggregates, difficulties related to the proper description of weak interactions resulting from the balancing of both dispersion and induction attraction, dominant at long-range, with exchange (size) repulsion, dominant at short-range,22 must be overcome. In the past decade, some of us developed a potential energy model for weak interactions where the combination of effective (i.e., including also less important terms) electrostatic and nonelectrostatic components allow us to describe consistently and accurately both neutral and ionic systems of increasing complexity. In particular, the potential model has been used successfully to describe hydrogen bonds (H-bonds) and it has been applied to investigate small clusters and liquid water,23−25 the bonds between water, hydrogen sulfide, ammonia and methane with benzene,26 and some small (NH3)2−5 clusters and liquid ammonia.27 Such potential model exploits an accurate description of the energy contribution arising from the balance of dispersion−induction attraction with exchange repulsion, globally represented by Vnel and formulated as an improved Lennard-Jones (ILJ) function.28 In such function a shape parameter (see next section), not present in the traditional LJ model, indirectly accounts for the role of less important effects and for the incomplete separability, at intermediate and short range, of the basic interaction components. The relevant parameters of the ILJ function, unlike what happens with other potential models, are not parametrized but predicted by using the values of molecular permanent charges, multipoles, and electronic polarizabilities.29 Moreover, the ILJ function can be combined with different representations of the electrostatic contributions (Vel) arising from the interaction of permanent charges and/or multipoles of the involved partners. In our preceding studies, though Vel has been calculated without including explicitly polarization effects, the model predictions have been successfully compared with both high quality experiments28 and accurate ab initio calculations.30,31 In the model, electronic polarization effects can be indirectly accounted for by modifying the value of the permanent multipole moments and, consequently, the charge distribution defining Vel and/or that of the additional shape parameter in the ILJ function (see next section). For instance, increases of about 13% and 36% of the dipole moment with respect to those corresponding to gas phase were found adequate to describe small clusters of water23−25,32 and ammonia,27 respectively. All these results encouraged us to use the same methodology to investigate hydrogen sulfide. In the present study, the hydrogen sulfide potential energy function is formulated following the same methodology used for water and ammonia. The proposed formulation is used in MD simulations to study small clusters of H2S and liquid hydrogen sulfide. In particular, the predictions of the model for the dimer and trimer, (H2S)2 and (H2S)3, are compared with results from ab initio calculations and with predictions from other potential models. Having assessed the reliability of the model, MD simulations considering an ensemble of 512 H2S molecules, at different conditions of pressure and temperature

2. COMPUTATIONAL METHODOLOGY The present study has been performed by means of MD simulations using the DL_POLY program,33 where the potential energy function describing the H2S−H2S intermolecular interaction and its derivatives has been implemented. In order to test the reliability of the potential model, ab initio calculations of the small (SH2)2,3 clusters have been carried out. In this section, the formulation of the potential energy function and details of both ab initio calculations and MD simulations are presented. 2.1. H2S−H2S Intermolecular Interaction. The total intermolecular interaction, Vint, is formulated as a sum of two effective and nearly independent components, which are representative of the leading terms and also indirectly including additional effects coming from less important contributions. Specifically, the interaction between permanent charge distributions, the so-called electrostatic energy (Vel), evaluated as a sum of Coulomb potentials between point of charges, is combined with Vnel, globally including the remaining nonelectrostatic components. In our potential model, a crucial aspect in the formulation of Vnel, which results from the combination of a size repulsion component with an induction and dispersion attraction one, is based on the choice of the interaction centers, localized on selected groups (bonds or other groups of atoms; see for instance refs 34 and 35) of the interacting molecules. Each interaction center has an assigned value of the polarizability, which is considered the key property to scale both attraction and repulsion in molecular systems,29 and it is decomposed as bond36 or other groups contributions (see for instance ref 37). Accordingly, the polarizability decomposition must be compatible with the total polarizability and, for a given molecule, the sum of the polarizability values assigned to each group of atoms must be equal to the total molecular polarizability. In the present study, following the same procedure used to describe the intermolecular interactions involving H2O or NH3 and considering also the limited dimension of H2S, a unique interaction center placed on the S atom is considered. The polarizability assigned to S is equal to 3.631 Å3,38 which is that of the H2S molecule. As a matter of fact, for small molecules, the consideration of only one interaction center, localized on the leading atom, for instance, the O atom of H2O 25,37 and the N atom of NH3,27 has been found to be adequate to describe the total system. Accordingly, for H2S−H2S, Vnel is defined as the interaction between two centers placed at an intermolecular distance and localized on the two S atoms. It is represented by means of an improved Lennard-Jones function (VILJ), where B

DOI: 10.1021/acs.jpca.5b11843 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A ⎡ ⎛ r0 ⎞n(r) n(r ) ⎛ r0 ⎞m⎤ m ⎜ ⎟ ⎜ ⎟ ⎥ − VILJ(r ) = ε⎢ ⎢⎣ n(r ) − m ⎝ r ⎠ n(r ) − m ⎝ r ⎠ ⎥⎦

phase (0.97 D).38 This allows to take indirectly into account for polarization effects without increasing the computation time. An increase of about a 67% (μ = 1.620 D) has been applied to study the (H2S)2 system. Such value is suggested from the values of the molecular polarizability of H2O, NH3, and H2S, equal to 1.501 Å3,38 2.16 Å3,40 and 3.631 Å3,38 respectively, and from the increases of approximately 13% and 36% of the gas phase dipole moments applied to investigate water and ammonia dimers, respectively. As it was made in our preceding studies of water and ammonia, higher increases of the dipole moment have been used to study (H2S)3 (μ = 1.750 D) and liquid hydrogen sulfide (μ = 1.843 D). The value of 1.843 D is in the range of other values used in different potential energy models. With the exception of the polarizable Drude potential model,19 in which the dipole moment of the molecule in gas phase equal to 0.98 D is used, other nonpolarizable models use values of 1.4 D,16,17 1.7 D,18 and 2.1 D.15 However, it must be stressed again that, in our model, different values of the dipole moment can be used in combination with different values of β.25 2.2. Ab Initio Calculations. Ab initio calculations on the monomer, dimer, and trimer have been performed to obtain their equilibrium structures and the binding energies of the (H2S)2 and the (H2S)3 systems. Geometry optimizations have been computed, departing from the MD equilibrium structures, using MP241,42 and aug-cc-pVDZ basis set using the Gaussian 09 suite.43 On the optimized structures, single point coupled cluster CCSD(T)44−46 calculations using the aug-cc-pVTZ basis set have been performed including the counterpoise correction47 for the basis set superposition error (BSSE) in the intermolecular interaction. 2.3. Molecular Dynamics Simulations. MD simulations for the small clusters of H2S have been performed, without imposing boundary conditions, by considering a microcanonical ensemble (NVE), for which the number of particles, N, the volume, V, and the total energy (kinetic + potential), E, are conserved along the trajectory. Due to the conservation of the total energy, for very low values of the kinetic energy (at very low temperatures) one can say that the potential energy, Ecfg, remains almost constant along the trajectory at a given temperature. In fact, the near constancy of the potential energy along a MD simulation can be observed by analyzing the energy values at each step (Ecfg,i) of the trajectories run at low temperatures. Accordingly, if the configuration of the cluster is close to the equilibrium one, the extrapolation of the potential energy to 0 K should give a value very similar to that of the binding energy. To avoid working with a geometry corresponding to a local minimum, we have performed calculations at increasing T so that the system can surmount the barrier between the local minimum and the stablest configuration. Then, by diminishing T and, therefore, the available energy, the system tends to the equilibrium configuration and the extrapolation of the Ecfg values to 0 K provides an estimation of the binding energy. Moreover, from the analysis of the trajectories performed at the lowest temperature for different systems,27,48 it has been observed that the extrapolated value is very similar to the lowest step Ecfg,i value attained along the trajectory, always corresponding to a step with a kinetic energy value close to zero. For liquid hydrogen sulfide, instead, we are interested in the calculation of the density values under different conditions of p and T. In this case, the possibility of the expansion or contraction of the simulation cell is necessary. Such study

(1)

with ⎛ r ⎞2 n(r ) = β + 4⎜ ⎟ ⎝ r0 ⎠

(2)

The relevant parameters of the ILJ function, r0 and ε, have been calculated following the guidelines given by Cambi et al.29 In particular, the r0 parameter for the H2S−H2S interaction, equal to 4.253 Å, has been calculated from the following equation, r0 = 2C

(αH2S)1/3 (αH2S)2γ

(3)

where αH2S is the polarizability of H2S. In eq 3, the C coefficient and the γ exponent are equal to 1.767 and 0.095, respectively, to obtain r0 in Å. Such values have been derived by averaging different values corresponding to well-known systems.29 On the other hand, ε, equal to 1.915 kJ mol−1 for the H2S−H2S interaction, has been calculated by considering the near proportionality between the C6,eff coefficient and r0 as ⎛ C6,eff ⎞ ε = 0.720⎜ 6 ⎟ ⎝ r0 ⎠

(4)

The C6,eff coefficient is an effective long-range interaction term, including dipole−dipole, dipole−multipole, and multipole−multipole interactions, which can be determined from scattering experiments.29 Such experiments have proved the dependence of C6,eff on the polarizability of the interacting particles. In fact, C6,eff, for the H2S−H2S interaction, can be determined from the following equation, K (Neff )1/2 (αH2S)3/2 (5) 2 where Neff is the number of effective electrons and K = 1513 to obtain C6,eff in kJ mol−1 Å6. Neff can be obtained from ab initio calculations by determining the relative perturbation induced by a point charge on the various atomic orbitals. The ILJ function, because of the additional parameter β, becomes more flexible and realistic with respect to the LJ one.28 Without modifying the relevant parameters of the potential (the well depth, ε, and the distance of the minimum of energy, r0), it can be applied to describe the interacting system in different environments. In fact, n(r) allows us to model the dependence on r of the steepness of the repulsion and also of the strength of the interaction as balancing of repulsion and attraction.39 This means that, without modifying ε and r0, the β modulation allows us not only to include minor effects from other components emerging at intermediate and short-range as the stabilization from charge transfer in the perturbative limit but also to partially compensate the use of different representations of the charge distributions defining Vel. As for water−water interactions,23−25,32 values of β = 6.7 and m = 6 have been considered. In the present study Vel has been evaluated by considering a three-point charge distribution on the H2S molecules, as it was performed for H2O−H2O23−25,32 and NH3−NH3.27 Also in the present evaluation of Vel in H2S clusters the dipole moment, μ, has been increased with respect to that of the monomer in gas C6,eff =

C

DOI: 10.1021/acs.jpca.5b11843 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A Table 1. Main Distances and Angles from the Dimer (H2S)2 and the Trimer (H2S)3a ab initio (H2S)2 (H2S)3

ILJ

rSH

rSS

α

θ

Ebind

rSH

rSS

α

θ

Ebind

2.745 2.716

4.079 4.008

173.3 161.0

95.2 86.6

−6.39 −19.87

2.720 2.683

3.989 3.900

176 154

105 103

−6.38 −19.65

a

For the trimer values are averaged. Paremeters listed are hydrogen bond length (rSH), distance between S atoms (rSS), SHS angle (α), angle between the donor SH axis and the C2 axis of the acceptor H2S molecule (θ), and binding energy Ebind. Distances, angles, and Ebind are given in Å, deg, and kJ mol−1, respectively.

requires, therefore, a NpT ensemble that mimics closely the laboratory conditions in a flask, where variations in the volume are retrieved when T and/or p change. Temperature and pressure are controlled using the Nosé−Hoover method49 and applying the Berendsen algorithm,50 respectively. In almost all simulations we have considered an ensemble of 512 molecules which, in agreement with the value of the density at the boiling point, have been randomly distributed in a cubic box of 31.25 Å of length. Moreover, in order to approximate the case of 512 H2S molecules to a large (infinite) system, cubic periodic boundary conditions have been imposed. MD simulations have been performed at increasing temperatures, in the 200−300 K range, and at a pressure of 1.013 bar. In addition, also investigated was the behavior of H2S at the critical point pressure (pc = 90 bar) and at temperatures in the 350−400 K range (including the critical value of 373.1 K) and at some selected thermodynamic states. The study has been completed with the investigation of the liquid/vapor interface where, since a high number of molecules is required to avoid system size effects,51 an ensemble of 2197 molecules has been used. For this system, the densities at two temperature/ambient pressure conditions have been computed. Simulations along 5 ns (including 1 ns of equilibration), using a cutoff radius of 12.5 Å and a time step of 1 fs, have been performed. Once the system is equilibrated (first 1 ns), the probed potential and kinetic energy values (Ecfg,i and Ek,i), the instantaneous temperature, Ti, and the evaluated atom positions have been collected along the remaining 4 ns of the trajectory. At the end of the simulation, the temperature and energy have been averaged to obtain the corresponding mean values, represented by Ecfg, Ek, and T. Bearing in mind that the geometry of H2S in the dimer is almost unchanged from that of the monomer,52 all the simulations have been performed by considering H2S as a rigid molecule.

at decreasing temperatures, which has been used previously to estimate binding energies for different aggregates, can be regarded as a simulated annealing minimization.53 At very low temperatures, even when several isomers exist, the system has not enough energy to isomerize, and a linear dependence of Ecfg on T can be expected. Sometimes the linear variation is maintained even at high temperatures, thus indicating either that isomerization does not occur or that the involved isomers have very similar energies.37,54 The Ecfg values obtained at different temperatures for (H2S)2 and (H2S)3, are represented in the top and lower panels of Figure 1, respectively.

Figure 1. Mean values of the potential (configuration) energies at the range of temperatures investigated. Top panel: (H2S)2. Lower panel: (H2S)3.

As can be observed from Figure 1, the potential model predicts binding energies of −6.38 kJ mol−1 and −19.65 kJ mol−1 for the dimer and the trimer, respectively. Such results are in excellent agreement with ab initio results. Moreover, the value of −6.38 kJ mol−1 for the dimerization energy is in very good agreement with the value of −6.53 kJ mol−1 reported by Kristóf and Liszi,17 our dimer being slightly more stable than that predicted by the polarizable Drude potential model, equal to −4.10 kJ mol−1.19 Other potential models predict a higher stability of the dimer, with dimerization energies equal to −7.37 kJ mol−1,18 −7.24 kJ mol−1,16 and −10.80 kJ mol−1.15 Reinforcing the validity of the binding energy estimations, the analysis of the instantaneous values of the potential energy, Ecfg,i, allows us to see that for the dimer the lowest energy value at the lowest temperature investigated is equal to −6.38 kJ mol−1 and for the trimer it is equal to −19.66 kJ mol−1, which are very similar to the extrapolated values shown in Figure 1. From the previous results, it emerges that the structures associated with the lowest Ecfg,i should be very similar to those

3. (H2S)2 AND (H2S)3 SYSTEMS The ab initio optimizations performed show that in the dimer, each molecule acts as H-donor or as H-acceptor only, with one of the H of the first molecule pointing to the S of the second one. Contrarily, in the trimer each molecule plays simultaneously the role of donor and acceptor, forming an equilateral triangle structure between three sulfurs and three intermolecular bonds. The relevant distances and angles of (H2S)2 and (H2S)3 together with the corresponding binding energies obtained from present ab initio calculations are given in the left side of Table 1. To estimate the binding energy, from MD simulations, trajectory calculations at temperatures ranging from 5 to 80 K have been carried out. Then, the mean potential energy values of the configurations attained along the trajectories, Ecfg(T), have been represented as a function of T and extrapolated to 0 K to make an estimation of Ebind.26,37 For weakly bound small aggregates, the process of running consecutive MD simulations D

DOI: 10.1021/acs.jpca.5b11843 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Figure 2. MD simulations on the (H2S)2 dimer performed at two different values of T: (a) evolution of the four intermolecular SH distances, rSH; (b) instantaneous values of the potential energy, Ecfg,i; (c) evolution of one of the four intermolecular SH distances, rSH; (d) instantaneous values of the potential energy, Ecfg,i. Temperatures are indicated in each panel.

The situation of no interexchange of the role as a donor or acceptor of the H atom is shown in the panel a of Figure 2, where the evolution of the four intermolecular SH distances, rSH, along the 4 ns of simulation obtained at 9 K is represented. Accordingly, no noticeable variations on the dimer configuration occurs at 9 K, and consequently, the rSH distances and the values of Ecfg,i remain almost constant along the trajectory, as can be seen in panels a and b of Figure 2. Panel a of Figure 2 also shows that the mean value of the shortest rSH distance fluctuates around 2.7 Å, this value being in very good agreement with that of 2.745 Å derived from ab initio calculations. In addition, it has been found a mean value of the distance between the S atoms equal to 4.0 Å, comparable with the value of 4.079 Å given in Table 1. Contrarily, when increasing T, the H2S molecules can interexchange their role of H-acceptor and donor along the trajectory, originating variations in the values of Ecfg,i. In fact, both Cs linear hydrogen bonded and C2v bifurcated dimer structures are observed along the trajectory. A detailed analysis of both the geometries and the energies, Ecfg,i, allows us to conclude that, in agreement with several ab initio calculations, the Cs structure is open and more stable than the C2v bifurcated one.56 At T = 40 K, noticeable fluctuations of one of the rSH distances (the remaining three rSH distances evolve in a similar way) are observed (see panel c of Figure 2). However, as shown in the panel d of Figure 2, the pronounced fluctuation of the Ecfg,i energy does not allow us to individuate the two structures. The regular fluctuation of Ecfg,i along the trajectory confirms that the isomerization process takes place between two equivalent isomers and that the C2v bifurcated structure should correspond to a transition state. Regarding the hydrogen sulfide trimer (SH2)3, similar to the structures found for water23 and ammonia,27 a triangular structure in which every molecule acts as a donor and acceptor is the most stable. Moreover, it is very compact and isomerizations are unlikely. In fact, only small variations on the geometry (which originate a gradual increase of the potential energy (less negative values) are observed along the trajectories, even at the highest T investigated.

at the equilibrium. Geometric characteristics of such equilibrium-like structures are given in the right side of Table 1 compared to the ab initio ones. With regard to the geometry, in agreement with experimental data52,55 and other theoretical calculations56,57 the minimum energy of the dimer corresponds to an almost linear hydrogen bond. Both distances and angles are in agreement with the values reported in the literature.57 In particular, for the dimer, our potential model predicts a rSS distance equal to 3.899 Å, very similar to the value of 3.90 Å reported by Kristóf and Liszi17 and slightly shorter than that calculated using the polarizable Drude potential model, equal to 4.06 Å.19 Values of 3.83 and 3.89 Å were obtained by Forester et al.16 and by Potoff and co-workers.18 As can be expected, from the value of the dimerization energy, the lowest distance, equal to 3.62 Å, was obtained using the model of Jorgensen.15 On another hand, the IR spectra of H2S molecules in Kr solids assigns several bands of the IR spectrum to the hydrogen bond stretching of the trimer.52 In the same study, using DFT, an optimized structure for the trimer in which the three S atoms approximately form an equilateral triangle is also obtained, in agreement with our calculations. The rSH and rSS distances for both dimer and trimer are very similar to those obtained from ab initio calculations, the highest difference being related to the θ angle. However, calculations show that geometries differing only in the θ angle exhibit very similar energies, according to the fluxional nature of the dimer. In fact, as it was observed for water and ammonia dimers, H2S− H2S is a fluxional H-bonded cluster with the propensity to easily interexchange the donor/acceptor nature of the two molecules. The equivalent isomers, only differing in the H-bond donor and acceptor roles played by the two molecules, are interconnected through very small barriers, and isomerization can take place even at very low temperatures. At T below 12 K no isomerizations have been observed. This means that, at these temperatures, along the trajectory the H2S molecules in the dimers do not interexchange their role of H-donor or acceptor along the trajectory. E

DOI: 10.1021/acs.jpca.5b11843 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

4. LIQUID HYDROGEN SULFIDE A first batch of trajectories, considering a NpT ensemble of 512 molecules, has been run at ambient pressure and temperatures in the 180−240 K range, which includes melting (187.4 K) and boiling (212.9 K) values. A very good agreement of our result of the enthalpy of vaporization with the experimental one of 18.68 kJ mol−1 58 has been also found. As a matter of fact, at 212.9 K and 1 atm of pressure our model gives an enthalpy of vaporization of 19.23 kJ mol−1, differing only about a 3% from the experimental value. All the potential models for hydrogen sulfide analyzed in ref 19 with the exception of that of Kristóf and Liszi17 show larger deviations from the experimental value of the enthalpy of vaporization. Moreover, the predicted value of the relative permittivity equal to 8.10, calculated using the standard expression given in ref 59, is also in agreement with the experimental value of 8.04.60 The evolution of density and diffusion coefficients with T has been analyzed, and as shown in Figure 3, a gradual variation on the density (d) and on the

As can be seen in Figure 4, results show a moderate density decrease and a diffusion coefficient increase with T up to

Figure 4. Density and diffusion coefficients at 90 bar (the critical pressure) of pressure as a function of T (around the critical temperature) in the top and lower panels, respectively.

around 375 K. The changes become much more pronounced when T is further increased, indicating that at 375 K we are close to the critical point. Moreover, at temperatures around 390 K, the system suddenly expands, which can be attributed to the liquid to gas phase transition. In order to further analyze the gas phase transition, MD simulations for a NpT ensemble of 2197 molecules have been performed at 1.013 bar of pressure and at the temperatures of 212.9 and 235 K. To prepare the initial configuration in agreement with the density at 212.9 K, 2197 molecules of H2S have been randomly distributed into a cubic box (Lx = Ly = Lz) of 50.785 Å of length. A NVE ensemble has been used to thermalize the system, and its final configuration, placed in the middle of a simulation cell with the z size increased (Lx = Ly= 50.785 Å and Lz = 150.785 Å), has been chosen as initial configuration to perform NpT MD simulations at T = 212.9 and 235 K and p = 1.013 bar. Such setting up allows us to obtain two vapor regions with a liquid slab in between.51 To analyze the results obtained, the simulation shell has been divided into slabs of 2.0 Å of thickness along the z axis and, every 50 time steps, the local density profile has been computed for each slab as explained in ref 62. Results are reported in Figure 5 where the plateau region corresponds to the density of a bulk hydrogen sulfide.

Figure 3. Density and diffusion coefficients at 1.013 bar of pressure as a function of T in the top and lower panels, respectively. The red mark in the top panel represents the experimental value of the density at the boiling point.

diffusion coefficients (D) is observed when the temperature is increased. At the boiling point a density of 0.946 g cm−3 has been computed, in very good agreement with the experimental value of 0.949 g cm−3. At the same thermodynamic conditions, we have calculated a diffusion coefficient equal to 4.3 × 10−9 m2 s−1, similar to the experimental values of 3.7 × 10−9 m2 s−1, and 4.6 × 10−9 m2 s−1 obtained at the temperatures of 206.5 and 223.2 K, respectively, along the coexistence curve.61 The behavior shown in Figure 3 can be explained because when T increases, the molecules tend to be placed far away, originating a slight increase of the volume and, accordingly, a decrease of the density. Conversely, due to the increase of T, the kinetic energy of the molecules is higher and their diffusion becomes easier, thus originating an increase of D. Moreover, a decrease of about a 13% in the configuration energy, ascribed to both changes on the nonelectrostatic and electrostatic interaction energy contributions, is observed, indicating that less attractive configurations of the molecules can be probed. The diffusion and density behaviors have also been investigated at conditions close to the critical ones. Several MD simulations at the critical pressure of 90 bar have been performed at increasing values of T, starting from T = 350 K.

Figure 5. Density profiles for hydrogen sulfide at 212.9 K (black markers) and 235 K (green markers). The red line corresponds to the fit of eq 6 at 212.9 K. F

DOI: 10.1021/acs.jpca.5b11843 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A The density profile can be fit with a hyperbolic tangent function:62,63 ⎡ A(z − z 0) ⎤ d(z) = 0.5(dl + d v ) − 0.5(dl − d v )tanh⎢ ⎥ ⎣ t inter ⎦ (6)

where dl and dv are the densities corresponding to the liquid and vapor phases, z0 is the position of the Gibbs dividing surface, and tinter is the interface thickness, defined as the thickness over which the density of hydrogen sulfide varies from a 90% to a 10% of dl. The position of the Gibbs surface, the interface thickness, and dl are given in Table 2.

Figure 6. S−S radial distribution function at 293.15 K and 30.40 bar. The g(rSS) function and n(rSS) are represented by continuous black and red lines, respectively. The dashed blue line defines the area used in the integration of g(rSS) to derive the coordination number (see text).

Table 2. Results of dl, z0, and tinter Derived from Fitting the Density Profiles with Equation 6 at Temperatures of 212.9 and 235 K T, K

dl, g cm−3

z0 , Å

tinter, Å

212.9 235.0

0.910 0.864

26.6 26.9

9.6 11.9

represented in Figure 6, where the coordination number is also indicated. However, it must be indicated that the integration of the first peak in the standard way gives a value of 11.9. Finally, the g(rSH) and g(rHH) functions have been also computed and compared with neutron diffraction65 experimental data. In Figure 7, the S−H and H−H radial distribution

The density of the vapor phase takes values close to zero. In fact, the values of dv obtained using different potential models21 are very small (of the order of 3 × 10−3 g cm−3), while those obtained for dl are slightly higher (about 0.89 g cm−3) than the values reported in Table 2. However, it must be taken into account that the values of Table 2 are derived from a fit of eq 6. As it can be observed by comparing the fitted value of 0.910 g cm−3 with the calculated value of 0.949 g cm−3 at 212.9 K, fitted values are underestimated. In spite of this, such values allow us to see that the percent of molecules involved in the phase transition increase with T. Results indicate that about a 10% of molecules pass from liquid to gas phase at 212.9 K, while at 235 K this percent increases to about a 18%. The capability of the potential model to reproduce the main structural characteristics of liquid hydrogen sulfide has also been investigated. The variation of the density and of the mean diffusion coefficients with T and the increase (less negative value) of the configuration energy suggest a change on the liquid structure. Such modification has been investigated by analyzing some radial distribution functions, g(r). Overall, small changes associated with the increase of T have been observed in the S−S radial distribution function, g(rSS), where it originates only a slightly broadening of the first peak. On another hand, Xray scattering data64 at 293.15 K and 30.40 bar have been compared to the g(rSS) function computed at the same p and T conditions. MD simulations predict a density of 14.5 × 10−3 molecules/Å3, in very good agreement with the experimental value of 14.0 × 10−3 molecules /Å3, while the first peak of the S−S in RDF is very similar to the experimental one, with the maximum at rSS = 4.0 ± 0.5 (0.1 Å below the experimental value).64 The main differences with X-ray diffraction data deal with the second peak which shows a double structure, associated with different microscopic configurations, not reproduced by MD simulations. Moreover, the integration of the first peak, following the same criteria that in ref 64, has been performed in a nonstandard way (see Figure 6) and we have calculated a coordination number, n(rSS), associated with the number of neighbors in the first shell, equal to 8, in good agreement with the value of 7.5 derived from X-ray experiments. The g(rSS) function (in black) and n(rSS) (in red) are

Figure 7. S−H (top panel) and H−H (lower panel) radial distribution functions at 208 K and 1.013 bar.

functions, obtained at the experimental thermodynamic conditions (T = 208 K, p = 1.013 bar), are reported in the top and lower panels, respectively. The g(rSH) and g(rHH) functions, in contrast with the equivalent ones for water,25 are quite broad. In fact, for water, the two radial distribution functions show two well-defined peaks (mainly in the O−H radial distribution function), such characteristic being absent in g(rSH) and in g(rHH) functions. In fact, regarding g(rSH) and in agreement with neutron diffraction experiments,65 the function shows a shoulder at about 2.8 Å and a sort of double peak centered around 4.5 Å. The position of the shoulder reflects the H-bond distance (rSH) reported for the dimer and the trimer in Table 1, while the large peak at 4.5 Å is compatible with the wide rSH distribution reported in Figure 2 that tends to broaden as T increases. As can be expected, in comparison with g(rSH), the g(rHH) function is less structured and no shoulder or double peaks have been observed. The analysis of the radial distribution functions and G

DOI: 10.1021/acs.jpca.5b11843 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

point, the values of the predicted density, the enthalpy of vaporization, and the relative permittivity are equal to 0.946 g cm−3, 19.23 kJ mol−1, and 8.10, respectively, and are in good agreement with the corresponding experimental values of 0.949 g cm−3, 18.68 kJ mol−1, and 8.04. Consistent with the experiments, MD simulations performed at the critical pressure and at a temperature very close to the critical value show a noticeable variation of both the density values and the mean diffusion coefficients. By increase of the size of the system, the density values at the interface liquid/vapor have been evaluated at two different temperatures, observing a decrease with T of the density of liquid at the interface. Structural changes when increasing T have been analyzed by considering the S−S, S−H, and H−H radial distribution functions. The results have been found to be in good agreement with X-ray64 and neutron diffraction experiments.65 It can be concluded that the potential model is able to describe in a consistent way the differences of the H-bonds networks of similar systems as those formed by water, ammonia, and hydrogen sulfide. Moreover, in view of the energy and structural predictions of the potential model it can be concluded that in spite of its simplicity, the present formulation of the effective H2S−H2S interaction potential opens the way toward the characterization of systems formed by hydrogen sulfide and several other molecules, and this could be useful in several fields concerning environment, safety, and health problems.

its comparison with those of water seem to indicate that liquid hydrogen sulfide is less hydrogen bonded. However, in agreement with the results of Santoli et al.,66 the characteristics of g(rSH) suggest that a weak degree of orientational correlations among first neighboring molecules may exist. Such results are very similar to those obtained for ammonia. In fact, by comparing the g(rXH) functions (X = O, N, S), it can be observed that the g(rOH) is much more structured than the g(rNH) and g(rSH) ones. While in the g(rOH) function two welldefined peaks appear, the g(rNH) and g(rSH) ones show a very similar broad and almost no structured peak, in agreement with the fact that the H-bond interaction in liquid ammonia67 is much weaker than in the case of water. The model predicts binding energies for the dimers of water, ammonia, and hydrogen sulfide to be equal to −21.99, −12.71, and −6.38 kJ mol−1, respectively, and bearing in mind that the three dimers are stabilized by the formation of one H-bond, more structured radial distribution functions can be expected for liquid water in comparison with ammonia and hydrogen sulfide. Therefore, considering the differences in the NH3−NH3 and the H2S−H2S interaction energies, liquid ammonia is a little bit more structured than hydrogen sulfide.

5. CONCLUSIONS Extensive MD simulations, performed by exploiting a new formulation of the intermolecular interaction, provide a proper and internally consistent rationalization of the basic properties of H2S in both gaseous and liquid phases and at the transition between them. Basically, the intermolecular interaction between H2S molecules has been constructed by assuming the separability of permanent (Vnel) and nonpermanent (Vel) interactions. Vnel has been formulated using an improved Lennard-Jones (ILJ) function and considering only one interaction center placed on the S atom. The parameters of the effective interaction potentials defining Vnel, that is, the well depth (ε) and the distance of the minimum energy (r0), calculated from the polarizability of the molecule, have been used to investigate systems of increasing complexity involving H2S. Polarization effects have been taken into account indirectly by increasing the value of the dipole moment of the molecule going from small clusters to liquid and modulating Vel. One of the main advantages of the model is that because of the modulation of β, the relevant parameters defining Vnel do not have to be changed when considering different environments like small clusters or condensed phases. To test the potential model predictions, the monomer, dimer, and trimer of hydrogen sulfide have been optimized at the MP2/aug-ccpVDZ level of calculation, and single point coupled cluster CCSD(T) using the aug-cc-pVTZ including the counterpoise correction for the basis set superposition error (BSSE) has been performed. MD simulations under different conditions have been carried out by considering rigid molecules. Results from MD calculations indicate that, in spite of the simplicity of the potential model, it is able to reproduce the main structural characteristics and the energetics of the small clusters. In fact, the predicted binding energies for the dimer and the trimer are in excellent agreement with the ab initio calculations, differing only by 0.02 kJ mol−1 and 0.20 kJ mol−1, respectively. Using a NpT ensemble of 512 molecules, MD simulations have been performed in a wide range of thermodynamic conditions. At ambient pressure, the dependence on T of the density and mean diffusion coefficients have been analyzed in an interval of T values around the boiling point (212.9 K). At the boiling

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpca.5b11843. Three DL_POLY files (used to perform the simulations at 212.9 K): CONFIG, CONTROL, FIELD (ZIP)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS A. Amat thanks the Italian Ministry of University and Research (MIUR) for the PRIN Contract PRIN-2010 20104XET32 ”DSSCX”. M. Alberti ́ and A. Aguilar are thankful for financial support from the Ministerio de Educación y Ciencia (Spain, Project CTQ2013-41307-P) and from the Generalitat de Catalunya (Grant 2009SGR-17). Also thanks are due to the Center de Supercomputació de Catalunya CESCA-C4 and Fundació Catalana per a la Recerca for the allocated supercomputing time. F. Pirani acknowledges financial support from the MIUR for the PRIN 2010−2011, Grant 2010 ERFKXL_002 PRIN.

■

REFERENCES

(1) Reiffenstein, R. J.; Hulbert, W. C.; Roth, S. H. Annu. Rev. Pharmacol. Toxicol. 1992, 32, 109−134. (2) O’Rourke, D.; Connolly, S. Annual Review of Environment and Resources 2003, 28, 587−617. (3) Predmore, B.; Lefer, D.; Gojon, G. Antioxid. Redox Signaling 2012, 17, 119−140. H

DOI: 10.1021/acs.jpca.5b11843 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A (4) Pryor, W.; Gojon, G.; Church, D. J. Org. Chem. 1978, 43, 793− 800. (5) Stoyanovsky, D.; Maeda, A.; Atkins, J.; Kagan, V. E. Anal. Chem. 2011, 83, 6432−6438. (6) Polhemus, D.; Calvert, J.; Butler, J.; Lefer, D. Scientifica 2014, 2014, 768607. (7) Zheng, Y.; Ji, X.; Ji, K.; Wang, B. Acta Pharm. Sin. B 2015, 5, 367−377. (8) López-Rendón, R.; Alejandre, J. J. Mex. Chem. Soc. 2008, 52, 88− 92. (9) Berendsen, H.; Postma, J.; van Gunsteren, W.; Hermans, J. Intermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, The Netherlands, 1981. (10) Berendsen, H.; Grigera, J.; Straatsma, T. J. Phys. Chem. 1987, 91, 6269−6271. (11) Jorgensen, W.; Chandrasekhar, J.; Madura, J.; Impey, R.; Klein, M. J. Chem. Phys. 1983, 79, 926−935. (12) Jorgensen, W.; Madura, J. Mol. Phys. 1985, 56, 1381−1392. (13) Mahoney, M.; Jorgensen, W. J. Chem. Phys. 2000, 112, 8910− 8922. (14) Wang, L.; Head-Gordon, T.; Ponder, J.; Ren, P.; Chodera, J.; Eastman, P.; Martinez, T.; Pande, V. J. Phys. Chem. B 2013, 117, 9956−9972. (15) Jorgensen, W. J. Phys. Chem. 1986, 90, 6379−6388. (16) Forester, T.; McDonald, I.; Klein, M. Chem. Phys. 1989, 129, 225−234. (17) Kristóf, T.; Liszi, J. J. Phys. Chem. B 1997, 101, 5480−5483. (18) Kamath, G.; Lubna, N.; Potoff, J. J. Chem. Phys. 2005, 123, 124505. (19) Riahi, S.; Rowley, C. J. Phys. Chem. B 2013, 117, 5222−5229. (20) Orabi, E.; Lamoureux, G. J. Chem. Theory Comput. 2014, 10, 3221−3235. (21) Shah, M.; Tsapatsis, M.; Siepmann, J. J. Phys. Chem. B 2015, 119, 7041−7052. (22) Albertí, M.; Castro, A.; Laganà, A.; Moix, M.; Pirani, F.; Cappelletti, D. Eur. Phys. J. D 2006, 38, 185−191. (23) Albertí, M.; Aguilar, A.; Bartolomei, M.; Cappelletti, D.; Laganà, A.; Lucas, J.; Pirani, F. Phys. Scr. 2008, 78, 058108. (24) Albertí, M.; Aguilar, A.; Cappelletti, D.; Laganà, A.; Pirani, F. Int. J. Mass Spectrom. 2009, 280, 50−56. (25) Faginas Lago, N.; Huarte-Larrañaga, F.; Albertí, M. Eur. Phys. J. D 2009, 55, 75−85. (26) Albertí, M.; Aguilar, A.; Huarte-Larrañaga, F.; Lucas, J.; Pirani, F. J. Phys. Chem. A 2014, 118, 1651−1662. (27) Albertí, M.; Amat, A.; Farrera, L.; Pirani, F. J. Mol. Liq. 2015, 212, 307−315. (28) Pirani, F.; Brizi, S.; Roncaratti, L.; Casavecchia, P.; Cappelletti, D.; Vecchiocattivi, F. Phys. Chem. Chem. Phys. 2008, 10, 5489−5503. (29) Cambi, R.; Cappelletti, D.; Liuti, G.; Pirani, F. J. Chem. Phys. 1991, 95, 1852−1861. (30) Albertí, M.; Aguilar, A.; Lucas, J.; Pirani, F.; Cappelletti, D.; Coletti, C.; Re, N. J. Phys. Chem. A 2006, 110, 9002−9010. (31) Albertí, M.; Aguilar, A.; Lucas, J.; Pirani, F.; Coletti, C.; Re, N. J. Phys. Chem. A 2009, 113, 14606−14614. (32) Albertí, M.; Aguilar, A.; Bartolomei, M.; Cappelletti, D.; Laganà, A.; Lucas, J.; Pirani, F. Lecture Notes in Computer Science 2008, 5072, 1026−1035. (33) Todorov, I.; Smith, W.; Trachenko, K.; Dove, M. J. Mater. Chem. 2006, 16, 1911−1918. (34) Albertí, M. J. Phys. Chem. A 2010, 114, 2266−2274. (35) Albertí, M.; Costantini, A.; Laganà, A.; Pirani, F. J. Phys. Chem. B 2012, 116, 4220−4227. (36) Pirani, F.; Cappelletti, D.; Liuti, G. Chem. Phys. Lett. 2001, 350, 286−296. (37) Albertí, M.; Faginas-Lago, N.; Laganà, A.; Pirani, F. Phys. Chem. Chem. Phys. 2011, 13, 8422−8432. (38) NIST Computational Chemistry Comparison and Benchmark Database. NIST Standard Reference Database Number 101. http:// cccbdb.nist.gov/.

(39) Pirani, F.; Albertí, M.; Castro, A.; Moix Teixidor, M.; Cappelletti, D. Chem. Phys. Lett. 2004, 394, 37−44. (40) Linstrom, P.; Mallard, W. NIST Chemistry WebBook, NIST Standard Reference Database Number 69; National Institute of Standards and Technology: Gaithersburg, MD; http://webbook.nist. gov. (41) Head-Gordon, M.; Pople, J.; Frisch, J. Chem. Phys. Lett. 1988, 153, 503−506. (42) Head-Gordon, M.; Head-Gordon, T. Chem. Phys. Lett. 1994, 220, 122−128. (43) Frisch, M. J.; et al. Gaussian 09, revision D.01; Gaussian Inc.: Wallingford, CT, 2009. (44) Purvis, G., III; Bartlett, A. J. Chem. Phys. 1982, 76, 1910−1918. (45) Scuseria, G.; Janssen, C.; Schaefer, H., III J. Chem. Phys. 1988, 89, 7382−7387. (46) Pople, J.; Head-Gordon, K.; Raghavachari, M. J. Chem. Phys. 1987, 87, 5968−5975. (47) Boys, S.; Bernardi, F. Mol. Phys. 1970, 19, 553−566. (48) Albertí, M.; Amat, A.; De Angelis, F.; Pirani, F. J. Phys. Chem. B 2013, 117, 7065−7076. (49) Nosé, S. A. Mol. Phys. 1984, 52, 255−268. (50) Berendsen, H.; Postma, J.; Vangunsteren, W.; Dinola, A.; Haak, J. J. Chem. Phys. 1984, 81, 3684−3690. (51) Trokhymchuk, A.; Alejandre, J. J. Chem. Phys. 1999, 111, 8510− 8523. (52) Tsujii, H.; Takizawa, K.; Koda, S. Chem. Phys. 2002, 285, 319− 326. (53) Kirkpatrick, S.; Gelatt, C., Jr.; Vecchi, M. Science 1983, 220, 671−680. (54) Albertí, M.; Amat, A.; Aguilar, A.; Huarte-Larrañaga, F.; Lucas, J.; Pirani, F. Theor. Chem. Acc. 2015, 134, 161(1)−161(12). (55) Odutola, J.; Viswanathan, R.; Dyke, T. J. Am. Chem. Soc. 1979, 101, 4787−4792. (56) Sodupe, M.; Oliva, A.; Bertran, J. J. Am. Chem. Soc. 1995, 117, 8416−8421. (57) Bhattacherjee, A.; Matsuda, Y.; Fujii, A.; Wategaonkar, S. ChemPhysChem 2013, 14, 905−914. (58) Clarke, E.; Glew, D. Can. J. Chem. 1970, 48, 764−775. (59) Yoshii, N.; Miura, S.; Okazaki, S. Chem. Phys. Lett. 2001, 345, 195−200. (60) Havrilaik, S.; Swenson, R.; Cole, R. J. Chem. Phys. 1955, 23, 134−135. (61) Dupré, F.; Piaggesi, D.; Ricci, F. Phys. Lett. A 1980, 80, 178− 180. (62) Taylor, R.; Shields, R. J. Chem. Phys. 2003, 119, 12569−12576. (63) Townsend, R.; Gryko, J.; Rice, S. J. Chem. Phys. 1985, 82, 4391− 4392. (64) Andreani, C.; Petrillo, C.; Rocca, D. Europhys. Lett. 1988, 5, 145−149. (65) Andreani, C.; Merlo, V.; Ricci, M.; Soper, A. Mol. Phys. 1991, 73, 407−415. (66) Santoli, G.; Bruni, F.; Ricci, F.; Ricci, M.; Soper, A. Mol. Phys. 1999, 97, 777−786. (67) Ricci, M.; Nardone, M.; Ricci, F.; Andreani, C.; Soper, A. J. Chem. Phys. 1995, 102, 7650−7655.

I

DOI: 10.1021/acs.jpca.5b11843 J. Phys. Chem. A XXXX, XXX, XXX−XXX