J. Phys. Chem. B 2008, 112, 7403–7409
7403
Molecular Dynamics Simulations of Heptyl Phosphonic Acid: A Potential Polymer Component for Fuel Cell Polymer Membrane Sudip Roy,* Tamer M. Ataol, and Florian Mu¨ller-Plathe Eduard-Zintl Institut fu¨r Anorganische and Physikalische Chemie, Technische UniVersita¨t Darmstadt, Petersenstr.-20, 64287 Darmstadt, Germany ReceiVed: July 20, 2007; ReVised Manuscript ReceiVed: April 7, 2008
Atomistic molecular dynamics simulations have been performed on heptyl phosphonic acid (HPA) to understand the dynamic hydrogen bonding network in the liquid phase. HPA is a phosphonic-acid functionalized alkane (heptane) and a model compound for one of the promising polymers for high temperature (>100 °C) fuel cell polymer electrolyte membranes. For the simulation, a force field for this molecule has been generated with the help of quantum chemical calculations and optimized by simplex algorithm. The force field has been validated against experimentally measured properties, for example, density and self-diffusion constant. From molecular dynamics simulations conducted at different temperatures, we have confirmed the hypothesis of dynamic hydrogen bond network formation in this material. Introduction Polymer electrolyte membrane (PEM) fuel cells are getting immense scientific attention due to their potential use in automobiles as a replacement of combustion engines. The polymer electrolyte membrane is a thin membrane, which separates the cathode and the anode in the fuel cell and helps the protons to be transferred from cathode to anode. Therefore, in a fuel cell the PEM is a crucial component, which determines the efficiency of the cell by its proton conductivity. It has to be chemically as well as mechanically stable against the fuel cell’s electrochemical environment. In modern PEM fuel cell technology, the polymer Nafion is widely used as membrane. Nafion is a good proton conductor under wet conditions. Experimental and theoretical studies have confirmed that the proton conductivity performance of Nafion is directly related to the degree of hydration of the membrane. Therefore, this material works efficiently only below the boiling temperature of water, because above it water starts to evaporate and the efficiency of the fuel cell decreases. To keep the high efficiency of the fuel cell at elevated temperatures, much effort has been focused on different strategies, for example, addition of high boiling-point proton conductive liquid such as imidazolium triflate1 or ionic liquids2 as well as the development of new membrane materials, which have an intrinsic capability of carrying protons at high temperatures. Other common fuel cells, apart from PEM, use liquid phosphoric acid as electrolyte. The operating temperature of phosphoric acid fuel cells is between 150 and 200 °C. The amphoteric liquid phosphoric acid shows high conductance with a high contribution from structure diffusion, that is, intermolecular proton transfer due to the hydrogen bond breaking and formation processes.3 So, on one hand we have heterocyclic compounds as carriers of protons at high temperatures and on the other hand we have liquid phosphoric acid. But both compounds have their own limitations as integrated components of a fuel cell. The moderate amphotericity of heterocyclic imidazole, pyrazole, and benzimadazole has fewer exchangeable protons, which participate in proton conductance and, hence, * To whom correspondence should be addressed. E-mail: s.roy@ theo.chemie.tu-darmstadt.de.
Figure 1. Hepthylphosphonic acid (HPA).
the conductivity is less than in the presence of water. Another problem of these materials is their adsorption on the platinum surface of electrode and a reduction of the electrode’s electrochemical efficiency. Therefore, there is a need for new materials to be investigated which can be used as PEM in high-temperature fuel cells. Among many other proposed materials, organic phosphonic acids and polymers made thereof as membrane material are attractive because of their high conductivity due to structural diffusion of hydrogen. Neat liquid phosphonic acid has a high proton conductivity due to its structure, diffusion and degree of self-dissociation. In the case of polymers containing phosphonic acid, a high proton conductivity could be achieved by introducing flexible spacers, for example, alkane or ethyleneoxide fragments, through which proton-carrying groups are attached to the polymer backbone. The flexibility of the long spacer groups facilitates the aggregation of the terminal protogenic groups and allows intermolecular hydrogen bonding. These spacer segments accelerate the hydrogen bond breaking and forming process by allowing better reorientation.4 Before polymers functionalized with phosphonic acid can be used in fuel cell as membranes, the mechanical, electrical (conductance) properties at high temperatures must be understood. Therefore, it is important to start calibrating the system from scratch using different levels of theory and experiment. The spacer concept, which has been shown to work well for heterocyclic compounds, may also play a role in proton conductance for phosphonic-acid-functionalized polymers, or monomer liquids.4 In this contribution, we have selected a model compound for our atomistic molecular dynamic simulation, which is phosphonic-acid-functionalized heptane or heptylphosphonic acid
10.1021/jp0757107 CCC: $40.75 2008 American Chemical Society Published on Web 05/31/2008
7404 J. Phys. Chem. B, Vol. 112, No. 25, 2008
Roy et al.
TABLE 1: Force Field Parameters atoms
ε (kJ/mol)
σ (nm)
atom (description)
charge (|e|)
bond
distance (nm)
H (bonded to C) H (bonded to O) All C O P
0.189 0.189 0.299 0.677 1.8
0.26582 0.0 0.32766 0.288 0.338
C (bonded to P) P H (bonded to O) O (bonded to H) H (bonded to charged C) O (doubly bonded) all others
-0.25 1.28 0.47 -0.73 0.135 -0.78 0
C-H C-C C-P PdO P-O O-H
0.109 0.1526 0.1795 0.15089 0.1557 0.094
bond angle
φ0 (degrees)
kθ (kJ/(mol rad2))
torsion
τ0[deg]-periodicity-kτ[kJ/(mol rad2)]
C-C-H H-C-H C-C-C P-C-H (1) P-C-H (2) C-P-O C-PdO O-P-O O-P ) O P-O-H
109.5 107.9 111.0 106.64 108 108.23 111.13 105.17 112.76 113.8
366.9 306.4 482.3 390.0 390.0 610 610 420 420 390
H-C-C-C C-C-C-C C-C-C-P C-C-P-O C-P-O-H
180s3-11.7 180s3-11.7 180s3-11.7 176s3-10.5 -142s1-19.4
(HPA, Figure 1) for which there are experimental reference data.5 The alkyl chain of this model compound introduces the spacer concept, for example, the heptylphosphonic acid can be grafted to a polymer backbone which may give mechanical stability to the membrane material. Theoretical elucidation of the proton transfer mechanism at the atomic level means that one has to model hydrogen bonding and its dynamics by bond breaking and bond forming processes. Ideally, one needs to carry out atomistic simulations at the ab initio level. But for the HPA fluid, ab initio molecular dynamics simulations would be prohibitively expensive. There are some attempts made on molecular phosphonic acid (H3PO3) with the density functional tight binding molecular dynamics method,6 which suggests a low energy barrier for proton transfer due to intermolecular hydrogen bonding. The ab initio based molecular modeling by Paddison et al.7 suggested that the energy barriers for proton transfer between pairs for methyl phosphonic acid and imidazole molecules are 37.2 and 120.1 KJ mol-1, respectively. Therefore, phosphonic acid has better prospects for fast hydrogen bonding dynamics. In this work, we have carried out classical molecular dynamics simulation of liquid heptyl phosphonic acid (HPA) structure depicted in Figure 1 to elucidate the mechanism of hydrogen bond network formation with the help of cluster analysis sampled over several nanoseconds. This type of analysis helps to understand the network formation processes, which may lead to a static or dynamic percolating networks inside the liquid sample and hence the enhanced proton conductivity. HPA is a new compound for atomistic molecular dynamics studies. There is no force field reported in the literature for this molecule which reproduces thermodynamical properties. Therefore, the major part of this paper is devoted to the force field development and validation. This is very important for future studies of phosphonic acid functionalized polymers with different spacer sizes. The paper is organized as follows. In the subsequent section, we describe the computational method used for the calculations. The section is followed by the results divided into two parts, first of which is devoted to the force field development and its validation and the other to the hydrogen bond network formation, that is, mainly its static and dynamic structural properties. Finally, the last section summarizes our findings.
Computational Details Molecular Dynamic. We performed classical molecular dynamics (MD) simulations of liquid HPA. All simulations are
Figure 2. C-C-P-O torsional potential fitting.
Figure 3. C-P-O-H torsional potential fitting.
Molecular Dynamics Simulations of Heptyl Phosphonic Acid
Figure 4. Density of liquid heptyl phosphonic acid at different temperatures.
J. Phys. Chem. B, Vol. 112, No. 25, 2008 7405 calculations, first- and second-neighbor intramolecular nonbonded interactions have been excluded. In eq 1, the qi are the partial charges of the individual atoms. For the initial guess of the charges, we have done quantum chemistry calculations with the Gaussian 03 code.12 For this purpose, geometry optimization of one HPA molecule has been done by using the Hartree-Fock method with the 6-31G* basis set, trying five randomly chosen different starting configurations. After the geometry optimization, we used CHELPG method13 to find the partial charges of the individual atoms in the molecule. These charges, however, were changed for the alkyl part and set to zero (except carbon attached to phosphorus). All the final parameters are given in Table 1. Besides nonbonded interactions, the bonded interactions are implemented in YASP code. The bond lengths (Table 1) were kept rigid during the simulation by the SHAKE procedure.14 The angle and torsional potentials are written in the following form
kφ (φ - φ0)2 2 angles
(3)
kτ [(1 - cos p(τ - τ0)] 2 torsions
(4)
Uangles ) Utorsions )
Figure 5. Potential energy of liquid heptyl phosphonic acid at different temperatures.
done by YASP,8 a molecular dynamics code developed in our group. In YASP, the form of the nonbonded potential is energy is
U)
[( ) ( ) ]
∑ 4εij ij
σij rij
12
-
σij rij
6
+
(
qiqj 1 εRF - 1 rij2 + 4πε0 rij 2εRF + 1 r3
cutoff
)
(1) where rij is the distance between the atoms i and j, ε0 is the vacuum permittivity, εRF ) 6, the reaction field dielectric constant. εij and σij are the Lennard-Jones parameters between atoms i and j and are obtained from the εii, εjj, σii, and σjj by the mixing rules
εij ) √εiiεjj
and
σij )
σii + σjj 2
(2)
The initial values of the Lennard-Jones parameters for the atoms of HPA have been taken from GROMOS9 and from the literature: εii and σii for the heptyl part were taken from ref 10, and εii and σjj of oxygen were taken from ref 11. Some of the Lennard-Jones parameters were changed during the simplex optimization process (details below) to fit the experimental density and diffusion coefficient data to the simulation. In all
∑
∑
In these equations kφ and kτ are angle and torsional force constants, respectively, for each different angle and dihedral in the molecule, and φ0 and τ0 are the equilibrium values; p is the periodicity of the torsional potential. The angle potentials of the alkane part were taken from ref 10 and from GROMOS.9 We have used the Hartee-Fock method with 6-31G* basis functions for the calculations of the torsional potentials for phosphonic acid part (last two torsions in Table 1). In these calculations, first the geometry of the molecule was optimized. Then, the torsion angle in question was changed in steps of 18 degrees. After each rotation the geometry was optimized again and the energy was recorded. The energy as a function of torsional angle was fitted by the torsional potential function of eq 4 and the required parameters were determined from this fitting. The fitting has been shown in Figure 2 and 3 for C-C-P-O and C-P-O-H torsions, respectively. The fitting has been simplified because of unnecessary higher periodicity in cosine potential which makes the calculations heavier particularly for simplex optimization procedure. After the fitting, the Lennard-Jones σ of carbon and hydrogen of alkyl chain has been optimized against the target density and diffusion. All force-field parameters are listed in Table 1. Simplex Optimization. The chosen set of parameters is optimized according to an iterative algorithm, that is, cycles of MD-simulation, evaluation of target properties, and eventually readjustment of parameters according to the downhill simplex method are performed repeatedly until the abort condition is reached.15 The readjustment of parameters is achieved by reflection, expansion, or contraction in parameter space.16 This optimization method is robust and applicable even if an analytical functional form is not available, as is the case here. For starting the simplex optimization procedure, initial values for the variables have to be provided manually. The values chosen should cover quite a large (but still reasonable) region in parameter space in order to avoid the results being trapped in a local minimum all over the time. After parametrization of force field constant pressure (NPT) simulations of a HPA liquid with 500 molecules have been performed at five different temperatures: 300, 350, 400, 450, and 500 K. The Berendsen thermostat17 was employed with a
7406 J. Phys. Chem. B, Vol. 112, No. 25, 2008
Roy et al.
TABLE 2: Cluster Evolution Patternsa clusters of Figure 12a
clusters of Figure 12b
clusters of Figure 12c
picoseconds snapshots
S
L
N
S
L
S
L
N
S
L
N
1ps 2ps 3ps 4ps
112 65 21 24
0 73 44 8
0 26 0 5
clusters 92 47 6 27
from Figure 12 0 0 45 0 40 1 0 21
97 6 17 30
0 97 6 3
0 6 17 16
102 47 20 25
0 54 27 14
0 1 0 9
a
N
clusters of Figure 12d
Size of cluster ) S. Number of members lost from previous snapshot ) L. Number of new members arrived after previous snapshot ) N.
Figure 6. Comparison of simulated and experimental densities at different temperatures. Experimental data are taken from ref 5.
Figure 8. Intermolecular radial distribution function of phosphorus atoms (solid curve) with integral (dotted curve).
Figure 7. Comparison of simulated and experimental self-diffusion coefficients at different temperatures. Experimental data are taken from ref 5.
Figure 9. Mixed intermolecular radial distribution function of doublebonded oxygen atoms and hydroxyl hydrogen atoms.
coupling time of 0.5 ps. The pressure was kept constant during the simulation at atmospheric pressure of 101.3 kPa by the Berendsen manostat17 with a coupling time of 5 ps. At each temperature, the liquid has been simulated for 5 ns. The cutoff and neighbor list cutoff for nonbonded interactions were taken as 0.9 and 1.0 nm, respectively. The system was equilibrated for 1 ns. After this, the density was constant (Figure 4). Following equilibration, a 4 ns simulation has been continued for the production simulation. The analysis of the simulation has been done over the recorded (every 1 ps) coordinates and velocities of the atoms. For analysis, the last 1 ns of the production simulation has been
used. The density (Figure 4) and total energy (Figure 5) for the last 2 ns of the simulation show equilibration of the density and the energy. Results and Discussion Force Field Validation. The force field given in Tables 1 and 2 has been obtained by changing the initial values of few parameters iteratively by simplex. The simplex algorithm implemented in YASP as published by Faller et al.15 have been applied to optimize the parameter Lennard-Jones σ of aliphatic carbon and hydrogen. These two parameters are optimized to
Molecular Dynamics Simulations of Heptyl Phosphonic Acid
Figure 10. Size distribution (unnormalized) for hydrogen-bonded clusters in heptyl phosphonic acid at 450 K.
get the density in accordance with experimental value. In each optimization cycle, at least three different simulations have been carried out: (1) A pre-equilibration of the starting configuration (which is the output coordinate of the previous simplex steps with most similar parameter, or already equilibrated system with initial parameter set) for 0.5 ns. (2) True relaxation of the system which was achieved previous step by running another 0.5 ns or until the equilibration criterion was met. The system was considered sufficiently equilibrated as soon as the slope of the linear regression of the pressure times the number of data points was smaller than their standard deviation.
J. Phys. Chem. B, Vol. 112, No. 25, 2008 7407 (3) Starting with the equilibrated system production run was performed for another 0.5 ns with time step 1 fs. After this sequence of simulations, the performance of the current parameter set has been calculated by comparing the obtained density to a target liquid density of 1.018 g/cm3 and the obtained self-diffusion constant to a target of 5.56 × 10-7 cm2/sec at 394 K temperature. These target values were obtained from the experimental work of Schuster et al.5 on HPA. With the resulting performance value, the simplex steps were executed, that is, a new set of parameters was generated and a new cycle was started. This procedure was repeated until the target density is reproduced within an error of 1% and target diffusion constant within an error of 10%. The aim of the optimization simulations was to fit the simulated data to the experimental data of HPA. Only the density and the diffusion coefficient have been used, since the heat of vaporization cannot be reliably measured for HPA, because of condensation occurring just a few degrees above the melting point, which is around 105 °C.5 The simulated densities along with experimental values are given in Figure 6. There is a good match (>1% deviation) of experimental and simulated data except around 325 K. Below 380 K, HPA freezes in experiment with a concomitant increase of the density. In the simulation, crystallization is not observed, but supercooling of the liquid. The force field is also optimized against the self-diffusion coefficients. Here, the simulated results are slightly lower than the experimental values (Figure 7). From the slope of the Arrhenius plot (Figure 7), the activation energy for diffusion can be obtained. It is 41 ( 1.5 kJ/mol from simulation, whereas the experimental value is 34.5 ( 0.4 kJ/mol. Interestingly, the simulation reproduces this energetic quantity quite well, although no energetic property, for example, heat of vaporization, was
Figure 11. A cluster of hydrogen-bonded heptyl phosphonic acid molecules. The formation of a channel of phosphonic-acid moieties is clearly visible.
7408 J. Phys. Chem. B, Vol. 112, No. 25, 2008
Roy et al.
Figure 12. Cluster evolution: (a) up left, (b) up right, (c) down left, (d) down right.
used in the force field parametrization. This gives additional confidence in the force field. Static Structural Properties. The proton conductivity of pure liquid HPA should be governed by the intermolecular hydrogen bonding. Hydrogen bonds facilitate the hopping of protons between phosphonic acid groups (structural diffusion). While proton-hopping is not included in our model, the formation and breakage of hydrogen bonds is. Therefore, MD simulations are reliable to elucidate the local structural and dynamical properties of also this aspect of the liquid. For HPA molecules doublebonded oxygen atoms and hydroxyl groups are taking part in the intermolecular hydrogen bonding. Thus, radial distribution functions between the atoms of phosphonic acid groups of different molecules give a qualitative understanding of hydrogen bonding and phosphonic acid aggregation. The RDFs for two different sets of atoms, first P atoms around P atoms and second H atoms of the hydroxyl group around the doubly bonded O atoms, are shown in Figures 8 and 9, respectively. In Figure 8, the first peak corresponds to the first solvation shell of the phosphonic acid group of one molecule. The integral over this peak until the first minimum (0.65 nm) gives the number of phosphonic acid groups forming this solvation shell. The integral of the RDF (Figure 8) shows that until the first minimum of the RDF there are approximately 6 P atoms around a given one. Therefore, there is a clustering of phosphonic acid groups around each other in the material. In Figure 9, the RDF shows only the intermolecular distances
between double bonded oxygen and hydrogen of hydroxyl group. There is a small peak observed at around 0.21 nm and the first minimum at 0.25 nm. This small peak suggests that irrespective of much favorable intramolecular interaction between double bonded oxygen and hydrogen of hydroxyl group there is still some intermolecular bridging at distance of 0.25 nm a typical hydrogen bond distance for such molecules. The second solvation shell of intermolecular distance for the same is also highly pronounced. Hydrogen Bond Percolation and Cluster Dynamics. The proton conduction benefits from the presence of large hydrogenbonded clusters of phosphonic acid groups as well as from a fast dynamical rearrangement of hydrogen bonded networks. Both aspects facilitate a static and/or dynamic percolation network through which proton transport can proceed. Two molecules are considered as clustered if the distance between their P atoms is below 0.51 nm. The comparison of the RDFs showed that below this P-P distance, the presence of the hydrogen bonding can be safely assumed. During the last 1 ns of the simulation, the coordinates of all atoms are recorded every 1 ps. From these coordinates the cluster sizes have been analyzed. Histograms of the cluster size versus the number of occurrences of the clusters are shown for 450 K in Figure 10 (in dotted line). To highlight the effect of hydrogen bonding in the formation of these clusters, the simulation has been repeated with all partial charges set to zero in NVT ensemble to keep the density fixed, while all other force field
Molecular Dynamics Simulations of Heptyl Phosphonic Acid
J. Phys. Chem. B, Vol. 112, No. 25, 2008 7409
parameters have been kept. This removes all hydrogen-bond interactions from the force field. The same cluster-size histogram at 450 K has been plotted in Figure 10 (solid line) from this noncharged simulation. The simulation with partial charges and hence with hydrogen bonds produces larger clusters than the noncharged simulation. In the noncharged model, the phosphonic acid groups have no special attraction and can be assumed to be more randomly distributed that in the model with charges. It is obvious that hydrogen bonding introduces a structuring into the liquid which leads to larger clusters of protogenic groups. In Figure 11, a randomly chosen big cluster snapshot is shown as an example. The phosphonic acid groups are bonded to each other by hydrogen bonds and they are forming channels. The clusters have random shapes and they do not show an overall pattern in general. However, it is possible to see pronounced channels formed by the phosphonic acid groups in all clusters. For larger clusters the channels become more apparent. As the clusters quickly rearrange (below), this leads to a dynamic percolating network of phosphonic acid functional groups inside the material. During the evolution of the system, clusters are breaking up and new ones are forming. The time required for this is on the order of picoseconds. Dynamically, the clusters are not stable and their molecules are exchanged continuously. Consequently, they are not long lasting structures but transient. The fast dynamics is evident in visualizations of clusters which are taken at 1-ps intervals (Figure 12). The number of cluster members lost and gained is indicated in the accompanying table (Table 2). It shows that any one phosphonic acid can find itself as member of several different clusters over a short time span. It can, thus, act as a carrier of a proton from one cluster to the next. This indicates that, while there is no infinite percolating cluster in the network at any one time, protons are effectively transported between smaller transient clusters. This would imply a combined transport mechanism of (i) structural diffusion (proton hopping) within an existing cluster plus (ii) transfer between different cluster by the continuous breaking and reformation of clusters. The simulations show that the heptyl phosphonic acid at the investigated temperatures >100 °C has a dynamics fast enough to support such a transport mechanism.
with their alkane tails. The cluster shapes are often found to be elongated with a phosphonic-acid channel at the center. The dynamics of the clusters is fast at the high temperatures investigated. Clusters break up and reform quickly, so that a given molecule can be part of several clusters over a short time span. Through this mechanism, protons may be transferred between different clusters, so that charge migrates through the liquid efficiently, although there is no infinite cluster with a percolating channel at any one time. This mechanism has possible implications also for phosphonic-acid-containing polymer phases. If the polymer’s architecture allows quick opening and closing of hydrogen bonds between phosphonic acid groups, the mechanism implied in this paper for the low molecular weight liquid heptyl phosphonic acid could still prevail, at least to some extent. This is the subject of ongoing investigations.
Conclusions We have parametrized an atomistic force field for heptyl phosphonic acid, a model compound for phosphonic-acidcarrying polymers, which could act as proton-conducting fuel cell membrane mataerials. The model shows good agreement with available experimental static and dynamic properties of liquid heptyl phosphonic acid at the temperatures relevant for fuel cell applications. It forms a solid basis for force fields for polymers such as vinyl phosphonic acid, which are yet to be developed. In the liquid, we find clustering of the phosphonic acid groups. Visual inspection of cluster geometries reveals similarities to amphiphilic systems. The molecules aggregate via hydrogen bonds among the polar head groups and push their alkane tails to the outside of the clusters, so that neighboring clusters touch
References and Notes (1) Doyle, M.; Choi, S.; Proulx, G. J. Electrochem. Soc. 2000, 147, 34. (2) Kreuer, K.; Paddison, S.; Spohr, E.; Schuster, M. Chem. ReV 2004, 104, 4637. (3) Sekhon, S.; Krishnan, P.; Singh, B.; Yamada, K.; Kim, C. Electrochim. Acta 2006, 52, 1639. (4) Steininger, H.; Schuster, M.; Kreuer, K.; Kaltbeitzel, A.; Bingo¨l, B.; Meyer, W.; Schauff, S.; Brunklaus, G.; Maier, J.; Spiess, H. Phys. Chem. Chem. Phys. 2007, 9, 1764. (5) Schuster, M.; Rager, T.; Noda, A.; Kreuer, k. D.; Maier, J. Fuel Cells 2005, 5, 355. (6) Joswig, J.-O.; Hazebroucq, S.; Seifert, G. THEOCHEM 2007, 816, 119–123. (7) Paddison, S.; Kreuer, K.; Maier, J. Phys. Chem. Chem. Phys. 2006, 8, 4530. (8) Mu¨ller-Plathe, F. Comput. Phys. Commun. 1993, 78, 77. (9) van Gunsteren, W. F.; Billeter, S. R.; Eising, A. A.; Hu¨nenberger, P. H.; Kru¨ger, P.; Mark, A. E.; Scott, W. R. P.; Tironi, I. G. Biomolecular Simulation: The GROMOS96 Manual and User Guide Zurich, Groingen, 1996. (10) Milano, G.; Mu¨ller-Plathe, F. J. Phys. Chem. B 2004, 108, 7415. (11) Mu¨ller-Plathe, F. Mol. Simul. 1996, 18, 133. (12) 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 B. 05; Gaussian, Inc.: Wallingford, CT, 2003. (13) Breneman, C. M.; Wiberg, K. B. J. Comput. Chem. 1990, 11, 361. (14) Ryckaert, J.; Ciccotti, G.; Berendsen, H. J. Comput. Phys. 1977, 23, 327. (15) Faller, R.; Schmitz, H.; Biermann, O.; Mu¨ller-Plathe, F. J. Comput. Chem. 1999, 20, 1009. (16) Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; Vetterling, W. T. Numerical Recipes in FORTRAN 77. Press Syndycate of the University of Cambridge, 1992 (17) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684.
JP0757107
