A molecular dynamics study on surface properties of supercooled water

616Science in China Series G: Physics, Mechanics & Astronomy 2006 Vol.49 No.5 616- -625DOI: 10.1007/s11433-006 2019-6A molecular dynamics study on surfaceproperties of supercooled waterL0 Yongjun & WEI BingboDepartment of Applied Physics, Northwestern Polytechnical University, Xi'an 710072, ChinaCorrespondence should be addressed to Wei Bingbo (email: bbwei@ nwpu.edu.cn)Received May 10, 2005; accepted July 14, 2006Abstract Molecular dynamics simulations were performed to study the surface proper-ties of water in a temperature range from 228 to 293 K by using the extended simple pointcharge (SPC/E) and four-site TIP4P potentials. The calculated surface tension increaseswith the decrease of temperature, and moreover the slopes of the surface ten-sion-temperature curves show a weak rise below 273 K, whereas no obvious anomaliesappear near 228 K, which accords with the previous experiments. Compared with themeasured values, the SPC/E potential shows a good agreement, and the TIP4P potentialunderestimates the surface tension. The main reason for that may be the reasonable de-scription of the surface structure of supercooled water for the SPC/E. When simulating theorientational distributions of water molecules near the surface, the SPC/E potential pro-duces higher ordering and larger surface potentials than the TIP4P potential.Keywords: surface tension, supercooled water, molecular dynamics simulation.Liquid-vapor interface has aroused wide research interest in experimental and theo-retical fields in recent several years because of its importance in physics, biology andenvironmental sciencel- 5. As a transition phase between the two bulk systems, liq-uid-vapor interface displays different physical and chemical properties. In particular, Wa-ter, the most popular liquid in nature, has especially large surface tension compared tosome simple liquids owing to the hydrogen-bonding structure at the surface. Therefore,the better understanding of thermodynamic and dynamic properties associated with liq-uid-vapor interface requires a molecular-level investigation of surface structure.Molecular dynamics (MD) simulation is one of the effective methods to study interfa-cial phenomena in molecular level. In the past decade, the surface properties of waterhave been subjected to extensive investigation by MD method using various interactionpotentials7- 9. Guissani and Guillot studied the liquid vanr cnexistence curve of SPC/Bwater and discussed the relationship between中国煤化工:rater moleculesand the slope of the liquid-vapor coexistence curHCNMHGtedthesufacetension and surface potential of water at 325 K for the TIP4P model8l. The correctedwww.scichina.com www.springerlink.comA molecular dynamics study on surface properties of supercooled water617value of the surface tension is 59 mNm , which is 13% smaller than the experimentalresults. And the calculated surface potential has reasonable magnitude but opposite sign.Alejandre et al.!9l investigated the liquid vapor equilibrium of water in the range of 316一573 K by using the SPC/E model with the Ewald sum. The simulated surface tensions arein good agreement with the experimental results, which suggests that the SPC/E model is amore effective potential function for studying the liquid vapor interface of water.It is noted that most studies on the surface properties of water concentrate on hightemperature region, and less attention is paid to supercooled water. Although Matsumotoand Kataokalu studied the surface properties of water in a wide temperature range from250 to 400 K with MD method, including surface tension, surface excess energy andorientational structure near the surface, the simulated surface tensions are 50% smallerthan the experimental values, and moreover, the specific surface entropy Ss(s。=-dy/dT， γis the surface tension, T the temperature) diminishes monotonouslywith the decrease of temperature in the supercooled region, which is contradictory to theexperimental results by Floriano et al.". Floriano et al.'measured the surface tensionof water using a small-sample capillary rise method down to 245.8 K and the resultsshowed that s is inclined to increase below 273 K. Therefore, it is necessary to investi-gate the surface tension of supercooled water in detail.In this paper, we report a study on the surface properties of supercooled water down to228 K with a classical MD method, and the behaviors of the surface tension and surfacestructure in the supercooled region are discussed.1Simulation methodDuring the simulations, the SPC/E and 4 site TIP4P potential functions are utilized.The interaction between two water molecules is a combination of Lennard-Jones 6-12potential function and Coulomb potential:U(n)=4ejt qi9j.(1])where r;j is the distance between atom i and j in two different molecules; qi and qj are thecharges of atom i and j; Ej and σj are Lennard-Jones parameters. The potential parame-ters for the SPC/E and TIP4P models are given in Table 1.Table 1 Molecular potental parameters for the SPC/B and TIP4P potential functios12,13)SPCETIP4ProH(A)0.9572∠HOH (deg)109.47104.52σ(A)3.1663.154E(kJ-mor)0.65020.6481go (e)-0.8746中国煤化工9(e)0.4238MHC NMH Goqm(e)roM (A)0.15618Science in China Series G: Physics, Mechanics & AstronomyThe simulation domain is schematically shown in Fig. 1. Periodical boundary condi-tions are applied in the three coordinate directions. The simulation system consists of 512.water molecules. The simulation is performed with a time step of 2.5 fs, and the simula-tion temperature ranges from 228 to 293 K. The short-range potential cutoff is 9.8 A. Ateach temperature, the bulk water system without liquid-vapor interface is equilibrated atfirst and then the z dimension of the simulation cell is increased to 100 A in order to formtwo liquid-vapor interfaces. The new system is equilibrated in a NVT ensemble rangingfrom 10 to 5x10' time steps with the decrease of simulation temperature. The productionruns of 10’ time steps are required to calculate the surface properties accurately. All ofhe codes run in HP-rx2600 system at Center for High Performance Computation ofNPU.VaporLiquidX业l。Fig.1. Simulation cell. L= 100 A for N = 512. The middle part is liquid water and vapor phase is fixed in eachside.In the MD simulations of dipolar fluids such as water, the long range electrostatic in-teractions play an important role and cannot be neglected. Therefore, the Coulomb poten-tial consists of two parts: short- range and long -range forces. The Ewald sum is an effec-tive way to calculate the ionic interaction and it can embody the long-range force in asimulation with periodic boundaries4. By means of the Ewald sum, the total Coulombpotential Ucoul is given byerf(KTiajb)UCoulrexp|(h) +;ZZZqaZqmhz0h'4x2ZTajtib(2)erf(KTrai)元2Zq -2222qw9nTaiba b#awhere V is the volume of simulation cell, h the reciprocal lattice vector, K(= 5.6/Lx) theconvergence factor, erfc(x) the complementary error function, S(h)= 22 Gia exp(ih.ra)中国煤化工and riajb =ra -rpb. In the simulations, a cutoff ispace, and theMYHCNMHGmaximum reciprocal lattice vectors are given by Ihmax H hmax =5. Because the size ofA molecular dynamics study on surface properties of supercooled water619simulation cell in the z: direction is lengthened, the maximum of Ihax | should be in-creased accordingly and is equal to 15 in this work.2 Results and discussion2.1 Density profilesIn order to observe the formation of liquid-vapor interface, the distribution of numberdensity along z coordinate is calculated. The simulation cell is divided into slabs with0.25 A thickness, and molecule number is statistically averaged in each slab. The densityprofiles at 233 and 293 K for the SPC/E and TIP4P potentials are shown in Fig. 2(a) and(b) respectively. There exists obviously a continuous density drop between the bulk liq-uid and vapor phases, indicating that the liquid-vapor interface forms. The density pro-files at every temperature are found to fluctuate around the bulk density. The amplitudeof this fluctuation reaches a maximum at the liquid/vapor interface and then falls intooscillations with a weak decay in liquid region. Moreover, the fluctuation becomes moreand more pronounced with the decrease of temperature, as shown in the insets of Fig. 2.The oscillatory density profile is a result of surface layering structure which is often ob-served in some liquids with low melting point and regarded as a simple geometrical re-1.2-(a)SPC/E0293K233 K- Fiting0.8-51.04+Normalg0.4Supercooled0.92L1.2-(6)TIP4P。293Kg 0.8--Fitting1.0420.960926古12(A)0.01015中国煤化工Z(A)THCNMH GFig. 2. Density profiles at 293 K (open circles) and 233 K (dots) for the SPC/E and TIP4P potentials. The solidlines are the fiting results according to the byperbolic tangent function. The insets show the enlargement of liquiddensity profiles.620Science in China Series G: Physics, Mechanics & Astronomysponse to a sharp interface with a short-range correlation. At the free liquid surface, thecapillary wave associated with the long-range correlation often suppresses the inho-mogeneity of density profile and damps the oscillation. Its effect on the suface layeringdepends on the transverse size. In simulations, the transverse size of the simulation cell isvery small and the effect of capillary waves is often ignored. In fact, according to theanalyses of Tarazona et al.!5 ' the layering structure is independent of temperature, andthe temperature dependence of the oscillatory density profiles shown in this simulation iscaused by the capillary wave essentially. The effect of temperature on the oscillatory am-plitude Am is described by means of the transverse sizeL: Am ~ L7(T.r) (7 is the decayexponent, a function of temperature and surface tension). Therefore, the increase of theinterface area (LxL) can damp the density oscillation to a certain extent, but it may notbe pronounced in view of the limitation of the typical scope in the present simulations.The density profiles are ftted by a hyperbolic tangent function:p(z)=-(Pr + Pv)-j(AL- Pv)tanb[(z- zo)1d],(3)where PL and Pv are bulk densities of liquid and vapor phase, Zo is the position of theGibbs surface, and d is a parameter of surface thickness. The fitted p (z) is shown as solidlines in Fig. 2.The surface thickness t, defined as the“10- 90 thickness"', is a function of d, t =2.1972d. For both the SPCE and TIP4P potentials, the surface thickness decreases withthe decrease of temperature as shown in Fig. 3. Moreover, in the simulation temperaturerange, the surface thickness for the TIP4P is about 5% larger than the SPC/E potential.According to Schwartz's measurements'o, the surface thickness of water reaches 7.2 Aat 293 K, much larger than the values of 2.3 and 2.4 A obtained in this work. The differ-ence between the simulations and experiments is attributed to the surface capillary wavesinvolved in the experimental measurements. The present results for the SPC/E model isalso smaller than that of Matsumoto et al. for the CC potentialo, which may relate to themolecular structures defined by these potential functions. The CC and TIP4P potentials2.5 t-▲SPC/E一●TIP4P三2.0-中国煤化工.5HCNMHG220240260280300T(K)Fig. 3.“10 -90 thickness" as a function of temperature for the SPC/E and TIP4P potential functions.A molecular dynamics study on surface properties of supercooled water621yield weaker interaction force, smaller liquid densities and higher vapor densities, result-ing in the larger surface thickness.2.2 Surface tensionSurface tension is one of the important parameters to characterize the surface structureand chemical activity and associates closely with some surface anomalies. In order tounderstand quantitatively the surface properties of water in the supercooled region, wesimulated the surface tension of water as a function of temperature from 228 to 293 K.The surface tension can be obtained by calculating the components of the pressuretensorl17:A/12\2(Px +P)-Pz(4)where Px Pyy and Pzz are the X, Y and Z element of the pressure tensor, andA = LLy issurface area. The element of the molecular pressure tensor isVPap=Z,m(V)a(V)p+ 2(q)a(faqi)p,(5)=l j>i a,b=lwhere N is the number of molecules; mi; V; are the molecular mass and the velocity of thecenter of mass; u the number of sites in a water molecule; r'ij the distance between themolecules i and j; andf the force between atom a in molecule i and atom b in moleculej.The electrostatic long range interactions make an important contribution to the surfacetension and are involved by using the Ewald sum in the simulations based on eq. (2). Thetail correction to the surface tension is also performed.The surface tensions of supercooled water as a function of temperature for the SPC/Eand TIP4P models are shown in Fig. 4. The open circles are the measured values byFloriano et al."; and the solid lines are the best fit for the simulated results. The ex-perimental data above 273 K have been measured extensively and can be accurately fit-ted by Vargafik equation[18!:「T。-T](T。-7γ=B1+b(6)1where Tc= 647.15 K, B= 0.2358 mNm~", h= 1.256 and b = -0.625. In general, the sur-face tension in the supercooled region is obtained approximately from the extrapolationof eq. (6), shown as the dashed line in Fig. 4. It can be seen that the absolute values of theslope of temperature surface tension curve diminish as the temperature decreases. How-ever, the experimental data given by Floriano display that the temperature dependence ofthe surface tension begins to increase below 273 K. For the present simulated results, thesurface tensions are in good agreement with the measured values in the supercooled statefor the SPC/E potential, and are somewhat smaller while T> 273 K. Once the tempera-ture decreases below 273 K, the temperature depe中国煤化Insion starts toslowly rise and the specific surface entropy, ss,MHCNMHGwith the de-crease of temperature, as plotted in Fig.4, which is more eviaent Ior tne 1 IP4P potential.This change tendency agrees qualitatively with the experimental results and is contrary to622Science in China Series G: Physics, Mechanics & AstronomyMatsumoto and Kataoka's results. It is suggested that there may be a second inflectionpoint in the surface tension temperature curve in the low temperature region as the firstone at about 200Cl9. Although the present results show that the temperature depend-ence of the surface tension is strengthened below 273 K, the increase appears so slowthat we cannot confirm the location of the inflection point from the relatively smallersimulation data. And the more extensive simulations in a wide temperature range, espe-cially at larger surpercooling, are required.100 r◆SPC/ETIP4P-。. Fittedo Floriano et al!". Vargaftik et al!8]80F”mpo,0-△-s60SupercooledNormalTm =273 K40 L220240260280300T(K) .Fig. 4. Surface tension of supercooled water. The solid diamonds and open triangles are calculated values with theSPC/E and TIP4P potentials, respectively. The open circles are measured results in ref. [11]. The solid curves are thebest fit to the simulated results. The dash line is the prediction of Vargaftik equation.In addition, no anomalous behavior of the surface tension like some bulk properties ofsupercooled water appears when the temperature approaches 228 K[20]. Florianol1l ex-plained that the divergence of the surface free energy per mole, γ213, at 228 K ismainly caused by the mole volume Vm, independent of x However, it is near 228 K thatthe anomalous behavior of γ2 3 appears, and if the behavior of the surface tension isunknown at 228 K, the above interpretation looks inadequate. The present simulationsprovide a complete picture of the surface tension near 228 K and can more rigorouslyelucidate this question.The comparison of the surface tension between the SPC/E and TIP4P potentials dem-onstrates that the SPC/E potential is a better model for reproducing the surface character-istics of water in the supercooled state. Considering a good agreement with the experi-mental values in Alejandre's simulationsl up to the triple point, the SPC/E potential canquantitatively represent the surface tension of water over a wide temperature range.中国煤化工2.3 Orientational structure and surface potentialMH.CNMHGThe orientation of water molcules is determinea oy two angies σ anu中according tothe definition of Matsumoto et al.0. θis the angle between the surface normal vector nzA molecular dynamics study on surface properties of supercooled water623and molecular dipole, and中is the rotational angle around the molecular dipole. It isfound that there are preferential orientation angles for SPC/E and TIP4P water moleculesnear the surface. The statistical averages of the two angles in the liquid side, <θ> and<办>, are 103° and 62° at 293 K for the SPC/E potential. In the vapor side, <的> and <φ>are 75° and 37° respectively. It can be seen that the dipole of water molecules in the liq-uid side points to the bulk liquid and one hydrogen atom of the molecule projects to-wards the liquid phase. On the other hand, the molecular dipole in the vapor side prefersto point to the bulk vapor, and one hydrogen atom projects towards the vapor phase.Moreover, these angles are independent of temperature in the simulations. For the TP4Ppotential, the values of <θ> and <φ> are smaller than that for the SPC/E potental, 95°and 54° in the liquid side, 84° and 38° in the vapor side respectively. Clearly, the degreeof orientation ordering of SPC/E water molecules near the surface is higher, which de-termines stronger surface polarity of water.The surface potential x is defined as the difference of the electrostatic potential be-tween the bulk liquid and vapor phases. It is calculated according to the following rela-tion:x=丛r ,dzp(z)(cos8),(7)where μ is the dipole moment of water molecules, equal to 2.351D for SPC/E and2. 180D for the TIP4P potential (1D = 3.335x10 30 Cm); E is the dielectric constant of thevacum. The calculated results are shown in Fig. 5. The surface potentials for the SPC/Eand TIP4P models are -0.58 V and - 0.39 V at 293 K, respectively. Compared with thevalue of -0.55 V reported by Zakharov et al.5l at 300 K using the TIP4P, the presentTIP4P results are large, which is relevant to the definition of surface potential. In Zak-harov et al's simulation', the surface potential consists of dipolar and quadrupolar con-tributions, whereas only dipolar contribution is taken into account in this work. The ab-solute values of the surface potential show an increasing tendency with the decrease oftemperature. If the surface potential is regarded as a linear function of temperature in thesimulation temperature range, its temperature coefficients (dydT) are equal to 13 and 12mV/K for the SPC/E and TIP4P potentials respectively. The values have opposite sign toSchiffrin's measurement!. In fact, the surface potential has great uncertainty in experi-ments and simulations, and even the results vary signs, which is induced by the measuremethod, potential function and calculation technique. As few experimental results can beobtained in the supercooled region to test the simulations, our aim is only to provide cal-culated values to predict the surface potential of supercooled water. In addition, the sur-face potential-temperature curves display a flat near 273 K as illustrated in Fig. 5.Whether this behavior corresponds to a transition of molecular structure near surface or isrelevant to the second inflection point of the surface tensinn is etil1 nnt clear. The abso-lute values of the surface potential for the SPC/EIYR中国煤化工than that forthe TIP4P, which is consistent with the above anaC N M H Gistribution ofsurface molecules and the result of the molecular structure determined by potential models.624Science in China Series G: Physics, Mechanics & Astronomy-▲SPCE-0.5- -●. TIP4P己-1.0-心▲一-1.5240260280300T(K)Fig.5. Surface potentials of supercooled water as a function of temperature. The triangles and squares are resultsfor the SPC/E and TIP4P potentials respectively.3 ConclusionsMD simulations are accomplished to study the surface properties of supercooled waterfrom 228 to 293 K. For the SPC/E potential, the simulated surface tension agrees wellwith the experimental results. Moreover, the temperature dependence of the surface ten-sion for both the SPC/E and TIP4P potentials tends to increase in the supercooled region,which is consistent with the experiments. The simulated orientational structure of watermolecules indicates that there exists the orientational ordering near the surface. The cal-culated surface potential displays negative value and positive temperature coefficient.Combined with the above results, the SPC/E potential reproduces accurately the surfacetension, causes stronger orientational ordering of water molecules near the surface, andgives a better description of the surface properties of supercooled water.AcknowledgementsThis work is financially supported by the National Natural Science Founda-tion of China (Grant Nos. 50121101, 50395105 and 50271058), and the Doctorate Foundation ofNorthwestern Polytechnical University. The authors would like to thank the Center for High Per-formance Computation of Northwestern Polytechnical University.ReferencesBasu J K, Hazra s, Sanyal M K. Growth mechanism of Langmuir Blodgett films. Phys Rev Lett, 1999, 82:4675 - 46782 Taylor R s, Shields R L. Molecular dynamics simulations of the ethanol liquid-vapor interface. J Chem Phys,2003, 119: 12569- 125763 Velev 0 D, Gurkov T D, Ivanov I B, et al. Abnormal thickness and stability of nonequilibrium liquid films.Phys Rev Lett, 1995, 75: 264- - 267中国煤化工4 Weng J G Park s, Lukes J R, et al. Molecular dynamics iYHCNMHGon liquid fims. JChem Phys, 2000, 113: 5917- 59235 Zakharov V V, Brodskaya E N, Laaksonen A. Surface tension of water droplets: A molecular dynamics study ofmodel and size dependencies. J Chem Phys, 1997, 107: 10675- 10683A molecular dynamics study on surface properties of supercooled water6 Wang J Z Cben M, Guo Z Y. A two-dimensional molecular dynamics simulation of liquid-vapor nucleation.Chin Sci Bull, 2003, 48(7): 623- 6267 Guissani Y, Guillot B. A computer simulation study of the liquid-vapor coexistence curve of water. J ChemPhys, 1993, 98: 8221 - 82358 Wilson M A, Porile A, Pratt L R. Surface potential of the water liquid -vapor interface. J Chem Phys, 1988,88: 3281一32859 Alejandre J, Tildesley D J, Chapela G A. Molecular dynamics simulation of the orthobaric densities and surfacetension of water. J Chem Phys, 1995, 102: 4574- 45830 Matsumoto M, Kataoka Y. Study on liquid- vapor interface of water ();: Simulational results of thermodynamicproperties and orientational structure. J Chem Phys, 1988, 88: 3233- 32451 Floriano M A, Angell C A. Surface tension and molar surface free energy and entropy of water to- -27.2C.JPhys Chem, 1990, 94: 4199- 420212 Jorgensen w L, Chandrasekhar I, Madura J D. Comparison of simple potential functions for simulating liquidwater. J Chem Phys, 1993, 79: 926- 93513 Berendsen HJ C, Grigera J R, Stratsma T P. The missing term in efctive pair potentials. J Phys Chem, 1987,91: 6269-627114 Arbuckle B W, Clancy P. Effects of the Ewald sum on the free energy of the extended simple point chargemodel for water. J Chem Phys, 2002, 116: 5090- 509815 Tarazona P, Chacon E, Reinaldo-Falagan M, et al. Layering structures at free liquid surfaces: TheFisher- Widom line and the capillary waves. J Chem Phys, 2002, 117: 3941- 395016 Schwartz D K, Scholssman E K, Kellogg G J, et al. Thermal diffuse X-ray-scattering studies of the water-vaporinterface. Phys Rev A, 1990, 41: 5687- 569017 NijmeijerN J P, Bakker A F, Bruin C, et al. A molecular dynamics simulation of the Lennard-Jones liq-uid-vapor interface. J Chem Phys, 1988, 89: 3789- 379218 Vargaftik N B, Volkov B N, Voljak L D. International tables of the surface tension of water. J Phys Chem RefData, 1983, 12:817- -82019 Pellicer J, Garcia-Morales V, Guanter L, et al. On the experimental values of the water surface tension used insome textbook. Am J Phys, 2002, 79: 705 - 70920 Leyendekkers J V, Hunter R J. Thermodynamic properties of water in the subcooled region. J Chem Phys, 1982,82: 1447- 145321 Shifrin D J. Real standard entropy of ions in water. Trans Faraday Soc, 1970, 66: 2464- -2468中国煤化工MYHCNMHG