7

Thermodynamics
Keith E. Gubbins
Department of Chemical Engineering
Cornell Universiiy
Ithaca, New York

I. Introduction
For the past 50 years, research in chemical engineering thermodynamics has
been largely concerned with problems related to the oil and petrochemical
industries. The primary methods of attack have been direct experiment and
macroscopic thermodynamic treatments of the data. Examples of the lat-
ter are empirical equations of state or activity coefficient expressions, group
contribution methods, and macroscopic corresponding states correlations.
These methods provide a convenient way to correlate large amounts of
experimental data, interpolate between different state conditions or between
different but structurally similar molecules, and make minor extrapolations
from experimentally measured conditions. They are most successful for fluids
of relatively simple molecules, such as the constituents of natural gas or
of low-molecular-weight hydrocarbon mixtures. However, they offer little
insight into the relation between the desired properties and the underlying
intermolecular forces or molecular structure and so are of little predictive
value. They will be much less useful for many of the new technologies, such
as the processing or design of electronic, photonic, or ceramic materials,
for biochemical processes, or for predicting the behavior of matter at and
near surfaces (e.g., in micelles, porous materials, or thin films).

125
Copyright 0 1991 by Audcmic Press. Inc.
ADVANCEBIN CHEMICAL ENGINERING. VOL. 16 All rights of reproduction in m y form rcmved.
126 Thermodynamics

A molecular approach based on statistical mechanics, combined where

possible with experiment and molecular simulation methods, will be a
necessary starting point for many of the most challenging problems of the
present and future. The rapidly increasing power of computer simulations
based on molecular modeling will play an important role in many appli-
cations, and an interdisciplinary approach will be essential for the new
technologies. In this paper, I shall focus on the present and future status
of the molecular theory and simulation approaches, rather than direct ex-
periment. This emphasis has been chosen because molecular theory and
simulation seem to me the areas that will move most rapidly in the next
few years; while experiment will continue to play a central role, the field
of thermodynamics has not been dominated by new experimental techniques
in the way that surface science, materials, and biotechnology have, for
example. I first provide a brief overview of the methods that are most likely
to prove useful for chemical engineering applications (Sect. 11), followed
by some examples of current problems (Sect. 111). Finally, some areas that
are likely to be important in the future are considered, together with im-
plications for the teaching of chemical engineering thermodynamics (Sect.
IV).

11. Modern Methods

Thermodynamic or transport properties can be studied by experiment, mac-
roscopic correlations, molecular theory, or computer simulation (Fig. 1).
In both the theory and simulation methods, one must specify an equation
for the intermolecular potential energy and then calculate the macroscopic
properties using the laws of statistical mechanics. The equations involve
N-body integrals (where N is the number of molecules) over functions of
the intermolecular potential energy which can be solved exactly in only a
few special cases (e.g., the ideal gas or gas adsorption at low pressure). The
theoretical methods therefore involve approximations to make the equations
tractable. The simulations, on the other hand, solve the full N-body equa-
tions by numerical methods; if sufficiently long computer runs are carried
out, the results should be exact for the model substance studied. A com-
parison of the simulation and theoretical results then provides a direct test
of the approximations used in the theory, without any uncertainty as to the
form of the potential. Thus, an important use of the simulations is to evaluate
new theories; those theories that meet this test can then be compared with
experimental data for real substances. If one neglects the test against simu-
lation and proceeds to compare theory directly with experiment, it is gen-
erally necessary to fit parameters in the potential (and often the theory) to
experimental data; this fitting often obscures defects in the theory, and much
Keith E. Gubbins 127

Theory t Potential Potential model

D MACROSCOPIC

CORRELATIONS
Correlation
I
EXPERIMENT

Figure 1. The four methods for studying physical properties and what may be learned
from a comparison of any two of them.

time can be wasted. The same difficulty arises with testing the macroscopic
correlations. An example is provided by the large amount of work done on
molecular theories of liquid mixtures (various I-, 2-, and 3-fluid confor-
mal solution theories, cell theories, random mixture theory, etc.) in the period
1936-1970, before the first simulations of mixtures were reported. Although
comparisons of theory and experiment had indicated good agreement, com-
parisons with computer simulations showed that most of the theories were
quite inaccurate, in some cases not falling on the same piece of graph paper
as the experimental results. A second important use of simulations is to evalu-
ate intermolecular potential models by making comparisons with experi-
mental data; in this case, the statistical mechanics is exact in both simulation
and experiment, so the only source of error is the potential model assumed
in the simulation. In some important applications, computer simulation may
provide the best way to obtain a detailed understanding of the molecular
behavior because suitable experiments cannot be devised at the present time;
examples are the breakdown of some classical thermodynamic equations
for small drops, the study of phase transitions in narrow pores, and the
dynamics and thermodynamics of protein folding.

A . Macroscopic Correlation Methods

These methods make use of available experimental knowledge combined
with correlation techniques that may or may not be partially based on
molecular concepts. They range from methods that are completely empirical
to ones in which there is significant input from statistical mechanics (see
Fig. 2). As an example of an approach that is almost fully empirical, one
can cite the many modifications of the Benedict-Webb-Rubin equation of
128 Thermodynamics

MACROSCOPIC APPROACH

< TEST VS. EXPERIMENT -P

I I I 1
UNBRIDLED BWR EOS LOCAL CORRESPONDING
EMPIRICISM COMPOSITIONS STATES
ACTIVlTY
COEFFICENT GROUP THEORY-BASED EOS
EQNS CONTRIBUTIONS

MOLECULAR APPROACH

-TEST VS. SIMULATION AND EXPERIMENT-

I
I I I I
CELL PERTURBATION DENSITY FULL
MODELS THEORY FUNCTIONAL AB
THEORY INlTIO
KINETIC
THEORY
Figure 2. Examples of macroscopic and statistical mechanical methods for studying
physical properties. The degree of molecular basis increases from left to right.

state, which involves expansions in density and temperature, with one or

more nonlinear terms often added for good measure. For a fluid of spheri-
cal Lennard-Jones or argon molecules, it is found that at least 34 adjust-
able constants are needed to obtain a fit to existing data over a reasonably
wide range of temperatures and pressures. To fit more complex fluids such
as water, many more terms are needed, and the equation still fails in the
critical region. Such methods are useful when one wants to fit a large amount
of data for a pure fluid; however, they cannot be extended to mixtures in
any systematic way and are unreliable even for phase regions slightly out-
side the range of fit.
Many of the macroscopic methods incorporate some features from mo-
lecular theory and are then usually more useful for prediction and extrapo-
Keith E. Gubbins 129

lation. Examples are equations of state of the modified van der Waals type
(Redlich-Kwong, etc.), local composition methods (which attempt to account
for the fact that the composition around a molecule of a particular species
differs from the mean composition of the mixture), group contributions
(which approximates the intermolecular forces as a sum of group-group
interactions), and corresponding states. A common feature of all of these
methods is that the desired macroscopic properties are not related explic-
itly to some expression for the intermolecular forces, in contrast to the mo-
lecular theories described in the following section. Because of this, it is not
possible to test the macroscopic methods against simulation results, but only
against actual experimental data. Since much fitting to this data has often
been involved, such comparisons are usually a weak test.
Before leaving these methods, we should note that it is often possible to
make a compromise between the molecular and macroscopic approaches by
starting from a sound molecular theory and adopting approximations to obtain
a semiempirical correlation that does not involve the intermolecular potential
or molecular correlation functions. Many examples occur in the chemical
engineering literature, and I mention only one, taken from recent work by
Bryan and Prausnitz [ l ] on developing an equation of state for polar flu-
ids. They write the equation for the compressibility factor Z = PVIRT in
the usual way as Z = Z,,f + Zpert, but in place of the hard-sphere reference
term they take Zref to be the equation of state for polar hard spheres, for
which an accurate statistical mechanical theory (i.e.. one that agrees closely
with the computer simulation results for polar hard spheres) exists. Such
a reference fluid is much closer to the real fluid of interest than a hard-sphere
one, so that the perturbation term Zpea is much smaller. Bryan and Prausnitz
used the mean field term of van der Waals, Zpert = -a/RTv, where a is the
van der Waals attraction constant and v is the molar volume. This approach
gives good results for polar fluids (see Fig. 3), in contrast to the many at-
tempts to doctor the perturbation term while retaining the traditional hard-
sphere reference.

B . Molecular Theory
The configurational part (the part involving the intermolecular forces) of
the Helmholtz free energy for a system of N molecules in volume V at
temperature T is given by

where drN = drldr2 * . drN, d d = dwldw2, U is the intermolecular po-

tential energy (which depends on all the rj and mi), and the integrations over
 1
130 Thermodynamics

7 I 1 1 8 1 1 I I I

6-
-Exp.
- .&C

5-
Acetonitrile 13
6
n -
z'a
-c
0-
-I - -
-2
I h i i'17 8 b i o i i i2
2 3
10001K-I
T
Figure 3. Calculated and observed vapor pressures for argon and several polar fluids
from an equation of state that uses a reference fluid of polar hard spheres, for
which an accurate statistical mechanical theory exists. Reprinted from Bryan
and Prausnitz [ 11 with the permission of Elsevier Science Publishers.

rj are over the volume V of the system. Here rj is the position of the cen-
ter of molecule i, and Oj (= 6 , p, for linear or 0, 6, for nonlinear mol- x
ecules) is its orientation relative to some space-fixed set of axes. Equation
(1) is valid for nonflexible molecules in which translational quantum ef-
fects are negligible. The integrations can be easily carried out for
noninteracting molecules (e.g., ideal gases), moderately dense gases, or
crystals, but for most other cases approximations are necessary. The prin-
cipal theories are [2] corresponding states theory, perturbation and cluster
expansions, density functional theory, lattice models, and integral equation
theory.
The first three of these theories are of particular interest in chemical en-
gineering. The molecular principle of corresponding states is based on the
idea that one can identify a group of substances, all of which obey a single
intermolecular potential law; they differ only in the values of the poten-
tial parameters. It provides the foundation for many existing correlations
of thermodynamic and transport properties, but these can be expected to
work only for groups of similar substances such as simple inorganics or low-
molecular-weight hydrocarbons.
Keith E. Gubbins 131

The perturbation theories [2, 31 go a step beyond corresponding states;

the properties (e.g., A,) of some substance with potential U are related to
those for a simpler reference substance with potential Uo by a perturbation
expansion (A, = A 0 + A1 + A2 + - The properties of the simple refer-
a).

ence fluid can be obtained from experimental data (or from simulation data
for model fluids such as hard spheres) or corresponding states correlations,
while the perturbation corrections are calculated from the statistical mechani-
cal expressions, which involve only reference fluid properties and the per-
turbing potential. Cluster expansions involve a series in molecular clusters
and are closely related to the perturbation theories; they have proved par-
ticularly useful for moderately dense gases, dilute solutions, hydrogen-bonded
liquids, and ionic solutions.
Density functional theories [2,4] are similar in spirit to the perturbation
theories, but are of particular value for problems involving nonuniform
systems in which the density (or molecular orientation) varies with posi-
tion (direction) in the system-surface phenomena, solidification and melting,
thin films, liquid crystals, polydisperse systems, and so on. In this approach,
one starts from the fact that the free energy density a ( r ) at some point r
in the system is a functional of the density profile p(r'). Using variational
methods, one finds the global minimum of the free energy density with
respect to the density profile and so determines the density profile itself.
The success of the method hinges on the accuracy of the expression used
for the free energy density functional a ( r ) .Usually this is written as a sum
of repulsive and attractive force contributions, the two terms being treated
by different approximations. Several versions of the theory exist, differing
in the approximations used for the repulsive term in a ( r ) , and these have
been successfully applied to problems in adsorption, micelles, fluids near
charged walls, and melting. This approach is likely to prove valuable for
the study of many of the interfacial problems to be met in the new tech-
nologies associated with electronic and microstructured materials, thin films,
etc.
Lattice models were used extensively to describe fluids from the 1930s
to the 1970s, but really describe solids rather than fluids, and they have been
superseded for most applications by the more sophisticated theories described
above. An exception is the study of the critical point, a singular point in
the phase diagram where conventional mean field and perturbation theo-
ries fail. The Ising model, a lattice theory that can be solved essentially
exactly, can be used successfully in this region and, together with some
related cell theories that can be mapped onto the king model (e.g., the deco-
rated lattice gas model), gives valuable information on the equation of state
in the critical region.
132 Thermodynamics

Integral equation methods provide another approach, but their use is limited
to potential models that are usually too simple for engineering use and are
moreover numerically difficult to solve. They are useful in providing equa-
tions of state for certain simple reference fluids (e.g., hard spheres, dipo-
lar hard spheres, charged hard spheres) that can then be used in the
perturbation theories or density functional theories.

C. Computer Simulation
In molecular simulation [ 5 , 61, one starts from a molecular model and an
equation for the intermolecular forces and calculates the macroscopic prop-
erties by a numerical solution of the equations. Two methods are in com-
mon use, the Monte Carlo and molecular dynamics techniques (Fig. 4).The
Monte Carlo (MC) method makes use of a random number generator to
“move” the molecules in a random fashion. Statistical mechanics tells us
that, for a fixed temperature and density, the probability of a particular
arrangement of the molecules is proportional to exp(-UlkT), where U is the
total energy of the collection of molecules, k is the Boltzmann constant, and
T is the temperature. In MC, the random moves are accepted or rejected
according to a recipe that ensures that the various molecular arrangements
that are generated appear with probabilities given by this law. After gen-
erating a long series of such arrangements, they can be averaged to obtain
the various equilibrium properties of the system of molecules.
These two techniques have several features in common. Accurate results
can be expected, provided that the simulation runs are carried on long enough
and that the number of molecules is large enough. In practice, the results
are limited by the speed and storage capacity of current supercomputers.
Typically, the number of molecules in the sample simulated can range up
to a few thousand or tens of thousands; for small molecules, the real time
simulated in MD is of the order of a nanosecond.
In order to minimize boundary effects in such small samples, it is cus-
tomary to use periodic boundary conditions; that is, the sample is surrounded
on all sides by replicas of itself, so that when a molecule moves through a
boundary and so out of the sample, it is automatically replaced by a mol-
ecule moving into the sample through the opposite face of the box. (Any-
one who has played Pacman, Asteroids, or similar video games is familiar
with this periodic boundaries trick.) Although molecular simulation can be
successfully applied to a wide range of problems, difficulties can arise with
some applications because of storage or speed limitations. Examples include
ionic fluids, such as plasmas and electrolytes, in which the range of the
intermolecular forces is very large, necessitating a large number of mol-
ecules. Difficulties also arise with substances in which long-range fluctuations
Keith E. Gubbins 133

N - -
100 10,000
Periodic Boundaries
Prescribed Intermolecular Potential

Monte Carlo Molecular Dynamics

Specify N,V, T Specify N, V, E

1
Generate random moves
1
Solve Newton's equations
5 = miq

Take averages Take averages

Obtain equilibrium properties Obtain equilibrium and

nonequilibrium properties
Figure 4. Two methods of molecular simulation. Typically, both treat a sample of N
molecules in a box of volume V (here shown as being of cubic shape). In the
Monte Carlo method it is common, but not necessary, to choose N, V, and
temperature T as the independent variables and to keep these fixed during the
simulation. Molecules are moved randomly, generating a new molecular
arrangement with the intermolecular potential energy U . These new arrange-
ments are accepted or rejected in such a way that those accepted occur with
the probability distribution that is required by the laws of statistical mech-
anics. In molecular dynamics the energy E, rather than the temperature T, is
fixed. The molecular positions r and velocities v (and also orientations and
angular velocities) for each molecule i are followed in time. In the molecular
dynamics (MD) method, the molecules are allowed to move naturally under
the influence of their own intermolecular forces. The positions and velocities
of each molecule are followed in time by solving Newton's equation of
motion (force equals mass times acceleration, a second-order differential
equation) using standard numerical methods. The macroscopic properties are
calculated by averaging the appropriate functions of molecular positions and
velocities over time.
134 Thermodynamics

occur; these are associated with substances very near critical points and those
exhibiting certain surface phenomena. Long-time phenomena are also apt
to make simulation difficult.
There are not only common features but also significant differences
between the MC and MD methods. MC is easy to program and can be easily
adapted to different conditions; it is adaptable, for example, to mixture studies
at constant pressure or to adsorption at constant chemical potential. MD is
more difficult to program and, in its conventional form, energy must be
conserved, providing a less convenient set of state variables to work with
for some applications. However, it is now possible to overcome this prob-
lem and to carry out MD calculations at constant temperature or pressure
[ 5 , 61. MD has two important advantages over MC; it can be used to study
time-dependent phenomena and transport processes, and the molecular
motions are natural and therefore can be observed and photographed eas-
ily using computer graphics.
Several specialized simulation techniques exist for particular applications
[ 5 ] : nonequilibrium MD for the study of transport properties and nonlinear
response, Brownian dynamics for the study of large molecules (e.g., pro-
teins) in solution, and quantum simulations for the study of non-classical
fluids and solids. Simulated annealing is a Monte Carlo technique for
optimization subject to a set of constraints and is finding widespread use
in design of chemical processes, circuit design (especially VLSI), image
processing, and protein engineering (see Sect. IV A).

D. Computer Graphics
A single molecular simulation provides a vast amount of detailed informa-
tion-typically lo8 coordinate positions and an equal number of orienta-
tions and linear and angular velocities. Even after averaging to obtain
macroscopic properties, the amount of molecular data (spatial, angular, and
time correlation functions, diffusion rates, and so on) is difficult to assimilate
from graphs or tables. Computer graphics provides a way to present such
large arrays of data and is particularly useful in visualizing such physical
processes as nucleation and phase separation, molecular motion near a sur-
face, adsorption and hysteresis in porous materials, and the folding of pro-
teins. In the simulation of very large molecules or materials, computer
graphics can help decide which structures or motions are important and which
functions would best characterize them. The ability to rotate the figure and
to zoom onto regions of particular interest is a great help in viewing such
processes. Although computer graphics is a rapidly advancing area, several
problems remain. The graphics workstations commonly available to uni-
versity research groups are still too slow to represent molecular motions
or rotations of the system in real time for most problems, and one must be
Keith E. Gubbins 135

satisfied to view a series of still shots, or make a movie of such a series to

see the motions. Moreover, the software needed for such displays must
usually be prepared locally and is not easily portable to different hardware.
E. Intermolecular Potentials
A key question about the use of any molecular theory or computer simu-
lation is whether the intermolecular potential model is sufficiently accu-
rate for the particular application of interest. For such simple fluids as argon
or methane, we have accurate pair potentials with which we can calculate
a wide variety of physical properties with good accuracy. For more com-
plex polyatomic molecules, two approaches exist. The first is a full ab initio
molecular orbital calculation based on a solution to the Schrodinger equa-
tion, and the second is the semiempirical method, in which a combination
of approximate quantum mechanical results and experimental data (second
virial coefficients, scattering, transport coefficients, solid properties, etc.)
is used to arrive at an approximate and simple expression.
The a b initio approach, while holding great promise for the future, suf-
fers from several difficulties at present. First, the calculations are very
lengthy, requiring long runs on very fast machines and massive data stor-
age. The computing power needed increases as N4 for Hartree-Fock cal-
culations (zeroth-order perturbation theory), where N is the number of
electrons per molecule; the rate of increase with N is even greater if elec-
tron correlations past the Hartree-Fock level (i.e., higher-order perturbation
terms) are included, up to N7 or N 8 , so that the calculations become rap-
idly more difficult as molecular size increases. Also, if the calculation is
done for just two molecules, the resulting pair potential takes no account
of the many-body forces that exist in dense gases or liquids; these are often
significant for dense fluids, particularly when induction forces (due to
molecular distortion in the presence of an electric field, e.g., due to a neigh-
boring polar molecule or ion) are present. These multibody forces can be
included in the ah initio calculations by increasing the number of molecules
present in the cluster, but at the cost of greatly increased computer time and
storage. Among the many pair potential models proposed for water are several
that come from ab inirio calculations. None of these are able to predict a
wide range of physical properties with the accuracy required for engineering
calculations at the present time. For example, the MCY (Matsuoka-Clementi-
Yoshimine) and CC (Caravetta-Clementi) potentials both predict pressures
at specified temperatures (or energy) and density (or volume), that are
considerably higher than the experimental values when used in computer
simulations, although they do predict correctly many of the anomalous prop-
erties of water [7]. For molecules with more electrons than water, the situation
is generally worse.
136 Thermodynamics

Because of the present difficulties facing the ah initio calculations, most

workers in this field use semiempirical expressions for the intermolecular
potentials. Among the older and more familiar ones are the Lennard-Jones
and Stockmayer models for spherical nonpolar and polar molecules, respec-
tively, and the Kihara model for molecules of nonspherical shape. More recent
developments have included group contribution models in which the pair
potential is made up of a sum of interactions between sites in the molecules;
the sites may be the nuclei of the atoms, groups of atoms (CH3, N 0 2 , etc.),
lone pair electrons, or bonds, and it is often assumed that such site inter-
actions are the same for different molecules having the same groups. Sev-
eral versions of this idea now exist and have been applied to large organic
molecules, particularly proteins and polypeptides. Examples include the TIPS
(transferable intermolecular potential functions) model, which works well
for water, alcohols, ethers, etc.; a modified form of TIPS is the OPLS
(optimized potentials for liquid simulations) model, in which the sites are
taken to be the nuclei or CH, groups [8]. Other examples of this approach,
developed mainly for conformational energy studies of proteins or polypep-
tides, but also applicable to a wide range of organics, include the ECEPP

Table 1. Liquids Simulated with OPLS Potential Functionsa

HCONH, 25 pyrrole 25
HCON(CH,), 25,100 pyridine 25
CH,CONHCH, 100 CH4 -161
CH,OH 25 C2H6 -89
C,H,OH 25 C3H8 -42, 25
,-C,H,OH 25 n-C4H 10 -0.5, 25
i-C,H,OH 25 I-C4H 25
t-C,H,OH 25 ,-'SH12 25
CH,SH 6 i-C5H12 25
C,H,SH 25 neo-CSHI , 25
(CH3)$ 25 C-CSHIO 25
C,H,SCH 25 n-C6H 14 25
(C,H,),S 25 CH,CH,CH=CH, 25
CH3SSCH3 25 I-CH,CH=CHCH, 25
(CH3)2O -25 c-CH,CH=CHCH, 25
C,H,OCH, 25 (CH,),C=CH, 25
(C2Hd20 25 benzene 25
THF 25 CH,CO,CH, 25
a From Jorgensen and Tirad+Rives [8].
Keith E. Gubbins 137

[9], CHARMM [lo], AMBER [ l l ] , and DISCOVER [I21 models. These

models differ in how they treat the intramolecular modes (some treat the
bond lengths as fixed, while others allow them to vary, for example) and
the kind of data (solid or liquid) used to fit parameters.
We briefly consider only the OPLS model, since it has been used to study
organic liquid properties. In OPLS, the sites are taken to be point charges
and/or Lennard-Jones centers placed on the nuclei or CH, groups, and the
total potential is a sum of Coulombic and Lennard-Jones terms. Such a
potential approximately reproduces the shape and the electrostatic and
dispersion forces. The necessary potential parameters for neutral molecules
are obtained by fitting directly to liquid property data, which requires
numerous simulations. For ions the parameters are derived by fitting to
structural and energetic results from ab initio calculations on ion-water
complexes. Some results of computer simulations using this potential model
for the 36 organic liquids listed in Table 1 are shown in Figs. 5 and 6 . The
model works quite well for a variety of aliphatic and aromatic hydrocar-

230

I95

160

u)
E! 125
0

90

55

20
20 55 90 125 160 195 230
EXPERIMENTAL
Figure 5. Comparisons of calculated and experimental volumes per molecule in A3 for
the liquids in Table 1 and TIP4P water. Calculated values are from computer
simulation using the OPLS potential method. Reprinted with permission from
W. L. Jorgensen and J. Tirado-Rives, J . Am. Chem. SOC.110, 1657 (1988) [8].
Copyright 1988 American Chemical Society.
138 Thermodynamics

15

12

v)
J 9
n
0

0
0 3 6 9 12 15 18
EXPER I M ENTAL
Figure 6. Comparison of calculated (computer simulation results for the OPLS model)
and experimental heats of vaporization in kcal mo1-l for the liquids in Table 1
and TIP4P water. Reprinted with permission from W. L. Jorgensen and J.
Tirado-Rives, J . Am. Chem. SOC.110, 1657 (1988) [8]. Copyright 1988
American Chemical Society.

bons, alcohols, amines, sulfur compounds, ethers, and so on, and has been
successfully applied to protein crystal structures and energy minimization
[81.
The models described so far use isotropic site-site interactions. They give
a good description of the molecular shape, but in most cases neglect the
rearrangement of the valence electrons that occurs on bonding. This shift
of the electrons into bonds, x orbitals, and lone pairs has a significant effect
on the intermolecular forces, and its description requires the use of aniso-
tropic site-site interactions [ 131. The electrostatic interactions between
molecules can be qu ite accurately described by using a series of sites within
the molecule (i.e., the nuclei of the atoms or CH, groups), each of which
interacts with sites on neighboring molecules with multipole forces (point
charge, dipole, quadrupole, and so on); an approach called distributed mul-
tipole analysis [13, 141 can be used to determine the best location of the
sites and the multipole moments needed from ab initio calculations. Such
Keith E. Gubbins 139

anisotropic site-site potentials have been successfully applied to a variety

of inorganic and aromatic molecules and give distinctly better results than
the isotropic site-site models. It is often possible to use fewer sites per
molecule, so that computational times may not be significantly greater than
with the simple isotropic site-site models. Also, the anisotropic atom-atom
potential can be more confidently transferred to other molecules where it
has a similar charge density.

111. Current Problems: Some Examples

Some current research areas using the methods described in the previous
section are shown in Table 2. Many of these involve molecular simulations
that exploit the limit of speed and storage of currently available
supercomputers [6]. In this section I shall consider two examples from among
these: (1) the determination of phase equilibria by computer simulation and
(2) the behavior of fluids in microporous materials.

A . Phase Equilibria of Complex Systems

An understanding of phase equilibria is of fundamental importance in a wide
variety of chemical processes, as well as in such diverse areas as biology,

Table 2. Some Current Research Areas

Continuous mixtures
Phase equilibria of mixtures
Polar and associating liquids
Electrolyte solutions
Solvation, folding of biological molecules
Gas hydrates (clathrates)
Polymers, advanced materials
Micelles, colloids, vesicles
Adsorption problems
Surfactants
Fluid behavior in porous materials
Nucleation
Thin films (Langmuir-Blodgett, etc.)
Chemical equilibria
Polydisperse fluids
Zeolites
140 Thermodynamics

geology, and the study of planetary atmospheres. Since the range of pos-
sible compositions, temperatures, and pressures that are met in practice is
enormous, it is not feasible to carry out experiments for more than a small
fraction of the systems of interest, and there is therefore much benefit to
be gained from developing computer simulation and theoretical prediction
methods. Some approximate estimates of the cost and time needed for such
calculations and experiments at the present time are shown in Table 3. Theo-
retical and empirical correlation methods are satisfactory only for rather
simple mixtures at the present time. For more complex fluids, simulation
or experiment is more reliable. Estimates of the cost and time for simula-
tions are strongly dependent on the complexity of the molecular model and
the accuracy desired and thus are difficult to make. At present, the simu-
lations are cheaper than experiments for simple mixtures but are still rela-
tively expensive for complex fluids, e.g., hydrogen-bonded ones; these costs
will decrease as faster machines become available. The simulations are in
general considerably faster than the experiments.
It is not straightforward to calculate phase equilibria or chemical poten-
tials in a simulation, and special techniques must be used. This has been
an active research area over the past few years [ 151. Several methods have
been proposed, and these can be divided into direct methods, in which the
coexistence properties of the phases are calculated directly, and indirect
methods, where the chemical potential is first calculated and then used to
determine the phase equilibrium conditions (Table 4).
The most straightforward direct method is to simulate a two-phase sys-
tem and allow it to equilibrate. This approach is valuable for studying the
properties of the interface but is less satisfactory for determining the prop-
erties of the bulk phases themselves because of slow diffusion across the
interface; the results are also sensitive to the interfacial area. An important
new development is the Gibbs ensemble Monte Carlo method [16], a di-
rect method that avoids the interfacial diffusion problem and is much faster
than simulating the two-phase system, especially for mixtures. The method
is illustrated in Fig. 7 and involves setting up two homogeneous phases
(I and 11) that are in thermodynamic equilibrium but not in physical con-
tact. Equilibration is achieved by allowing changes in the volumes and
number of molecules in each phase (keeping the total volume and number
of molecules for the two-phase system constant), together with the usual
Monte Carlo moves of molecules in each box to obtain equilibrium. Some
typical results for mixtures are shown in Fig. 8. Agreement with results from
other (indirect) methods is good, and the Gibbs method offers a great
improvement in speed because no interface is involved. Typically, for a binary
mixture of 600 molecules, the time required is five CPU minutes per mil-
Keith E. Gubbins 141

Table 3. Approximate Costs of Binary Vapor-Liquid Equilibria

Determinationa
Method cost Time
Redlich-Kwong eq. $10 0.1 hour
Perturbation theory $1~1000 1 hour

a The figures for the Redlich-Kwong equation and experiment are

taken from B. Moser and H. Kistenmacher, Fluid Phase Equilibria,
34, 195 (1987). The estimated costs and times for computations are
for a currently available supercomputer and are for a complete bi-
nary vapor- liquid phase diagram with about four or five isotherms.

lion molecular configurations on a Cray XMP supercomputer. The Gibbs

method can also be applied to equilibria in membranes or pores.
The major indirect methods [ 151 are ( 1 ) the test particle and grand ca-
nonical ensemble methods, which yield a chemical potential @) value in
a single simulation at the state point of interest, and (2) thermodynamic
integration, where one carries out a series of simulations at different state
points (e.g., different densities along an isotherm) and uses the standard
relations of classical thermodynamics to integrate along a path to get ,u at
the point of interest. The test particle method involves using an imaginary
test molecule (a ghost-like observer molecule) to measure the intermolecular
potential Ut due to the other (real) molecules; an exact equation relates p
to an integral over exp(-Ut/kT), so that attractive molecular configurations
(negative Ut) will be the dominant contributions to p. In the grand canonical
method, one uses as independent variables p, V , and T and calculates the
density in the simulation. This method is useful near critical points, e.g.,
in studies of supercritical extraction (an example is shown in Table 5 ) , since

Table 4. Methods for Determining Phase Equilibria by

Computer Simulationo

Direct Methods Simulate 2-phase system

Gibbs ensemble Monte Carlo
Indirect Methods Test particle method
Grand canonical ensemble method
Biased sampling methods
Thermodynamic integration
a From Abraham [6].
142 Thermodynamics

E' + AE', E' + AE",

N',V' N",V"
I
0
0 0 0 0
0
0 .
O b
0 0
0

\KIc
0 . 0 0

El, E' + AE', E" + AE:

N',V' N', V' + AV N", V"- AV

E'+AE', E'+AE",
N1+ I , V ' N"- I,V"
Figure 7. Possible steps in the Gibbs method for simulating the properties of fluids. The
schematic illustrates the initial system configuration and three steps: (a)
particle displacement, (b) volume change, and ( c )particle transfer. Variables
are defined as in Fig. 4.Reprinted with permission from W. L. Jorgensen and
J. Tirado-Rives, J . Am. Chern. SOC.110, 1657 (1988) [8]. Copyright 1988
American Chemical Society.

i t permits large density fluctuations in the fluid. The test particle, grand
canonical, and Gibbs ensemble methods work well at moderate densities
but become difficult to use at high densities, e.g., a dense liquid near its
triple point or a solid. The reason is most easily seen for the test particle
method, since the random insertion of a test molecule in a dense fluid will
almost certainly result in molecular overlap and consequently a large positive
value of Ut and a very small contribution to the integral that determines
p.The same problem occurs in the Gibbs and grand canonical methods with
the addition of new molecules to the system. What one needs to do is develop
some way of guiding the test particles or molecules toward any hole that
may be present in the fluid or solid. Several "biased sampling" methods have
Keith E. Gubbins 143

I I I I I
0 0.2 0.4 0.6 0.8 1.0
x2

(b)
1201 I I I I 1

EX PER1MENT
I00
P/bar
80

60

40

20

0
0 0.2 0.4 0.6 0.8 1.0
x2
Figure 8. (a) Vapor-liquid coexistence curves for a mixture of Lennard-Jones molecules
with parameters chosen to approximate acetone ( I ) - carbon dioxide (2):
0 = Gibbs method at constant pressure; A = test particle method at constant
volume. Horizontal and vertical lines are error bars. (b) Experimental (x) and
empirical equation of state (-) results for acetone-carbon dioxide mixtures.
Reprinted with the permission of Taylor & Francis Ltd. from Panagiotopoulos
et al. [ 161; and with permission from A. Z. Panagiotopoulos, U. W. Suter, and
R. C. Reid, Ind. Eng. Chem. Fundam. 25,525 (1986). Copyright 1986
American Chemical Society.
144 Thermodynamics

Table 5. Effect of Water Addition on the Solubility of Naphthalene (1) in Compressed

Carbon Dioxide (2) at 340 K from grand canonical MC simulations

Bl -5.0 -5.0 -5.0

B2 3.1 3.1 3.1
Nti 2 0 0 1 2
PO3 0.15fO.01 0.16 f 0.01 0.36 f 0.25
4 1 ’ 0.17f0.03 0.23 f 0.04 10.6 k 8.4
42’ 44.7 f 3.2 46.4 f 3.5 92.8 f 30.5
y,=N,lN 0.0038 0.0049 0.10
Here B, =p,JkT + I n N i , where pi,.is the residual (nonideal gas) part of the chemical
potential.
M. Nouacer and K. S. Shing, Molec. Simul. 2 , 5 5 (1989).

been developed for doing this; they extend the density range in which these
methods can be used. However, they still usually fail for solids or dense
fluids of highly nonspherical molecules (e.g., liquid crystals). For these more
difficult situations, one must resort to thermodynamic integration.

B. Fluid Behavior in Micropores

In microporous materials such as activated carbons, silica gel, porous glasses
and ceramics, clays, and zeolites, adsorbed fluids show many unusual and
interesting properties: ( 1 ) preferential adsorption of certain fluid compo-
nents, (2) hysteresis effects, and (3) a variety of new or unusual phase
transitions, including capillary condensation, wetting and prewetting phe-
nomena, layering transitions, and two-dimensional solid melting. For very
small pores, the behavior becomes characteristic of one-dimensional (cy-
lindrical geometry-no phase transitions) or two-dimensional (parallel plate
geometry) systems. These phenomena are poorly understood but play an
important role in separation and purification processes, catalysis, and many
geophysical and biological processes. They arise because of the strong in-
termolecular forces between the pore walls and the fluid molecules, which
dominate the fluid behavior for small pores; the fluid thus behaves quite
differently from the bulk. They are difficult to study experimentally because
the characteristics of most porous materials are not well defined; computer
simulation and molecular theory have therefore played an important role
in understanding them over the last two or three years. The long-term aim
of this work is to design adsorbents and membranes that will give improved
performance for specific applications in chemical processing.
Keith E. Gubbins 145

GCMC SIMULATION OF METHANE IN CARBON PORE

I I 1 I
0.8
kT/E = 0.5
R 7u Adsorption
N/V9
0 = rnol. dio. o Desorption
0.6

0.4

L
1st. LAYER
0.2

d
o k I I I I

Adsorption isotherm for a pure Lennard-Jones fluid in a cylindrical pore of

radius R = 70, from grand canonical Monte Carlo simulations. The system is a
crude model that approximates a methane-like fluid in a carbon pore at a
temperature in the liquid range; N is the number of fluid molecules, V is the
pore volume,pconis the configurational part of the chemical potential, and E
is the potential well depth for the fluid molecules. Increasing pconcorresponds
to increasing the pressure. On raising the pressure the system passes through
three layering transitions and finally condenses to liquid. The third layering
transition is believed to be a metastable state, and it as well as liquefaction is
accompanied by hysteresis. At low pressures the first layer is a two-dimen-
sional liquid, and this solidifies to a two-dimensional solid at somewhat
higher pressures (at pCon/&- -1 1). Reprinted with permission from B. K.
Peterson, G. S. Heffelfinger, F.van Swol, and K. E. Gubbins, J . Chem. Phys.
93,679 (1990).

The principal tools have been density functional theory and computer
simulation, especially grand canonical Monte Carlo and molecular dynamics
[17-191. Typical phase diagrams for a simple Lennard-Jones fluid and for
a binary mixture of Lennard-Jones fluids confined within cylindrical pores
of various diameters are shown in Figs. 9 and 10, respectively. Also shown
in Fig. 10 is the vapor-liquid phase diagram for the bulk fluid (i.e., a pore
of infinite radius). In these examples, the walls are inert and exert only weak
forces on the molecules, which themselves interact weakly. Nevertheless,
146 Thermodynamics

0 0.2 0.4 0.6 0.8 1.0

'A Or YA
Figure 10. Vapor-liquid equilibria for an argon-krypton mixture (modeled as a Lennard-
Jones mixture) for the bulk fluid (R* = m) and for a cylindrical pore of radius
R* = R/o,, = 2.5. The dotted and dashed lines are from a crude form of
density functional theory (the local density approximation, LDA). The points
and solid lines are molecular dynamics results for the pore. Reprinted with
permission from W. L. Jorgensen and J. Tirado-Rives, J . Am. Chem. SOC.
110, 1657 (1988) [8]. Copyright 1988 American Chemical Society.

condensation occurs at pressures far below the vapor pressure of the bulk
fluid (capillary condensation), and the critical temperature and pressure are
reduced substantially. For pores of intermediate size, e.g., a radius of 7.00
(0 is the molecular diameter), there are so-called layering transitions; these
are first-order phase transitions between adsorbed layers of one molecular
thickness and two, between two molecular layers and three, and so on. Such
transitions depend strongly on the strength of the wall forces, the temperature,
and the pore radius; for activated carbon pores, which interact more strongly
with the fluid molecules, the layering transitions would be more pronounced
than those shown in Fig. 9. For very small pores, they are inhibited by
molecular packing effects, while in very large pores they have less effect
on the overall phase diagram because of the large amount of bulk fluid
Keith E. Gubbins 147

present. When the pore radius becomes very small, a cylindrical pore ap-
proaches a one-dimensional limit in which the fluid molecules can only move
along a line. Exact statistical mechanics tells us that such a system cannot
show any phase transitions. This seems to occur at a radius between l a and
20. For mixtures, confinement within a pore will lead to large shifts in relative
adsorption and volatility, in addition to these effects (see Fig. 10). Much
remains to be done to understand these complex phase diagrams and ad-
sorption effects in terms of the underlying intermolecular forces, particu-
larly for the more complex fluids and porous materials of technological
interest.

IV. The Future

A. Research
Chemical engineering thermodynamics has been undergoing a transition
from an empirical, experimental field to one in which soundIy based mo-
lecular ideas play an increasingly prominent role. This trend has been made
possible by the rapid development of computers and of computer simula-
tion and molecular theory. Computers are expected to increase in power by
a factor of lo3 to lo6 over the next decade; such developments are bound
to have a major effect on the way we approach thermodynamics research.
The semiempirical method will continue to remain useful for many of the
traditional applications for thermodynamics research and will be most suc-
cessful when it has some molecular basis. However, the macroscopic ap-
proach will not be sufficient to meet the challenge of such new areas as
biotechnology, microstructured materials, and technologies that exploit the
unique properties of surfaces [20]. These areas will call for an interdisci-
plinary, molecular approach which combines the methods of simulation,
theory, and experiment. Among the areas that are likely to provide chal-
lenges and opportunities for thermodynamicists in the near future are the
following:

1. Computer-Aided Materials Design

We are moving toward an era in which simulation and molecular theory will
play an important role in the design of new polymers, composites, ceram-
ics, and electronic and photonic materials [21, 221. Much of the theoretical
and simulation work to date has asked questions of the form: what are the
properties of this particular model substance? We need to invert this ques-
tion and ask: here is a specific need, for example a difficult or expensive
separation-how can we use our theory and simulation techniques to design
a material or process to best meet this need [22]? One example would be
148 Thermodynamics

the design of improved porous materials as adsorbents, their composition

and surface characteristics being guided by a combination of simulation,
theory, and experiment. Molecular simulation techniques are already at a
stage where they can provide an understanding of the relation between struc-
ture and properties of these materials and can point the way to new designs.
The design of polymers having specific desired properties is also foresee-
able in the next decade [23]. The engineering of “designer” polymer mol-
ecules is in some ways more difficult than that of biological molecules (see
below), since in the latter case the structure of a single molecule or the in-
teraction of a pair of molecules (e.g., an enzyme and an inhibitor) is often
the main concern, whereas for polymers the bulk properties (mechanical,
thermal, etc.) are of major interest, so that it is necessary to simulate a sample
containing a reasonable number of molecules. In many of the materials ap-
plications, it will be necessary to combine molecular simulation methods with
quantum mechanical calculations of electronic structure [24].
I will give two examples drawn from recent work. The first is catalyst
design. It is now possible, using ab initio methods (generalized valence bond
with configuration interaction calculations), to estimate the energetics of
simple surface reactions to good accuracy [25]. Such calculations can help
to guide experiments on new catalysts and promoters, as well as in discrimi-
nation among various reaction mechanisms. These methods have recently
been used to study the epoxidation of olefins on silver catalysts [26]. The
next step in such work will be to use the energetics derived from the ab
initio calculations as input to molecular dynamics stimulations to study
reactions on catalytic surfaces [27]. With such methods it should eventu-
ally be possible to study the rates of elementary catalytic reactions. To do
this successfully, one must learn how to deal with the different time scales
involved in such reactions, a problem that will require some novel devel-
opments in the simulation techniques. The second example is the search for
high-temperature superconductors. This field has been almost entirely em-
pirical, but recently developed theories [27], though not yet fully tested,
may point the way to better materials and greatly speed up the search.
2. Surface Phenomena
Many new technologies rely on the unusual properties of interfaces-
Langmuir-Blodgett and other films, micelles, vesicles, small liquid drops,
and so on. Classical thermodynamics is often inadequate as a basis for treating
such systems because of their “smallness,” and experimental probes of the
interface are limited, especially for fluid systems. Computer simulation can
play an important role here, both in understanding the role of intermolecular
forces in obtaining desired properties and, in combination with experiment,
in designing better materials and processes [ 6 , 281.
Keith E. Gubbins 149

3. Biotechnology and Biomedicine

Thermodynamic stability plays a major role in the protein folding prob-
lem, in predicting the three-dimensional structure of enzymes from their
chemical composition, in predicting the limits of protein stability under a
variety of conditions, and in designing enzyme inhibitors [6, 21, 291. Quite
good potential energy functions now exist (OPLS, CHARMM, ECEPP, etc.-
see Sect. I1,E) for the configurational and intermolecular energies of many
biopolymers, including polypeptides, proteins, and nucleic acids, and these
can be used to predict the energies of the vast number of possible molecu-
lar conformations of these molecules [30]. A currently active research area
is the determination of the low-potential-energy conformations (i.e., local
and global minima) of such large flexible molecules, using one of the
potential energy functions mentioned above together with Monte Carlo or
some other systematic trial method; such calculations are also important for
biopolymers in solution, where it is thefree energy minima that are important.
These calculations will be particularly valuable when combined with high-
resolution proton NMR measurements on biopolymer molecules in solution.
Such measurements can accurately determine the distances between hydrogen
nuclei [31] and, when coupled with simulation methods, can yield the mo-
lecular structure in solution. It is now possible, using free energy pertur-
bation methods, to calculate changes in solvation energy and binding
constants for molecules of different structures [21, 32, 331. This technique
is beginning to play an important role in the design of new drugs which
are enzyme inhibitors. The free energy perturbation method involves a Monte
Carlo or molecular dynamics simulation with the potential drug molecule
bound to the active site on the enzyme (generally it is necessary only to
simulate a substructure of the enzyme); if the molecule is strongly bound,
it prevents the normal action of the enzyme [34]. Computer simulation of
such binding activity can be used to help design the drug molecule and to
guide experiments. A somewhat different method, simulated annealing [35],
can also prove useful in determining the most likely conformations of bio-
logical molecules [36]. The system is first simulated at a high temperature,
where thermodynamic equilibrium is easily attained, and then slowly cooled
(“annealed”); this helps to avoid becoming trapped in local minima.
Such calculations should eventually be able to provide a better understand-
ing of the thermodynamic fundamentals of biochemical separation meth-
ods, particularly for such processes as membrane and chromatographic
techniques; a clearer understanding of the causes of molecular clustering
will be important in this area. In the next five to ten years, these methods
should be much more sophisticated and will play a major role in the de-
sign of drugs, affinity agents, and proteins that fold into desired patterns.
150 Thermodynamics

Somewhat later it should be possible to design molecules with desired

biological properties.
4. Phase Equilibria and Fluid Properties
For the traditional oil and chemical industries, many of the most important
thermodynamic problems will continue to involve calculations of fluid phase
equilibria and physical properties of multicomponent fluid mixtures. Much
improved theories or theoretically based correlations are needed for mix-
tures containing molecules of disparate types-hydrocarbons with associ-
ating molecules and ions, for example. This is an area where we already
have a good understanding of the physical phenomena involved, but because
of the size of these industries small improvements in our predictive abil-
ity may lead to large economic gains. The computer simulation approach
is still some way from providing these small improvements because of the
slowness of current computers and the sensitivity of the results to small errors
in the intermolecular potentials. The development of the Gibbs method for
studying multicomponent phase equilibria [ 161 is a major step forward. For
simulation calculations to reach the stage where they can supplement or
replace experiments in a major way, we will need increases of computing
speed of lo2 or more coupled with better potentials, probably of the
aniostropic site-site form discussed in Sect. I1,E. On the experimental side,
there have been significant improvements in techniques [37], particularly
for measurements at high pressures and over wide ranges of temperature.
Phase equilibrium experiments are time-consuming, particularly because of
the time needed for equilibration; techniques for reducing the sample size,
while keeping surface effects minimal, are needed to improve this situa-
tion.
To summarize, molecular theory and computer simulation, in combina-
tion with experiment, will play a central role in many of the new areas,
particularly in helping to produce “designer microstructures” and “designer
molecules.” This will be particularly noticeable in many of the materials
areas in the next five years, for example, the design of porous materials,
polymers, and catalysts for simple reactions, the design of proteins and drugs
for particular purposes, and the molecular engineering of some thin films.
In these areas, classical or approximate quantum mechanical methods are
used, so that relatively modest improvements in computing speed and in
the accuracy of intermolecular potentials will have considerable impact.
However, the use of the full a b initio approach, involving a rigorous quanta1
treatment, seems to be much further away for problems of interest to chemical
engineers. Because the computing requirements rise so rapidly with the size
of the molecular species (typically, as the fifth or sixth power of the num-
ber of electrons), such calculations require much more computing power;
Keith E. Gubbins 151

in conventional simulations, where the molecules obey classical mechan-

ics, the computing power needed rises only as N 2 , where N is the number
of molecules. Nevertheless, quantum Monte Carlo methods based on
Feynman path integral techniques [38] are now being used to treat electron
transfer processes and in the long term should have an important impact
on electrode processes, semiconductor properties, electron transfer reactions,
and the use of lasers to transform matter. A spinoff of the molecular simu-
lation work has been the development of simulated annealing [35], widely
used in optimization of process and electrical circuit design but of future
interest in some of the new research areas; a key advantage of this tech-
nique is the avoidance of local minima in the function to be minimized (e.g.,
the free energy).
Chemical engineers working in these newer areas will need to become
versed in new fields and in many cases to work with scientists and engi-
neers from other disciplines. To deal successfully with biological systems
requires a group familiar with the intricacies of cell biology, biophysics,
and biochemistry, for example. Among the major international conferences
in applied thermodynamics are the International Conferences on Fluid Phase
Equilibria and Physical Properties for Chemical Process Design [39], held
every three years and alternating between Europe and North America, and
the IUPAC Conference on Chemical Thermodynamics [40].These confer-
ences will need to give more emphasis to these new areas and approaches
than they have in the past.

B . Teaching
These developments will call for a restructuring and rethinking of the teaching
of chemical engineering thermodynamics. Many widely used texts concentrate
exclusively on the classial approach, the word molecule never appearing
in the course. Statistical mechanics and quantum mechanics are usually taught
as part of the chemistry or physics sequence, but the student rarely if ever
sees any applications of this material in chemical engineering courses.
Implicitly, these subjects are treated like Greek mythology-part of the
student’s general education, but not of great importance in chemical engi-
neering. Such an approach may prevent our graduates from involvement in
important technologies of the future. We need to introduce examples into
our existing chemical engineering courses that apply the fundamental knowl-
edge of these topics that the student has gained in chemistry or physics. A
problem at present is the lack of suitable textbooks at the undergraduate
level. Books on statistical mechanics written by chemists usually do not
progress past applications to the ideal gas and other systems of noninteracting
particles, so that the student gets little feel for the relevance of the subject
152 Thermodynamics

to modern areas of technology. The challenge is to produce a book that is

usable at the undergraduate level but that has well thought out examples,
including application of the principles to nonconventional processes such
as thin films, materials, surface phenomena, and so on. It would be particu-
larly valuable to have a book that covers both classical and statistical ther-
modynamics at the junior or senior level. There is also a need for a graduate
level textbook on computer simulation methods. A possible approach to the
need for interdisciplinary education would be to develop more MS course
programs that cross traditional departmental lines. Obvious examples would
be degrees in chemical engineering and microelectronics and in chemical
engineering and biology.

Acknowledgments
I am grateful to many colleagues for helpful discussions and for putting up with many ques-
tions on areas outside my own direct experience. In particular, I thank F. H. Arnold, B. J.
Berne, J. C. G. Calado, E. A. Carter, P. T. Cummings, W. A. Goddard, W. L. Jorgensen, F.
Kohler, J. A. McCammon, A. Z. Panagiotopoulos, N. Quirke, H. A. Scheraga, W. C. Still,
and D. N. Theodorou. Part of this work was supported by grants from the National Science
Foundation and the Gas Research Institute.

References
1. Bryan, P. F., and Prausnitz, J. M., Fluid Phase Equilibria 38,201 (1987).
2. See, for example: Hansen, J. P., and McDonald, I. R., Theory of Simple Liquids, 2nd
Ed. Academic Press, London, 1986; Gray, C. G., and Gubbins, K. E., Theory of Mo-
lecular Fluids. Clarendon Press, Oxford, 1984; Lee, L. L., Molecular Thermodynam-
ics of Nonideal Fluids. Buttenvorths, Boston, 1987.
3. Chapman, W. G., Gubbins, K. E., J o s h , C. G., and Gray, C. G., Pure Appl. Chem. 59,
53 (1987).
4. Evans, R., Adv. Phys. 28,43 (1979).
5 . Allen, M. P., and Tildesley, D. J., Computer Simulation of Liquids. Clarendon Press,
Oxford, 1987.
6. Abraham, F. F., Adv. Phys. 35, 1 (1986).
7. Kataoka, Y., J. Chem. Phys. 87, 589 (1987).
8. Jorgensen, W. L., and Tirado-Rives, J., J. Am. Chem. SOC. 110, 1657 (1988), and ref-
erences therein.
9. Momany, F. A., McGuire, R. F., Burgess, A. W., and Scheraga, H. A., J. Phys. Chem.
79,2361 (1975); Nemethy, G., Pottle, M. S., and Scheraga, H. A., J. Phys. Chem. 87,
1883 (1983); Sippl, M. J., Nemethy, G., and Scheraga, H. A., J. Phys. Chem. 88,623 1
(1984).
10. Brooks, B. R., Bruccoleri, R.E., Olafson, B. D., States, D. J., Swaminathan, S., and
Karplus, M. J., J. Comp. Chem. 4, 187 (1983). See also Brooks et al. [30].
1 I. Weiner, S. J., Kollman, P. A., Nguyen, D. T., and Chase, D. A., J. Phys. Chem. 7,230
(1986).
Keith E. Gubbins 153

12. Dauber-Osguthorpe. P., Roberts, V. A., Osguthorpe, D. J., Wolff, J., Genest, M., and
Hagler, A. T., Proteins 4, 31 (1988).
13. Price, S. L., Mol. Simul. 1, 135 (1988), and references therein.
14. Stone, A. J., and Alderton, M., Mol. Phys. 56, 1047 (1988); Price, S. L., Stone, A. J.,
and Alderton, M., Mol. Phys. 52,987 (1984).
15. See Gubbins, K. E., Mol. Simul. 2,223 (1989), and references therein.
16. Panagiotopoulos, A. Z., Mol. Phys. 61.8 13 (1987); Panagiotopoulos, A. Z., Quirke, N.,
Stapleton, M., andTildesley, D. J., Mol. Phys. 63,527 (1988).
17. See, for example: Peterson, B. K., Gubbins, K. E., Heffelfinger, G. S., Marini Bettolo
Marconi, U., and van Swol, F., J. Chem. Phys. 88, 6487 (1988); Heffelfinger, G. S.,
Tan, Z., Gubbins, K. E., Marini Bettolo Marconi, U., and van Swol, F., Mol. Simul. 2,
393 (1989), and references therein.
18. Magda, J. J., Tirrell, M., and Davis, H. T., J. Chem. Phys. 83, 1888 (1985); Bitsanis,
I., Magda, J. J., Tirrell, M., and Davis, H. T., J. Chem. Phys. 87, 1733 (1987).
19. Walton, J. P. R. B., and Quirke, N. P., Mol. Simul. 2, 361 (1989).
20. National Research Council, Committee on Chemical Engineering Frontiers: Research
Needs and Opportunities. Frontiers in Chemical Engineering. Research Needs and
Opporrunities. National Academy Press, Washington, D. C., 1988; Krantz, W. B.,
Wasan, D. T., and Nerad, P. V., eds., Interfacial Phenomena in the New and Emerging
Technologies, Proceedings of workshop held at University of Colorado, May 29-3 1,
1986, National Science Foundation, Division of Engineering, Washington, D. C., 1987.
21. McCammon, J. A., Science 238,486 (1987).
22. See for example: Proceedings of the Conference on Industrial Applications of Molecu-
lar Simulation, Mol. Simul., Vol. 2 and 3 (1989).
23. For some recent work in this area see: Theodorou, D. N., and Suter, U. W., Macromol-
ecules 19, 139,379 (1986); Mansfield, K. F., and Theodorou, D. N., “Atomistic Simu-
lation of Glassy Polymer Surfaces and Glassy Polymer Solid Interfaces,” Annual
AIChE Meeting, Washington, D. C., Nov. 1988.
24. For an application to the properties of crystalline silicon see: Carr, R., and Parrinello,
M., Phys. Rev. Lett. 55,2471 (1985). For such treatments for a variety of materials, see:
Proceedings of the CCPS-CCP9 Conference on Computer Modeling of New Materi-
als, University of Bristol, UK, Jan. 4-6, 1989; Mol. Simul. 4, Nos. 1-3 (1989).
25. Carter, E. A., and Goddard, W. A., J. Chem. Phys. 88,3132 (1988); Carter, E. A., and
Goddard, W. A., J. Cutul. 112.80 (1988); Carter, E. A., and Goddard, W. A., Surface
Sci. 209, 243 (1988).
26. Carter, E. A., private communication (1988).
27. Guo, Y., Langlois, J.-M., and Goddard, W. A., Science 239,896 (1988); Chen, G., and
Goddard, W. A., Science 239,899 (1988). See also Science 242,31 (1988). This theory
explains high-T superconductivity in terms of magnetic interactions of electrons and
pairs of copper atoms, which leads to electron pairing; it predicts that the highest critical
temperature for the copper oxide superconductors under development is likely to be
about 225 K, about 100 K higher than it is now. The theory may also point the way to
improved superconductors based on materials other than copper oxide. For a discus-
sion of other theories, see: Emery, V. J., Physics Today Jan. 1989, pp. 5-26; Little, W.
A., Science 242, 1390 (1988).
154 Thermodynamics

28. For an example of the use of computer simulation to study micelles, see: Woods, M.
C., Haile, J. M., and O’Connell, J. P., J . Phys. Chem. 90, 1875 (1986). See also Davis,
H. T., in Perspectives in Chemical Engineering: Research and Education (C.K. Colton,
ed.), p. 169. Academic Press, San Diego, Calif., 1991 (Adv. Chem. Eng. 16).
29. Bash, P. A., Singh, U. C., Langridge, R., and Kollman, P. A., Science 236,564 (1987);
Kollman, P. A,, Ann. Rev. Phys. Chem. 38,303 (1987).
30. For a comparison of these methods, see: Hall, D., and Pavitt, N., J . Comp. Chem. 5,
441 (1984); also, Brooks, C. L., 111, Karplus, M., and Pettitt, B. M.,Adv. Chem. Phys.
71 (1988).
3 1. Wuthrich, K., NMR of Proteins and Nucleic Acids. Wiley, New York, 1986.
32. van Gunsteren, W. F., and Berendsen, H. J. C., J . Cornput.-Aided Molec. Des. 1, 171
(1987).
33. McCamrnon, J. A., and Harvey, S . C., Dynamics of Proteins and Nucleic Acids. Cam-
bridge University Press, New York, 1987.
34. Wong, C. F., and McCamrnon, J. A., J. Amer. Chem. SOC. 108,3830 (1986).
35. Kirkpatrick, S., Gelatt, C. D., and Vecchi, M. P., Science 220,67 1 (1983).
36. See, for example: Bringer, A. T., Kuriyan, J., and Karplus, M., Science 235,458 (1987).
37. See Section 7 of the Proceedings of the 9th IUPAC Conference on Chemical Thermo-
dynamics, Lisbon, July 1986. Published in: Pure Appl. Chem. 59 (1987).
38. Berne, B. J., and Thirumalai, D., Ann. Rev. Phys. Chem. 37,401 (1986).
39. The last of these conferences was held in Banff, Alberta, in May 1989.The conference
proceedings appeared in Fluid Phase Equilibria 52 (1 989).
40. The most recent of these conferences was held in Prague, Aug. 29-Sept. 2, 1988. The
proceedings appeared in Pure Appl. Chem. 61 (1989).