MM Simulations: Length Dependence of

J. Phys. Chem. B 2008, 112, 7061–7070

7061

Fragmentation-Based QM/MM Simulations: Length Dependence of Chain Dynamics and Hydrogen Bonding of Polyethylene Oxide and Polyethylene in Aqueous Solutions Hui Li, Wei Li, Shuhua Li, and Jing Ma* School of Chemistry and Chemical Engineering, Institute of Theoretical and Computational Chemistry, Key Laboratory of Mesoscopic Chemistry of MOE, Nanjing UniVersity, Nanjing 210093, P. R. China ReceiVed: January 24, 2008; ReVised Manuscript ReceiVed: March 14, 2008

Molecular fragmentation quantum mechanics (QM) calculations have been combined with molecular mechanics (MM) to construct the fragmentation QM/MM method for simulations of dilute solutions of macromolecules. We adopt the electrostatics embedding QM/MM model, where the low-cost generalized energy-based fragmentation calculations are employed for the QM part. Conformation energy calculations, geometry optimizations, and Born-Oppenheimer molecular dynamics simulations of poly(ethylene oxide), PEOn (n ) 6-20), and polyethylene, PEn (n ) 9-30), in aqueous solution have been performed within the framework of both fragmentation and conventional QM/MM methods. The intermolecular hydrogen bonding and chain configurations obtained from the fragmentation QM/MM simulations are consistent with the conventional QM/MM method. The length dependence of chain conformations and dynamics of PEO and PE oligomers in aqueous solutions is also investigated through the fragmentation QM/MM molecular dynamics simulations. 1. Introduction With recent advances in computer science, theoretical methodologies, and computational technologies, various levels of theoretical methods have been widely used in studying macromolecules at different levels, from quantum mechanics (QM) calculations on electronic structures and atomistic molecular mechanics (MM) simulations of chain and packing structures to coarse-grain modeling of mesoscopic systems. Although the MM-based classical molecular dynamics (MD) and Monte Carlo simulations are frequently applied to macromolecule solutions, the limits of force fields (such as the poor transferability of the force-field parameters and the many-body interactions) are distinct. Thus, it is necessary to study macromolecule solutions at higher level. For simulations of solutions, there are mainly two types of solvent models, the implicit continuum dielectric approach and the explicit solvent model.1–4 In the computationally economical continuum models, some special short-range interaction, such as hydrogen bonds, cannot be described well. On the other hand, it is impossible to use the explicit model (with thousands of solvent molecules) to include long-range interactions in full QM calculation because of tremendous computational cost. Thus, several kinds of hybrid models5 have been developed for describing both specific and bulk effects. Another efficient way to describe the solvent effect is the hybrid QM/MM method, which was originally proposed by Warshel and Levitt.6 In the QM/MM framework, the solute is treated as the QM part, and the solvent molecules are taken as the MM part. Many works on QM/MM methods have been published in recent years.7 However, the size of the QM part in conventional QM/MM is still restrained by the exponential scale of computational cost with the system size in conventional QM methods. Thus, for the macromolecule solutions, if we want to describe the largesized solute in terms of QM calculations, an efficient low-cost approach for quantum chemical calculations is required. * Corresponding author. E-mail: [email protected].

Molecular fragmentation methods for electronic structure calculations have been developed for years,8–20 and with such methods, it is possible to perform ab initio calculations for large systems containing hundreds and even thousands of atoms. The success of fragmentation-based methods lies in the locality and transferability of properties of a subsystem in molecules. Belonging to this category, the energy-based approaches have advantage of easy implementation at various theoretical levels for geometry optimization and calculations of various properties of large-sized systems.15–18 The basic idea was independently proposed by Li et al.15 and Collins’s group.16 Recently, Gadre et al. also formulated their previous molecular tailoring approach17b,c into a similar energy-based strategy.17a In order to treat the polar and charged macromolecules, the generalized energy-based fragmentation (GEBF) approach15d employs background point charges to consider the long-range interaction in the calculation of each fragment. This approach can reproduce the ground-state energies of highly polar and charged systems within a few millihartrees (0.63 kcal/mol). In addition, some electronic properties, such as atomic charges, dipole moments, and static polarizability, can also be precisely predicted.15d On the other hand, within the framework of the density matrix fragmentation scheme, Gadre et al.17c and Mezey et al.10e,f have given satisfactory results of electrostatics potential, electron density, and dipole moments, and recently, Mezey et al.10f have further improved the accuracy of these one-electron properties by employing point charges to mimic the distant parts of macromolecules. Such a strategy was also applied in some other works by Sakai et al.,19 Jiang et al.,15f and Truhlar’s group.18 In fact, the introduction of background point charges in the fragmentation quantum chemical method to consider the longrange electrostatics interaction is the same idea as that behind electrostatics embedding model21 in QM/MM method. In this way, the polarization effect from the MM part can be partially considered. In the present work, we implement a hybrid method by combining GEBF calculations with MM to study dilute solutions of macromolecules. The combined strategy of fragmentation

10.1021/jp800777e CCC: $40.75  2008 American Chemical Society Published on Web 05/20/2008

7062 J. Phys. Chem. B, Vol. 112, No. 23, 2008

Li et al. of the MM part treated by force field. For the studied aqueous solutions of oligomer (in which no covalent bond exists between the QM and MM parts), EQM-MM is the interaction between the elec QM and MM parts, consisting of electrostatics (EQM-MM ) and VdW van der Waals (EQM-MM) interactions. Equation 1 can be written as Frag elec VdW ETOL ) EQM + EMM + EQM-MM + EQM-MM

method,15d

Within the framework of the GEBF energy of the system can be expressed as Figure 1. Schematic description of fragmentation QM/MM models, where the QM (oligomer) and MM (solvent molecules) parts are demonstrated by Corey-Pauling-Koltun and stick models, respectively. The red region in the center of the QM part is a fragment in molecular fragmentation, and the green regions are considered together with the central fragment in a subsystem for QM treatment. The yellow regions, which are far from the fragment, are taken as background point charges in fragmentation QM calculations.

quantum chemistry and MM, called fragmentation QM/MM, in the electrostatics embedding model is applied to calculations of conformation energies, geometry optimizations, and BornOppenheimer molecular dynamics (BOMD). Poly(ethylene oxide), PEO, is a kind of water-soluble polymer with extensive industrial applications, including cosmetics, pharmaceutical formulations, adhesives, lubricants, films, and so on.22 Polyethylene (PE) is the most important polymer in the plastic industry.23 Several MD and Monte Carlo simulations within the framework of MM have been carried out to investigate the configuration and other physical characteristics of PEO24 and PE25 in various media. Therefore, dilute PEO and PE solutions are fundamental systems to test the performance of our fragmentation QM/MM method. Chain conformations and hydrogen-bonding interactions between solutes and solvents are studied in comparison with the conventional QM/MM method. Furthermore, oligomers with different chain lengths are analyzed on the basis of our fragmentation QM/MM MD simulations. The article is organized as follows. In Section 2, we will give a brief introduction of the fragmentation QM/MM Hamiltonian, as well as the details of geometry optimizations and MD simulations. In Section 3, we will compare the fragmentation QM/MM results of conformation energies, geometry optimizations, and BOMD simulations for PEO and PE with those of conventional QM/MM method. Finally, the length dependence of configurations and chain dynamics of PEO and PE have been investigated. 2. Fragmentation-Based QM/MM Model 2.1. EnergiesandForces.Asintheelectrostaticsembedding7o,21 model, the partial point charges of the MM part are included in the Hamiltonian of the QM part when the GEBF calculation15d is performed. The partition of a long oligomer into fragments in the QM part is illustrated in Figure 1. The subsystem, constructed by the fragment (red region) and its nearby fragments (green regions), is calculated by QM, and the distant QM fragments (yellow regions) and MM parts are considered as background point charges in the QM calculation of each subsystem. The GEBF QM/MM Hamiltonian can be expressed as Frag ETOL ) EQM + EMM + EQM-MM Frag EQM

(1)

where is the energy of the QM part determined by using the energy-based fragmentation method and EMM is the energy

Nsub Frag EQM )

∑

(∑ ) ∑

Natom

Nsub

CiE˜i -

i)1

Ci - 1

a>b

i

(2)

the total QM

qaqb rab

(3)

where Nsub is the number of subsystem, and E˜i and Ci are the energy (including the self-energy of the point charges) and coefficient of the ith subsystem, respectively. The details were given in ref 15d. In combination with the MM part in the GEBF QM/MM method, we can obtain E˜i′, the total energy of the ith subsystem embedded in background point charges on those distant atoms in the QM part and solvent molecules in the MM part. Finally, the total energy with the fragmentation QM/MM is obtained.

(

Nsub

ETOL )

Nsub

)

QM+MM

∑ CiE˜i′ - ∑ Ci - 1 ∑ i)1

i)1

a>b

qaqb MM qaqb + rab rab a>b

∑

VdW EMM + EQM-MM (4)

where in the second term, the electrostatics interaction between atomic charges covers the whole system including both QM and MM parts. The electrostatics interaction within the MM part is also counted in EMM and should be removed, as shown in the third term. It should be mentioned that if the QM and MM parts are linked by some covalent bonds, the fragmentation QM/MM scheme will become complicated, deserving our efforts in a future work. The forces in the fragmentation QM/MM model are divided into two groups: forces involved in the QM part and forces involved in the MM part. Forces on the QM part can be derived as

f QM R ) gR -

∂ VdW E ∂rR QM-MM

(5)

where gR is the force on the Rth atom, which is obtained from the fragmentation QM calculation, and rR is the position of the Rth atom in the QM part. By using QM charges, which are taken from natural population analysis,26 the force on the MM part is approximately written as

f MM ≈a

∂ VdW (Eelec + EMM + EQM-MM ) ∂ra QM-MM

(6)

where Eelec QM-MM is the electrostatics interaction between the QM and MM parts, calculated by using point charges, and ra is the position of the ath atom within the MM part. 2.2. Geometry Optimization. For a dilute solution of a long oligomer with thousands of atoms, it is not practical to use conventional quasi-Newton algorithms of optimization, which need at least Nf (Nf + 1)/2 (where Nf is the number of freedoms) storage of the Hessians. Instead, the limited-memory BroydenFlether-Goldfarb-Shanno method,27,28 which utilizes the Hessian information accumulated from the gradients over the history of optimization and scales as O(Nf), is adopted in geometry optimization. Because of the approximations made in the force

Fragmentation-Based QM/MM Simulations

J. Phys. Chem. B, Vol. 112, No. 23, 2008 7063

TABLE 1: Studied Models of Aqueous PEOn and PEn Solutions for Fragmentation QM/MM Simulations oligomer

number of solvent molecules

volume (Å)

density (g/cm3)

PEO6 PEO8 PEO10 PEO15 PEO20 PE9 PE12 PE15 PE21 PE30

1339 1555 1912 1847 3595 1249 1885 2346 3285 5175

37.6 × 37.9 × 31.4 39.9 × 43.2 × 31.8 45.4 × 45.8 × 31.6 31.7 × 64.5 × 31.3 57.6 × 66.7 × 31.4 51.2 × 30.2 × 29.8 44.4 × 50.7 × 27.7 48.9 × 57.6 × 26.7 57.8 × 69.8 × 29.8 70.5 × 87.5 × 27.6

0.90 0.86 0.88 0.88 0.90 0.82 0.91 0.91 0.82 0.92

calculation, the forces do not correspond to exact gradients of the potential energy surface. This is not a serious problem when the gradients are large, but when the system gets close to the minimum, where the forces are very small, the error will make the fragmentation QM/MM forces hard to converge. Furthermore, for such large systems containing thousands of atoms and innumerable minima, it is meaningless to locate an exact minimum; therefore, we give the criterion of convergence of maximum force less than 0.002 Eh/Å, and the energy change is smaller than 10-4 Eh (0.063 kcal/mol). Of course, the accuracy of the gradient increases with the increase in the size of fragment, and more strict criteria could be used. 2.3. BOMD Simulations. In MD simulations, the approximations made in calculations of forces give rise to slight errors of pairwise forces of the QM/MM method. Although these errors do not affect geometry optimization significantly, they can be accumulated and cause problems in Hamiltonian energy conservation in MD simulations (Supporting Information, Figure S1), as also noticed in ref 29. There have been some attempts to develop efficient methods to improve the energy conservation in QM/MM MD simulations, such as smooth functions for modifying the forces and energies.30–32 In our simulations, we use the Nose-Hoover-chain (NHC)33,34 method for controlling temperature in the BOMD simulations in canonical NVT ensemble so that the energy drift can be adsorbed by the fictitious heat bath and a constant temperature can be held. The NHC method is implemented in the Velocity Verlet34 method to control the temperature at 300.0 K, and the initial temperature is also set at 300.0 K. The total simulation period is 5.0 ps with time steps of 1.0 fs, and the trajectories are recorded every 20 steps. The conformation distributions and pair distribution functions discussed in the following section are sampled in the last 4.0 ps of the simulation time. 2.4. Computational Details of PEO and PE Models. Different lengths of PEO [H(CH2CH2O)nH, n ) 6, 8, 10, 15, 20] and PE [H(CH2CH2)nH, n ) 9, 12, 15, 21, 30] have been studied. Among them, PEO10 and PE9 are chosen to validate our fragmentation-based method by making comparison with conventional QM/MM and QM methods. As shown in Figure 1, the oligomer is divided into several fragments. Here, we chose every three successive heavy atoms along the main chains as a fragment. Then, the nearest neighboring fragments of this central fragment are included to form a work subsystem for QM calculations. In the next section, it will be shown that such a fragmentation scheme can give satisfactory results relative to the conventional QM method. Furthermore, the accuracy of the fragmentation-based QM calculations can be further improved by using longer fragments and hence larger subsystems, but the computational time will increase accordingly.15 Because we attempt to apply the GEBF QM/MM to MD simulations of

Figure 2. Force deviations (∆F) in unit of mEh/Å between the fragmentation and conventional QM/MM methods in PEO10 and PE9 solutions. The energy-based fragmentation QM calculations are carried out on small fragments (A and B) and larger fragments (C and D). The fragmentation schemes are also given in each figure, where the black balls denote O atoms in PEOn and the red bars denote the cut positions.

solutions of longer oligomers with up to 30 repeat units, we use the modest fragmentation scheme. Helix conformations of PEO are used as the initial geometries, and linear conformations are used for PE. Orthorhombic boxes of solvent with 12.0 Å buffers are added for a solute molecule to build the solvent shells, and the details are shown in Table 1. The initial positions of water molecules are randomly produced by the Amber package.35 The single-point charge36 model is employed for water solvents, where the partial charges of O and H are -0.82 and 0.41, respectively. A 12-6 Lennard-Jones potential is used for describing the van der Waals interactions. Parameters of consistent valence force field (CVFF)37 are employed, and all the force-field parameters are listed in the Supporting Information, Table S2. To validate the results of the fragmentation QM/ MM method, the conventional QM/MM is also carried out by using the Gaussian03 package38 and the same force-field parameters. 3. Results and Discussion First of all, the performance and applicability of the present fragmentation QM/MM method are tested in conformation energies, geometry optimizations, and BOMD simulations, in comparison with the conventional QM/MM and full QM results. And then, on the basis of fragmentation-based QM/MM simulations, the hydrogen bonding and chain configurations of PEO and PE as a function of chain length in aqueous solutions are investigated. 3.1. Conformation Energies. Two kinds of fragmentation schemes of PEO10 and PE9 are shown in Figure 2. The singlepoint energy calculations of PEO10 and PE9 (in the initial configurations) are performed in the gas phase, with the results shown in Table 2. The comparisons are also made between the fragmentation and conventional QM/MM models in Table 2. For PEO10 and PE9, single-point energies obtained by energybased fragmentation QM calculations have very small errors, no larger than 2.0 mEh (1.25 kcal/mol) at various levels with different basis sets. In the fragmentation-based QM/MM framework, it can be found that although thousands of molecules are added into the system, the accuracy of total energy is still


Li et al.

TABLE 2: Total Energies in Single-Point Calculations of PEO10 and PE9 Models in Gas Phase (full QM) and Aqueous Solutions (QM/MM), with Both the Fragmentationa and the Conventional Methods HF/sto-3g

B3LYP/3-21g

MP2/6-31g*

frag. (au)a conv. (au) ∆E (mH)

-1510.95493 (-1510.95506) -1510.95505 -0.12 (0.01)

PEO10 in gas phase -1530.95923 (-1530.95842) -1530.95843 0.80 (-0.01)

-1534.55310 (-1534.55299) -1534.55338 -0.28 (-0.39)


-1527.04082 (-1527.04087) -1527.04096 -0.14 (-0.09)

PEO10 in aqueous solution -1547.04809 (-1547.04733) -1547.04755 0.54 (-0.22)

-1550.64233 (-1550.64208) -1550.64229 0.04 (-0.21)


-695.46494 (-695.46493) -695.46492 0.02 (0.01)

PE9 in gas phase -704.93460 (-704.93284) -704.93250 2.10 (0.34)

-706.05444 (-706.05334) -706.05308 1.46 (0.26)


-683.23442 (-683.23440) -683.23440 0.02 (0.00)

PE9 in aqueous solution -692.70781 (-692.70599) -692.70569 2.12 (0.30)

-693.82707 (-693.82567) -693.82544 1.63 (0.23)

a Two kinds of fragmentation schemes are adopted in QM parts, as shown in Figure 2. The results in the parentheses are obtained from the fragmentation schemes with larger fragments (Figure 2C,D).

TABLE 3: Structures and Energies of PEO10 and PE9 in Aqueous Solution Optimized by Both Fragmentation and Conventional QM/MM Methodsa aqueous solution frag.

Figure 3. Atomic charge deviations (∆q) between the fragmentation (with schemes in Figure 2A,B) and conventional QM/MM methods of QM parts of PEO10 and PE9 backbones in aqueous solutions.

similar to that of the fragmentation QM model. Furthermore, the accuracy of energy increases with the increase of the fragment size, as shown in Table 2. Energy-based fragmentation QM/MM calculations also give a satisfactory description of forces. The force deviations between fragmentation (with the modest fragment size) and conventional QM/MM in PEO10 and PE9 solutions are shown in Figure 2A,B, respectively. It can be seen that the maximum error is within 1.0 × 10-3 Eh/Å. Obviously, the accuracy of the forces increases when the fragment size increases (see Figure 2C, D), but the needed computational effort also increases. To make a compromise between the accuracy and the computational cost, we select the fragmentation schemes shown in Figure 2A,B for the following geometry optimizations and MD simulations. The deviations in the QM part and the linkage region between the QM and the MM parts are larger than those in the MM part far from the QM part. It is not surprising because of the approximation made in the energy and force calculations in the fragmentation QM model. Furthermore, atomic partial charges of QM parts in PEO10 and PE9 obtained from the fragmentation and conventional QM/ MM calculations are also compared in Figure 3. The atomic charges derived from the GEBF method are also consistent with the conventional method (less than 0.01 e). The fragmentationbased QM/MM method can give a reasonable description of electrostatics interaction between the QM and the MM parts.

C-C (Å) 1.527 C-O (Å) 1.447 C-H (Å) 1.082 C-C-O (deg) 108.6 C-O-C (deg) 115.6 H-C-H (deg) 108.8 energy (au) -1560.3674

conv. P(EO)10 1.525 1.446 1.082 108.0 115.8 108.7 -1560.38457

gas phaseb crystalc,d 1.524 1.436 1.082 108.2 116.8 108.6

1.54 1.43 1.09 110 112 109.5

1.541 1.087 112.5 107.5

1.53 1.069 112.0 107.0

P(E)9 C-C (Å) C-H (Å) C-C-C (deg) H-C-H (deg) energy (au)

1.540 1.088 112.3 107.0 -724.40171

1.540 1.088 112.2 107.2 -724.41256

a Data from gas-phase full quantum chemical calculation and X-ray crystal experiment are also listed for comparison. b Full quantum optimized structures on HF/3-21G level. c X-ray diffraction of PEO crystal.39 d X-ray scattering result of PE at 4 K.40

In comparison with the time consumed in QM calculations, the CPU time for MM calculations can be ignored. Hence, the computational cost of the fragmentation QM/MM is dominated by the time consumed in the fragmentation-based QM calculation, and the effort needed for the GEBF method has been discussed a lot in previous works.15 It has been demonstrated that the fragmentation method can achieve the linear scaling with the number of basis sets for macromolecules.8–20 3.2. Geometry Optimizations. In order to save CPU time, HF/3-21G level is used for all the geometry optimizations and MD simulations. The same initial structures (helix for PEO and linear for PE) are used as starting points for both fragmentation and conventional QM/MM methods. The optimized structures and corresponding energies of PEO10 and PE9 are shown in Table 3 and Figure 4. Although there are innumerous minima in such systems, with flat potential surfaces, the optimized backbones of oligomers from the fragmentation and conventional QM/MM methods superpose well (Figure 4a). The average bond lengths and bond angles obtained by the frag-


Figure 4. Comparison between the fragmentation and conventional QM/MM optimizations. (a) Superposition of the optimized oligomer backbones from the two methods and (b) first solvation shells (r < 3.0 Å) of PEO10 (left) and PE9 (right) solutions.

mentation-based QM/MM method are consistent with those obtained with the conventional method, as shown in Table 3. Full QM optimizations of PEO10 and PE9 in the gas phase have also been performed to make a comparison with QM/MM methods.41 Compared with the PEO oligomer in the gas phase, the C-O bond length in aqueous solution increases obviously because of the polarization effect, but the covalent C-C bonds are hardly affected by the adjacent water molecules. The hydrogen bonding interactions between water molecules and the oxygen atoms in PEO make the C-O bonds weaker than those in the gas phase and crystal phase (where the solvent-water molecules are absent). Significant interaction between water molecules and oxygen atoms in PEO will be further quantitatively described by the MD simulations in Section 3.4. 3.3. BOMD Simulations. BOMD simulations have been performed for PEO10 and PE9 with both the fragmentation and the conventional QM/MM methods. As mentioned in Section 2, the Nose-Hoover-chain method has been employed to adsorb the drifted Hamiltonian energy in NVT ensembles. Because the C-O bonds in PEO are more polar than the corresponding C-C bonds in PE, the drift of Hamiltonian energy in PEO is slightly more notable than in PE systems, as shown in the Supporting Information, Figure S1. It implies that the error in GEBF gradients is mostly involved in the point-charge approximation for the distant parts from the central fragments. The temperatures and energies of PEO10 and PE9 in the framework of fragmentation-based QM/MM are comparable with those in the conventional method, as shown in Figure 5. Like most MD simulations, the starting points are at the minima of the simulated systems at 0 K; thus, the potential energies and total energies in NVT ensembles keep on rising in the first picosecond. After that period, the potential energies decrease slightly to reach the most probable configurations, and subsequently, the configurations and energies frustrate around the equilibrium positions. The superposed trajectories of the oligomer chains are shown in Figure 6, from which one can see that the snapshots from fragmentation and conventional QM/MM MD simulations are similar to each other. Usually, the local conformations of oligomer skeleton are classified into trans (T) and gauche (G/ G′, two gauche conformations with opposite dihedral angles) on the basis of the torsion angles, Φ, of the backbone atoms.

J. Phys. Chem. B, Vol. 112, No. 23, 2008 7065 Distributions of torsion angles (Φ) of the oligomer chains are shown in Figure 7. It is clearly seen that both fragmentation and conventional QM/MM simulations give nearly the same populations of T and G/G′ local conformations with Φ of around 0° and 110°/-110°, respectively. Therefore, the more economical fragmentation-based QM/MM can be used for MD simulations of aqueous solutions of longer oligomers (e.g., PEO20 and PE30, consisting of about 180 atoms in QM parts), without much loss of accuracy. 3.4. Intermolecular Hydrogen Bonds. In aqueous solutions, the hydrogen bonds, such as O-H · · · O and C-H · · · O, have been a long-standing topic in both experimental and theoretical works.42 Here, hydrogen bonding also plays an important role in the chain configurations in aqueous PEO and PE solutions. In the present work, we focus our attention to the hydrogen bonds between oligomers and solvent-water molecules. In the analysis of optimized structures of PEO10 and PE9 solutions, the hydrogen-bond criterion for the O · · · H distance is set as being less than 3.0 Å, and the O-H · · · O (or C-H · · · O) is larger than 110°, as suggested in other works.42f From the snapshots of the first solvation shells of PEO10 in Figure 4b and the optimized structure parameters listed in Table 4, significant O-H · · · O hydrogen bonds have been observed with a short O · · · H distance, and the coordinate number of the O-H · · · O hydrogen bond is one or two. Relatively weak C-H · · · O hydrogen bonds also exist in both PEO10 and PE9 in the optimized structures at 0 K, and their O · · · H distances are in good agreement with the neutron diffraction data of (NC)Csp3-H · · · O and (ClC)Csp3-H · · · O in solid states.42f A more detailed picture of hydrogen bonding can be obtained by the survey of intermolecular pair distribution functions from MD simulations. For this purpose, PEO oxygen/water hydrogen (OPEO-HW), PEO hydrogen/water oxygen (HPEO-OW), and PE hydrogen/water oxygen (HPE-OW) intermolecular pair distribution functions are shown in Figure 8. The fragmentation and conventional QM/MM methods give nearly identical curves for OPEO-HW, HPEO-OW, and HPE-OW pair distribution functions. In the OPEO-HW pair distribution function, the first aqueous shell of oxygen atoms in PEO is distinctly denoted by the peak at 1.7 Å. In contrast to the strong O-H · · · O hydrogen bonds, C-H · · · O hydrogen bonds are absent in the MD simulations under 300 K, as shown by the loose aqueous shell around PE9 (Figure 4b) and the flat shape of HPE-OW pair distribution curves without any specific peaks. We also show in Figure 9 the distributions of O-H · · · O (when H · · · O < 2.5 Å) and C-H · · · O (when H · · · O < 3.0 Å) angles in the vicinity of solute molecules. The sharp peak at 170° in the distribution curve of O-H · · · O angles again denotes the existence of specific O-H · · · O hydrogen bonds in aqueous PEO10 solution. Full QM calculations of supermolecular solvent cluster models, Models 1-3 (Supporting Information, Figure S2) with polarized continuum model (PCM)43–45 have also been performed for comparison. The QM calculations at HF/3-21G* level give 1.85 Å and 2.30 Å for H · · · O distances in O-H · · · O and C-H · · · O, respectively. Correspondingly, the O-H · · · O hydrogen bonds are nearly linear, and C-H · · · O bends slightly. The optimized structures from both fragmentation and conventional QM/MM methods give satisfactory descriptions of orientations and saturations of both O-H · · · O and C-H · · · O hydrogen bonds. In contrast, these two characteristics of intermolecular hydrogen bonding cannot be described well in most force fields. In force fields such as CVFF, hydrogenbonding interaction is usually considered as the combination of two-body electrostatics interaction and modified VdW


Li et al.

Figure 5. Comparison of temperatures (left) and energies (right) in BOMD simulations of PEO10 and PE9 at 300 K. In (a) and (b), the black and red lines denote temperatures in fragmentation and conventional QM/MM methods, respectively. In (c) and (d), the black lines and red lines denote potential energies and total energies (potential energy and kinetic energy) from fragmentation QM/MM method, respectively, and the green and blue lines are potential and total energies from the conventional method, respectively.

TABLE 4: Comparison of Intermolecular Hydrogen Bonds in PEO10 and PE9 Optimized by Fragmentation and Conventional QM/MM Methods and Parameters from Full QM Calculations on Cluster Models Embedded in PCM and Experiments for Solid States explicit QM/MM model conv.

PCMa

exptl.b

1.83 159.3 2.54 137.4

1.85 176.4

1.98

2.50 141.3

2.29 178.4

frag.

Figure 6. Superposed trajectories of MD simulations with an interval of 200 fs for PEO10 (left) and PE9 (right) by using the fragmentation (up) and conventional (down) QM/MM methods. The water molecules are omitted for clarity.

P(EO)10

OPEO · · · HW (Å) O-H · · · O (deg) HPEO · · · OW (Å) C-H · · · O (deg)

1.80 160.0 2.50 136.5

HPE · · · OW (Å) C-H · · · O (deg)

2.49 142.76

2.47

P(E)9 2.59

a Optimized geometries of models 1-3 optimized at HF/3-21G level with PCM model are shown in the Supporting Information, Figure S2. b Neutron diffraction data in solid states.42f

Figure 7. Distributions of overall dihedral angles (Φ) along the main chains of PEO10 (left) and PE9 (right) simulated by fragmentation and conventional QM/MM methods in NVT ensembles.

interaction. The atomic partial charges are usually obtained from QM calculation or experiments and fixed in the durations of MD simulations. Then, the LennardJones potential is modified to consider the hydrogen-bonding interaction. Therefore, the simultaneous polarization effects are neglected in classic MM, hindering its application to properly describe the hydrogen

bonding. In the present electrostatics embedding model of QM/ MM, the MM part appears as point charges in QM Hamiltonian, approximately accounting for the polarization effects from solvent molecules. However, new problems will arise in consideration of the mismatch of the MM parameters of partial charges and LennardJones potential with the QM electrostatics energy. For example, the QM/MM calculations of PEO aqueous solution yield QM charges of O, C, and H atoms of about -0.7, 0.12, and -0.23, respectively, much larger than the default CVFF parameters for O, C, and H partial charges (-0.3, -0.05, and 0.1, respectively). Therefore, the electrostatics interaction between QM and MM parts may overestimate the O-H · · · O hydrogen-bonding interaction. This also gives us a new hypothesis to develop new force-field parameters for QM/MM methods. 3.5. Length Dependence of Chain Dynamics and Hydrogen Bonding. In order to study the influence of chain length on the inter- and intramolecular interactions, we have performed the GEBF QM/MM MD simulations on PEO and PE with


Figure 8. Pair distribution functions of PEO10 and PE9 solutions in simulations by both fragmentation and conventional QM/MM methods at 300 K.

Figure 9. Distributions of O-H · · · OPEO and C-H · · · OW angles (θ) in the vicinity of solutes (H · · · O within 2.5 Å for O-H · · · OPEO and H · · · O within 3.0 Å for C-H · · · OW) in the simulations by fragmentation and conventional QM/MM methods at 300 K.

Figure 10. Distributions of the overall dihedral angles (Φ) along the main chains of PEO and PE oligomers with different lengths in the MD simulation (NVT, 300 K).

various lengths in aqueous solutions. The heavy computational load of QM/MM in treatment of PEO and PE oligomers with more than 100 backbone atoms limits the MD simulation within several picoseconds. Thus, we concentrate on the microscopic characteristics such as the length dependence of local conformations and motions of chains in the initial solvation stage. Local Conformations. The chain conformation is usually described by the dihedral angles along the main chain. Probability distributions of overall dihedral angles along the main chains of PEO and PE in the MD simulations are shown in Figure 10. There is little length dependence in the distributions of T and G/G′ conformations in both PEO and PE solutions. In PEO, the population of G conformations is larger than that of G′. However, in PE, the two peaks of G conformations have similar heights. The trans conformation is much more favorable than G in PE, and the ratio of G:T:G′ is about 1:5:1. The G conformations in PEO have more populations than those in PE.

J. Phys. Chem. B, Vol. 112, No. 23, 2008 7067

Figure 11. Distributions of OCCO and COCC dihedral angles in PEO10 during different time periods of simulations.

Figure 12. End-to-end distances of PEO (left) and PE (right) evolved with simulation time (NVT, 300 K).

In order to understand the dihedral angle distribution in PEO, detailed analysis is performed on the conformations of PEO10 chains. There are two different kinds of dihedral angles along the main chain of PEO, OCCO and COCC. The distributions of these two kinds of dihedrals in PEO10 updated along with the simulation time are shown in Figure 11. For OCCO, peaks of the two G conformations, G and G′, both increase and reach the same area, from the first to the last picosecond in the MD simulation. Unlike the distribution of CCCC dihedral angles in PE, G conformation is much preferable for OCCO dihedrals, in agreement with previous MD investigations within the framework of MM.24a For COCC dihedral angles, the distribution curve is quite different from that for OCCO, only the peak of the G conformation is distinct, and the peak of G′ is inconspicuous. This is closely related to the initial geometry of PEO. When the MD simulation is started, the initial G conformations of helix PEO chain transfer to trans conformations rapidly; at the same time, G′ conformations also appear. It can be predicted that the distribution of COCC dihedral angles will become a symmetrical shape such as the CCCC distribution in PE when the simulation time is long enough. Chain Dynamics. We also investigate the time dependence of end-to-end distance of PEO and PE, as shown in Figure 12. By starting from the helix and linear initial conformations for PEO and PE, all the end-to-end distances gradually decrease during the simulations. Hence, we define an average curling rate, Vc, as

VC )

d0 - dt t

(7)

where d0 and dt are the end-to-end distances at the initial (t ) 0) and the tth time step, respectively. The curling rates of PEO and PE with different lengths are shown in Figure 13. Both PEO and PE have a similar tendency for the evolvement of VC


Figure 13. Average curling rates of the evolvement of end-to-end distance of PEO and PE oligomer chains (NVT, 300 K).

Figure 14. Pair distribution functions of OPEO-HW in simulations of PEO with various lengths in aqueous solutions (NVT, 300 K).

versus the chain length, and the curling rate for PE is larger than that for PEO. When the number of atoms in the main chain is smaller than 20, the curling rate is very slow. As the chain length increases, the curling rate increases rapidly, implying that longer oligomer chains have larger chain flexibility and hence prefer to take a curl conformation. This is consistent with our knowledge of long oligomers and polymers. However, when the chain becomes longer than 30 atoms, VC increases much more slowly and even decreases when the chain contains more than 50 backbone atoms. This suggests that the curling occurs in some local regions in the whole polymer chain. Thus, when the chain length is larger than the local curling region, the curling rate will not increase linearly. Hydrogen Bonding. The pair distribution functions of OPEO-HW in various aqueous PEO solutions are nearly the same, as shown in Figure 14. This means that in the case of PEO solutions, the intermolecular hydrogen-bonding interaction has little length dependence in the initial salvation stage. 4. Conclusions We have implemented a hybrid QM/MM method combining the GEBF quantum chemical method with MM. The singlepoint calculations, geometry optimizations, and BOMD simulations of PEO and PE in dilute aqueous solutions have been performed to test the performance of the present GEBF QM/ MM method. The fragmentation QM/MM method can give satisfactory energies, forces, geometries, electrostatics potentials, and other characteristics, in comparison with the conventional QM/MM method. The present fragmentation QM/MM method has some advantages over conventional MM models. First, the polarization

Li et al. of solvents is included in the QM/MM method. The electrostatics embedding model considers the polarization by containing the solvent atomic charge in QM Hamiltonian. Second, the hydrogen-bonding interaction is well described in the fragmentation-based QM/MM method. Within the conventional MM, the hydrogen bonding is usually described by modified VdW and electrostatics interactions, which are of two-body nature. However, hydrogen bonding is a kind of many-body interaction. Thus, the two-body nonbonded interaction potential cannot give a reasonable description of some important characteristics of hydrogen bonds such as orientations and saturations, which are well described in our QM/MM simulations. Third, the QM/MM method does not rely too much on the transformability of forcefield parameters in various chemical environments. The present GEBF QM/MM simulations provide us detailed information on the chain dynamics and length dependence of local conformations and hydrogen bonding in PEO and PE solutions. In comparison with conventional QM/MM methods, the fragmentation-based method combines advantages of both the QM/ MM method and the molecular fragmentation method. The GEBF method can efficiently decrease the computational cost of QM calculations of long oligomers with up to 30 repeat units, and the QM/MM method can treat the interactions between the QM part and complex environments. The GEBF QM/MM method gives a new insight into studying macromolecule solutions in an efficient way. With this method, it is also possible to perform ab initio MD simulations of large systems such as polymers, biomolecules, and even bulk materials. The present work is just a demonstration of the fragmentationbased QM/MM method at the preliminary stage. The problem of maintaining Hamiltonian energy conservation requires more advanced strategy to precisely deal with the interaction between QM and MM parts. The overestimated O-H · · · O hydrogenbonding interaction calls for the reparameterization of some force-field parameters in the polarizable system. Additional functions such as periodic boundary condition will be also implemented to make the method applicable to more realistic models. We believe that, with the development of computer science, the low-scaling QM method combined with MM will become a practical tool in the simulations of macromolecule solutions in the future. Acknowledgment. This work was supported by the National Natural Science Foundation of China (Grant no. 20433020 and 20573050), the Ministry of Education (NCET-05-0442), and the National Basic Research Program (Grant no. 2004CB719901). Supporting Information Available: Hamiltonian energies in P(EO)10 and P(E)9 solutions during the MD simulations (Figure S1); selected clusters of full quantum chemistry optimized at HF/3-21G and MP2/6-31G* levels with PCM model (Figure S2); initial geometries of PEO10 and PE9 (QM parts) (Table S1); and force-field parameters in QM and MM parts (Table S2). This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Tomasi, J.; Persico, M. Chem. ReV. 1994, 94, 2027–2094 (and references therein) (2) Cramer, C. J.; Truhlar, D. G. Chem. ReV 1999, 99, 2161–2200 (and references therein) (3) Reichardt, C. Chem. ReV. 1994, 94, 2319–2358 (and references therein) (4) Orozco, M.; Luque, F. J. Chem. ReV. 2000, 100, 4187–4225. (and references therein) (5) See, for example: (a) Alema`n, C. Chem. Phys. 1999, 244, 151– 162. (b) Alema`n, C.; Galembeck, S. E. Chem. Phys. 1998, 232, 151–159.

Fragmentation-Based QM/MM Simulations (c) Alema`n, C. Chem. Phys. Lett. 1999, 302, 461–470. (d) Ekanayake, K. S.; Lebreton, P. R. J. Comput. Chem. 2006, 27, 277–286. (e) Kong, S.; Evanseck, J. D. J. Am. Chem. Soc. 2000, 122, 10418–10427. (f) Kiruba, G. S. M.; Wong, M. W. J. Org. Chem. 2003, 68, 2874–2881. (g) Mennucci, B.; Martı´nez, J. M. J. Phys. Chem. B 2005, 109, 9830–9838. (h) Bandyopadhyay, P.; Gordon, M. S. J. Chem. Phys. 2000, 113, 1104–1109. (i) Meng, S.; Ma, J.; Jiang, Y. J. Phys. Chem. B 2007, 111, 4128–4136. (j) Meng, S.; Ma, J. J. Phys. Chem. B 2008, ASAP (6) Warshel, A.; Levitt, M. J. Mol. Biol. 1976, 103, 227–249. (7) See, for example: (a) Singh, U. C.; Kollman, P. A. J. Comput. Chem. 1986, 7, 718–730. (b) Dewar, M. J. S.; Zoebisch, E. G.; Healy, E. A.; Stewart, J. J. P. J. Am. Chem. Soc. 1985, 107, 3902–3909. (c) Field, M. J.; Bash, P. A.; Karplus, M. J. Comput. Chem. 1990, 11, 700– 733. (d) Gao, J.; Xia, X. Science 1992, 258, 631–635. (e) Stanton, R. V.; Hartsough, D. S.; Merz, K. M., Jr. J. Phys. Chem. 1993, 97, 11868– 11870. (f) Thery, V.; Rinaldi, D.; Rivail, J.-L. J. Comput. Chem. 1994, 15, 269–282. (g) Maseras, F.; Morokuma, K. J. Comput. Chem. 1995, 16, 1170–1179. (h) Thompson, M. A. J. Phys. Chem. 1995, 99, 4794– 4804. (i) Eurenius, K. P.; Chatfield, D. C.; Brooks, B. R.; Hodoscek, M. Int. J. Quantum Chem. 1996, 60, 1189–1200. (j) Bakowies, D.; Thiel, W. J. Phys. Chem. 1996, 100, 10580–10594. (k) Bersuker, I. B.; Leong, M. K.; Boggs, J. E.; Pearlman, R. S. Int. J. Quantum Chem. 1997, 63, 1051–1063. (l) Gao, J.; Amara, P.; Alhambra, C.; Field, M. J. J. Phys. Chem. A 1998, 102, 4714–4721. (m) Lyne, P. D.; Hodoscek, M.; Karplus, M. J. Phys. Chem. A 1999, 103, 3462–3471. (n) Eichinger, M.; Tavan, P.; Hutter, J.; Parrinello, M. J. Chem. Phys. 1999, 110, 10452–10467. (o) Zhang, Y.; Lee, T-S.; Yang, W. J. Chem. Phys. 1999, 110, 46–54. (p) Zhang, Y.; Liu, H.; Yang, W. J. Chem. Phys. 2000, 112, 3483–3492. (q) Field, M. J.; Albe, M.; Bret, C.; Proust-De Martin, F.; Thomas, A. J. Comput. Chem. 2000, 21, 1088–1100. (r) Yarne, D. A.; Tuckerman, M. E.; Martyna, G. J. J. Chem. Phys. 2001, 115, 3531–3539. (s) Laio, A.; VandeVondele, J.; Rothlisberger, U. J. Chem. Phys. 2002, 116, 6941– 6947. (t) Crespo, A.; Scherlis, D. A.; Marti, M. A.; Ordejon, P.; Roitberg, A. E.; Estrin, D. A. J. Phys. Chem. B 2003, 107, 13728–13736. (u) Sherwood, P.; et al. J. Mol. Struct. 2003, 632, 1–28. (v) Ferre, N.; Olivucci, M. J. Am. Chem. Soc. 2003, 125, 6868–6869. (w) Loferer, M. J.; Loeffler, H. H.; Liedl, K. R. J. Comput. Chem. 2003, 24, 1240–1249. (x) Nam, K.; Gao, J.; York, D. M. J. Chem. Theory Comput. 2005, 1, 2–13. (y) Gao, J.; Ma, S.; Major, D. T.; Nam, K.; Pu, J.; Truhlar, D. G. Chem. ReV. 2006, 106, 3188–3209. (8) (a) Yang, W. Phys. ReV. Lett. 1991, 66, 1438–1441. (b) Yang, W. Phys. ReV. A 1991, 44, 7823–7826. (c) Lee, C.; Yang, W. J. Chem. Phys. 1992, 96, 2408–2411. (d) Yang, W.; Lee, T.-S. J. Chem. Phys. 1995, 103, 5674–5678. (e) Zhao, Q.; Yang, W. J. Chem. Phys. 1995, 102, 9598–9603. (f) Lee, T-S.; York, D. M.; Yang, W. J. Chem. Phys. 1996, 105, 2744– 2750. (g) York, D. M.; Lee, T-S.; Yang, W. Phys. ReV. Lett. 1998, 80, 5011–5014. (9) Dixon, S. L.; Merz, K. M., Jr. J. Chem. Phys. 1996, 104, 6643– 6649. (b) Dixon, S. L.; Merz, K. M., Jr. J. Chem. Phys. 1997, 107, 879– 893. (c) Monard, G.; Bernal-Uruchurtu, M. I.; van der Vaart, A.; Merz, K. M., Jr.; Ruiz-Lo´pez, M. F. J. Phys. Chem. A 2005, 109, 3425–3432. (10) (a) Mezey, P. G. J. Math. Chem. 1995, 18, 141–168. (b) Mezey, P. G. AdV. Quantum Chem. 1996, 27, 163–222. (c) Mezey, P. G. Int. J. Quantum Chem. 1997, 63, 39–48. (d) Mezey, P. G. Int. ReV. Phys. Chem. 1997, 16, 361–388. (e) Exner, T. E.; Mezey, P. G. J. Comput. Chem. 2003, 24, 1980–1986. (f) Exner, T. E.; Mezey, P. G. J. Phys. Chem. A 2004, 108, 4301–4309. (11) (a) Zhang, D. W.; Zhang, J. Z. H. J. Chem. Phys. 2003, 119, 3599– 3605. (b) Chen, X. H.; Zhang, D. W.; Zhang, J. Z. H. J. Chem. Phys. 2004, 120, 839–844. (c) Chen, X. H.; Zhang, J. Z. H. J. Chem. Phys. 2004, 120, 11386–11391. (d) Xiang, Y.; Zhang, D. W.; Zhang, J. Z. H. J. Comput. Chem. 2004, 25, 1431–1437. (e) Zhang, D. W.; Xiang, Y.; Zhang, J. Z. H. J. Phys. Chem. B 2003, 107, 12039–12041. (f) Zhang, D. W.; Chen, X. H.; Zhang, J. Z. H. J. Comput. Chem. 2003, 24, 1846–1852. (g) Chen, X. H.; Zhang, Y.; Zhang, J. Z. H. J. Chem. Phys. 2005, 122, 184105–184109. (h) Mei, Y.; Zhang, D. W.; Zhang, J. Z. H. J. Phys. Chem. A 2005, 109, 2–5. (i) Chen, X. H.; Zhang, J. Z. H. J. Chem. Phys. 2006, 125, 044903–044909. (12) Gu, F. L.; Aoki, Y.; Korchowiec, J.; Imamura, A.; Kirtman, B. J. Chem. Phys. 2004, 121, 10385–10391. (13) (a) Kitaura, K.; Ikeo, E.; Asada, T.; Nakano, T.; Uebayasi, M. Chem. Phys. Lett. 1999, 313, 701–706. (b) Kitaura, K.; Sugiki, S.-I.; Nakano, T.; Komeiji, Y.; Uebayasi, M. Chem. Phys. Lett. 2001, 336, 163–170. (c) Fedorov, D. G.; Kitaura, K. Chem. Phys. Lett. 2004, 389, 129–134. (d) Fedorov, D. G.; Kitaura, K. J. Chem. Phys. 2004, 120, 6832–6840. (e) Fedorov, D. G.; Kitaura, K. J. Chem. Phys. 2005, 122, 054108–054117. (f) Fedorov, D. G.; Kitaura, K. J. Chem. Phys. 2005, 123, 134103–134113. (14) Hirata, S.; Valiev, M.; Dupuis, M.; Xantheas, S. S.; Sugiki, S.; Sekino, H. Mol. Phys. 2005, 103, 2255–2265. (15) (a) Li, W.; Li, S. J. Chem. Phys. 2004, 121, 6649–6657. (b) Li, S.; Li, W.; Fang, T. J. Am. Chem. Soc. 2005, 127, 7215–7226. (c) Li, W.; Fang, T.; Li, S. J. Chem. Phys. 2006, 124, 154102–154107. (d) Li, W.; Li, S.; Jiang, Y. J. Phys. Chem. A 2007, 111, 2193–2199. (e) Cao, H.; Fang, T.; Li, S.; Ma, J. Macromolecules 2007, 40, 4363–4369. (f) Li, S.; Li, W.;

J. Phys. Chem. B, Vol. 112, No. 23, 2008 7069 Fang, T.; Ma, J.; Jiang, Y. Low Scaling Quantum Chemical (LSQC) program, version 1.1; Nanjing University: Nanjing, 2006. (g) Jiang, N.; Ma, J.; Jiang, Y. J. Chem. Phys. 2006, 124, 114112–114120. (16) (a) Deev, V.; Collins, M. A. J. Chem. Phys. 2005, 122, 154102– 154113. (b) Collins, M. A.; Deev, V. A. J. Chem. Phys. 2006, 125, 104104– 104118. (17) (a) Ganesh, V.; Dongare, R. K.; Balanarayan, P.; Gadre, S. R. J. Chem. Phys. 2006, 125, 104109–104118. (b) Gadre, S. R.; Shirsat, R. N.; Limaye, A. C. J. Phys. Chem. 1994, 98, 9165–9169. (c) Babu, K.; Gadre, S. R. J. Comput. Chem. 2003, 24, 484–495. (18) (a) Dahlke, E. E.; Truhlar, D. G. J. Chem. Theory Comput. 2007, 3, 46–53. (b) Dahlke, E. E.; Truhlar, D. G. J. Chem. Theory Comput. 2008, 4, 1–6. (c) Dahlke, E. E.; Leverentz, H. R.; Truhlar, D. G. J. Chem. Theory Comput. 2008, 4, 33–41. (19) (a) Morita, S.; Sakai, S. J. Comput. Chem. 2001, 22, 1107–1112. (b) Sakai, S.; Morita, S. J. Phys. Chem. A 2005, 109, 8424–8429. (20) Bettens, R. P. A.; Lee, A. M. J. Phys. Chem. A 2006, 110, 8777– 8785. (21) Sherwood, P. In Modern and Algorithm of Quantum chemistry, 2nd ed.; Grotendorst, J. Ed.; NIC: Ju¨lich, 2000; Vol. 3, p 285. (22) Body, R. W.; Kyllingstad, V. L. Encyclopedia of Polymer Science and Engineering, 2nd ed.; Wiley-Interscience: New York, 1986; Vol. 6; p 225. (23) Aggarwal, S. L.; Sweeting, O. J. Chem. ReV. 1957, 57, 665–742. (24) See, for example: (a) Rigby, D.; Sun, H.; Eichinger, B. E. Polym. Int. 1997, 44, 311–330. (b) Smith, G. D.; Bedrov, D.; Borodin, O. Phys. ReV. Lett. 2000, 85, 5583–5586. (c) Ahlstro¨m, P.; Borodin, O.; Wahnstro¨m, G.; Wensink, E. J. W.; Carlsson, P.; Smith, G. D. J. Chem. Phys. 2000, 112, 10669–10679. (d) Smith, G. D.; Bedrov, D.; Borodin, O. J. Am. Chem. Soc. 2000, 122, 9548–9549. (e) Bedrov, D.; Smith, G. D. Langmuir 2006, 22, 6189–6194. (f) Sasanuma, Y.; Ohta, H.; Touma, I.; Matoba, H.; Hayashi, Y.; Kaito, A. Macromolecules 2002, 35, 3748–3761. (g) Siqueira, L. J. A.; Ribeiro, M. C. C. J. Chem. Phys. 2006, 125, 214903–214910. (25) See, for example: (a) Fujiwara, S.; Sato, T. J. Chem. Phys. 1997, 107, 613–622. (b) Harmandaris, V. A.; Mavrantzas, V. G.; Theodorou, D. N. Macromolecules 1998, 31, 7934–7943. (c) Jin, Y.; Boyd, R. H. J. Chem. Phys. 1998, 108, 9912–9923. (d) Liu, C.; Muthukumar, M. J. Chem. Phys. 1998, 109, 2536–2542. (e) Pu¨tz, M.; Kremer, K.; Grest, G. S. Europhys. Lett. 2000, 49, 735–741. (f) Frankland, S. J. V.; Caglar, A.; Brenner, D. W.; Griebel, M. J. Phys. Chem. B 2002, 106, 3046–3048. (g) Hirvi, J. T.; Pakkanen, T. A. J. Chem. Phys. 2006, 125, 144712–144722. (h) Knyazev, V. J. Phys. Chem. A 2007, 111, 3875–3883. (26) Reed, A. E.; Curtiss, L. A.; Weinhold, F. Chem. ReV. 1988, 88, 899–926. (27) Nocedal, J. Math. Comput. 1980, 35, 773–782. (28) Liu, D. C.; Nocedal, J. Math. Prog. 1989, 45, 503–528. (29) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1987; p 99. (30) Kerdcharoen, T.; Liedl, K. R.; Rode, B. M. Chem. Phys. 1996, 211, 313–323. (31) Morokuma, K.; Kerdcharoen, T. Chem. Phys. Lett. 2002, 355, 257– 262. (32) Hofer, T. S.; Pribil, A. B.; Randolf, B. R.; Rode, B. M. J. Am. Chem. Soc. 2005, 127, 14231. (33) Martyna, G. J.; Klein, M. L.; Tuckerman, M. J. Chem. Phys. 1992, 97, 2635–2643. (34) Jang, S.; Voth, G. A. J. Chem. Phys. 1997, 107, 9514–9527. (35) Case, D. A.; Darden, T. A.; Cheatham, T. E., III. AMBER8; University of California: San Francisco, 2004. (36) Berendsen, H. J.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. In Intermolecular Forces; Pullman, B. Eds.; Reidel: Dordrecht, 1981, p331. (37) Dauber-Osguthorpe, P.; Roberts, V. A.; Osguthorpe, D. J.; Wolff, J.; Genest, M.; Hagler, A. T. Proteins: Struct., Funct., Genet. 1988, 4, 31–47. (38) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. 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 D.01; Gaussian, Inc.: Wallingford, CT, 2004. (39) Takahashi, Y.; Tadokoro, H. Macromolecules 1973, 6, 672–675.

7070 J. Phys. Chem. B, Vol. 112, No. 23, 2008 (40) Shukla, A.; Dolg, M.; Stoll, H. Chem. Phys. Lett. 1998, 294, 126– 134. (and references therein) (41) The seemingly good agreement between the gas-phase QM results and the crystal structures, shown in Table 3, may benefit from an error cancelation, because the relatively low theory level of HF/3-21G is employed here. (42) See, for example: Kollman, P. A.; Allen, L. C. Chem. ReV. 1972, 72, 283–325. (b) Berkovitch-Yellin, Z.; Leiserowitz, L. J. Am. Chem. Soc. 1982, 104, 4052–4064. (c) Taylor, R.; Kennard, O. Acc. Chem. Res. 1984, 17, 320–326. (d) Steiner, T.; Desiraju, G. R. Chem. Commun. 1998, 891–

Li et al. 892. (e) Novoa, J. J.; Lafuente, P.; Mota, F. Chem. Phys. Lett. 1998, 290, 519–525. (f) Steiner, T. Angew. Chem., Int. Ed. 2002, 41, 48–76. (43) Cance`s, M. T.; Mennucci, B.; Tomasi, J. J. Chem. Phys. 1997, 107, 3032–3041. (44) Cossi, M.; Barone, V.; Mennucci, B.; Tomasi, J. Chem. Phys. Lett. 1998, 286, 253–260. (45) Mennucci, B.; Tomasi, J. J. Chem. Phys. 1997, 106, 5151–5158.

JP800777E

MM Simulations: Length Dependence of

Recommend Documents