The article information
 Miao Han, Ma Hongru
 缪晗, 马红孺
 Confinement Induced Ordering in Fluid of Hard Ellipsoids
 受限导致硬椭球液体的取向有序
 Chinese Journal of Chemical Physics, 2016, 29(2): 212218
 化学物理学报, 2016, 29(2): 212218
 http://dx.doi.org/10.1063/16740068/29/cjcp1506119

Article history
 Dated: Received on June 7, 2015
 Accepted on October 18, 2015
b. School of Mechanical Engineering, and Key Laboratory of Arti cial Structures and Quantum Control (Ministry of Education), Shanghai Jiao Tong University, Shanghai 200240, China
Hardsphere (HS) fluid [1, 2, 3, 4, 5, 6, 7] is an ideal theoretical model that can be used to describe the fluid phase of matter,and it is also the simplest one that predicts a phase transition between a fluid phase and a solid phase. The HS model has been studied extensively for more than four decades,and thus deep understanding the various aspects of the model has been well established. Among these,one of the most important aspects is about the impact of the hardwall confinement on the properties of the HS fluid. For example,Götzelmann et al statistically computed the direct correlation function and the surface tension of uniform hardsphere mixtures within the frame of the weighted density functional theory [8]. Later on,Roth et al [9] studied the structures of polydisperse hardsphere systems by using modified fundamental measure theory (FMT) originally proposed by Rosenfeld [10],and predicted the existence of the oscillatory size segregation. On the other hand,Lin et al experimentally investigated the hindered diffusion of an isolated Brownian particle system confined within two flat walls and implied that ignoring the pressure and velocity fields led to underestimation of the influence of the walls on the dragging force [11].
Since 1960s an increasing interest has been focused on the study of isotropically interacting hardbody systems,which has led to a lot of achievements such as applications of the scaling particle theory [12] and a vast knowledge of correlation function [13],whereas these studies are still based on single direct assumptions. In many real systems the effects caused by the anisotropic shape of particles cannot be ignored [14, 15, 16, 17, 18, 19]. As a matter of fact,many anisotropic systems,such as rods [20],spherocylinders [21],and plates [22] have been proposed and examined in great detail. For example,Pfleiderer et al [23] performed a numerical study on the stretchedfcc phase of hard ellipsoids of revolution and identified the SM2 phase for the range of aspect ratio 3.0$$6.0. Recently,Cohen et al [24] identified the structure of a simple dense fluid of ellipsoids using the direct confocal imaging technique that was dictated by the interplay between rotations and translations caused by the anisotropic shape. Cleaver et al [25] simulated thin lyotropic liquid crystal films for which the anchoring strength and orientation were varied. Teixeira et al [26]carried out a simulation study of the effects of nanoconfinement on a system of hard Gaussian overlap particles. These studies lead to the conclusion [27, 28, 29, 30, 31, 32] that various particle shape and confinement have remarkable effects on the configurations of fluid.
To date,however,systems composed of monodisperse ellipsoids confined between two parallel hard walls have not been fully studied,and hence deserve further investigations [33, 34, 35]. Specifically,when the system is strongly confined in at least one direction,such as confined between two closely apposing walls or in a narrow cylindrical pore [36, 37],the confining perturbations may have significant effects on the structures of the anisotropic systems. Furthermore,hard ellipsoid as a typical anisotropic geometry can represent solid shapes from a hard disk to a needle by simply changing the aspect ratio,and thus the fluid of ellipsoids serves as an ideal model systems of anisotropic rigid particles. In order to deeply understand the effects of the rigid confinements on the various properties of an anisotropic ellipsoid system,we focus on the fluid of ellipsoids confined in one of the most representative confinements, i.e. two apposing flat walls that are readily accessed in experiments.
In this work,we systematically study the fluid of 3D hard prolate ellipsoids with various densities and aspect ratios confined between two impenetrable walls with varying distances using Monte Carlo simulations.
Ⅱ. THEORY AND MODELWe consider a model fluid consisting of $N$ identical hard ellipsoids confined between two parallel impenetrable walls with the normal direction along the $z$ axis located at the positions of $z$=0 and $z$=$L_z$,respectively. Reflective boundary condition is utilized at the $z$ direction for the surfaces of the walls,while periodic boundary conditions are imposed on the $x$ and $y$ directions that have lengths of $L_x$ and $L_y$,respectively. Accordingly,the bulk density is determined as $\rho_{\rm{b}}$=$N/L_xL_yL_z$. The geometry of the ellipsoid is specified by the lengths of three principal semi axes $a$,$b$,and $c$,where $b$ is set to be equal to $c$ and $a$ is chosen as the longest semi axis in our simulations,and thereby the aspect ratio of the ellipsoid is simply characterized by $a/b$. Three unit vectors $\mathbf{l}$,$\mathbf{m}$,and $\mathbf{n}$ along the principal axes are implemented to describe the orientation of an ellipsoid, and hence the vector $\mathbf{r}$ indicating the position of an arbitrary point on the surface of the ellipsoid satisfies
$ \mathbf{r}\cdot\left(\frac{\mathbf{ll}}{a^2}+\frac{\mathbf{mm}}{b^2}+\frac{\mathbf{nn}}{c^2}\right)\cdot\mathbf{r}=1 $  (1) 
For a fluid system consisting of anisotropic particles,there are two interchangeable types of entropy in the simulation system [38],i.e. the translational entropy $E_{\rm{t}}$ which denotes how many states the structure possesses,and the rotational entropy $E_{\rm{r}}$ which is related to the number of distinct states that the particles can be arranged in the same structure. The free energy $F$ in lowdensity systems is dominated by the translational entropy $E_{\rm{t}}$ $(F$$\rightarrow$$E_{t})$,and thus in order to minimize $F$,the translational entropy $E_{\rm{t}}$ should be maximized. In other words,the value of $E_{\rm{t}}$ is much larger than that of $E_{\rm{r}}$ in a low density system,and therefore a change in $E_{\rm{r}}$ can rarely influence $E_{\rm{t}}$. In a marked contrast,a change in $E_{\rm{r}}$ can drastically affect $E_{\rm{t}}$ in a high density system,in order to maximize $E_{\rm{t}}$ the anisotropic particles are aligned along a predominant direction facilitating the transformation of $E_{\rm{r}}$ into $E_{\rm{t}}$. Moreover,the conversion of $E_{\rm{r}}$ into $E_t$ results in the loss of itself,which is proceeded via one of the following two possible mechanism: (i) $L_{\rm{p}}$ is realized by the formation of an ordered nematic phase with the long axis of ellipsoids aligning parallel to the walls,the depletion potential of the walls induces more ellipsoids to be aligned parallelly to the walls [39],(ii) $L_{\rm{v}}$ is proceeded by the formation of another ordered nematic phase with the long axis of ellipsoids aligning vertical to the walls [40]. The depletion potential of ellipsoids that results in the alignment of more ellipsoids to be vertical to the walls. So for the sake of maximizing $E_{\rm{t}}$ we need to compare $L_{\rm{p}}$ with $L_{\rm{v}}$ to determine which mechanism is dominant.
If $L_{\rm{p}}$ is more dominant than $L_{\rm{v}}$,the long axes of ellipsoids prefer to be ordered parallelly to the wall,indicating a stronger directing effect from the depletion potential of the walls. Otherwise the long axes of ellipsoids prefer to be aligned normal to the walls,indicating a stronger directing effect from the depletion potential of the ellipsoids. If $L_{\rm{p}}$ and $L_{\rm{v}}$ are comparable in magnitude,the alignments of the ellipsoids can be either parallel or vertical to the walls.
In this anisotropic system,one interesting quantity is the reduced excess adsorption [41, 42],which is defined as
$ \Gamma=\intop \textrm{d}z\left[\rho^{*}\left(z\right)\rho^{*}_{\rm{b}}\right] $  (2) 
where the absolute value of $\Gamma$ represents the adsorption capability of the walls to the ellipsoids, $z$ is the centerofmass of the ellipsoid, $\rho^{*}(z)$ is the local reduced number density,and $\rho^{*}_{\rm{b}}$ is the reduced number density in the bulk. The reduced number density is defined as the number density multiplied by the volume $V^{*}$ of the cube that circumscribes an ellipsoid,i.e. $\rho^{*}$=$\rho V^{*}$. $V^{*}$ is set to be the unit of volume in the simulations to guarantee that the unit of length in every simulation is adjusted with the aspect ratio to conserve $V^{*}$. The local orientation order parameter $S$,defined as the average of the largest eigenvalue of the $Q$ tensor [43] is
$ Q_{\alpha\beta}=\frac{1}{N}\sum_{i=1}^{N}\left(\frac{3}{2}u_{\alpha}^{i}u_{\beta}^{i}\frac{1}{2}\delta_{\alpha\beta}\right) $  (3) 
which is another critical quantity in the anisotropic system that quantifies the orientational order. Here $u_{\alpha}^{i}\left(\alpha=x,y,z\right)$ is the Cartesian component of the longest axial unit vector of ellipsoid $i$. The value of the order parameter is close to zero for an isotropic phase,or close to $1$ for a nematic phase in a finite system [35].
The $NVT$ Monte Carlo simulation is performed to obtain the equilibrium properties [44]. Hardcore interactions are implemented between ellipsoids as well as between ellipsoids and the walls,i.e. the interaction is zero when there is no any overlap between ellipsoids,or between the ellipsoids and the walls,otherwise the interaction is set to be infinite.
In the simulations,different sampling algorithms introduce errors to different degrees. In molecular dynamics simulations,a basic task is to estimate whether the results can be validated or not [45, 46]. Instead,in Monte carlo simulations [47],an importance$$sampling algorithm introduced by Metropolis et al is often adopted,in which every accessible point in configurational space is considered to be reached in a finite number of Monte Carlo steps from any given point. Moreover,many efficient Monte Carlo simulations have not been proven to be ergodic [44]. In order to lower the systematic errors in the biased sampling algorithms, we carry out the simulations with a large number of ellipsoidal particles,$N$=1000. All results are obtained for prolate ellipsoids with aspect ratios $a/b$ of 1.0$$3.5. Each simulation consists of a succession of trial moves where the average acceptance of translational/orientational moves is about 25%$$30%. Each trial move of a chosen ellipsoid consists of a small random translational displacement and a small angle of random rotation. If the trial move results in any overlap between the moved ellipsoid with others,it is rejected,and otherwise,it is accepted. It takes $10^{6}$ MC steps per particle for equilibration and 2$\times$10$^{6}$ MC steps for data collection.
Ⅲ. RESULTS AND DISCUSSIONFigure 1 shows the local density profiles of the ellipsoids with respect to the distance $z$ to the wall at $z$=0 for four different aspect ratios $a/b$=1.5,2.4,2.7,and 3.0,and a fixed bulk density of $\rho_\textrm{b}$=0.3. There are two obvious features in the four density profiles. The first one is that the density vanishes when approaching $z$=$b$,which is dictated by the excludedvolume effect of the hard ellipsoid. The other is that there is a remarkable peak between $z$=$b$ and $z$=$a$ in each density profile,but the peak width and height are sensitive to the aspect ratio,i.e. that the peak is widened while the peak height is lowered as the aspect ratio increases. The presence of the peak indicates the layering arrangement of particles induced by the adsorption of the walls,and the variation tendency is dictated by the orientational configurations of the ellipsoids from parallel to normal to the walls that weakening the layering effect. Accordingly,the magnitude of oscillation is also reduced by increasing the anisotropy of ellipsoids. For this low density system of $\rho_\textrm{b}$=0.3,the bulk state is an isotropic fluid. Limited local densities when $z$ is close to $b$ for all four aspect ratios suggest that few ellipsoids are aligned with their long axes parallel to the walls.
Based on these density distributions,the surface excess adsorption is calculated for various aspect ratios and plotted in Fig. 2. We find that the reduced excess adsorption $\Gamma$ increases (or its magnitude decreases) with the aspect ratio in Fig. 2,and we ascribe this to the orientational ordering effect of the two walls to the ellipsoids. Specifically, when the aspect ratio $a/b$=1.0,which leads to the isotropic hard sphere system,our result is similar to that reported by Roth et al [34]. In general,the reduced surface excess adsorption is rather weak at the considered range of aspect ratios in Fig. 2,which simply indicates that the effect of the wall confinement is weak in the anisotropic fluid of ellipsoids at the low density. The weak excess adsorption is attributed to the factor that the translational entropy $E_{\rm{t}}$ is expected to be more dominant over the rotational entropy $E_{\rm{r}}$ at $\rho_{\textrm{b}}$=0.3 [48],and thereby the ellipsoids prefer to be more uniformly distributed in the entire available space.
The excess adsorption curves with respect to the aspect ratio of ellipsoids for two higher bulk density of $\rho_{\rm{b}}$=0.4 and $\rho_{\rm{b}}$=0.6 are present in Fig. 3. For $\rho_{\rm{b}}$=0.4,similar to that of $\rho_{\rm{b}}$=0.3,the excess adsorption monotonically increases as the aspect ratio is increased. Surprisingly,when the bulk density is increased to $\rho_{\rm{b}}$=0.6,the excess adsorption does exhibit a monotonic behavior as the aspect ratio varies from that at two lower bulk densities of $\rho_{\rm{b}}$=0.3 and $\rho_{\rm{b}}$=0.4. Indeed,it increases at the range of small aspect ratios until it reaches a maximum point at $a/b$$\approx$2.9,and then it abruptly drops off. At the aspect ratio of $a/b$=3.0,the magnitude of the excess adsorption $\Gamma$ has dropped down to $\Gamma$$\approx$$$0.076.
In order to understand this underlying mechanism in the nonmonotonic variation behavior of the excess adsorption with respect to the aspect ratio,we probe into the internal configurations of the fluid at the high bulk density of $\rho_{\rm{b}}$=0.6 for various aspect ratios,i.e. the inhomogeneous density profiles induced by the perturbation of the walls. Figure 4 shows four typical density distributions obtained for distinct aspect ratios located around the maximal point of the excess adsorption,$a/b$=2.7,2.8,2.9,and 3.0. Generally,each density distribution exhibits two notable peaks. One peak is located at the position very close to $z$=$b$,which is contributed by the ellipsoids aligned parallelly to the wall. The other peak is located at the position slightly further than $z$=$a$,and hence it is resulted in by the aggregation of the ellipsoids with the along axes normal to the wall. Importantly,when the aspect ratio increases from $a/b$=2.7 smaller than the critical value of 2.9 to be 3.0, the magnitudes of the two peaks change drastically. On the one side,the two peaks are maintained until they merges into a unique peak in the limit case of isotropic hard spheres as the aspect ratio is decreased. On the other side, the first peak is reduced indicating that less ellipsoids are aligned with the long axes parallel to the wall,whereas the second peak is strengthened indicating that more ellipsoids with the long axes normal to the wall are aggregated to the surface with a distance of around the length of the semimajor axis to the wall. In particular,for $a/b$=2.8,the first peak is only slightly weaker than the second one,while for $a/b$=2.9,the first peak has changed to be significantly weaker than the second one. For further increased $a/b$=3.0,the first peak has become very weak compared with the second peak. The varying tendency of the magnitudes of the two peaks when the aspect ratio is increased from $a/b$=2.7 going through the critical value of 2.9 to 3.0 leads to an obvious conclusion that the aggregation of ellipsoids near the wall is predominantly governed by the ellipsoids aligned with the long axes normal to the walls when the aspect ratio exceeds 2.9 in the highdensity ellipsoid fluid confined between two walls.
In order to understand the variation of the peaks further, we estimate the probability of the ellipsoids aligned in the angle interval $P$($\theta$$ < $$\theta_{0}$),where $\theta$ is the angle between the long axis of the ellipsoid and the $z$ direction,i.e.,$\theta$=0 means that the long axis of the ellipsoid is along the $z$ axis while $\theta$=$\pi/2$ means that the ellipsoid is aligned in the plane of the wall. Here we choose $\theta_{0}$=$\pi/5$ as a specific value to represent the case of small angles,and the data are collected in the slit area with the distance to the wall from $b$ to $1.4a$. The probability of $P$($\theta$$ < $$\theta_{0}$), which measures the aligning degree of the ellipsoids with the long axes normal to the walls,is plotted in Fig. 5. We find that the probability increases slowly as the aspect ratio is till $\alpha$=2.9,and then the probability increases rapidly. This indicates that more and more ellipsoids tend to be aligned perpendicularly to the wall when the aspect ratio exceeds 2.9,which results in the abrupt decrease (or increase in the magnitude) of the excess adsorption observed in Fig. 3.
At a dense density $\rho_{\rm{b}}$=0.6,the rotational entropy $E_{\rm{r}}$ is not negligible [34],and can be transformed into $E_{\rm{t}}$. Thus,the directed ordering of the ellipsoids near the walls is dictated by the delicate balance between $L_{\rm{p}}$ and $L_{\rm{v}}$. For aspect ratios smaller than $2.9$,Fig. 4 implies that the competition between $L_{\rm{p}}$ and $L_{\rm{v}}$ is too fierce and $L_{\rm{p}}$ is comparable to $L_{\rm{p}}$,and thereby the ellipsoids are aligned either parallelly or vertically to the walls.
However,for aspect ratios larger than 2.9,there exists a single dominant peak because $L_{\rm{v}}$ is much larger than $L_{\rm{p}}$,indicating that the vertical configurations are dominant. In fact,Fig. 4 also suggests that the depletion potential between the ellipsoids is more pronounced that the depletion potential between the ellipsoids and the walls,and as a result the ellipsoids prefer to stand on the walls.
Inspired by the variation of the ordering configurations near two walls sensitive to the aspect ratio,we examine the order parameter $S$ of the hard ellipsoids within the wall confinement for two aspect ratios $a/b$=2.8 and 3.0,which are smaller and larger than the critical value of 2.9,respectively. The results of orientational order parameters with respect to the bulk density in Fig. 6 reveal that the variation of orientational order parameters are also sensitive to the aspect ratio near the critical value of $a/b$=2.9. For $a/b$=2.8 smaller than the critical value,the orientational order parameter maintains roughly a constant low value of $0.1$ indicating a near isotropic state. However,for $a/b$=3.0 larger than the critical value,$S$ increases almost linearly as the bulk density increasing. For example,for $\rho_{\rm{b}}$=0.7,$S$ is increased to be 0.55,which indicates an isotropicnematic transition near the walls. This observation supports the results further in Fig. 4. For instance,for $a/b$=2.8 at $\rho_{\rm{b}}$=0.7,the low orientational order suggests the the ellipsoids do not have a preferred alignment between parallel and vertical to the walls,indicating that $L_v$ and $L_p$ are comparable. Whereas,when the aspect ratio becomes $a/b$=3.0,the ellipsoids are directed to be ordered by the walls,and accordingly,$L_v$ becomes more dominant than $L_p$.
In order to distinguish the spontaneous ordering in the bulk anisotropic system from the induced ordering by the wall confinement,we compare their orientational order parameters of varying bulk densities for a specific aspect ratio of $a/b$=3.0. In Fig. 7,we can see that $S$ in the bulk systems only slightly increases as the bulk density increases and sustains a small value indicating that the bulk system of the ellipsoids with $a/b$=3.0 is truly isotropic. Furthermore, it has been known that the transition from isotropic to nematic in the bulk system occurs above the bulk density of $\rho_{\rm{b}}$=0.7 [49, 50, 51, 52, 53]. Therefore we are able to conclude that the anisotropic ordering behavior in the confined system arises from the confinement effect of the walls.
Ⅳ. CONCLUSIONWe study the structure of a hard ellipsoid fluid confined between two hard walls by using the $NVT$ Monte Carlo simulation. Our results suggest that the reduced surface excess adsorption is a sensitive parameter to identify the transition between an isotropic state and a nematic phase. Specifically,for the high bulk density $\rho_{\rm{b}}$=0.6,but below the bulk transition density from isotropic to nematic, the reduced excess adsorption surprisingly exhibits a nonmonotonic change as the aspect ratio varies,i.e. that it increases (or decreases in magnitude) at small aspect ratios until reaches a maximum at a critical aspect ratio of $a/b$=2.9,and then decreases abruptly. This finding is confirmed by the analysis of different relevant quantities,the density configurations and the orientational order parameter for various aspect ratios. Based on these consistent results we conclude that when the aspect ratio $a/b$$\geq$2.9 the confinement effect of the walls induces an ellipsoid fluid to be ordered with the the long axes of ellipsoids aligned normally to the walls,while the walls have minor directing effect on the ordering of the ellipsoids when $a/b$$ < $2.9. This will place an useful,tight theoretical constraint on the investigation of an isotropicnematic phase transition in the future.
Ⅴ. ACKNOWLEDGEMENTSThis work is supported by the National Natural Science Foundation of China (No.10874111,No.11304169,and No.11174196).
[1]  D. E. Sullivan and G. Stell, J. Chem. Phys. 69, 5450 (1978). 
[2]  J. C. Neu, Phys. Rev. Lett. 82, 1072 (1999). 
[3]  C. N. Patra and S. K. Ghost, J. Chem. Phys. 106, 2752 (1997). 
[4]  R. D. Groot, N. M. Faber, and J. P. van der Eerdenn, Mol. Phys. 62, 861 (1987). 
[5]  V. M. Bedanov and F. M. Peeters, Phys. Rev. B 49, 2667 (1994). 
[6]  V. A. Schweigert and F. M. Peeters, Phys. Rev. B 51, 7700 (1995). 
[7]  R. Zangi, J. Phys.: Condens. Matter 16, S5371 (2004). 
[8]  B. Götzelmann, A. Haase, and S. Dietrich, Phys. Rev. E. 53, 3456 (1996). 
[9]  R. Roth, R. Evans, A. Lang, and G. Kahl, J. Phys.: Condens. Matter 14, 12063 (2002). 
[10]  Y. Rosenfeld, Phys. Rev. Lett. 63, 980 (1989). 
[11]  B. H. Lin, J. Yu, and S. A. Rice, Phys. Rev. E 62, 3909 (2000). 
[12]  H. Reiss, H. L. Frisch, E. Helfand, and J. L. Lebowitz, J. Chem. Phys. 32, 119 (1960). 
[13]  P. Tarazona, Phys. Rev. A 31, 2672 (1985). 
[14]  A. Donev, J. Burton, F. H. Stillinger, and S. Torquato, Phys. Rev. B 73, 054109 (2006). 
[15]  T. Schilling, S. Pronk, B. M. Mulder, and D. Frenkel, Phys. Rev. E 71, 036138 (2005). 
[16]  R. Blaak, B. M. Mulder, and D. Frenkel, J. Chem. Phys. 120, 5486 (2004). 
[17]  Y. Li, H. Miao, H. R. Ma, and J. Z. Y. Chen, Soft Matter 9, 11461 (2013). 
[18]  C. Avendano and F. A. Escobedo, Soft Matter 8, 4675 (2012). 
[19]  K. Wojciechowski and D. Frenkel, Comp. Met. Sci. Technol. 10, 235 (2004). 
[20]  M. A. Bates and D. Frenkel, J. Chem. Phys. 112, 10034 (2000). 
[21]  P. G. Bolhuis and D. Frenkel, J. Chem. Phys. 106, 666 (1997). 
[22]  A. G. Moreira and R. R. Netz, Euro. Phys. J. E 8, 33 (2002). 
[23]  P. P eiderer and T. Schilling, Phys. Rev. E 75, 020402(R) (2007). 
[24]  A. P. Cohen, E. Janai, E. Mogilko, A. B. Scho eld, and E. Sloutskin, Phys. Rev. Lett. 107, 238301 (2011). 
[25]  D. J. Cleaver and P. I. C. Teixeira, Chem. Phys. Lett. 338, 1 (2001). 
[26]  P. I. C. Teixeira, F. Barmes, F. C. AnquetilDeck, and D. J. Cleaver, Phys. Rev. E 79, 011709 (2009). 
[27]  F. Barmes and D. J. Cleaver, Phys. Rev. E 71, 021705 (2005). 
[28]  F. Barmes and D. J. Cleaver, Phys. Rev. E 69, 061705 (2004). 
[29]  M. Moradi, R. J. Wheatley, and A. Avazpour, Phys. Rev. E 72, 061706 (2005). 
[30]  P. I. C. Teixeira, A. Chrzanowska, G. D. Wall, and D. J. Cleaver, Mol. Phys. 99, 889 (2001). 
[31]  G. W. Wall and D. J. Cleaver, Mol. Phys. 101, 1105 (2003). 
[32]  G. W. Wall and D. J. Cleaver, Phys. Rev. E 56, 4306 (1997). 
[33]  W. H. Li and H. R. Ma, J. Chem. Phys. 119, 585 (2003). 
[34]  D. Frenkel, B. M. Mulder, and J. P. McTague, Phys. Rev. Lett. 52, 287 (1984). 
[35]  J. Talbot, D. Kivelson, M. P. Allen, G. T. Evans, and D. Frenkel, J. Chem. Phys. 92, 3048 (1990). 
[36]  P. A. Kralchevsky and K. Nagayama, Adv. Colloid Interface Sci. 85, 145 (2000). 
[37]  M. Breuer, J. Bernsdorf, T. Zeiser, and F. Durst, Internat. J. Heta Fluid 21, 186 (2000). 
[38]  M. S. Searle and D. H. Williams, J. Am. Chem. Soc. 114, 10690 (1992). 
[39]  H. Miao, Y. Li, and H. R. Ma, J. Chem. Phys. 140, 154904 (2014). 
[40]  B. Götzelmann, R. Roth, S. Dietrich, M. Dijkstra, and R. Evans, Europhys. Lett. 47, 398 (1999). 
[41]  R. Roth and S. Dietrich, Phys. Rev. E 62, 6926 (2000). 
[42]  J. R. Henderson and F. van Swol, Mol. Phys. 51, 991 (1984). 
[43]  U. Fabbri and C. Zannoni, Mol. Phys. 58, 763 (1996). 
[44]  D. Frenkel and B. Smit, Understanding Molecular Simulation, New York: Academic Press, A Division of Harcourt Inc., (2002). 
[45]  W. F. van Gunsteren and A. E. Mark, J. Chem. Phys. 108, 6109 (1998). 
[46]  W. F. van Gunsteren, J. Dolenc, and A. E. Mark, Curr. Opin. Struct. Biol. 18, 149 (2008). 
[47]  X. G. Hong and K. Tamura, Chin. Phys. Lett. 20, 1315 (2003). 
[48]  D. Frenkel and B. M. Mulder, Mol. Phys. 55, 1171 (1985). 
[49]  X. Zheng, W. Iglesias, and P. Pal yMuhoray, Phys. Rev. E 79, 057702 (2009). 
[50]  M. Baus, J. L. Colot, X. G. Wu, and H. Xu, Phys. Rev. Lett. 59, 2184 (1987). 
[51]  R. Holyst and A. Poniewierski, Mol. Phys. 68, 381 (1989). 
[52]  U. P. Singh and Y. Singh, Phys. Rev. A 33, 2725 (1986). 
[53]  J. F. Marko, Phys. Rev. Lett. 60, 325 (1988). 
b. 上海交通大学机械与动力工程学院, 上海 200240