Hard-sphere mixture

The previous example showed the case of a 1-component system. Let's instead look at mixtures here. Consider a 3:1 mixture of hard spheres with sizes 0.5 and 1.0. We solve the system with the same three steps as before.

First, we define the potential. Here, we must use HardSpheres, which takes a vector containing the diameters of each species and applies an additive mixing rule. For more information see HardSpheres.

using OrnsteinZernike
D = [0.5, 1.0]
potential = HardSpheres(D)

HardSpheres([0.5 0.75; 0.75 1.0])

Secondly, we define the system. In this example, the total number density is $\rho = 1.6$. For mixtures, the system expects a vector of individual densities. Those are computed by multiplying the total density with the concentration fractions.

ρ_total = 1.6
ρ = ρ_total*[0.75, 0.25] ## it is a 3:1 system
dims = 3 # we consider a 3D system
kBT = 1.0 # thermal energy
system = SimpleLiquid(dims, ρ, kBT, potential)

3 dimensional SimpleLiquid:
 ρ = [1.2000000000000002 0.0; 0.0 0.4]
 kBT = 1.0
 potential = HardSpheres{StaticArraysCore.SMatrix{2, 2, Float64, 4}}([0.5 0.75; 0.75 1.0])

Thirdly, we define the closure

closure = PercusYevick()

PercusYevick()

And now we solve the system.

sol = solve(system, closure)

OZSolution{Vector{Float64}, Array{Float64, 3}}:
 r = [0.005 0.015 … 9.985 9.995000000000001]'
 k = [0.15707963267948966 0.47123889803846897 … 313.6880264609408 314.00218572629984]'
 cr = [-4.907370880024967 -4.922251537243953; -4.828585215242844 -4.922236155034784; … ; 0.0 0.0; 0.0 0.0;;; -4.922251537243951 -16.552339427948283; -4.922236155034795 -16.35006544342643; … ; 0.0 0.0; 0.0 0.0]
 gr = [-3.766986722553156e-11 -8.103739901343943e-12; -3.6840308581531644e-11 -8.1574746957358e-12; … ; 1.0000000000179485 1.0000000000919336; 1.0000000000179385 1.0000000000921834;;; -8.109513061071993e-12 -5.509726008767757e-11; -8.167244658352502e-12 -5.3434590085998934e-11; … ; 1.0000000000919336 0.9999999999571715; 1.0000000000921834 0.9999999999578801]
 ck = [-1.3544005039517346 -5.733054985781598; -1.3483234766315801 -5.676806415958998; … ; 0.00018955687242657574 -0.0002998037504976306; 0.00019414850674488198 -0.0003170267388477564;;; -5.733054985781601 -21.713019192358065; -5.676806415959002 -21.37210713033137; … ; -0.00029980375049763085 0.0004441194941040604; -0.0003170267388477567 0.0004918339066380198]
 Sk = [0.7527530004381633 -0.1782336315822241; 0.7514659792847361 -0.17869921868599145; … ; 0.7501706723784687 -8.997757675530539e-5; 0.7501748105808242 -9.514890848941872e-5;;; -0.17823363158222422 0.06801395551266855; -0.1786992186859916 0.06867600175438471; … ; -8.997757675530547e-5 0.2500444306326604; -9.514890848941882e-5 0.25004920513686496]

For mixtures, the fields sol.gr, sol.cr, sol.ck, and sol.Sk are now three-dimensional arrays with shape (Nr, Ns, Ns). For example, $g_{12}(r_6)$ is stored in sol.gr[6,1,2].

We just solved the system using the default iterative solver NgIteration introduced by Ng. However, in this specific case, an exact solution is implemented. To use this, we specify the method Exact().

method = Exact()
sol_exact = solve(system, closure, method)

OZSolution{Vector{Float64}, Array{Float64, 3}}:
 r = [0.005 0.015 … 9.985 9.995000000000001]'
 k = [0.15707963267948966 0.47123889803846897 … 313.6880264609408 314.00218572629984]'
 cr = [-4.883337720022288 -4.921630466346176; -4.820547994642327 -4.921609849449453; … ; 0.0 0.0; 0.0 0.0;;; -4.921630466346171 -16.489170751785; -4.92160984944945 -16.327959345117435; … ; 0.0 0.0; 0.0 0.0]
 gr = [0.0 0.0; 0.0 0.0; … ; 1.0000000000179539 1.000000000091713; 1.0000000000179388 1.0000000000919926;;; 0.0 0.0; 0.0 0.0; … ; 1.000000000091713 0.9999999999570701; 1.0000000000919926 0.9999999999576932]
 ck = [-1.3542110635234053 -5.732061421503805; -1.3481337052508677 -5.675818584421081; … ; 0.00012083435991868302 -0.0001908870980911698; 0.0001236511084374651 -0.00020183403704973214;;; -5.732061421503805 -21.708321868707827; -5.675818584421081 -21.36747031400456; … ; -0.0001908870980911698 0.0002822216330594651; -0.00020183403704973214 0.00031293547727240245]
 Sk = [0.7527369373510386 -0.059411172495853294; 0.7514498662444521 -0.0595663338703082; … ; 0.7501087798181567 -1.9093633885573983e-5; 0.7501113171843109 -2.0188926517075634e-5;;; -0.5347005524626797 0.06801981064775306; -0.5360970048327746 0.06868189601974367; … ; -0.00017184270497016586 0.25002822972379546; -0.0001817003386536808 0.25003130235574933]

Let's plot the resulting $g(r)$.

using Plots
p = plot(sol.r, sol.gr[:, 1, 1], xlims=(0,5), xlabel="r", ylabel="g(r)", lw=4, label="g11(r) iterative")
plot!(sol.r, sol.gr[:, 1, 2], xlabel="r", ylabel="g(r)", lw=4, label="g12(r) iterative")
plot!(sol.r, sol.gr[:, 2, 2], xlabel="r", ylabel="g(r)", lw=4, label="g22(r) iterative")
plot!(sol_exact.r, sol_exact.gr[:, 1, 1], lw=2, c=:black, label="exact")
plot!(sol_exact.r, sol_exact.gr[:, 1, 2], lw=2, c=:black, label=nothing)
plot!(sol_exact.r, sol_exact.gr[:, 2, 2], lw=2, c=:black, label=nothing)
p