Research: Science and Education
Modeling the Hydrogen Bond within Molecular Dynamics Peter Lykos Chemistry Division, Illinois Institute of Technology, Chicago, IL 60616;
[email protected] The kinetic molecular theory of gases is the traditional first application of classical mechanics in general chemistry. The model assumes a dilute gas to be made up of a large number of particles moving randomly in a container with no interparticle forces operating. Furthermore all collisions of the particles with the container walls are assumed to be elastic. The result of that analysis, based on Newton’s laws of motion, is compared with the equation of state of an ideal gas, that is, with the empirical Charles and Boyle’s laws to which is added Avogadro’s hypothesis. Via that comparison we find the average kinetic energy of many particles comprising a gas at equilibrium to be proportional to the absolute temperature of that gas. Thus we do introduce the concept of temperature into a mechanical many-particle model of matter. Real gases do not follow the equation of state for the perfect gas well at low temperatures and high pressures. Indeed at sufficiently low temperatures, on compression, real gases undergo a change in state to a liquid. Further the resulting liquid is essentially incompressible. A Nobel Prize was awarded to van der Waals for extending the equation of state for an ideal gas. He recognized that there must be forces of attraction operating between particles for condensation to take place and also that the particles must have finite extension in space or, put another way, experience strong repulsive forces when in close proximity. Thus an important refinement was introduced using bulk variables but reasoning from a molecular model of gases. In the process, the van der Waals equation of state introduced two parameters that, unlike the ideal gas model, particularized it to the specific gas of interest. The advent of the electronic digital computer made it feasible to extend the kinetic molecular theory of gases. Numerical solution of Newton’s laws applied to conservative systems made up of many particles, including provision for intermolecular forces, led to the development of molecular dynamics. Initially that was done for spherically symmetrical rare gas atoms such as argon. But once that method was well developed, Rahman and Stillinger (R-S) extended it to their landmark and comprehensive treatment of liquid water (1). Water’s physical properties1 required the introduction of the concept of the hydrogen bond (H-bond) to rationalize the trends observed for similar compounds made up of elements in the periodic chart in the column beginning with oxygen (similarly for nitrogen and fluorine). Thus, even for a fairly simple model of water, provision had to be made for the H-bond, a bond with a dissociation energy about one twentieth the dissociation energy of ordinary covalent bonds. The H-bond of course is essential to us as living beings. The “rungs” in the twisted ladder model of DNA are just representations of the hydrogen bonds that hold the two strands together. However a simple yet revealing and easily accessible quantitative mechanical model of the ubiquitous hydrogen bond is not customarily developed in any of the usual underwww.JCE.DivCHED.org
•
graduate courses in biology, chemistry, or physics. Thus we took the initiative to extract from the R-S liquid water model how they accommodated hydrogen bonding in the context of molecular dynamics and to incorporate that treatment in undergraduate physical chemistry together with the entire subject of molecular mechanics applied to liquid water. In that manner we hoped to reveal the anatomy of a hydrogen bond. In class we review the general chemistry treatment of the kinetic molecular theory of gases, its features listed in a column, and then go on to systematically refine that model, in a parallel column, until we arrive at the R-S model for liquid water. We then append their treatment of the water dimer as a separate example of their two-body interaction. As an extension of the R-S article, we assign as a problem for the class calculation of the two body potential for the vibration of the two rigid water molecules comprising the water dimer model using the given R-S potential. The Rahman and Stillinger Molecular Dynamics Model of Liquid Water The R-S model for liquid water, limited in scope by the size and speed of the IBM 360兾75 available at Argonne National Laboratory in the late 1960s, had these specifications: The container selected is a mathematical cube containing 216 (6 × 6 × 6) molecules with a volume corresponding to the actual density of liquid water at room temperature. Each of the faces, sides, and corners of that cube are shared with other similar cubes forming a stack of 27 cubes (3 × 3 × 3 or three cubes on each edge) such that periodic boundary conditions could be imposed. What that means is that as a molecule approaches one of the walls in the central cube, instead of colliding with the wall and bouncing back it is allowed to pass through the wall and its mirror image molecule enters the central cube through the opposite wall. That feature was introduced to offset the very large ratio of surface to volume that would have dominated the trajectories of the molecules. The geometry of the water molecule, as well as the isotopes of hydrogen and oxygen, need to be specified so that the principal components of the molecular moment of inertia could be determined (needed for the rotational kinetic energy expression) as well as computation of the intermolecular potential as a function of intermolecular distance and relative orientation. As the model is approximate, R-S assumed that it was sufficient to take the water molecule as a regular tetrahedron with the oxygen at the center, a proton at each of two vertices, and a lone pair of electrons at each of the remaining two vertices. Net partial charges of equal magnitude (0.19 e), with appropriate sign at each vertex are selected, presumably to replicate the permanent electric dipole moment of the water molecule. Of course as an additional exercise this may be checked by actual computation of the electric dipole moment from the model.
Vol. 81 No. 1 January 2004
•
Journal of Chemical Education
147
Research: Science and Education
Twelve parameters need to be specified for each molecule at the time the simulation is started: three Cartesian coordinates and their time derivatives for linear velocity; and three Eulerian angles and their time derivatives for angular velocity. The molecules are assumed to be rigid rotors. The potential is assumed to be a sum over all unique water molecule pairs of an effective two-body potential. The two body potential itself is a sum of two contributions from the interacting pair of molecules as follows: (i)
(ii)
Lennard-Jones (L-J) potential taken to be the neon–neon L-J potential (1) as Ne is the isoelectronic united atom limit for water (as are each in the series methane, ammonia, water, hydrogen fluoride, and of course neon itself ). The L-J potential accounts primarily for the repulsion at small intermolecular distances (hence the incompressibility of liquid water). The electrostatic potential arising from the Coulomb attraction and repulsion between the four partial charges on one rigid molecule and the four partial charges on the other molecule for a total of 16 attractions and repulsions (1). This component would be a point dipole– point dipole interaction at large distances but not at the intermolecular distances encountered here. R-S introduced as a precaution a multiplicative switch factor, changing smoothly from zero at small distances to one at larger distances. That was done because use of a finite increment of time to numerically solve the second-order ordinary differential equations might project a next step to nearly superpose two point charges and thus lead to an enormous nonphysical Coulombic interaction. More specifically, the switching function is zero from zero distance to 2.0379 Å and becomes one at R = 3.1877 Å.2 Between those two values of R the function is S-shaped rising from zero (initially with zero slope) to one, again with zero slope at that point. The students were asked to verify the switch function properties of being continuous as well as having a continuous first derivative needed for a conservative system.
The Hydrogen Bond According to R-S The minimum energy configuration for the two rigid water molecule hydrogen-bonded dimer is displayed in Figure 1. The dimer has a distance of 0.76 Å between the hydrogen of one water molecule and the lone pair of the other of the water molecules (1). As is displayed in Figure 1 each of the four vertices in a tetrahedral water molecule is 1 Å from the oxygen and each bears a charge of appropriate sign
ⴚ
H ⴙ ⴙ H
O 1.00 Å
0.76 Å
ⴚ
O 1.00 Å
ⴚⴚ
ⴙ ⴙ H H
Figure 1. The water dimer at its equilibrium configuration according to the R-S model.
148
Journal of Chemical Education
•
of +0.19 e at each proton and ᎑0.19 e at each lone pair. Thus the distance between the two oxygen atoms is 1.00 + 0.76 + 1.00 Å. There are three equivalent planes of symmetry in the water dimer model, each of which contain the pair of oxygen atoms, and one hydrogen atom and one lone pair between them, all in a straight line. For example, the plane of symmetry in the plane of the paper has a hydrogen atom above the first oxygen that is contributing a hydrogen atom to the H-bond and a lone pair above the second oxygen that is contributing a lone pair to the H-bond. Where physical chemistry has calculus-based physics as a prerequisite it is reasonable to pose the following problem to physical chemistry students: For the water dimer described in Figure 1: (a) Find the four unique intercharge (intermolecular) distances at 2.76 Å. (b) Increase (stretch) and decrease (compress) the Hbond by 0.1 Å and find the corresponding unique intercharge distances between the two sets of four charges. (c) Find the values for the total Coulomb energy at each of the three intermolecular distances (0.75, 0.76, 0.77 Å). (d) Find the values of the switching function at the three intermolecular distances. (e) Find the value of the L-J energy at each of the three distances. (f )
Find the potential energy of the dimer at the three distances. Note the increase in potential energy when the H-bond is both stretched and compressed from the equilibrium geometry.
(g) Assuming the potential is a quadratic in the displacement find the force constant for vibration of two rigid water molecules as a pseudo diatomic molecule. (h)
Check the value of the equilibrium H-bond distance.
(i)
To more fully explicate the several contributions to the total potential energy as a function of separation, examine separately the effects each of the parts.
The potential energy at the equilibrium configuration gives 6.5 kcal兾mole-pair for the dissociation energy of the dimer. Of course that result, based on use of classical mechanics, does not include the zero-point energy correction one would get using quantum mechanics. A more recent and more sophisticated determination of the water pair potential energy per se, estimated from spectroscopic data, gives 5.40 kcal兾mole-pair and 1.74 kcal for the zero-point correction for a net dissociation energy of 3.66 kcal兾mole-pair (2). Thus the relatively crude model used by R-S captures the dissociation energy of the water dimer (and hence the energy of that H-bond) reasonably well. The question remains as to why the H-bond is linear in this model? There are two aspects to that question; namely, rotation about the linear H-bond and bending of the H-bond in any of the planes of symmetry. In both aspects the L-J potential energy is independent of both rotation and bend-
Vol. 81 No. 1 January 2004
•
www.JCE.DivCHED.org
Research: Science and Education
ing. Also note that the R-S model of the dimer is isomorphic with ethane. In ethane the set of three protons at one carbon is staggered with respect to the other set of three protons at the other carbon to minimize the Coulomb repulsion of the positive charges. In the case of the R-S model, the point charges are eclipsed as the alignment of positive charges at the protons (subtended by one oxygen) with the negative charges due to the lone pairs (subtended by the second oxygen) maximizes the Coulomb attraction. The electrostatic interaction is such that rotation around the H-bond increases the distance between opposite charges and decreases the distance between like charges so the Coulombic energy increases. Similarly bending in one of the symmetry planes decreases the distance between one pair of opposite charges but increases the distance between two pairs of opposite charges leading to a net increase in Coulomb energy. The switching function acts in such a manner as to exaggerate slightly those trends. Thus the linear configuration, with eclipsed opposite point charges at the given equilibrium inter-oxygen distance, is the lowest energy configuration. For those people interested in an introductory treatment of the solution of the differential equations in the context of molecular dynamics, a good resource is the Saiz and Tarazona article, Molecular Dynamics and the Water Molecule (3). Conclusion In chemistry, experiment and computer-based modeling are progressing hand in hand to elucidate the structure and other properties of molecules, as well as the properties of bulk matter modeled as assemblages of molecules near equilibrium or even significantly away from equilibrium. Where chemical bonds are made or broken, or where other quantum effects are operating, one needs to apply quantum mechanics. However, where the molecules remain intact, often classical mechanics is sufficient and much simpler to apply to both complicated individual molecules with many internal degrees of freedom as well as to assemblages of molecules such as in liquid water. In this brief introduction to the seminal Rahman and Stillinger article, only the model is explored in the context of using it as a vehicle to introduce undergraduates in physical chemistry to the power and ease of application of classical mechanics to bulk matter from a molecular perspective. Rahman and Stillinger, however, went on to analyze their results from their extensive and carefully done simulation. They took great care to explicate thoroughly and exhaustively, not only the features of the model itself, but also the range of liquid water physical properties that could be assessed with good fidelity. Subsequently they published three additional articles examining the effect on their model of varying the pressure, varying the temperature, and varying the time increment used in solution of the ordinary second-order differential equations (4–6). Thus the article used as the basis of this work, together with the immediately following three sequels, provide a pow-
www.JCE.DivCHED.org
•
erful tutorial framework within which students can be shown the power of classical physics when joined with the modern digital computer and applied to the familiar and, even today, challenging sample of matter—liquid water. Today derivative versions of these methods are being used widely, especially in molecular biochemistry. A recent example appeared in Science (7). A nearly complete treatment of hydrogen bonds appears in Jeffrey’s remarkable book on hydrogen bonding (8). A very recent article (9), “Ultrafast HydrogenBond Dynamics in the Infrared Spectroscopy of Water”, describes an analysis of femtosecond infrared spectroscopic studies of the O⫺H-stretching frequency of HOD in D2O that reveals rearrangements of the hydrogen-bonded network in water. That distinguishes intra hydrogen bond changes that happen over shorter times (170 femtoseconds) from collective structural changes that happen over longer times (1.2 picoseconds). Finally we should note the vibrational potential of the O⫺H stretch shows extreme anharmonicity arising from the H-bond interaction in O⫺H…H as reported in Science (10). Note 1. The Web site http://www.sbu.ac.uk/water (accessed Aug, 2003), provided and maintained by Martin Chaplin, contains a rich and varied display of data and information about water from many perspectives. The site includes experimental information, critiques of various models, and an extensive (over 500 references) bibliography. 2. As the 1971 Rahman and Stillinger article may not be conveniently available to some readers, a referee has urged us to write down the explicit expression for the “switching function” as it is the only formula that is not common knowledge. The expression is, S(rij) = 0
for 0 < rij < RL
S(rij) = (rij – RL)2(3RU – RL – 2rij)兾(RU – RL)3 for (RL < rij < RU) S(rij) = 1
for RU < rij
where RL = 2.0379 Å and RU = 3.1877 Å. One may ask students to verify that the switching function varies continuously and differentially between zero and one as rij grows from zero to infinity.
Literature Cited 1. Rahman, A.; Stillinger, F. J. Chem. Phys. 1971, 55, 3336. 2. Fellers, R.; Leforestier, S. C.; Braly, L. B.; Brown, M. G.; Saykally, R. J. Science 1999, 284, 945. 3. Saiz, E.; Tarazona, M. P. J. Chem. Educ. 1997, 74, 1350. 4. Stillinger, F.; Rahman, A. J. Chem. Phys. 1972, 57, 1281. 5. Stillinger, F.; Rahman, A. J. Chem. Phys. 1974, 60, 1545. 6. Stillinger, F.; Rahman, A. J. Chem. Phys. 1974, 61, 4973. 7. de Groot, B. L.; Grubmueller, H. Science 2001, 294, 2353. 8. Jeffrey, G. A. An Introduction to Hydrogen Bonding; Oxford University Press: Oxford, United Kingdom, 1997. 9. Fecko, C. J.; Eaves, J. D.; Loparo, J. J.; Tokmakoff, A.; Geissler, P. L. Science 2003, 301, 1698. 10. Bakker, H. J.; Nienhuys, H.-K. Science 2002, 297, 587.
Vol. 81 No. 1 January 2004
•
Journal of Chemical Education
149
