Excess volume of homonuclear diatomic mixtures

Journal of Molecular Liquids 125 (2006) 22 – 28 www.elsevier.com/locate/molliq

Excess volume of homonuclear diatomic mixtures Pavel Mora´vek *, Jirˇ´ı Kolafa, Toma´x Hujo, Stanislav Labı´k Department of Physical Chemistry, Institute of Chemical Technology, Technicka´ 5, 16628 Prague 6, Czech Republic Available online 27 December 2005

Abstract A new fast and robust numerical method of solving the Ornstein – Zernike integral equation is described and tested on the binary mixture of hard diatomic molecules. The method combines direct iterations and the Newton-GMRES algorithm. Excess volumes of a mixture of hard homonuclear dumbbells with reduced elongations L* = 0.3 and 0.6 are calculated and compared to computer simulations and several equations of state including the virial expansion up to B 5. D 2005 Elsevier B.V. All rights reserved. Keywords: Homonuclear diatomic; Excess volume; Ornstein – Zernike; Newton-GMRES; NITSOL; Virial coefficient

1. Introduction The numerical solution of the Ornstein– Zernike integral equation has been studied for many years. Several classic methods are known. The simplest but inefficient direct iteration (Picard) method with damping [1] can be improved by variable relaxation parameter [2]. Considerable improvement for simple liquids is achieved by the hybrid Picard/Newton– Raphson (NR) methods as Gillan [3] and Labı´k et al. [4] presented. The solution of the OZ equation for molecular fluids and their mixtures is more complicated. Some calculations based on Picard/NR algorithm were done by Labı´k et al. [5], using the earlier work of Lado [6] for this application. The hybrid Picard/NR algorithm, which eliminates the slow convergence of Picard iterations and the necessity of the large memory to keep many variables for the NR method, can be improved by replacing the Newton – Raphson algorithm by the Newton-GMRES (Generalized Minimal RESidual method) algorithm. Mathematical basics of GMRES are described in Ref. [7]. Various implementations of this algorithm suitable for solving the Ornstein – Zernike equation have been compared by Booth et al. [8]. In this work we have used the Newton-GMRES algorithm implemented in the public domain Krylov solver NITSOL [9] combined with the classic Picard iterations to obtain thermo* Corresponding author. E-mail address: [email protected] (P. Mora´vek). URL: http://www.vscht.cz/fch/en/people/ (P. Mora´vek, J. Kolafa, T. Hujo, S. Labı´k). 0167-7322/$ - see front matter D 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.molliq.2005.11.013

dynamical properties of the binary mixture of hard homonuclear diatomics with reduced elongations L* = 0.3 and 0.6 using six different closures of the Ornstein– Zernike equation. The kernel of the algorithm is based on the work of Lado [6]. It consists in expanding the correlation functions in spherical harmonics and writing the Ornstein – Zernike equation as a set of nonlinear equations in Fourier space. The results have been compared in terms of the excess volume with the classic equation of state proposed by Boublı´k [10], two modern ones (Largo and Solana [11] and Boublı´k [12]), and Monte Carlo simulations in the NPT ensemble, and the virial expansion up to the fifth order. 2. Methods 2.1. The Ornstein –Zernike equation for a mixture of linear molecules The OZ equation [13] for linear molecules is given by Z q cð12Þ ¼ cð13Þ½cð32Þ þ cð32Þdð3Þ; ð1Þ 4k where q is the number density, c(12) the direct correlation function, c(12) = h(12)  c(12) the series function, where h(12) = g(12)  1 is the total correlation function, and d(13) denotes integration over all positions and orientations of molecule 3. The OZ equation must be coupled with a closure, which may be expressed in the form cð12Þ ¼ exp½  buð12Þ þ cð12Þ þ Bð12Þ  cð12Þ  1:

ð2Þ

P. Mora´vek et al. / Journal of Molecular Liquids 125 (2006) 22 – 28

The function B(12) is known as the bridge function, u(12) is the intermolecular potential and b = 1 / kT. For better numerical solution, and extended to a binary mixture, the OZ equation may be expressed as Z   qa caa ð13Þ cab ð32Þ þ cab ð32Þ dð3Þ cab ð12Þ ¼ 4k Z   qb cab ð13Þ cbb ð32Þ þ cbb ð32Þ dð3Þ; þ ð3Þ 4k where indices a and b label the components. In order to decrease the number of variables, any correlation function X(12) can be expanded in normalized spherical harmonics [6], X ð12Þ ¼ 4p

V X V n X X

Xl1 l2 m ðrÞYl1 m¯ ðh1 ; /1 ÞYl2 m¯ ðh2 ; /2 Þ;

l1 ¼0 l2 ¼0 m¼n

ð4Þ where n = min(l 1, l 2), and m ¯ =  m. The harmonic coefficients derived from Eq. (4) are Z Z Z 1 1 1 k Xl1 l2 m ðrÞ dcosh1 dcosh2 d/12 X ð12ÞYl1 m¯ ðh1 ; /1 Þ 2 1 1 k  Yl2 m ðh2 ; /2 Þ: ð5Þ The OZ equation for linear molecules in the Fourier space is given by [6] Z q c˜ ð12Þ ¼ c˜ ð13Þ½c˜ ð32Þ þ c˜ ð32Þdð3Þ; ð6Þ 4k where the harmonic coefficients of c˜ and c˜ from Eq. (4) are written as c˜ ð12Þ ¼ 4k

V X

minX ðl1 ;l2 Þ

c˜ l1 l2 m ðk ÞYl1 m ðh1 ; /1 ÞYl2 m¯ ðh2 ; /2 Þ;

l1 l2 ¼0 minðl1 ;l2 Þ

ð7Þ

c˜ ð12Þ ¼ 4k

V X

minX ðl1 ;l2 Þ

c˜ l1 l2 m ðk ÞYl1 m ðh1 ; /1 ÞYl2 m¯ ðh2 ; /2 Þ:

l1 ;l2 ¼0 minðl1 ;l2 Þ

ð8Þ From the last three expressions one can derive a relation suitable for our algorithm c˜ l1 l2 m ðk Þ ¼ ð  Þm q

V X

  c˜ l1 l3 m ðk Þ c˜ l3 l2 m ðk Þ þ c˜ l3 l2 m ðk Þ :

ð9Þ

l3 ¼m

Finally, after extending Eq. (9) to a binary mixture, m c˜ ab l1 l2 m ðk Þ ¼ ð  Þ q

V X

V X

Let us summarize the methods of Lado [6] and LabNk et al. [5]. The variables r and k are discretized at values r i = iDr, k j = jDk, where i, j = 1,. . .N. The OZ equation then becomes a set of nonlinear equations. The algorithm is based on the successive solution of harmonic coefficients defined by Eq. (5). Both methods start from an initial estimate of c˜ l 1l 2m(k j ) and each iteration consists of the following steps: & Transformation to the configuration space, c˜ l 1l 2m (k j ) Y c l 1l 2m (r i ); & Application of closure, c l 1l 2m (r i ) Y c l 1l 2m (r i ); & Transformation to the Fourier space, c l 1l 2m (r i )c˜l 1l 2m (k i ); & Using the OZ equation, c˜l 1l 2m (k i ) Y c˜l 1l 2m (k j ). The Hankel transformation is used for the direct and inverse transformation of c˜ l 1l 2m (k j ) and c l 1l 2m (r i ). This part is described by Lado [6]. In detail, harmonic coefficients are transformed from the intermolecular coordinate frame to the space-fixed coordinate frame. Then the Hankel (Fourier –Bessel) transformation is performed and obtained coefficients are transformed back to the intermolecular coordinate frame. The OZ equation is computed as the set of three equations (10) for indices aa, ab and bb (two pure fluids and a mixture). The final results of c˜ are taken as the initial estimate for the next iteration cycle. A relaxation parameter is incorporated to guarantee convergence at higher densities. The application of closure is the most time-consuming step in the algorithm described above because it involves numerical calculations of a large set of nonlinear equations. The number of equations depends on N and the number of harmonic coefficients. For example, for the expansion in 125 harmonic coefficients (to l 1l 2m = 444) of binary mixture, and each coefficient adopting about 1000 discretized values, we arrive at a set of about 125,000 linear equations. This set of equations may be solved by the Newton– Raphson method but it is efficient only for small sets. For this reason, the NR method is used only for calculating the most significant values of the correlation functions. The remaining values are calculated by direct iterations. Our improvement of this method is based on new mathematical algorithms for solution of large sets of nonlinear equations. Booth et al. [8] successfully used the NewtonGMRES algorithm implemented in nonlinear Krylov solvers NKSOL [14] and NITSOL [9]. This algorithm is based on the classic Newton iterative algorithm, where trial steps are obtained by a Krylov subspace method GMRES. A nonlinear system defined by ð11Þ

can be solved by the Newton method in an iteration i as

l3 ¼m

þ ð  Þm q

2.2. Hybrid Picard iteration and Newton-GMRES algorithm

F ðxÞ ¼ 0;

h i ab ab ˜ ˜ ð k Þ c ð k Þ þ c ð k Þ c˜ ab l1 l3 m l3 l2 m l3 l2 m

23

F ðxi Þ ¼  J ðxi ÞI Dxi :

h i ˜ ab ˜ ab c˜ ab l1 l3 m ðk Þ c l 3 l 2 m ðk Þ þ c l3 l2 m ðk Þ :

l3 ¼m

ð10Þ

ð12Þ

Common methods need whole huge Jacobian matrix N  N. The advantage of methods based on the Krylov algorithm consists in using only a product of Jacobian and the vector of

24

P. Mora´vek et al. / Journal of Molecular Liquids 125 (2006) 22 – 28

residuals J ˙ Dx i . This product is calculated numerically by the equation J ðxi ÞIDxi ¼

F ðxi þ eDxi Þ  F ðxi Þ ; e

ð13Þ

where e is a small number. The next problem in classic Newton method is the determination of the relaxation parameter g i in each step to provide convergence of the solution xiþ1 ¼ xi þ gi Dxi :

ð14Þ

NITSOL offers several terms for the relaxation parameter. We used the default term jjFcur jj ¼

jjFprev þ J prev IDxjj : jjFprev jj

ð15Þ

Another advantage of the NITSOL algorithm is possibility to calculate all values of discretized correlation functions of mixture within this algorithm, not the most significant ones by the Newton-GMRES method and the remaining by Picard iterations. The time consumption is not very different and there is no need for determination of how many values will be calculated in the Newton-GMRES algorithm. There are several closure approximations implemented in our program: classic closures of Percus and Yevick [15], HNC [16], Martynov and Sarkisov [17], Rogers and Young [18] and two modern ones—the closure of Duh and Henderson [19] and the modified Verlet closure [20]. Most calculations were performed with 125 harmonic coefficients (to l 1l 2m = 444). Comparison with calculations using only 34 coefficients (to 222) in Table 1 shows that 125 coefficients are sufficient for both elongations. Although the excess volume is a difference of close numbers, the errors caused by truncation are systematic and should not affect the results dramatically. Problems could arise for much longer elongations. 2.3. Excess volumes Various thermodynamics properties are of interest in the mixture of diatomics but a fast method to obtain them is needed. We would like to compare the speed and accuracy with MC calculations and values computed from equations of state. Therefore we chose a ‘‘difficult’’ thermodynamic property, the excess volume, because it is very close to zero for the hard diatomic model. Table 1 Compressibility factor from OZ algorithm truncated to coefficients 222 and 444 (q = 0.5, Duh – Henderson closure) as a function of composition x 1 x1

Z 222

Z 444

0.0 0.2 0.4 0.6 0.8 1.0

12.78 10.66 9.06 7.86 6.97 6.27

12.83 10.69 9.07 7.87 6.98 6.29

The calculation proceeds as follows. First, the contact values of the pair distribution functions were calculated from our program. In the next step the compressibility factor Z is given by the sum. 2 N XX 3 Z ¼1þ k xi xj bgij;c rij;c x 1 ;x 2 ; ð16Þ 3 V i j where index c means the contact value and expression b>x 1x 2 is the unweighted average over all orientations of molecules 1 and 2. The excess volume of the mixture is given by V E ¼ V  x1 V1  x2 V2 ;

ð17Þ

where V is the volume of the mixture and Vi volume of pure fluid of diatomics of species i at the same reduced pressure bp. Throughout the paper we use units reduced by sphere diameter and the volumes are expressed per one particle, i.e., V = Vrealr  3 / N, where Vreal is the volume of the box with N molecules, and r is the diameter of the spheres composing dumbbells. Volumes are calculated form equation bpV ¼ Z ð1=V ; xi Þ;

ð18Þ

where the function Z is given by values obtained by Eq. (16) and x 1 is the mole fraction of species 1. It appears that the excess volume of our mixture is an almost symmetric function in x 1 so that it can be correlated by the formula V E ¼  Ax1 ð1  x1 Þ:

ð19Þ

Values of A will therefore facilitate comparison of different methods. 2.4. Equations of state Three equations of state were compared with our results. The first two [10,12] are based on the SPT theory. Hard convex bodies are described by the nonsphericity parameters a¼

rs ; 3qv

b¼

qs2 9qv2

ð20Þ

where r, s and v are geometric functionals of the mixture given by X y¼q xi Yi ; Yi ¼ Ri ; Si ; Vi ; Qi ; X q¼q xi Qi ; Qi ¼ R2i ; ð21Þ and Vi and S i are the volume and surface of the dumbbell i and R i is the mean curvature integral divided by 4k. EOSs are defined by Z 1989 ¼

1 3ag þ 1  g ð1  gÞ2 þ

g2 ½ð49a  31Þ  gð11a  7Þ  g2 ð25a  21Þ ð 1  gÞ 3 ð22Þ

P. Mora´vek et al. / Journal of Molecular Liquids 125 (2006) 22 – 28

25

0

-0.015

-0.02

VE

-0.01

VE

-0.025

βp=0.5

-0.02

βp=2

-0.03

-0.03

β p=2 -0.035 0.2

0.3

0.4

0.5

0.6

0.7

0.8

βp=4 0

0.2

0.4

Fig. 1. Excess volume calculated from various closures of the OZ equation. Solid line: the Percus – Yevick’s closure, dashed line: HNC, dotted line: the modified Verlet’s, dash-dotted line: the Martynov – Sarkisov’s, double-dashed line: Rogers – Young’s, thick solid line: the Duh – Henderson’s closure.

0.8

1

Fig. 2. Pressure dependence of the excess volume. Results from the Ornstein – Zernike integral equation with Duh and Henderson’s closure. Solid line: bp = 4, dashed line: bp = 2, dotted line: bp = 5.

were comparable, subsequent calculations at bp = 0.5 and 2 were performed only for 500 particles.

and Z 2002 ¼

1 3ag g2 ½3bð1  2gÞ þ 5ag þ þ ; 1  g ð 1  gÞ 2 ð 1  gÞ 3

2.6. Virial expansion ð23Þ The virial expansion for a binary mixture reads as

where g is the packing fraction and the superscripts of Z refer to the year. The third EOS was derived by Largo and Solana [11]. This EOS is based on effective parameters g ef, v ef and a ef. These parameters represent the dumbbell real shape better than the Boublı´k parameters designed for general HCB. The effective molecular volume is given by [23] vef ¼

0.6

x1

x1

  kr3 1 1 þ 3L4  L43  3h4h ; 2 6

ð24Þ

where the reduced elongation L* = L / r, h* = h / r = (1  L*2 / 4)1 / 2 and h = arcsin(L* / 2). The form of the EOS is  HS  ef   LS ¼ 1 þ aef Zmix mix Zmix gmix  1 :

Z ¼1þ

V X

Bn ðx1 ÞV 1n ;

ð26Þ

n¼2

where the mixing rule for B i is n  X n nk k x x B : k 1 2 nk;k k¼0

Bn ðx1 Þ ¼

ð27Þ

The mixture virial coefficients B n k,k to the order n = 5 were calculated by the method used by Enciso et al. [22] for the mixture of hard spheres. The virial coefficients are represented as sums of Ree –Hoover diagrams (called also modified stars, complement blocks, etc.). In the MC integration a linear chain,

ð25Þ 0

2.5. Monte Carlo simulation

-0.001

VE

The MC simulations provide ‘‘experimental’’ results for model systems. We used the standard cubic periodic simulation box and NPT ensemble [21]. The acceptance ratios were adjusted to 0.25 for displacements, and about 0.4 for rotations and volume changes. The calculations started from a random configuration in a very large box and were followed by a sufficient period of equilibration until measurements were started. The NPT ensemble for hard body systems is controlled by a single parameter bp and gives the averaged volume. Both pure compounds as well as the mixtures at a range of compositions were simulated and the excess volumes calculated from Eq. (17). The MC calculations were first performed for 1000 and 500 particles at the reduced pressure bp = 4. Because the results

βp=0.5

-0.002

0

0.2

0.4

x1

0.6

0.8

1

Fig. 3. Excess volume at reduced pressure bp = 5. Thick dashed line: data from virial coefficients, thick solid line with standard errors: MC calculations, dashdotted line: Boublı´k’s old EOS, dotted line: Boublı´k’s new one.

26

P. Mora´vek et al. / Journal of Molecular Liquids 125 (2006) 22 – 28 0.004 βp=0.5

0

βp=2

0.002 βp=4

0

VE

VE

-0.0005

-0.002 -0.004

-0.001

-0.006 βp=2

βp=0.001

-0.008 0

0.2

0.4

x1

0.6

0.8

1

0

0.2

0.4

x1

0.6

0.8

1

Fig. 4. Excess volume at reduced pressure bp = 2. Thick dashed line: data from virial coefficients, thick solid line with standard errors: MC calculations, dashdotted line: Boublı´k’s old EOS, dotted line: Boublı´k’s new one.

Fig. 6. Pressure dependence of the excess volume. Results from the Largo and Solana EOS. Solid line: bp = 4, dashed line: bp = 2, dotted line: bp = 0.5, dashdotted line: bp = 0.001.

serving as the spanning diagram, is generated by reptation. In one reptation step a molecule (with probability 1 / 2 of species 1 and with probability 1 / 2 of species 2) is added to the chain head so that it overlaps with the head molecule; this molecule becomes a new head while one molecule is removed from the tail. The neighboring molecules in the spanning chain are thus connected by f-bonds. Then the remaining f- and e-bonds are calculated and the counter of the corresponding Ree – Hoover diagram is advanced by 1. From the counters and from the known value of the spanning chain the Ree – Hoover integrals are calculated. We generated 400 cycles by 109 MC steps. The statistical errors were calculated by the block method because the individual MC steps are correlated. The excess volumes are calculated as in Section 3 with Eq. (16) replaced by (26).

The GMRES-based algorithm worked about 15 times faster than direct iterations [1] or their improved version [2], in addition, it is much more stable even at high densities and does not require adjustment of parameters ‘‘by hand’’. One iteration lasted 1.8 s on a workstation with the Athlon64 3000+ processor. There are typically several tens of iterations needed for one density (see below for details). The resulting time is several hours if Eq. (18) is solved for a set of compositions. The results from various closures for bp are compared in Fig. 1. It is seen that the two simplest closures, Percus – Yevick’s and HNC, bound other more sophisticated closures. The best closure is one from the authors Duh and Henderson [19] which seems to be slightly better than Verlet’s modified closure [20] for hard diatomics. Closures were compared in the previous work [24]. The value of parameter A, see (19), is around  0.1. It was found that the excess volume goes to more negative values

3. Results and discussion In the first part of this work we calculated the excess volumes from the Ornestein – Zernike equation. Several classic and modern closures were implemented in this algorithm.

VE

0.0003

0

-0.0003

-0.0006 βp=4 0

0.2

0.4

x1

0.6

0.8

1

Fig. 5. Excess volume at reduced pressure bp = 4. Thick dashed line: data from virial coefficients, thick solid line with standard errors: MC calculations, dashdotted line: Boublı´k’s old EOS, dotted line: Boublı´k’s new one.

Table 2 Reduced virial coefficients of the mixture of homodumbbells k

l

Bkl r  3(k+l  1)

2 1 0

0 1 2

3.100038(7) 3.615854(9) 4.198210(12)

3 2 1 0

0 1 2 3

5.95090(6) 7.23093(7) 8.78144(7) 10.65826(9)

4 3 2 1 0

0 1 2 3 4

8.29383(21) 10.20819(26) 12.56452(29) 15.46531(32) 19.03672(39)

5 4 3 2 1 0

0 1 2 3 4 5

9.5875(6) 11.8374(8) 14.6146(9) 18.0425(9) 22.2731(12) 27.4936(15)

P. Mora´vek et al. / Journal of Molecular Liquids 125 (2006) 22 – 28

27

0.0004

0.001 2

0

0 3

VE

VE

-0.0004

-0.001

4 -0.0008 5 -0.0012

-0.002

βp=2

βp=0.5 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0

1

0.2

0.4

x1

Fig. 7. Excess volume calculated from virial coefficients to various order at bp = 0.5. Dashed line: only the second order, dotted line: to the third order, dash-dotted line: to the fourth order, double dashed line: to the fifth order. Solid line with error bars is the fit to MC data (cf. Fig. 8).

with the increasing pressure, see Fig. 2. We can also mention that the number of calculations is rapidly growing with more elongated molecules with L* close to 1. There are necessary 94 iterations for L* = 0.3, 149 iterations for L* = 0.6 and 1756 iterations for L* = 1 in case of pure diatomics and reduced density equal to 0.4. The results from the OZ equation were compared with equations of state. We used two modern equations and one classic one but the results are very different from the first method, see Figs. 3 –5. The old Boublı´k’s EOS [10] gives consistently small and negative excess volumes while the new version [12] predicts positive values at higher pressures, so does the Largo and Solana EOS [11], see Fig. 6. An independent route to the excess volume, though limited to small pressures, is the virial expansion. Our MC results on the mixture virial coefficients are collected in Table 2. The predicted excess volumes for bp = 0.5, as calculated from the volume virial expansion (26) truncated at n = 2, 3, 4, 5, are shown in Fig. 7. They seem to slowly converge to values within standard errors the same as MC results shown in the same figure. The convergence of the virial expansion is much better for bp < 0.2 and becomes unusable for bp > 1, see Table 3. The low-density limit with the second virial only gives A = 2B 11  B 02 - B 20 =  0.06654(2). These excess volumes are by one order of magnitude greater than the excess volumes at bp å 1. We have two notes about the virial expansion. First, we have used the volume virial expansion (series in q). Using the Table 3 Excess volume expressed using parameter A (Eq. (19)) calculated from truncated virial expansion up to B n n

bp = 0.01

bp = 0.1

bp = 0.5

bp = 2

bp = 4

2 3 4 5

0.05234 0.05907 0.05932 0.05933

0.01111 0.02352 0.02718 0.02784

 0.00224 0.00214 0.00472 0.00577

0.002511 0.000900 0.000218 0.000609

 0.001968  0.000928  0.000134 0.000075

0.6

0.8

1

x1

Fig. 8. Raw and correlated MC data on the excess volume for bP = 2 and N = 500. Dots with error bars: MC data with estimated standard deviations, parabola with error bars: model (19) with standard deviation.

pressure-explicit virial expansion (series in bp) may be considered as a natural choice because it can be easily converted to an explicit (albeit rather complex) formula for the excess volume as a function of bp. However, the coefficients in this expansion are large in absolute value and have alternating signs. The expansion apparently diverges for bp > 0.05. Second, a widely used trick to improve the convergence of virial expansion is its resummation via Pade´ approximant, i.e., a ratio of two polynomials with the same density expansion. The two most natural choices for the expansion to q 5 are P(3, 2) (numerator to q 3, denominator to q 2) and P(2, 3); both give unusable results of the excess volume. Only some Pade´ approximants of more truncated expansion, for instance P(2, 2), give results comparable to the truncated virial series. The last method giving the excess volumes is the NPT MC simulation. The results for N = 500 and bp = 2 are shown in Fig. 8. It is seen that the accuracy of MC calculations is not sufficient to determine the skewness of the curve. Therefore in subsequent analysis of MC data we assumed that the excess volumes are given by formula (19). This is equivalent to fitting the series of MC data at the same bp and different compositions (including the pure compounds) to a quadratic function; the coefficient at x 12 gives directly parameter A. This fit, along with the error estimate (combined both from the MC-based error estimates of individual data and the error of

Table 4 Direct MC simulations of the excess volume given by parameter A, Eq. (19) N

bp

A

500 500 500 1000

4 2 0.5 4

0.0022(5) 0.0044(7) 0.0079(12) 0.0023(5)

N is the number of particles and bp is the reduced pressure. Estimated standard errors are given in parentheses.

28

P. Mora´vek et al. / Journal of Molecular Liquids 125 (2006) 22 – 28

the fit) is also shown in Fig. 8. Values of parameter A are given in Table 4. The simulation-based excess volumes are small, negative, and decreasing (in absolute value) as the pressure rises. 4. Conclusions In the first part of this work, a novel numerical procedure for fast solution of the Ornstein – Zernike equation was successfully implemented and tested. The solution is efficient and numerically stable even at high pressures. It is about fifteen times faster in comparison with direct iterations. Since it is based on a public domain software and its usage is simple and independent on the particular form of the equation, it can be recommended. In the second part, excess volumes of a mixture of homodumbbells were studied. The results are subject of discrepancies. MC simulations indicate that the excess volumes are small (compared to the low-pressure limit) and decreasing with increasing pressure. At low pressure the MC results agree with the volume-explicit virial expansion which cannot be however used at higher pressures. As regards the equations of state, the equation in best agreement with these finding is the 15 years old Boublı´k’s equation [10], other equations lead to values of wrong sign at high pressures. On the other hand, the Ornstein –Zernike equation results are not satisfactory because they are too large (in absolute value) and exhibit wrong pressure dependence. As the Ornstein– Zernike equation reproduces the third virial coefficient, it may be interesting to compare the results to the virial expansion truncated at n = 3. Here we find that the volumeexplicit virial expansion works better than the Ornstein –Zernike equation, but the pressure-explicit virial expansion much worse. This discrepancy may be then caused by neglected q 4 and higher terms (if not by errors in the q 2 and q 3 terms caused by insufficient number of harmonic expansion coefficients).

Acknowledgment This work was supported by the Ministry of Education, Youth and Sports of the Czech Republic under project LN00A032 (Structure and Dynamics of Complex Molecular Systems and Biomolecules). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24]

F. Lado, J. Chem. Phys. 47 (1967) 4828. J. Kolafa, Mol. Phys. 75 (1992) 577. M.J. Gillan, Mol. Phys. 38 (1979) 1781. S. Labı´k, A. Malijevsky´, P. Vonˇka, Mol. Phys. 56 (1985) 709. S. Labı´k, R. Pospı´xil, A. Malijevsky´, W.R. Smith, J. Comp. Phys. 115 (1994) 12. F. Lado, Mol. Phys. 47 (1982) 283. H.F. Walker, L. Zhou, Numer. Linear Algebra Appl. 1 (1994) 571. M.J. Booth, A.G. Schlijper, L.E. Scales, A.D.J. Haymet, Comp. Phys. Commun. 119 (1999) 122. M. Pernice, H.F. Walker, J. Sci. Comput. 19 (1998) 302. T. Boublı´k, Mol. Phys. 68 (1989) 191. J. Largo, J.R. Solana, Mol. Phys. 96 (1999) 1367. T. Boublı´k, Mol. Phys. 100 (2002) 3443. J.P. Hansen, I.R. McDonald, Theory of Simple Liquids, Academic Press, 1976. P.N. Brown, Y. Saad, J. Sci. Comput. 11 (1990) 450. J.K. Percus, G.J. Yevick, Phys. Rev. 110 (1958) 1. M.J. Van Leeuwen, J. Groenweld, J. De Boer, Physica 25 (1959) 729. G.A. Martynov, G. Sarkisov, Mol. Phys. 49 (1983) 1495. F.J. Rogers, D.A. Young, Phys. Rev., A 30 (1984) 999. D.M. Duh, D. Henderson, J. Chem. Phys. 104 (1996) 6742. S. Labı´k, A. Malijevsky´, W.R. Smith, Mol. Phys. 69 (1990) 649. M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Clarendon Press, Oxford, 1986. E. Enciso, N.G. Almarza, M.A. Gonza´lez, F.J. Bermejo, Mol. Phys. 100 (2002) 1941. J.L.F. Abascal, J. Lago, J. Mol. Liq. 30 (1985) 133. P. Mora´vek, Diploma Thesis, ICT Prague, 2002.