Computational Studies of Structure and Dynamics of Clathrate

Ind. Eng. Chem. Res. 2007, 46, 131-142

131

Computational Studies of Structure and Dynamics of Clathrate Inhibitor Monomers in Solution Diego A. Go´ mez Gualdro´ n, Santiago Aparicio-Martı´nez,† and Perla B. Balbuena* Artie McFerrin Department of Chemical Engineering, Texas A&M UniVersity, College Station, Texas 77843

Theoretical studies of five monomer constituents of gas clathrate inhibitorssvinylpyrrolidone, vinylvalerolactam, L-proline (LP), 1-formylpyrrolidine, and (dimethylaminoethyl)methacrylate (DMAEMA)sare performed in the presence of water, with the solvent treated both explicitly and as a continuum. The binding energy increases as the inhibitor is more polarized, but no clear relationship is found between polarization and free energy of solvation values. Comparison of the binding energies of the inhibitor-water complexes and analysis of the hydration shells around the inhibitor sites indicate that short-range interactions influence the clathrate inhibitor effect more than long-range ones do. The balance of hydrophobic and hydrophilic moieties is found to be a determinant factor in the behavior of the inhibitor. In contrast to the other monomers, analyses of the DMAEMA and LP behaviors show that the N site can be associated with the hydrophilic moiety. Besides, DMAEMA has the best hydrophobic solvation around its methyl groups, resembling clathrate formation around methane. The time evolution of the mean square displacement and the power spectrum calculated from velocity autocorrelation functions reveal slight changes in the dynamic behavior of water that are not strong enough to be considered decisive in the inhibition efficiency. 1. Introduction Gas clathrates are icelike structures occurring in gas-water systems, where water forms a crystal trapping small molecules of a gas such as methane or ethane.1 High pressure or low temperature, as often found during natural gas exploration and processing and in transportation pipelines and processing equipment, favors the formation of clathrates, causing flow assurance problems and safety risks; therefore, there is great interest in using clathrate inhibitors to avoid the blockage of pipelines.2 The most used kind of clathrate inhibitor is the thermodynamic one, which shifts the boundaries of clathrate stability in the phase diagram for the system gas-waterclathrate toward more extreme conditions of temperature and pressure, away from the operation conditions of the systems of interest. A common thermodynamic clathrate inhibitor is methanol used in concentrations around 40% by weight. However, such an amount of methanol is expensive and also harmful for the environment.3 A promising approach to the clathrate problem is the use of low-dosage “kinetic” inhibitors,4 which consist mostly of organic polymers whose structural units are molecules containing both oxygen and nitrogen. These inhibitors act through a kinetic mechanism in which the onset of the hydrate formation is delayed. Previous theoretical studies using molecular dynamics (MD) and Monte Carlo simulations showed that monomers are sufficient to describe the inhibitor effect,5 illustrating the adsorption of the inhibitor on large clathrate surfaces; most of these works have been related to poly(vinyl pyrrolidone) (PVP),6 poly(vinyl caprolactam) (PVcap), and their monomers or similar molecules.7,8 Aparicio-Martinez et al.9 studied the properties of some linear and cyclic amides in solution and their possible use as kinetic inhibitors. In this work we are interested in analyzing the behavior of monomers of some experimentally tested kinetic inhibitors. We * To whom correspondence should be addressed. Tel.: 979-8453375. Fax: 979-845-6446. E-mail: [email protected]. † Permanent address: Department of Chemistry, University of Burgos, 09001 Burgos, Spain.

focus on the interactions of the inhibitors with water, using density functional theory (DFT) and classical MD simulations. (Dimethylaminoethyl) methacrylate (DMAEMA) and vinylpyrrolidone (VP) are constituents of the commercial inhibitor VC713, which has proved to be more efficient than PVP, one of the first kinetic inhibitors available. L-Proline (LP) is an amino acid involved with antifreeze proteins10 and has been tested experimentally together with 1-formylpyrrolidine (1FP) and vinylvalerolactam (Vval).11 Previous MD-based reports about the changes of water structure induced by VP had some limitations because of the use of some constraints on the VP molecule or the use of implicit atoms. We use an all-atom force field to investigate monomers different than VP and attempt a comparison among them. Special emphasis is given to the interaction of hydrophilic and hydrophobic moieties with water as well as to the alteration of dynamic properties of water, analyzing how these are related to the inhibition mechanism. DFT calculations are used to investigate the effect of the electronic structure on the behavior of the inhibitor molecules. The solvent is simulated as a continuum by use of the self-consistent reaction field (SCRF) method,12 and short-range interactions are studied by analyzing monomer-water complexes via DFT and SCRF methods, whereas the classical MD approach is used to incorporate a much larger number of water molecules describing short- and long-range interactions at a less expensive computational cost. 2. Computational Methods 2.1. Density Functional Theory Calculations. Gaussian 0313 was utilized for all the ab initio computations. Full optimizations were carried out for the molecules of interest, using DFT with the three-parameter Becke gradient-corrected exchange functional14 in conjunction with the Lee-Yang-Parr correlation functional (B3LYP).15,16 In a few cases, for instance in some LP conformations, the less expensive Hartree-Fock (HF) method was required because of convergence problems with the DFT calculations. Due to strong electrostatic interactions present either under the SCRF or due to hydrogen bonds with

10.1021/ie061123g CCC: $37.00 © 2007 American Chemical Society Published on Web 12/08/2006

132

Ind. Eng. Chem. Res., Vol. 46, No. 1, 2007

Figure 1. Top and side views of optimized gas-phase structures for the five monomers; the carbonyl groups always are in the same plane with several other atoms in each molecule. The nomenclature for the main carbon atoms is also displayed. The following color code is used: blue, N atoms; red, O atoms; dark gray, C atoms; light gray, H atoms. Greek letters are used to designate the position of the C atoms starting with R as the first C after the C attached to the carbonyl group, except in 1FP where we use the IUPAC nomenclature (numbers instead of Greek letters), because in this case the first C is attached to a N atom of the ring.

explicit water molecules, it is important to use a large and flexible basis set that describes properly electrons close and far from the nuclei; here the 6-311++g** basis set was utilized to include diffuse functions and polarization effects. Atomic charges were calculated to fit the electrostatic potential17 according to the Merz-Singh-Kollman (MK) scheme;18 these site charges were used in classical MD simulations. For aqueous-phase calculations, the SCRF approach was utilized with the solvent treated as a continuum. The polarizable continuum model (PCM)12 was used to obtain quantitative estimations of the different electrostatic and nonelectrostatic contributions to the total free energy of solvation. In the PCM approach the solute is located in a cavity defined as interlocking van der Waals spheres centered at the nuclei, and a reaction field is added to simulate the solvent effect. In this work, the united atom model is used to build the cavity, a value of 1.2 was used to scale all the radii, and 60 tesserae were used to divide the spherical surfaces; the contribution of electrostatic and nonelectrostatic terms used to calculate free energies is described in our previous work.9 To study short-range interactions, 1:1 and 2:1 complexes of water and monomer were optimized; several initial configurations were tested to evaluate possible interactions of water with the hydrophobic and hydrophilic moieties of each monomer. The continuum approach was used again to calculate in solution the best complexes obtained in gas phase. The binding energies were calculated and only the largest one for each inhibitor was corrected by use of the basis superposition error (BSSE) and the counterpoise procedure proposed by Boys and Bernardi,19 since the main purpose is comparative and the procedure appears not to alter the order of going from the inhibitor with the highest binding energy to the one with the lowest energy for any kind of interaction. 2.2. Classical Molecular Dynamics. MD simulations in the canonical NVT ensemble20 were carried out with the DL_POLY

program.21 Simulations were performed at 298.15 K with the temperature controlled through the Nose´-Hoover thermostat22 with temperature scaled every 15 fs. No significant differences were found in the structural results obtained in the range of 268.15-298.15 K. The equations of motion were solved by use of the Verlet Leapfrog20 integration algorithm with a Verlet width shell of 1.0 Å. The molecular geometries were constrained according to the SHAKE algorithm.23 Long-range electrostatic interactions are treated with the smooth particle mesh Ewald method.24 The simulation cell is a cubic box of ∼18.4 Å side containing one inhibitor monomer and 200 water molecules, yielding a density ∼0.99 g cm-3 for the whole system; periodic boundary conditions are applied in the three directions to simulate an infinite system. An additional simulation performed as a control was carried out without the monomer. All the simulations were done with a cutoff radius of 9.0 Å for the nonbonded interactions. The simulation length was 600 ps, with a time step of 1 fs, providing the first 200 ps for equilibration and the remaining 400 ps for data collection purposes. To describe the inter- and intramolecular interactions involving the monomers, the Dreiding force field25 was used. The electrostatic interactions mainly describe hydrogen bonds. The charges assigned to each atom are those that fit the electrostatic potential according to the results of the DFT calculations, by use of the MK18 procedure, which has proven useful for similar interactions; besides, no evidence was found that parameters for van der Waals energies need to be refitted for a specific electrostatic model.26 Such short-range intermolecular interactions are described by the 12-6 Lennard-Jones model, with parameters in the original Dreiding force field reference.25 For the cross-pair interactions the mixing rules of LorentzBerthelot20 were used. In the selection of atom type from the Dreiding force field, no explicit hydrogen atoms were chosen to participate in hydrogen bonds, that is, there is no specific force field defining such interaction; thus, we model the

Ind. Eng. Chem. Res., Vol. 46, No. 1, 2007 133

Figure 2. Optimized structures for L-proline: (left) in aqueous phase where the zwitterion is the stable form, with the hydrogen atom migrated from the hydroxyl group to the N atom; (right) in gas phase. The energy difference ∆E is defined as Eaqueous - Egasphase. The migrating hydrogen is inside the red ellipse. Color codes for the atoms are as given for Figure 1.

hydrogen bond as completely of electrostatic nature. The potential used to mimic water was the extended simple point charge model SPC/E,27 in which the molecule is assumed to be rigid and three point charges are considered, one for each of the atoms. This model performs very well in reproducing the structural and thermodynamic properties and the dynamics of water.28 3. Results and Discussion 3.1. DFT Calculations: 3.1.1. Structural Properties. The DFT-optimized structures of the inhibitor monomers are shown in Figure 1. An important feature common to the molecules studied is the presence of regions of hydrophobic and hydrophilic behavior, whose balance is believed to generate unusual structuring properties on water. For example, inhibitors such as VP7 have ring carbons that together with methyl and vinyl groups belong to the hydrophobic moiety, whereas the carbonyl group common to all the monomers belongs to the hydrophilic moiety; however, the role played by N of the amine group for 1FP, DMAEMA, and LP and that of the amide group for VP and Vval is not well defined. For this latter group, the stiffness of the C-N bond due to its partial double-bond character makes possible the existence of cis-trans isomerism.9 The preference of that kind of compounds for the trans conformation both in gas phase and in aqueous solution is confirmed by our calculations (Table 1). 1FP and DMAEMA do not present such isomerism because the symmetry of the ring in 1FP and that of the two methyl groups in DMAEMA does not allow the existence of a geometric isomer caused by torsion of the C-N bond, although the torsional barrier is found to be on the order of 25 kcal/mol, revealing the double-bond character of that bond. All the monomers studied are cyclic except the ester derivative DMAEMA, which appears to be in syn conformation with the first carbon of the dimethylaminoethyl group eclipsing the Cd O bond. The cis conformer (with the dimethylaminoethyl substituent group in cis position with respect to the oxygen of the carbonyl) is the only DMAEMA conformer for which optimized solutions were found in gas or aqueous phase, since attempts to obtain the trans structure ended up in the cis conformation because of the great steric effect when the first carbon of the dimethylaminoethyl group eclipses the methyl substituent of the acrylic group. The five-member rings of 1FP, VP, and LP are in envelope conformation, whereas in Vval, the six-member ring is in chairlike conformation, the carbons in the ring lie in staggered conformation, alleviating steric effects, and hydrogen atoms are arranged in alternating axial-equatorial directions. For VP and Vval, the N atom and the carbonyl and vinyl groups lie in the same plane, and a similar structure is found for the N, the R-carbon, and the carbonyl and hydroxyl groups in L-proline;

Table 1. Energy Differencesa in Gas and Aqueous Phase between Trans and Cis Conformers for VP and Vval VP ∆Etc

Vval

gas

sol

gas

sol

-3.52

-3.21

-4.08

-3.58

∆Etc )Etrans - Ecis in kilocalories/mole. The difference increases with ring size and is smaller in solution than in gas phase. a

in DMAEMA almost all of the heavy atoms lie in the same plane excepting the N and its two methyl substituents. On the other hand, 1FP is planar except for the carbon-3 that gives the envelope conformation to the ring. In fact all of the molecules present a plane containing a hydrophilic carbonyl group and a hydrophobic group out of the plane. The spatial arrangement of the atoms of the inhibitor determines the accessibility of the active sites and likely the orientation and other features of the inhibitor when it is adsorbed on the clathrate surface. The calculated atomic charges for the structures in Figure 1 are provided as Supporting Information. Except for LP, which as an amino acid shows the conversion to zwitterionic conformation, no drastic changes are observed between gas and water solution structures for a given molecule; therefore for the sake of simplicity, only the gas-phase structures are displayed in Figure 1. LP has two conformations (Figure 2), each one having a well-defined minimum in gas and aqueous phase, respectively: the normal conformation is the most stable for gas phase and the zwitterionic for the aqueous phase. In solution, the hydrogen migrates from the hydroxyl group to the N, leaving two oxygen atoms with high possibilities of interaction, possibly in a resonant structure where both oxygen atoms are indistinct from each other. This migration enhances the negative charge of the hydroxyl hydrogen and leaves the two hydrogen atoms attached to N, with a high positive charge. However, because the migrating proton holds a positive charge in the vicinity of the oxygen in cis position, a stronger solventmonomer interaction is found involving the oxygen in trans position with the nitrogen receiving the proton. For LP in aqueous phase, the formation of hydrogen bonds is likely to be found, as in the proline crystalline structure, where a three-dimensional hydrogen-bonding network linked via N-H‚ ‚‚O occurs.29 The same kind of interaction is likely to take place with the proton received by the nitrogen acting as a bridge between N and the oxygen atom of water molecules. Some general tendencies are found for all monomers. Table 2 summarizes the changes of the most important monomer geometric parameters in gas and aqueous phase. Slight changes are observed when the solvent effect on the monomers are incorporated by use of a polarizable continuum model (PCM); the angle O-C(O)-Y, where Y is O for DMAEMA and LP and N for the remaining monomers, increases for LP, 1FP, and DMAEMA, diminishing charge repulsion between negatively

134


Table 2. Summary of the Monomers’ Main Geometric Parameters monomer

phase

C(O)-O, Å

C(O)-Y.a Å

C(O)-CR, Å

N-Cx,b Å

N-Cz,c Å

O-C(O)-Y,a deg

DMAEMA

gas solution gas solution gas solution gas solution gas solution gas solution gas solution

1.21 1.22 1.20 1.26 1.22 1.24 1.21 1.23 1.22 1.23 1.21 1.23 1.22 1.23

1.35 1.34 1.34 1.25 1.35 1.34 1.39 1.37 1.39 1.38 1.39 1.37 1.39 1.38

1.50 1.50 1.54 1.55

1.47 1.46 1.49 1.51 1.47 1.48 1.48 1.47 1.49 1.48 1.46 1.47 1.47 1.48

1.46 1.46 1.49 1.52 1.47 1.48 1.40 1.40 1.41 1.41 1.39 1.39 1.40 1.40

123.17 123.21 123.35 128.08 124.90 125.16 125.82 125.53 122.22 121.90 125.71 124.95 121.71 121.55

LP 1FP VP cis Vval cis VP trans Vval trans

1.53 1.52 1.53 1.52 1.52 1.52 1.52 1.52

a Y represents N for VP, Vval, and 1FP and O for DMAEMA and LP. b x represents δ for LP and Vval, γ for VP, 1 for 1FPand the carbon of bond ethyl-N for DMAEMA. c z represents R for LP, 4 for 1FP, the methyl substituent for DMAEMA, and vinyl for Vval and VP.

charged atoms O and Y. All of the monomers show an elongation of the bond length between carbon and oxygen of the carbonyl group, owing to the interaction of the carbonyl group with the self-consistent reaction field that simulates an aqueous environment around the monomer. The carbonyl group is a site of high electronic density and the enhanced separation between the oxygen atoms and the molecule induces charge separation and therefore favors an increase of the dipole moment in aqueous solution. Similar elongation is found for the C-H distances and for the H-C-H angles, where the interaction with the solvent gives rise to an enhancement of the C-H distances and of the H-C-H angle due to repulsion between positively charged hydrogen atoms. On the other hand, the bond length between C of the carbonyl group and the next most electronegative atom attached to it is shorter (i.e., N in LP, VP, and Vval and O in LP and DMAEMA), revealing the presence of a stronger bond. If we take bond length as an indicator of bond strength, it can be inferred that the strong interaction with the solvent weakens the carbonyl bond, strengthening the bond between carbon and another atom. The carbonyl group, which has the strongest interaction with the surrounding water due to the strong polarization that enables strong hydrogen bonds, is responsible for the most important intermolecular interactions with the solvent; this is in agreement with former studies on VP interactions with a clathrate surface,5,8 where the carbonyl group was found to be the adsorption site on the clathrate crystal. Therefore, it seems likely that future clathrate inhibitors still will have the carbonyl group in their structure. Even though a high participation of the N atom in the intermolecular interactions could be expectedsbecause of its lone pairsits contribution to the inhibitor effect is less important, since the N position is not easily accessible or at least is clouded by the presence of other groups. For instance, in DMAEMA a “hydrophobic shield” is made up around the N position by the presence of the methyl groups; and in VP and Vval the performance of N as hydrogen-bond acceptor is eclipsed by the closeness of better hydrogen-bond acceptors such as oxygen interacting with the solvent. Yet DMAEMA is a more flexible molecule than the cyclic Vval and VP, and a rearrangement of the methyl groups will allow the interaction of water molecules with the N position, even when the interaction is not as strong as it would be if the access to the nitrogen atom were free. Anyway, DMAEMA has an advantage over VP and Vval because, even if small, the interaction with N exists, and that feature can be involved with the better performance of VC-713 as hydrate inhibitor. In this sense LP would also have that

Figure 3. ESP mapped on an electronic density surface isovalue of 0.0004 au; the color scale (blue, negative; red, positive) is the same for aqueous (left) and gas phase (right). A drastic change of the ESP for L-proline is due to conversion to the zwitterionic conformation. 1FP is poorly polarized in gas phase, but in aqueous phase the negative ESP on the carbonyl group is well defined. The ESP of VP and Vval (not shown) follow the same pattern of those of DMAEMA in both phases, with a negative ESP located on the carbonyl group slightly enhanced in aqueous with respect to gas phase.

advantage, but since LP is a rather rigid molecule, its inhibitor properties are not a priori clear, especially if compared with DMAEMA, and maybe the adsorption orientation may play against the interaction with water that has not crystallized yet. 3.1.2. Electronic Distribution. The electrostatic potentials (ESP) for 1FP, LP, and DMAEMA shown in Figure 3 illustrate drastic changes in some of the molecules (LP, 1FP) and negligible changes in others (DMAEMA) in going from gas phase to solution. As already mentioned, most potential clathrate polymer inhibitors share the characteristic of having monomers containing both the N atom and the carbonyl (CdO) group. Although the role of the carbonyl group is rather well understood, the role of N is not as well clarified. The monomers studied in this work have the carbonyl in different functional groups, leading

Ind. Eng. Chem. Res., Vol. 46, No. 1, 2007 135 Table 3. Dipole Moments µ for All Monomers in Gas and Aqueous Phasea LP

1FP Vval cis VP cis Vval trans VP trans DMAEMA

gas 5.98 4.76 sol 14.89 6.51 Dm 8.91 1.75

4.48 6.37 1.89

4.34 6.29 1.95

4.15 5.75 1.60

3.97 5.63 1.66

2.43 3.29 0.86

a ∆µ is defined as the difference µ sol - µgas. Units for dipoles are debyesD.

to different intramolecular interactions that could affect the carbonyl group; these different functions are aldehyde for 1FP, carboxylic acid for LP, ester for DMAEMA, and amide for VP and Vval. The application of the SCRF method affects the electronic distribution of each monomer; however, this is not clearly observed in the shape and location of the molecular orbitals, being almost unaffected on going from gas phase to aqueous phase (simulated by the SCRF); in general, the nitrogen and oxygen pz atomic orbitals make the largest contribution to the highest occupied molecular orbital (HOMO) and the carbon and hydrogen p and s atomic orbitals in the hydrophobic moieties make the largest contribution to the lowest unoccupited molecular orbital (LUMO). The absence of changes in the orbitals suggests no change in reactivity, but that is not our concern in this work. Yet changes in the values of dipole moment and in the electrostatic potential convoluting the molecule can be observed. In gas phase, a highly negative ESP surrounds the oxygen atom of the nucleophilic carbonyl group corresponding to its electrostatic negative charge, whereas on the aliphatic ring or on the aliphatic groups conforming the hydrophobic moiety, the ESP is just slightly positive and covers a much larger zone than the negative one (Figure 3). The difference between the ESP value on the N atom and the almost neutral value found on the hydrophobic moiety is small, supporting the idea that N does not participate in strong hydrogen bonds since these are basically electrostatic interactions; thus N appears to belong to the hydrophobic moiety rather than to the hydrophilic one. Nevertheless, for LP in gas phase a small zone of highly positive ESP below the N position is observed. In the aqueous phase, a larger zone of negative electrostatic potential can be observed for most monomers around the carbonyl group, even when there is no dramatic change in the structure of the molecule (Figure 3b,c). If the change in structure is severe, as in LP, the new ESP surface shows even more contrasting zones of positive and negative ESP (Figure 3a). That contrast indicates the displacement of electrons toward the red negative region, increasing the possibility of strong hydrogen bonding in that zone and probably favoring a good inhibitor performance since it is apparent that the interaction of the inhibitors and the methane-water-clathrate medium is of electrostatic nature, and charge distribution may be an influencing factor in adsorption phenomena due to favorable hydrogen bonding. However, at this point it is not clear yet whether, in a system like that already mentioned, the inhibitor would prefer to interact with the waters forming clathrate or with the waters still free in solution. For L-proline in gas phase, the oxygen of the hydroxyl group does not present as strong polarization as the carbonyl group, but once in solution where the zwitterionic conformation is stable, both O atoms are resonant and similar negative ESP is observed in Figure 3a. Although according to this resonance the strength of hydrogen bonding should be equivalent in both oxygen atoms, due to the lack of symmetry of the molecule there is a region of negative as well as of positive electrostatic potential in the neighborhood of one oxygen atom corresponding

Table 4. Changes of Energy (∆E) and Free Energy of Solvation (∆Gsol) for All of the Monomers, and Contributions to the Free Energy of Electrostatic and Nonelectrostatic Termsa monomer

∆E

∆Gelec

∆Gnonelec

∆Gsol

LPb 1FP cis-VP trans-VP cis -Vval trans-Vval DMAEMA

-19.18 -9.11 -8.86 -8.55 -8.11 -7.61 -8.15

-37.17 -9.49 -9.10 -8.80 -8.32 -7.82 -8.49

0.49 0.31 1.53 1.69 1.21 1.36 3.63

-17.27 -9.18 -7.57 -7.11 -7.11 -6.46 -4.86

a Energy units are kilocalories/mole. b For LP, ∆G sol is the sum of electrostatic and nonelectrostatic terms plus a correction to account for the change from zwitterion to gas-phase structure.

to the nearby N, such that a water molecule approaching this oxygen is affected by both potentials and the strength of a potential hydrogen bond is diminished. Nevertheless, L-proline is expected to have more affinity with the aqueous environment since it has double the number of favorable sites where strong interaction with water molecules may take place. 1FP evolves from a state possessing an electron density surface with a very similar ESP along the whole surface to one where localization of electrons on the oxygen of the aldehyde group is very clear, showing the solvent effect on this molecule. For DMAEMA the oxygen defining the ester function increases only slightly its negative polarity (Figure 3c). The ESP behavior for VP and Vval is similar to that of DMAEMA, having the negative potential on the carbonyl group but extending itself slightly toward the vinyl group. As seen in Figure 3, except for LP, molecules show a hydrophobic moiety characterized by an almost neutral potential on a large part of their surface, where little or no interaction with water is expected; thus the ratio of hydrophobic to hydrophilic moiety could be important in the inhibition process, with the hydrophilic one adsorbing on the crystal and the hydrophobic one occluding the pass of water toward the crystal. Because LP is much more polarized and has one more site to interact with water molecules, it is interesting to investigate the relative importance of adsorption on crystal with respect to the ability of a monomer for attracting free water molecules; this may be done by analyzing the behavior of LP in comparison to that of the other monomers. The tendency to increase polarization when interacting with the SCRF is common to all the monomers. The LP amino acid presents the highest polarization in aqueous solution and the linear molecule DMAEMA has the weakest polarization, as measured by the values of dipole moment in gas and aqueous phase (Table 3). LP shows a much greater change in polarization than the other molecules (8.91 D). Even though in gas phase it is still the most polarized molecule, its dipole moment is on the same order as those of the other monomers, but in aqueous phase the dipole moment of LP is more than double that of the closest dipole, which corresponds to 1FP. Clearly that is owed to the drastic change in structure of LP, since the other monomers do not suffer such dramatic change, instead bearing a ∆µ in the range of 0.9-1.8 D. It also can be observed that the same order occurring in gas phase for the dipole moment, LP > 1FP > cis-Vval> cis-VP > trans-Vval > trans-VP > DMAEMA, is kept in solution phase. 3.1.3. Energy and Stability. The structures of all monomers show a decrease in energy when passing from gaseous to solution phase, indicating a favorable interaction with the aqueous environment in the latter phase conferring more stability to the molecules. This result is favorable to the molecule

136


Figure 4. Energy differences for changes in conformation and environment for LP at HF level. ∆Es are calculated as EB - EA for the process A f B.

functioning as a clathrate inhibitor, assuming that the clathrate inhibitor acts by either adsorbing on the water-made crystal surface or having a structural effect on free water molecules. Table 4 summarizes changes in total and Gibbs free energies together with the main contribution to the free energies of solvation; the value of ∆E means Esol - Egas. LP showed the best solvation, in agreement with the fact that LP is highly hygroscopic and with high polarization. Although water affinity is a requirement, and a high free energy of solvation could be a signal of this, it appears not to be enough to determine good performance of a clathrate inhibitor. For instance, although VP has better solvation than DMAEMA, VC-713 has proven to be a better inhibitor than PVP. Nevertheless, this affirmation should be taken with caution since VC-713 also has VP in its structure and the versatility of three different groups could result in an enhancement of the inhibitor performance. Notice that the cis conformer has better solvation than the trans one for both VP and Vval. Although trans conformers are more stable both in gas phase and in solution, the cis one is better solvated because it has higher polarization and thus a stronger interaction with the reaction field; the dipole moment is a determinant factor to calculate the dominant electrostatic contribution to ∆Gsol. However, VP is better solvated despite having consistently lower polarization than Vval, with the corresponding conformer cis or trans in both phases; perhaps the difference is not so significant to counteract the fact that in Vval the hydrophobic moiety is larger. However the relationship is not clear. Also, a correspondence of ∆E with the free energy of solvation can be observed, at least for the cyclic molecules. However, DMAEMA appears not to follow the tendency indicating that the larger the value of ∆E, the better the solvation, apparently because of the much larger contribution of nonelectrostatic terms compared with the other monomers, due to a significant difference between the work done to build the PCM model cavity in comparison with the van der Waals terms. Notice that DMAEMA is the one that occupies a larger volume and has the larger ratio of hydrophobic to hydrophilic moiety. The PCM model calculates the free energy of solvation of a given molecule with the assumption that a molecule shows the

same stable geometry in gas phase aas in solution; for LP, that would be misleading since LP exists in gas phase in normal conformation and in aqueous phase in zwitterionic conformation. Because the zwitterion is not a well-defined minimum in gas phase and the normal gas-phase conformer is not the most stable in solution, converged structures for both zwitterion and normal conformation in both phases were performed at the lessexpensive Hartree-Fock (HF) level. A scheme with the changes in energy associated with changes of structure and phase for LP at HF level is illustrated in Figure 4. The ∆E values calculated for LP at the HF level are close to those at the B3LYP level, as seen by comparison of the ∆Es between normal LP gas phase and aqueous phase for B3LYP (Table 4) and HF (diagonal arrow in Figure 4). In order to calculate ∆Gsol for LP, the energy involved in the transition zwitterion-normal conformation in gas phase was subtracted from the value of ∆Gsol given by the PCM model, as an approximation to find a more realistic value of free energy of solvation (Table 4). Notice that according to the HF calculation the energetic barrier between the two conformations of LP diminishes considerably in solution. The nature of short-range interactions is studied by analyzing the geometries and energetic behavior of water-monomer structures and the stability of the interacting atomic pair. Full optimizations of these complexes were carried out at the B3LYP/ 6-311++g** level. Different initial configurations were tested in order to determine possible interactions of water with monomer sites, including the hydrophobic moiety. The local minimum structures displayed in Figure 5 show that the most important water-monomer interaction is via the carbonyl group, which is the one showing the largest binding energy. The distances involved in the interaction of water and the oxygen atom of the carbonyl group (Figure 5) are typical of those for hydrogen bonds. As a result of complexation, increases in the bond lengths O(w)-H(w) for water and C(O)-O of the carbonyl group are found. However, no correlation was found between the increase of these bond lengths or the hydrogen-bond length and the stability of the complexes. Moreover, in aqueous phase the PCM method predicts shorter distances than those found in gas phase; however, they show consistently lower binding energies. The stability of the complexes follows the order LP


Figure 5. Gas-phase structures corresponding to the highest binding energy complexes found for each monomer. Gas-phase interaction distances (indicated in the figures) are within the range of hydrogen bonds. Binding energy is calculated by ∆E ) Ecomplex - (Ewater + Emonomer), where each term is calculated in gas phase or solution for ∆Egas and for ∆Esol, respectively. The geometry of the complex does not change much in aqueous phase except for a slight diminution (0.5-1 Å) in the hydrogen-bond length. (a) DMAEMA; (b) cis-Vval; (c) trans-Vval; (d) 1FP; (e) LP; (f) cis-VP; (g) trans-VP.

> 1FP > cis-Vval > cis-VP > trans-Vval > trans-VP > DMAEMA, when calculated in vacuum (gas phase). For the cyclic molecules, the angle C(O)-O-H(w) is found to be within the range of 115-130° favoring the interaction between O(w) and the H-Cβ for LP and H-CR for the other cyclic monomers. This interaction certainly stabilizes the structure since other interactions involving just O‚‚‚H(w) are less energetically favorable. The interaction of water and the carbonyl of the noncyclic DMAEMA is linear with the C(O)-O bond and H(w)-O(w) in the same direction. The interaction distances of water with the hydrophobic moiety are much larger (values in the range from 2.5 to 2.8 Å) than those of hydrogen bonds, and the associated binding energy is very low (-1.1 to -1.8 kcal/mol difference). A favorable N-water interaction was found in LP, which is the only monomer with a nonneutral electrostatic potential around N and located far enough from the carbonyl group, but the corresponding gas-phase binding energy of -4.88 kcal/mol is lower than that obtained when H(w)-O is the dominant interaction. VP and Vval cis conformers were also found forming stable complexes with water (binding energies -3.5 and -4.0 kcal/ mol, respectively), but the position and orientation of water does not make clear whether the interaction is with the positively charged N or with the slightly negative vinyl group, possibly making a non conventional hydrogen bond with the π bond of the vinyl group as proton acceptor.30 This feature of VP and Vval could help to support the argument that pendent groups in clathrate inhibitors might play an important role in the inhibition, based on the observed interaction with lateral groups. This is another example of better interaction of water with the cis than with the trans conformer. Notice that although VP is better solvated than Vval according to free energy calculations with the PCM model, Vval shows better binding energies for these water-monomer complexes than the corresponding conformers in VP.

The gas-phase lowest energy complexes were also optimized in solution, using the PCM approach to build the cavity and the SCRF approach to simulate the solvent, thus assessing the relative role of long-range interactions. A diminution in the binding energy is observed for the complexes of the cyclic monomers, especially in trans-Vval, which has a binding energy close to zero. The optimization for DMAEMA-water had convergence problems; therefore we report the energetic results of a single-point calculation on the gas-phase complex geometry, which shows an increase of the binding energy in the presence of the self-consistent reaction field. This result is at least qualitatively correct when it is considered that a geometry optimization would reduce the energy of the complex even further. The interaction distances in solution for all cases are found to be shorter than in gas phase, but the increase in the O(w)-H(w) bond length of the proton donor group in water and that involving the acceptor O group in the monomers is found to be of the same order as in gas phase. Complexes of two water molecules and a monomer were also optimized to investigate the effect of the increase of the number of water molecules on the interaction properties. Although the carbonyl group was expected to hold hydrogen bonds with the two waters, in all cases the most stable structure obtained was that of water involved in hydrogen bonding both with the other water and with the carbonyl group (Figure 6). Yet it is not clear if this kind of behavior is useful or not for a prospective clathrate inhibitor, since it is necessary to establish if this predicts the molecules to be “distracted” from forming the clathrate crystal or if it is a harbinger of strong adsorption on the clathrate surface or if it promotes the crystal growing. In any case, it could be taken as an indication of the energy of adsorption of the inhibitor on the crystal, based on the fact that the adsorption on the crystal is given by the interaction between one pendent hydrogen in the clathrate molecules and a carbonyl group of the inhibitor. The structural effect shown in Figure 6a

138


Figure 6. (H2O)2-monomer complexes: (a) one water forming hydrogen bonds with the carbonyl group and with the other water; (b) two water molecules forming hydrogen bonds with the carbonyl group. In all cases structure a is more stable than b.

could be also an indication of the effect of some clathrate inhibitors such as PVP and VC-713, which are able to change the type of crystallization structure31 of the hydrate, perhaps by distorting the arrangement tendency or growing habit of the water molecules.31 Besides, the favorability toward the structure type a could be interpreted that even though the monomers interact strongly with water and eventually the monomer could attach some waters to itself, the preferred interaction of water molecules is with themselves and therefore favoring the formation of crystal, which is expected because the effect of the inhibitor is kinetic and not thermodynamic. Thus, the main role of the carbonyl group should be the adsorption on crystal surface rather than its interaction with free water. For Vval we could not obtain a complex showing two hydrogen bonds with the carbonyl group (as in b), which suggests that, if adsorbed on the clathrate surface, Vval might not have a hydrogen bond with free water, or vice versa, if adsorbed to free water, it might not adsorb to the clathrate surface. An additional effect of the inhibitor in the polymer form may be given by the interaction of some nonadsorbed carbonyl groups with free water molecules delaying their diffusion toward the crystal surface. 3.2. Classical Molecular Dynamics: 3.2.1. Structure. In the computational details we mentioned that although different temperatures from 298 to 268 K were used to run simulations, no significant differences in the radial distribution functions (RDFs) were found except for the sharpness of the peaks, which increases at lower temperatures (clearly due to the decrease in mobility of water), but the peaks’ location and the number of water molecules in each shell of solvation remained unaltered, so we considered that the study of the arrangement of water at room temperature is adequate to get an insight into the behavior of the system water-inhibitor even at low temperatures where clathrate is formed. Hence this section is based upon simulations at 298 K. The RDFs revealing interactions between the carbonyl O and water are displayed in Figure 7. All the monomers show a first peak for the pair O‚‚‚O(w) at a distance of 2.6 Å and a first peak for the pair O‚‚‚H(w) at 1.6 Å; the latter accounts for hydrogen-bonding interactions with a distance slightly smaller than that given by DFT calculations (1.8 Å). The number of water molecules in the first coordination shell for DMAEMA is 2.3, for Vval 2.6, for VP 2.8, and for 1FP and LP 3 and 3.4, respectively; this number is commensurate with the monomer solubility as given by the solvation free energy. The increase in ring size seems to diminish the effectiveness of the carbonyl group interacting with water; according to the coordination numbers found around the O site in VP and Vval, DMAEMA is the least effective monomer for interacting with a given number of waters via the carbonyl group. This is in agreement with the fact that DMAEMA is the only molecule among the studied ones where the carbonyl group is less accessible due to steric effects arising from the relative

proximity (2.5 Å) of the hydrophobic end carbon of the acrylic group and the initial carbon of the ethyl group. Though cis-VP and cis-Vval have the end of the vinyl group at a similar distance, for them this carbon shows an electronegative potential that rather helps to increase slightly the carbonyl capability of interaction, at least in relation to DMAEMA. A second wide peak in the RDFs is encountered at 4.2 Å, delimited for two valleys at 3.2 and 5.4 Å; the depth of the valleys and the definition of the peaks is again correlated to the solubility of each monomer being stronger for LP. The presence of the valleys shows the existence of a second shell of solvation where there are 13.2 waters for DMAEMA and Vval, 13.4 for VP, 15 for 1FP, and 15.1 for LP. A third shell of solvation can be observed with 38.7 waters for 1FP, 38.9 for LP, 42.4 for DMAEMA, 54.7 for VP, and 55.4 for Vval; that is to say, the third solvation shell for VP and Vval is broader than the corresponding ones in the other monomers, and in 1FP the appearance of the shell is not clearly defined. Interactions with the ester oxygen of DMAEMA are negligible on the basis of its RDF, which does not present prominent peaks; moreover, in Figure 7 a low number of water molecules is involved in interaction with this group in comparison to that corresponding to the carbonyl groups in each monomer. On the other hand, interactions with the hydroxyl oxygen of LP are of the same nature as those with the carbonyl one, differing only in intensity as evidenced by its RDF. This shows that although interaction with the LP hydroxyl oxygen is less strong than that with the carbonyl one, it is strong enough to provide an advantage for LP with respect to the other molecules concerning the number of water molecules strongly interacting with it. The RDF corresponding to the interactions with the N atom are shown in Figure 8. Weak interactions are observed in DMAEMA, where a well-defined peak with location characteristic of a hydrogen bond but with a low intensity is found; the low coordination value indicates that occasionally no water is found within this solvation shell. That is perhaps due to the difficult accessibility to the N atom, where a water molecule has to overcome steric effects of the methyl groups to reach the nitrogen position. For DMAEMA the peak for N-H(w) is located at shorter distances than for N-O(w), indicating the interaction with the water bond O(w)-H(w) as proton donor, whereas for LP the peak for N-O(w) is the one located at shorter distances; thus the water molecule is the proton acceptor, although not exactly to the N but rather to the H previously attached when the zwitterion was formed. It is because of this fact that the location of this first peak for N-O(w) in LP is located farther than N-H(w) even though the interaction through the hydrogen bond is stronger in the former (shell coordination number 4.6). In the second shell, the difference still holds with 27 and 32 waters for LP and DMAEMA, respectively. Figure 8 also illustrates that for 1FP, VP, and Vval the interactions for the pair N-O(w) or N-H(w) are not as clear as for LP and DMAEMA, and the shapes of the RDFs are more similar to the RDFs around atoms belonging to a hydrophobic moiety as C-3. Also it can be observed from the plot of the coordination number versus shell radius how the interaction increases according to the order Vval < VP < 1FP. The interaction with the hydrophobic moiety is also observed from the RDF peak appearing in the interval 4-5 Å, revealing 1012 waters solvating around the atom C-3. The carbon atom type C-3 was chosen to describe this interaction since these atoms are located in the rings or in the hydrophobic groups, whereas the C-2 type was not chosen because in most of the monomers studied this RDF would just reflect the solvation around oxygen


Figure 7. Radial distribution functions and water coordination numbers versus distance (r, given in angstroms) from O: black, DMAEMA; red, LP; blue, 1FP; green, VP; orange, Vval. The interaction O-O(w) is shown in red, and O-H(w), in black. The intensity of the peaks and the number of waters solvated around the O carbonyl at a given distance increases as ∆Gsol does according to the ab initio calculations (Table 4). The number of water molecules around the ester oxygen in DMAEMA and hydroxyl oxygen in LP are plotted as dashed lines.

since in most monomers it is too close to the oxygen of the carbonyl group, with the shape of the pair correlation being affected by the strong solvation around oxygen. RDFs and solvation numbers around C-3 are shown in Figure 9. Slight broadening in the RDF for the pair C-3‚‚‚O(w) is observed for the cyclic molecules at 3.8 and 4.6 Å; the intensity of the broadening decreases in the order LP > 1FP> VP> Vval as the hydrophobic character of the simulated inhibitor decreases. For DMAEMA, it is possible to observe from the C-3‚ ‚‚O(w) RDF a group of water molecules around the methyl groups at a distance of 4 Å, which can be considered a welldefined solvation shell separated from bulk water whose peaks appear at 6 Å from the methyl groups. When the solvation of water around the hydrophobic moiety is compared with that around methane in solution, it is found that the solvation around the methyl groups in DMAEMA

resembles more the RDF of water and methane (cage effect)32 than the solvation around any of the rings does; this could be taken as an indicator of hydrophobic effect around these groups being similar to the one occurring around methane, and according to that it would be easier to form clathrate-like structures around DMAEMA than around the other inhibitors. This could be an advantage for hydrophobic solvation as shown by the plot of the integral of the RDF around C-3 (Figure 9), where it is observed that the solvation around the hydrophobic moiety of DMAEMA is slightly more effective than for the other monomers. This is probably a factor that makes VC-713 more effective in the inhibition than PVP,31,33 since the arrangement of water resembling clathrate around the methyl groups of DMAEMA, where anyway the cage cannot be completed because of the size of the inhibitor, can be a delaying factor in the successful formation of clathrates. Moreover, this kind of

140


Figure 8. Radial distribution functions of O(w) and H(w) around N in each monomer and the number of water molecules versus distance from N: black, DMAEMA; red, LP; blue, 1FP; green, VP; orange, Vval.

Figure 9. Radial distribution function for the pairs C_3sO(w) and C_3sO(w). At the right, solvation number versus distance from C_3. Black, DMAEMA; red, LP; blue, 1FP; green, VP; orange, Vval. Distance (r) is given in angstroms.

effect could be still in effect for the other inhibitors, as reflected by the RDFs showing a hydrogen-bond network between waters around the hydrophobic moiety. 3.2.2. Dynamics. The water molecule has the largest part of its mass concentrated in the oxygen atom; thus the analysis of its dynamics can be done by means of aanalysis of the motion of this atom. In Figure 10 it is apparent that there is a slight change in the dynamics of water as reflected for the time evolution of the mean square displacement (MSD). The diffusion coefficients were calculated by linear regression in the interval of 8-20 ps, and the values are shown in Table 5. It can be seen that the diffusion coefficient for the water-inhibitor system is larger than the diffusion coefficient for pure water except for the case of water-1FP. This result is remarkable since, owing to the strong interactions of the monomers with free water molecules, a slight diminution in the global average motion of water would be expected. This apparent contradiction is understood if these interactions are thought of as having a dynamic nature instead

of a static one. For instance, there is a competition of the residence times of water molecules forming H-bonds with the inhibitor versus those for water-water H-bonds; this competition accounts for both the structural and dynamic results. The power spectra for the water molecules are shown in Figure 11. For SPC/E water, the peak at 40 cm-1 corresponds to the bending interaction O(w)‚‚‚O(w)‚‚‚O(w) and the shoulder at 240 cm-1 corresponds to O(w)‚‚‚O(w) bond stretching. The modification on the power spectra for O(w) can be a qualitative estimative of how much the tetrahedral structure of water is modified by addition of the inhibitor. The curves show a fitting to a sixth-grade polynomial to smooth the curves, which reveals the very small differences among the systems. The pure water spectrum shows a slightly higher peak for the vibrations O(w)‚ ‚‚O(w)‚‚‚O(w) than when the inhibitor is present; the location of the peak is not modified by the introduction of inhibitor to the system. Although this reduction could be not very significant as to indicate that is a determinant factor in the inhibition, it can be noticed that VP modifies this characteristic of water to


Figure 10. Mean square displacement for the water-inhibitor system and for pure water. The figure at the left shows the msd for the first 0.6 ps, and that at the right, for the last 12 ps. Black, DMAEMA;-red, LP; blue, 1FP; green, VP; orange, Vval; cyan, pure water.

Figure 11. (Left) Power spectra for the different water-inhibitor systems and pure water; (right) velocity autocorrelation functions from which the power spectra are calculated. Black, DMAEMA; -red, LP; blue, 1FP; green, VP; orange, Vval; cyan, pure water. Table 5. Diffusion Coefficientsa Calculated for Water-Inhibitor Systems from MSD Plots with Two Different Intervals system

D × 109 (8-20 ps)

D × 109 (0.8-20 ps)

VP-water Vval-water DMAEMA-water LP-water 1FP-water pure water

2.29 2.24 2.10 2.07 1.96 1.98

2.29 2.24 2.14 2.09 2.02 1.97

a

D values are given in square meters per second.

a lesser extent than the other inhibitors; Vval and 1FP alter the behavior of water in a similar way but not as much as LP and DMAEMA do. As with the MSD, the variation is very small since the interaction of the monomer with the water solution seems to be very local, not transcending to long range, and the average dynamics of water molecule is not greatly affected; notice, however, that the only two inhibitors, DMAEMA and LP, that present alternate sites for hydrogen bonding other than the carbonyl group are the ones mostly affecting the structure of water. VP, Vval, and 1FP, whose only strong active site is the carbonyl according to DFT calculations and the RDFs, affect

water behavior to an extent that can be considered correlated to the binding energy for the hydrogen bond between the carbonyl group and water. 4. Summary and Conclusions DFT and classical MD simulations suggest that the group determining the interaction with water, and in consequence with a clathrate, is the carbonyl group, LP was found to be the one with the strongest interactions via hydrogen bond. Short-range interactions are dominant, and evidence of longrange structures of water was not found in spite of the presence of more than one solvation shell found around the carbonyl group. Since experimental data indicate that VC-713 is a better inhibitor than PVP, we conclude that an analysis of the free energy of solvation is not enough to predict the good performance of an inhibitor. The nitrogen site was found to belong to the hydrophilic moiety only for LP and DMAEMA, giving an additional advantage for these two molecules; this advantage is even more accented in DMAEMA because of its flexibility, which allows higher versatility of interactions.

142


The better binding energy for water-water than for waterinhibitor confirms that there is no thermodynamic effect for the kind of inhibitors in this study. The increase in the size ratio of hydrophobic moiety/hydrophilic moiety seems to weaken the strength of hydrogen bonding of the carbonyl group with waters. Vval, which was found to be able to engage in just one hydrogen bond, is thought to lose its ability to interact with water once adsorbed on the crystal. On the contrary, LP has the possibility of using one of its oxygen atoms to adsorb on the crystal surface and the other to interact with free water. Although some alterations of the dynamics of water were shown, those appear to be just a consequence of short-range interactions, not altering much the global behavior of water. Methyl groups in DMAEMA behave as methane in the presence of water, resembling the formation of clathrate around, but as the cavity cannot be completed because of the volume of the DMAEMA molecule where the methyl is attached, this could be a factor in retarding the formation of crystal. This analysis suggests that a good inhibitor must have a mix of different features given by the various monomers in this study; therefore it would be interesting, for example, to combine the strength of adsorption of LP with the flexibility and alternative kind of interaction offered by DMAEMA.

N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, Revision C.02; Gaussian, Inc.: Wallingford, CT, 2004. (14) Becke, A. D. Density-functional Exchange-energy approximation with asymptotic behavior. Phys. ReV. A 1988, 38, 3098. (15) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys. ReV. B 1988, 37, 785-789. (16) Becke, A. D. Density-functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98, 5648-5652. (17) Singh, U. C.; Kollman, P. A. An approach to computing electrostatic charges for molecules. J. Comput. Chem. 1984, 5, 129-145. (18) Besler, B. H.; Merz, K. M.; Kollman, P. A. Atomic charges derived from semiempirical methods. J. Comput. Chem. 1990, 11, 431-439. (19) Boys, S. F.; Bernardi, F. Mol. Phys. 1970, 19, 553.

Acknowledgment This work is supported by the Gas Processors Suppliers Association. Supporting Information Available: Calculated atomic charges for the structures in Figure 1. This material is available free of charge via the Internet at http://pubs.acs.org. Literature Cited (1) Sloan, E. D. Clathrate Hydrates of Natural Gases; Marcel Dekker: New York, 1998. (2) Behar, E.; Delion, A. S.; Sugier, A.; Thomas, M. Plugging control of production facilities by hydrates. Ann. N.Y. Acad. Sci. 1994, 715, 94105. (3) Anderson, F. E.; Prausnitz, J. M. Inhibition of gas hydrates by methanol. AIChE J. 1986, 32, 1321-1333. (4) Makogon, Y. F.; Makogon, T. Y.; Holditch, S. A. Kinetics and mechanisms of gas hydrate formation and dissociation with inhibitors. Ann. N.Y. Acad. Sci. 2000, 912, 777-796. (5) Carver, T. J.; Drew, M. G. B.; Rodger, P. M. Inhibition of Crystal Growth in Methane Hydrate. J. Chem. Soc., Faraday Trans. 1995, 91, 3449-3460. (6) Sloan, E. D. U.S. Patent No. 5,432,292, Colorado School of Mines, 1995; Vol. 432 (292). (7) Carver, T. J.; Drew, M. G. B.; Rodger, P. M. Molecular dynamics calculations of N-methylpyrrolidone in liquid water. Phys. Chem. Chem. Phys. 1999, 1, 1807-1816. (8) Kvamme, B.; Huseby, G.; Kristian, O. Molecular dynamics simulations of PVP kinetic inhibitor in liquid water and hydrate/liquid water systems. Mol. Phys. 1997, 90, 979-991. (9) Aparicio-Martinez, S.; Hall, K. R.; Balbuena, P. B. Theoretical Study on the Properties of Linear and Cyclic Amides in Gas Phase and Water Solution. J. Phys. Chem. A 2006, 110, 9183-9193. (10) Wu, Y.; Banoub, J.; Goddard, S. V.; Kao, M. H.; Fletcher, G. L. Antifreeze glycoproteins: relationship between molecular weight, thermal hysteresis and the inhibition of leakage from liposomes during thermotropic phase transition. Comp. Biochem. Physiol. Part B: Biochem. Mol. Biol. 2001, 128, 265-273. (11) Jussaume, L.; Canselier, J.-P.; Monfort, J.-P.; Rodger, P. M. Molecular Modelling of Additives on type II gas hydrates surfaces: Development of a criterion efficiency. AIChE Top. Conf. 1998. (12) Miertus, S.; Scrocco, E.; Tomasi, J. Chem. Phys. 1981, 55, 117. (13) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, Jr., J. A.; Vreven, T.; Kudin, K.

(20) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: Oxford, U.K., 1990. (21) Smith, W.; Forester, T. R. DL_POLY; Daresbury Laboratory: Daresbury, U.K., 1996. (22) Hoover, W. Canonical dynamics: equilibrium phase-space distributions. Phys. ReV. A 1985, 31, 1695. (23) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. Numerical integration of the cartesian equations of motion of a system with constraints: Molecular dynamics of n-alkanes. J. Comput. Phys. 1977, 23, 327-341. (24) Essmann, U.; Perera, L.; Berkowitz, M. L. A smooth particle mesh Ewald method. J. Chem. Phys. 1995, 103, 8577-8593. (25) Mayo, S. L.; Olafson, B. D.; Goddard, W. A. Dreiding: A generic force field for molecular simulations. J. Phys. Chem. 1990, 94, 88978909. (26) Leusen, F. A study of different approaches to the electrostatic interaction. organic crystals. Phys. Chem. Chem. Phys. 2003, 5, 49234931. (27) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The missing term in effective pair potentials. J. Phys. Chem. 1987, 91, 6269-6271. (28) Guillot, B. A reappraisal of what we have learnt during three decades of computer simulations in water. J. Mol. Liq. 2002, 101, 219260. (29) Janczak, J.; Luger, P. L-Proline monohydrate at 100 K. Acta Cryst. Sect. C: Crystal Struct. Commun. 1997, 53, 1954-1956. (30) Chowdhury, H.; Rahaman, S. H.; Ghosh, R.; Sarkar, S. K.; Fun, H.-K.; Ghosh, B. K. Synthesis and characterization of copper(II)-bis(pyridin2-carbonyl)amine (pyca)-pseudohalide complexes: 2D superstructure formation through CH..O hydrogen bond and [pi]..[pi] interaction in [Cu(pyca)Cl(NCO)]. J. Mol. Struct. 2006 (in press). (31) Koh, C. A.; Westacott, R. E.; Zhang, W.; Hirachand, K.; Creek, J. L.; Soper, A. K. Mechanisms of gas hydrate formation and inhibition. Fluid Phase Equilib. 2002, 194-197, 143-151. (32) Ben-Naim, A. Hydrophobic Interactions; Plenum: New York, 1980. (33) Lederhos, J. P.; Long, J. P.; Sum, A.; Christiansen, R. L.; Sloan, J. E. D. Effective kinetic inhibitors for natural gas hydrates. Chem. Eng. Sci. 1996, 51, 1221-1229.

ReceiVed for reView August 24, 2006 ReVised manuscript receiVed October 17, 2006 Accepted October 23, 2006 IE061123G

Computational Studies of Structure and Dynamics of Clathrate

Recommend Documents