2
S i m u l a t i n g L i q u i d Water near M i n e r a l Surfaces: Current Methods and Limitations David J. Mulla Department of Agronomy and Soils, Washington State University, Pullman, WA 99164-6420 It is important to propose molecular and theoretical models to describe the forces, energy, structure and dynamics of water near mineral surfaces. Our understanding of experimental results concerning hydration forces, the hydrophobic effect, swelling, reaction kinetics and adsorption mechanisms in aqueous colloidal systems is rapidly advancing as a result of recent Monte Carlo (MC) and molecular dynamics (MD) models for water properties near model surfaces. This paper reviews the basic MC and MD simulation techniques, compares and contrasts the merits and limitations of various models for water-water interactions and surfacewater interactions, and proposes an interaction potential model which would be useful in simulating water near hydrophilic surfaces. In addition, results from selected MC and MD simulations of water near hydrophobic surfaces are discussed in relation to experimental results, to theories of the double layer, and to structural forces in interfacial systems. Recent evidence (1) suggests that reactions at the mineral/liquid interface were involved in the beginnings of life on Earth. Not surprisingly, the nature and properties of mineral/water interfaces are of interest to physicists, chemists, physical chemists, applied mathematicians, colloid scientists, geochemists, soil scientists and civil engineers. Of particular interest is an increased understanding of the role of water in colloidal swelling, solute hydration, reaction kinetics, adsorption mechanisms, and ion exchange. The theoretical study (2,3) of this interface is made inherently difficult by virtue of the complex, many-body nature of the interaction potentials and forces involving surfaces, counterions, and water. Hence, many models of the interfacial region explicitly specify the forces between colloidal particles or between solutes, but few account for the many-body interaction forces of the solvent. 0097-6156/ 86/ 0323-0020S06.00/ 0 © 1986 American Chemical Society
2.
MULLA
Simulating
Liquid
Water near Mineral
Surfaces
21
E x p e r i m e n t a l s t u d i e s o f the thermodynamic, s p e c t r o s c o p i c and t r a n s p o r t p r o p e r t i e s o f m i n e r a l / w a t e r i n t e r f a c e s have been e x t e n s i v e , a l b e i t c o n f l i c t i n g a t times ( 4 - 1 0 ) . Ambiguous terms such as " h y d r a t i o n f o r c e s " , " h y d r o p h o b i c i n t e r a c t i o n s " , and " s t r u c t u r e d w a t e r " have a r i s e n to d e s c r i b e i n t e r f a c i a l p r o p e r t i e s which have been d i f f i c u l t to q u a n t i f y and e x p l a i n . A d e t a i l e d s t a t i s t i c a l - m e c h a n i c a l d e s c r i p t i o n of the f o r c e s , e n e r g i e s and p r o p e r t i e s of water a t m i n e r a l surfaces i s c l e a r l y d e s i r a b l e . M o l e c u l a r p r e d i c t i o n s of the p r o p e r t i e s o f i n t e r f a c i a l systems are now becoming p o s s i b l e as a r e s u l t o f r a p i d advances i n l i q u i d s t a t e chemical p h y s i c s and computer t e c h n o l o g y . The o b j e c t i v e s of t h i s paper are 1) to r e v i e w the general approaches and models used i n Monte C a r l o (MC) and m o l e c u l a r dynamics (MD) s i m u l a t i o n s of i n t e r f a c i a l s y s t e m s , 2) to d e s c r i b e and d i s c u s s r e s u l t s from s e l e c t e d s i m u l a t i o n s t u d i e s of i n t e r f a c i a l w a t e r , and 3) to d i s c u s s the major l i m i t a t i o n s of these t e c h n i q u e s and to o f f e r s u g g e s t i o n s f o r overcoming them. General S i m u l a t i o n Approaches In most MC (11,12) and MD (12,13) p a r t i c l e s are p l a c e d i n a c e l l of i n t e r a c t i o n p o t e n t i a l energy (U ) p o t e n t i a l s (U-jj) between p a r t i c l e s N
u
N
=
S
I
s t u d i e s , a s m a l l number (N) of f i x e d volume (V) and the t o t a l from a l l p a i r w i s e i n t e r a c t i o n i and j i s c a l c u l a t e d :
Ï J
Uij
(l)
P a r t i c l e i n t e r a c t i o n s are not computed beyond a c u t o f f r a d i u s o f from f o u r t o e i g h t Angstroms t o improve c o m p u t a t i o n a l e f f i c i e n c y by n e g l e c t i n g l o n g - r a n g e i n t e r a c t i o n s which c o n t r i b u t e l i t t l e t o the o v e r a l l s t r u c t u r e of the f l u i d . P e r i o d i c boundary c o n d i t i o n s are imposed by f i l l i n g the space around the b a s i c c e l l w i t h image c e l l s t r a n s l a t e d by m u l t i p l e s o f the u n i t l e n g t h o f the b a s i c c e l l . Thus, p a r t i c l e s near the c e l l b o u n d a r i e s of the b a s i c c e l l i n t e r a c t w i t h image p a r t i c l e s , not w i t h empty c a v i t i e s . This technique prevents s p u r i o u s edge e f f e c t s from a f f e c t i n g the r e s u l t s o f the s i m u l a t i o n . P e r i o d i c c e l l b o u n d a r i e s a l s o a l l o w a s m a l l sample o f p a r t i c l e s to e x h i b i t p r o p e r t i e s c h a r a c t e r i s t i c of a much l a r g e r sample s i z e . A l t h o u g h the method o f moving p a r t i c l e s t o new l o c a t i o n s and o f o b t a i n i n g e q u i l i b r i u m c o n f i g u r a t i o n s f o r the p a r t i c l e s d i f f e r s f o r the MC and MD methods, i n both t e c h n i q u e s the p o s i t i o n s and o r i e n t a t i o n s of thousands o f c o n f i g u r a t i o n s are generated and used t o c a l c u l a t e average p r o p e r t i e s o f the s y s t e m . In the MC method, ensemble average p r o p e r t i e s can be computed. These may i n c l u d e s t r u c t u r a l and thermodynamic p r o p e r t i e s such as d e n s i t y , d i p o l e d i r e c t i o n c o s i n e , hydrogen bond e n e r g i e s and number, r a d i a l d i s t r i b u t i o n f u n c t i o n s , i n t e r n a l e n e r g y , heat c a p a c i t y and i n t e r n a l p r e s s u r e . In the MD method, time average p r o p e r t i e s which are e i t h e r s t r u c t u r a l or dynamic can be computed. For i n s t a n c e , these p r o p e r t i e s may i n c l u d e d e n s i t y , d i p o l e d i r e c t i o n c o s i n e , hydrogen bond e n e r g i e s and number, r a d i a l d i s t r i b u t i o n f u n c t i o n s , d i p o l e r e l a x a t i o n t i m e , and s e l f - d i f f u s i o n coefficients. Thus, the key d i f f e r e n c e between the two t e c h n i q u e s i s t h a t the MC method g e n e r a l l y a l l o w s o n l y s t a t i c p r o p e r t i e s t o be
22
G E O C H E M I C A L PROCESSES AT M I N E R A L SURFACES
e v a l u a t e d , w h i l e t h e MD method g e n e r a l l y dependent phenomena to be s t u d i e d .
a l l o w s both s t a t i c and time
Monte C a r l o Methods. A l t h o u g h s e v e r a l s t a t i s t i c a l mechanical ensembles may be s t u d i e d u s i n g MC methods ( 2 , 1 2 , 1 4 ) , t h e c a n o n i c a l ensemble has been the most f r e q u e n t l y used ensemble f o r s t u d i e s o f i n t e r f a c i a l systems. In the c a n o n i c a l ensemble, the number o f m o l e c u l e s ( N ) , c e l l volume (V) and temperature (T) a r e f i x e d . Hence, the c a n o n i c a l ensemble i s denoted by t h e symbols NVT. The c h o i c e o f ensemble determines which thermodynamic p r o p e r t i e s can be computed. In the NVT ensemble one cannot compute the chemical p o t e n t i a l or e n t r o p y o f t h e system; two p r o p e r t i e s which are o f c r i t i c a l importance f o r i n t e r f a c i a l systems. The c h o i c e o f an ensemble a l s o determines the s a m p l i n g a l g o r i t h m used t o generate m o l e c u l a r c o n f i g u r a t i o n s from random moves o f the m o l e c u l e s . In MC methods the u l t i m a t e o b j e c t i v e i s t o e v a l u a t e m a c r o s c o p i c p r o p e r t i e s from i n f o r m a t i o n about m o l e c u l a r p o s i t i o n s generated over phase s p a c e . To e v a l u a t e average m a c r o s c o p i c p r o p e r t i e s , p, i n t h e c a n o n i c a l ensemble from s t a t i s t i c a l m e c h a n i c s , the f o l l o w i n g expression i s used: ρ = Σ
ρ
[ exp(-U (q)/kT)]/Q(N,V,T) Pq
N
(2)
where p i s t h e v a l u e o f the m a c r o s c o p i c p r o p e r t y ρ i n the q t h c o n f i g u r a t i o n , k i s B o l t z m a n n ' s c o n s t a n t , Q(N,V,T) i s the c a n o n i c a l ensemble p a r t i t i o n f u n c t i o n , and t h e summation runs over a l l q e q u i l i b r i u m molecular configurations. Thus, E q u a t i o n 2 suggests t h a t to e v a l u a t e p r o p e r t i e s i n t h e c a n o n i c a l ensemble MC method, t h e p r o b a b i l i t y w i t h which any m o l e c u l a r c o n f i g u r a t i o n o c c u r s s h o u l d be proportional to exp(-U (q)/kT). The s p e c i f i c a l g o r i t h m f o r g e n e r a t i n g new c o n f i g u r a t i o n s t h a t s a t i s f y t h i s requirement i n v o l v e s i ) s e l e c t i n g a m o l e c u l e a t random, i i ) s e l e c t i n g C a r t e s i a n c e n t e r - o f - m a s s d i s p l a c e m e n t c o o r d i n a t e s randomly over an i n t e r v a l which i s not g r e a t e r than h a l f t h e c e l l l e n g t h , i i i ) s e l e c t i n g a r o t a t i o n angle a t random, and i v ) c a l c u l a t i n g t h e new energy, U ( q + l ) , o f t h e new c o n f i g u r a t i o n generated by moves i ) - i i i ) . The f i n a l s t e p i n v o l v e s d e c i d i n g whether t o a c c e p t or r e j e c t the new, but random, configuration. T h i s i s done by g e n e r a t i n g a random number between z e r o and u n i t y and comparing t h e random number t o t h e q u a n t i t y : q
N
N
exp[-(U (q+l) N
- U (q))/kT] N
(3)
I f t h e random number i s l e s s than or equal t o t h e q u a n t i t y i n E q u a t i o n 3, then the new move i s a c c e p t e d ; o t h e r w i s e the move i s r e j e c t e d . T h i s acceptance c r i t e r i a i s o f t e n made even more s t r i n g e n t by r e q u i r i n g t h a t as few as 50% o f the moves s a t i s f y i n g E q u a t i o n 3 a r e actually selected. Note t h a t these procedures always f a v o r moves which l e a d t o reduced t o t a l i n t e r a c t i o n e n e r g i e s . M o l e c u l a r Dynamics Methods. In c o n t r a s t t o the MC method, both k i n e t i c and s t r u c t u r a l p r o p e r t i e s o f a m o l e c u l a r system can be e v a l u a t e d from MD s t u d i e s . These p r o p e r t i e s are e v a l u a t e d as averages over c o n f i g u r a t i o n s generated d u r i n g t i m e . In m i c r o c a n o n i c a l ensemble s t u d i e s w i t h t h e MD method, the p r o p e r t i e s which a r e c o n t r o l l e d
2.
MULLA
Simulating
Liquid
Water near Mineral
23
Surfaces
i n c l u d e Ν, V and t o t a l energy, E. T o t a l energy i s computed from the sum of t o t a l k i n e t i c energy and t o t a l p o t e n t i a l energy. The p o t e n t i a l energy i s e v a l u a t e d from an e x p r e s s i o n i n v o l v i n g a summation over the i n t e r a c t i o n p o t e n t i a l s between i n d i v i d u a l p a r t i c l e s i and j , U-jj, as g i v e n i n E q u a t i o n 1. Temperature i s not f i x e d i n the m i c r o c a n o n i c a l MD method, s i n c e i t v a r i e s w i t h f l u c t u a t i o n s i n t o t a l k i n e t i c energy. To determine the movement of m o l e c u l e s , the f o l l o w i n g a l g o r i t h m (15) i s o f t e n u s e d . The f o r c e a c t i n g on the i t h atom i n a m o l e c u l e (F^) i determined from the s p a t i a l d e r i v a t i v e of the t o t a l i n t e r a c t i o n p o t e n t i a l energy of t h a t p a r t i c l e : s
Pi = -V Ej U i j
(4)
A c e n t e r e d f i n i t e d i f f e r e n c e scheme i s used to c a l c u l a t e p o s i t i o n of the i t h atom, x-j, a s h o r t time At i n the f u t u r e : Xj(t
+ At)
= -x^t
- At)
+ 2x-j(t) + ( A t / ^ ) ^ 2
the
(5)
where a l l t h a t i s r e q u i r e d to determine t h i s new p o s i t i o n are the p r e s e n t and p a s t l o c a t i o n s , and the f o r c e from E q u a t i o n 4. Typically, to ensure n u m e r i c a l s t a b i l i t y of the a l g o r i t h m , the magnitude of At i s on the o r d e r of 1 0 " seconds or l e s s . The v e l o c i t y of the atom, v-j, i s determined by the e x p r e s s i o n : 1 5
V i
(t)
= [l/UAtOHx^t
+ At)
- Xi(t
- At)]
(6)
Hence, the MD method i n v o l v e s n u m e r i c a l i n t e g r a t i o n of the e q u a t i o n s of motion f o r a l l p a r t i c l e s each time a new c o n f i g u r a t i o n i s g e n e r a t e d , w h i l e the MC method o n l y i n v o l v e s movement of one random p a r t i c l e f o r each new c o n f i g u r a t i o n . I t s h o u l d be noted t h a t many o t h e r n u m e r i c a l a l g o r i t h m s are a v a i l a b l e f o r the MD method ( 1 3 ) , and t h a t maximum a c c u r a c y r e s u l t s from the use of a l g o r i t h m s t h a t i n c l u d e h i g h e r powers of At than are g i v e n i n E q u a t i o n 5. Interaction
Potentials
E q u a t i o n s 3-4 show t h a t the form o f the i n t e r a c t i o n p o t e n t i a l s used i n s i m u l a t i n g i n t e r f a c i a l water i s c r i t i c a l . Of i n t e r e s t f o r i n t e r f a c i a l systems are both the i n t e r a c t i o n p o t e n t i a l between water m o l e c u l e s and t h a t between the s u r f a c e and a water m o l e c u l e . The f i r s t MC (16) and MD (Γ7) s t u d i e s were used to s i m u l a t e the p r o p e r t i e s of s i n g l e p a r t i c l e f l u i d s . A l t h o u g h the b a s i c MC (11,12) and MD (12,13) methods have changed l i t t l e s i n c e the e a r l i e s t s i m u l a t i o n s , the systems s i m u l a t e d have c o n t i n u a l l y i n c r e a s e d i n complexity. The a b i l i t y to s i m u l a t e complex i n t e r f a c i a l systems has r e s u l t e d p a r t l y from improvements i n s i m u l a t i o n a l g o r i t h m s (15,18) or i n the i n t e r a c t i o n p o t e n t i a l s used to model s o l i d s u r f a c e s ( 1 9 ) . The major r e a s o n , however, f o r t h i s a b i l i t y has r e s u l t e d from the i n c r e a s i n g s o p h i s t i c a t i o n of the i n t e r a c t i o n p o t e n t i a l s used to model liquid-liquid interactions. These advances have i n v o l v e d the use of the f o l l o w i n g p o t e n t i a l s : Lennard-Jones 12-6 ( 2 0 ) , Rowlinson ( 2 1 ) , BNS
24
G E O C H E M I C A L P R O C E S S E S AT M I N E R A L S U R F A C E S
(22), (28).
ST2 ( 2 3 ) ,
MCY ( 2 4 ) ,
CF ( 2 5 ) ,
PE ( 2 6 ) ,
TIP4P ( 2 7 ) ,
and MCY+CC+DC
Water P o t e n t i a l s . The ST2 ( 2 3 ) , MCY ( 2 4 ) , and CF (25) p o t e n t i a l s are c o m p u t a t i o n a l l y t r a c t a b l e and a c c u r a t e models f o r two-body w a t e r - w a t e r i n t e r a c t i o n p o t e n t i a l s . The ST2, MCY and CF models have f i v e , f o u r , and t h r e e i n t e r a c t i o n s i t e s and have f o u r , t h r e e and t h r e e charge centers, respectively. N e i t h e r the ST2 nor the MCY p o t e n t i a l s a l l o w OH or HH d i s t a n c e s t o v a r y , whereas bond l e n g t h s are f l e x i b l e w i t h the CF model. W h i l e both the ST2 and CF p o t e n t i a l s are e m p i r i c a l models, the MCY p o t e n t i a l i s d e r i v e d from ab i n i t i o c o n f i g u r a t i o n i n t e r a c t i o n m o l e c u l a r o r b i t a l methods (24) u s i n g many g e o m e t r i c a l arrangements of water d i m e r s . The MCY+CC+DC w a t e r - w a t e r p o t e n t i a l (28) i s a r e c e n t m o d i f i c a t i o n of the MCY p o t e n t i a l which a l l o w s f o u r body i n t e r a c t i o n s to be e v a l u a t e d . In comparison to the two-body p o t e n t i a l s d e s c r i b e d above, the MCY+CC+DC p o t e n t i a l r e q u i r e s a supercomputer or a r r a y p r o c e s s o r i n o r d e r t o be c o m p u t a t i o n a l l y f e a s i b l e . T h e r e f o r e , the ST2, MCY and CF p o t e n t i a l s are g e n e r a l l y more economical t o use than the MCY+CC+DC p o t e n t i a l . A comparison of the b u l k water p r o p e r t i e s p r e d i c t e d by the ST2, MCY, and CF models i n s i m u l a t i o n s i s g i v e n i n T a b l e I. These d a t a were o b t a i n e d from ( 2 ) , u n l e s s o t h e r w i s e n o t e d . T a b l e I. Comparison of water p r o p e r t i e s f o r the ST2, MCY and CF s i m u l a t i o n models and b u l k water a t a p p r o x i m a t e l y 298 K. Property
ST2
MCY
CF
b u l k water
-U (kJ/mol) C (J/K/mol) μ (Debye u n i t s ) PV/NkT (V/N=l) D ( 1 0 " m /sec)
34 71 2.35 0.09 3.1 (30)
28.5 79 2.26 8.5 (29) 2.3 (29)
33
34 75 1.86 0.05 2.85
v
9
2
--
1.86 0.1 1.10
(31)
(32)
R e s u l t s i n T a b l e I i l l u s t r a t e some o f the s t r e n g t h s and weaknesses o f the ST2, MCY and CF models. A l l models, except the MCY model, a c c u r a t e l y p r e d i c t the i n t e r n a l e n e r g y , - U . C o n s t a n t volume heat c a p a c i t y , C , i s a c c u r a t e l y p r e d i c t e d by each model f o r which data i s a v a i l a b l e . The ST2 and MCY models o v e r p r e d i c t the d i p o l e moment, μ, w h i l e the CF model p r e d i c t i o n i s i d e n t i c a l w i t h the v a l u e for bulk water. The r a t i o PV/NkT a t a l i q u i d d e n s i t y of u n i t y i s tremendously i n e r r o r f o r the MCY model, w h i l e both the ST2 and CF models p r e d i c t i o n s are r e a s o n a b l e . T h i s l a r g e e r r o r u s i n g the MCY model s u g g e s t s t h a t i t w i l l n o t , i n g e n e r a l , s i m u l a t e thermodynamic p r o p e r t i e s o f water a c c u r a t e l y ( 2 9 ) . V a l u e s o f the s e l f - d i f f u s i o n c o e f f i c i e n t , D, f o r each o f the water models e x c e p t the CF model agree f a i r l y w e l l w i t h the v a l u e f o r b u l k w a t e r . v
In s i m u l a t i n g i n t e r f a c i a l w a t e r , i t i s i m p o r t a n t t o use a model f o r w a t e r - w a t e r i n t e r a c t i o n s which y i e l d s a c c u r a t e r e s u l t s i n s i m u l a t i o n s of bulk water. Each of the models d i s c u s s e d here have o b v i o u s advantages and d i s a d v a n t a g e s . The CF model i s g e n e r a l l y more
2.
MULLA
Simulating
Liquid
Water near Mineral
25
Surfaces
a c c u r a t e i n p r e d i c t i n g b u l k water p r o p e r t i e s than the o t h e r models i n Table I. Two drawbacks o f the ST2 model are i t s r i g i d n e s s and o v e r l y t e t r a h e d r a l geometry. The MCY p o t e n t i a l may l e a d to s p u r i o u s r e s u l t s for i n t e r f a c i a l water, since i t generates excessive i n t e r n a l pressures. Surface P o t e n t i a l s . C o n s i d e r the form of the s u r f a c e - w a t e r i n t e r a c t i o n p o t e n t i a l f o r an i n t e r f a c i a l system w i t h a hydrophobic surface. The oxygen atom of any water m o l e c u l e i s a c t e d upon by an e x p l i c i t l y uncharged s u r f a c e d i r e c t l y below i t v i a the Lennard-Jones p o t e n t i a l (U|_j): U (Rij) = A[(a/R L J
i : j
)
a
- (a/R
) ]
(7)
b
i : j
where R^j i s the d i s t a n c e between the j t h s u r f a c e atom and the oxygen atom on the i t h water m o l e c u l e , and A and σ are parameters which s p e c i f y the depth of the p o t e n t i a l energy w e l l and the d i s t a n c e at which i t s v a l u e f i r s t e q u a l s z e r o , r e s p e c t i v e l y . The exponents a and b s p e c i f y the power law f o r the r e p u l s i v e and a t t r a c t i v e components of the Lennard-Jones p o t e n t i a l , r e s p e c t i v e l y . Three commonly used p a i r s of v a l u e s f o r a and b are 12 and 6, 9 and 3 , or 4 and 2, which produce the Lennard-Jones 12-6 ( 3 3 , 3 4 ) , 9-3 ( 3 5 ) , and 4-2 (36) p o t e n t i a l s , respectively. T y p i c a l v a l u e s f o r parameters of the l a t t e r L e n n a r d Jones p o t e n t i a l s are r e p o r t e d i n Table I I . In g e n e r a l , the depth of the p o t e n t i a l w e l l f o r these p o t e n t i a l s (about - 0 . 5 kcal/mole) i s t y p i c a l o f the energy between hydrophobic s u r f a c e s and p h y s i s o r b e d noble gases. Table I I . T y p i c a l v a l u e s f o r parameters of the Lennard-Jones 1 2 - 6 , 9 - 3 , and 4-2 p o t e n t i a l s .
Lennard-Jones potential 12-6 9-3 4-2
A (kcal/mole ) 0.303 1.202 1.728
a
b
σ (Angstroms)
12 9 4
6 3 2
3.1 2.5 2.0
The above forms f o r the Lennard-Jones s u r f a c e -water i n t e r a c t i o n p o t e n t i a l have been used as models of hydrophobic s u r f a c e s such as p y r o p h y l l i t e , g r a p h i t e , or p a r a f f i n . I f the i n t e n t i o n o f the s t u d y , however, i s to understand i n t e r f a c i a l p r o c e s s e s a t m i n e r a l s u r f a c e s r e p r e s e n t a t i v e of s m e c t i t e s or m i c a , e x p l i c i t e l e c t r o s t a t i c i n t e r a c t i o n s betweeen water m o l e c u l e s and l o c a l i z e d charges a t the s u r f a c e become i m p o r t a n t . Two methods f o r i n c l u d i n g e x p l i c i t e l e c t r o s t a t i c i n t e r a c t i o n s are proposed. In the f i r s t , and more d i f f i c u l t a p p r o a c h , one would need to conduct e x t e n s i v e quantum mechanical c a l c u l a t i o n s o f the p o t e n t i a l energy v a r i a t i o n between a model s u r f a c e and one a d j a c e n t water molecule u s i n g thousands o f d i f f e r e n t g e o m e t r i c a l o r i e n t a t i o n s . This approach has been used i n a l i m i t e d f a s h i o n to s t u d y the i n t e r a c t i o n p o t e n t i a l between water and s u r f a c e Si-OH groups on a l u m i n o s i 1 i c a t e s , s i l i c a t e s and z e o l i t e s ( 3 7 - 3 9 ) .
26
G E O C H E M I C A L PROCESSES AT M I N E R A L SURFACES
A s i m p l e r approach i s t o use an e m p i r i c a l model c o n s i s t i n g o f a Lennard-Jones p o t e n t i a l p l u s a Coulombic term t o e x p l i c i t l y a c c o u n t f o r charges on t h e s u r f a c e . To c a l c u l a t e the magnitude o f t h i s Coulombic c h a r g e , c o n s i d e r t h e p h y s i c a l p r o p e r t i e s o f s m e c t i t e and mica a l u m i n o s i 1 i c a t e s . I t i s known t h a t t h e u n i t c e l l o f s m e c t i t e and mica s u r f a c e s has dimensions o f about 46 square Angstroms ( 3 4 ) . U s i n g t h i s v a l u e f o r t h e a r e a o f a u n i t c e l l and the d a t a i n Table I I I , t h e charge per u n i t c e l l and d e l o c a l i z e d charge on s u r f a c e oxygens o f t y p i c a l s m e c t i t e and mica m i n e r a l s can be computed. Table
Physical
I I I . C a l c u l a t i o n o f t h e number o f charges per u n i t c e l l on t y p i c a l s m e c t i t e and mica s u r f a c e s .
Property
Smectite
s u r f a c e a r e a (m /g) charge d e n s i t y ( e . s . u . / m ) c a t i o n exchange c a p a c i t y (meq/g) charge d e n s i t y (number/unit c e l l ) d e l o c a l i z e d charge per s u r f a c e oxygen 2
2
750 4 χ 10 1.038 0.38 0.06
Mica 100 10 χ 1 0 0.346 0.96 0.16
8
8
To s i m u l a t e s m e c t i t e o r mica m i n e r a l s , a t o t a l o f about 0.4 and 1 e x p l i c i t n e g a t i v e c h a r g e s , r e s p e c t i v e l y , need t o be a s s i g n e d t o each u n i t c e l l on t h e s u r f a c e . T h i s charge s h o u l d be d e l o c a l i z e d over about s i x oxygen atoms s u r r o u n d i n g t h e d i t r i g o n a l c a v i t i e s o f t h e s m e c t i t e and mica s u r f a c e s , s i n c e these charges o r i g i n a t e from o c t a h e d r a l o r t e t r a h e d r a l s i t e s w i t h i n t h e c r y s t a l and n o t from t h e s u r f a c e atoms. A proposed form f o r t h e w a t e r - s u r f a c e i n t e r a c t i o n p o t e n t i a l , U ^ , s u i t a b l e f o r s i m u l a t i o n s o f s m e c t i t e o r mica s u r f a c e s i n t e r a c t i n g w i t h t h e ST2 model o f water i s : 4 U
WS( ijAj) R
= A[(a/R
i ; j
)
a
-
(o/R
i : j
) ] + S(RT j ) b
Σ
(q qj) / d a
a j
(8)
a=l where R ^ j , A , a , b and σ a r e as d e f i n e d i n E q u a t i o n 7 and T a b l e I I , d j i s t n e d i s t a n c e between t h e j t h s u r f a c e atom and t h e a t h charge ( h a v i n g charge q ) on t h e water m o l e c u l e , qj i s t h e d e l o c a l i z e d charge on t h e j t h s u r f a c e atom from T a b l e I I I , and S(R^j ) i s t h e s w i t c h i n g f u n c t i o n o f t h e ST2 water p o t e n t i a l ( 2 3 ) . The magnitude o f charge on each o f t h e f o u r p o i n t charges f o r ST2 water i n E q u a t i o n 8 i s 0.2357. a
a
A p l o t o f t h e Lennard-Jones 9-3 form o f E q u a t i o n s 7 and 8 f o r ST2 water i n t e r a c t i n g w i t h s m e c t i t e and mica s u r f a c e s i s shown i n F i g u r e 1. V a l u e s f o r the parameters used i n F i g u r e 1 a r e g i v e n i n T a b l e s II and I I I , and i n r e f e r e n c e ( 2 3 ) . The water m o l e c u l e i s o r i e n t e d so t h a t i t s p r o t o n s f a c e t h e s u r f a c e and i t s lone p a i r e l e c t r o n s f a c e away from t h e s u r f a c e , and t h e p r o t o n s a r e e q u i d i s t a n t from t h e s u r f a c e . Note t h a t t h e depth o f t h e p o t e n t i a l w e l l i n F i g u r e 1 f o r i n t e r a c t i o n s w i t h t h e s m e c t i t e s u r f a c e and mica s u r f a c e a r e
MULLA
Simulating
-4
Liquid
ι 0
.
Water near Mineral
ι
2
.
ι
4
.
ι
6
Surfaces
.
ι
8
,
I
10
DISTRNCE (Angstroms) F i g u r e 1. Comparison o f ST2 w a t e r - s u r f a c e i n t e r a c t i o n s computed from E q u a t i o n s 7 o r 8 u s i n g parameters f o r the Lennard-Jones 9-3 p o t e n t i a l i n Table II and the d e l o c a l i z e d charge magnitude f o r s m e c t i t e and mica s u r f a c e s i n Table III.
28
G E O C H E M I C A L PROCESSES AT M I N E R A L SURFACES
about - 1 . 5 and - 3 . 5 k c a l / m o l e , r e s p e c t i v e l y . These i n t e r a c t i o n e n e r g i e s a r e s i m i l a r i n magnitude t o weak hydrogen bond e n e r g i e s . An i n t e r a c t i o n p o t e n t i a l between t h e s u r f a c e and i o n s may a l s o be needed i n s i m u l a t i n g c o u n t e r i o n d i f f u s i o n f o r t h e s m e c t i t e and mica s u r f a c e models. The form o f such an i n t e r a c t i o n p o t e n t i a l remains t o be d e t e r m i n e d . T h i s may n o t pose a s i g n i f i c a n t p r o b l e m , s i n c e r e c e n t e v i d e n c e (40) suggests t h a t over 98% o f t h e c a t i o n s near s m e c t i t e s u r f a c e s l i e w i t h i n t h e shear p l a n e . For s p e c i f i c a l l y adsorbed c a t i o n s such as p o t a s s i u m or c a l c i u m , t h e s u r f a c e - i o n i n t e r a c t i o n s can a l s o be n e g l e c t e d i f i t i s assumed t h a t c a t i o n d i f f u s i o n c o n t r i b u t e s l i t t l e t o t h e water s t r u c t u r e . In s i m u l a t i n g t h e i n t e r a c t i o n p o t e n t i a l between c o u n t e r i o n s and i n t e r f a c i a l w a t e r , a w a t e r - i o n i n t e r a c t i o n p o t e n t i a l s i m i l a r t o those a l r e a d y developed f o r MD s i m u l a t i o n s (41-43) c o u l d be s p e c i f i e d . Simulations of I n t e r f a c i a l
Water
S e v e r a l MC and MD s t u d i e s o f i n t e r f a c i a l water near hydrophobic s u r f a c e s have been r e p o r t e d ( 3 3 - 3 6 , 4 4 - 4 8 ) . Both o f t h e MC s t u d i e s ( 3 5 , 4 5 ) , as w e l l as the f o u r MD s t u d i e s ( 3 3 , 3 4 , 3 6 , 4 7 ) r e p o r t i n g d e t a i l e d o b s e r v a t i o n s o f i n t e r f a c i a l water are d i s c u s s e d h e r e . This comparison w i l l show t h a t c h o i c e o f t h e w a t e r - w a t e r p o t e n t i a l i s c r i t i c a l f o r such s t u d i e s . I t w i l l a l s o i l l u s t r a t e t h e wide range o f i n t e r f a c i a l p r o p e r t i e s which can be s t u d i e d u s i n g computer simulations. R e s u l t s from t h e e a r l y p i o n e e r i n g MC s t u d i e s f o r i n t e r f a c i a l water a r e summarized i n T a b l e I V . Table
IV.
R e s u l t s from Monte C a r l o S i m u l a t i o n s o f Water Hydrophobic S u r f a c e s .
# molecules c e l l dimensions (nm^) c e l l d e n s i t y (g/cc) temperature (K) # configurations water p o t e n t i a l surface potential range o f d e n s i t y o s c i l l a t i o n s (g/cc) d e n s i t y t r e n d towards s u r f a c e s hydrogen bonding t r e n d towards s u r f a c e s i n t e r n a l energy t r e n d towards s u r f a c e s preferred dipole orientation near s u r f a c e s
Near
Reference 35
Reference 45
216 7.127 0.906 298 2.5 χ 1 0 Row!inson L - J 9-3 1.1 t o 1.5 decreases decreases decreases
150 4.5 0.997 300 1 χ 10 MCY hard w a l l 0.5 t o 2 . 3 increases increases increases
yes
yes
6
4
These r e s u l t s i n d i c a t e t h a t , compared t o b u l k w a t e r , i n t e r f a c i a l water e x h i b i t s unique o s c i l l a t i o n s i n d e n s i t y w i t h d i s t a n c e from t h e s u r f a c e and p r e f e r e n t i a l d i p o l a r o r i e n t a t i o n s . Both s i m u l a t i o n s r e p o r t d e n s i t y v a l u e s which a r e u n r e a s o n a b l e . P a r t o f t h i s problem a r i s e s from a t t e m p t i n g t o f i x t h e water d e n s i t y based on t h e average c e l l volume and t h e number o f water m o l e c u l e s ; an approach which
2.
MULLA
Simulating
Liquid
Water near Mineral
29
Surfaces
o v e r l o o k s the f a c t t h a t the c e l l volume near the s u r f a c e s i s g e n e r a l l y f r e e of water m o l e c u l e s due to r e p u l s i v e f o r c e s a r i s i n g from the surfaces. For the study u s i n g the MCY p o t e n t i a l ( 4 5 ) , a more s e r i o u s problem i n v o l v e s the e x c e s s i v e i n t e r n a l p r e s s u r e s generated by the MCY p o t e n t i a l , which l e a d to e x c e s s i v e d e n s i t y o s c i l l a t i o n s and i n c r e a s i n g d e n s i t y near the s u r f a c e s . C o n s i s t e n c y between the two s i m u l a t i o n s i s poor; whereas the f i r s t (35) p r e d i c t s decreased d e n s i t y , hydrogen bonding and i n t e r n a l energy near the s u r f a c e s , the second (45) r e p o r t s e x a c t l y the o p p o s i t e t r e n d s . These d i f f e r e n c e s are a l l p r o b a b l y r e l a t e d to the use of d i f f e r e n t w a t e r - w a t e r p o t e n t i a l s . R e s u l t s o f s e l e c t e d MD s t u d i e s o f i n Table V. Table V.
i n t e r f a c i a l water are
R e s u l t s from S e l e c t e d M o l e c u l a r Dynamics S t u d i e s o f Near Hydrophobic S u r f a c e s .
# molecules c e l l volume (ntrr*) c e l l d e n s i t y (g/cc) temperature (K) t r a j e c t o r y time (ps) water p o t e n t i a l surface p o t e n t i a l range of d e n s i t y o s c i l l a t i o n s (g/cc) density trend towards s u r f a c e s hydrogen bonding t r e n d towards s u r f a c e s preferred dipolar orientations self-diffusion coeff. near s u r f a c e s (m /s) near midplane (m /s) d i p o l e r e l a x a t i o n time near s u r f a c e s (ps) near midplane (ps) 2
2
Water
(36)
Reference (34)
216 8.787 0.74 287 20 ST2 L - J 12-6
150 5.198 0.87 304 14 MCY L - J 4-2
256 9.156 0.84 286 0.75 ST2 L - J 12-6
0.9 t o
0.5 to
Reference (47)
Reference
150 4.5 1.0 301 25 ST2 hard w a l l 0.9 to
Reference
(33)
1.0
reported
1.0
3.2
0.8 to
1.1
decreases
decreases
increases
decreases
—
—
—
decreases
yes
yes
yes
yes
3.3 χ 1 0 " 4.2 χ 1 0 " 3.1 2.1
9 9
4.8 χ 1 0 " 3.3 χ 1 0 " — —
9 9
3.1 χ 1 0 " 3.7 χ 1 0 " 2.3 2.0
9 9
2.1 χ 1 0 " 2.7 χ 1 0 " — —
The r e s u l t s i n Table V i l l u s t r a t e t h a t MD s t u d i e s , compared t o the MC r e s u l t s i n Table IV, f a c i l i t a t e the i n v e s t i g a t i o n o f t r a n s p o r t and time-dependent p r o p e r t i e s . A l s o , they show t h a t use of the MCY p o t e n t i a l l e a d s to v e r y l a r g e d e n s i t y o s c i l l a t i o n s and i n c r e a s i n g water d e n s i t y near the s u r f a c e s . T h i s appears t o be a s e r i o u s drawback t o the use of the MCY p o t e n t i a l i n s i m u l a t i o n s of i n t e r f a c i a l water. R e s u l t s from the i n v e s t i g a t i o n s u s i n g the ST2 p o t e n t i a l show t h a t i n t e r f a c i a l water d e n s i t y i s a p p r o x i m a t e l y 1.0 g/cc, w i t h a tendency f o r decreased d e n s i t y and hydrogen bonding near the s u r f a c e s . As i n the MC s i m u l a t i o n s , o r i e n t a t i o n s o f the water d i p o l e moment are a f f e c t e d by the presence of a s o l i d / l i q u i d i n t e r f a c e , and an
9 9
30
G E O C H E M I C A L P R O C E S S E S AT M I N E R A L S U R F A C E S
a p p r e c i a b l e decrease i n d i p o l e r e l a x a t i o n and the water c o e f f i c i e n t are u s u a l l y observed near the s u r f a c e s .
self-diffusion
Comparison w i t h E x p e r i m e n t . How do the r e s u l t s of these s i m u l a t i o n s compare w i t h e x p e r i m e n t a l r e s u l t s on c o r r e s p o n d i n g s y s t e m s , and what do they i n f e r about i n t e r p r e t i n g the " h y d r o p h o b i c e f f e c t " , " h y d r a t i o n f o r c e s " , and the " s t r u c t u r e " of i n t e r f a c i a l water? S t r u c t u r a l l y , the i n t e r f a c i a l water e x h i b i t s c l e a r d e n s i t y o s c i l l a t i o n s which extend t o a t l e a s t 15 Angstroms from the s u r f a c e s ( 3 4 ) . Since t h i s s i m u l a t i o n i n v o l v e d a h y d r o p h o b i c , n e u t r a l s u r f a c e , these e f f e c t s are d i r e c t l y a t t r i b u t a b l e t o the presence of the s u r f a c e s , and not to an e f f e c t o f charged c o u n t e r i o n s . F u r t h e r m o r e , s i n c e no l o n g - r a n g e changes i n hydrogen bonding p a t t e r n s were observed due to t h i s s t r u c t u r a l r e o r d e r i n g ( 3 4 ) , the MD r e s u l t s suggest t h a t the hydrophobic e f f e c t i s due t o e n t r o p y changes i n the i n t e r f a c i a l l i q u i d r a t h e r than to l o n g range bonding e f f e c t s . The presence of d e n s i t y o s c i l l a t i o n s i n s i m u l a t i o n r e s u l t s of i n t e r f a c i a l water has s t i m u l a t e d e x p e r i m e n t a l s t u d i e s which were e x p l i c i t l y designed to d e t e c t t h e i r p r e s e n c e . Structural forces a s s o c i a t e d w i t h d e n s i t y o s c i l l a t i o n s i n water next t o mica s u r f a c e s have r e c e n t l y been measured e x p e r i m e n t a l l y ( 4 9 ) . A l t h o u g h , mica i s not a hydrophobic s u r f a c e , i t s h o u l d be p o i n t e d out t h a t t h e r e i s t h e o r e t i c a l b a s i s f o r s u g g e s t i n g t h a t m o l e c u l a r l a y e r i n g near s u r f a c e s and the accompanying o s c i l l a t o r y f o r c e s are r e s p o n s i b l e f o r both the " h y d r a t i o n " and " h y d r o p h o b i c " e f f e c t s ( 2 ) . Whether t h i s f o r c e i s a t t r a c t i v e (hydrophobic e f f e c t ) or r e p u l s i v e ( h y d r a t i o n e f f e c t ) depends upon how the d e n s i t y o s c i l l a t i o n s f i t i n t o the r e g i o n between the s u r f a c e s . Much MC work i s needed w i t h both hydrophobic and h y d r o p h i l i c s u r f a c e s u s i n g the grand c a n o n i c a l ensemble to determine the chemical p o t e n t i a l and e n t r o p y of the i n t e r f a c i a l water a t v a r i o u s s u r f a c e s e p a r a t i o n s i n o r d e r t o b e t t e r understand the magnitude of these e f f e c t s . T h e o r e t i c a l e x p l a n a t i o n s f o r the " h y d r a t i o n f o r c e " (50) and the " h y d r o p h o b i c e f f e c t " (51) o f t e n i n v o l v e an a n a l y s i s of the f o r c e s emanating from o r i e n t e d m o l e c u l a r d i p o l e or quadrupole moments. All o f the computer s i m u l a t i o n s f o r i n t e r f a c i a l water d i s c u s s e d above found s i g n i f i c a n t e v i d e n c e f o r p r e f e r r e d d i p o l a r o r i e n t a t i o n s near the surfaces. In most o f the s t u d i e s , a l l o f which used n o n - p o l a r i z a b l e , p a i r w i s e i n t e r a c t i n g water models, the tendency f o r p r e f e r r e d d i p o l a r o r i e n t a t i o n s diminishes in a continuous fashion with i n c r e a s i n g d i s t a n c e from the s u r f a c e s and was n e g l i g i b l e a t a d i s t a n c e of from ten to f i f t e e n Angstroms from the s u r f a c e s . When r e a l i s t i c water p o t e n t i a l s i n c o r p o r a t i n g c o o p e r a t i v e e f f e c t s become a v a i l a b l e , t h i s e f f e c t can be expected to become even more s i g n i f i c a n t . A c c o r d i n g t o the Kirkwood t h e o r y of p o l a r d i e l e c t r i c s , s i m p l e r e l a t i o n s (23) between m o l e c u l a r d i p o l e moment v e c t o r s and the meansquare t o t a l d i p o l e moment o f water c l u s t e r s can be used to compute the s t a t i c d i e l e c t r i c c o n s t a n t of w a t e r . As the n o r m a l i z e d meansquare t o t a l d i p o l e moment i n c r e a s e s towards u n i t y , t h e o r y p r e d i c t s d e c r e a s e s i n the s t a t i c d i e l e c t r i c c o n s t a n t . S i n c e MD r e s u l t s i n d i c a t e t h a t the mean-square t o t a l d i p o l e moment of i n t e r f a c i a l water i s g r e a t e r than t h a t f o r b u l k water ( 4 8 ) , the s t a t i c d i e l e c t r i c
2.
MULLA
Simulating
Liquid
Water near Mineral
31
Surfaces
c o n s t a n t of i n t e r f a c i a l water s h o u l d be lower than t h a t i n b u l k w a t e r . T h i s c o n c l u s i o n i s meant to be o n l y q u a l i t a t i v e , s i n c e c a l c u l a t i o n s of the s t a t i c d i e l e c t r i c c o n s t a n t u s i n g the Kirkwood t h e o r y may be i n c o n s i d e r a b l e e r r o r i f the e f f e c t s of the r e a c t i o n f i e l d are not accounted f o r ( 5 2 ) . The r e a c t i o n f i e l d i s an e l e c t r i c a l f i e l d w i t h i n the c u t o f f r a d i u s which c o n t r i b u t e s t o the d i e l e c t r i c c o n s t a n t . It r e s u l t s from the p o l a r i z i n g e f f e c t s o f the d i p o l e moments o u t s i d e the c u t o f f r a d i u s , and i t s e f f e c t can be i n c l u d e d ( w i t h c o n s i d e r a b l e e f f o r t ) i f a c c u r a t e computations of the d i e l e c t r i c c o n s t a n t are desired (53). The heterogeneous d i e l e c t r i c p r o p e r t i e s i n the l i q u i d medium near s u r f a c e s have i m p o r t a n t i m p l i c a t i o n s f o r t h e o r i e s of the s t r u c t u r e of the i n t e r f a c i a l r e g i o n . A r e c e n t t h e o r e t i c a l s t u d y of the e f f e c t of decreased medium d i e l e c t r i c c o n s t a n t near s u r f a c e s (54) shows t h a t i t i s associated with s i g n i f i c a n t reductions in surface potential computed from double l a y e r t h e o r y . A note of c a u t i o n to t h e o r e t i c i a n s is in order. The m o l e c u l a r s i m u l a t i o n s i n d i c a t e t h a t the s t a t i c d i e l e c t r i c p r o p e r t i e s of the i n t e r f a c i a l water decay g r a d u a l l y w i t h i n c r e a s i n g d i s t a n c e from the s u r f a c e . Hence, the use o f d i s c r e t e m i x t u r e models (54,55) i n which water near the s u r f a c e s i s d i v i d e d i n t o two z o n e s ; one having p r o p e r t i e s c h a r a c t e r i s t i c of water i n the f i r s t adsorbed l a y e r and the second h a v i n g b u l k p r o p e r t i e s are not l i k e l y to represent actual s u r f a c i a l c o n d i t i o n s . Similar caution s h o u l d be used i n a d s o r p t i o n models t h a t assume d i s c r e t e m o l e c u l a r l a y e r s of s u r f a c e complexes i n the i n t e r f a c i a l r e g i o n . R e a c t i o n k i n e t i c s and many t r a n s p o r t p r o p e r t i e s i n l i q u i d s are c o n t r o l l e d by r a t e s of d i f f u s i o n . P r o c e s s e s t h a t may be c o n t r o l l e d by d i f f u s i o n i n c l u d e , f o r example, r a t e s of l i g a n d exchange from t r a n s i t i o n metal i o n s ( 5 6 ) , r e a c t i o n s i n v o l v i n g p r o t o n t r a n s f e r ( 5 7 ) , and exchange r e a c t i o n s near m i n e r a l s u r f a c e s ( 5 8 ) . R e s u l t s i n Table V (from s t u d i e s 3 4 , 3 6 , 4 7 ) i n d i c a t e t h a t the v a l u e f o r the s e l f - d i f f u s i o n c o e f f i c i e n t of water near m i n e r a l s u r f a c e s i s c o n s i s t e n t l y about 80% lower than the v a l u e near the m i d p l a n e . Q u a n t i t a t i v e l y , the v a l u e s f o r D near the midplane are h i g h e r than v a l u e s r e p o r t e d f o r b u l k water at comparable temperatures ( 5 5 ) . For i n s t a n c e , the v a l u e s o f D i n b u l k water a t 285 and 300 Κ are about 1.8 and 2.9 χ 1 0 ~ m / s , respectively. The r e s u l t s from one study ( 3 3 ) , appear to be e x c e s s i v e l y h i g h , c o n s i d e r i n g the temperature of the s i m u l a t i o n , and are a l s o i n c o n s i s t e n t w i t h e x p e r i m e n t a l r e s u l t s i n t h a t they p r e d i c t d i f f u s i o n r a t e s which are g r e a t e r near the s u r f a c e s than near the midplane. E x c l u d i n g the l a t t e r r e s u l t s , v a l u e s f o r D from the MD s i m u l a t i o n s appear to q u a l i t a t i v e l y obey the expected t r e n d f o r i n c r e a s i n g v a l u e s of D w i t h i n c r e a s i n g t e m p e r a t u r e . F u r t h e r m o r e , the decreases i n v a l u e s f o r s e l f - d i f f u s i o n c o e f f i c i e n t near the s u r f a c e s are q u a l i t a t i v e l y c o n s i s t e n t w i t h e x p e r i m e n t a l measurements of decreased water m o b i l i t y near n e u t r a l s i l i c a t e s u r f a c e s ( 5 9 - 6 2 ) . 9
2
Another t r a n s p o r t p r o p e r t y of i n t e r f a c i a l water which can be s t u d i e d by MD t e c h n i q u e s i s the d i p o l e r e l a x a t i o n t i m e . This property i s computed from the d i p o l e moment c o r r e l a t i o n f u n c t i o n , which measures the r a t e a t which d i p o l e moment a u t o c o r r e l a t i o n i s l o s t due to r o t a t i o n a l motions i n time ( 6 3 ) . L a r g e r v a l u e s f o r the d i p o l e r e l a x a t i o n time i n d i c a t e slower r o t a t i o n a l motions o f the d i p o l e
32
G E O C H E M I C A L P R O C E S S E S AT M I N E R A L S U R F A C E S
moment. Both of the MD s t u d i e s r e p o r t i n g v a l u e s f o r the d i p o l e r e l a x a t i o n time i n Table V i n d i c a t e t h a t r e l a x a t i o n times are l a r g e r f o r water near the s u r f a c e s than f o r water near the m i d p l a n e . The MCY p o t e n t i a l , however, i s not as s e n s i t i v e t o changes i n r e l a x a t i o n time as i s the ST2 p o t e n t i a l . These r e s u l t s can be i n t e r p r e t e d to mean t h a t m o l e c u l e s near the s u r f a c e s e x p e r i e n c e h i n d e r e d r o t a t i o n a l movement as compared t o water near the m i d p l a n e . Experimental e v i d e n c e from near i n f r a r e d ( 6 ) , e l e c t r o n s p i n resonance ( 5 9 ) , and n u c l e a r magnetic resonance (60) s t u d i e s s u p p o r t the MD s i m u l a t i o n r e s u l t s , i n t h a t they i n d i c a t e h i n d e r e d r o t a t i o n a l motion of water near uncharged s i l i c a t e s u r f a c e s . Summary Monte C a r l o and M o l e c u l a r Dynamics s i m u l a t i o n s of water near hydrophobic s u r f a c e s have y i e l d e d a w e a l t h of i n f o r m a t i o n about the s t r u c t u r e , thermodynamics and t r a n s p o r t p r o p e r t i e s of i n t e r f a c i a l water. In p a r t i c u l a r , they have demonstrated the presence of m o l e c u l a r l a y e r i n g and d e n s i t y o s c i l l a t i o n s which extend many Angstroms away from the s u r f a c e s . These o s c i l l a t i o n s have r e c e n t l y been v e r i f i e d e x p e r i m e n t a l l y . Ordered d i p o l a r o r i e n t a t i o n s and reduced d i p o l e r e l a x a t i o n times are observed i n most of the s i m u l a t i o n s , i n d i c a t i n g t h a t i n t e r f a c i a l water i s not a u n i f o r m d i e l e c t r i c continuum. Reduced d i p o l e r e l a x a t i o n times near the s u r f a c e s i n d i c a t e t h a t i n t e r f a c i a l water e x p e r i e n c e s h i n d e r e d rotation. The m a j o r i t y of s i m u l a t i o n r e s u l t s i n d i c a t e t h a t water near hydrophobic s u r f a c e s e x h i b i t s fewer hydrogen bonds than water near the midplane. S e v e r a l m e r i t s and s t r e n g t h s of m o l e c u l a r s i m u l a t i o n s of i n t e r f a c i a l water are a p p a r e n t . S i n c e these methods y i e l d s t r u c t u r a l , thermodynamic and t r a n s p o r t p r o p e r t i e s f o r b u l k water which are i n good agreement w i t h many e x p e r i m e n t a l l y measured p r o p e r t i e s of b u l k water over a wide temperature range, they seem t o o f f e r a p r o m i s i n g approach f o r s t u d y i n g i n t e r f a c i a l w a t e r . I n t e r f a c i a l systems are g e n e r a l l y composed of s e v e r a l components and are d i f f i c u l t t o characterize. The n a t u r e of m o l e c u l a r s i m u l a t i o n s a l l o w s the system b e i n g a n a l y z e d to be e x a c t l y s p e c i f i e d i n terms of the types o f components, t h e i r i n t e r a c t i o n p o t e n t i a l s , the i n i t i a l atomic or m o l e c u l a r l o c a t i o n s and the types o f boundary c o n d i t i o n s imposed. Thus, e f f e c t s o f the s u r f a c e s can be s t u d i e d i n d e t a i l , s e p a r a t e l y from e f f e c t s o f c o u n t e r i o n s or s o l u t e s . In a d d i t i o n , i n d i v i d u a l l a y e r s o f i n t e r f a c i a l water can be a n a l y z e d as a f u n c t i o n of d i s t a n c e from the s u r f a c e and d i r e c t i o n a l a n i s o t r o p y i n v a r i o u s p r o p e r t i e s can be s t u d i e d . F i n a l l y , one computer experiment can o f t e n y i e l d i n f o r m a t i o n on s e v e r a l water p r o p e r t i e s , some o f which would be t i m e consuming or even i m p o s s i b l e t o o b t a i n by e x p e r i m e n t a t i o n . Examples of i n t e r f a c i a l water p r o p e r t i e s which can be computed v i a the MD s i m u l a t i o n s but not v i a experiment i n c l u d e the number of hydrogen bonds per m o l e c u l e , v e l o c i t y a u t o c o r r e l a t i o n f u n c t i o n s , and r a d i a l distribution functions. S e v e r a l weaknesses and d i s a d v a n t a g e s o f the computer s i m u l a t i o n methods can a l s o be mentioned. Foremost among these l i m i t a t i o n s i s the f a c t t h a t none o f the commonly used models f o r water i n t e r a c t i o n s
2.
MULLA
Simulating
Liquid
Water near Mineral
Surfaces
33
account f o r t h r e e - or f o u r - b o d y i n t e r a c t i o n s . Y e t , the complex b e h a v i o u r of i n t e r f a c i a l water i s a p p a r e n t l y d i c t a t e d by these cooperative e f f e c t s (26). S i g n i f i c a n t d i f f e r e n c e s e x i s t i n the i n t e r n a l energy, hydrogen bonding and d i f f u s i o n r a t e s o f i n t e r f a c i a l water as p r e d i c t e d by the ST2 and MCY p o t e n t i a l s , w i t h r e s u l t s f o r the former model being more c o n s i s t e n t w i t h e x p e r i m e n t a l r e s u l t s than those f o r the l a t t e r . These d i f f e r e n c e s p r o b a b l y r e l a t e to excess i n t e r n a l p r e s s u r e s generated by the MCY p o t e n t i a l . A second problem i n v o l v e s u n c e r t a i n t y about what d i s t a n c e o f s e p a r a t i o n s h o u l d be imposed between s u r f a c e s f o r a p a r t i c u l a r c h o i c e o f the number o f water m o l e c u l e s . I f e x c l u d e d volume near the s u r f a c e s i s not accounted f o r , the r e s u l t i n g average water d e n s i t y and i t s o s c i l l a t i o n s w i l l be e x c e s s i v e . F i n a l l y , d i e l e c t r i c p r o p e r t i e s of i n t e r f a c i a l water are not e a s i l y q u a n t i f i e d when r i g i d , p a i r w i s e i n t e r a c t i n g water p o t e n t i a l s are u s e d . T h i s l i m i t a t i o n can be p a r t i a l l y overcome by a c c o u n t i n g f o r r e a c t i o n f i e l d e f f e c t s , but even t h e n , water i n t e r a c t i o n p o t e n t i a l s i n c u r r e n t use are not e f f e c t i v e i n modeling the d i e l e c t r i c p r o p e r t i e s o f water ( 5 3 ) . On the w h o l e , the advantages and s t r e n g t h s of MC and MD s i m u l a t i o n s of i n t e r f a c i a l water outweigh t h e i r d i s a d v a n t a g e s and weaknesses. Even i f q u a n t i t a t i v e p r e d i c t i o n o f i n t e r f a c i a l water p r o p e r t i e s i s not p o s s i b l e i n some c a s e s , a knowledge of q u a l i t a t i v e t r e n d s as a f u n c t i o n o f d i s t a n c e from the s u r f a c e s or r e l a t i v e t o r e s u l t s from s i m u l a t i o n s o f b u l k water are o f t e n e x t r e m e l y illuminating. What i s the l i k e l y f u t u r e use o f MC and MD t e c h n i q u e s f o r s t u d y i n g i n t e r f a c i a l systems? S e v e r a l p r o m i s i n g approaches are possible. C o n t i n u e d i n v e s t i g a t i o n o f double l a y e r p r o p e r t i e s , " h y d r a t i o n f o r c e s " , " h y d r o p h o b i c e f f e c t s " , and " s t r u c t u r e d w a t e r " are c l e a r l y a w a i t i n g the development of improved models f o r w a t e r - w a t e r , s o l u t e - w a t e r , s u r f a c e - w a t e r , and s u r f a c e - s o l u t e p o t e n t i a l s . S i m u l a t i o n s o f o r g a n i c s near s u r f a c e s are c l e a r l y p o s s i b l e s i n c e p o t e n t i a l s d e s c r i b i n g w a t e r - o r g a n i c i n t e r a c t i o n s are p r e s e n t l y available (64,65). S t u d i e s o f r e a c t i o n k i n e t i c s and l i g a n d exchange p r o c e s s e s near s u r f a c e s are c l e a r l y p o s s i b l e u s i n g the m o l e c u l a r t i m e s c a l e g e n e r a l i z e d Langevin e q u a t i o n approach ( 6 6 ) . Gaseous a d s o r p t i o n on metal ( 6 7 ) , hydrophobic (68) and o t h e r s i m p l e s u r f a c e s have been e x t e n s i v e l y s t u d i e d ( 2 ) , but s i m i l a r approaches u s i n g models f o r c l a y and z e o l i t e s u r f a c e s are a l s o p o s s i b l e . Mechanisms o f c r y s t a l growth and d e f e c t s t r u c t u r e s i n v i t r e o u s s i l i c a (69) and g l a s s (70) have been s t u d i e d , and s i m i l a r s t u d i e s on a l u m i n o s i l i c a t e m i n e r a l s (even under c o n d i t i o n s of h i g h temperature and p r e s s u r e ) are p o s s i b l e . Finally, new t h e o r e t i c a l developments are a l l o w i n g thermodynamic p r o p e r t i e s (71,72) and n o n - e q u i l i b r i u m c o n d i t i o n s (73) to be s t u d i e d w i t h MD methods. In s h o r t , MC and MD s t u d i e s of i n t e r f a c i a l systems are s t i l l in t h e i r infancy.
34
GEOCHEMICAL PROCESSES AT MINERAL SURFACES
L i s t o f Symbols a A b C d j D Ε v
a
k m-j Ν Pq ρ Ρ qj q Q R-jj S(R.jj) Τ U U-jj U|_j U|\j U 5 v.j V x-j μ σ At a
W
parameter f o r t h e Lennard-Jones p o t e n t i a l parameter f o r t h e Lennard-Jones p o t e n t i a l parameter f o r t h e Lennard-Jones p o t e n t i a l c o n s t a n t volume heat c a p a c i t y d i s t a n c e between s u r f a c e and a t h charge on water m o l e c u l e self-diffusion coefficient t o t a l energy f o r c e a c t i n g on i t h atom Boltzmann's constant mass o f i t h atom number o f m o l e c u l e s value of pth property of qth configuration ensemble average v a l u e o f p r o p e r t y ρ i n t e r n a l pressure charge on j t h s u r f a c e atom charge on a t h charge c e n t e r o f water m o l e c u l e c a n o n i c a l ensemble p a r t i t i o n f u n c t i o n d i s t a n c e between i t h and j t h atoms s w i t c h i n g f u n c t i o n o f ST2 water p o t e n t i a l a b s o l u t e temperature i n t e r n a l energy pairwise interaction potential Lennard-Jones i n t e r a c t i o n p o t e n t i a l t o t a l i n t e r a c t i o n p o t e n t i a l energy p o t e n t i a l energy f o r water i n t e r a c t i n g w i t h a charged s u r f a c e v e l o c i t y o f i t h atom volume p o s i t i o n o f i t h atom d i p o l e moment parameter f o r t h e Lennard-Jones p o t e n t i a l time s t e p f o r MD a l g o r i t h m
Literature Cited 1. 2.
Cairns-Smith, A. G. Scientific American 1985; 252(6), 90-100. Nicholson, D.; Parsonage, N. G. "Computer Simulation and the Statistical Mechanics of Adsorption"; Academic: New York, 1982; Chap. 4,6. 3. van Megen, W.; Snook, I. Adv. Colloid Interface Sci. 1984; 21, 119-194. 4. Low, P. F. Adv. Agronomy 1961, 13, 269-327. 5. Graham, J. Rev. Pure Appl. Chem. 1964; 14, 81-89. 6. Klier, K.; Zettlemoyer, A. C. J. Colloid Interface Sci. 1977; 58(2), 216-229. 7. Low, P. F. Soil Sci. Soc. Am. J. 1979; 43(5), 651-658. 8. Derjaguin, Β. V.; Churaev, Ν. V. In "Progress in Surface and Membrane Science"; Cadenhead, D. Α.; Danielli, J. F., Eds.; Academic: New York, 1981; Vol. 14, pp 69-130. 9. Sposito, G.; Prost, R. Chem. Rev. 1982; 82(6), 553-573. 10. Void, R. D.; Void, M. J. "Colloid and Interface Chemistry"; Addison-Wesley: Reading, Mass., 1983; Chap. 7,9.
2. MULLA 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. 34. 35. 36. 37. 38. 39. 40.
Simulating Liquid Water near Mineral Surfaces
35
Valleau, J. P.; Whittington, S. G. In "Modern Theoretical Chemistry"; Berne, B. J., Ed.; Plenum: New York, 1977; Vol. 5, Chap. 4. Barker, J. A.; Henderson, D. Rev. Mod. Phys. 1976; 48(4), 587671. Wood, D. W. In "Water: A Comprehensive Treatise"; Franks, F., Ed.; Plenum: New York, 1979; Vol. 6, Chap. 6. Snook, I. K.; van Megen, W. J. Chem. Phys. 1980; 72(5), 2907-2913. Verlet, L. Phys. Rev. 1967; 159(1), 98-103. Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H. J. Chem. Phys. 1953; 21(6), 1087-1092. Alder, B. J.; Wainwright, T. E. J. Chem. Phys. 1957; 27, 12071209. Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J. Comp. Phys. 1977; 23, 327-341. Steele, W. A. Surface Sci. 1973; 36, 317-352. Rahman, A. Phys. Rev. 1964; 136(2A), A405-A411. Barker, J. Α.; Watts, R. O. Chem. Phys. Lett. 1969; 3(3), 144145. Rahman, Α.; Stillinger, F. H. J. Chem. Phys. 1971; 55(7), 33363359. Stillinger, F. H., Rahman, A. J. Chem. Phys. 1974; 60(4), 15451557. Matsuoka, O.; Clementi, E.; Yoshimine, M. J. Chem. Phys. 1976; 64, 1351-1361. Stillinger, F. H.; Rahman, A. J. Chem. Phys. 1978; 68(2), 666-670. Barnes, P.; Finney, J. L.; Nicholas, J.; Quinn, J. E. Nature 1979; 282, 459-464. Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983; 79, 926-935. Detrich, J.; Corongiu, G.; Clementi, E. IBM Res. Rept. 1984; KGN3, 1-11. Impey, R. W.; Madden, P. Α.; McDonald, I. R. Mol. Phys. 1982; 46(3), 513-539. van Gunsteren, W. F.; Berendsen, H. J. C.; Rullman, J. A. C. Faraday Disc. 1978; 66, 58-70. Stillinger, F. H.; Rahman, A. J. Chem. Phys. 1978; 68(2), 666-670. Rahman, Α.; Stillinger, F. H. J. Chem. Phys. 1971; 55(7), 33363359. Sonnenschein, R.; Heinzinger, K. Chem. Phys. Lett. 1983; 102(6), 550-554. Mulla, D. J.; Cushman, J. H.; Low, P. F. Water Resour. Res. 1984; 20(5), 619-628. Christou, Ν. I.; Whitehouse, J. S.; Nicholson, D.; Parsonage, N. G. Faraday Symp. Chem. Soc. 1981; 16, 139-149. Barabino, G.; Gavotti, C.; Marchesi, M. Chem. Phys. Lett. 1984; 104(5), 478-484. Hobza, P.; Sauer, J.: Morgeneyer, C.; Hurych, J.; Zahradnik, R. J. Phys. Chem. 1981; 85, 4061-4067. Geerlings, P.; Tariel, N.; Botrel, Α.; Lissillour, R.; Mortier, W. J. J. Phys. Chem. 1984; 88, 5752-5759. Sauer, J . ; Zahradnik, R. Intl. J. Quant. Chem. 1984; 26, 793-822. Low, P. F. Soil Sci. Soc. Am. J. 1981; 45(6), 1074-1078.
36
GEOCHEMICAL PROCESSES AT MINERAL SURFACES
41.
Probst, M. M.; Radnai, T.; Heinzinger, K.; Bopp, P.; Rode, B. M. J. Phys. Chem. 1985; 89, 753-759. 42. Heinzinger, K. Pure Appl. Chem. 1985; 57(8), 1031-1042. 43. Bounds, D. G. Mol. Phys. 1985; 54(6), 1335-1355. 44. Gruen, D. W. R.; Marcelja, S.; Pailthorpe, B. A. Chem. Phys. Lett. 1981; 82(2), 315-320. 45. Jonsson, B. Chem. Phys. Lett. 1981; 82(3), 520-525. 46. Lee, C. Y.; McCammon, J. Α.; Rossky, P. J. J. Chem. Phys. 1984; 80(9), 4448-4445. 47. Marchesi, M. Chem. Phys. Lett. 1983; 97(2), 224-230. 48. Mulla, D. J.; Low, P. F.; Cushman, J. H.; Diestler, D. J. Colloid Interface Sci. 1984; 100(2), 576-580. 49. Pashley, R. M.; Israelachvi1i, J. N. J. Colloid Interface Sci. 1984; 101(2), 511-523. 50. Gruen, D. W. R.; Marcelja, S. J. Chem. Soc. Faraday Trans. 2 1983; 79, 225-242. 51. Goldman, S. J. Chem. Phys. 1981; 75(8), 4064-4076. 52. Steinhauser, O. Mol. Phys. 1982; 45(2), 335-348. 53. Neumann, M. J. Chem. Phys. 1985; 82(12), 5663-5672. 54. Helmy, A. K.; Natale, I. M. Clays Clay Min. 1985; 33(4), 329-332. 55. Fripiat, J. J.; Letellier, M.; Levitz, P. Phil. Trans. R. Soc. London 1984; A311, 287-299. 56. Margerum, D. W.; Cayley, G. R.; Weatherburn, D. C.; Pagenkopf, G. K. In "Coordination Chemistry"; Martell, A. E., Ed.; ACS Monograph 174, American Chemical Society: Washington, D. C., 1978; Vol. 2, Chap. 1. 57. Eigen, M. Angew. Chemie 1964; 3(1), 1-72. 58. Breen, C.; Adams, J.; Riekel, C. Clays Clay Min. 1985; 33(4), 275284. 59. Bassetti, V.; Burlamacchi, L.; Martini, G. J. Am. Chem. Soc. 1979; 101(19), 5471-5477. 60. Cruz, M. I.; Letellier, M.; Fripiat, J. J. J. Chem. Phys. 1978; 69(5), 2018-2027. 61. Martini, G. J. Colloid Interface Sci. 1981; 80(1), 39-48. 62. Roberts, Ν. K.; Zundel, G. J. Phys. Chem. 1980, 84, 3655-3660. 63. Mulla, D. J.; Low, P. F. J. Colloid Interface Sci. 1983; 95(1), 51-60. 64. Remerie, K.; van Gunsteren, W. F.; Postma, J. P. M.; Berendsen, H. J. C.; Engberts, J. B. F. N. Mol. Phys. 1984; 53(6), 1517-1526. 65. Kuharski, R. Α.; Rossky, P. J. J. Am. Chem. Soc. 1984; 106, 57865793. 66. Adelman, S. A. J. Phys. Chem. 1985; 89, 2213-2221. 67. Broughton, J. Q. Surface Sci. 1980; 91, 91-112. 68. Talbot, J.; Tildesley, D. J.; Steele, W. A. Mol. Phys. 1984; 51(6), 1331-1356. 69. Garofalini, S. H. J. Non-Crystalline Solids 1984; 63, 337-345. 70. Soules, T. F. J. Non-Crystalline Solids 1982; 49, 29-52. 71. Nose', S. Mol. Phys. 1984; 52(2), 255-268. 72. Ray, J. R.; Rahman, A. J. Chem. Phys. 1984; 80(9), 4423-4428. 73. Hoover, W. G.; Moran, B.; Haile, J. M. J. Stat. Phys. 1984; 37(1/2), 109-121. RECEIVED
June 18, 1986
