- Email: [email protected]
Intrinsic frame transport for a model of nematic liquid crystal
PHYSICA ELSEVIER
Physica A 240 (1997) 173-187
Intrinsic frame transport for a model of nematic liquid crystal S. C o z z i n i a'*, L . F . R u l l a, G . C i c c o t t i b, G . V . P a o l i n i c'1 a Dpto. Fisica At6miea Molecular y Nuclear Aptdo. 1065, Universidad de Sevilla, 41080 Sevilla, Spain b Dipartimento di Fisica, Universitgt La Sapienza Roma, Piazzale Aldo Moro 2. 00185 Roma, Italy e Laboratoire des Solides Irradibs, Ecole Polytechnique, 91128 Palaiseau Cedex, France
Abstract
We present a computer simulation study of the dynamical properties of a nematic liquid crystal model. The diffusional motion of the nematic director is taken into account in our calculations in order to give a proper estimate of the transport coefficients. Differently from other groups we do not attempt to stabilize the director through rigid constraints or applied external fields. We instead define an intrinsic frame which moves along with the director at each step of the simulation. The transport coefficients computed in the intrinsic frame are then compared against the ones calculated in the fixed laboratory frame, to show the inadequacy of the latter for systems with less than 500 molecules. Using this general scheme on the Gay-Berne liquid crystal model, we evidence the natural motion of the director and attempt to quantify its intrinsic time scale and size dependence. Through extended simulations of systems of different size we calculate the diffusion and viscosity coefficients of this model and compare our results with values previously obtained with fixed director.
1. I n t r o d u c t i o n
In the nematic phase o f a liquid crystal all the molecules are aligned, in average, along a common direction, identified by a unit vector n called director. As a consequence o f the broken symmetry along n transport properties in the nematic phase are considerably more complicated than in the isotropic one. While in the isotropic phase the second rank diffusion and heat conductivity tensors have only one independent component, this number doubles in the nematic phase. Cross coupling terms which are forbidden in the isotropic phase raise the number o f independent components o f the fourth rank viscosity tensor from three to seven [1]. * Corresponding author. E-mail: [email protected].
I Present address: Department of Chemistry, University of Liverpool, Liverpool L69 3BX, U.K. 0378-4371/97/$17.00 Copyright (~) 1997 Elsevier Science B.V. All rights reserved P l l S0378-4371 (97)001 39-8
S. Cozzini et al./Physica A 240 (1997) 173-187
174
The transport coefficients of a nematic system model can be computed by standard Equilibrium Molecular Dynamics (EMD) by evaluating the appropriate Green-Kubo (GK) integrals along the simulation trajectory [2]. GK relations for nematic liquid crystals are straightforward in the case of the diffusion and heat conductivity coefficients, while the derivation of formulae for the viscosity coefficients requires some care. GK relations for viscosity coefficients were derived more than 20 years ago by Forster [3] using projection operator techniques; a more recent formulation was proposed by Sarman and Evans [4] by means of linear response theory. However, it must be noted that the two formulations are not completely equivalent and a satisfactory full formalization of the GK expressions for all the viscosity coefficients is still lacking. From a theoretical point of view, all GK relations are expressed in a coordinate system in which one of the axes is aligned with the director, whose orientation is considered constant in time. This assumption raises a delicate problem in a computer experiment. In the absence of an extemal field the director orientation is a quasiconserved degree of freedom, and is expected to diffuse slowly on the unit sphere. This diffusion can be considerable in samples of typical simulation size (few hundred molecules) and therefore needs to be taken into account. Indeed, neglecting this effect can lead to large errors in the estimation of the transport coefficients. So far, this problem has been circumvented by fixing the director with appropriate constraints [4-6] or by applying external fields to stabilize its orientation [7]. The main drawback of these approaches is that they tend to bias the dynamics of the system. In this work, we propose a new approach which leaves the system unperturbed, by performing our measurements in the director's intrinsic reference frame. This frame is defined at each instant by the position of n. All the physical quantities of interest are calculated with respect to this moving frame. Our investigation is twofold. First, we study the size dependence of the director motion and check for an effective separation in time scale between the director's diffusion and molecular relaxation. Second, we test our approach by calculating the diffusion and the viscosity coefficients for a model system (the Gay-Berne model) and comparing our results with previous studies of similar models. The paper is organized as follows: in Section 2 we introduce the relevant theoretical formulae and summarize the GK relations employed. In Section 3 we describe our computer experiment, while in Section 4 the results are presented and discussed. Finally in Section 5 we draw some conclusions.
2. Formalism and kinematics of the director
The degree of molecular alignment in a liquid crystal is measured by the order parameter S, which is defined as the largest eigenvalue of the second rank-order tensor [8]:
Q=
E 1
UiUi
-
-
1 ,
(1)
S. Cozzini et al./Physica A 240 (1997) 173-187
175
where 1 is the unit tensor, N is the number o f molecules and ui is the unit vector parallel to the axis o f molecule i. The eigenvector corresponding to S is defined as n, the director. The angular velocity o f the director is f2 - n x ft. This definition rests on the fact that for uniaxial molecules the linear velocity is orthogonal to n. To define the intrinsic system associated with the moving director we observe that, by definition, ( Q - S 1 ) - n = 0 . This means that every row or column of the Q - S1 matrix lies in a plane perpendicular to the director. Therefore, choosing the first column o f this matrix we obtain a vector m orthogonal to n. The vector product between m and n defines the remaining axis. It is worth noting that this formulation of the intrinsic frame guarantees its continuity in time as a consequence o f the continuity of Q. In a liquid crystal there are two independent components o f the second rank diffusion tensor D~I3: DIIII, which relates forces and fluxes parallel to the director, and D x ± relating forces and fluxes in the plane perpendicular to the director. It is easy to derive the Green-Kubo relations for these coefficients: Ollll ---
/0
dt (vll(t)Vll(0)),
D±± =
F
dt (Vx(t)v±(O)),
(2)
where vii and Vx are the parallel and perpendicular component of the linear velocity of each molecule with respect to n. Here, for simplicity, we drop out the average over all the molecules o f the sample. The viscosity coefficients can be defined using the phenomenological relations existing between the pressure tensor and the strain rate. To derive them, we note that the entropy production a for a nematic liquid crystal subject to a shear stress can be expressed, using the same notation o f Ref. [5], as a= -l[ps
: (Vn)S + p a . (~7 x u - 2f2) + ½(tr(P) - Peq)V" n],
(3)
where T is the temperature, P is the pressure tensor, u is the velocity field and f2 the angular velocity of the director. The notation A s = l [ A + A r] - 1 tr(A) stands for the symmetric traceless part o f a second rank tensor while A a = ½[A - A r] is its ant±symmetric part. In the above expression the thermodynamical fluxes can be identified as ps, pa and tr(P) - Peq, conjugated, respectively, with the three forces ( V u ) s, V x u - 212 and V . u . The condition of uniaxal symmetry yields the following phenomenological relations between fluxes and forces [4,6,9]: ps
=
_2~/(Vu)S _ 2~ 1[(nn)S. (Vu)~]~ _ 2~3(nn)~[(nn)S : (Vu)~] _ ((nn)SX7. + 2f/2[(nn) ~. ~. (½V × u - f2)] ~ ,
u
(4)
2P" = - ~l [32-1 - (nn)S] - (½V × u - f2) - ~2e : [(nn) s. (Vu)S],
(5)
½[tr(P) - Peq] -- -r/v x~7"u - K(nn)~ "(Vu) ~ .
(6)
Here e is the third-order Levi Civita tensor. The bulk viscosity ~z, the twist viscosity 71 and the shear viscosity ~/ must be positive because they relate tensorial quantities o f the same order and would be the only non-zero terms in an isotropic system. The anisotropic contribution is given by the two coefficients ~1 and 43 which relate the same
S. Cozzini et aL /Physica A 240 (1997) 173-187
176
quantities like r/ but in this case expressed using the order tensor nn [9] to take into account the lower symmetry. All the other terms are cross couplings between tensors of different parity and ranks and can have either sign. The coefficient ~ relates the symmetric traceless part of the pressure tensor to the divergence of the velocity field while q2 relates the same quantity to the angular velocity of the director relative to the rotating fluid. These two coefficients are linked to the coefficients x and 72 by the following Onsager relations: ~2 z ~72/2 and ~ =- K. The number of independent variables is therefore seven. GK expressions for five of these coefficients are derived in Ref. [4] and we will use that formulation throughout this work. In a director-based coordinate system with the third axis along the director, the GK expressions for the five viscosity coefficients are given by
2tlpt2tlQ__.._p r/= ½(2q3131 + 111212)
3~bno '
(7)
2 tlpI2tlK2P
~1 =2(t/3j31 -
q1212) -
700
(8) ,
7]3 -- 3(r13333 .2_ 2r/ill 1 ) _ 2173131 _ r/1212 +
2
-
2 tlPOrlOP ~
,
(10)
~boo ' ~bo~ .;4
(9)
( 11 ) '
where the r//jkt,qSno, q~o and q~Pa are all time-correlation fimction integrals (TCFI) and are given by the following relations: tlijkI = flV c~O0=flV
tl~a = flV
(Pi;(S)PkSl(O))eqds,
/o
(Q2(s)Q2(O))eqds=flV
(P~j(s)fa2(O))eq ds = flV
(12)
io /o
(Ql(S)Ql(O)}eqds,
(13)
(p2s3(s)f~l (0))e0 ds,
(14)
tl~pa=flV fo ~ (p~3(s)f21(O))~qds=flV fo ~c (P~ll(S)Q2(O))eqds= ~1 .
(15)
Here P,] ( ~ ) are elements of the symmetric traceless (antisymmetric) pressure tensor, ~ is the component of the angular velocity of the director, fl = 1/ksT and V is the volume of the simulation cell. Some symmetry relations hold among different r/0kt [4]. We note that the TCFI in Eq. (15) is not an experimental number but a sum rule [4]. Finally, a word of caution must be spent. The above formulae are derived under the assumption that the director remains fixed during the whole simulation. In that sense
177
S. Cozzini et aLIPhysica A 240 (1997) 173 187
a director-based coordinate system is an inertial frame, a condition not necessarily met during simulation. This effectively poses a problem only if the time scale over which the displacement of the director is appreciable is comparable with the decay time of the correlation functions we want to measure. We will discuss this issue further in the light of the results presented in Section 4.
3. Model and computer experiment In this section we briefly discuss the Gay Berne (GB) model used and the simulations performed. The GB model [10] has been widely used as a semi-realistic model for liquid crystals in recent years. It can be thought of as a fluid of pointqike particles, each carrying a unit vector u which mimics the long molecular axis. The particles interact via an anisotropic Lennard-Jones potential, which depends on the relative orientation and position of a pair of molecules. The explicit interaction energy of two molecules i and j is given by
Uij = 4g(t ij, hi, uj)
E(
(70
rij -
o(~g, ui, uj) + (70
12
(7O
r6 - cr(~ii, ui, uj) + (70
' (16)
where r~ is the relative distance between the centers of mass of the molecules i and j and bj = ri;/ro is a unit vector along the direction of the intramolecular distance: r U- - r i - rj. The range parameter a(~,y,u/,nj) is expressed as
{ a(~ij, ui, u/)=6o
Z [(ro"ni%F-ij'Uj)2 1-~
[
l+)~(ui.uj)
+
(Fij.ui--ro'nj)2J } -1':2 l - z(ui.uj)
(17)
where Z = ( K2 - 1)/(K2 + 1) and K is the length to breadth ratio for a prolate ellipsoid of revolution. The strength of the interaction e, is a product of two terms:
(18)
g(r.ij, ui, nj)--~ F,O[gl(ni, nj)] v )< [~2(Ftj, ni, nj)] 'u
with < (u~, uj) = [1 -
z2(ui • uj):]-1/2
~:2(~0,ui, u j ) =
1- ~
(19)
and
1 - z'(ni'u/)
J
Here •' = [(K') 1/" - 1]/[(~') l/" + 1] and s:' is the ratio of the potential energy minima of the side-by-side and the end-to-end configuration. The phase diagram of this fluid has been extensively studied by de Miguel et al. [11,12].
S. Cozzini et aL /Physica A 240 (1997) 173-187
178
Recently, some transport coefficients were measured for the GB-models by different methods. Sarman and Evans used Equilibrium and Non-Equilibrium Molecular Dynamics simulations to measure diffusion coefficients, heat fluxes [13] and viscosities [4] for a purely repulsive GB system. De Miguel et al. [14] measured the diffusion coefficients by EMD while Smondyrev et al. performed EMD calculations [14] of viscosity coefficients. We therefore chose this model for our computer experiment in order to check our results against a relatively wide spectrum of data obtained with different techniques. Table 1 List of the runs. Here T is the averaged temperature and the S the averaged order parameter Name
Size
Run length (~)
T
S
T PL B LC XX YY
178 178 256 256 500 500
3250 4000 4000 4000 3000 3000
1.242 ± 0.009 1.260 ± 0.010 1.186 4- 0.009 1.244 4- 0.010 1.259 i 0.006 1.260 ± 0.006
0.77 4- 0.02 0.77 4- 0.02 0.774-0.02 0.76 4- 0.02 0.75 4- 0.02 0.75 ± 0.02
i
178
0
................................................................
-1
-2
-3
(a)
I
i
I
I
-4
-2
2
4
Fig. 1. The director motion for systems of different size. See text for explanation.
- -
P I
3 I
'(p~nupuo3) 1 ':~!2, (o) 0 ]
3-
'¢8-
I
t-
00£ i
i
i
£
i
(q) #
3
I
0
I
17-
3l
{
g-
3"
993 I
6LI
ZgI-EZI (Z661) OPE V va!g,fVd/"l n la lUlZZOD "~
180
S. Cozzini et al./Physica A 240 (1997) 173-187
MEAN SQUARED 0.008
J
0.006
- ....
DISPLACEMENT i
OF THE DIRECTOR F
N=178 N=256
r
j
"
o
"E 0.004 I
¢-
0.002
0.000
0.0
2000
4000
6000
8000
10000
time ( in timesteps) Fig. 2. The director MSD for systems of different size.
In all our simulations we use the standard anisotropy parameter Z = 3 and Z~-- 5 and exponents p = 2 and v = 1. The simulation results are given in length, energy, mass and time units of a0, e0, m and z--ao(m/eo) 1/2. Three systems of different sizes were investigated: 178, 256 and 500 molecules, all at the same reduced density (p* --0.345) and the same starting reduced temperature T * = 1.25. The stability of the nematic phase at this state point is well known [ll,12]. For the largest systems (256 and 500 molecules) the cutoff distance beyond which the potential was set equal to zero was 4.0a0. For the smallest system (N = 178 molecules) the correspondingly smaller box size dictated a slightly shorter spherical cutoff distance r~* = 3.8a0. Consequently the smallest system does not correspond exactly to the larger ones, but the nematic phase is still stable and the thermodynamical properties of the three systems are virtually the same. We performed constant energy molecular dynamics simulations with the usual periodic boundary conditions. The integration algorithm used is a standard velocity Verlet for the translational degrees of freedom while for the rotational degrees we adapted the standard algorithm for uniaxial molecules [15] to the velocity Verlet scheme. The time step was 0.001z and reduced moment of inertia was set to unity. For each size we performed two separate runs starting from independently equilibrated
S. Cozzini et al./Pkysica A 240 (1997) 173 187
181
Diffusion at short times 0.050
i
i
a--od
par(int)
~-----~ d p a r (fix) *-----~ d p a r (fix) o
--~ d p a r ( i n t )
0.040
¢-1 "O ~) ¢.) -1
0.030
0.020
I
0.0100. 0
50.0
I
100.0 150.0 time ( reduced units)
I
I
200.0
250.0
Fig. 3. The diffusion coefficients measured over a short run for N = 178. Here Dper and Dp,r indicate, respectively, D ± ± and DII!I, while int and fix means intrinsic and fixed frame. See text for discussion.
configurations with different director orientation. Table 1 summarizes the simulations performed. The length of the runs is necessarily large in the attempt to get good convergence of the results. Along each trajectory we calculated the TCFI over a correlation time of 5~. The correlations were also computed in a director-based coordinate system fixed throughout the simulation and defined by the position of the director at the beginning of the run. In order to monitor the director's dynamics its orientation was stored every ten time steps (0.01~) for each run.
4. Results The first aspect we are interested in is the dependence of the motion of the director on the size of the system. A representation of this motion is obtained by describing the position of n(t) via the angles 0 and 4) formed with respect to the fixed frame associated to n(t = 0). In Fig. 1 the quantities tg(0)cos(~b) and tg(0)sin(~b) are plotted at each recorded time step. The starting point of the trajectory is by definition at the origin. The director's drift can be better appreciated because of the factor tg(0).
S. Cozzini et al./Physica A 240 (1997) 173-187
182
For instance, an azimuthal angle 0 - - 9 0 ° corresponds to infinite displacement. Fig. 1 clearly indicates a size dependence of the motion of the director. Figs. l(a) and (b) correspond to N = 178 and N = 2 5 6 , respectively, and Fig. l(c) refers to N - - 5 0 0 . More quantitative information is contained in the mean square displacement (MSD) of the director, reported in Fig. 2, where again the size dependence is apparent. We are currently running some tests with larger systems in order to see if a general relationship between the director's diffusion and the system size can be inferred. Let us discuss now the results found for the diffusion coefficients. Fig. 3 shows the values of the two diffusion coefficients DII II and D ± ± calculated in the two different reference frames. The length of the trajectory over which the TCF1 is calculated is reported in the abscissas. By increasing the statistics we expect the average to stabilize around the final value. This is indeed what we observe for the two diffusion coefficients FlintII is twice as large as Flint ~ x x in agreement measured in the intrinsic frame. In this case ~11 with previous results [14,13]. The behaviour of the coefficients measured in the fixed frame is however different. At short times the diffusion coefficients measured in the two frames coincide within the error. By increasing the statistics DilfixII and D ~ eventually deviate from their initial values. This result is a clear signature of the motion of the
N=178
0.050
i
i
d_par z...---~ d p a r *---.--~ d _ p e r o o d_per
(Int) (fix) (fix) (Int)
0.040
¢.-
o "1o
0.030
0.020
n
0'01500.0 (a)
i
L
1500.0
2500.0
3500.0
time ( reduced units)
Fig. 4. The diffusion coefficientscalculated along the whole runs for systemsof different size: (a) 178 atoms, (b) 256, (c) 500. Labels are as in previous figure.
S. Cozzini et al./Physica A 240 (1997) 173-187
183
N=256 0.050
i
i
o o = ...... *------~ ~
d__par (Int) d _ p a r (fix) d _ p e r (fix) d _ p e r (Int)
0.040
¢9 t"0
0.030
¢O
0.020
0.010 0,0
' 1000.0
2000.0
3000.0
time ( reduced units)
(b)
N=500 0.050
i
o
o d _ p a r (int)
~*----~ d _ p a r (fix) c--
•.-> d _ p e r (fix) ~ d _ p e r (Int)
0.040
ro
0.030
o
0.020
0.010 0.0
(c)
r 1000.0
,
2000.0 time ( reduced units)
Fig. 4 (continued).
3000.0
184
S. Cozzini et al./Physica A 240 (1997) 173 187
director. After about 100r the initial director orientation with respect to which we are measuring our diffusion coefficients (fixed frame) is no longer the symmetry axis of our system. Fig. 4 shows the same quantities measured for longer runs of about 3000r for the three systems studied. The behaviour exhibited over short trajectories is here more pronounced due to the longer times. In all the three cases the diffusion coefficients drift away from their initial values, the more rapidly the smaller the size of the system. This behaviour is not surprising and is due to the fact that the direction along which we are measuring our coefficients is losing correlation with the symmetry axis of the system. It could be expected that the diffusion coefficients measured parallel and perpendicularly to this arbitrary direction would eventually tend to roughly the same value. Indeed this seems the case for N = 500. For the smaller samples the values of the coefficients vary without any regular pattern probably indicating some kind of dependence between the motion of the director and the molecular dynamical time. Finally, we present our results for the viscosity coefficients. We start by noting that, since viscosity coefficients are defined using collective (as opposed to singleparticle) quantities, very long simulations are needed to obtain convergence. In this sense a proper choice of the coordinate system plays a fundamental role. The validity of the intrinsic frame against the fixed one is proved in Fig. 5, where the values for the TCFI of Eq. (15) are shown as before for runs of increasing length. The integral calculated in the intrinsic frame converges to the correct analytic result, while the one computed in the fixed frame underestimates the TCFI by over 20%. In order to compute the viscosity coefficients we need to evaluate correlations coupling: (a) pressure tensor and angular velocity components (mixed term); (b) homologous components of the angular velocity and (c) different elements of the pressure tensor. These three types of correlations exhibit different decay times, determining the different levels of accuracy with which the coefficients can be estimated within a given trajectory. Correlations of type (a) and (b) decay rapidly: this ensures a clear separation of time scale between the director motion and the corresponding molecular properties. For these integrals, Table 2 reports our results at different system sizes. The agreement with the results of Ref. [4] is completely satisfactory in the first column; a small size dependence is apparent in the last two columns. The third type of TCFI are the most difficult to calculate due to very long relaxation time. In order to get a converged value of the TCFI we integrated the correlation function up to 15r which is probably greater of the recurrence time of our systems: in any case we just integrated till we reach a plateau. This time is considerably large and perhaps begins to be comparable with the time scale of the director motion, especially for the smaller systems. For this reason we are not confident about the results obtained for the two smaller samples. All the necessary symmetry relations are badly satisfied for TCFI obtained for N = 178 and N = 256 samples. This seems to indicate that the converge is difficult to reach. On the contrary, the values of the TCFI for N = 500 are more stable and all the symmetry relations hold within the statistical error. However, our values are quite different from those reported in Ref. [4] as indicated in Table 3. The reason of this discrepancy could be referred to the difference between the GB potential used in Ref. [4] (purely
S. Cozzini et al./Physica A 240 (1997) 173-187
185
N=500
0.54
I
I
I
I
intrinsic fixed 0.52
0.$
0.48
0.46 0.44 0.42
0.4
0.38 500
0
1000
1500 time (reduced units)
2500
2000
3000
Fig. 5. The values of the TCFI q~pa calculated in the intrinsic and fixed frame for N = 500. Table 2 Values of TCFI for systems of different sizes: the values are obtained as averages over the different components and over different sub-runs Size
qSQn
178 256 500 Ref. [4]
0.09 ± 0.09 ± 0.09 ± 0.09 ±
* q~?Pa 0.02 0.01 0.01 0.01
0.486 + 0.490 ± 0.495 ± 0.500 ±
* qP~2 0.008 0.008 0.008 0.001
0.61 ± 0.64 ± 0.65 ± 0.67 ±
0.04 0.03 0.03 0.001
repulsive) and the one we use (the full attractive GB potential). Finally, it is interesting to compare our results with the work o f S m o n d y r e v et al. [7]. In this work the authors investigated our same GB system b y stabilizing the director with an external field. T h e y use the formulation o f Forster [3] for the calculation o f the viscosity coefficients and they express their final results in terms o f Miesowicz viscosity coefficients. The following set o f equations connects the Miesowicz expressions to our coefficients [4]: r/~ = r / +
~ , / 6 + q2 + 71/4,
q~ =r/+
01/6 - 02 -~- ;)1/4 ,
r/~ = t / - 0 , / 3 .
(21)
186
S. Cozzini et al. /Physica A 240 (1997) 173-187 Table 3 Values for integral of type q~/J,/6 Size
?]3333
?]3131
?]111I
500 Ref. [4]
4.0 dz 0.3 3.2 4- 0.2
0.95 ~ 0.08 0.58 ~ 0.04
4.3 4- 0.4 3.3 i 0.1
Table 4 Comparison with Ref. [7] Viscosity
This work (500)
qM
0.2 ± 0.1
qm
14.3 4- 0.8
qM
3.3 4- 0.5
71
5.6 ± 0.6
Ref [7] 2.5 12 3.1 10
Table 4 compares the two sets of results. The values for Ref. [7] are given for a state point very close to ours. We include in this table also the values for the twist viscosity ~l- It must be noted that the authors exchanged the definitions for the coefficients ~/~ and ~/~: here we consider the inverse definitions and thus we exchanged the values. We note that the agreement with two coefficients is quite good but the first Miesowicz coefficient and the twist viscosity are very different. The reason of such different behaviour is not clear. However, we note that in Ref. [7] the statistics over which the coefficients are measured is considerably less than the statistics we produced in order to compute the quantities of interest.
5. Conclusion In this work we presented an intrinsic method to calculate the dynamical properties of a nematic liquid crystal avoiding the use of any external perturbation on the system. We checked that this approach gives the good physical information by comparing our results with those obtained using always a fixed framed with different methods to stabilize the director. The agreement obtained indicates that our method can be used to obtain intrinsic dynamical properties for nematic liquid crystals without perturbing the director motion. The measurements performed in the presence of a pertubation are only correct for the extrapolated values at zero perturbation. In fact, we know that in thermodynamical limit no perturbation is needed to keep the director fixed. The validity of our approach rests on the assumption that there is an effective time scale separation between the director motion and the molecular motion so that the intrinsic properties are still a good formulation for the transport. We found that this assumption is not completely statisfied for the smaller samples where the time over which the director
S. Cozzini et al./Physica A 240 (1997) 173- 187
187
changes significantly becomes comparables with the decay time o f some correlation function.
Acknowledgements We thank E. de Miguel, M.P. Allen and D. Frenkel for useful discussions. The work is partially supported by EEC network contracts number ERBCHRXCT930282 "Organization and dynamics o f molecules in ordered phases" and number ERBCHRXCT930351 "Molecular Dynamics and Monte Carlo simulations o f quantum and classical systems" and by project PB94-1442 from "Ministerio de Educaci6n y Ciencia" D G I C Y T (Spain). The supercomputer center Cineca (Italy) is also acknowledged for the allocation o f computer time trought Icarus project.
References [1] S.R. de Groot, P. Mazur, Nonequilibrium Thermodynamics, Dover, New York, 1984. [2] D.J. Evans, G.P. Morris, Statistical Mechanics of Nonequilibrium Liquids, Academic Press, London, 1990. [3] D. Forster, Ann. Phys. 85 (1974) 505-534. [4] S. Sarman, D.J. Evans, J. Chem. Phys. 99 (1993) 9021. [5] S. Sarman, J. Chem. Phys. 101 (1994) 480. [6] S. Sarman, J. Chem. Phys. 103 (1995) 393. [7] A.M. Smondyrev, G.B. Loriot, R.A. Pelcovits,Phys. Rev. Lett. 75 (1995) 2340. [8] T.C. Lubensky, Phys. Rev. A 2 (1970) 2497. [9] S. Hess, J. Non Equillibrium Therm. I1 (1986) 175. [10] J.G. Gay, B.J. Berne, J. Chem. Phys. 74 (1981) 3316. [11] E. de Miguel, L.F. Rull, M.K. Chalam, K.E. Gubbins, F. van Swol, Molec. Phys. 74 (1991) 405. [12] E. de Miguel, L.F. Rull, M.K. Chalam, K.E. Gubbins, Molec. Phys. 71 (1990) 1223. [13] S. San-nan, D.J. Evans, J. Chem. Phys. 99 (1993) 620. [14] E. de Miguel, L.F. Rull, K.E. Gubbins, Phys. Rev. A 45 (1992) 3813. [15] M.P. Allen, D.J. Tildesley, Computer Simulations of Liquids, Clarendon, Oxford, 1987.
Physica A 240 (1997) 173-187
Intrinsic frame transport for a model of nematic liquid crystal S. C o z z i n i a'*, L . F . R u l l a, G . C i c c o t t i b, G . V . P a o l i n i c'1 a Dpto. Fisica At6miea Molecular y Nuclear Aptdo. 1065, Universidad de Sevilla, 41080 Sevilla, Spain b Dipartimento di Fisica, Universitgt La Sapienza Roma, Piazzale Aldo Moro 2. 00185 Roma, Italy e Laboratoire des Solides Irradibs, Ecole Polytechnique, 91128 Palaiseau Cedex, France
Abstract
We present a computer simulation study of the dynamical properties of a nematic liquid crystal model. The diffusional motion of the nematic director is taken into account in our calculations in order to give a proper estimate of the transport coefficients. Differently from other groups we do not attempt to stabilize the director through rigid constraints or applied external fields. We instead define an intrinsic frame which moves along with the director at each step of the simulation. The transport coefficients computed in the intrinsic frame are then compared against the ones calculated in the fixed laboratory frame, to show the inadequacy of the latter for systems with less than 500 molecules. Using this general scheme on the Gay-Berne liquid crystal model, we evidence the natural motion of the director and attempt to quantify its intrinsic time scale and size dependence. Through extended simulations of systems of different size we calculate the diffusion and viscosity coefficients of this model and compare our results with values previously obtained with fixed director.
1. I n t r o d u c t i o n
In the nematic phase o f a liquid crystal all the molecules are aligned, in average, along a common direction, identified by a unit vector n called director. As a consequence o f the broken symmetry along n transport properties in the nematic phase are considerably more complicated than in the isotropic one. While in the isotropic phase the second rank diffusion and heat conductivity tensors have only one independent component, this number doubles in the nematic phase. Cross coupling terms which are forbidden in the isotropic phase raise the number o f independent components o f the fourth rank viscosity tensor from three to seven [1]. * Corresponding author. E-mail: [email protected].
I Present address: Department of Chemistry, University of Liverpool, Liverpool L69 3BX, U.K. 0378-4371/97/$17.00 Copyright (~) 1997 Elsevier Science B.V. All rights reserved P l l S0378-4371 (97)001 39-8
S. Cozzini et al./Physica A 240 (1997) 173-187
174
The transport coefficients of a nematic system model can be computed by standard Equilibrium Molecular Dynamics (EMD) by evaluating the appropriate Green-Kubo (GK) integrals along the simulation trajectory [2]. GK relations for nematic liquid crystals are straightforward in the case of the diffusion and heat conductivity coefficients, while the derivation of formulae for the viscosity coefficients requires some care. GK relations for viscosity coefficients were derived more than 20 years ago by Forster [3] using projection operator techniques; a more recent formulation was proposed by Sarman and Evans [4] by means of linear response theory. However, it must be noted that the two formulations are not completely equivalent and a satisfactory full formalization of the GK expressions for all the viscosity coefficients is still lacking. From a theoretical point of view, all GK relations are expressed in a coordinate system in which one of the axes is aligned with the director, whose orientation is considered constant in time. This assumption raises a delicate problem in a computer experiment. In the absence of an extemal field the director orientation is a quasiconserved degree of freedom, and is expected to diffuse slowly on the unit sphere. This diffusion can be considerable in samples of typical simulation size (few hundred molecules) and therefore needs to be taken into account. Indeed, neglecting this effect can lead to large errors in the estimation of the transport coefficients. So far, this problem has been circumvented by fixing the director with appropriate constraints [4-6] or by applying external fields to stabilize its orientation [7]. The main drawback of these approaches is that they tend to bias the dynamics of the system. In this work, we propose a new approach which leaves the system unperturbed, by performing our measurements in the director's intrinsic reference frame. This frame is defined at each instant by the position of n. All the physical quantities of interest are calculated with respect to this moving frame. Our investigation is twofold. First, we study the size dependence of the director motion and check for an effective separation in time scale between the director's diffusion and molecular relaxation. Second, we test our approach by calculating the diffusion and the viscosity coefficients for a model system (the Gay-Berne model) and comparing our results with previous studies of similar models. The paper is organized as follows: in Section 2 we introduce the relevant theoretical formulae and summarize the GK relations employed. In Section 3 we describe our computer experiment, while in Section 4 the results are presented and discussed. Finally in Section 5 we draw some conclusions.
2. Formalism and kinematics of the director
The degree of molecular alignment in a liquid crystal is measured by the order parameter S, which is defined as the largest eigenvalue of the second rank-order tensor [8]:
Q=
E 1
UiUi
-
-
1 ,
(1)
S. Cozzini et al./Physica A 240 (1997) 173-187
175
where 1 is the unit tensor, N is the number o f molecules and ui is the unit vector parallel to the axis o f molecule i. The eigenvector corresponding to S is defined as n, the director. The angular velocity o f the director is f2 - n x ft. This definition rests on the fact that for uniaxial molecules the linear velocity is orthogonal to n. To define the intrinsic system associated with the moving director we observe that, by definition, ( Q - S 1 ) - n = 0 . This means that every row or column of the Q - S1 matrix lies in a plane perpendicular to the director. Therefore, choosing the first column o f this matrix we obtain a vector m orthogonal to n. The vector product between m and n defines the remaining axis. It is worth noting that this formulation of the intrinsic frame guarantees its continuity in time as a consequence o f the continuity of Q. In a liquid crystal there are two independent components o f the second rank diffusion tensor D~I3: DIIII, which relates forces and fluxes parallel to the director, and D x ± relating forces and fluxes in the plane perpendicular to the director. It is easy to derive the Green-Kubo relations for these coefficients: Ollll ---
/0
dt (vll(t)Vll(0)),
D±± =
F
dt (Vx(t)v±(O)),
(2)
where vii and Vx are the parallel and perpendicular component of the linear velocity of each molecule with respect to n. Here, for simplicity, we drop out the average over all the molecules o f the sample. The viscosity coefficients can be defined using the phenomenological relations existing between the pressure tensor and the strain rate. To derive them, we note that the entropy production a for a nematic liquid crystal subject to a shear stress can be expressed, using the same notation o f Ref. [5], as a= -l[ps
: (Vn)S + p a . (~7 x u - 2f2) + ½(tr(P) - Peq)V" n],
(3)
where T is the temperature, P is the pressure tensor, u is the velocity field and f2 the angular velocity of the director. The notation A s = l [ A + A r] - 1 tr(A) stands for the symmetric traceless part o f a second rank tensor while A a = ½[A - A r] is its ant±symmetric part. In the above expression the thermodynamical fluxes can be identified as ps, pa and tr(P) - Peq, conjugated, respectively, with the three forces ( V u ) s, V x u - 212 and V . u . The condition of uniaxal symmetry yields the following phenomenological relations between fluxes and forces [4,6,9]: ps
=
_2~/(Vu)S _ 2~ 1[(nn)S. (Vu)~]~ _ 2~3(nn)~[(nn)S : (Vu)~] _ ((nn)SX7. + 2f/2[(nn) ~. ~. (½V × u - f2)] ~ ,
u
(4)
2P" = - ~l [32-1 - (nn)S] - (½V × u - f2) - ~2e : [(nn) s. (Vu)S],
(5)
½[tr(P) - Peq] -- -r/v x~7"u - K(nn)~ "(Vu) ~ .
(6)
Here e is the third-order Levi Civita tensor. The bulk viscosity ~z, the twist viscosity 71 and the shear viscosity ~/ must be positive because they relate tensorial quantities o f the same order and would be the only non-zero terms in an isotropic system. The anisotropic contribution is given by the two coefficients ~1 and 43 which relate the same
S. Cozzini et aL /Physica A 240 (1997) 173-187
176
quantities like r/ but in this case expressed using the order tensor nn [9] to take into account the lower symmetry. All the other terms are cross couplings between tensors of different parity and ranks and can have either sign. The coefficient ~ relates the symmetric traceless part of the pressure tensor to the divergence of the velocity field while q2 relates the same quantity to the angular velocity of the director relative to the rotating fluid. These two coefficients are linked to the coefficients x and 72 by the following Onsager relations: ~2 z ~72/2 and ~ =- K. The number of independent variables is therefore seven. GK expressions for five of these coefficients are derived in Ref. [4] and we will use that formulation throughout this work. In a director-based coordinate system with the third axis along the director, the GK expressions for the five viscosity coefficients are given by
2tlpt2tlQ__.._p r/= ½(2q3131 + 111212)
3~bno '
(7)
2 tlpI2tlK2P
~1 =2(t/3j31 -
q1212) -
700
(8) ,
7]3 -- 3(r13333 .2_ 2r/ill 1 ) _ 2173131 _ r/1212 +
2
-
2 tlPOrlOP ~
,
(10)
~boo ' ~bo~ .;4
(9)
( 11 ) '
where the r//jkt,qSno, q~o and q~Pa are all time-correlation fimction integrals (TCFI) and are given by the following relations: tlijkI = flV c~O0=flV
tl~a = flV
(Pi;(S)PkSl(O))eqds,
/o
(Q2(s)Q2(O))eqds=flV
(P~j(s)fa2(O))eq ds = flV
(12)
io /o
(Ql(S)Ql(O)}eqds,
(13)
(p2s3(s)f~l (0))e0 ds,
(14)
tl~pa=flV fo ~ (p~3(s)f21(O))~qds=flV fo ~c (P~ll(S)Q2(O))eqds= ~1 .
(15)
Here P,] ( ~ ) are elements of the symmetric traceless (antisymmetric) pressure tensor, ~ is the component of the angular velocity of the director, fl = 1/ksT and V is the volume of the simulation cell. Some symmetry relations hold among different r/0kt [4]. We note that the TCFI in Eq. (15) is not an experimental number but a sum rule [4]. Finally, a word of caution must be spent. The above formulae are derived under the assumption that the director remains fixed during the whole simulation. In that sense
177
S. Cozzini et aLIPhysica A 240 (1997) 173 187
a director-based coordinate system is an inertial frame, a condition not necessarily met during simulation. This effectively poses a problem only if the time scale over which the displacement of the director is appreciable is comparable with the decay time of the correlation functions we want to measure. We will discuss this issue further in the light of the results presented in Section 4.
3. Model and computer experiment In this section we briefly discuss the Gay Berne (GB) model used and the simulations performed. The GB model [10] has been widely used as a semi-realistic model for liquid crystals in recent years. It can be thought of as a fluid of pointqike particles, each carrying a unit vector u which mimics the long molecular axis. The particles interact via an anisotropic Lennard-Jones potential, which depends on the relative orientation and position of a pair of molecules. The explicit interaction energy of two molecules i and j is given by
Uij = 4g(t ij, hi, uj)
E(
(70
rij -
o(~g, ui, uj) + (70
12
(7O
r6 - cr(~ii, ui, uj) + (70
' (16)
where r~ is the relative distance between the centers of mass of the molecules i and j and bj = ri;/ro is a unit vector along the direction of the intramolecular distance: r U- - r i - rj. The range parameter a(~,y,u/,nj) is expressed as
{ a(~ij, ui, u/)=6o
Z [(ro"ni%F-ij'Uj)2 1-~
[
l+)~(ui.uj)
+
(Fij.ui--ro'nj)2J } -1':2 l - z(ui.uj)
(17)
where Z = ( K2 - 1)/(K2 + 1) and K is the length to breadth ratio for a prolate ellipsoid of revolution. The strength of the interaction e, is a product of two terms:
(18)
g(r.ij, ui, nj)--~ F,O[gl(ni, nj)] v )< [~2(Ftj, ni, nj)] 'u
with < (u~, uj) = [1 -
z2(ui • uj):]-1/2
~:2(~0,ui, u j ) =
1- ~
(19)
and
1 - z'(ni'u/)
J
Here •' = [(K') 1/" - 1]/[(~') l/" + 1] and s:' is the ratio of the potential energy minima of the side-by-side and the end-to-end configuration. The phase diagram of this fluid has been extensively studied by de Miguel et al. [11,12].
S. Cozzini et aL /Physica A 240 (1997) 173-187
178
Recently, some transport coefficients were measured for the GB-models by different methods. Sarman and Evans used Equilibrium and Non-Equilibrium Molecular Dynamics simulations to measure diffusion coefficients, heat fluxes [13] and viscosities [4] for a purely repulsive GB system. De Miguel et al. [14] measured the diffusion coefficients by EMD while Smondyrev et al. performed EMD calculations [14] of viscosity coefficients. We therefore chose this model for our computer experiment in order to check our results against a relatively wide spectrum of data obtained with different techniques. Table 1 List of the runs. Here T is the averaged temperature and the S the averaged order parameter Name
Size
Run length (~)
T
S
T PL B LC XX YY
178 178 256 256 500 500
3250 4000 4000 4000 3000 3000
1.242 ± 0.009 1.260 ± 0.010 1.186 4- 0.009 1.244 4- 0.010 1.259 i 0.006 1.260 ± 0.006
0.77 4- 0.02 0.77 4- 0.02 0.774-0.02 0.76 4- 0.02 0.75 4- 0.02 0.75 ± 0.02
i
178
0
................................................................
-1
-2
-3
(a)
I
i
I
I
-4
-2
2
4
Fig. 1. The director motion for systems of different size. See text for explanation.
- -
P I
3 I
'(p~nupuo3) 1 ':~!2, (o) 0 ]
3-
'¢8-
I
t-
00£ i
i
i
£
i
(q) #
3
I
0
I
17-
3l
{
g-
3"
993 I
6LI
ZgI-EZI (Z661) OPE V va!g,fVd/"l n la lUlZZOD "~
180
S. Cozzini et al./Physica A 240 (1997) 173-187
MEAN SQUARED 0.008
J
0.006
- ....
DISPLACEMENT i
OF THE DIRECTOR F
N=178 N=256
r
j
"
o
"E 0.004 I
¢-
0.002
0.000
0.0
2000
4000
6000
8000
10000
time ( in timesteps) Fig. 2. The director MSD for systems of different size.
In all our simulations we use the standard anisotropy parameter Z = 3 and Z~-- 5 and exponents p = 2 and v = 1. The simulation results are given in length, energy, mass and time units of a0, e0, m and z--ao(m/eo) 1/2. Three systems of different sizes were investigated: 178, 256 and 500 molecules, all at the same reduced density (p* --0.345) and the same starting reduced temperature T * = 1.25. The stability of the nematic phase at this state point is well known [ll,12]. For the largest systems (256 and 500 molecules) the cutoff distance beyond which the potential was set equal to zero was 4.0a0. For the smallest system (N = 178 molecules) the correspondingly smaller box size dictated a slightly shorter spherical cutoff distance r~* = 3.8a0. Consequently the smallest system does not correspond exactly to the larger ones, but the nematic phase is still stable and the thermodynamical properties of the three systems are virtually the same. We performed constant energy molecular dynamics simulations with the usual periodic boundary conditions. The integration algorithm used is a standard velocity Verlet for the translational degrees of freedom while for the rotational degrees we adapted the standard algorithm for uniaxial molecules [15] to the velocity Verlet scheme. The time step was 0.001z and reduced moment of inertia was set to unity. For each size we performed two separate runs starting from independently equilibrated
S. Cozzini et al./Pkysica A 240 (1997) 173 187
181
Diffusion at short times 0.050
i
i
a--od
par(int)
~-----~ d p a r (fix) *-----~ d p a r (fix) o
--~ d p a r ( i n t )
0.040
¢-1 "O ~) ¢.) -1
0.030
0.020
I
0.0100. 0
50.0
I
100.0 150.0 time ( reduced units)
I
I
200.0
250.0
Fig. 3. The diffusion coefficients measured over a short run for N = 178. Here Dper and Dp,r indicate, respectively, D ± ± and DII!I, while int and fix means intrinsic and fixed frame. See text for discussion.
configurations with different director orientation. Table 1 summarizes the simulations performed. The length of the runs is necessarily large in the attempt to get good convergence of the results. Along each trajectory we calculated the TCFI over a correlation time of 5~. The correlations were also computed in a director-based coordinate system fixed throughout the simulation and defined by the position of the director at the beginning of the run. In order to monitor the director's dynamics its orientation was stored every ten time steps (0.01~) for each run.
4. Results The first aspect we are interested in is the dependence of the motion of the director on the size of the system. A representation of this motion is obtained by describing the position of n(t) via the angles 0 and 4) formed with respect to the fixed frame associated to n(t = 0). In Fig. 1 the quantities tg(0)cos(~b) and tg(0)sin(~b) are plotted at each recorded time step. The starting point of the trajectory is by definition at the origin. The director's drift can be better appreciated because of the factor tg(0).
S. Cozzini et al./Physica A 240 (1997) 173-187
182
For instance, an azimuthal angle 0 - - 9 0 ° corresponds to infinite displacement. Fig. 1 clearly indicates a size dependence of the motion of the director. Figs. l(a) and (b) correspond to N = 178 and N = 2 5 6 , respectively, and Fig. l(c) refers to N - - 5 0 0 . More quantitative information is contained in the mean square displacement (MSD) of the director, reported in Fig. 2, where again the size dependence is apparent. We are currently running some tests with larger systems in order to see if a general relationship between the director's diffusion and the system size can be inferred. Let us discuss now the results found for the diffusion coefficients. Fig. 3 shows the values of the two diffusion coefficients DII II and D ± ± calculated in the two different reference frames. The length of the trajectory over which the TCF1 is calculated is reported in the abscissas. By increasing the statistics we expect the average to stabilize around the final value. This is indeed what we observe for the two diffusion coefficients FlintII is twice as large as Flint ~ x x in agreement measured in the intrinsic frame. In this case ~11 with previous results [14,13]. The behaviour of the coefficients measured in the fixed frame is however different. At short times the diffusion coefficients measured in the two frames coincide within the error. By increasing the statistics DilfixII and D ~ eventually deviate from their initial values. This result is a clear signature of the motion of the
N=178
0.050
i
i
d_par z...---~ d p a r *---.--~ d _ p e r o o d_per
(Int) (fix) (fix) (Int)
0.040
¢.-
o "1o
0.030
0.020
n
0'01500.0 (a)
i
L
1500.0
2500.0
3500.0
time ( reduced units)
Fig. 4. The diffusion coefficientscalculated along the whole runs for systemsof different size: (a) 178 atoms, (b) 256, (c) 500. Labels are as in previous figure.
S. Cozzini et al./Physica A 240 (1997) 173-187
183
N=256 0.050
i
i
o o = ...... *------~ ~
d__par (Int) d _ p a r (fix) d _ p e r (fix) d _ p e r (Int)
0.040
¢9 t"0
0.030
¢O
0.020
0.010 0,0
' 1000.0
2000.0
3000.0
time ( reduced units)
(b)
N=500 0.050
i
o
o d _ p a r (int)
~*----~ d _ p a r (fix) c--
•.-> d _ p e r (fix) ~ d _ p e r (Int)
0.040
ro
0.030
o
0.020
0.010 0.0
(c)
r 1000.0
,
2000.0 time ( reduced units)
Fig. 4 (continued).
3000.0
184
S. Cozzini et al./Physica A 240 (1997) 173 187
director. After about 100r the initial director orientation with respect to which we are measuring our diffusion coefficients (fixed frame) is no longer the symmetry axis of our system. Fig. 4 shows the same quantities measured for longer runs of about 3000r for the three systems studied. The behaviour exhibited over short trajectories is here more pronounced due to the longer times. In all the three cases the diffusion coefficients drift away from their initial values, the more rapidly the smaller the size of the system. This behaviour is not surprising and is due to the fact that the direction along which we are measuring our coefficients is losing correlation with the symmetry axis of the system. It could be expected that the diffusion coefficients measured parallel and perpendicularly to this arbitrary direction would eventually tend to roughly the same value. Indeed this seems the case for N = 500. For the smaller samples the values of the coefficients vary without any regular pattern probably indicating some kind of dependence between the motion of the director and the molecular dynamical time. Finally, we present our results for the viscosity coefficients. We start by noting that, since viscosity coefficients are defined using collective (as opposed to singleparticle) quantities, very long simulations are needed to obtain convergence. In this sense a proper choice of the coordinate system plays a fundamental role. The validity of the intrinsic frame against the fixed one is proved in Fig. 5, where the values for the TCFI of Eq. (15) are shown as before for runs of increasing length. The integral calculated in the intrinsic frame converges to the correct analytic result, while the one computed in the fixed frame underestimates the TCFI by over 20%. In order to compute the viscosity coefficients we need to evaluate correlations coupling: (a) pressure tensor and angular velocity components (mixed term); (b) homologous components of the angular velocity and (c) different elements of the pressure tensor. These three types of correlations exhibit different decay times, determining the different levels of accuracy with which the coefficients can be estimated within a given trajectory. Correlations of type (a) and (b) decay rapidly: this ensures a clear separation of time scale between the director motion and the corresponding molecular properties. For these integrals, Table 2 reports our results at different system sizes. The agreement with the results of Ref. [4] is completely satisfactory in the first column; a small size dependence is apparent in the last two columns. The third type of TCFI are the most difficult to calculate due to very long relaxation time. In order to get a converged value of the TCFI we integrated the correlation function up to 15r which is probably greater of the recurrence time of our systems: in any case we just integrated till we reach a plateau. This time is considerably large and perhaps begins to be comparable with the time scale of the director motion, especially for the smaller systems. For this reason we are not confident about the results obtained for the two smaller samples. All the necessary symmetry relations are badly satisfied for TCFI obtained for N = 178 and N = 256 samples. This seems to indicate that the converge is difficult to reach. On the contrary, the values of the TCFI for N = 500 are more stable and all the symmetry relations hold within the statistical error. However, our values are quite different from those reported in Ref. [4] as indicated in Table 3. The reason of this discrepancy could be referred to the difference between the GB potential used in Ref. [4] (purely
S. Cozzini et al./Physica A 240 (1997) 173-187
185
N=500
0.54
I
I
I
I
intrinsic fixed 0.52
0.$
0.48
0.46 0.44 0.42
0.4
0.38 500
0
1000
1500 time (reduced units)
2500
2000
3000
Fig. 5. The values of the TCFI q~pa calculated in the intrinsic and fixed frame for N = 500. Table 2 Values of TCFI for systems of different sizes: the values are obtained as averages over the different components and over different sub-runs Size
qSQn
178 256 500 Ref. [4]
0.09 ± 0.09 ± 0.09 ± 0.09 ±
* q~?Pa 0.02 0.01 0.01 0.01
0.486 + 0.490 ± 0.495 ± 0.500 ±
* qP~2 0.008 0.008 0.008 0.001
0.61 ± 0.64 ± 0.65 ± 0.67 ±
0.04 0.03 0.03 0.001
repulsive) and the one we use (the full attractive GB potential). Finally, it is interesting to compare our results with the work o f S m o n d y r e v et al. [7]. In this work the authors investigated our same GB system b y stabilizing the director with an external field. T h e y use the formulation o f Forster [3] for the calculation o f the viscosity coefficients and they express their final results in terms o f Miesowicz viscosity coefficients. The following set o f equations connects the Miesowicz expressions to our coefficients [4]: r/~ = r / +
~ , / 6 + q2 + 71/4,
q~ =r/+
01/6 - 02 -~- ;)1/4 ,
r/~ = t / - 0 , / 3 .
(21)
186
S. Cozzini et al. /Physica A 240 (1997) 173-187 Table 3 Values for integral of type q~/J,/6 Size
?]3333
?]3131
?]111I
500 Ref. [4]
4.0 dz 0.3 3.2 4- 0.2
0.95 ~ 0.08 0.58 ~ 0.04
4.3 4- 0.4 3.3 i 0.1
Table 4 Comparison with Ref. [7] Viscosity
This work (500)
qM
0.2 ± 0.1
qm
14.3 4- 0.8
qM
3.3 4- 0.5
71
5.6 ± 0.6
Ref [7] 2.5 12 3.1 10
Table 4 compares the two sets of results. The values for Ref. [7] are given for a state point very close to ours. We include in this table also the values for the twist viscosity ~l- It must be noted that the authors exchanged the definitions for the coefficients ~/~ and ~/~: here we consider the inverse definitions and thus we exchanged the values. We note that the agreement with two coefficients is quite good but the first Miesowicz coefficient and the twist viscosity are very different. The reason of such different behaviour is not clear. However, we note that in Ref. [7] the statistics over which the coefficients are measured is considerably less than the statistics we produced in order to compute the quantities of interest.
5. Conclusion In this work we presented an intrinsic method to calculate the dynamical properties of a nematic liquid crystal avoiding the use of any external perturbation on the system. We checked that this approach gives the good physical information by comparing our results with those obtained using always a fixed framed with different methods to stabilize the director. The agreement obtained indicates that our method can be used to obtain intrinsic dynamical properties for nematic liquid crystals without perturbing the director motion. The measurements performed in the presence of a pertubation are only correct for the extrapolated values at zero perturbation. In fact, we know that in thermodynamical limit no perturbation is needed to keep the director fixed. The validity of our approach rests on the assumption that there is an effective time scale separation between the director motion and the molecular motion so that the intrinsic properties are still a good formulation for the transport. We found that this assumption is not completely statisfied for the smaller samples where the time over which the director
S. Cozzini et al./Physica A 240 (1997) 173- 187
187
changes significantly becomes comparables with the decay time o f some correlation function.
Acknowledgements We thank E. de Miguel, M.P. Allen and D. Frenkel for useful discussions. The work is partially supported by EEC network contracts number ERBCHRXCT930282 "Organization and dynamics o f molecules in ordered phases" and number ERBCHRXCT930351 "Molecular Dynamics and Monte Carlo simulations o f quantum and classical systems" and by project PB94-1442 from "Ministerio de Educaci6n y Ciencia" D G I C Y T (Spain). The supercomputer center Cineca (Italy) is also acknowledged for the allocation o f computer time trought Icarus project.
References [1] S.R. de Groot, P. Mazur, Nonequilibrium Thermodynamics, Dover, New York, 1984. [2] D.J. Evans, G.P. Morris, Statistical Mechanics of Nonequilibrium Liquids, Academic Press, London, 1990. [3] D. Forster, Ann. Phys. 85 (1974) 505-534. [4] S. Sarman, D.J. Evans, J. Chem. Phys. 99 (1993) 9021. [5] S. Sarman, J. Chem. Phys. 101 (1994) 480. [6] S. Sarman, J. Chem. Phys. 103 (1995) 393. [7] A.M. Smondyrev, G.B. Loriot, R.A. Pelcovits,Phys. Rev. Lett. 75 (1995) 2340. [8] T.C. Lubensky, Phys. Rev. A 2 (1970) 2497. [9] S. Hess, J. Non Equillibrium Therm. I1 (1986) 175. [10] J.G. Gay, B.J. Berne, J. Chem. Phys. 74 (1981) 3316. [11] E. de Miguel, L.F. Rull, M.K. Chalam, K.E. Gubbins, F. van Swol, Molec. Phys. 74 (1991) 405. [12] E. de Miguel, L.F. Rull, M.K. Chalam, K.E. Gubbins, Molec. Phys. 71 (1990) 1223. [13] S. San-nan, D.J. Evans, J. Chem. Phys. 99 (1993) 620. [14] E. de Miguel, L.F. Rull, K.E. Gubbins, Phys. Rev. A 45 (1992) 3813. [15] M.P. Allen, D.J. Tildesley, Computer Simulations of Liquids, Clarendon, Oxford, 1987.