7 Computer Simulation of the Liquid-Vapor Surface of Molecular Fluids S. M. THOMPSON and K. E. GUBBINS
Downloaded by PURDUE UNIV on July 7, 2016 | http://pubs.acs.org Publication Date: June 1, 1978 | doi: 10.1021/bk-1978-0086.ch007
School of Chemical Engineering, Cornell University, Ithaca, NY 14853
The r e s u l t s described here are from the initial stages o f a program designed to study the s t r u c t u r e and p r o p e r t i e s of liquidvapor surfaces of systems composed o f polyatomic s p e c i e s . Both computer s i m u l a t i o n and p e r t u r b a t i o n theory methods are the t o o l s used to o b t a i n these r e s u l t s . Here we report molecular dynamics (MD) computer simulations o f two homonuclear diatomic l i q u i d s , each at a temperature in the r e g i o n of the triple p o i n t . The l i q u i d - v a p o r s u r f a c e of systems of monatomic molecules has been the subject o f considerable i n v e s t i g a t i o n (1-7) i n the past few years, but as f a r as we are aware, t h i s is the first work of t h i s type on molecular s p e c i e s . This work is initially d i r e c t e d towards o b t a i n i n g both e q u i l i b r i u m and non-equilibrium p r o p e r t i e s in the inhomogeneous s u r f a c e zone f o r one component systems. E q u i l i b r i u m p r o p e r t i e s i n c l u d e the d e n s i t y - o r i e n t a t i o n profile, which provides i n f o r mation on p r e f e r r e d o r i e n t a t i o n s ( i f any) i n the surface zone, and s u r f a c e tensions and energies. Non-equilibrium p r o p e r t i e s i n c l u d e t r a n s l a t i o n a l and r e - o r i e n t a t i o n a l v e l o c i t y a u t o c o r r e l a t i o n functions and t h e i r a s s o c i a t e d memory f u n c t i o n s , l e a d i n g to information on the d i f f u s i o n o f molecules both perpendicular to and p a r a l l e l to the plane o f the i n t e r f a c e . Here only e q u i l i b rium p r o p e r t i e s are presented. Future extensions i n c l u d e the study o f binary mixtures of molecules of v a r y i n g complexities and the behaviour of r e l a t i v e l y massy s u r f a c t a n t molecules i n the surface region. I.
The Simulation
Each system c o n s i s t s o f 216 molecules confined to square prisms o f dimensions x , V L ( = X ) and Z ( > X L ) with the usual p e r i o d i c boundary c o n d i t i o n s i n the (xy) plane of the s u r f a c e . The l i q u i d was confined to the centre o f the c e l l with vapor phases a t e i t h e r end. The two surfaces are s t a b l e without the i n c l u s i o n o f e x t e r n a l forces (2,_3>4_»j3) . I n i t i a l l y the center of mass i s f i x e d i n the center ( z = z / 2 ) o f the c e l l . The L
l
L
C 0 M
L
0-8412-0463-2/78/47-086-076$05.00/0 © 1978 American Chemical Society Lykos; Computer Modeling of Matter ACS Symposium Series; American Chemical Society: Washington, DC, 1978.
7.
THOMPSON AND
Liquid-V apor Surface of Molecular Fluids
GUBBINS
77
r
t r a n s l a t i o n a l and angular v e l o c i t i e s o f the molecules a r e r a n domly chosen subject to the c o n s t r a i n t s that the centre o f mass has no r e s u l t a n t t r a n s l a t i o n a l or angular v e l o c i t y and that the d e s i r e d k i n e t i c energy i s obtained, t h i s being d i v i d e d between t r a n s l a t i o n and r o t a t i o n a l according to the e q u i p a r t i t i o n law. Random o r i e n t a t i o n s were assigned. The molecular centers o f mass are assigned coordinates on an f . c . c . l a t t i c e s t r u c t u r e o f the appropriate l i q u i d d e n s i t y , while the 'vapor phases a r e 'empty . As the s i m u l a t i o n proceeds, the l a t t i c e s t r u c t u r e melts and a vapor phase develops. When the r e s u l t a n t density p r o f i l e has become s t a b l e , t h i s i n i t i a l e q u i l i b r a t i o n period ( t y p i c a l l y 10* i n t e g r a t i o n steps) i s r e j e c t e d and the run proper begins. A vapor molecule which attempts to e x i t the c e l l by means o f one of the end w a l l s i s returned to the c e l l by simply bouncing i t o f f the w a l l . The low vapor d e n s i t y ensures that t h i s procedure r e s u l t s i n t o t a l l y n e g l i g i b l e d i s t o r t i o n s to the surface s t r u c ture. The c e l l s were each 19 molecular diameters (a) i n height ( Z L = 19a), and the width o f the c e l l (XL) was computed as twice the p o t e n t i a l c u t - o f f d i s t a n c e . This l a t t e r quantity was 2.5 molecular diameters plus the molecular elongation (see below). Thus the c e l l widths were 5.6584 and 6.2160 molecular diameters f o r N and C l r e s p e c t i v e l y . Each f l u i d was modelled by means o f a s i t e - s i t e or atomatom p o t e n t i a l (8). Each molecule i s assumed to c o n s i s t o f two i n t e r a c t i o n centers (commonly assumed p o s i t i o n a l l y c o i n c i d e n t with the atomic n u c l e i ) . The intermolecular p o t e n t i a l i s then the sum o f four s h i f t e d Lennard-Jones (12,6) i n t e r a c t i o n s (see f i g u r e 1) 1
Downloaded by PURDUE UNIV on July 7, 2016 | http://pubs.acs.org Publication Date: June 1, 1978 | doi: 10.1021/bk-1978-0086.ch007
1
2
2
4 u(r u)ia>2) = 12
I ^JS^K^ K l
^
=
where \jS
- \j
M
M
~ LJ u
( r
c>
+
U
i j
(
r
c
)
(
r
c
"
r
)
2
and ^ ( r )
- 4e[(a/r)
x a
6
- (a/r) ]
(3)
Where r i s the p o t e n t i a l c u t - o f f d i s t a n c e , the prime denoting d i f f e r e n t i a t i o n , z i s the w e l l depth and a the molecular c o l l i s i o n diameter o f the Lennard-Jones i n t e r a c t i o n , co^ = { S ^ ^ } are the p o l a r angles s p e c i f y i n g the o r i e n t a t i o n o f molecule i , and = i - 2- The f i n a l term i n (2) ensures that the p o t e n t i a l and i t s d e r i v a t i v e a r e continuous a t the c u t - o f f d i s t a n c e , l e a d i n g to b e t t e r dynamics. The four atom-atom distances are given by: c
12
Lykos; Computer Modeling of Matter ACS Symposium Series; American Chemical Society: Washington, DC, 1978.
78
COMPUTER
r! - r f 2
r
2
= r
i
3
= ri
2
X
2
12
2
2r
2
2
= r
MATTER
2
+ £| + 2 r ( £ i C O S 6 i + £ cos 6 ) + 2*ifc f(wiw )
2
2
OF
2*f + 2£iri (cos 6 i - cos 0 ) - 2£if(u)i )
2
2
r
2
MODELING
2
i 2
2
2
2
2
(£ cos 9 i + £iCos 0 ) + 2£i£ f (o)ia) ) 2
2
2
2
2£| - 2 r i £ ( c o s 0 i - cos 0 ) - 2fc|f (u>i(i>) 2
2
2
2
(4) with f(u)iU) ) = cos 0 i cos 0 2
+ s i n 0i s i n 0
2
2
cos at 172.0 K po
Lykos; Computer Modeling of Matter ACS Symposium Series; American Chemical Society: Washington, DC, 1978.
84
COMPUTER MODELING
Table I
Averages Obtained
During
OF
MATTER
Simulations
N2_
Cl
2
40.00 x 10
3
Number of time steps
42.75
Temperature, K
66.52 + 2.01
172.0
0.608 + 0.024
0.526 ± 0.023
Surface thickness , d/a
1.79
1.14
L o c a t i o n o f Gibbs d i v i d i n g surface, z / a
4.20
4.00
Surface t e n s i o n ,
11.98
L i q u i d d e n s i t y , p*^
s
i
m
X
10
3
± 5.0
Q
Downloaded by PURDUE UNIV on July 7, 2016 | http://pubs.acs.org Publication Date: June 1, 1978 | doi: 10.1021/bk-1978-0086.ch007
JnT
x
2
10
+ 0.44
37.35 ± 0.83
+ 3
Experimental c o e x i s t i n g l i q u i d density, P* (12) L,exp —
0.664 a t 66.5 K
0.55 at 172.0 K
Experimental s u r f a c e t e n s i o n , Jm" x 1 0 (13)
11.41 at 66.5 K
39.06 at 172.0 K
2
3
z ^ 4a and z ^ 5a i s genuine s t r u c t u r e , being free o f these features. The corresponding curve f o r N i s not p l o t t e d because of poor s t a t i s t i c s f o r t h i s n e a r - s p h e r i c a l molecule. A pro nounced tendency to adopt p r e f e r r e d o r i e n t a t i o n s i s i n d i c a t e d , t h i s tendency being height dependent. In the l i q u i d phase at z - z = 0.85a (corresponding to the maximum i n p ( z ) ) the mole cules have a tendency to o r i e n t with t h e i r axes v e r t i c a l ( 9 = 0 ) , while at the Gibbs s u r f a c e ( c l o s e to the minimum i n p ( z ) ) a r e v e r s a l takes place and the molecules bend to o r i e n t with t h e i r axes p a r a l l e l to the i n t e r f a c i a l plane. This r e s u l t i s q u a l i t a t i v e l y i n agreement with the p r e d i c t i o n s o f f i r s t order p e r t u r b a t i o n theory f o r a n i s o t r o p i c overlap forces (14). A p e r t u r b a t i o n treatment using the f u l l atom-atom p o t e n t i a l i s under way. The p a i r p o t e n t i a l i s c u r r e n t l y being modified by the i n c l u s i o n of a quadrupole-quadrupole p o t e n t i a l , and the e f f e c t of s u r f a c e area on s u r f a c e thickness found i n (7) f o r monatomic molecules i s being i n v e s t i g a t e d . 2
Q
2
2
Acknowledgment It i s a pleasure to thank the N a t i o n a l Science Foundation and the Petroleum Research Fund, administered by the American Chemical S o c i e t y , f o r grants i n support of t h i s work, and S o h a i l Murad f o r a copy of h i s s i m u l a t i o n program f o r hydrogen c h l o r i d e .
Lykos; Computer Modeling of Matter ACS Symposium Series; American Chemical Society: Washington, DC, 1978.
7. THOMPSON AND GUBBiNs
Liquid-Vapor Surface of Molecular85Fluids
Abstract An application of the molecular dynamics method to simulate the liquid-vapor surface of molecular fluids is described. A predictor-corrector algorithm is used to solve the equations of translational and rotational motion, where the orientations of molecules are expressed in quaternions. The method is illustrat ed with simulations of 216 homonuclear (N and Cl ) diatomic molecules. Properties calculated include surface tensions and density-orientation profiles. 2
2
Downloaded by PURDUE UNIV on July 7, 2016 | http://pubs.acs.org Publication Date: June 1, 1978 | doi: 10.1021/bk-1978-0086.ch007
Literature Cited 1. Chapela, G. Α., Saville, G. and Rowlinson, J. S., Faraday Disc. Chem. Soc. (1975) 59, 22. 2. Lee, J. Κ., Barker, J. A. and Pound, G. M., J. Chem. Phys. (1974) 60, 1976. 3. Abraham, F. F., Schreiber, D. E. and Barker, J. Α., J. Chem. Phys. (1975) 62, 1958. 4. Opitz, A. C. L., Phys. Letters A (1974) 47, 439. 5. Liu, K. S., J. Chem. Phys. (1974) 60, 4226. 6. Rao, M. and Levesque, D., J. Chem. Phys. (1976) 65, 3233. 7. Chapela, G. Α., Saville, G., Thompson, S. M. and Rowlinson, J. S., Trans. Far. Soc. II (1977) 73, 1133. 8. Sweet, J. R. and Steele, W. Α., J. Chem. Phys. (1967) 47, 3029. 9. Cheung, P. S. Y. and Powles, J. G., Mol. Phys. (1975) 30, 921. 10. Singer, Κ., Taylor, A. and Singer, J. V. L., Mol. Phys. (1977) 33, 1757. 11. Evans, D. J. and Murad, S., Mol. Phys. (1977) 34, 327. 12. Vargaftik, Ν. Β., Tables on the Thermodphysical Properties of Liquids and Gases, 2nd ed., Halsted Press Div., Wiley, New York (1975). 13. Jasper, J. J., J. Phys. Chem. Ref. Data (1972) 1, 841. 14. Haile, J. M., Gubbins, Κ. E. and Gray, C. G., J. Chem. Phys. (1976) 64, 1852. RECEIVED
August 15, 1978.
Lykos; Computer Modeling of Matter ACS Symposium Series; American Chemical Society: Washington, DC, 1978.
