I. The quantum mechanical description of atoms & molecules

Newton's mechanics provide an adequate description of large bodies, but small, rapidly moving particles such as electrons require a quantum mechanical treatment. The need for the quantum mechanics was revealed by development of specroscopic methods in the late 19^th and early 20^th centuries. In biophysics, our need for this description of matter also arises mainly in the description of macromolecules provided by spectroscopic experiments. Exact quantum mechanical descriptions of macromolecules such as proteins require a great deal of mathematical sophistication. However, as biophysicists, we usually need understand only the basic principles of the quantum view of matter in order to understand the nature of the chemical bond and to interpret spectroscopic observations on macromolecues. Here we review these principles and show how quantum mechanical treatments of four fairly simple systems allow interpretation of the properties of even complex macromolecules.

A. Fundamental assumptions or postulates

1. electromagnetic radiation has a particle nature; it's energy is quantized

- Max Planck's explanation of black body radiation spectrum and Albert Einstein's explanation of the photoelectric effect

E-C 10-4

2. small particles have a wave-like nature: de Broglie's PhD thesis: l = h/p

- this explained why electrons were diffracted like x-rays or light and accounted for the quantization of angular momentum assumed by Niels Bohr's view of the hydrogen atom as having "orbitals" of quantized energy (E=-e²/2a₀n²; a₀=h²/4p²me² = Bohr radius = .053nm).

E-C 10-6

- a corrolary of this postulate is Heisenberg's proposal that the exact position and momentum of a particle could not be specified simultaneously:

DpDr » DEDt » h

- another corrolary is that the Schrödinger equation can describe the behavior of particles interms of wave mechanics:

The term in brackets is the Hamiltonian (or energy) operator (H), E is its' expectation value (the energy) obtained by operating on the wave function y. This leads to the operator form of the Schrödinger equation: Hy = Ey. The first term of H is the kinetic energy operator (p²/2m) while the second term is the potential energy operator. In general, the Schrödinger equation has an infinate set of solutions, y_i, called eigen functions of the Hamiltonian operator, each one with its corresponding expectation value or eigen value of the Hamiltonian.

B. Some Key Simple Models: For the most part, the motions of even complex macromolecules (N atoms) are described in terms of three simple types of motion: translation in three orthogonal directions of the molecular center of mass, rotation of the three principle molecular moments of inertia, and 3N-6 internal vibrations of atoms determined by covalent bond strengths, molecular geometry, and atomic masses. Here, we give the quantum mechanical description of the three simple models describing these motions.

Particle in a box: d²y/dx² {-h²/(8p²m)} = Ey

y_n = (2/a)^½ sin(npx/a)

E_n = n²h²/8ma², n=1,2,3....

E-C 10-10

E-C 10-11

Harmonic oscillator or the weight on a spring. This is one of the most important physical models in biophysics. We will examine this model and its' classical description in detail in the Spectroscopy module. For now, we simply state the quantum mechanical description. The energy levels are:

E = (n+½)hu, with u = (k/m)^½; note the ground state energy, ½hu.

E-C10-12

E-C10-13 E-CT10-1

Rigid rotator. The rotation of any polyatomic molecule can be treated in terms of rotations about the three principle axies of an unsymmetric top (three unequal principle moments of inertia: I₁, I₂, I₃). The largest possible moment of inertia of a molecule is I₁, the smallest is I₃, I₂ is the moment about the axis perpendicular to the axies of the other two principle moments. The point of intersection turns out to be the center of mass. Unfortunately, the general quantum mechanical treatment of rotation of an unsymmetrical top is quite complex. Even the treatment of a symmetric top (two equal moments) is complex (e.g., see asymmetric diatomic rigid rotator in E + C). Fortunately, the important features of the quantum mechanical treatment of rotation are contained in the treatment of rotation of a mass in a plane about an axis perpendicular to the plane (Atkins, 408-411). This rotator moves in one circular dimension sweeping out an angle f, i.e., it embodies the simplest type of angular motion. Its' quantum mechanical description is:

d²y/dx² = -2EI/h² y

y_m = e^imf/(2p)^½

E = m² · h²/2I, m = ±0, ±1, ±2, ....

For convenience, the moment of inertia is often expressed in temperature or frequency units:

Q_r = h²/2Ik, where k is the Boltzman constant. E = m² · kQ_r. Note that Q_r has units of ^°K.

B_r = h/4pI has units of frequency. E = m² · hB_r

Two observations about the properties of the single dimensional quantum mechanical rotator are worth noting:

1. The energy levels of the rotator are degenerate, there being two wave functions (corresponding to clockwise and counterclockwise rotation) for each energy level. This degeneracy becomes much more complex in the case of the symmetric and asymmetric tops, but the principle is the same. It will be important to remember this degeneracy when we consider the partition function of the simple, one-dimensional rotator.

2. For most moderately large and heavy molecules (certainly for all macromolecules but probably not for H₂), Q_r <<< T (the experimental temperature). Thus, molecular rotations are nearly always in the classical limit (DE <<< kT), meaning that molecules are distributed according to the Boltzmann distribution over their molecular energy levels. Thus, a complex quantum mechanical treatment is not necessary when calculating the rotational contribution to the molecular partition function (see Section IV).

C. Hydrogen atomic orbitals. This is the most important simple quantum mechanical system in chemistry. It is the basis of most interpretations of IR, visible, or UV spectroscopic experiments and provides the framework for our understanding of atomic structure and chemical bonding. The simple spherical symmetry of this problem makes possible an analytical solution, the wave functions of which are outlined in Table 10-3 from E & C.

The energy of the system is determined by the principle quantum number, n:

E-C 10-17; T10-3

E_n=-{2(pZe)²m}/(nh)², where

m=the reduced mass » m_e.

The energy levels are degenerate, with the degenerate states being described by the azymuthal (l, angular momentum), the magnetic (m, direction of ang. mom.), and the spin (m_s) quantum numbers. The spin quantum number does not result from the simple treatment in E & C but is due to a relativistic treatment of Dirac. The degeneracy due to m and m_s is broken in

a magnetic field.

The azymuthal quantum numbers are usually given th familiar designations s, p, d, f.

The familiar shapes of the first few H-atom wave functions are given in E & C Fig 10-18. The squares of these functions (which look much like these functions) are the probability distibution functions of electrons around the nucleus. By comparison to the Bohr description, these are called orbitals, although this is a poor nomenclature.

10-18 & EC T10-3

The H atom is the only atom described exactly by quantum theory! The theory of atomic structure approximates the effect of filled inner orbitals by shielding functions and then uses H-atom wave functions to describe the behavior of outer or valence electrons. This results in the familiar

rules for the order of filling of valence orbitals, which provides the organizing principle for the periodic table of the elements. Recall that completed orbitals always have paired electrons. The Pauli exclusion principle states that in allowed atomic states no more than one electron can have the same set of four quantum numbers. Thus, paired electrons must have m_s = +1 and -1.

D. H₂⁺ molecule; LCAO MO's. The H-atom wavefunctions also serve as the basis for the quantum mechanical description of chemical bonding, i.e., the description of molecules.

The Hamitonian operator belongs to a class of operators called Hermitean operators. Hermitean operators have several useful properties, including that their eigenfunctions are orthogonal and form a complete set of basis functions (see Bohm for a complete discussion). Orthogonal means:

where "*" designates the complex conjugate of a wavefunction. "Completeness" means that any other wavefunction, Y, can be represented as a sum of contributions from the members of an infinite complete set of wavefunctions of the Hamiltonian of the H-atom:

This means that any molecular orbital (MO, wavefunction for a collection of atoms) can be represented as a linear combination of atomic orbitals (LCAO) of the H-atom. The appropriate linear combination is found using the variation principle, which states that the expectation value of the Hamiltonian (energy) of the true molecular ground state wavefunction £ the energy of any other approximate molecular wavefunction. The energy (for that matter, the expectation value of any operator) is written as:

In practice, then, the j^th molecular orbital wave function (Y_j) is defined by finding, using linear algebra methods, the coefficients, a_ij, that minimize the energy of the molecular orbital ground state. E & C discuss in detail the rules for forming MO's and give examples in terms of diatomic molecules, which can form bonding or anti-bonding MO's. (This chapter might be useful to read to remind yourself of the origen of all the material you learned in freshman chemistry, although we will not be concerned with these details here.) For any molecular geometry, in principle, a set of MO's and E's can be found in this way, and electrons can be assigned to these by procedures similar to those used for assigning electrons to AO's. In most cases, unless there is substantial delocalization of electrons between several atoms (e,g,, aromatics or H-bonded systems of water molecules), it is not necessary to carry out such extensive calculations, and chemists often content themselves with considering the distibution of electrons between pairs of atoms. Before this is done, chemists have often found it necessary to take appropriate linear combinations of atomic orbitals to form hybrid atomic orbitals that can account for the observed geometries of covalent bonding patterns (e.g., sp₃ orbitals of oxygen account for the tetrahedral distribution of electrons around oxygen in water).

E & C 11-8

From the viewpoint of biophysics, the LCAO-MO desciption will be useful mainly when we come to discussing the spectroscopy of macromolecules. For now, it is worth noting how covalent bonds between atoms and the electric field of a molecule can be described in terms of probability distributions of electrons. Covalent bonds reflect the high probabilty of electrons being located between ajacent atomic nucleii in bonding orbitals, resulting in lower potential energy, as well as the lower kinetic energy associated with electrons being distributed over larger distances as compared to their atomic or anti-bonding orbitals (see E & C figs 11-2 and 11-1).

E&C 11-2; 11-1

Another molecular property that must be understood from the point of view of the quantum mechanical description is the molecular charge distribution. This is needed to describe the interaction of molecules with other molecules (see non-covalent forces below) or with electromagnetic radiation, as we will discuss in the module on spectroscopy. As we will see below, the simplest description of a general charge distribution is in terms of a dipole moment (E&C fig 11-13).

Quantum mechanically, the dipole moment (m) is given as the expectation value of the dipole moment opperator:

, where r_i is the distance of the i^th electron from the charge center of the molecule. Anything that alters the quantum state of the molecule will also alter the molecular dipole moment (induce a dipole). This is the basis of the description of UV, visible, and IR molecular spectroscopy to be given in Biophys 146.

II. Covalent forces between atoms

A. Born-Oppenheimer Approximation In treating the H₂⁺ molecule, it is always assumed that the motion of the nucleii can be ignored when treating the motion of electrons, since nucleii are roughly 2000 times more massive than electrons. This has the effect of ignoring coupling between nuclear and electronic motions. As a corrolary of this approximation, the potential energy of a diatomic covalent bond can be obtained by calculating the MO energy as a function of interatomic distance. This interatomic potential is thus seen a provided by the probabilistic distribution of electrons in bonding MO, with this distribution being able to adjust instantaneously to alterations in interatomic distance. In practice, these potentials are not obtained by MO calculations but derive instead from IR spectroscopy, since the frequncies of bond stretching vibrations provide information on the shape of the potential.

B. The Morse Potential & the anharmonic oscillator Experimentally, it is known that the interatomic potential is anharmonic (not having the same potential energy function as the harmonic oscillator, ½kx²). This is evident from the facts that energy levels are not all equally spaced (IR spectra are temperature dependent) and that covalent bonds do break if sufficiently excited by thermal energy. Of the many empiracle functionalities that have been proposed for this potential, the Morse functionality is most widely used, and is commonly written as:

Solution of Schrödinger's equation for an anharmonic potential of this form yields energy levels:

This functional form has V(0)=0 while it should be a negative number. This is easily accomplished by subtracting D_e from the functional form shown. This also makes V(r)-> 0 as r->¥.

Note that this approaches the harmonic potential as D_e approaches infinity.

III. Non-covalent forces between and within molecules The forces between molecules or atoms are largely electrical in nature, although our description of them only sometimes reveals this.

A. Potentials & forces; electrostatic interactions; point charges. In general, the force of interaction between two molecules will be a function of the distance between them and a complex function of their orientation relative to each other. Their potential energy of interaction is the integral of the force of interaction over the distance or relative orientation variables that describe the force. If we ignore for now the orientation problem, then

Note that the reference zero for potential energy of interaction is taken at r = ¥. For orientation-dependent potentials, there is an additional torque tending to rotate the molecules. This problem becomes much more complex, but can be treated by Lagrange's equations (see Chapter 1.4 of Hirschfelder et al.) to obtain the complex motion of two such approaching bodies.

Polar molecules can often be approximated as a collection of point charges. In this case, Coulomb's law gives the force between two point charges and the force per unit charge on q₁ (E₁ = field at charge 1 due to charge 2) as

, where the proportionality constant = k = 1/(4pe₀), e₀ being the permitivity of a vacuum (8.8·10^-12 s²A²kg^-1m^-3). Note that the force and the field are related vector quantities, where e₁₂ is the unit vector along r from 2 to 1. Note that the Coulomb force between point charges is strictly distance and not orientation dependent. However, the force between two collections of point charges will depend on all inter-charge distances and unit vectors, as such, will depend on the orientation between these fixed charge collections.

By contrast to the force and electrostatic field, the potential energy of charge 1 in the field of charge 2 (or, visa versa of 2 in the field of 1) is a scalar quantity defined as the work needed to move charge 1 from infinity to a distance r form charge 2.

The electrostatic potential (f) is defined as the potential energy that a unit charge would have if brought from infinity to some point (in our case a distance r from charge 2).

B. A collection of charges, the dipole approximation. The electrostatic potential (f) due to any distribution of charges at a distance R from some point in the distribution (see Fig 6-7 of Feynman-II) is given by:

The sum over 1/r_i is a complex function of the geometry and orientation of the charge distribution of the molecule under consideration. Since r_i » R - d_i·e_r (e_r = unit vector from origin to P along R), and d_i << R, the quantity 1/r_i can be expanded as a Taylor series about 1/R in powers of d_i/R. This gives the expression:

If we define the dipole moment as m = S_iq_id_i, then the field due to the whole molecule is approximated by:

The ignored higher order terms (called quadrapole, octopole, etc.) become important at very small distances or when the net charge on the molecule (Q) is zero and the dipole moment is zero or very small. Note that the dipole moment is a vector quantity, i.e., it has both a direction in the molecule and a magnitude. Thus, the potential between two ideal point dipoles is a complex function of orientation (see chapter 12.1 of Hirschfelder et al. to obtain the expression given below).

If this complex function is averaged over all orientations in space, an average, temperature-dependent dipole-dipole potential is obtained (see chapter 13.5 of Hirschfelder et al.).

These expressions are for ideal point dipoles; even more complex expressions hold for real dipoles (charges separated by some distance l) especially as r_ab becomes comparable to l. Obviously, the dipole approximation is useful only when molecules are reasonably well separated.

C. Dispersion & repulsion; Lennard-Jones & other potentials; H-bonds. Even if a molecule does not have a net charge or a dipole moment (i.e., a non-polar molecule), it can interact with other charged, polar, or non-polar molecules through induced dipoles. The dipole induced in a molecule that experiences an external field (E) is m = a E, where the tensor a is the molecular polarizability. The origen of polarizabiltiy lies in electron redistribution in response to a field, as

we will describe in detail in the module on spectroscopy. As two electron distributions approach each other, they perturb each others' wavefunctions causing redistributions of electrons in each molecule. The resulting induced dipole-induced dipole interaction can be described and calculated classically or quantum mechanically, and is thereby shown to vary as 1/r⁶. However, the magnitude of the potential of interaction (called the London Dispersion Potential) is most often estimated experimentally from transport coefficients and second virial coefficients of dilute gases. This interaction is usually assumed to be independent of relative molecular orientation, but, of course, need not be, since the polarizability

tensor has a distinct orientation in a molecule.

Atoms and molecules also experience strong repulsive interactions as they approach to the point that their electronic wave functions begin to intermingle. When this occurs, MO's can form. If unoccupied valence orbitals are available, the occupancy of bonding MO's will exceed that of anti-bonding MO's and a covalent bond will form. If not (e.g., two He atoms), the two AO's will overlap, resulting in a probability that two electrons in the same orbital have the same spin quantum number, which violates the Pauli exclusion principle. The result is a very sharp (roughly exponential) and strong repulsive potential (see E & C fig 11-17). This has been approximated in several ways. One is in terms of hard spheres with radii given by the van der Waals radii compiled by Linus Pauling (E & C table 11-5) obtained from crystal data and gas phase virial coefficients. The second is to assume a 1/rⁿ functionality. With n=12, this potential combined with that for the London Dispersion interaction results in the familiar Lennard-Jones 6-12 potential:

, where s is the r value for which V=0 (approximates the van der Waals radius) and e is the depth of the potential at its minimum. This is one of several potentials that have been suggested to account for the interactions between non-polar molecules (others include the Sutherland and Buckinham potentials). The combination of repulsive and non-polar attractive interactions between molecules is often referred to as van der Waals interactions.

For polar molecules, this is clearly inadequate. The Stockmayer potential (Hirschfelder et al., p23) adds to a Lennard-Jones potential a term for orientation-dependent interactions of point dipoles; real dipoles would require an even more complex additional term. Dipole-induced dipole interactions are included in the 1/r⁶ term, as are orientation-averaged dipole-dipole interactions, which make the 1/r⁶ term temperature-dependent. Other types of interactions (see Atkins table 22.3) can be included by adding additional terms in 1/rⁿ.

A significant complication occurs when H atoms occur between two electo-negative atoms (A××H-B). When the A××H distance can be shown to be < R_A + R_H - 0.2A (R = van der Waals radius), a hydrogen-bond is said to exist. In some cases, H may be shared equally between A and B (e.g., ice) but usually one interaction is greater. The energy of this interaction (~5Kcal/mol) is intermediate between that of covalent (~100-200Kcal/mol) or ionic (~50Kcal/mol) and that of dispersion or dipolar interactions (~0.5-1Kcal/mol). The nature of the H-bond has been widely argued, but some features seem clear:

1] they involve redistribution of electrons within as well as transfer between A and B.

2] they need not be linear; bent bonds as well as bifurcated (A₁,A₂××H-B) happen.

3] in most cases, A contains lone-pair electrons, which may be responsible for the H-bond.

D. Empirical force fields. It should be evident from the descriptions of intra- and inter-molecular forces given so far that a complete and accurate description of the structure and motions of a macromolecule would be an enormously difficult and complicated undertaking. Bonded interactions and even some non-bonded interactions (e.g., 1-2-3-4 torsional interactions, H-bonds, etc.) require ab initio LCAO-MO calculations for an accurate description, and our discussions have demonstrated that even these involve certain assumptions. We have also seen how even classical treatments of non-bonded interactions between polar molecules in condensed systems and moderate temperatures can become enormously complicted due to orientational effects. So how might one go about predicting the structure and dynamics of a macromolecule? The answer is "only by using a very approximate description of the intra- and inter-molecular potentials designed specifically for the problem at hand." This is called an empiracle potential. The article by Bowen and Allinger describes the art and science of defining and parameterizing such a potential, the MM3 potential and computational method. Several points about such potentials and their use in molecular dynamics calculations can be made:

1] potential parameters are not general, i.e., a potential optimized to predict the properties of one type of molecule or group will not necessarily work for a different molecule or group.

2] intra-molecular vibrations are often treated very approximately as bond stretches, bends, and torsions. Coupling between these motions is sometimes included to match spectroscopic data. All vibrations are treated using approximate harmonic or slightly off-harmanic potentials. While this is computationally convenient, it is theoretically incorrect, since it does not properly keep track of the internal degrees of freedom of the molecule. We will see how to do this correctly during our description of vibrational spectroscopy.

3] non-bonded interactions are treated by a van der Waals interaction (a Lennard-Jones or some similar potential) in combination with point charge or dipole approximations to the electrostatic contribution. It is these terms that will be very system dependent.

E. Structure calculations: energy minimization. The earliest attempts to predict the structure of biological macromolecules used empirical force fields to calculate the potential energy (V{r_i}) of a given arrangement or "configuration" of atoms in a macromolecule and then standard minimization routines to minimize this function, rationalizing that the most likely structure would be the one with the lowest internal energy. Advantages: very low computation time µ N. Disadvantages: ignores the average or dynamic nature of macromolecule configuration; => "local minima".

IV. Microscopic energy distribution => macroscopic thermodynamic properties

A. The Boltzmann distribution and the molecular partition function. Here we examine how we make a link between the microscopic or molecular description of macromolecules and their macroscopic or thermodynamic behavior (E&C fig 14-1). This is done using the tools and concepts of statistical mechanics. The usefulness of this subject in not just in the formalism of the partition function but also and especially in the conceptual basis that it gives to thermodynamics and to the farily new field of molecular modelling.

If we have N independent and distinguishable molecules, each of which is characterized by a quantum state y_i and energy e_i, there are usually many ways in which these molecules can be distributed over the available energy states (n₁, in a state with e₁, n₂ in e₂, etc., describesd as {n_i}). Each is referred to as a configuation of the system of molecules. It is evident (see A fig 19-1) that

some configurations will have many more arrangements that others. For very large N, one distribution is represented in much greater numbers than all the others.

Boltzmann in Austria and Maxwell in England showed at the end of the 19^th century that the most likely configuation could be found by maximizing the function W_{ni} subject to the condition that the total energy and number of molecules in the system remained constant.

The resulting configuration is characterized by a distribution of molecules into states such that the probability that a molecule is in state i is

This is called the Boltzmann distribution. Distributing molecules according to this distribution maximizes the randomness of the system of molecules, since this distribution allows for distinguishable molecules to be distributed over available energy states in the greatest number of ways, subject to the stated constraints. In this sence, the observation that systems of molecules obey the Boltzmann distribution is a statement of the second law of thermodynamics. The probability that a molecule exists in a high energy state increases with temperature (E&C fig 14-4).

This has two physical consequences. First, molecules in very high energy states are highly unlikely but not forbidden. Second, increasing E of a system must increase the temperature of the system so as to make high energy states more probable.

It is convenient to define the molecular partition function as:

where the sum over i is over quantum states while the sum over j is over energy levels of degeneracy w_j. Physically, the molecular partition function is the average number of states available to a molecule at temperature T. As T->0, q-> degeneracy of the ground state level; as T->¥, q-> total number of states, whether this is finite or infinite. This is evident in the temperature dependence of the partition function for a molecule free to distribute over two equally spaced quantum states (A fig 19.3).

B. The configurational energy and entropy of a system of Boltzman molecules. It is now easy to calculate the total configurational energy (due to distribution over available states) as:

Since the Boltzmann distⁿ maximizes the number of energy state arrangements in a particular configuration, it is easy to see why Boltzmann postulated that the configurational entropy of such a system of independent, distinguishable molecules is:

The justification of this fundamental postulate is given in Atkins (p679), and the relation of W to p_i comes from the derivation of the Boltzmann distribution (A667,eq3). Since

, the Helmholtz free energy is:

These expressions are the basis for the relationship between the statistical description of systems of molecules and the macroscopic formalism of thermodynamics.

C. The cannonical or system partition function and free energy. The molecular partition function was defined under conditions such that the system was defined in terms of

the thermodynamic variables N, E, and V. It is much more common and convenient experimentally for system variables to be N, V, and T. To facilitate description of macroscopic systems in these terms, it has been found convenient to define an ensemble of thermodynamic systems, each with the same N, V, and T but with different energies, E_i. This is illustrated in A fig 19.13 for 20 systems; in general, we envision a very large number of systems. By analogy with the description of the molecular partition function, there are many ways in which these systems can be assigned to the energies, E_i. For the dominant distribution, the probability, P_i, that a system will have energy E_i is:

Q is called the cannonical partition function and the ensemble of systems is called a cannonical ensemble. There are two important points that must be made about the cannonical partition function. The first is that, if this partition function can be defined on the basis of some molecular description of the system, then all the thermodynamics of the system can be obtained from it. Table 20.1 (A) gives the expressions that accomplish this.

In these expressions, the reference state for which internal energy = U(0) is taken as the internal energy at T = 0°K. To this must be added the internal energy due to thermal distribution over available energy levels. Recall that the Helmholtz free energy is that portion of the internal energy of a system available to do work

(Arbeit). At the reference state (T = 0°K), A(0) is totally in the form of internal energy, U(0). At temperature = T, A(T) is decreased by the energy devoted to increasing the configurational entropy of the system, energy unavailable to do work. Also recall that a spontaneous change must have DA£0. Thus, for a change in a system to occur spontaneously, it must result in an increase in the system partition function, i.e., in the average number of states that contribute to the average system energy. The entropy, as defined by Boltzmann, is seen to be related to the total internal energy change (rel. to T=0°K) minus that portion avalable for work.

Also given are the expressions for enthalpy (H) and Gibbs free energy (G), which is the internal energy available to do work after accounting for the effect of pressure-volume work. Although A is the more fundamentally important quantity, G is of greater practical significance since most experiments are performed at constant p rather than constant V.

The cannonical partition function can be expressed in terms of the product of molecular partition functions for independent molecules, molecules whose energy level distributions do not depend on the energy level distributions of other molecules in the system. Thus, the partition function for a system of distinguishable molecules, N₁ of type 1, N₂ of type 2, etc., is

Hill ("Statistical Thermodynamics") derives this expression for the partition function for an Einstein crystal by showing that the sum over all the possible states of independent molecules is equal to the product of the sums for individual molecules (see Prob 4). The fact that all molecules reside at particular lattice sites makes the molecules distinguishable. This is exactly the result that we obtained when we derived the Boltzmann distribution. Note, however, that molecules in an ideal solution or in an ideal gas are not distinguishable, so that we count too many configurations by using the above expression. Thus, the partition function for an ideal gas containing a mixture of independent, indistinguishable molecules (N_i of each kind) becomes:

Once the molecular partition function of each type of molecule is known, the system partition function can be written and leads immediately to the system thermodynamic properties. Now, the issue becomes how to write the molecular partition function.

D. Macromolecules: independent translations, rotations, and vibrations. The problem of writing the molecular partition function simplifies greatly if the quantum states associated with different modes of motion of the molecule are independent. Then, the total energy is e = e_elec + e_vib + e_trans + e_rot, and the partition function becomes: q = q_elecq_vibq_transq_rot (again, sum = product of sums). At normal temperatures, you know now that only the ground state electronic level is occupied, so q_elec = w₀e^-e0/kBT, where w₀ is the degeneracy and e₀ is the energy of the ground state. The partition functions for each of the nuclear motions can be obtained from the known quantum mechanical description of the harmonic oscillator, the particle in a box, and the rigid rotator:

Note that this form of q_rot is the form given in texts; my form contains another factor of p as derived in your problem set. Note also that q_rot and q_trans are each written for 3 degrees of freedom, while q_vib is written for 3N-6 independent vibrational degrees of freedom. A linear molecule will have only 2 rotational degrees of freedom and 3N-5 vibratonal ones. These contributions to the molecular partition function are all simple expressions in the classical limit, and it is thus easy to calculate from them the contributions of each molecular degree of motional freedom to the system thermodynamics. In your problem set, you demonstrate that each rotational and translational degree of freedom contributes ½k_BT and each vibrational one k_BT (½k_BT from potential energy and ½k_BT from kinetic energy) to the system internal energy. This is the very important and useful principle of equipartition of thermal energy.

Now, we see the elegance of statistical mechanics; the thermodynamic properties of even the most complex macromolecule can be and are described in terms of very simple models for the various motions of all the atoms of the macromolecule (see prob 5).

E. Fluctuations One of the major insights from the statistical treatment of molecular events is that even unlikely events are possible. Thus, a system with N, V, and T fixed (described by the cannonical partition function) has a probable energy E_av, but deviations from this value do occur. For this system, there is a distribution of energies with a characteristic width given by the root mean square deviation from the mean (s_E, Hill, p35):

Note that fluctuations from the average energy are proportional to the heat capacity (related to the second temperature derivative of the partition function) and increase with the square of the temperature.

Similarly, for an open system described by fixed T, V, and m (chemical potential), we can define a grand cannonical ensemble and a grand cannonical partition function that allows us to calculate not only the most probable energy but also the most probable number of particles, N_av. This partition function and the thermodynamic functions of such a system are described in detail in Hill's book, if you are interested. For the present, we wish only to note that there is in such a system

a similar probability distribution in N with a characteristic width of s_N:

Note that the fluctuations in number density in such an open system increase with temperature but decrease with volume and are proportional to the isothermal compressibility, k. By the proper choice of system and partition function, the fluctuations in any dependent thermodynamic quantity can be described.

F. Structure Calaculations: Monte Carlo Methods One of the currently popular methods for estimating the structure of biological macromolecules from first principles uses the principles of thermodynamics to estimate the probability that a given structure or configuration contributes to the average structure. An empirical force field (Section IIID) is used to calculate the potential energy of a given configuration (V{r_i}), and this n^th configuration is essentially "counted" with a probability proportional to the Boltzmann factor, e^-V/kT. Advantages: computation time µ N·# configs; Boltzmann weighted average structure recognizes macromolecules are not rigid; "local minima" and energy barriers are overcome. Disadvantages: difficult to avoid "unreasonable" configurations with stochastic or random generation procedures; kinetic energy contributions not included.

V. Water: the unique biological environment For the most part, liquid water is the environment in which biological macromolecules must perform their functions. Unfortunately, liquid water is one of the most difficult substances to understand on a molecular level; there have been volumes written on the subject of water and ice! This is first because of the highly polar and tetrahedral structure of the water molecule, second because of water's ability to form H-bonds, and finally because water is a small molecule whose dynamics are difficult to treat.

A. The water molecule & hydrogen bonds The open and tetrahedral structure of the ice I lattice (normal ice - see figure from Pauling) implies that the 2s and 2p orbitals of the water molecule oxygen must hybridize to form the MO description of H-bonded water molecules in ice. This was first recognized by John Pople in 1951 and is illustrated in figure 1.5 from Eisenberg & Kauzmann.

A model for the structure of ice shows H atoms located in discrete positions between oxygens. Actually, Pauling recognized in 1935 that the residual configurational entropy (k_Bln3/2) displayed by ice at 0°K could be explained by H atoms distribution between two possible O atoms. This also explained the very high conductivity of ice. This H redistribution attributes to partial covalent (i.e., electron sharing)and non-additive nature of the interaction between water molecules.

Another feature of the ice structure worth noting is the amount of unoccupied space which it contains due to the tetrahedral nature of the O atom arrangement. Water is nearly unique in the fact that the solid form is less dense than the liquid (ice floats on water). Indeed, the density of water continues to increase with increasing temperature until 4°C, after which it decreases as for any normal liquid.

These features of ice structure help us understand the nature of the interactions between water molecules and the resulting unique structure and solvent properties of liquid water.

B. Structure of liquid water Water has a number of unusual properties, as summarized in the review article by Stillinger (1980). Most remarkable are:

i. d_liq > d_solid.

ii. d_liq has a maximum at 4°C.

iii. anomalously high T_M, T_B, and T_C.

iv. high H⁺ and OH^- mobility in both the solid and fluid states.

v. large "configuational" heat capacity (see fig 4.10 from Eisenberg & Kauzmann).

The "configurational" heat capacity might better be called the "structural" heat capacity. It is that part of the heat capacity that arises from a continuous change in structure of a phase, rather than from excitation of mechanical degrees of freedom (translation, rotation, vibration). Eisenberg & Kauzmann (1969) attribute nearly half of water's unusually large heat capacity to configurational heat capacity (see Figure 4.10). Unlike the heat capacity of normal liquids, the heat capacity of water decreases with temperature due to a decrease in this configurational term (see fig).

The latter three observations suggest that a great deal of the H-bonded structure of ice persists in liquid water, although the first two observations do indicate that this open structure is broken down relative to ice and continues to break down with increasing temperature. Several models of water structure have been proposed to account for these properties. In general, these have fallen into two classes: continuum theories and discrete

cluster theories. The former took several forms, although all involved some continuous, temperature-induced breakdown of the ice lattice. Discrete species models were first suggested by Felix Frank, who proposed that water could be viewed as an equilibrium mixture of small, ice-like clusters and monomers. Nemethy and Scheraga adopted this view in developing their once-famous cluster thermodynamic treatment of water and aqueous solutions. Theories based on both models could account qualitatively for the unusual properties of water and some (e.g., Nemethy & Scheraga) could even account for the unusual thermodynamics of the liquid. None, however, attempted to account for the dynamic properties (spectra, dielectric relaxation, etc.) The arguement between these two camps was largely put to rest by the pioneering molecular dynamics treatment of Stillinger, which we will consider in detail. Stillinger used an empiracle force field that accounted for the non-additive nature of water H-bonds, i.e., the fact that the potential energy of small cyclic H-bonded arrangements is greater than predicted by the sum of pair-wise interactions of water molecules (Lentz & Scheraga, 1973). Stillinger accounted for the structural, dynamic, and thermodynamic properties of liquid water to a reasonable approximation. The picture that emerged was a composite of the continuum and discrete views: small polyhedral clusters of H-bonded water molecules tend to persist on the picosecond time scale, but these have distorted H-bonds and not the perfect tetrahedral structure of ice. At the same time, a class of "non-bonded" (less perfectly H-bonded) water molecules exists (see Fig. 3 from Rahman & Stillinger, 1973). Thus, the

highly polar, roughly tetrahedral nature of the water molecule results in an open, H-bonded structure for liquid water but with a great deal more molecular disorder than found in ice. It is interesting that an improved treatment of water based on the cluster model and also accounting for the non-additive nature of water H-bonds and the dynamic properites of clusters of water molecules simultaneously gave nearly the same picture of liquid water (Lentz et al., 1974).

C. The hydrophobic effect. Solutions of non-polar molecules in water also show unusual properties (see Tables taken from Nemethy & Scheraga).

This figure illustrates a water-hydrocarbon clathrate structure. Shown is the crystal structure for tert.-butyl amine clathrate.

Nearly all of the unfavorable free energy associated with these solutions comes from the entropy rather than the enthalpy term (Table 1). In addition, the partial molar volumes of non-polar molecules are often much lower in water than in non-polar solvents (Table 2). Finally, non-polar compounds are known to form crystalline hyrates whose cystal structures reveal cages of water molecules surrounding the non-polar molecule. These are called clathrate structures. On the basis of these observations, Nemethy & Scheraga, after the initial suggestion of Frank & Wen, proposed that non-polar solutes actually interact favorably with the open structure of liquid water because they are able to fit into the large free volumes in water. This has the effect, however, of making the fully H-bonded water species energetically even more stable, thus enhancing the ice-like arrangements in water. This results in the large negative entropy change associated with hydration of non-polar solutes. Whether or not one accepts the simple cluster model of Nemethy & Scheraga, the fact is that non-polar solutes appear to enhance the structure of already structured liquid water. This is the basis of the hydrophobic interaction. It comes about because water next to a non-polar solute is structured or entropically unfavorable relative to bulk water. There is, therefore, a "thermodynamic potential" driving the association of non-polar solutes so as to decrease the water-solute interface:

This reaction has a large positive DS and a small (positive or negative) DH and is thus entropy driven. Note that this can not be termed a "force" or "potential", since there is no distance dependence, it depends on the comparison of two thermodynamic states. There is recently confusion on this issue, as Helm et al. (1992) report a "hydrophobic force" between bilayers mounted on mica sheets. If Helm et al. are correct, we would need to redefine what is meant by a "hydrophobic interaction." However, Tsao et al. (1993) have disputed this claim and propose that Helm's force was due to dipole-dipole interactions between membrane micro-domains.

VI. Molecular mechanics & introduction to molecular dynamics

A. The molecular dynamics calculation.

System of N atoms (masses = m_i) with 3N Cartesian coordinates expressed as a vector (). The dynamic behavior of this system is calculated by applying Newton's equations of motion:

1] Þ V_pot() Þ the internal force on particle i due to each other particle j and the acceleration of particle i due to each other particle j (a_ij):

2] Assuming constant acceleration over time step dt, we obtain the change in velocity of particle i under the influence of particle j:

dv_ij = a_ij dt

3] Numerical integration of the equations of motion yields a matrix of values dx_ij, the i^th row of which gives the incremental change in the postion of atom i due to all the other atoms in the system:

4] Summing the elements of each row yields a new set of coordinates, , and these yield a new V().

Now, repeat steps 1 - 4. By saving the _n after each step, we have a complete record of the positions of all the particles in the system at t = ndt.

This calculation is carried out on a relatively small group of molecules in a small, properly scaled volume and then this volume is repeated such that the left boundary of the cell to the right of the central cell is identical to the left boundary of the central cell (see figure). This is called periodic boundary conditions. To further limit the time required for calculations, certain atomic pairs beyond a cutoff sphere are ignored.

We have made two implicit assumptions in this procedure that deserve mention. First, we have assumed that all molecules in our system behave according to classical mechanics. This is reasonable for most proteins but might not be reasonable for systems that involve electron tunneling or other quantum mechanical systems.

Second, we have assumed that we have a function V() that represents the potential energy of all the atoms in a system under the influence of all the others. As we have already seen, the approximations made in developing an empiracle force field potential function are considerable, making this one of the most severe limitations of the method. This accounts for the many force fields/MD programs available (e.g., AMBER, GROMOS, CHARMM, CEDAR, DISCOVER).

B. The time scale of the molecular dynamics simulation or "experiment". This depends both on the process you wish to understand and on the time scale of the fastest motion in the system. The time step must be at least one-tenth the period of the fastest motion that you wish to include in your calculation (e.g., in the attached figure, this would be that of a bond vibration, ca., 10^-14 s, and the step would be 10^-15sec or a femptosecond). Note that it requires a thousand steps with such a step size to simulate a system for even 1psec!

C. The temperature of the molecular dynamics "experiment" is fixed by the initial velocities and positions given all the atoms, since these independent variables determine the kinetic and potential energies of each particle. Since the total energy of the system of particles must remain unchanged during the simulation, an "equilibration period" of several hundreds or thousands of time steps must be followed to allow the proper balance between kinetic and potential energy to be reached, and then the temperature can be obtained from the calculated average kinetic energy and the classical equipartition of energy principle that we encountered in our discussions of statistical mechanics:

E_kin = S_i ½m_iv_i² = 3/2Nk_BT

If the temperature of the "equilibrated" system is too high, the system is "cooled" by slowing all the particles, reequilibrating, then redetermining the temperature. This process can be continued until a close to desired temperature is found.

D. Uses of the molecular dynamics or molecular mechanics method. Several important applications of this method have been developed and are discussed in some detail in the review by Karplus & Petsko (1990). These include:

1. Estimating the structure and dynamics of proteins in solution. Of course, this is practical only if one already has an estimate of the structure (from X-ray and/or NMR data) and wants to know how the protein responds to a ligand or some other perturbation. Popular variations on this theme include focussing only on a portion of the protein (usually the binding pocket), calculating only V() as a ligand is manually positioned in the binding pocket then minimizing this function to estimate the energy of binding, etc. These and related methods have become very popular in the pharmaceutucal industry for the rational design of drugs.

2. Free energy changes associated with residue mutations or ligand changes. This is a useful tool in estimating the importance of a particular residue or of a ligand to the thermodynamic stability of the native state of a protein in solution. This is discussed in the excerpt from Karplus & Petsko on the next page.

3. Refinement of crystal structures and solution structures from X-ray and NMR data. This is accomplished by making the diffraction pattern R factor or the squared deviation from observed NOE distances as terms in the empiracle force field.

E. Evaluation and Prospects for Improvements. Shortcomings of the method derive largely from shortcomings of the force fields or from assumptions made in order to accomodate limited availabilty of computational time. Efforts to fix major shortcomings include the following:

1. the solvation problem: This is especially important for biomolecules. Karplus favors using statistical or averaging approximations such as solution of the Poisson-Boltzmann equation or use of

long-term electrostatic approximations (see below). Others argue that the highly polar and structured nature of water make such approaches unlikely to account for the structure of solvent water near a mac

romolecule. They advocate use of specific potentials designed to account for the properties of water and still be compatible with the force fields used for proteins or other biomolecules (see the SPC potential discussed after Stillenger's water).

2. environment sensitive molecular polarity: Real molecules (e.g., water) have different dipole moments in condensed phases than in the vapor phase (from where most experimental estimates come). Efforts to develop force fields containing molecular polarizabilities may help with this problem but are probably a long way from yielding useful results.

3. long-range forces: Electrostatic interactions are often quite long-range. Limiting these to a small cut-off sphere is expected to result in significant systematic errors. This is especially true in trying to predict hydration effects. Fast computational methods ("fast mutipole" methods, Ewald sums) are being developed to address this issue.

4. limitations on computer time: Since the time required to carry out MD simulations goes up as N² (N = # of atoms), these calculations quickly become impractical or impossible for biological macromolecules including solvent waters. Obviously, bigger and faster computers will eventually solve this problem, but people are working all the time to make the computations more efficient. Current efforts include i) mutiple step times (short for rapidly varying contributions, long for slowly varying ones) and ii) continuum models for water contributions.

VII. Stillinger's water

A. The Force Field. In the early 1970's, Ben Naim and Stillinger proposed a simple potential for the interaction between water molecules that was designed to take into account the tetrahedral nature of water, provide H-bond-like interaction energies and expected linear H-bond geometries, and produce a dipole moment in the range thought to occur in liquid water. This was accomplished by arranging point charges at the verticies of a tetrahedron, including a Lennard-Jones 6-12 potential, and including a switching function that turned off the electrostatic interaction when oxygens approached within predetermined separation limits (see diagram from Stillinger & Rahman, 1974). The minimum energy configuration for a pair of waters possessing this inter-molecular potential is consistent with reasonable expectations (see Figure from Stillinger & Rahman, 1974; note that the original potential used was that of Ben Naim & Stillinger = BNS, but that this was later replaced with a slightly less tetrahedral model = ST2).

B. The Simulation. Using this force field, Stillinger and his collaborators carried out a series of MD simulations of liquid water and aqueous solutions based on a system containing 216 water nolecules

(a 6x6x6 cube of water). Structural properties: The BNS force field was only moderately successful in reproducing the complex X-ray scattering pattern of water; it overestimated the intensity of the main peak and only very weakly produced the low angle shoulder on this peak, a feature unique to liquid water (Fig 7, Rahman & Stillinger, 1971). Since the main peak corresponds to tetrahedral O - O pair distributions, it was concluded that the potential favored too much tetrahedral arrangements of oxygens. This led to the ST2 potential, which has less of a tetrahedral tendency (Fig 3, Stillinger & Rahman, 1974). Simulation with this potential led to a much better fit of the X-ray scattering and even a reasonable description of neutron diffraction of D₂O (sensitive to D - D pair distributions) (see Figs 9 & 10 from Stillinger & Rahman, 1974; next page). The BNS potential also underestimated the dielectric constant, so that the point charges were moved farther from each other in the ST2 model, resulting in an excellent fit to the steady state dielectric constant.

Dynamic properties: Both the BNS and ST2 simulations led to reasonable estimates of the self-diffusion constant of water molecules (4·10^-5 cm²/sec; observed is 3·10^-5). D_self is likely most sensitive to the size of a molecule and to the interaction energy between molecules, both of which were similar for and probably reasonably estimated by the two potentials. The dipole-dipole relaxation time (t_d) was also reasonably estimated (8·10^-12 < t_d <11·10^-12; observed 7·10^-12). Thus, the shape of the orientational potential between water molecules is reasonably approximated by the tetrahedral point charge model.

Thermodynamic properties: U/N reflects the arbitrary

fixed total energy for each simulation run and can not be compared to experiment without a reference state. However, it's temperature derivative, C_V, is a critical thermodynamic peoperty of water. The simulation with ST2 overestimates C_V severely (35 vs 18 cal/mol·°K observed), although the simulation does predict the proper decrease with increasing temperature. Apparently, ST2

does not produce a realistic estimate of the energy difference between an H-bonded versus a non-H-bonded state, which is not unexpected given the clear quantum mechanical and non-additive nature of the water H-bond (Lentz & Scheraga, 1973).

This description illustrates the process by which comparison of simulation results with experiments can be used to obtain improved estimates of empiracle force fields, which then yield what are presumed to be better descriptions of the dynamic and structural properties of complex molecular systems.

C. Liquid Water. MD yields a view of water as consisting of a random, 3-D network of H-bonded molecules. The network has a local preference for tetrahedral geometry but deviates substantially from this over longer distances due to a large

number of strained or "broken" H-bonds. These strained or broken bonds produce a highly

dynamic or random network that accounts for the entropy of liquid water relative to ice. When ice melts, these broken or

strained bonds allow water

molecules to occupy less space by assuming more compact

(relative to ice) polyhedral network patterns (see Fig 3, p25). Conversion to more compact but more strained (higher E) patterns continues with increasing temperature, leading to the increase in molar

volume and large configurational C_V. Of course, the larger intermolcular distances between non-bonded molecules will lead to a competing volume expansion, accounting for the minimum in molar volume (experi

mentally at 4°C, at 27°C in the ST2 simulation).

D. The Hydrophobic Effect; Lennard-Jones Solutes: Geiger, Rahman & Stillinger used the ST2 potential in 1979 to simulate hydration of non-polar solutes (essentially, He or Ne atoms that replaced two water molecules in their water simulation). The Ne atoms had two major effects. First, water molecules in the first shell next to Ne reoriented so that their H-bond directions were all pointed toward other water molecules and never toward the Ne. Each Ne was surrounded by roughly 8 nearest neighbor and 6 next-nearest neighbor waters (see Fig 7, Geiger et al.). The reorientational correlation times of molecules in this 12 member first shell increased 1.2 to 1.7 fold relative to correlation times for molecules in the bulk. This corresponds roughly to a distorted clathrate-like structure similar to those seen in crytalline hydrates of organic molecules. Second, the Ne-Ne distance oscillates around about 3A for roughly 2.5psec, then switches to about 6A. Analysis of H-bonding patterns showed that this corresponded to a shift from two Ne encased in a single cage to each Ne having its own cage (i.e., a single layer of water being shared by the two Ne atoms - see Fig 13). In either case, the presence of the non-polar Ne solute induces a cage of more ordered and H-bonded water in its neighborhood. Although it was not possible to estimate the entropy change associated with the presence of Ne, it is clear it will be positive. The partial molar C_V was estimated from the simulation as 37 - 67 eu/mol solute (observed for Ne is 30-65). While the double cage structure is probably unique to the small Ne atom; larger hydrophobic solutes might be expected to induce similar ordered layers of water near their surfaces. Thus, this simulation supports the general view of the hydrophobic interaction first suggested by Frank & Wen and developed by Nemethy & Scheraga.

E. Simpler Water Model Needed for Protein Simulations. While the Stillinger ST2 potential does a

credible job reproducing the properties of water (a job for which it was designed), it is too complex

for use in protein simulations. Berendsen et al. (1981) has developed a Single Point Charge (SPC) potential that, while not reproducing the properties of water as completely as does ST2, does do a credible job and involves only 3 instead of 5 point charges, leading to a great reduction in calculation times.