VAPOR-LIQUID EQUILIBRIUM CALCULATIONS WITH A NEW STEP POTENTIAL
MODEL AND EBULLIOMETRIC MEASUREMENTS
By
Abu Mokhtarul Hassan

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Chemical Engineering

2011

ABSTRACT
VAPOR-LIQUID EQUILIBRIUM CALCULATIONS WITH A NEW STEP POTENTIAL
MODEL AND EBULLIOMETRIC MEASUREMENTS
By
Abu Mokhtarul Hassan
Robust thermodynamic property prediction is vital to the design of chemical process
technologies. SPEAD (Step Potential Equilibria and Dynamics) is a rapid method of vaporliquid equilibrium property calculations using discontinuous molecular dynamics (DMD)
simulations in conjunction with thermodynamic perturbation theory (TPT). Though the potential
predicts vapor pressures of several compounds reasonably well, the liquid densities seem to be
systematically overpredicted. Equilibrium densities provide a vital test for the accuracy and
transferability of a potential, since thermodynamic properties are sensitive to small changes in
densities.
A new four step potential model called StePPE (Step Potentials for Phase Equilibria) is
developed, based on the theory of SPEAD, which reproduces the liquid densities significantly
better. The improvement in vapor pressure predictions for the compounds in this study is
moderate, but the accuracy is improved for a wider range of reduced temperatures, and is better
especially at lower reduced temperatures in comparison to SPEAD. Furthermore, the
transferability of parameters is observed to be better than SPEAD for a broader chain length
within a homologous series and from one homologous series to another, evidenced by a
reduction in the number of adjustable parameters for most of the compounds studied. A guideline
for obtaining the first guess of the diameters of the functional groups is developed based on the

ii

van der Waal’s volumes of the sites. The optimized site diameters are observed to be close to
97% of the first guess.
To provide vapor pressure data for compounds for which data is limited or unavailable in
literature, new ebulliometers are designed, which can accurately measure the boiling
temperatures of pure compounds using a small amount of sample. Activating the heated surface
with powdered glass and the use of Cottrell type tube were observed to be effective in removing
superheat. Condensate hold-up in these ebulliometers is estimated in order to use these
ebulliometers for the measurement of bubble points of binary mixtures, without measuring any
of the equilibrium compositions.

iii

Copyright by
ABU MOKHTARUL HASSAN
2011

iv

To my mother Haidery Begum, who fought against all odds to fulfill one of her very few dreams
of getting me an education

v

ACKNOWLEDGEMENTS

I thank Allah for any benefit this work can bring to anyone, and seek his forgiveness for
the shortcomings in me, from which I have tried to save this work.
Any accomplishment in this work would not have been possible without the support of
some talented individuals, and this short acknowledgement is not enough to appreciate their
efforts towards completion of this dissertation. I would first of all, like to thank my advisor Dr.
Carl Lira, for having the confidence in me to do this work. I still remember the first day when I
talked to him about possible projects in his group, when he showed me a beautiful ethyl lactate
molecule on his computer screen and asked me if I was interested in simulating properties of
molecules using molecular simulations. Though I was quite impressed with the visuals, it wasn’t
enough to dispel my fear of anything that had to do with programming triple integrals and
quadruple summations. I however started working on it somewhat, just to get a ‘flavor’, and
about fifty run-time errors later, I realized I had underestimated Dr. Lira’s patience more, than I
had underestimated my abilities. I remember the many discussions we had that lasted four hours
and several cappuccinos, but in the end they were all worth it. I am thankful to him for his help,
not only as a good research supervisor, but also for being a good mentor and friend outside
research.
I am thankful to my committee members Dr. Nikolai Priezjev, Dr. Micahel Feig and Dr.
Dennis Miller (not in the order of people who asked me the most difficult questions). While Dr.
Feig’s quantum level questions helped me grasp some of the minutest details in our model, Dr.
Priezjev’s molecular level questions helped me understand some of the shortcomings of our
model and think about improvements. Dr. Miller’s uncanny knack of extracting out the smallest

vi

details from a pile of data is almost superhuman. Without the thought provoking comments and
questions from these individuals, I wouldn’t have had a complete understanding of this work.
My lab mates in Lira/Miller group deserve special thanks for developing a workplace
with learning and nurturing environment. Our post docs Dr. Lars Peerebom, Dr. Aspi Kolah and
Dr. Ambareesh Murkute have been very helpful in making me learn to be competitive yet
collaborative. Venkata Sai has been a good friend with whom I have often shared my ups and
downs during my studies. Alvaro Orjuela, Anne Lown, Xi Hong, Arathi Santhankrishnam and
Aron Oberg have been very good lab mates. Undergrads Joe Zalokar, Pete Rossman, Dushyant
Barpaga, Jonathan Evans and Alex Smith have been a crucial help in both, my experimental and
computational endeavors.
I am also thankful to my friends and teachers from undergrad. Tabish Maqbool deserves
special words of appreciation for being a constant support throughout my career and a true friend
who I could always rely on. I will also always remember my days with Yusuf Ansari and Zaid
Bin Khalid who have been best of friends. My undergrad teacher Miss Seemi Rafique is one of
the important individuals who helped in my early days when I had been struggling, not really
knowing whether I wanted to go for higher studies. Dr. Hamid Ali’s constructive criticism is
something I will always appreciate.
My days in East Lansing would not have been as enjoyable without some important
friends. I have enjoyed and company and learned something from all of them. My earlier
roommate Farhan Ahmad’s passion for research is unmatched. It’s been hard for me to persuade
him to talk about anything but research, though I’m somewhat at a loss trying to figure out his
penchant, bordering on monomania, for a certain chicken curry, cooked exactly in a certain way
with precise amounts of certain ingredients! Rehan Baqri’s confidence in talking about subjects

vii

like chemical engineering and world politics, and anything that lies in between amazes me. He’s
a neuroscientist by the way. Basir Ahmad’s simplicity and willingness to help, sometimes by
putting himself in trouble, is something really scarce these days. I thank all my other friends in
the Islamic Center of East Lansing and our little basketball group in Spartan Village for making
my days in East Lansing memorable.
I would also like to acknowledge the help of the staff in the CHEMS office. Joe Ann,
Lauren Brown, Jennifer Peterman, Nikole Shook, Donna Fernandez and Kim McClung, have
been our hands and feet, and I don’t even want to imagine the mayhem the department would be
in without them.
Finally, and importantly, I thank my family for their selfless support. Ma and baba, may
Allah always take care of you as you took care of me when I was little. Your unflinching faith in
me and unconditional support is something I cannot compensate for, no matter what. The
sacrifices of my brothers Zafar and Hassan to make me reach where I am today is something I
will be indebted to all my life. I will always cherish the boundless love and affection of my
sisters Ishrat and Zakia. I cannot imagine going through this journey without the love and
support of my dear wife Zahida, who has always stood beside me in sharing my difficulties and
celebrating my achievements. I thank Allah for blessing me with this wonderful family.
I cannot count the bounties Allah has bestowed upon me, and the people mentioned here
are my greatest blessings. As I embark upon a new journey, I ask Allah for his guidance.

viii

TABLE OF CONTENTS

LIST OF TABLES......................................................................................................................... xi
LIST OF FIGURES ..................................................................................................................... xiii
1. INTRODUCTION ...................................................................................................................... 1
1.1 Research Overview ............................................................................................................... 1
1.2 Objectives and overview of thesis ........................................................................................ 2
1.3 References............................................................................................................................. 7
2. OVERVIEW OF VAPOR-LIQUID EQUILIBRIUM CALCULATION WITH SPEAD
METHOD ....................................................................................................................................... 8
2.1 Introduction........................................................................................................................... 8
2.2 DMD simulations for thermodynamic properties ............................................................... 10
2.3 SPEAD background ............................................................................................................ 11
2.4 References........................................................................................................................... 25
3. PREDICTION OF LIQUID DENSITIES AND VAPOR PRESSURES WITH SPEAD:
MOTIVATION FOR A NEW STEP POTENTIAL MODEL...................................................... 28
3.1 Summary of results with SPEAD: model parameters and the liquid density errors ........... 28
3.2 Improvements: the StePPE model ...................................................................................... 30
3.3 Summary of comparison of errors with SPEAD and StePPE............................................. 37
3.4 References........................................................................................................................... 40
4. LIQUID DENSITIES AND VAPOR PRESSURES OF PURE COMPOUNDS WITH
StePPE MODEL ........................................................................................................................... 41
4.1 Site diameter and step potential optimization: summary.................................................... 41
4.2 Discussion of results with StePPE parameters and comparison with other UA models .... 43
4.3 References.......................................................................................................................... 70
5. MEASUREMENTS OF VAPOR PRESSURES AND BUBBLE POINTS WITH
EBULLIOMETRY ....................................................................................................................... 74
5.1 Introduction......................................................................................................................... 74
5.2 Ebulliometer designs for pure compounds ......................................................................... 75
5.3 Ebulliometer for measuring T-x of a binary liquid............................................................. 79
5.4 References........................................................................................................................... 84

ix

6. CONCLUSIONS AND RECOMMENDATIONS FOR FUTURE WORK............................. 86
6.1. VLE calculations with StePPE........................................................................................... 86
6.2 Ebulliometric measurements............................................................................................... 96
6.3 References........................................................................................................................... 99
APPENDICES ............................................................................................................................ 102
A.1 Angle distribution in butyl propionate ............................................................................. 102
A2. Bond lengths .................................................................................................................... 103
A3. Finding the first guess of site sizes from excluded volumes............................................ 104
A.4. Liquid density and vapor pressure plots ......................................................................... 106
A.5. Tail corrections ............................................................................................................... 117
A.6. Ratio of the number of intramolecular attractive pairs to total attractive pairs .............. 119

x

LIST OF TABLES

Table 3.1: Summary of errors in vapor pressure and liquid density predictions with
SPEAD and StePPE ...................................................................................................................... 37
Table 4.1: Summary of optimized parameters for sites used in this study. The sites within
parentheses in the first columns are the sites to which the group in that column is bonded,
or a special case. For example ‘CH3 (-CH2-/-CH3)’ represents a methyl group, which is
bonded to another methyl or methylene group. The ‘>’ or ‘<’ signs mean two sigma
bonds and the ‘=’ sign denotes a double bond. For example, the carbonyl carbon, which
has two sigma bonds and a double bond is represented by ‘>C=’. θ is the bond angle. The
term ‘σ’ denotes the site diameter and the term ‘k’ in well energies is the Boltzmann
constant. ........................................................................................................................................ 42
sat

min

Table 4.2: Vapor pressures and liquid densities of n-alkanes with StePPE. P , Tr ,
max
and Tr
in this and the subsequent tables mean vapor pressure, minimum reduced
temperature and maximum reduced temperature.......................................................................... 44
Table 4.3: Vapor pressures and liquid densities of branched alkanes with StePPE ..................... 48
Table 4.4: Vapor pressures and liquid densities of alkenes with StePPE..................................... 55
Table 4.5: Vapor pressures and liquid densities of primary alcohols with StePPE...................... 57
Table 4.6: Vapor pressure and liquid density deviations of secondary alcohols with
StePPE........................................................................................................................................... 61
Table 4.7: Vapor pressure and liquid density deviations of ethers with StePPE.......................... 62
Table 4.8: Vapor pressure and liquid density deviations of ketones with StePPE ....................... 65
Table 4.9: Deviations in liquid densities and vapor pressures of esters with StePPE .................. 68
Table 6.1: Comparison of errors of secondary alcohols with two different sizes of the
methine group ............................................................................................................................... 89

xi

Table 6.2: Errors in liquid densities and vapor pressures of n-alkanes with StePPE with
tail corrections............................................................................................................................... 91
Table 6.3: Comparison of the radius of gyration of n-alkanes with the most likely end-toend distance in StePPE simulations .............................................................................................. 96
Table A2.1 Bond lengths used in StePPE model........................................................................ 103

xii

LIST OF FIGURES

Figure 1.1: Comparison of methods for performing molecular simulations .................................. 3
Figure 2.1: Illustration of bonds and pseudobonds for three adjacent sites in an alkane
represented as overlapping spheres. The solid lines a and b represent the covalent bonds
between the adjacent sites (1-2, 2-3) and the dashed line c represent the pseudobond
between sites 1 and 3 to constrain the angle θ .............................................................................. 11
Figure 2.2: Illustration of the multistep potential used in SPEAD. The innermost well
represents a covalent bond and the second inner well represents pseudobond to constrain
bond angles ................................................................................................................................... 12
Figure 2.3: Finding the pseudobond length (given by the dashed lines) for cis/trans
interactions in cis/trans butene...................................................................................................... 13
Fig 2.4: SPEAD simulation steps to get the EoS .......................................................................... 21
Fig 2.5: Flowsheet of equilibrium pressure and density calculations. Symbol meaningsPr : Reduced pressure, Tr : Reduced temperature, ω: Compressibility factor, b : Core
volume (volume/molecule), η : Packing fraction (V: vapor, L: liquid), f : Fugacity, A :
Helmholtz energy departure.......................................................................................................... 22
Figure 3.1: Published results of vapor pressure and liquid densities using SPEAD..................... 28
Figure 3.2: Density predictions of normal alkanes with Akron SPEAD. The dots are
experimental points and the lines are model predictions .............................................................. 29
Figure 3.3: Optimized site diameters (ordinate) compared to diameters obtained from van
der Waal's volumes (abscissa) in Å; (a) StePPE, based on 20% volume correction per
bond, (b) TraPPE, based on 24% volume correction per bond..................................................... 33
Figure 3.4: Core volumes predicted with StePPE plotted against van der Waal’s volumes
of alkanes. The x-axis denotes the number of C atoms in the alkane chain ................................. 34

xiii

Figure 3.5: Comparison of A2 using polynomial function as in Akron- SPEAD and
exponential function as in MSU-SPEAD for trioxane.................................................................. 35
Figure 3.6: Liquid density deviations of selected compounds with SPEAD (filled
symbols). Symbols are: diamonds- alkanes, circles- 1-alkenes, triangles- ketones,
squares- acetates. The open symbols are with StePPE ................................................................. 38
Figure 4.1: Liquid densities from left to right of (a) n-propane, n-pentane and n-octane,
and (b) n-dodecane and n-hexadecane. The symbol meanings are same in (a) and (b) ............... 45
Figure 4.2: Vapor pressures, from right to left, of, n-propane, n-pentane, n-octane, ndodecane and n-hexadecane. For interpretation of the references to color in this and all
other figures, the reader is referred to the electronic version of this dissertation......................... 46
Figure 4.3: (a) Liquid densities of (left to right) isopentane, neopentane and 2,2dimethylhexane; (b) saturation pressures of (right to left) isopentane, neopentane and 2,2dimethylhexane. The legend meanings are same in (a) and (b).................................................... 49
Figure 4.4: Liquid densities (a) and saturation pressures (b) of 3,4-dimethylhexane and
2,5-dimethylhexane. The lines are from correlations in DIPPR, filled symbols are from
TraPPE and the open symbols are from StePPE........................................................................... 50
Figure 4.5: Experimental liquid densities of C9 isomers and 2,2,3,3- tetramethylbutane.
The curves are from equations in DIPPR 801 database................................................................ 51
Figure 4.6: (a) liquid densities (left to right) of propene, 1-butene, 1-pentene and 1octene; (b) vapor pressures of (right to left), propene, 1-butene, 1-pentene and 1-octene.
The symbol meanings are same in (a) and (b) .............................................................................. 53
Figure 4.7: liquid densities (a) and vapor pressures (b) of butene isomers. The meanings
of the legends are same in (a) and (b) ........................................................................................... 54
Figure 4.8: Liquid densities (a) and vapor pressures (b) of (right to left) propanol, 1pentanol and 1-octanol. In (a) the dotted, solid and dashed lines are for propanol, 1pentanol and 1-octanol respectively.............................................................................................. 58
Figure 4.9: Liquid densities of selected 1-alcohols with SPEAD and StePPE. Dashed
lines are with SPEAD and the solid lines are with StePPE .......................................................... 59

xiv

Figure 4.10: Liquid densities of methanol (a) and ethanol (b) with SPEAD and StePPE.
The symbol meanings are same in (a) and (b) .............................................................................. 59
Figure 4.11: Liquid densities and vapor pressures of 2-butanol (circles) and 2-propanol
(squares) with TraPPE and SPEAD. The closed symbols are with TraPPE and the open
symbols are from StePPE.............................................................................................................. 60
Figure 4.12: (a) liquid densities of (left to right) dimethyl ether, dipropyl ether and dipentyl ether and di-nonyl ether; vapor pressures of (right to left) dimethyl ether, dipropyl
ether and di-pentyl ether and di-nonyl ether ................................................................................. 63
Figure 4.13: Liquid densities (a) and vapor pressures (b) of acetone, 2-pentanone and 2octanone with SPEAD, TraPPE and StePPE. In (a), the filled symbols are from TraPPE
and open symbols are from StePPE. In (b), from top to bottom: acetone, 2-pentanone, 2octanone ........................................................................................................................................ 66
Figure 5.1: Ebulliometer based on Hoover’s design for vapor pressure measurement ................ 76
Fig 5.2: Cottrell tube ebulliometer. In (a), the entire assembly is shown, (b) shows the
Cottrell tube insert with a funnel shaped bottom to catch the vapors, which get carried to
the equilibrium chamber with the boiling liquid. The contents get squirted on the
temperature probe (not shown) in this chamber. The height of the outside tube is about 6”
and the diameter is about 1”.......................................................................................................... 77
Figure 5.3: Ebulliometer with Washburn pump ........................................................................... 77
Figure 5.4 Stirred round bottom flask ebulliometer for vapor pressure measurements................ 78
sat

Figure 5.5: Measurement of P of monoethyl succinate with Cottrell, Hoover and stirred
flask ebulliometers ........................................................................................................................ 78

Figure 5.6: Measurements of water P

sat

with Cottrell and Washburn ebulliometers ................... 79

Fig 6.1: Liquid densities of neopentane and 2,2-dimethylhexane with StePPE. The
symbols are from DIPPR801 correlations .................................................................................... 88

xv

Fig 6.2: Vapor pressure and liquid density deviations with StePPE, using ‘34’ and ‘77’
intramolecular interactions; x-axis denotes the total number of carbon atoms in the ester
chain.............................................................................................................................................. 92

Fig 6.3: Ratio of intramolecular attractive pairs (NI) to total number of attractive pairs
(NT) with x. The legends denote the total number sites in the molecule ...................................... 93
Figure 5.4: Vapor pressures of ethers with StePPE with and without correction for partial
charges .......................................................................................................................................... 95
Figure A.1: Distribution of some angles in butyl propionate. The legends denote the three
sites that form the angle, and the site within the parentheses is the central angle. The
vertical lines denote the equilibrium values (109.5° and 120°) ................................................ 102
Fig A3.1: Schematic for finding site diameters from given excluded volumes.......................... 104
Figure A.4. 1: Liquid densities (a) and vapor pressures (b) of n-alkanes ................................... 106
Fig A.4.2: Liquid densities of branched alkanes (a); isobutane (b) ............................................ 107
Fig A.4.3: Vapor pressures of branched alkanes (a); isobutane (b)............................................ 108
Figure A.4.4: Liquid densities (a) and vapor pressures (b) of alkenes ....................................... 109
Figure A.4.5: Liquid densities (a) and vapor pressures (b) of primary alcohols. The
symbol ‘1olCx’ denotes a primary alcohol with ‘x’ carbons in the backbone. For
example, 1olC8 in 1-octanol ....................................................................................................... 110
Figure A.4.6: Liquid densities (a) and vapor pressures (b) of secondary alcohols..................... 111
Figure A.4.7: Liquid densities (a) and vapor pressures (b) of ethers.......................................... 112
Figure A.4.8: Liquid densities (a) and vapor pressures (b) of ketones ....................................... 113
Figure A.4.9: Liquid densities (a) and vapor pressures (b) of acetates....................................... 114

xvi

Figure A.4.10: Liquid densities (a) and vapor pressures (b) of methyl esters ............................ 115
Figure A.4.11: Liquid densities (a) and vapor pressures (b) of esters higher than acetates ....... 116

xvii

CHAPTER I
1. INTRODUCTION

1.1 Research Overview
Phase equilibria plays a critical role in a variety of technical and industrial applications
ranging from combustion of fuels inside an engine to process technology development in the
chemical industry. From applying adhesives and coatings, to simply melting ice, phase behavior
is also at the core of many day-to-day applications. Physical property calculations are at the core
of process engineering design because the properties determine the sizes for equipment and
utility demands, which in turn govern economic profitability of processes.
This thesis encompasses both, experimental and modeling aspects of fluid properties,
primarily focusing on predictions and measurement of vapor pressure. Recent advancements in
computational efficiency have enabled the application of advanced molecular simulation
techniques to generate phase equilibria of relatively complex systems, using atomistic or united
atom (UA) models. Vapor-liquid coexistence properties, which are available from experiments
for a variety of systems, can be used to parameterize, or test the accuracy and transferability of
the force-fields developed for the different simulation schemes. Liquid density has been the
preferred property for force-field optimization and testing, since system thermodynamic
properties are highly sensitive to density. Vapor pressure is another key equilibrium property of
interest for the study of potentials, which is also critical in industrial applications such as
chemical separation[1] and liquid fuel combustion[2]. The prediction of accurate vapor pressures
can be leveraged to secondary properties by integration, differentiation, or other correlation. For
example, the vapor pressure governs the flash point of a fuel[3] and is related to its surface
tension and heat of vaporization[4].
1

This thesis also presents the development and testing of various ebulliometers. To
measure vapor-liquid equilibrium (VLE) experimentally, either static equilibrium cells, or
dynamic equilibrium circulation stills are used. The static method is more accurate, but it is slow,
and the auxiliary equipment needed makes it expensive[5]. On the other hand, ebulliometry,
which is the most widely used dynamic method, is fast, inexpensive, and has been improved over
time to give extremely accurate results[6-8]. Ebulliometry works on the principle that at a given
temperature, a liquid will boil at its vapor pressure. In this method, rather than measuring the
bubble pressure, the bubble temperature of the fluid is measured at a fixed pressure.
In this work, we have used ebulliometry as our method of measuring VLE and have
developed an ebulliometer to be used for vapor pressures of pure compounds in a broad pressure
range of 30 mbar to atmospheric. The focus is also to develop ebulliometers versatile enough for
both, pure compounds and binary mixtures.

1.2 Objectives and overview of thesis
The objective of this work is to develop a tool for engineering predictions of equilibrium
properties. A spectrum of simulation methods are available spanning from quantum mechanics to
empirical group contribution methods. As shown in Figure 1.1, the level of confidence increases
with the level of detail, but the computational costs for quantum mechanical calculations of
systems of molecules are rarely justified. On the other hand, group contribution methods for
predictions of properties such as vapor pressure typically require estimation of critical properties
and then apply a corresponding states model to estimate the properties. Such methods have been
shown to lead to large and unpredictable errors. On the other hand, group contribution methods
such as UNIFAC use group-based energies and are fairly reliable for predicting activity
coefficients in mixtures. For the purposes here, it is desirable to include more detailed molecular
2

information, but also to compromise on the level of detail to provide a robust predictive method.
The existing model called SPEAD (Step Potentials for Equilibria And Dynamics) has been
shown to provide rapid accurate estimates for vapor pressure. The method uses realistic
simulation of the molecular geometry, but simplifies the pair potentials and the intramolecular
vibrations and torsions. The method is also limited because it does not provide an direct accurate
prediction of the critical point due to the limitations of the truncated series expansion used to
represent the Helmholtz energy. However, the method is accurate enough in the engineering
ranges needed for fuel properties that reexamination of the capabilities is justified. The method
appears to provide a good balance of realistic simplifications and molecular detail to provide
good engineering estimates of vapor pressure and density.

Quantum
Mechanics

All atom
molecular
simulation

United atom
molecular
simulation

Empirical
Corresponding
States
Prediction

Level of detail and computational cost
Speed of prediction
Figure 1.1: Comparison of methods for performing molecular simulations

The purpose of this study is to develop a transferable four step potential model similar
SPEAD for accurately predicting the equilibrium densities, which gives the vapor- liquid
coexistence curves (VLCC), vapor pressures and heats of vaporization of fluids. This study will
address the inadequacies in the SPEAD model and provide a fundamental guideline for site
diameter selection based on van der Waal’s volumes. All the improvements are aimed at

3

improving transferability, and wherever possible, reducing the number of parameters. Though
the other UA models reproduce the liquid densities well, vapor pressure, which is a relatively
important fluid property in the chemical/petrochemical industry, is not reproduced accurately
with these models, especially for longer molecules. This model would thus be a much faster
alternative compared to traditional UA models for accurate property predictions. The goals of the
work are to provide predictions for ketones, esters, ethers, and alcohols which are all potential
biofuel molecules. Further, we would like to extend the work to higher molecular weights typical
of diesel-type fuels.
In order to complement the model, ebulliometers for accurate vapor pressure
measurement would be designed, to provide data for pure compounds as well as binary mixtures,
which are currently unavailable in literature. The purpose of the new designs is to eliminate
superheat effectively yet keeping the operation simple and inexpensive.
The following chapters discuss the development of the potential model and the
ebulliometry apparatus in order to meet the goals stated above. Conclusions and
recommendations for future work are provided in the final chapter. The following paragraphs
give a brief overview of the chapters in this thesis.
Chapter II evaluates simulation-perturbation technique of SPEAD and shows systematic
errors in liquid density predictions. Citing the importance of reproducing the liquid densities
accurately, improvements to the model are discussed and then a revised model called StePPE
(Step Potentials for Phase Equilibria) is developed in Chapter III. Vapor liquid equilibrium
calculations for pure compounds with the new model are discussed in Chapter IV. Comparisons
of the results for individual homologous series, wherever possible, are made with the TraPPE,
NERD and SPEAD models. Next, the thesis discusses the development of ebulliometers for

4

vapor pressures measurements of pure compounds and bubble pressures of binary mixtures in
Chapter V. Finally, the conclusions and recommendations for future work are presented in
Chapter VI. Some of the supporting material is presented in the Appendix at the end.

5

REFERENCES

6

1.3 References
1.

Prausnitz, J.M.; Lichtenthaler, R.N.; de Azevedo, E.G., "Molecular Thermodynamics of
Fluid-Phase Equilibria". 3rd ed. 1999, Prentice-Hall, Inc.: Upper Saddle, New Jersey.

2.

Faeth, G.M., "Evaporation and Combustion of Sprays". Progress in Energy and
Combustion Science, 1983, 9, 1-76.

3.

Perry, R.H.; Green, D.W., "Perry's Chemical Engineers' Handbook". 1997.

4.

Pollara, L.Z., "Surface Tension and Vapor Pressure". J. Phys. Chem., 1942, 46, 11631167.

5.

Hala, E.; Pick, J.; Fried, V.; Vilim, O., "Vapor- Liquid Equilibrium". 1958, Pergammon
Press: New York.

6.

Rogalski, M.; Malanowski, S., "Ebulliometers Modified for the Accurate Determination
of Vapor-Liquid Equilibrium". Fluid Phase Equilibria, 1980, 5, 97-112.

7.

Chernyak, Y.; Clements, J.H., "Vapor Pressure and Liquid Heat Capacity of Alkylene
Carbonates". Journal of Chemical & Engineering Data, 2004, 49, 1180-1184.

8.

Raal, J.D.; Gadodia, V.; Ramjugernath, D.; Jalari, R., "New Developments in Differential
Ebulliometry: Experimental and Theoretical". Journal of Molecular Liquids, 2006, 125,
45 – 57.

7

CHAPTER II
2. OVERVIEW OF VAPOR-LIQUID EQUILIBRIUM CALCULATION WITH SPEAD
METHOD

2.1 Introduction
To obtain phase equilibrium properties of fluids and mixtures, the industry has relied on
semi empirical equations of state (EoS) like the Peng-Robinson[1] or SAFT[2], or macroscopic
group contribution methods, like UNIFAC[3]. Both methods become increasingly inaccurate
with increasing temperatures and pressures as well as with molecular complexity[4-7]. An
alternative is to use site-site interaction potentials (considering either an atom or a group of
atoms as a “site”), and use molecular simulations for equilibrium thermodynamic properties.
Such formalism would include well established bond lengths, angles and torsions in a molecule,
thereby incorporating more realistic structural contributions to the properties being calculated.
References [4-7] demonstrate the advantages molecular simulation models have over the
equations of state or the group contribution methods. The key, however, is the development of a
site-site potential model (or a force-field), which is both accurate, and “transferable” from
molecule to molecule. That is, the parameters optimized for a particular site-site interaction is
usable in any molecule in which that site is present.
Several UA models with 12-6 Lennard-Jones potential, like OPLS[8], TraPPE[9] and
NERD[10] have been used successfully for generating vapor-liquid coexistence curves and
critical properties. While all of these models reproduce saturation liquid densities and critical
temperatures well for the chain lengths for which they were optimized, the vapor pressure
estimation lacks accuracy. OPLS, which was optimized to reproduce saturation liquid densities
and heats of vaporization of short chain molecules at ambient conditions, overpredicts liquid

8

densities at elevated temperatures, and underestimates vapor pressures even at ambient
conditions,[9, 11] and the errors increase with increasing chain length. TraPPE predicts the liquid
densities better than OPLS for relatively longer chains (up to C10), but higher errors have been
reported for longer chains[10], where NERD seems to perform better. Vapor pressures are
consistently overpredicted with both the models, but the errors are smaller with NERD[12].
All atom (AA) simulations would be more accurate, especially at higher densities, but
their CPU times are roughly an order of magnitude higher compared to UA simulation models. A
discussion of various all atom simulations for phase behavior for n-alkanes can be found in the
study of Chen et al[13]. An intermediate solution is the anisotropic models[14], which are
computationally less demanding than the AA simulations.
The SPEAD method is another UA model which can predict vapor pressures and heats of
vaporization approximately 50 times faster than conventional full potential molecular
simulations[15]. It is a combination of discontinuous molecular dynamics (DMD) simulation and
thermodynamic perturbation theory (TPT), which gives an equation of state for VLE
calculations. SPEAD estimates vapor pressures of hydrocarbons including alkanes[16],
aromatics[17], low molecular weight ethers and primary alcohols[16], and acetates[18] with
errors of less than 10% of the experimental values for a reduced temperature range of 0.45 to 1.0.
However, as discussed later in this study, a systematic overprediction in liquid densities is
observed with increasing chain length. Notwithstanding the inadequacies, for calculations away
from the critical point, SPEAD offers a great industrial tool. If the force-field can be improved to
estimate liquid densities better, SPEAD would be a more robust alternative to equations of state
and group contribution methods on one hand, and a much faster alternative to the other UA
models with similar accuracy, on the other.

9

2.2 DMD simulations for thermodynamic properties
In molecular dynamics (MD) simulations, a system of particles that obey Newton’s Laws
of motion is created and the configurational energy of the system is given by a defined potential
function (also known as a force-field). Based on initial positions and momenta, the trajectory of
the particles can be created by integrating the Newton’s laws of motion numerically, and by
using standard statistical mechanical analysis tools, the required time averaged property of the
system can be calculated from the trajectory. The kinetic energy of the particles defines the
system temperature, and the forces on the particles are calculated using the force-field. For
Lennard- Jones type continuous potential functions, the choice of a suitable time step is
important. Also, to avoid ‘heating’ of the system during the simulation, a suitable numerical
“thermostat” has to be employed. An overview of the method can be found in standard
simulation textbooks[19, 20].
Obtaining the system trajectory is much simpler in discontinuous MD simulations. The
system acts as hard billiard balls which undergo completely elastic collisions. There is no timestep to be chosen since no velocity change occurs until a collision takes place, and since the
collisions are completely elastic, the positions and velocities of the particles post collision can be
calculated by conservation of energy and momenta considerations. Alder and Wainwright[21]
devised a method for calculating the trajectory of square well spheres. In addition to hard
collisions, well interactions were considered, and the dependency of velocities and positions on
the well depth was obtained. Rappaport[22] extended the study to chains formed by tethering the
spheres by using infinite potential wells to represent bonds. For a more in-depth analysis of the
continuous and discontinuous MD simulation methods, Maginn and Elliott’s review[23] on MD
is suggested.

10

The DMD simulation method in SPEAD develops from Rappaport’s treatment of
tethered hard sphere polymer chains[24]. The step-potentials are added post simulation using
TPT, which allows equilibrium property calculations. The time required is an order of magnitude
smaller compared to MD simulations with full potentials, as only hard interactions are accounted
for. The following sections discuss the SPEAD theory in a greater detail.

2.3 SPEAD background
SPEAD model combines DMD simulation of a “reference” system with thermodynamic
perturbation theory, which incorporates the disperse interactions. The reference fluid is
composed of hard-sphere chains, with each hard sphere representing a united atom site. For
example, the alkanes would be formed by methyl (CH3) and methylene (CH2) sites, with the
distance between each site equivalent to C-C bond length in literature. The bonds angles are
constrained by using a pseudobond between alternative sites (Fig 2.1). The bonded spheres
within a molecule overlap due to the fact that the bond length is smaller than the site diameters.
For non-bonded interactions however, the sites act as hard-spheres.

Distance (nm x 10)

4
3

2

2
1

a

0

1

θ b
c

3

-1
-2
-3

-2

0
2
Distance (nm x 10)

4

Figure 2.1: Illustration of bonds and pseudobonds for three adjacent sites in an alkane
represented as overlapping spheres. The solid lines a and b represent the covalent bonds
between the adjacent sites (1-2, 2-3) and the dashed line c represent the pseudobond
between sites 1 and 3 to constrain the angle θ
11

1. SPEAD force- field
Infinite potential wells are used to restrict bond lengths and bond angles within a
molecule (Fig 2.2). The innermost infinite well represents a covalent bond, and the second
infinite well is used for a pseudobond between first and third neighbors in order to restrict the
bond angles. Both the wells are centered corresponding to realistic equilibrium values of bond
lengths and angles, and allowed to vary within 5% of those values. For example, for a C-C bond
would be represented by a well centered at 0.154nm and a C-C-C bond angle would be depicted
by a well centered at a value corresponding to ~110 degrees. A third well is used for additional
bonds for 1-4 interactions to constrain cis/trans configuration. The center of the well for the
cis/trans configuration is obtained from simple geometry (Fig 2.3). The non-bonded interactions
during the simulation are repulsive only.

Figure 2.2: Illustration of the multistep potential used in SPEAD. The innermost well
represents a covalent bond and the second inner well represents pseudobond to constrain
bond angles
12

1.34 Å
1.54 Å

120°

120°
1.54 Å

1.54 Å

1.54 Å
120°

1.34 Å

Figure 2.3: Finding the pseudobond length (given by the dashed lines) for cis/trans
interactions in cis/trans butene

The effect of the non-bonded disperse interactions is added post-simulation as
perturbation. A step-wise approximation to the Lennard-Jones model is used, which mimics the
-6

r behavior of potential expected theoretically (Fig 2.2). Such an approximation can provide
properties that are similar to the continuous potential for spheres[25]. While an infinite number
of small steps may be the best approximation of the continuous potential, recent efforts focus on
optimization of a four-step potential. Following the reference simulation, the inner and outer well
depths corresponding to the four steps are adjusted, and the two central wells are interpolated.
The four step-wells are located at r/σ = 1.2, 1.5, 1.8, 2.0, where σ is the site diameter, as shown
in Figure 2.2. Hydrogen bonding is represented as an additional perturbation using Wertheim’s
theory as discussed in the next section.

2. Reference Simulation
The reference system is composed of hard sphere chains with bond lengths and bond
angles constrained at literature values as described in the previous section. First the core volume
(volume occupied by one molecule, b) is found by simulating a single molecule. For n number of
unit cells, 4n molecules are then arranged in an fcc lattice and the simulation is started. For
molecules with more than 20 sites, 256 molecules are used and for smaller molecules, 108

13

molecules are used for the simulation. The simulation box is further divided into a number of
cells for event tracking, such that for a particular site, only the sites in the adjacent cells are
considered for an event. The events that are kept track of are bond “collisions”, cell crossings,
and site-site collisions. For efficiency, the algorithm for tracking the events is based on
Rappaport’s event tree formalism[26, 27]. The system equilibrates for a given number of
intermolecular hard collisions, before the production run starts when data is collected. Typically,
the number of collisions for equilibration is set to about 500,000. The total number of collisions
range from about 6000,000 to 10,000,000, the higher number generally being used for bigger
molecules. The actual simulation time ranges between a few hundred to a few thousand
picoseconds depending on the size of the molecule and density of the system. Low density
(packing factors < 0.01) runs with small molecules take the highest amount of wall time to run.
For example, the CPU time using a 64 bit processor for a molecule with ten sites is about two
hours at packing fractions of 0.001 and about one hour at a packing fraction of 0.56. During the
simulation, there is no fluctuation in the velocities since the attractive interactions are not present
to accelerate the particles. The velocities are checked periodically to confirm the temperature
remains constant. The simulation is repeated for about 21 dimensionless packing fractions
(represented by ‘η’ which is equal to bρ, where b is the molecular core volume and ρ is the
density corresponding to the reciprocal units of b) ranging from vapor-like to liquid- like
densities (bρ ranging from 0.001 to 0.54). During the simulation, the number counts for pair
types in well limits are collected and averaged in order to later apply the TPT, as discussed in the
section below. Typically 300 samples are averaged over the course of the run. The parameters
for interactions of unlike sites are represented by Lorenz-Berthelot combining rules,

σ ij =

1
(σ i + σ j ) and ε ij = (ε i ε j )12 , where i and j are different site types. While σ ij
2
14

is

used in the reference simulation for the hard-site collisions,

ε ij is used during the application of

TPT, where the attractive wells are introduced.
Pair types are distinct from pairs in following discussions. For example, consider a
system of propane molecules. Using the united atom approach, the propane molecule is
represented by two CH3 groups connected by a CH2 group. When two propane molecules
interact, there are nine site pairs that interact. However, there are only three pair types (CH3 +
CH3, CH2 + CH2, CH2 + CH3). When accumulating interactions, the two approaches share
many similarities, and can be proven to be equivalent. In the following discussion, sums are
written over pair types.
3. Perturbation theory in SPEAD
TPT in SPEAD is based on Barker and Henderson’s formalism[28], which builds on the
commonly known high temperature series expansion of the configurational partition function
developed by Zwanzig[29]. The basic idea is to divide the potential energy of the system of
interest into a reference (unperturbed) and a perturbed part, and mathematically express the
Helmholtz free energy as a series expansion in inverse temperature. The different terms in the
expansion can be found using ensemble averages over the reference system. The derivation is
shown below.
From statistical mechanics:
A = −kT ln Q

(1)

where, k is the Boltzmann constant, A is the extensive Helmholtz energy and Q is the canonical
partition function, or the configurational integral. Zwanzig divided the potential energy of the
system to an unperturbed part (Uo) and a perturbation contribution (UP), as:

15

U = U0 + UP

(2)

This facilitates writing the configurational integral in (1) as:

 U 
Q = Q0 exp − P 
 kT 

(3)
O

The angular brackets denote an ensemble average and the subscript “o” implies an average over
the unperturbed system (dropped in following equations for simplicity) and

Q0 denotes the

configurational integral of the unperturbed system.
For an ensemble of monatomic molecules with a single square well on each site, Barker
and Henderson simplified equation (3) to:

 1
Q = Q0 exp −
 kT

∑
i


 1
N i ε i  exp −

 kT

In equation (4), “i” denotes a pair type,

Ni

εi

∑ (N
i

i


− Ni ) i 
ε


(4)

denotes the potential well depth for pair type i and

is the number of pairs of type i within the well, as determined by evaluating snapshots of

the simulation. Expanding the exponentials in equation (4) and taking the logarithm (see
equation (1)), we get:

A
A
1
= 0 +
NkT NkT kT
where,

∑
i

Ni εi −

∑∑ ( N N
2(kT )
1

i

2

i

j

− N i N i ε iε j

)

(5)

j

A0 is the Helmholtz energy of the reference system, N is the total number of molecules

and i and j denote pair types. Extending the analysis to four step wells, the Helmholtz free energy
departure can be expressed as,

16

A − Aig A0 − Aig
=
+
NkT
NkT
1
1
N im ε im −
∑∑
NkT i m
2 N (kT ) 2

∑∑∑∑ ( N
i

j

m

im

)

N jn − N im N jn ε imε jn

n

(6)
In equation (6) the subscript i, j are for pair types, while m, n are for wells, such that

N im

would be the average number of pairs of type i in well m. For example, in alkanes, there are three
pair types, so i and j would vary from 1 to 3 and m and n would vary from 1 to 4 for the wells.
The well energies for unlike sites in a pair type (the CH3-CH2 pair in alkanes) would be given
by the Lorenz-Berthelot combining rule, as stated above. An excellent review on the progress of
TPT for thermodynamic properties is presented in the review by Zhou and Solana[30].
The second and third terms on the right-hand side of the equation are the commonly
known first and second order perturbation contributions to the Helmholtz free energy. When
hydrogen bonding is present, an association contribution to the Helmholtz energy is added to (6).
The first and second order perturbation sums are abbreviated by A1 and A2 (except T dependence
is kept explicit), and the hydrogen bonding perturbation is added so that (6) can be written as:

A0 − A ig A1 A2 A assoc
A − A ig
=
+
+
+
NkT
NkT
T T 2 NkT

(7)

The term A1 represents the average configurational energy of the system due to
perturbation, while A2 represents the fluctuation of that energy around the mean value. The
A

assoc

term is included if hydrogen-bonding sites are present. Care is taken in implementation of

17

the sums to avoid double-counting in the programming of the sums. From classical
thermodynamics:

( A − A ig ) T ,V
RT

ρ

=∫

( Z − 1)dρ

ρ

0

(8)

Therefore, the corresponding equation of state (Z = f(P,V,T)) can be represented as:

Z = 1+
Since η = bρ ,


ρ d
( A − A ig ) T ,V 
RT  dρ



(9)

ρ
η
. Therefore:
=
dρ d η

 d ( A0 − Aig ) 
dA 1
dA
 +η 1   +η 2
Z = 1 + η
 


dη
dη  T 
dη



 1 
Assoc
 2 +Z
T 

(10)
or introducing shorthand variables:

1
 1
Z = Z 0 + Z1   + Z 2  2
T 
T


Assoc
+Z


(11)

The term Z0 in equation (11) denotes the compressibility factor of the unperturbed fluid, which
is used to calculate the Helmholtz energy departure of the reference fluid. From (10), one can
deduce:

(

 d A0 − A ig
Z 0 = 1 + η

dη


)  ⇒ (Z




0

− 1)

18

dη

η

(

= d A0 − A ig

)

(12)

Therefore the Z0 terms are collected from the simulations at all the packing fractions and a
polynomial function is fit to the values as shown below, which can be used in equation (12) to
calculate the Helmholtz energy departure

(A

0

− Aig

) of the reference fluid.

1 + a1η + a2η 2 + a3η 3
Z0 =
(1 − η )3

(13)

To represent hydrogen bonding, Wertheim’s theory[31] is used and the contribution of
association has to be included in the Helmholtz energy and the equation of state. For a molecule
with M different associating sites, the contribution of association to the Helmholtz energy is
given by[32]:


A assoc
Xi 1
i
= ∑ ln X −
+ M
NkT
2  2
i 
where,

Xi

(14)

is the mole fraction of hydrogen bonding sites that are not bonded. For one donor

and one acceptor site, M = 2, so that the above equation can be written as,

A assoc
XA
XD
A
D
= ln X −
+ ln X −
+1
NkT
2
2

(15)

The superscripts A and D denote acceptor and donor in the above equation. Furthermore, since
the current implementation of SPEAD does not differentiate between donors and
acceptors,

X D = X A , therefore, equation (15) can be simplified as:
A assoc
= 2 ln( X A ) + (1 − X A )
NkT

(16)

The corresponding contribution to the compressibility factor can be given by:

19

Z

assoc

 1 + η −η 2 / 2 
= −(1 − X A )
 (1 − η )(1 − η / 2) 




(17)

A

where, η is the packing fraction, and X is calculated by:

X

A

=

− 1 + 1 + 4α
2α

(18)

In the above equation, α , and is given by[33]:

α = ρ∆AD

(19)

In equation (19), ρ is the molar density and

∆AD

is the “strength” of association as it includes

the range and depth of the association square well and is dependent upon the hard-sphere radial
hs

distribution function (g(r) )[2, 33]:

(

)

∆ AD = g (r ) hs K AD exp(ε AD / kT ) − 1
In the above equation, K
AD

association, and ε

AD

(20)

is the bond volume, depicted by the width of the square well for

is the bond energy. Using equations (19) and (20) and the Carnahan and

Starling’s equation for hard sphere radial distribution function:

g (r ) hs =

(1 − η / 2)

(1 − η )3

(21)

we get:

α=

ρ (1 − η / 2) AD
K [exp(ε AD / kT ) − 1]
3
(1 − η )

The bond volume is calculated from an empirical correlation:

20

(22)

K

AD


n
= κ 1 + λ  hb
 n
 t







(23)

where nhb is the number of hydrogen-bonding sites in the molecule, nt is total number of sites,
and κ and λ are empirical constants for a given type of hydrogen bonding moiety. The empirical
relation (23) causes the perturbation to be smaller in larger molecules.

4. Summary of steps in the SPEAD method
Once the compressibility factor for the reference fluid (Z0) is calculated from the
reference simulation, TPT can be applied to include the contributions from the first and second
order perturbation terms as well as association, to complete the calculation of the compressibility
factor for the intended fluid. A flowsheet of the steps is shown below.

Simulate repulsive chain at specified packing fractions,
Collect radial distribution functions

Reference
Simulation

Apply well depths (ε1, ε4) to calculate A1, A2, to get
simulated Helmholtz energies at each packing fraction

Including
TPT

Fit A1, A2 ‘data’ with analytical functions

 d ( A0 − Aig ) 
dA  1 
dA  1 
 +η 1   +η 2  2 
Z = 1+η


dη
dη  T 
dη  T 


Fig 2.4: SPEAD simulation steps to get the EoS

21

Fitting

EoS

5. Vapor pressure and orthobaric densities from SPEAD EoS
To perform the vapor pressure calculations, the compressibility factor is used to
determine the Gibbs energy departure (a measure of fugacity). Solving for densities that provide
equal Gibbs departures for the vapor and liquid phases gives the vapor pressure for the given
temperature, and the densities at the saturation pressure correspond to the equilibrium liquid and
vapor densities. The temperatures for the vapor pressure and density calculations are obtained
from a data file that includes the experimental vapor pressures and liquid densities for the given
compound. The algorithm is shown in the flowsheet below.

Start

Input T

sat
r

log P


7
1
= (1 + ω )1 −
 T
3
r







End

P

sat

Guess
new
sat
P

η V = bP sat / RT
η L = 0.65

P calc (η L ,η V ) = ZRTη / b

sat

=P

Yes

No
L V

f /f -1 < 0.001

Guess new η
No

P

calc

sat

/P

Yes

-1< 0.01

f (η L ,η V ) = P sat exp( A + Z − 1 − log Z )

η V = η V ,η L = η L
Fig 2.5: Flowsheet of equilibrium pressure and density calculations. Symbol meanings- Pr :
Reduced pressure, Tr : Reduced temperature, ω: Compressibility factor, b : Core volume
(volume/molecule), η : Packing fraction (V: vapor, L: liquid), f : Fugacity, A : Helmholtz
energy departure

22

To obtain the optimum values of the step wells, vapor pressure and density calculations
are performed by changing the well parameters globally of a selected set of compounds (training
set) for which experimental data is available. The well-parameters that give the minimum global
error for the training set give the optimized values of the well energies. The parameters can then
be used for calculating the densities and vapor pressures of an unknown system.

23

REFERENCES

24

2.4 References
1.

Prausnitz, J.M.; Lichtenthaler, R.N.; de Azevedo, E.G., "Molecular Thermodynamics of
Fluid-Phase Equilibria". 3rd ed. 1999, Prentice-Hall, Inc.: Upper Saddle, New Jersey.

2.

Chapman, W.G.; Gubbins, K.E.; Jackson, G.; Radosz, M., "New Reference Equation of
State for Associating Liquids". Industrial & Engineering Chemistry Research, 1990, 29,
1709-1721.

3.

Fredenslund, A.; Jones, R.L.; Prausnitz, J.M., "Group-Contribution Estimation of
Activity Coefficients in Nonideal Liquid Mixtures". AIChE Journal, 1975, 21, 10861099.

4.

Spyriouni, T.; Economou, I.G.; Theodorou, D.N., "Phase Equilibria of Mixtures
Containing Chain Molecules Predicted Through a Novel Simulation Scheme". Physical
Review Letters, 1998, 80, 4466-4469.

5.

Nath, S.K.; Escobedo, F.A.; de Pablo, J.J.; Patramai, I., "Simulation of Vapor-Liquid
Equilibria for Alkane Mixtures". Industrial & Engineering Chemistry Research 1998, 37,
3195-3202.

6.

Bureau, N.; Jose, J.; Mokbel, I.; de Hemptinne, J.-C., "Vapour Pressure Measurements
and Prediction for Heavy Esters". Journal of Chemical Thermodynanmics, 2001, 33,
1485-1498.

7.

Asher, W.E.; Pankow, J.F.; Erdakos, G.B.; Seinfeld, J.H., "Estimating the Vapor
Pressures of Multi-Functional Oxygen-Containing Organic Compounds Using Group
Contribution Methods". Atmospheric Environment, 2002, 36, 1483–1498.

8.

Jorgensen, W.L.; Madura, J.D.; Swenson, C.J., "Optimized Intermolecular Potential
Functions for Liquid Hydrocarbons". Journal of the American Chemical Society, 1984,
106, 6638-6646.

9.

Martin, M.G.; Siepmann, J.I., "Transferable Potentials for Phase Equilibria. 1. UnitedAtom Description of n-Alkanes". Journal of Physical Chemistry B, 1998, 102, 25692577.

10.

Nath, S.K.; Escobedo , F.A.; de Pablo, J.J., "On the Simulation of Vapor–Liquid
Equilibria for Alkanes". Journal of Chemical Physics, 1998, 108, 9905-9911.

25

11.

Delhommelle, J.; Tschirwitz, C.; Ungerer, P.; Granucci, G.; Millie´, P.; Pattou, D.; Fuchs,
A.H., Journal of Physical Chemistry B, 2000, 104, 4745-4753.

12.

Errington, J.R.; Panagiotopoulos, A.Z., "A New Intermolecular Potential Model for the nAlkane Homologous Series". Journal of Physical Chemistry, 1999, 103, 6314-6322.

13.

Chen, B.; Martin, M.G.; Siepmann, J.I., "Thermodynamic Properties of the Williams,
OPLS-AA, and MMFF94 All-Atom Force Fields for Normal Alkanes". The Journal of
Physical Chemistry, 1998, 102, 2578-2586.

14.

Toxvaerd, S., "Molecular Dynamics Calculation of the Equation of State of Alkanes".
Journal of Chemical Physics, 1990, 93, 4290-4295.

15.

Case, F.; Chaka, A.; Friend, D.G.; Frurip, D.; Golab, J.; Gordon, P.; Johnson, R.; Kolar,
P.; Moore, J.; Mountain, R.D.; Olson, J.; Ross, R.; Schiller, M., "The second industrial
fluid properties simulation challenge". Fluid Phase Equilibria, 2005, 236, 1-14.

16.

Unlu, O.; Gray, N.H.; Gerek, Z.N.; Elliott, J.R., "Transferable step potentials for the
straight-chain alkanes, alkenes, alkynes, ethers, and alcohols". Ind. Eng. Chem. Res.,
2004, 43, 1788-1793.

17.

Gray, N.H.; Gerek, Z.N.; Elliott, J.R., "Molecular Modeling of Isomer Effects in
Naphthenic and Aromatic Hydrocarbons". Fluid Phase Equilibria, 2005, 228, 147-153.

18.

Baskaya, F.S.; Gray, N.H.; Gerek, Z.N.; Elliott, J.R., "Transferable Step Potentials for
Amines, Amides, Acetates, and Ketones". Fluid Phase Equilibria, 2005, 236 42–52.

19.

Frenkel, D.; Smit, B., "Understanding Molecular Simulation". 2nd ed. 1996, Academic
Press: San Diego.

20.

Allen, P.; Tildesley, D.J., "Computer Simulation of Liquid". 1987, Oxford Science
Publication: Oxford.

21.

Alder, B.J.; Wainwright, T.E., "Studies in Molecular Dynamics. I. General method".
Journal of the American Chemical Physics, 1959, 31, 459- 436.

22.

Rappaport, D.C., "Molecular Dynamics Simulation of Polymer Chains with Excluded
Volume". Journal of Physics A: Mathematical and General, 1978, 11, L213-L217.

26

23.

Maginn, E.J.; Elliott, J.R., "Historical Perspective and Current Outlook for Molecular
Dynamics As a Chemical Engineering Tool". Industrial & Engineering Chemistry
Research 2010, 49, 3059-3078.

24.

Rapaport, D.C., "Molecular Dynamics Study of a Polymer Chain in Solution ". Journal
of Chemical Physics, 1979, 71, 3299-3303.

25.

Chapela, G.A.; Scriven, L.E.; Davis, H.T., "Molecular Dynamics for Discontinuous
Potential. IV. Lennard Jonesium". J. Chem. Phys., 1989, 91, 4307.

26.

Rapaport, D.C., "The Event Scheduling Problem in Molecular Dynamic Simulation".
Journal of Computational Physics, 1980, 34, 184-201.

27.

Rappaport, D.C., "The Art of Molecular Dynamics Simulation". 2nd ed. 2004,
Cambridge University Press: Cambridge, UK.

28.

Barker, J.A.; Henderson, D., "Perturbation Theory and Equation of State for Fluids: The
Square-Well Potential". The Journal of Chemical Physics, 1967, 47, 2856-2861.

29.

Zwanzig, R.W., "High- Temperature Equation of State by a Perturbation Method. I.
Nonpolar Gases". The Journal Of Chemical Physics, 1954, 22, 1420-1426.

30.

Zhou, S.; Solana, J.R., "Progress in the Perturbation Approach in Fluid and FluidRelated Theories". Chemical Reviews, 2009, 109, 2829-2858.

31.

Wertheim, M.S., "Fluids with highly directional attractive forces. II. Thermodynamic
perturbation theory and integral equations". J. Stat. Phys., 1984, 35, 35-47.

32.

Wertheim, M.S., "Fluids with Highly Directional Attractive Forces. IV. Equilibrium
Polymerization". Jornal of Statistical Physics, 1986, 42, 477-492.

33.

Elliott, J.R., "Efficient Implementation of Wertheim's Theory for Multicomponent
Mixtures of Polysegented Species". Industrial & Engineering Chemistry Research, 1996,
35, 1624-1629.

27

CHAPTER III
3. PREDICTION OF LIQUID DENSITIES AND VAPOR PRESSURES WITH SPEAD:
MOTIVATION FOR A NEW STEP POTENTIAL MODEL

3.1 Summary of results with SPEAD: model parameters and the liquid density errors
The SPEAD force-field is a major improvement in comparison to TraPPE for vapor
pressure predictions away from the critical point, especially at lower reduced temperatures[1].
For the classes of compounds shown in Figure 3.1, the average absolute deviations (AAD) in
vapor pressure predictions with SPEAD are within 16% of the experimental data.

18
16
14

Vapor pressure
Liquid density

%AAD

12
10
8
6
4
2
0
alkanes

alcohols

ethers

acetates

ketones

Amines

Figure 3.1: Published results of vapor pressure and liquid densities using SPEAD

However, this improvement comes at the expense of accuracy in the liquid density
predictions. Except for the ethers, the density deviations of all the other compounds in Figure 3.1
are approximately between 2.5- 6%. It is desirable that the density deviations fall within 1-2% of
the experimental values, since small changes in densities can greatly affect other thermodynamic

28

properties. Figure 3.2 below shows the simulated liquid densities of alkanes with the SPEAD
model. With increasing molecular weight, the densities are increasingly overpredicted. A similar
trend is observed in the density errors of ketones, acetates, α-olefins, ketones and acetates (Fig
3.6). This trend is not so clear in some of the compounds which have been optimized only for
smaller chains, for example, ethers[1]. Intuitively, this could be a result of “overoptimization” of
the non-alkane parameters in these compounds, in order to compensate for the errors being
carried over from the alkane parameters. These errors in the densities are observed despite the
use of extra parameters for some of the functional groups, like alcohols[1], ketones[2] and
acetates[2].

700

600

T(K)

500
400
300
200
100

0

0.1

0.2

0.3

0.5
0.4
ρ (g/mL)

0.6

0.7

0.8

0.9

Figure 3.2: Density predictions of normal alkanes with Akron SPEAD. The dots are
experimental points and the lines are model predictions

29

Since thermodynamic properties are sensitive to small changes in densities, it is
imperative for any new potential to reproduce densities well. Care also needs to be taken to
ensure transferability of the parameters, which can be achieved by including a broader range of
molecular weights in order to incorporate the smallest possible contribution from a functional
group.

3.2 Improvements: the StePPE model
To remedy the discrepancies in SPEAD discussed above, some of the model parameters,
such as site sizes and intramolecular interaction parameters, needed readjustment in order to
study their effect on the vapor pressure and density predictions. The site sizes can be changed
manually, but to change any other simulation parameter, access to the source code is required.
The SPEAD software is a compiled version instead, owned by ChemStations Inc. of Houston,
Texas, and the source code is unavailable. Therefore, with the aim of developing a robust step
potential model, using DMD and TPT, which could be used as a fast tool in a variety of
industrial applications, a new source code was developed. To avoid confusion with the SPEAD
model and parameters, we refer to this new model as StePPE, which stands for Step Potentials
for Phase Equilibria, though it uses much of the same methodology as SPEAD. As the name
suggests, the form of the force-field is the same (see “SPEAD theory” in the introduction part),
but some changes were necessary in order to reproduce the liquid densities more accurately,
without compromising the speed or accuracy in the vapor pressure predictions observed with the
SPEAD model. The changes are discussed in the following sections.

30

1. Hard-site diameters from excluded volumes
Intuitively, the density errors in alkanes can be addressed by revisiting the hard-site
diameters used for the methyl and methylene groups, since the errors in other compounds may be
due to poor transferability of the alkane parameters. The site diameters in SPEAD were found by
fitting the experimental vapor pressures. In the following paragraphs, we discuss our method of
site-size selection using excluded volume considerations.
Hard site diameters can be found by mapping of the explicit atom potentials to UA
potentials and then mapping it to a hard site potential, as described by McCoy[3]. A coarser
method is to use a combination of the tabulated excluded volumes of the individual atoms in a
particular group, and find the diameter of an equivalent sphere. Both methods have been shown
to be in close agreement for the diameter of methane[4]. If the volume lost due to overlaps is
accounted for correctly, the latter method can work reasonably well for other sites as well.
In our model, we use the core volumes of groups, as tabulated by Bondi[5]. Bondi
determines the group core volumes by summing the contributions of spheres for the individual
atoms. In this work, we find an equivalent sphere to represent the group core volumes as a single
united atom group sphere and not spheres for the individual atoms, and we then calculate the
united atom site diameter. Volume lost due to overlap by each covalent bond is calculated by
simple geometry (see Appendix A.2) and added to the excluded volume of the group. The
diameter of an equivalent sphere with this corrected volume gives the first guess of the site
diameter. The volume correction would obviously depend on the adjacent sites and bond lengths,
which would be tedious to tabulate for every site. Just for the quaternary carbon in branched
alkanes, a total of 24 different sizes would have to be tabulated!
We calculated the volume corrections using all possible site combinations in an n-alkane,
namely, CH3-CH3, CH2-CH2 and CH3-CH2. Using an analysis based on the bond lengths and
31

the excluded volumes as described by Bondi[5] and Slonimskii et al[6], we found that the
volume lost for each covalent bond in each case was close to 20%. We used that correction for
sites in all other functional groups to provide a first guess at the site diameters. For example, the
first guess for the diameter of the quaternary carbon was obtained from adding an 80%
correction, 20% for each covalent bond. We refer to the diameters calculated by the
approximation as σvdw. Starting from this diameter, we found the optimized diameters to be
slightly smaller than σvdw. Figure 3.3 (a) and (b) show the site sizes based on our analysis and
the optimized sizes in StePPE and TraPPE respectively. In both cases, the final optimized site
3

sizes follow the sizes obtained from van der Waal’s volumes closely, the sp CH in TraPPE
being the only outlier that doesn’t follow that trend. The sensitivity of the site diameters on the
3

volume correction increases with the number of covalent bonds. For example, the sp C in the
branched alkanes has four covalent bonds, so will be affected most by a change in the volume
3

correction. The methyl group would be least sensitive. The unusually big size of the sp C in
TraPPE (6.4 Å) can be explained based on this reasoning.
The site sizes in the UA models in general, are obtained from fitting available
experimental thermodynamic data, such as vapor-liquid coexistence properties or critical
properties. In the continuous models, sigma and epsilon are optimized simultaneously. Changing
either parameter would require the simulation to be run again. In SPEAD, the well parameters
are optimized post simulation, but changing the site diameter would require the simulation to be
repeated. Therefore, whether it is the computationally more expensive continuous models, or the
four step TPT model, a methodology to get a good first guess for the size diameter would be
extremely helpful in minimizing the CPU time for optimization.

32

3.6

7.0

(a)

>C=
CH2=
CH
CH3
CH2
C
-CH=

>C=
CH2=
CH
CH3
CH2
C
-CH=

6.0

σOptimized

σOptimized

3.8

5.0

(b)

3.4
4.0

σOptimized =

0.97
3.0

3.2
3.2

3.4

3.6

3.8

3.0

4.0

5.0

6.0

σvdW

σvdW

Figure 3.3: Optimized site diameters (ordinate) compared to diameters obtained from van
der Waal's volumes (abscissa) in Å; (a) StePPE, based on 20% volume correction per bond,
(b) TraPPE, based on 24% volume correction per bond

The molecular core volumes are determined during a run. Prior to a run, the sites are
placed at the bond distances given in Table A2.1 in the appendix. The core volume of the free
molecule is determined by then vibrating the molecule for 3000-5000 ps and taking 100
snapshots. To prevent translation of the center of mass, the first and last sites in the structure file
are accelerated towards each other. The degrees of freedom are reduced by three and when
relating the kinetic energy to temperature, the simulation is started at 300K. After an
equilibration time (10000 + 1000×Number of molecules per site), 100 snapshots are taken at
intervals of 30 - 60 ps depending on the size of the molecule. For each snapshot, a grid of points
placed 1/20 Å apart is placed over the molecule. The core volume of the molecule is determined
by calculating the fraction of points inside the molecule, which corresponds to the fraction of the
cell occupied by the core volume. Different spacings of grid points were tested, and the spacing
of 1/20 Å was determined to be suitable. Due to flexibility of the molecule, the core volume

33

7.0

varies amongst the snapshots. Use of a finer grid took a longer time, and did not result in
significantly different core volumes. Further, the free molecule core volume is slightly different
than the core volume determined during a simulation. During the simulation, collisions with
other molecules affect the core volume and the trend is that the core volume decreases slightly
(0.2 %) at high density when the molecules are tightly packed. When each density is run, the
core volumes of the first 100 molecules in the simulation box are sampled at the end of the run,
and the average core volume is determined. Then the core volumes from the series of densities
simulated are averaged to determine the reported core volume. As an example of the results for
the core volumes, the core volumes of n-alkanes are plotted against van der Waals volumes
tabulated in DIPPR in Figure 3.4. The core volumes in StePPE are smaller than the van der
Waal’s volumes since the optimized sizes of the CH3 and the CH2 groups are slightly smaller
than that predicted by the excluded volumes.

vdw Volume (cm3/mol)

300
vdW vol Lit

250

Core vol StePPE

200
150
100
50
0
0

10

20

30

Nc
Figure 3.4: Core volumes predicted with StePPE plotted against van der Waal’s volumes of
alkanes. The x-axis denotes the number of C atoms in the alkane chain
34

2. Exponential function for A2
The equation of state is obtained from derivatives of the first and second order
perturbation terms. Therefore, a good fit to the simulation data for the perturbation terms is
essential for the new equation of state to be accurate. For the second order perturbation term, we
use a new exponential function given by:

A2 = ( a21η + a22 ) exp(− a23η a24 ) − a25
where, a2is are constants to be found from the best fit. It doesn’t differ much from the published
SPEAD A2 polynomial function for lower and packing factors near 0.35 – 0.4, but in some
cases, fits better at moderate densities (Figure 3.5).

20000
0

A2

-20000
-40000
Simulation 'data'
-60000

Polynomial Function
Exponential Function

-80000
-100000
0

0.2

0.4

0.6

0.8

η
Figure 3.5: Comparison of A2 using polynomial function as in Akron- SPEAD and
exponential function as in MSU-SPEAD for trioxane

35

1

Also, the fit with the exponential function would be more reasonable if the runs were to be
extended to higher densities, though not expected at the current level of detail (note the
decreasing A2 in Figure 3.5 for the polynomial function). Overall, the function is used to
interpolate and differentiate the simulation results; we don’t expect a dramatic change in the
overall accuracy of the model with this change, but the exponential function would give a more
realistic fit overall. The function was under consideration since we had observed an unrealistic
“dip” in the prediction of A2 at higher densities seen in the earliest SPEAD studies[7]. Though
the functional form for A2 has changed in SPEAD too, we deemed it reasonable to use the
exponential function in our model.

3. Intramolecular interaction parameters
In SPEAD, any two sites separated by less than 6 bonds interact with special reduced
diameters (σIntra). The main purpose of the intra diameters is to let sites flip past each other so as
to make the dihedral angle flexible, which in turn imparts flexibility to the chain. In our model,
we use intra diameters only for a 1, 4 interactions. Also, the pairwise attractive potential starts at
a separation of 3 bonds in our model, as opposed to 7 bonds in SPEAD. The purpose of this
change was to incorporate all the intramolecular attractions possible. Bond lengths are taken
largely from the NIST Computational Chemistry Comparison and Benchmark DataBase[8]
Average values are used for each class of bond. No special bond lengths were used for specific
molecules. In SPEAD, a large table of bond lengths is used to characterize the second-neighbor
distances for constraining bond angles. In StePPE, we have instead specified bond angles for
3

2

given site types. For example, bond angles for sp molecules are typically 109°, but for sp they

36

are set to 120°. The lengths of pseudobonds are calculated at the start of a simulation run, but
matching the bond lengths with the site angles. The pseudobonds are allowed to vibrate ±5%.

3.3 Summary of comparison of errors with SPEAD and StePPE
The effect of the improvements discussed above is manifest in the results obtained with
StePPE. Not only is the systematic error in the densities removed (Fig 3.6), we have been able to
eliminate a significant number of parameters for a number of functional groups, as will be
discussed in the third chapter. An overall comparison of errors in vapor pressures and liquid
densities with SPEAD and StePPE is shown in Table 3.1

Table 3.1: Summary of errors in vapor pressure and liquid density predictions with
SPEAD and StePPE
sat

Compound
n-Alkanes
n-Alkenes
n-Alcohols
Branched alkanes
Secondary alcohols
Ethers
Ketones
Esters

P Error
(%AAD)
SPEAD StePPE
4.54
6.15
7.19
10.35
9.01
8.13
9.55
9.39
5.4
8.56
12.2
8.1
13.6
10.87

ρ Error
(%AAD)
SPEAD StePPE
3.1
0.76
1.61
0.84
2.49
0.65
1.1
0.58
1.83
0.49
3.4
0.37
5.5
0.93

Nc range
(#compounds simulated)
SPEAD
StePPE
1 – 16 (12)
2 – 24 (16)
3 – 8 (20)
3 – 10 (15)
1 – 10 (10)
1 – 20 (17)
3 – 10 (24)
3 – 9 (10)
2 – 6 (6)
2 – 18 (14)
3 – 9 (15)
3 – 9 (16)
3 – 12 (9)
3 – 19 (25)

Some other aspects have been the subject of recent studies, where we have tried to
analyze the effect of tail corrections and include the effect of electrostatic charges in the model.
These topics are discussed in the “Future work” section in CHAPTER VI, where the
implementation procedure and some preliminary results are discussed.

37

% Error bias

10

5

0

-5
0

5
10
15
# Carbon atoms

20

Figure 3.6: Liquid density deviations of selected compounds with SPEAD (filled symbols).
Symbols are: diamonds- alkanes[9], circles- 1-alkenes[1], triangles- ketones[2], squaresacetates[2]. The open symbols are with StePPE

38

REFERENCES

39

3.4 References
1.

Unlu, O.; Gray, N.H.; Gerek, Z.N.; Elliott, J.R., "Transferable Step Potentials for the
Straight-Chain Alkanes, Alkenes, Alkynes, Ethers, and Alcohols". Industrial &
Engineering Chemistry Research 2004, 43, 1788-1793.

2.

Baskaya, F.S.; Gray, N.H.; Gerek, Z.N.; Elliott, J.R., "Transferable Step Potentials for
Amines, Amides, Acetates, and Ketones". Fluid Phase Equilibria, 2005, 236 42–52.

3.

McCoy, J.D.; Curro, G.J., "Mapping Explicit Atom onto United Atom Potentials".
Macromolecules, 1998, 31, 9362-9368.

4.

McCoy, J.D.; Mateas, S.; Zorlu, M.; Curro, J.G., "The Role of Excluded Volume in
Polyethylene Intermolecular Potentials". Journal of Chemical Physics, 1995, 102.

5.

Bondi, A., "van der Waal's Volumes and Radii". The Journal of Physical Chemistry,
1964, 68, 441-451.

6.

Slonimskii, G.L.; Askadskii, A.A.; Kitaigordoskii, A.I., "The Packing of Polymer
Molecules". Polymer Science USSR, 1970, 12, 556-577.

7.

Cui, J.Y.; Elliott, J.R., "Phase Diagrams for a Multistep Potential Model of n-Alkanes by
Discontinuous Molecular Dynamics and Thermodynamic Perturbation Theory". Journal
of Chemical Physics, 2002, 116, 8625-8631.

8.

http://cccbdb.nist.gov/.

9.

Elliott, J.R.; Gray, N.H., "Asymptotic Trends in Thermodynamic Perturbation Theory".
Journal of Chemical Physics, 2005, 123, 184902-1-184902-6.

40

CHAPTER IV
4. LIQUID DENSITIES AND VAPOR PRESSURES OF PURE COMPOUNDS WITH
StePPE MODEL

4.1 Site diameter and step potential optimization: summary
Liquid density and vapor pressure data from DIPPR 801 were used to optimize the site
parameters in this work. Comparisons with SPEAD[1], TraPPE[2] and NERD[3] models are
made, wherever published data was found with those models. Experimental data obtained from
the DIPPR database was used for the optimization, but for simplicity, when comparisons are
made with other models, we use the correlations fit to the accepted experimental data available in
DIPPR, except when the comparison is only with SPEAD. Only the experimental data that was
categorized as “accepted” in the DIPPR database was used, and smoothed data was avoided. The
StePPE fits for complete set of molecules in a homologous series are shown in separate plots in
Appendix A.3.
Three parameters were optimized for each site, namely the hard site diameter (σ) and the
inner (ε1) and the outer (ε4) well depths using the HEEDS optimization software. The
optimization method used was proprietary method called SHERPA[4], which stands for
Simultaneous Hybrid Exploration that is Robust, Progressive and Adaptive. It is a hybrid method
which uses multiple search methods simultaneously using the best attributes of each method. It
uses a combination of global and local search methods, and the parameters of the methods used
are tuned internally. The number of evaluations to perform is defined by the user. Typically, the
SHERPA optimization was followed by a quadratic gradient search to further optimize the
parameters.

41

Table 4.1: Summary of optimized parameters for sites used in this study. The sites within
parentheses in the first columns are the sites to which the group in that column is bonded,
or a special case. For example ‘CH3 (-CH2-/-CH3)’ represents a methyl group, which is
bonded to another methyl or methylene group. The ‘>’ or ‘<’ signs mean two sigma bonds
and the ‘=’ sign denotes a double bond. For example, the carbonyl carbon, which has two
sigma bonds and a double bond is represented by ‘>C=’. θ is the bond angle. The term ‘σ’
denotes the site diameter and the term ‘k’ in well energies is the Boltzmann constant.

σ

σ Intra

(Å)
Site
(Å)
CH3 (-CH2-/-CH3)
3.68
2.63
CH3 (>CH-/>C<)
3.68
2.63
CH3 (>CH- in isobutane only)
3.68
2.63
CH3 (-O-)
3.68
2.63
CH3 (>C=)
3.68
2.63
CH3 (-OH in methanol only)
3.68
2.63
CH3 (>C= in acetone only)
3.68
2.63
CH3 (2-butanol only)
3.68
2.63
CH2 (-CH3/-CH2-)
3.68
2.63
CH2 (-CH-)
3.68
2.63
CH2 (-O-)
3.68
2.63
CH2 (>C=)
3.68
2.63
CH2 (-OH in ethanol only)
3.68
2.63
CH2 (-OH)
3.68
2.63
CH (in a branched alkane)
3.67
2.67
CH (-OH in secondary alcohols)
3.67
2.69
C (-CH3/-CH2-/>CH-/>C<)
3.66
2.63
CH2= (-CH=)
3.52
2.52
CH= (in α-olefins and α,ω-diolefins)
3.4
2.4
CH= (in cis/trans olefins only)
3.4
2.4
>C= (in ketone)
3.31
2.92
>C= (in ester)
3.31
2.92
-OH (-CH2- in primary alcohols)
2.8
1.8
-OH (-CH- in secondary alcohols)
2.8
1.8
-O- (-CH2-/CH3 in ester or ether)
2.47
1.83
O= (>C= in ketone or ester)
2.77
1.83
Association parameters

λ

κ

OH (primary or secondary alcohols) 0.000496 0.994347

ε 1/k
(K)
93.4
54.7
70.4
115.9
115.78
103.3
100.0
93
35.4
35.0
54.04
69.8
44.7
36
15.3
20.3
4.0
78.8
56.9
57.3
45.0
40.4
140
116
55.4
108.9

Angle
(θ°)
109.5
109.5
109.5
109.5
109.5
109.5
109.5
109.5
109.5
109.5
109.5
109.5
109.5
109.5
109.5
109.5
109.5
125
125
125
120
120
109.5
109.5
109.5
120

ε 4/k
(K)
13.0
27.3
20.4
6
11
34.5
26.5
16
17.2
15.7
8.1
5.8
27.0
25.2
9.1
8.9
0.7
6.9
12.8
8.9
10
2.0
34
6.2
37.3
12.8

ε AD(kcal/mol)
4.69423

The two intermediate well depths were found by interpolating between the inner and the
outer well. Wherever possible, a broad range of molecular weights was used in the training set
for each functional group, to account for the smallest contributions from a group, and hence to

42

improve the transferability of the parameters. Table 4.1 shows the optimized parameter values of
all the sites that were optimized in this study.

4.2 Discussion of results with StePPE parameters and comparison with other UA models
1. n- Alkanes
n-Alkanes, composed of methyl (-CH3) and methylene (-CH2-) groups, are saturated

hydrocarbons used in the chemical and petrochemical industry as starting materials and solvents,
and are also a critical source of energy for industrial, transport, as well as household needs. They
are also of vital importance to academic research of organic compounds, since the functional
groups in alkanes are also present in other homologous series, like unsaturated hydrocarbons,
alcohols, acids, aldehydes, ketones and esters. Because of their importance as such, and
ubiquitous presence in organic compounds, large amount of reliable thermodynamic data for
alkanes is available in literature. Hence, they are always used as the first class of compounds to
test the accuracy and transferability of potential models.
In the UA models, two sites with appropriate bonding parameters are required to define
the alkanes. In our model, the CH3 and CH2 parameters were obtained by fitting the
experimental liquid densities and vapor pressures of ten alkanes. The parameters were then used
to predict the liquid densities and vapor pressures of six other alkanes, shown in Table 4.1. The
overall deviation for liquid densities is less than 1%, and for vapor pressures is 6.1%. The
reduced temperature range shown in Table 4.2 and the following tables is only for accepted
experimental data in DIPPR that was used for the optimizations. The comparisons mostly range
to wider temperature ranges.
Among the continuous force-fields for alkanes, OPLS[5] , which was parameterized for
liquid densities and heats of vaporization of short chain alkanes near ambient conditions, has

43

been observed to show inaccuracies at elevated temperatures and pressures for both, liquid
densities[2, 6] and vapor pressures[2]. TraPPE –UA[2] predicts the liquid densities extremely
well (AAD < 1%) for the alkanes ranging between ethane and n-decane, but higher errors have
been reported[3] for longer normal alkanes. The reported TraPPE vapor pressures on the other
hand, show an average deviation of about 26%[1]. The errors appear to be higher for the shorter
n-alkanes and at lower reduced temperatures. The NERD[7] force-field predicts the liquid
density of alkanes longer than n-dodecane better than TraPPE, but vapor pressure calculations
for the long alkanes are less accurate[8].
sat

Table 4.2: Vapor pressures and liquid densities of n-alkanes with StePPE. P

, Tr
max
Tr
in this and the subsequent tables mean vapor pressure, minimum reduced
temperature and maximum reduced temperature
sat

Compound
nC2
nC3*2
nC4*
nC5*
nC6*
nC8
nC9
nC10*
nC12*
nC13
nC14
nC16*
nC18*
nC19*
nC22
nC24*
Overall

1

P
%AAD %Bias
2.27
0.28
9.97
-9.97
14.94
-14.94
8.47
-8.47
7.07
-7.07
6.49
-6.45
4.23
-4.16
4.57
-3.47
4.45
-2.31
4.2
1.02
4.72
0.77
6.85
3.38
7.27
4.55
7.62
2.11
10.19
1.71
9.26
-6.89
6.14

ρ
%AAD %Bias
3.48
-1.6
1.86
-1.82
0.94
-0.16
1.04
-0.36
0.69
-0.02
0.82
0.33
0.6
0.2
0.81
0.46
0.77
0.5
0.93
0.72
0.21
0
0.19
0.13
0.32
0.12
0.44
0.16
0.15
-0.08
0.09
-0.09

-2.71

0.75

1

Tr

min

0.37
0.24
0.33
0.31
0.38
0.4
0.38
0.39
0.4
0.4
0.4
0.4
0.4
0.4
0.4
0.43

Tr

min

, and

max

0.9
0.97
0.89
0.89
0.9
0.88
0.9
0.9
0.87
0.85
0.81
0.82
0.84
0.84
0.85
0.73

0.07

The ‘n’ in ‘nCx’ implies a normal alkane and ‘x’ denotes the chain length. For example, nC8 is
n-octane
2
The compounds with ‘*’ in this and the following tables were used in the training set
44

Unlu and Elliott[9] first reported vapor pressures and liquid densities of 12 n-alkanes
with step potentials between methane and n-hexadecane in 2004. With the optimized parameters,
the deviations in vapor pressure and liquid density were 4.3% and 4.0% respectively. They report
an improvement in the vapor pressure prediction by an order of magnitude compared to TraPPE.
Also, the SPEAD predictions range to lower reduced temperatures in comparison. It should be
noted however, that the adjustable parameters for n-alkanes in TraPPE are four, against six in
SPEAD. The density deviations in SPEAD on the other hand, showed a systematic
overprediction with increasing chain length (see Fig 3.2)). In a later publication[10], the
parameters were refined to include larger molecular weight compounds, at the expense of ethane.
With the new parameters, we see an improvement in the vapor pressure predictions of the
intermediate alkanes- n-heptane to n-dodecane, but the densities still remain overpredicted (Fig
3.5).

600

800

(a)

(b)

DIPPR
TraPPE
SPEAD
StePPE

700

500

600

T(K)

T(K)

400
500

300
400
200

300
200

100
0.3

0.4

0.5

0.6

0.7

0.8

0.3

0.4

0.5

0.6

0.7

0.8

ρ (g/mL)

ρ (g/mL)

Figure 4.1: Liquid densities from left to right of (a) n-propane, n-pentane and n-octane,
and (b) n-dodecane and n-hexadecane. The symbol meanings are same in (a) and (b)
45

The StePPE liquid densities are significantly improved in comparison to SPEAD, and the
vapor pressures are particularly improved for C16 and larger alkanes. Figures 4.1 and 4.2 show a
comparison of the liquid densities and vapor pressures for selected alkanes with TraPPE,
SPEAD, and StePPE. For molecules larger than propane, the liquid densities are overpredicted
with SPEAD and the error is maximum for n-hexadecane. With StePPE, not only is the
systematic error in the liquid densities eliminated for all but the smallest compounds, the
predictions range to lower temperatures. The densities for n-dodecane with TraPPE and StePPE
are very similar, but for n-hexadecane, StePPE seems to be more accurate than TraPPE at lower
reduced temperatures. The vapor pressures for the entire range of compounds with StePPE seem
to agree well with the experimental data.

10
DIPPR
SPEAD
StePPE
TraPPE

1

P(MPa)

0.1
0.01
0.001
0.0001
0.00001
0.000001
0.001

0.002

0.003

0.004

0.005

0.006

0.007

0.008

-1

1/T(K )
Figure 4.2: Vapor pressures, from right to left, of, n-propane, n-pentane, n-octane, ndodecane and n-hexadecane. For interpretation of the references to color in this and all
other figures, the reader is referred to the electronic version of this dissertation.

46

In StePPE, the sizes of CH3 and CH2 were kept same based on the excluded volume
considerations discussed above. Therefore, we have one less adjustable parameter compared to
SPEAD for n-alkanes. Also, unlike OPLS and NERD, no special parameters were used for the
smaller compounds.

2. Branched alkanes
Branched alkanes are important for the petroleum refining and petrochemical industries,
as they are good antiknock agents for gasoline engines, and have good lubrication properties
compared to their straight chain counterparts. A good design of the processes technologies for
refining and separation of branched alkanes requires reliable phase equilibrium data. The
branched alkane isomers also provide an attractive set of compounds to test a force-field for
property differences based on differences in molecular arrangements.
Of the force-fields developed so far, not many have been optimized for branched alkanes.
SPEAD parameters for branched alkanes have not been published. The TraPPE[11, 12] force –
field reproduces the liquid densities of short branched alkanes well (Fig 4.3 (a)), but the vapor
pressures are overpredicted, especially for the smaller molecules and at lower reduced
temperatures (Fig 4.3 (b)). The %AAD for the vapor pressures of the compounds in their study is
about 23%. Also, the liquid densities for compounds with less than six carbon atoms appear to be
underpredicted, while for the bigger compounds, an opposite trend is observed. Rather than
changing the site diameters, using a different potential depth for the side methyl group could be
helpful in rectifying the errors, but might necessitate the use of special parameters for some of
the smaller compounds, an approach used in the NERD[13-15] force-field. In comparison to
TraPPE, the NERD force-field has been tested for a wider range of branched alkanes and can
distinguish between structural isomers, but special parameters have been used for the methyl

47

groups in isobutane, isopentane, neopentane and neohexane. The vapor pressures of the branched
alkanes with NERD, except for a small set of compounds, are still to be tested. StePPE uses a
similar approach to NERD, but with special parameters only for isobutane.

Table 4.3: Vapor pressures and liquid densities of branched alkanes with StePPE
sat

3

Compound
C3-2-Me*
C3-2.2-diMe*
C4-2-Me*
C5-2-Me
C5-3-Me
C4-2.3-diMe
C4-2.2.3-triMe
C6-2-Me
C6-3-Me
C5-3-Et*
C6-2.5-diMe
C5-2.2.4-triMe
C6-3-Et
C6-2.4-diMe
C5-2.3.4-triMe
C7-2-Me
C6-2.2-diMe*
C7-4-Me*
C6-3.4-diMe*
C6-2.2.5-triMe
C7-2.6-diMe
C7-3-Et*
C9-5-Me*
C9-2-Me
Overall

P
%AAD %Bias
0.85
0.21
2.57
-1.08
5.07
-5.07
3.12
3.1
18.45
-18.45
5.05
4.53
9.77
5.94
2.37
1.56
10.74
-10.74
6.18
-0.64
4.79
4.79
28.36
28.36
7.05
-0.91
6.18
5.3
26.01
26.01
3.3
1.11
6.17
-0.42
4.96
-2.81
9.62
-7.11
9.43
-2.33
9.79
9.79
7.48
-7.37
5.37
-4.49
4.86
4.86
9.55

ρ
%AAD %Bias
0.79
-0.76
1.62
1.62
1.46
-1.25
1.05
-1.05
1.52
1.52
0.49
-0.49
2.61
2.61
0.51
-0.2
0.6
0.38
1.7
1.7
0.52
0.05
1.08
0.38
0.49
0.46
1.24
1.19
1.46
-1.46
0.74
-0.25
2.6
2.6
0.61
0.32
1.57
1.53
2.42
2.42
0.49
-0.49
1.18
1.18
0.65
0.65
0.16
-0.16

3.80

1.10

3

Tr

min

0.28
0.59
0.32
0.28
0.3
0.29
0.47
0.29
0.32
0.32
0.33
0.3
0.48
0.49
0.29
0.29
0.31
0.31
0.48
0.33
0.33
0.31
0.34
0.33

Tr

max

0.89
0.89
0.86
0.88
0.88
0.89
0.89
0.89
0.89
0.89
0.89
0.89
0.89
0.89
0.89
0.89
0.89
0.88
0.89
0.89
0.89
0.88
0.89
0.89

0.35

The first number after ‘C’ denotes the number of carbon atoms in the backbone, the numbers
between the dashes denote the location of side groups, and the suffixes denote the type of side
group. ‘Me’ stands for a methyl group and ‘Et’ stands for an ethyl group. For example, ‘C6-2.2diMe’ means, 2,2-dimethylhexane.

48

600

10

(a)

(b)

DIPPR
StePPE
TraPPE

500

Psat (MPa)

1

T (K)

400

300

0.1

0.01
200

100
0.35

0.45

0.55

0.65

0.75

0.85

0.001
0.001

0.002

0.003

0.004

0.005

0.006

1/T (K-1)

ρ (g/mL)

Figure 4.3: (a) Liquid densities of (left to right) isopentane, neopentane and 2,2dimethylhexane; (b) saturation pressures of (right to left) isopentane, neopentane and 2,2dimethylhexane. The legend meanings are same in (a) and (b)

The TraPPE parameters, though somewhat inadequate for a wide range of compounds,
are able to distinguish isomers with structural differences, without any special parameters (Fig
4.4). The only difference in the compounds in Figure 4.4 is in the intramolecular interactions due
to the different arrangement of the CH and CH3 groups. StePPE also makes this distinction
reasonably well, but positive deviations are observed in the liquid density predictions of
compounds with a quaternary carbon (Table 4.3). Treating the methyl side chains attached to a
quaternary carbon explicitly as in NERD, might alleviate the issue somewhat, but we might have
to compromise the accuracy in the vapor pressures of some of these compounds which seem to
be slightly underpredicted. Despite the discrepancies, StePPE reproduces the liquid densities and
vapor pressures of a wide range of branched alkanes very well.

49

550

10

(a)

(b)
3,4-dimethylhexane
2,5-dimethylhexane

Psat (MPa)

T(K)

450

350

1

0.1
250

150
0.3

0.4

0.5

0.6

0.7

0.8

0.01
0.0015

0.002

0.0025

0.003

0.0035

1/T (K-1)

ρ (g/mL)

Figure 4.4: Liquid densities (a) and saturation pressures (b) of 3,4-dimethylhexane and 2,5dimethylhexane. The lines are from correlations in DIPPR, filled symbols are from TraPPE
and the open symbols are from StePPE

Optimizing the branched alkanes for a wide range of compounds may not be that trivial
after all, because of the multiple factors that need to be taken into account to incorporate the
effect of branching on the densities as well as the shapes of the VLCC. For example, with
increased branching, naively, we would expect the densities to become smaller because of poor
packing and reduced dispersive attractions due to screening. But with much branching, we see an
opposite trend. A look at the experimental C9 isomer densities confirms this behavior (Fig 4.5).

50

600

2,2,4,4- tetramethylpentane
2,2,5- trimethylhexane
n-nonane
2,2,3,3-tetramethylbutane
2,6-dimethylheptane

550

T(K)

500

450

400

350
0.4

0.45

0.5

0.55

0.6

0.65

0.7

ρ (g/mL)

Figure 4.5: Experimental liquid densities of C9 isomers and 2,2,3,3- tetramethylbutane.
The curves are from equations in DIPPR 801 database

3. n- Alkenes
Because of the unsaturated double bonds present in them, alkenes are at the core of a
variety of useful chemical reactions as starting materials and intermediates. They are therefore,
aptly called the building blocks of the petrochemical and polymer industries. Reliable
thermodynamic data is crucial for industrial processes involving alkenes.
To represent the alkenes, two new sites (-CH= and CH2=) were introduced in the StePPE
2

database. The sp CH parameters were kept the same for the α-olefins, dienes but for cis/trans
2

alkenes, the well depths of the cis/trans sp CH were slightly different. Of the other available
force-fields for olefins, OPLS[5], Spyriouni et al’s model (SET)[16] and TraPPE[17] use only
two unique sites for the olefins. The OPLS parameters were optimized to reproduce

51

thermodynamic and conformational properties of 1 and 2-butenes at near ambient conditions,
and are not transferable to other alkenes[17]. The SET model seems to underpredict the liquid
densities of alkenes smaller than n-butene and overpredict the liquid densities of larger
alkenes[17, 18] and is increasingly less accurate than TraPPE and NERD for longer chain
lengths. The vapor pressure predictions with SET seem to be more accurate than TraPPE and
NERD. Like the SET model, NERD too was parameterized for α-olefins only, and represents the
vapor liquid coexistence densities of alkenes from propene to 1-octene well. However, NERD
2

uses special parameters for the methyl group in propene and the sp CH2 in ethylene. A
comparison of liquid density and vapor pressure predictions for α-olefins with TraPPE, SPEAD,
NERD and StePPE is shown in Figure 4.6.
The SPEAD model reproduces the liquid densities of α and higher olefins and dienes of
chain lengths upto C8 with an overall error of about 1.6%, using four new sites for alkenes- sp

2

2

CH2 for propylene and other alkenes, and sp CH for α-olefins, and dienes[1]. However, the
error in the α-olefins seems to increase with chain length from within 1% in 1-butene to 4.4% in
1-octene. The predictions with TraPPE and StePPE are in close agreement with experimental
data for the alkenes shown. The liquid densities with NERD however seem to be somewhat
undepredicted for the smaller alkenes at low reduced temperatures. The vapor pressures on the
other hand, are consistently overpredicted with TraPPE, while both the step potential models
agree well with experiments.

52

600

(a)

500

10
SPEAD
StePPE
NERD
TraPPE
DIPPR

(b)

1

T(K)

P(MPa)

400

300

0.1

0.01

200

100
0.4

0.5

0.6

0.7

0.8

ρ (g/mL)

0.001
0.001

0.003

0.005
-1

1/T(K )

Figure 4.6: (a) liquid densities (left to right) of propene, 1-butene, 1-pentene and 1-octene;
(b) vapor pressures of (right to left), propene, 1-butene, 1-pentene and 1-octene. The
symbol meanings are same in (a) and (b)

The liquid densities of the α, ω- diolefins seem to agree well with experiments, except 1,
4- pentadiene (Table 4.4). However, there have been very few studies, and in a very limited
temperature range, for densities of 1, 4- pentadiene, which show considerable scatter. Vapor
pressure measurements are also very rare in literature for 1, 4- pentadiene. In case for 1, 9decadiene, density data is limited and scattered too, but its vapor pressure data seems to be more
reliable.

53

10

350

1

300

0.1

T(K)

P(MPa)

400

250

cis-2-butene
cis-2-butene-StePPE
trans-2-butene
trans-2-butene-StePPE
1-butene
1-butene-StePPE

0.01

0.001

200

(b)

(a)
150
0.4

0.5

0.6

0.7

0.8

0.0001
0.002
0.9

0.003

0.004

0.005

0.006

0.007

-1

ρ (g/mL)

1/T(K )

Figure 4.7: liquid densities (a) and vapor pressures (b) of butene isomers. The meanings of
the legends are same in (a) and (b)

With the current set of parameters, the errors in cis/trans isomers are low but, it appears
that a very clear distinction between the isomers is not made for the liquid densities (Fig 4.7 (a)).
The liquid density predictions for cis-2-butene and trans-2-butene almost overlap each other. But
the trends in vapor pressures are distinguishable for the isomers (Fig 4.7 (b)). TraPPE parameters
for the alkenes are able to make that distinction without special cis/trans parameters, but have not
been tested for isomers of alkenes larger than butene. It is not clear whether SPEAD uses
2

separate parameters for the sp CH2 and CH groups for cis/trans compounds. Perhaps, different
sets of parameters for the cis/trans groups may be required to correctly differentiate the trends
between those isomers in StePPE, which would be a subject of future study.

54

Table 4.4: Vapor pressures and liquid densities of alkenes with StePPE
sat

P
4
Compound
%AAD %Bias
C3-1-ene*
10.93
-10.93
C4-1-ene
5.09
-3.40
C4-1.3-diene* 27.14
27.14
C4-2c-ene*
12.94
-12.94
C4-2t-ene*
5.13
2.04
C5-1-ene*
3.08
-1.66
C5-2t-ene
14.64
14.64
C6-1.5-diene* 12.40
10.73
C6-1-ene
3.85
-3.11
C7-1-ene*
3.53
-3.24
C7-2c-ene
6.21
-4.62
C8-1-ene
3.95
-2.41
C8-2t-ene
13.93
13.51
C5-1.4-diene
20.01
20.01
C10-1.9-diene
8.98
-7.80
Overall

10.35

ρ
%AAD %Bias
0.47
-0.10
0.50
0.41
0.71
0.36
0.59
0.53
1.50
-1.48
0.83
0.57
1.35
-1.35
0.40
0.07
1.05
1.04
0.26
-0.04
0.12
0.04
0.51
0.43
1.03
-1.03
2.70
2.70
0.54
0.27

2.38

0.81

Tr

min

0.28
0.30
0.39
0.45
0.43
0.33
0.47
0.52
0.30
0.48
0.48
0.31
0.45
0.65
0.33

Tr

max

0.89
0.89
0.89
0.88
0.89
0.90
0.72
0.90
0.90
0.73
0.73
0.87
0.88
0.88
0.58

0.02

4. Primary alcohols
The current implementation of StePPE does not include electrostatic interactions due to
partial charges, which are expected to be present in compounds with highly electronegative sites
like oxygen. In alcohols, partial charges would be present on the hydroxyl group and the αcarbon, which is implemented in many of the continuous potential models which include
electrostatic interactions in their force-field, such as OPLS[19], TraPPE[20], and NERD[21]. In
StePPE, a hydrogen bonding square well similar to SAFT[22] is used to account for association,
and its effect is added to the EoS as discussed in the ‘SPEAD theory’ section in CHAPTER III .

4

The prefix ‘Cx’ denotes the chain length, the number(s) in between the dashed denote the
position(s) of the double bond and the letters ‘c’ or‘t’ denote a cis or a trans case, and the suffix
‘ene’ denotes one double bond and ‘diene’ denotes two double bonds. For example, C8-2t-ene
denotes trans-2-octene
55

Also, the alkyl groups bonded to the hydroxyl group are optimized to give well parameters
different than those obtained from fitting the alkanes to compensate partially for the effect of the
charges. Unlike the UA models listed above, where the hydrogen atom in the hydroxyl group is
treated explicitly, OH in SPEAD and StePPE is one group. Special parameters are used for
methyl group in methanol in all the models. SPEAD also uses additional special parameters for
ethanol and propanol compared to the rest of the primary alcohols studied. Though SPEAD
reports a higher number of compounds simulated, the liquid densities are not reproduced well for
the longer alcohols (Fig 4.8 (a), 4.9). The vapor pressures deviations however, are generally
small with SPEAD (Fig 4.8 (b)). The OPLS parameters developed for the short chain alcohols
show poor transferability for vapor pressures and liquid densities for longer alcohols or elevated
temperatures [23]. The liquid densities with TraPPE and NERD are in good agreement with the
experimental values for methanol to 1-octanol, the longest n-alcohol studied in TraPPE and
NERD. The vapor pressure predictions with TraPPE seem to show a positive deviation (Fig 4.8
(b)). With NERD, the vapor pressure predictions need verification. With the availability of more
reliable data for alcohols larger than 1-octanol, the transferability of the parameters developed
for TraPPE and NERD force-fields can be tested in future.
The errors for the liquid densities and vapor pressures with StePPE are summarized in
Table 4.5. For the eleven compounds in the training set and six in the validation set, the overall
deviation in vapor pressure is about 8%, which is slightly better than SPEAD, where this
deviation was 9% for a similar range of reduced temperatures. The deviations in the liquid
densities however, are a more significant improvement, where the corresponding value is 0.72%,
compared to 2.49% in SPEAD. Also, the StePPE parameters are tested to longer chains. The
biggest alcohol tested in SPEAD is 1-hexanol, while in StePPE, it is 1-eicosanol. The liquid

56

density and vapor pressure data is limited for the longest alcohols, but within the temperature
range where data is available, the agreement is excellent with StePPE. The only anomalies are
the deviations in liquid density for methanol and vapor pressure of 1-octadecanol (Table 4.4).
Except methanol, the individual density errors are less than 1% with StePPE, whereas in
SPEAD, half of the compounds studied show a density deviation of above 2.5%. Moreover, the
density deviations with SPEAD tend to increase with increasing chain length (Fig 4.8 (a) and Fig
4.9), indicating that the –CH2– sizes might have been too small. The discussed improvement in
StePPE is achieved with three less site types (corresponding to eight adjustable parameters- six
well parameters and two hydrogen bonding parameters) compared to SPEAD, which uses special
parameters for the alkyl group in propanol in addition to the hydroxyl group in ethanol and
propanol.

Table 4.5: Vapor pressures and liquid densities of primary alcohols with StePPE
sat

Compound
Methanol*
Ethanol*
Propanol*
1-butanol
1-pentanol*
1-octanol*
1-nonanol
1-decanol
1-dodecanol*
1-tridecanol
1-pentadecanol
1-hexadecanol
1-heptadecanol
1-octadecanol
1-nonadecanol
1-eicosanol*

P
%AAD %Bias
4.25
-1.04
1.39
-0.77
4.76
3.72
9.73
-9.73
12.93
12.93
3.33
3.10
4.15
2.01
4.45
-2.29
8.47
-6.97
9.96
2.80
7.38
4.72
5.98
-1.92
12.32
3.58
26.66
-26.66
13.92
-10.21
11.77
-7.36

%AAD
2.48
0.60
0.86
0.48
0.65
0.26
0.19
0.19
0.20
0.24
0.27
0.19
0.43
0.49
0.44
0.21

%Bias
2.21
0.27
-0.44
0.28
-0.52
-0.15
0.10
-0.01
-0.06
-0.14
-0.27
0.09
-0.43
0.34
-0.44
-0.21

Overall

7.95

0.72

0.05

-1.61

ρ

57

Tr

min

0.40
0.43
0.34
0.34
0.43
0.40
0.41
0.41
0.41
0.41
0.44
0.42
0.43
0.42
0.44
0.43

Tr

max

0.88
0.89
0.90
0.89
0.89
0.84
0.90
0.80
0.86
0.90
0.89
0.74
0.90
0.73
0.89
0.87

10

(a)

600

(b)

1
0.1

P(MPa)

T(K)

500

SPEAD
StePPE
TraPPE
DIPPR

400

0.001

propanol - SPEAD
Propanol - StePPE
Propanol - TraPPE
Pentanol - SPEAD
Pentanol - StePPE
Pentanol - TraPPE
Octanol - SPEAD
Octanol - StePPE
Octanol -TraPPE

300

0.0001
0.00001

200
0.4

0.5

0.6

0.01

0.7

0.8

0.000001
0.9
0.0015

0.0025

0.0035

0.0045

-1

ρ (g/mL)

1/T(K )

Figure 4.8: Liquid densities (a) and vapor pressures (b) of (right to left) propanol, 1pentanol and 1-octanol. In (a) the dotted, solid and dashed lines are for propanol, 1pentanol and 1-octanol respectively

Methanol is the only compound with a better liquid density prediction with SPEAD at
lower reduced temperatures (Fig 4.10 (a)). It should be noted however, that only the CH3 wells
are special parameters in StePPE for methanol, compared to the well parameters of both, the
CH3 and OH groups, as well as the hydrogen bonding parameters for the OH in SPEAD. We
believe, considering only the temperature dependent hydrogen bonding parameters separately for
methanol might be enough to improve the density predictions.

58

650
600

ethanol
butanol
heptanol
nonanol

550

T(K)

500
450
400
350
300
250
200
0.5

0.6

0.7

0.8

0.9

ρ (g/mL)
Figure 4.9: Liquid densities of selected 1-alcohols with SPEAD and StePPE. Dashed lines
are with SPEAD and the solid lines are with StePPE

600

600
DIPPR
StePPE
SPEAD

(a)
500

(b)
500

400

T(K)

T(K)

400

300

300

200

200

100

100

0

0

0.4

0.6

0.8

1

ρ (g/mL)

0.4

0.6

0.8

1.0

ρ (g/mL)

Figure 4.10: Liquid densities of methanol (a) and ethanol (b) with SPEAD and StePPE. The
symbol meanings are same in (a) and (b)
59

5. Secondary alcohols
New sites used in defining the branched alcohols were the side hydroxyl group and the
CH group to which it is bonded. The size and the association parameters for the hydroxyl group
were borrowed from the optimizations of primary alcohols. The optimized wells of the hydroxyl
group were shallower than the hydroxyl in the n-alcohols, possibly because of some steric
hindrance, as it is in the middle of the chain. A similar observation was made for the side methyl
groups in branched alkanes, which are less exposed than the terminal methyl groups in long
alkanes. The CH3 potential wells were optimized separately for 2-propanol.
Very few force-fields have tested their optimized parameters for branched alcohols.
Figure 4.11 shows the liquid densities and vapor pressures of 2-propanol and 2-butanol with
TraPPE[20] and StePPE. The liquid densities of both alcohols are overpredicted by about 1%
with StePPE, while with TraPPE, the predictions appear to be underpredicted at higher reduced
temperatures. The vapor pressures appear to be somewhat overpredicted with both models.
500

10

(a)

2-propanol
2-butanol
1

400

0.1

T(K)

P(MPa)

450

(b)

350

300

0.01

0.001

250
0.5

0.6

0.7

0.8

0.9

0.0001
0.0015

0.0025

0.0035

0.0045

-1

ρ (g/mL)

1/T(K )

Figure 4.11: Liquid densities and vapor pressures of 2-butanol (circles) and 2-propanol
(squares) with TraPPE and SPEAD. The closed symbols are with TraPPE and the open
symbols are from StePPE.

60

Table 4.6: Vapor pressure and liquid density deviations of secondary alcohols with StePPE
sat

P
5
Compound %AAD %Bias
C3-2ol*
18.04
9.2
C4-2ol
10.6
10.6
C5-2ol*
10.36
10.36
C6-2ol
5.71
0.83
C7-2ol*
4.52
2.6
C8-2ol
5.63
0.4
C9-2ol*
5.49
-4.24
C5-3ol
7.99
-7.43
C6-3ol*
17.36
-17.29
C7-3ol
18.03
-16.39
Overall

.

9.39

2.65

ρ
%AAD %Bias
1.16
1.16
0.72
0.72
0.46
0.42
1.34
1.34
0.67
0.67
0.98
0.98
0.8
0.8
0.23
-0.23
0.16
0.1
0.56
0.54
0.58

Tr

min

0.54
0.52
0.42
0.41
0.4
0.4
0.41
0.44
0.45
0.4

Tr

max

0.9
0.9
0.7
0.89
0.71
0.76
0.62
0.73
0.71
0.71

0.52

Table 4.6 shows the liquid density and vapor pressure deviations with StePPE for

secondary alcohols with chain lengths from C3 to C9. For the 2-alcohols, the liquid density
deviations seem to be positive and increasing with chain length, which can be a result of the CH2
diameter being small. But there isn’t a clear trend in the density deviations of the 3-alcohols. On
the other hand, the vapor pressures of the 3-alcohols appear to be underpredicted, and the error
increases with chain size, suggesting that the potential well for CH2 may be too deep. However,
for the 2-alcohols, the vapor pressure deviations seem to be positive and appear to decrease with
chain length (if we ignore 2-hexanol, for which very few reliable data points are available in
DIPPR), which suggests the -OH or CH potential might be somewhat shallow. Therefore, there
is no clear trend about the deviations, but with the availability of more reliable data for longer
branched alcohols with the –OH at different locations, a more detailed analysis can be made. A
study with an increased diameter of the CH group is shown in the future work section in
CHAPTER VI.
5

In ‘Cx-yol’, x denotes the backbone chain length, and y denotes the position of the –OH group
in the chain. For example, C5-3ol would mean 3-pentanol.
61

6. Ethers
Ethers are important industrial solvents and also have medicinal uses as sedatives. Ethers
are also used as octane boosters (methyl-tertiary-butyl ether) for gasoline engines and cetane
boosters (diethyl ether) for diesel engines. Thermodynamic data of only the smallest aliphatic
ethers is available. In this study, the vapor pressures and liquid densities of ethers ranging from
dimethyl ether to dinonyl ether are presented.

Table 4.7: Vapor pressure and liquid density deviations of ethers with StePPE
sat

6

Compound
C1OC1*
C1OC3
C2OC2*
C1OC4
C2OC3
C1OC5*
C3OC3*
C4OC2
C2OC6*
C4OC4
C5OC5
C6OC6
C8OC8*
C9OC9
Overall

P

ρ

%AAD %Bias
8.7
8.7
2.89
1.92
8.33
8.33
8.96
-8.23
6.06
6.04
15.93
-13.7
9.18
5.45
8.71
-5.8
5.09
4.95
2.47
0.65
12.67
-7.54
9.07
9.07
17.85
17.85
24.47
20.81
8.56

4.19

%AAD %Bias
1.4
-1.4
1.26
-1.26
0.59
-0.49
0.13
-0.08
0.61
0.32
0.23
-0.22
1.02
-1.02
0.21
-0.17
0.19
-0.13
0.21
0.14
0.21
0.04
0.34
0.26
0.28
0.23
0.32
0.32
0.49

Tr

min

0.45
0.4
0.43
0.31
0.49
0.32
0.38
0.38
0.39
0.47
0.36
0.45
0.41
0.39

Tr

max

0.88
0.67
0.89
0.86
0.89
0.86
0.9
0.9
0.85
0.71
0.89
0.75
0.63
0.9

-0.27

3

The new sites to define ethers were the sp oxygen and the alkyl group adjacent to it. The
sizes of these alkyl groups were kept the same as in alkanes, and only the well depths were
varied, for similar reasons as in the alcohols. The parameters for the rest of the alkyl groups were
used from the optimization of the alkanes. The overall deviation in the liquid density prediction
6

In ‘CxOCy’, x and y denote the chain length of the alkyl group on each side of the oxygen. For
example C2OC3 is ethyl propyl ether
62

is less than 0.5% (Table 4.7) for the fourteen compounds studied. The largest ether reported in
both, TraPPE[24] and SPEAD[1] was dipropyl ether, for which the density deviations appear to
be better with StePPE for a broader temperature range (Fig 4.12). For the other compounds
common in the three models, the predictions are better with TraPPE and StePPE. The biggest
compound in StePPE was dinonyl ether, for which the density is overpredicted at temperatures
close to critical. Tuning the parameters further, for a better fit for dinonyl ether, cannot be
accomplished without increasing the errors in the smaller compounds, for example, dimethyl
ether, whose density is somewhat underpredicted with the current set of parameters. To our
knowledge however, no other model has reported such accuracy for such a wide range of ethers.

750

(a)

650

100
Trappe
SPEAD
StePPE
DIPPR

(b)

1

P (MPa)

T (K)

550
450

0.01

0.0001

350
0.000001

250
150
0.45

0.55

0.65

0.75

0.00000001
0.85
0.001

ρ (g/mL)

0.002

0.003

0.004

0.005

0.006

-1

1/T (K )

Figure 4.12: (a) liquid densities of (left to right) dimethyl ether, dipropyl ether and dipentyl ether and di-nonyl ether; vapor pressures of (right to left) dimethyl ether, dipropyl
ether and di-pentyl ether and di-nonyl ether

63

The overall deviation in the vapor pressure for the compounds studied with StePPE is
8.6%, which is significantly better than TraPPE for the relatively smaller compounds reported in
reference 24. It appears that the TraPPE parameters for the ethers get better with increasing chain
length. For example, the average deviation in the vapor pressure of dipropyl ether was about
40%, but for dipentyl ether, it was within 10%. Both, SPEAD and StePPE agree closely with the
experimental vapor pressures for the common ethers. Compared to the liquid densities, the
agreement is relatively better for vapor pressures for dinonyl ether with StePPE.

6. Ketones
To represent the ketones, two new sites- carbonyl oxygen (O=) and carbonyl carbon
(>C=), were added to the StePPE database. The parameters for the alkyl groups were obtained
from alkane optimizations, except for the alkyl groups attached to the carbonyl carbon, whose
well depths were readjusted. The ketone oxygen in SPEAD[25] has an additional adjustable
parameter, because it is treated as a hydrogen bonding site. Additionally, the wells and the
diameter of the CH2 next to the carbonyl carbon have been reoptimized in SPEAD in order to fit
the liquid densities better. As discussed earlier in this work, readjusting the CH2 size would also
help eliminate the systematic trend in the deviations of alkane densities. The overall average
error in liquid densities with StePPE is below 0.5%, and all of the individual average deviations
are below 1% (Table 4.8) from acetone to 5-nonanone, which is somewhat better than
TraPPE[24] for the common compounds in both studies, except at temperatures very close to
critical. The StePPE predictions however, span to much lower reduced temperatures compared to
TraPPE. The maximum errors observed in the reduced temperature range of about 0.3 to 0.9 for
the sixteen ketones are within 1.5%.

64

Table 4.8: Vapor pressure and liquid density deviations of ketones with StePPE
sat

7

Compound
C3.2one*
C4.2one*
C5.2one
C5.3one*
C6.2one*
C6.3one
C7.2one
C7.3one*
C7.4one
C8.2one*
C8.4one
C8.3one
C9.2one
C9.4one
C9.5one*
Overall

P

ρ

%AAD %Bias
2.5
-2.1
9.86
4.7
6.71
1.29
7.52
5.81
5.93
-5.39
4.88
1.8
9.42
-9.36
3.42
0.14
7.7
0.56
9.88
-9.8
18.48
-11.18
20.47
-19.58
10.3
-10.24
15.84
-2.06
5.9
0.01
8.12

-2

%AAD %Bias
1.42
-1.86
0.53
0.53
0.19
-0.11
0.28
0.28
0.53
0.53
0.28
-0.23
0.38
0.38
0.41
0.41
0.31
-0.29
0.66
0.66
0.3
-0.14
0.26
0.11
0.65
0.65
0.25
-0.25
0.3
0.22
0.37

Tr

min

0.36
0.36
0.36
0.49
0.48
0.5
0.41
0.47
0.4
0.4
0.47
0.45
0.42
0.43
0.44

Tr

max

0.90
0.89
0.89
0.89
0.73
0.72
0.74
0.69
0.88
0.76
0.71
0.71
0.77
0.71
0.76

0.18

The overall error in vapor pressures with StePPE and SPEAD are better than TraPPE.
However, it is important to note that StePPE has three less adjustable parameters in comparison
to SPEAD. In both step potential models, the CH3 in acetone was optimized separately, so it has
different well parameters compared to the CH3 adjacent to >C= in other ketones. If we do not
consider these special parameters for acetone, the error in vapor pressure is about 42%, athough
the density error is still less than 1%. The maximum deviation in vapor pressures with StePPE is
observed with 3-octanone, 4-octanone and 4-noanone, all of which have very scattered
experimental vapor pressure data in DIPPR. It should be noted that we used the experimental
points, rather than the fitted correlations, to tune the parameters in StePPE.

7

In ‘Cx.yone’ x denotes the chain length and y denotes the position of the ketone group. For
example, C6.3one is 3-hexanone
65

600

(a)
500

10
Acetone
2-pentanone
2-octanone

(b)

SPEAD
StePPE
TraPPE
DIPPR

1
0.1

P (MPa)

T(K)

400

300

0.01

0.001
200

0.0001

100
0.5

0.6

0.7

0.8

0.9

0.00001
0.001

0.002

0.003

0.004

0.005

0.006

-1

ρ (g/mL)

1/T(K )

Figure 4.13: Liquid densities (a) and vapor pressures (b) of acetone, 2-pentanone and 2octanone with SPEAD, TraPPE and StePPE. In (a), the filled symbols are from TraPPE
and open symbols are from StePPE. In (b), from top to bottom: acetone, 2-pentanone, 2octanone

8. Esters
Esters are important compounds in chemical industry with a plethora of applications, like,
fragrances, food additives and preservatives, adhesives and polymers. Esters are also the
important constituents in biodiesel[26] fuel blends, owing to their relatively desirable
combustion properties[27]. They are formed by replacing the hydroxyl group (–OH) in an acid
with an alkoxy group (–OR), which comes from an alcohol, or by a transesterification
method[28], where a triglyceride is converted to monoglycerides.
2

The only parameters optimized for esters were the outer well depths of the carbonyl (sp )
carbon (>C=). Its diameter, along with the parameters for the adjacent alkyl groups and the
2

carbonyl (sp ) oxygen (O=), was borrowed from the ketone optimizations. The parameters of the

66

3

alkoxy (sp ) oxygen (–O–) and the alkyl group to which it is attached (–R– (–O–)) are obtained
2

from the optimization of the ethers. The sp carbon in the ketones is bonded to alkyl groups,
3

whereas in esters, one of the bonds is with the sp oxygen. One would expect, as in the alcohols
and ethers, that the well depth for >C= would be higher in the ester in comparison to ketones.
However, the optimized well depths show an opposite trend, based on which, it might be
expected that the well depths for the carbonyl carbon in the ketones might have been too deep.
To test this hypothesis, we tried to reoptimize the ketones with smaller well size for the carbonyl
carbon, but we were unable to get a combination of wells for the carbon and oxygen, which gave
errors similar to or better than the current parameters for those groups.
The overall deviation in the liquid density predictions for a set of 26 esters was found to
be 0.8% and the corresponding deviations in the vapor pressure predictions was 10%. The set
includes multiple acetates, propionates and butanoates, along with methyl esters that are relevant
to biofuels, namely methyl decanoate to methyl stearate. We believe no other model has
considered such a wide range of esters, though the vapor pressure error for some acetates and
propionates appears to be high (between 15% and 21%), which is still better than SPEAD and
the group contribution methods[29]. Three of the nine acetates reported in SPEAD[25] show a
deviation of more than 25% in the vapor pressure prediction. Decylacetate, which is the biggest
acetate in the study, shows a deviation of about 27%. The liquid density predictions appear to be
consistently better with StePPE (Table 4.7). It should also be noted however, that in SPEAD, the
parameters for >C=, O=, –O–, and the methyl group attached to >C= are optimized for the
esters. Therefore, the number of sites for which parameters are adjusted in SPEAD is four,
compared to one in StePPE. The parameters for the group present in esters have also been

67

reported for OPLS[30], TraPPE[31] and NERD[32] force-fields. The OPLS parameters were
developed for methyl acetate and are inaccurate for liquid densities of ethyl acetae[33]. The
NERD parameters have been used to study triglycerides and very limited data is published for
monoesters. Only TraPPE reports results for a small set, but the state points are not published, so
we’re unable to make a comparison.

Table 4.9: Deviations in liquid densities and vapor pressures of esters with StePPE
sat

P
8
Compound %AAD
C2ateC1*
9.38
C2ateC2*
14.35
C3ateC1*
7.92
C2ateC3
7.46
C4ateC1
7
C3ateC2*
13.09
C2ateC4
2.71
C3ateC3
9.35
C4ateC2
11.67
C3ateC4
5.82
C6ateC1
16.33
C2ateC5*
11.56
C4ateC3*
8.52
C2ateC6
21.79
C4ateC4
18.36
C2ateC7
5.21
C5ateC4
16.31
C8ateC1*
15.08
C8ateC2
7.43
C2ateC8*
14.84
C10ateC1
17.5
C2ateC10
18.41
C12ateC1*
17.24
C14ateC1*
5.14
C16ateC1
5.19
C18ateC1*
4.14
Overall

10.03

ρ

%Bias
9.32
14.35
3.94
7.23
-3.75
13.09
-1.03
5.39
-1.72
3.64
-16.33
-10.44
6.33
-21.79
-5.38
-5.07
-9.85
-15.08
-2.88
-14.84
-17.13
-18.41
-17.24
-4.12
-2.93
-0.49
-2.78

%AAD %Bias
1.31
-1.28
1.62
-1.62
1.22
-1.22
1.04
-0.95
1.05
-1.05
1.22
-1.22
0.41
-0.41
1.28
-1.28
1.1
-1.1
0.49
-0.49
0.2
0.16
0.17
-0.08
0.9
-0.9
0.19
0.15
0.42
-0.31
0.35
0.35
0.3
0.11
0.54
0.54
0.07
-0.01
0.27
0.27
0.68
0.68
0.54
0.54
0.84
0.84
0.82
0.82
0.72
0.72
0.64
0.64
0.8

8

Tr

min

0.39
0.52
0.48
0.48
0.47
0.48
0.47
0.46
0.45
0.46
0.78
0.46
0.46
0.44
0.41
0.43
0.42
0.70
0.45
0.38
0.57
0.41
0.57
0.61
0.61

Tr

max

0.89
0.89
0.89
0.90
0.89
0.88
0.71
0.74
0.74
0.73
1.00
0.81
0.75
0.62
0.89
0.73
0.81
1.00
0.74
0.86
0.87
0.68
1.00
1.00
1.00

-0.55

‘CxateCy’ is ROOR’. Cx (R) is from acid, Cy (R’) is from alcohol part in the saturated ester.
For example, C3ateC2 means propanoic acid ethyl ester
68

REFERENCES

69

4.3 References
1.

Unlu, O.; Gray, N.H.; Gerek, Z.N.; Elliott, J.R., "Transferable Step Potentials for the
Straight-Chain Alkanes, Alkenes, Alkynes, Ethers, and Alcohols". Industrial &
Engineering Chemistry Research 2004, 43, 1788-1793.

2.

Martin, M.G.; Siepmann, J.I., "Transferable Potentials for Phase Equilibria. 1. UnitedAtom Description of n-Alkanes". Journal of Physical Chemistry B, 1998, 102, 25692577.

3.

Nath, S.K.; Escobedo , F.A.; de Pablo, J.J., "On the Simulation of Vapor–Liquid
Equilibria for Alkanes". Journal of Chemical Physics, 1998, 108, 9905-9911.

4.

SHERPA – An Efficient and Robust Optimization/Search Algorithm, Red Cedar
Technologies: www.redcedartech.com.

5.

Jorgensen, W.L.; Madura, J.D.; Swenson, C.J., "Optimized Intermolecular Potential
Functions for Liquid Hydrocarbons". Journal of the American Chemical Society, 1984,
106, 6638-6646.

6.

Benet, J.; MacDowell, L.G.; Menduin, C., "Liquid-Vapor Phase Equilibria and Surface
Tension of Ethane as Predicted By the Trappe and OPLS Models†". Journal of Chemical
& Engineering Data, 2010, 55, 5465–5470.

7.

Nath, S.K.; Escobedo, F.A.; de Pablo, J.J.; Patramai, I., "Simulation of Vapor-Liquid
Equilibria for Alkane Mixtures". Industrial & Engineering Chemistry Research 1998, 37,
3195-3202.

8.

Errington, J.R.; Panagiotopoulos, A.Z., "A New Intermolecular Potential Model for the nAlkane Homologous Series". Journal of Physical Chemistry, 1999, 103, 6314-6322.

9.

Cui, J.Y.; Elliott, J.R., "Phase Diagrams for a Multistep Potential Model of n-Alkanes by
Discontinuous Molecular Dynamics and Thermodynamic Perturbation Theory". Journal
of Chemical Physics, 2002, 116, 8625-8631.

10.

Elliott, J.R.; Gray, N.H., "Asymptotic Trends in Thermodynamic Perturbation Theory".
Journal of Chemical Physics, 2005, 123, 184902-1-184902-6.

70

11.

Martin, M.G.; Siepmann, J.I., "Novel Configurational-Bias Monte Carlo Method for
Branched Molecules. Transferable Potentials for Phase Equilibria. 2. United-Atom
Description of Branched Alkanes". Journal of Physical Chemistry B, 1999, 103, 45084517.

12.

Zhuravlev, N.D.; Martin, M.G.; Siepmann, J.I., "Vapor-Liquid Phase Equilibria of
Triacontane Isomers: Deviations from Principle of Corresponding States". Fluid Phase
Equilib., 2002, 202, 307-324.

13.

Nath, S.K.; de Pablo, J.J., "Simulation of Vapor-Liquid Equilibria for Branched
Alkanes". Molecular Physics, 2000, 98, 231-238.

14.

Nath, S.K.; Khare, R., "New Forcefield Parameters for Branched Hydrocarbons". Journal
of Chemical Physics, 2001, 115, 10837-10845.

15.

Wescott, J.T.; Kung, P.; Nath, S.K., "Vapour–Liquid Coexistence Properties and Critical
Points of Two Branched Alkanes Series". Fluid Phase Equilibria, 2003, 208, 123–139.

16.

Spyriouni, T.; Economou, I.G.; Theodorou, D.N., "Molecular Simulation of α-Olefins
Using a New United-Atom Potential Model: Vapor-Liquid Equilibria of Pure Compounds
and Mixtures". Journal of the American Chemical Society, 1999, 121, 3407-3413

17.

Wick, C.D.; Martin, M.G.; Siepmann, J.I., "Transferable Potentials for Phase Equilibria.
4. United-Atom Description of Linear and Branched Alkenes and alkylbenzenes".
Journal of Physical Chemistry B, 2000, 104, 8008-8016.

18.

Nath, S.K.; Banaszak, B.J.; de Pablo, J.J., "A New United Atom Force Field for αOlefins". Journal of Chemical Physics, 2001, 114, 3612-3617.

19.

Jorgensen, W.L., "Optimized Intermolecular Potential Functions for Liquid Alcohols".
Journal of Physical Chemistry, 1986, 90, 1276-1284.

20.

Chen, B.; Potoff, J.J.; Siepmann, J.I., "Monte Carlo Calculations for Alcohols and Their
Mixtures with Alkanes. Transferable Potentials for Phase Equilibria. 5. United-Atom
Description of Primary, Secondary, and Tertiary Alcohols". The Journal of Physical
Chemistry, 2001, 105, 3093-3104.

71

21.

Khare, R.; Sum, A.K.; Nath, S.K.; de Pablo, J.J., "Simulation of Vapor-Liquid Phase
Equilibria of Primary Alcohols and Alcohol-Alkane Mixtures". Journal of Physical
Chemistry, 2004, 108, 10071-10076.

22.

Liu, J.; Bowmann II, T.L.; Elliott, J.R., "Discontinuous Molecular Dynamics Simulation
of Hydrogen-Bonding Systems". Ind. Eng. Chem. Res., 1994, 33, 657-964.

23.

van Leeuwen, M.E., "Prediction of the Vapour- Liquid Coexistence Curve of Alkanols by
Molecular Simulation". Molecular Physics, 1996, 87, 87-101.

24.

Stubbs, J.M.; Potoff , J.J.; Siepmann, J.I., "Transferable Potentials for Phase Equilibria. 6.
United-Atom Description for Ethers, Glycols, Ketones, and aldehydes". Journal of
Physical Chemistry B, 2004, 108, 17596-17605.

25.

Baskaya, F.S.; Gray, N.H.; Gerek, Z.N.; Elliott, J.R., "Transferable step potentials for
amines, amides, acetates, and ketones". Fluid Phase Equilib., 2005, 236, 42-52.

26.

Demirbas, A., "Biodiesel from Vegetable Oils via Transesterification in Supercritical
Methanol". Energy Conversion and Management, 2002, 43, 2349-2356

27.

Dagaut, P.; Gail, S., "Chemical Kinetic Study of the Effect of a Biofuel Additive on JetA1 Combustion". Journal of Physical Chemistry A, 2007, 111, 3992-4000.

28.

Demirbas, A., "Biodiesel Fuels from Vegetable Oils via Catalytic and Non-Catalytic
Supercritical Alcohol Transesterifications And Other Methods: A Survey". Energy
Conversion and Management, 2003, 44, 2093–2109.

29.

Bureau, N.; Jose, J.; Mokbel, I.; de Hemptinne, J.-C., "Vapour Pressure Measurements
and Prediction for Heavy Esters". J. Chem. Thermodyn., 2001, 33, 1485-1498.

30.

Briggs, J.M.; Nguyen, T.B.; Jorgensen, W.L., "Monte Carlo Slmulations of Liquid Acetic
Acid and Methyl Acetate with the OPLS Potential Functions". The Journal of Physical
Chemistry, 1991, 95, 3315-3322.

31.

Kamath, G.; Robinson, J.; Potoff, J.J., "Application of Trappe-UA Force Field for
Determination of Vapor-Liquid Equilibria of Carboxylate Esters". Fluid Phase
Equilibria, 2006, 240, 46-55.

72

32.

Sum, A.K.; Biddy, M.J.; de Pablo, J.J., "Predictive Molecular Model for the
Thermodynamic and Transport Properties of Triacylglycerols". Journal of Physical
Chemistry, 2003, 107, 14443-14451.

33.

Kamath, G.; Cao, F.; Potoff, J.J., "An Improved Force Field for the Prediction of the
Vapor-Liquid Equilibria for Carboxylic Acids". Journal of Physical Chemistry, 2004,
108, 14130-14136.

73

CHAPTER V
5. MEASUREMENTS OF VAPOR PRESSURES AND BUBBLE POINTS WITH
EBULLIOMETRY

5.1 Introduction
Vapor pressures of pure compounds can be measures experimentally by using static or
dynamic methods. In the static method, the liquid is degassed to eliminate dissolved gases, so
that the vapor space in the equilibrium chamber is filled only with the vapors of the liquid of
interest and the pressure of the system is measured to give the vapor pressure. In dynamic
methods, the system pressure is fixed and the liquid is heated to bring it to boil. The
corresponding pressure gives the vapor pressure of the system at the measured boiling
temperature. Though static methods are more accurate, they are tedious to degass and the
auxiliary equipment needed is expensive[1]. Ebulliometry[2, 3] on the other hand, which is the
most widely used dynamic method, is fast and has undergone considerable developments in
recent years to give accuracies[4-6] close to static methods.
One of the major issues in ebulliometry is effective removal of superheat[3]. It can be
removed with the help of Cottrell[7], or Washburn[8] type tubes, which carry slugs of liquid and
vapor to an equilibrium chamber where the temperature is measured. On the way to this
chamber, the superheat in the liquid is used to vaporize extra liquid in the Cottrell pump, so that
the temperature being measured is the true equilibrium temperature. Another way to alleviate the
problem is to “activate” the heated surface with powdered glass[2] and use stirring to facilitate
nucleation points for the vapor and provide good mixing. At low pressures getting rid of the
superheat becomes difficult, because even a small liquid head inside the boiling chamber can
create large enough superheat to cause ‘bumping’, which can flash the entire sample at once and

74

lead to cooling of the equilibrium chamber. Measuring the equilibrium temperature under such
circumstances may require combining various methods to reduce the superheat.
In this study, we constructed ebulliometers based on previous models with two important
design considerations, effective removal of superheat at low pressures and use of small amount
of sample, since many of the compounds in high purity form are very expensive. At several
instances during the optimization of parameters in the StePPE model, we encountered scattered
or limited vapor pressure data, which needed measurements. Another important objective of this
study is to make the ebulliometer developed for pure compounds usable for measurements of
bubble points of binary mixtures. Adding the capability to measure bubble points of binary
mixtures would be helpful in achieving that goal and would also be helpful in validating the
StePPE model for binary mixtures in future.

5.2 Ebulliometer designs for pure compounds
The first ebulliometer studied in this work was based on Hoover’s[9] design, which uses
a tungsten wire to heat the liquid (Fig 5.1). The total height of the apparatus is about 8 inches.
The sample holder has an inside diameter of about 1 cm and requires about 10 mL of sample. It
works well when the wire is heated using a flame ( Fig 5.5) at pressures as low as 30 mbar,
below which, removing the superheat becomes difficult and bumping is observed. However,
owing to safety concerns with the use of flames in a lab that handles large amounts of flammable
liquids, this design was abandoned.

75

to vacuum
thermal well for
temperature sensor
cooling jacket
sample holder
condensate return
tungsten wire for
heating

Figure 5.1: Ebulliometer based on Hoover’s[9] design for vapor pressure measurement

Three other designs based on Cottrell’s (Fig 5.2), Washburn’s (Fig 5.3) and stirred flask
(Fig 5.4) designs were tested, and they all work well in a pressure range of approximately 50 –
1000 mbar (Fig 5.5, 5.6). With the Cottrell tube apparatus, foaming was observed at pressures
below 100 mbar for esters. Boiling chips, which were used in that design in order to facilitate
nucleation of vapors, might have added impurities like moisture to the sample. In the Washburn
and the stirred flask design, the heated surface is activated using powdered glass, in order to
provide nucleation points for vaporization, so that boiling chips can be avoided. In the stirred
tank apparatus, stirring is used in addition to the activation to help create good and rapid mixing,
which is the key to effective superheat removal. Since the mixing is good in this design, the
equilibrium temperature is measured by inserting the temperature probe in the liquid itself, rather
than in a separate chamber, as in the Cottrell or Washburn design.

76

b

a

Equilibration chamber

Cottrell pump

Fig 5.2: Cottrell tube ebulliometer. In (a), the entire assembly is shown, (b) shows the
Cottrell tube insert with a funnel shaped bottom to catch the vapors, which get carried to
the equilibrium chamber with the boiling liquid. The contents get squirted on the
temperature probe (not shown) in this chamber. The height of the outside tube is about 6”
and the diameter is about 1”

Insert with
Washburn pump

Activated surface

Figure 5.3: Ebulliometer with Washburn pump

77

Condenser tube
RTD probe
Temperature
readout
Ground glass
particles

Sample holder

heating mantle

Figure 5.4 Stirred round bottom flask ebulliometer for vapor pressure measurements

250
Hoover
Cottrell
ASTM
NIST Antoine Eqn

200

P(mbar)

150

100

50

0
400

420

440

460

480

T(K)
sat

Figure 5.5: Measurement of P
flask ebulliometers

of monoethyl succinate with Cottrell, Hoover and stirred

78

120.00
NIST
Cottrell
Washburn

100.00

T (°C)

80.00
60.00
40.00
20.00
0.00
0

200

400

600

800

P (mmHg)
sat

Figure 5.6: Measurements of water P

with Cottrell and Washburn ebulliometers

5.3 Ebulliometer for measuring T-x of a binary liquid
In addition to errors arising from superheat, compositional errors are also of concern
when measuring the boiling temperature of a binary mixture. The major source of compositional
error is the condensate hold-up on the walls of the ebulliometer. A relatively smaller error is
incorporated due to evaporation of the liquid to fill the vapor space, which can be neglected if the
relative volatilities of the components in the binary mixture do not differ significantly. Olson[10]
gives a guideline for the use of different ebulliometric methods based on the difference in the
boiling points of the compounds of interest, and suggests the use of static methods in case this
difference is above about 120 K, where the errors due to evaporation and condensed liquid holdup would be high.

79

Based on simple material balances, the condensed liquid hold-up can be calculated using
Olson’s[10] or Rogalski’s[11] method, which differ slightly in the approach, but practically give
the same result. For one mole of a sample feed, with V moles of vapor and L moles of liquid:
V + L =1

(1)

Applying material balance on component i, we get:

V . y i + L.xi = 1.qi

(2)

where, yi, xi and qi are the mole fractions of i in the vapor phase, the liquid phase and the feed,
respectively. Using the equilibrium ratio Ki = yi/xi, equation 2 can be written as:

q i = V .K i xi + L.xi

(3)

or:

xi =

qi
qi / L
=
V .K i + L (V / L).K i + 1

(4)

Defining vapor-to-liquid mole ratio as f = V/L, (4) can be simplified to:

xi =

q i (1 + f )
f .K i + 1

(5)

(∵ L + V = 1 ⇒ 1 / L = 1 + V / L = 1 + f )

or:

f =

q i − xi
y i − qi

(6)

The fraction f depends on the construction of the ebulliometer and the rate of heat supplied. In
other words, for a given ebulliometer and an approximate heating rate, f should be independent
of the system. Generally, the heating rate is monitored using a drop counter in the condensate

80

return line of the ebulliometer. Once f is known, equation (5) can be used to calculate xi
iteratively for a given feed composition (qi) for an unknown system, without having to measure
the equilibrium liquid or vapor compositions, given by:

xin +1 =

qi (1 + f )
f .( y in / xin ) + 1

(7)

where, n is the iteration number. The first guess of xi can be the feed composition (qi) and the
corresponding vapor composition (yi) can be found by using any activity model, using the
following equation:

γ i x i Pi sat = y i P
In equation (8),

(8)

γ i is the activity coefficient of component i, Pi sat

is the vapor pressure of

pure component i at the measured bubble temperature and P is the bubble pressure.
This approach has been used is many studies[4, 12, 13], where, first, f is determined using
equation (9) either by measuring both, the liquid and vapor compositions, or by measuring the
liquid composition and using an activity model to find the vapor composition, then equations (7)
and (8) are used to get the T-x-y data of a system for a given pressure.
The flask ebulliometer (Fig 5.4) uses a small amount of sample, is the simplest among the
ones studied and easiest to operate. Also, it has the smallest amount of vapor space. Therefore
corrections due to evaporation would be negligible, if it were to be used to measure the bubble
pressures of binary mixtures. The correction due to the condensate hold-up would be relatively
small too, but important enough to be evaluated. However, neither the equilibrium liquid nor

81

vapor compositions can be measured easily with this ebulliometer. If a thermodynamically stable
and well known binary system is used, the factor f in equation (6) can be calculated by
measuring the T, P and qi only. The key is measuring these quantities accurately. We will test
this ebulliometer with methyl palmitate – methyl myristate system, for which P-T-x data are
available[14]. If f stays stable, the ebulliometer can be used to measure the P-T-x of unknown
systems.

82

REFERENCES

83

5.4 References
1.

Hala, E.; Pick, J.; Fried, V.; Vilim, O., "Vapor- Liquid Equilibrium". 1958, Pergammon
Press: New York.

2.

Swietoslawski, W., "Ebulliometric Measurements". 1945, Reinhold Publishing: New
York.

3.

Ambrose, D., "Vapor Pressures", in Experimental Thermodynamics. 1975, Butterworths:
London. 607-656.

4.

Rogalski, M.; Malanowski, S., "Ebulliometers Modified for the Accurate Determination
of Vapor-Liquid Equilibrium". Fluid Phase Equilibria, 1980, 5, 97-112.

5.

Chernyak, Y.; Clements, J.H., "Vapor Pressure and Liquid Heat Capacity of Alkylene
Carbonates". Journal of Chemical & Engineering Data, 2004, 49, 1180-1184.

6.

Raal, J.D.; Gadodia, V.; Ramjugernath, D.; Jalari, R., "New Developments in Differential
Ebulliometry: Experimental and Theoretical". Journal of Molecular Liquids, 2006, 125,
45 – 57.

7.

Cottrell, F.G., "On the Determination of Boiling Points of Solutions". Journal of the
American Chemical Society, 1919, 41, 721-729.

8.

Washburn, E.W.; Read, J.W., "“Concentrated” Solutions. IV. The General Boiling Point
Law". Journal of the American Chemical Society 1919, 41, 729-741.

9.

Hoover, S.R.; John, H.; Mellon, E.F., "Vapor Pressure Ebulliometer for Milliliter
Samples". Analytical Chemistry 1953, 25, 1940-1941.

10.

Olson, J.D., "Measurement of Vapor-Liquid Equilibria by Ebulliometry". Fluid Phase
Equilibria, 1989, 52, 209-218.

11.

Rogalski, M.; Rybakiewicz, K.; Malanowski, S., "Rapid and Accurate Method for
Determination of Vapour-Liquid Equilibrium". Berichte der Bunsen-Gesellschaft für
Physikalische Chemie, 1977, 81, 1070-1073.

84

12.

Zhenghong, G.; Zhonglong, M.; Qin, C.L.W., "Isothermal Vapor–Liquid Equilibria of
Benzene-Dimethylformamide and Toluene-Dimethylformamide Systems Using a Small
Modified Ebulliometer". Fluid Phase Equilibria, 2000, 173, 253–262.

13.

Kim, I.; Svendsen, H.F.; Børresen, E., "Ebulliometric Determination of Vapor-Liquid
Equilibria for Pure Water, Monoethanolamine, N-Methyldiethanolamine, 3(Methylamino)-Propylamine, and Their Binary and Ternary Solutions". Journal of
Chemical & Engineering Data, 2008, 53, 2521–2531.

14.

Rose, A.; Supina, W.R., "Vapor Pressure and Vapor-Liquid Equid Equilibrium Data for
Methyl Esters of the Common Saturated Normal Fatty Acids". Journal of Chemical &
Engineering Data, 1961, 6, 173-179.

85

CHAPTER VI
6. CONCLUSIONS AND RECOMMENDATIONS FOR FUTURE WORK

This chapter is divided in two sections. The first part deals with the StePPE simulations
and the second part is about ebulliometric measurements of vapor pressures.

6.1. VLE calculations with StePPE
1. Conclusions
Robust thermodynamic property prediction is at the core of developing good process
technologies in the chemical industry. Molecular simulations with UA Lennard-Jones
potential[1-4], have come a long way in accurately predicting equilibrium liquid densities, vapor
pressures, and heats of vaporization for pure compounds and mixtures. Recently, molecular
dynamics simulations with step potentials[5], in conjunction with TPT has proven to be a fast
alternative to continuous potential simulations. The step potential model, called SPEAD, though
reproduces vapor pressures relatively better than other methods, the liquid density predictions
lack accuracy. Vapor –liquid coexistence densities provide a stringent test for accuracy of a
force-field, since thermodynamic properties are sensitive to changes in densities. A new UA step
potential model, named StePPE, has been developed aimed at addressing the discrepancies in the
SPEAD model, so that it can be used not only for vapor pressures, but also for accurate
predictions of liquid densities. A new guideline has been outlined for selecting the site sizes,
based on Bondi’s[6] tabulated excluded volumes of the sites. Compared to SPEAD, the number
of adjustable parameters has been reduced in StePPE for the alkenes, alcohols, ketones and
esters. With the changes, a significant improvement in the liquid density predictions is achieved,
while the accuracy in vapor pressure predictions is better than macroscopic group contribution
86

methods, especially for some of the complex compounds studied. The improvements have also
helped ameliorate the transferability of the parameters to a broad range of compounds, and a
broad chain length range within a homologous series, as discussed in CHAPTER IV.
Comparisons have been made with the more detailed continuous potential models. A good
agreement is observed for the liquid densities compared to the accurate continuous potential
models, and the vapor pressure predictions are generally more accurate with StePPE.

2. Future work
Many of pure compounds like, alkynes, aldehydes, naphthenes, aromatics and amines
may not need any changes to the model parameters or the theory. The new parameters for those
compounds should be optimized first, as experimental data for many of them is readily available.
For the pure compounds reported in this study, though the StePPE model reproduces the
vapor pressures and liquid densities very well, some limitations were observed in the predictions
for branched compounds and esters. Some improvements can be made for the charge neutral
compounds by revisiting some of the model parameters, while improving the predictions of
compounds with charges might involve some improvement to the theory itself.

a) Branched compounds and site sizes
The agreement with the experimental liquid densities and vapor pressures for branched
compounds with StePPE is generally excellent for a broad range of compounds and with fewer
special parameters compared to some other models. However, a closer look at the density
predictions reveals a systematic overprediction for compounds with quaternary carbon (Fig 6.1).
One remedy would be to treat the side methyl groups attached to the quaternary carbon with

87

reduced well parameters. Alternatively, the size of the quaternary carbon can be increased. The
experimental liquid densities of branched alkanes decrease with branching, which suggests that
among structural isomers, the core volume should increase with branching. The current sizes
however, show an opposite trend, and the expected error is apparently compensated for by the
potential wells of the CH and CH3s groups.
500
StePPE
2,2-dimethylhexane

450

neopentane
400

T (K)

350
300
250
200
150
100
0.45

0.55

0.65

0.75

0.85

ρ (g/mL)
Fig 6.1: Liquid densities of neopentane and 2,2-dimethylhexane with StePPE. The symbols
are from DIPPR 801 correlations

A similar observation was made for the secondary alcohols. An optimization with an
increased size of the methine group was done. The results show a better agreement with the
increased CH size than what is obtained with the current set of parameters for the liquid densities
as well as the vapor pressures (Table 6.1). Additionally, unlike the current implementation, no
special parameters were used for 2-methanol.

88

From a fundamental viewpoint, increasing the site diameters appears more justifiable
than changing the well parameters to overcome the errors, which can be achieved by setting the
volume correction factor, as discussed in CHAPTER III, to above 20%. A reasonable value
would be 20.5%, which would affect the CH3 and CH2 diameters only by a small amount.

Table 6.1: Comparison of errors of secondary alcohols with two different sizes of the
methine group
sat

P

σCH = 3.9Å
Compound %AAD %Bias
2-propanol
8.68
8.60
2-butanol
3.20
3.07
2-pentanol
3.20
2.94
2-hexanol
7.10
-6.86
2-heptanol
7.68
-7.63
2-octanol
10.45
-9.73
2-nonanol
14.98
-14.98
3-pentanol
9.76
9.34
3-hexanol
9.73
-0.72
3-heptanol
9.17
-1.45

ρ

σCH = 3.9Å
%AAD %Bias
9.20
9.20
12.53
12.53
14.96
14.96
7.80
4.99
7.46
6.73
6.55
4.65
3.43
-0.07
4.69
-2.50
12.25
-10.64
13.67
-9.69

σCH = 3.9Å

σCH = 3.9Å

%AAD %Bias %AAD %Bias
1.41
-1.39
1.84
1.84
0.26
0.10
1.25
1.25
0.35
-0.25
0.79
0.77
0.70
0.70
1.64
1.64
0.27
0.21
0.94
0.94
0.55
0.55
1.25
1.25
0.46
0.46
1.04
1.04
0.31
-0.31
0.17
0.16
0.15
-0.02
0.37
0.34
0.45
0.42
0.75
0.75

b) Including the analytical tail correction
For a Lennard– Jones particle, the potential well extends to infinity, while the simulation
box is limited. The general practice therefore is to truncate the potential function at a cut-off
distance (rc), and add a tail correction to compensate for the truncation error, such that the
energy due to a particle can be given by,

E
NkT

=

ρ

rc

u ( r ) g (r ) 4πr
2kT ∫

2

dr +

0

ρ

∞

u (r ) g (r )4πr
2kT ∫
rc

89

2

dr

where, N is the total number of particles, ρ is the average number density of particles, g (r ) is
the radial distribution function (rdf) at a distance r , and u (r ) is the potential function. The
second integral is the tail correction, which can be calculated analytically. Assuming g (r ) = 1
beyond the cutoff separation, the contribution of the tail to the configurational energy per particle
can be given by,

ρ

∞

E tail
∫ u(r )4πr dr = NkT
2kT rc
2

Generally, rc is chosen to be around 2.5σ, where the potential is about 1.6% of the
minimum. Though seemingly negligible, the tail correction due to the potential beyond rc can be
significant at liquid-like densities[7]. For the step potentials, with rc = 2σ, used in SPEAD and
StePPE, the tail correction for a pair i is given by (for a complete derivation, see Appendix A.2),

Eitail
32 ρiσ i ε 4i
=−
NkT
3kT
3

where,

ε 4i is the fourth well, ρ i

is the number density of pair i and

σ i is the adjusted

diameter for the pair. The contribution from each pair can be added to the Helmholtz energy
equation from TPT.
To replicate a Lennard-Jones type potential more realistically, an optimization with tail
corrections for the alkanes was done, keeping the fourth well close to about 6% of the first,
(which is the potential at 2σ for Lennard- Jones potential). The preliminary results are shown in
Table 6.2. Though the overall agreements for the vapor pressures and the liquid densities are

90

good, the density deviations seem to increase with increasing molecular weight. This might need
a reoptimization of the CH2 site size. The analysis can then be extended to ether, ketones and
esters.

Table 6.2: Errors in liquid densities and vapor pressures of n-alkanes with StePPE with tail
corrections
Compound
propane
n-butane
n-octane
n-decane
n-tetradecane
Overall

sat

% AAD P
6.48
11.14
9.04
8.05
11.19
9.18

Optimized well parameters
Site
ε1 (K)
CH3
86.94
CH2
44

sat

% Bias P
-2.82
-7.87
4.26
3.41
1.46
1.69

%AAD ρ
1.03
0.71
0.53
1.42
2.45
1.34

%Bias ρ
-1.03
0.71
0.53
1.42
2.45
1.07

ε4 (K)
5.4
4.4

c) Intramolecular attraction parameters
Reduced diameters are used for non-bonded interactions between third neighbors in order
to impart flexibility to the chain. Using reduced diameters for up to six neighbors is not expected
to affect the simulation results noticeably. The non- bonded intramolecular attractions on the
other hand, are expected to affect the results more. In StePPE, the attractive interactions start
with neighbors separated by three bonds (collectively, we call it a ‘34’ interaction, meaning that
the attractive interactions start at the third bond and the full-diameter repulsions start at the
fourth bond). In order to explore this effect, ester simulations were run with ‘77’ intramolecular
interactions and the results were compared with those with the ‘34’ interactions (Fig 6.2). The
errors with ‘77’ interactions show an increasing trend with increasing number of carbon atoms in
the chain, while the trend is more scattered with the ‘34’ interactions, though the vapor pressures
91

appear to be overpredicted and the densities are underpredicted for the smallest esters with the
‘34’ interactions.

(%Bias) / (%Bias)Max

1.5
1.0
0.5
Psat 77
Psat 34
Density 77
Density 34

0.0
-0.5
-1.0
-1.5
0

5

10

15

20

NC
Fig 6.2: Vapor pressure and liquid density deviations with StePPE, using ‘34’ and ‘77’
intramolecular interactions; x-axis denotes the total number of carbon atoms in the ester
chain

Intuitively, increasing the number counts for the intramolecular attractions would help
reduce the errors in the densities and vapor pressures of the small esters, without affecting the
predictions of the bigger molecules as much. For a system of M molecules, with Ns number of
sites per molecule, the ratio r of the intramolecular attractive interactions to the total attractive
interaction can be given by (see Appendix A.3),
x −1

N s ( N s − 1) − 2∑ ( N s − i )
i =1
x −1

N s (MN s − 1) − 2∑ ( N s − i )
i =1

92

where, x is the bond distance where the intramolecular attractions begin to be counted. Currently

x is set at 3. Fig 5.3 shows a plot of the ratio given by the above equation with changing x for
different molecular sizes. The plot is not completely quantitative since the equation developed
above includes all pairs whether they fall inside or out of the pair-potential. In a real situation,
the ratio should scale down by some factor. The smaller molecules would be affected the most
by changing the position since they have steeper slopes. In order to increase the intramolecular
attractions, the ratio r has to be increased. But the current parameters are set at the minimum
bond distance possible for the intramolecular attractions. It is to be noted that the esters use the
optimized parameters from alkanes, ethers and ketones, all optimized with ‘34’ interactions. So,
the parameters of these functional groups would have to be optimized first with the desired value
of x, in order to get a more complete picture of the effect of x on the results.

0.008
0.007

NI / NT

0.006
20
15
10
7
5

0.005
0.004
0.003
0.002
0.001
0.000
3

4

5

6

XATT
Fig 6.3: Ratio of intramolecular attractive pairs (NI) to total number of attractive pairs
(NT) with x. The legends denote the total number sites in the molecule
d) Including electrostatic interactions in TPT

93

The pair-potential in StePPE is truncated at 2σ. The potential due to electrostatic
interactions dies off much slowly and is multiple orders of magnitude higher than the dispersion
potential in the range 0-2σ. For systems where partial charges are expected, it is therefore
important to include electrostatic interactions in the force-field. Since the current StePPE
implementation does not include electrostatic interactions, adjusting the step wells based on the
type of interaction can somewhat compensate the effect of the partial charges. For example, in
dimethyl ether, the oxygen is expected to have a partial negative charge and the CH3 group is
expected to have a partial positive charge. For non-bonded interactions between oxygen and a
CH3 group, the potential well can be made deeper to compensate for the attraction between these
two sites, and likewise, the potential wells for the oxygen-oxygen, or CH3-CH3 interactions can
be made shallower as shown below,
0
ε ij ,m = ε ij ,m − K

where,

0
ε ij ,m

qi q j

σ ij ,m

is the original well depth of well m of pair i-j,

q i , q j are the charges, σ ij ,m

denotes the location of well m for pair i-j and K is a scaling factor.
Figure 6.4 shows a comparison of results for ethers based on this approach with the
original StePPE results. Charges were borrowed from the OPLS force-field[8]. These results are
preliminary and more detailed study needs to be done in order to incorporate the effect of
charges in the code. For better results, the factor K in the above equation should be made an
adjustable parameter or a more sophisticated method should be developed. Parameters obtained
from ethers and ketones can then be used for esters in order to test the transferability.

94

25

K=0
K = 30

%AAD Psat

20

15

10

5

C
9

C
9O

C
8

C
6

C
8O

C
5

C
6O

C
5O

C
4

C
6

C
4O

C
2O

C
2

C
4O

C
3

C
5

C
3O

C
1O

C
3

C
4

C
2O

C
2

C
1O

C
2O

C
3

C
1O

C
1O

C
1

0

Figure 5.4: Vapor pressures of ethers with StePPE with and without correction for partial
charges

e) Radius of Gyration
Another interesting study would involve comparison of the radius of gyration for StePPE
simulations with the reported radius of gyration. Below are some comparisons of the radius of
gyration reported by DIPPR to the most likely end-to-end distance from the intramolecular
distribution function for the n-alkanes. The values appear to be consistently lower with StePPE.
This is however a crude method of obtaining the radius of gyration using only the end-to-end
distance, and more precise method is needed to calculate it during the simulations.

95

Table 6.3: Comparison of the radius of gyration of n-alkanes with the most likely end-toend distance in StePPE simulations
RG (Å)
Compound StePPE DIPPR
butane
1.93
2.88
pentane
2.15
3.34
hexane
2.7
3.77
octane
3.69
4.55
nonane
4.68
4.81
decane
4.13
5.15

6.2 Ebulliometric measurements

1. Conclusions and recommendations
To measure vapor pressures of pure components and bubble pressures of binary mixtures,
ebulliometry provides a fast and reasonably accurate method, compared to static methods.
Commercially available ebulliometers are generally complicated and require sample sizes in the
range of about 150 – 200 mL[9, 10]. Several other ebulliometers have been used for laboratory
measurements with smaller samples[11].
For the measurement of pure component vapor pressures, our designs with a Cottrell[12]
or Washburn[13] type pump were found to be effective in removing superheat. Similar accuracy
was achieved with our stirred flask type design, which is the simplest among the three. The
design is a scaled-down form of the stirred-type ebulliometer discussed in the ASTM E1719
method. The heated surface of the flask was activated with powdered glass to facilitate the
nucleation of bubbles. The activation coupled with vigorous stirring was found to be very
effective in removing the superheat from the boiling liquid, and temperature measurement was

96

possible by putting the probe in the liquid itself rather than in a separate equilibrium chamber, as
in Cottrell’s or Washburn’s design.
Based on an analysis similar to the stirred flask design for binary mixtures, the
condensate hold –up in the Washburn apparatus can be estimated. And used for the measurement
of bubble pressures of unknown systems.

97

REFERENCES

98

6.3 References
1.

Siepmann, J.I.; Karaborni, S.; Smit, B., "Simulating the Critical Behaviour of Complex
Fluids". Nature 1993, 365, 330- 332.

2.

Jorgensen, W.L.; Madura, J.D.; Swenson, C.J., "Optimized Intermolecular Potential
Functions for Liquid Hydrocarbons". Journal of the American Chemical Society, 1984,
106, 6638-6646.

3.

Martin, M.G.; Siepmann, J.I., "Transferable Potentials for Phase Equilibria. 1. UnitedAtom Description of n-Alkanes". Journal of Physical Chemistry B, 1998, 102, 25692577.

4.

Nath, S.K.; Escobedo , F.A.; de Pablo, J.J., "On the Simulation of Vapor–Liquid
Equilibria for Alkanes". Journal of Chemical Physics, 1998, 108, 9905-9911.

5.

Unlu, O.; Gray, N.H.; Gerek, Z.N.; Elliott, J.R., "Transferable Step Potentials for the
Straight-Chain Alkanes, Alkenes, Alkynes, Ethers, and Alcohols". Industrial &
Engineering Chemistry Research 2004, 43, 1788-1793.

6.

Bondi, A., "van der Waal's Volumes and Radii". The Journal of Physical Chemistry,
1964, 68, 441-451.

7.

Frenkel, D.; Smit, B., "Understanding Molecular Simulation". 2nd ed. 1996, Academic
Press: San Diego.

8.

Briggs, J.M.; Matsui, T.; Jorgensen, W.L., "Monte Carlo Simulations of Liquid Alkyl
Ethers with the OPLS Potential Functions". Journal of Computational Chemistry, 1990,
11, 958–971.

9.

Olson, J.D., "Measurement of Vapor-Liquid Equilibria by Ebulliometry". Fluid Phase
Equilibria, 1989, 52, 209-218.

10.

ASTM Standard E1719, Standard Test Method for Vapor Pressure of Liquids by
Ebulliometry. 2005, www.astm.org: West Conshohocken, PA.

11.

Swietoslawski, W., "Ebulliometric Measurements". 1945, Reinhold Publishing: New
York.

99

12.

Cottrell, F.G., "On the Determination of Boiling Points of Solutions". Journal of the
American Chemical Society, 1919, 41, 721-729.

13.

Washburn, E.W.; Read, J.W., "“Concentrated” Solutions. IV. The General Boiling Point
Law". Journal of the American Chemical Society 1919, 41, 729-741.

100

APPENDICES

101

APPENDICES

A.1 Angle distribution in butyl propionate
To demonstrate that the bonds in the DMD are flexible, the following figure shows the
bond distribution from 10 snapshots from the simulation of butyl propionate. The curves are not
smooth due to the small sample.

0.08

CH2-(CH2)-CH2
=O-(>C=)-O
CH2-(>C=)-O=)
>C=(-O-)-CH2
CH2-(>C=)-O

0.07
0.06
0.05
0.04
0.03
0.02
0.01
0
80

90

100

110

120

130

140

Angle (degrees)

Figure A.1: Distribution of some angles in butyl propionate. The legends denote the three
sites that form the angle, and the site within the parentheses is the central angle. The
vertical lines denote the equilibrium values (109.5° and 120°)

102

A2. Bond lengths
Table A2.1 Bond lengths used in StePPE model
Length(Å) Remarks
1.54
Any bond in a saturated alkane
1.54
Bond between carbonyl carbon in an
ester or ketone and a CH3 or a CH2
group
>C=(ester)
O-(ester)
1.43
>C=(ester/ketone) O=(ester/ketone) 1.21
CH
C/CH/CH2/CH3 1.54
CH=
CH=/CH=c/CH=t 1.54
CH=c and CH=t denote a CH= in a
cis or trans alkene
CH=c
CH=c
1.34
CH=t
CH=t
1.34
CH=
CH2
1.54
CH=
CH2=
1.34
CH(OHs)
CH2
1.54
CH=c
CH2
1.54
CH=t
CH2
1.54
CH(OHs)
CH3
1.54
CH(OHs) denotes a CH in a
secondary alcohol
CH=c
CH3
1.54
CH=t
CH3
1.54
CH=
CH3
1.54
CH(OHs)
OHs
1.43
CH2
C/CH/CH2/CH3 1.54
CH2(OH)
CH2
1.54
CH2(r)
CH2(r)
1.54
Bonds between CH2 groups in ring
compounds
CH2e(OH)
CH3
1.54
CH2
O-(ester/ether)
1.43
CH2e(OH)
OH
1.43
Only in ethanol
CH2(OH)
OH
1.43
CH3
C/CH/CH2/CH3 1.54
CH3
O-(ether/ester)
1.43
CH3(OH)
OH
1.43
Only in methanol
SiteA
SiteB
C
C/CH/CH2/CH3
>C=(ester/ketone) CH2/CH3

103

A3. Finding the first guess of site sizes from excluded volumes

R2

R1
m
l
h1

h2

Fig A3.1: Schematic for finding site diameters from given excluded volumes

In Fig A3.1, we use a propane molecule to illustrate the calculation of the volume
correction factor for each covalent bond. The outer two spheres are CH3 groups and the middle
sphere is CH2. The excluded volumes of these sites and the bond length (l in the figure above)
are known. With use of simple geometry, R1 and R2 can be calculated. Once the diameters are
known, the volume lost by each site due to overlap can be calculated. The method is described
below.
The volume of a section of a sphere of height h and radius R can be given by:

V =

π
3

h 2 (3R − h

)

We can use the auxiliary parameter m to calculate R1 and R2 from the given excluded volumes

V1 and V2 using the following equations:

104

2
R12 + l 2 − R2
m=
2l
h1 = R1 − m

h2 = l + R2 − m

π
4
V1 = πR13 − 2 h12 (3R1 − h1 )
3
3
V2 =

π
3

h22 (3R2 − h2 )

105

A.4. Liquid density and vapor pressure plots
In all the plots in this section, the experimental data is shown as points and the StePPE
predictions are shown as solid lines. The lines and symbols are of the same color for a molecule.
700

nC2
nC3
nC4
nC5
nC6
nC8
nC9
nC10
nC12
nC13
nC14
nC16
nC18
nC19
nC22
nC24

(a)
600

500

T(K)

400

300

200

100

0
0.35

0.4

0.45

0.5

0.55

0.6

0.65

0.7

0.75

0.8

rho(g/mL)

10

(b)

1
0.1

Psat(MPa)

0.01
0.001
0.0001
1e-05
1e-06
1e-07

nC2
nC3
nC4
nC5
nC6
nC8
nC9
nC10
nC12
nC13
nC14
nC16
nC18
nC19
nC22
nC24

1e-08
1e-09
0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009

0.01

1/T(K)

Figure A.4. 1: Liquid densities (a) and vapor pressures (b) of n-alkanes

106

0.011

550

(a)
500
450
400

T(K)

350
300
250
200
150
100
0.45

0.5

0.55

0.6

0.65

0.7

0.75

rho(g/mL)
C3-2.2-diMe
C4-2-Me
C4-2.2.3-triMe
C4-2.3-diMe
C5-2-Me
C5-3-Me
C5-3-Et
C5-2.2.4-triMe

C5-2.3.4-triMe
C6-2-Me
C6-3-Me
C6-2.5-diMe
C6-2.4-diMe
C6-2.2-diMe
C6-3.4-diMe
C6-3-Et

380

C7-2-Me
C7-4-Me
C7-2.6-diMe
C7-3-Et
C9-2-Me
C9-5-Me

isobutane

360

(b)

340

T(K)

320
300
280
260
240
220
200
0.44 0.48 0.52 0.56

0.6

0.64 0.68

rho(g/mL)

Fig A.4.2: Liquid densities of branched alkanes (a); isobutane (b)

107

0.8

0.85

10

(a)
1
0.1
0.01

Psat(MPa)

0.001
0.0001
1e-05
1e-06
1e-07
1e-08
1e-09
0.001

0.002

0.003

0.004

0.005

0.006

0.007

1/T(K)
C3-2.2-diMe
C4-2-Me
C4-2.2.3-triMe
C4-2.3-diMe
C5-2-Me
C5-3-Me
C5-3-Et
C5-2.2.4-triMe

C5-2.3.4-triMe
C6-2-Me
C6-3-Me
C6-2.5-diMe
C6-2.4-diMe
C6-2.2-diMe
C6-3.4-diMe
C6-3-Et

C7-2-Me
C7-4-Me
C7-2.6-diMe
C7-3-Et
C9-2-Me
C9-5-Me

10

(b)

Psat(MPa)

1
0.1
0.01
0.001
0.0001

1e-05
0.002 0.003 0.004 0.005 0.006 0.007
1/T(K)

Fig A.4.3: Vapor pressures of branched alkanes (a); isobutane (b)

108

0.008

550

C3-1-ene
C4-1-ene
C4-2c-ene
C4-2t-ene
C4-1.3-diene
C5-1-ene
C5-2t-ene
C5-1.4-diene
C6-1-ene
C6-1.5-diene
C7-1-ene
C7-2c-ene
C8-1-ene
C8-2t-ene
C10-1.9-diene

(a)
500
450

T(K)

400
350
300
250
200
150
0.45

0.5

0.55

0.6

0.65

0.7

0.75

0.8

0.85

rho(g/mL)
10

C3-1-ene
C4-1-ene
C4-2c-ene
C4-2t-ene
C4-1.3-diene
C5-1-ene
C5-2t-ene
C5-1.4-diene
C6-1-ene
C6-1.5-diene
C7-1-ene
C7-2c-ene
C8-1-ene
C8-2t-ene
C10-1.9-diene

(b)

1
0.1

Psat(MPa)

0.01
0.001
0.0001
1e-05
1e-06
1e-07
1e-08
1e-09
0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009

0.01

1/T(K)

Figure A.4.4: Liquid densities (a) and vapor pressures (b) of alkenes

109

600

1olC3
1olC4
1olC5
1olC7
1olC8
1olC9
1olC10
1olC12
1olC13
1olC15
1olC16
1olC17
1olC18
1olC19
1olC20

(b)
550
500

T(K)

450
400
350
300
250
200
150
0.55

0.6

0.65

0.7

0.75

0.8

0.85

0.9

0.95

rho(g/mL)

10

1olC3
1olC4
1olC5
1olC7
1olC8
1olC9
1olC10
1olC12
1olC13
1olC15
1olC16
1olC17
1olC18
1olC19
1olC20

1
0.1

Psat(MPa)

0.01
0.001
0.0001
1e-05
1e-06
1e-07
0.001

0.0015

0.002

0.0025

0.003

0.0035

0.004

1/T(K)

Figure A.4.5: Liquid densities (a) and vapor pressures (b) of primary alcohols. The symbol
‘1olCx’ denotes a primary alcohol with ‘x’ carbons in the backbone. For example, 1olC8 in
1-octanol

110

500

C4-2ol
C5-2ol
C5-3ol
C6-2ol
C6-3ol
C7-2ol
C7-3ol
C8-2ol
C9-2ol

(a)
450

T(K)

400

350

300

250

200
0.55

0.6

0.65

0.7

0.75

0.8

0.85

0.9

rho(g/mL)
10

(b)
1
0.1

Psat(MPa)

0.01

C4-2ol
C5-2ol
C5-3ol
C6-2ol
C6-3ol
C7-2ol
C7-3ol
C8-2ol
C9-2ol

0.001
0.0001
1e-05
1e-06
1e-07
0.0015

0.002

0.0025

0.003

0.0035

0.004

0.0045

1/T(K)

Figure A.4.6: Liquid densities (a) and vapor pressures (b) of secondary alcohols

111

600

(a)
550
500

T(K)

450
400
350

C1OC1
C1OC3
C1OC5
C2OC2
C2OC3
C2OC6
C3OC3
C4OC2
C4OC4
C5OC5
C6OC6
C8OC8
C9OC9

300
250
200
0.5

0.55

0.6

0.65

0.7

0.75

0.8

0.85

rho(g/mL)
10

(b)

1
0.1

Psat(MPa)

0.01
0.001
0.0001
1e-05
1e-06

C1OC1
C1OC3
C1OC5
C2OC2
C2OC3
C2OC6
C3OC3
C4OC2
C4OC4
C5OC5
C6OC6
C8OC8
C9OC9

1e-07
1e-08
1e-09
0.0015 0.002 0.0025 0.003 0.0035 0.004 0.0045 0.005 0.0055 0.006
1/T(K)

Figure A.4.7: Liquid densities (a) and vapor pressures (b) of ethers

112

450

C4.2one
C5.2one
C6.2one
C7.2one
C8.2one
C9.2one
C5.3one
C6.3one
C7.3one
C7.4one
C8.4one
C9.4one
C9.5one

(a)
400

T(K)

350

300

250

200
0.68

0.7

0.72

0.74

0.76

0.78

0.8

0.82

0.84

0.86

0.88

0.9

rho(g/mL)
10

(b)
1

Psat(MPa)

0.1

0.01

0.001

C4.2one
C5.2one
C6.2one
C7.2one
C8.2one
C9.2one
C5.3one
C6.3one
C7.3one
C7.4one
C8.4one
C9.4one
C9.5one

0.0001

1e-05

1e-06
0.0015

0.002

0.0025

0.003

0.0035

0.004

0.0045

0.005

1/T(K)

Figure A.4.8: Liquid densities (a) and vapor pressures (b) of ketones

113

0.0055

500

C2ateC1
C2ateC2
C2ateC3
C2ateC4
C2ateC5
C2ateC6
C2ateC7
C2ateC8
C2ateC10

(a)
450

400

T(K)

350

300

250

200

150
0.6

0.65

0.7

0.75

0.8

0.85

0.9

0.95

1

1.05

1.1

rho(g/mL)
10

(b)
1

Psat(MPa)

0.1

C2ateC1
C2ateC2
C2ateC3
C2ateC4
C2ateC5
C2ateC6
C2ateC7
C2ateC8
C2ateC10

0.01

0.001

0.0001

1e-05
0.002

0.0025

0.003

0.0035

0.004

1/T(K)

Figure A.4.9: Liquid densities (a) and vapor pressures (b) of acetates

114

0.0045

380

C3ateC1
C4ateC1
C6ateC1
C8ateC1
C10ateC1
C12ateC1
C14ateC1
C16ateC1
C18ateC1

(a)
360

T(K)

340

320

300

280

260
0.78

0.8

0.82

0.84

0.86

0.88

0.9

0.92

0.94

rho(g/mL)
10

(b)
1
0.1

Psat(MPa)

0.01
0.001

C3ateC1
C4ateC1
C6ateC1
C8ateC1
C10ateC1
C12ateC1
C14ateC1
C16ateC1
C18ateC1

0.0001
1e-05
1e-06
1e-07
1e-08
0.0015

0.002

0.0025

0.003

0.0035

0.004

1/T(K)

Figure A.4.10: Liquid densities (a) and vapor pressures (b) of methyl esters

115

500

C3ateC2
C3ateC3
C4ateC2
C3ateC4
C4ateC3
C4ateC4
C5ateC4
C8ateC2

(a)
450

T(K)

400

350

300

250
0.6

0.65

0.7

0.75

0.8

0.85

0.9

0.95

rho(g/mL)
10

(b)
1

Psat(MPa)

0.1

0.01

C3ateC2
C3ateC3
C4ateC2
C3ateC4
C4ateC3
C4ateC4
C5ateC4
C8ateC2

0.001

0.0001

1e-05

1e-06
0.0015

0.002

0.0025

0.003

0.0035

0.004

1/T(K)

Figure A.4.11: Liquid densities (a) and vapor pressures (b) of esters higher than acetates

116

A.5. Tail corrections
The average potential energy due to particle i is given by (from the text “Understanding
Molecular Simulation” by Frenkel and Smit):

E
NkT

=

∞

ρ
2kT

u (r ) g (r )4πr 2 dr
∫

(A5.1)

0

where, ρ is the average number density of particles, g (r ) is the rdf at a distance r , u (r ) is the
potential function and N is the total number of particles. For a L-J particle, the potential well
extends to infinity, while the simulation box is limited. The general practice is to truncate the
potential function at a cut-off distance (rc), and add a tail correction to compensate for the
truncation error, such that:

E
NkT

=

ρ
2kT

rc

∫ u(r ) g (r )4πr dr +
2

0

ρ
2kT

∞

u ( r ) g (r )4πr 2 dr
∫

rc

(A5.2)
The second integral in equation A5.2 is the tail correction, which, on assuming g (r ) = 1 beyond
the cutoff separation, can be written as:

∞

E tail
ρ
=
u (r )4πr 2 dr
∫
NkT 2kT rc

(A5.3)

Generally, rc is chosen to be around 2.5σ, where the potential is about 1.6% of the minimum.
Though seemingly negligible, the tail correction due to the potential beyond rc can be significant
at liquid-like densities.

117

In SPEAD, Helmholtz energy is calculated using the configurational energy of the
system, using a cut-off distance of 2.0 σ. Tail corrections have not been accounted for in the
previous SPEAD publications, and the results nevertheless, appear to be reasonably good,
probably because the fourth well is an adjustable parameter, which in most cases is deeper than
about 6% (as in a Lennard-Jones potential at 2.0 σ) of the first well.
The integral for the tail correction is solved analytically, by using an approximation of
the attractive potential beyond 2σ. Ignoring the repulsive contribution, the potential can be
6

assumed to be of the form –B/r , and knowing that, at rc = 2σ, the potential is ε4, B is found to
6

be 64ε4σ , which would give:

64ε 4σ 6
E tail
ρ ∞
4πr 2
=−
dr
6
∫
NkT
kT rc
r

(A5.4)

Since we are calculating the Helmholtz energy by summing over pair type, the factor 2 is
dropped from the denominator. Integrating, we find for a particular pair type i:

tail
i

E
=−
NkT

32 ρ iσ i ε 4i
3

3kT

(A5.5)

The Helmholtz energy term can be modified as follows to include the tail correction:

A− A
NkT

ig

E+E
A −A
= 0
−
NkT
NkT
ig

tail

(

)

 E + E tail 2 − E + E tail

−
2 Nk 2T 2

where, E is the potential energy due to the wells. Since:

118

2





(A5.6)

E tail =

Eitail
∑
pairs

which is constant, it cancels out in the energy fluctuation term, and can simply be added to the
energy term such that:

E + E tail
NkT

N i πσ i3ε 4i
A1 − ∑
3kV
pairs
=
T

(A5.7)

Where, Ni is the number of pairs of types i and V is the volume of the simulation box.

A.6. Ratio of the number of intramolecular attractive pairs to total attractive pairs
This analysis is only for straight chain molecules. The derivation would be more tedious
for branched compounds.
Consider the following system:
Total number of molecules = M
Number of sites in a molecule = Ns
Therefore, total intramolecular interactions in a molecule =

N sC2 =

Ns!
N (N − 1)
= s s
2( N s − 2)!
2

Similarly, the total number of attractive pairs in the simulation box can be given by:

MN s (MN s − 1)
2

119

Let the position where intramolecular attractions are counted from = x, which also denotes the
minimum bond separation for the interaction to be counted. Therefore, number of intramolecular
interaction to be ignored for attractions per molecule can be given by:
x −1

(Ns -1) + (Ns-2) +······+ (Ns-(x-1)) =

∑ (N

s

− i)

i =1

The first term denotes number of pair interactions separated by 1 bond, the second term, by 2
bonds and so on. The last term is for the interactions separated by (x-1) bonds, since the
interaction separated by x bonds or more are included. So, the total intramolecular attractive
pairs in the simulation box (NI) can be given by:

 N s ( N s − 1) x −1

NI = M 
− ∑ ( N s − i )
2
i =1


Similarly, the total number of attractive pairs in the simulation box (NT) can be given by:

 N s (MN s − 1) x −1

− ∑ ( N s − i )
NT = M 
2
i =1


Therefore the ratio NI/NT can be given by:

 N s (MN s − 1) x −1

− ∑ ( N s − i)

2
NI
i =1


=
N T  N s (MN s − 1) x −1

− ∑ ( N s − i)

2
i =1



120