From 5e9d5cafa7d874a46b11b9986e39ddfdc1c6fb09 Mon Sep 17 00:00:00 2001 From: Ilian Pihlajamaa Date: Mon, 23 Oct 2023 09:32:28 +0200 Subject: [PATCH] add closures and checks --- Project.toml | 5 + checks/check_BBclosure.jl | 91 ++++++++++++ checks/check_CarbajalTinoko.jl | 34 +++++ checks/check_Lee.jl | 62 ++++++++ checks/check_compressibility.jl | 34 +++++ checks/check_differentiation.jl | 39 +++++ checks/check_self_consistency.jl | 242 +++++++++++++++++++++---------- src/Closures.jl | 140 ++++++++++++++++-- src/OrnsteinZernike.jl | 7 +- src/Thermodynamics.jl | 17 ++- 10 files changed, 570 insertions(+), 101 deletions(-) create mode 100644 checks/check_BBclosure.jl create mode 100644 checks/check_CarbajalTinoko.jl create mode 100644 checks/check_Lee.jl create mode 100644 checks/check_compressibility.jl create mode 100644 checks/check_differentiation.jl diff --git a/Project.toml b/Project.toml index 699ff60..158d8fd 100644 --- a/Project.toml +++ b/Project.toml @@ -4,9 +4,14 @@ authors = ["Ilian Pihlajamaa "] version = "0.1.1" [deps] +Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94" FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" Hankel = "74863788-d124-456e-a676-9b76578dd39e" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" +Optim = "429524aa-4258-5aef-a3af-852621145aeb" +Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" diff --git a/checks/check_BBclosure.jl b/checks/check_BBclosure.jl new file mode 100644 index 0000000..631de5a --- /dev/null +++ b/checks/check_BBclosure.jl @@ -0,0 +1,91 @@ +import Pkg; Pkg.activate(".") +using Revise +using OrnsteinZernike, Plots, Dierckx +import Roots, ForwardDiff + +function find_self_consistent_solution(ρ, kBT, M, dr, dims, pot) + + function pressure(ρ, α) + method = NgIteration(M=M, dr=dr, verbose=false) + system = SimpleLiquid(dims, ρ, kBT, pot) + sol = solve(system, BomontBretonnet(α), method) + p = compute_virial_pressure(sol, system) + return p + end + + function find_inconsistency(ρ, α) + system1 = SimpleLiquid(dims, ρ, kBT, pot) + method = NgIteration(M=M, dr=dr, verbose=false) + sol1 = solve(system1, BomontBretonnet(α), method) + + dpdρ = @time ForwardDiff.derivative(ρ -> pressure(ρ, α), ρ) + + @time begin + p1 = compute_virial_pressure(sol1, system1) + dρ = sqrt(eps(ρ)) + system2 = SimpleLiquid(dims, ρ+dρ, kBT, pot) + sol2 = solve(system2, BomontBretonnet(α), method) + p2 = compute_virial_pressure(sol2, system2) + dpdρ2 = (p2-p1)/dρ + end + @time begin + dρ = sqrt(eps(ρ)) + system2 = SimpleLiquid(dims, ρ+dρ, kBT, pot) + sol2 = solve(system2, BomontBretonnet(α), method) + p2 = compute_virial_pressure(sol2, system2) + system3 = SimpleLiquid(dims, ρ-dρ, kBT, pot) + sol3 = solve(system3, BomontBretonnet(α), method) + p3 = compute_virial_pressure(sol3, system3) + dpdρ3 = (p2-p3)/dρ/2 + end + @show dpdρ, dpdρ2, dpdρ3 + + χ = compute_compressibility(sol1, system1) + inconsistency = dpdρ/kBT - 1/(ρ*kBT*χ) + + return inconsistency + end + + func = α -> find_inconsistency(ρ, α) + α = Roots.find_zero(func, (0.0,1.0), Roots.Bisection(), atol=0.0001) + system = SimpleLiquid(dims, ρ, kBT, pot) + method = NgIteration(M=M, dr=dr, verbose=false) + sol = solve(system, BomontBretonnet(α), method) + return system, sol, α +end +println("Hard Spheres") +for ρstar = [0.3, 0.5, 0.7, 0.8, 0.9] + ρ = ρstar + M = 1000 + dr = 10.0/M + kBT = 1.0 + dims = 3 + pot = HardSpheres(1.0) + system, sol, α = find_self_consistent_solution(ρ, kBT, M, dr, dims, pot) + + P = compute_virial_pressure(sol, system)/ρ/kBT - 1 + gmax = maximum(sol.gr) + B = (OrnsteinZernike.bridge_function(BomontBretonnet(α), 0.0, 0.0, sol.gr .- 1 .- sol.cr)) + plot(sol.r, B) |> display + println("At ρ = $(ρstar), we find α = $(trunc(α,digits=4)), and βp/ρ - 1 = $(trunc(P,digits=4)).") + println("B(0) = $(B[1])") +end + +for closure in [PercusYevick, MartynovSarkisov] + for ρ in [0.3, 0.5, 0.7, 0.8, 0.9] + M = 1000000 + dr = 10.0/M + kBT = 1.0 + dims = 3 + pot = HardSpheres(1.0) + system = SimpleLiquid(dims, ρ, kBT, pot) + method = NgIteration(M=M, dr=dr, verbose=false) + sol = solve(system, closure(), method) + B = (OrnsteinZernike.bridge_function(closure(), 0.0, 0.0, sol.gr .- 1 .- sol.cr)) + + println("for closure $(closure) at ρ = $(ρ), we have B(0) = $(Spline1D(sol.r,B)(0.0))") + end +end + + + diff --git a/checks/check_CarbajalTinoko.jl b/checks/check_CarbajalTinoko.jl new file mode 100644 index 0000000..bf2b062 --- /dev/null +++ b/checks/check_CarbajalTinoko.jl @@ -0,0 +1,34 @@ +using OrnsteinZernike +M = 10000 +ρ = 0.9 +kBT = 1.0 +dims = 3 +dr = 10.0/M +pot = HardSpheres(1.0) +system = SimpleLiquid(dims, ρ, kBT, pot) +closure = CarbajalTinoko(0.3) +method = NgIteration() +sol = solve(system, closure, method) + +using Roots +λ = 0.3 +r, γ = 2.180762339605223, 0.0007099178871052923 +e = ifelse(λ > 0, 3 + λ, 3exp(λ*r)) + +f = function (b) + ω = γ + b + @show ω + if abs(ω) < 0.0001 + bfunc = e*(-(ω^2/6)+ω^4/360) + else + y = exp(ω) + bfunc = e*((2-ω)*y - 2 - ω)/(y - 1) + end + obj = b - bfunc + @show r, γ, obj + return obj +end +b = find_zero(f, γ) + +using Plots +plot(-10:0.1:10.0, f.(-10:0.1:10.0)) \ No newline at end of file diff --git a/checks/check_Lee.jl b/checks/check_Lee.jl new file mode 100644 index 0000000..b4e9edc --- /dev/null +++ b/checks/check_Lee.jl @@ -0,0 +1,62 @@ +using NLsolve, OrnsteinZernike + +function find_self_consistent_solution_Lee(ρ, kBT, M, dr, dims, pot) + function find_inconsistency(ρ, params) + ζ, ϕ, α = params + @show ζ, ϕ, α + closure = Lee(ζ, ϕ, α, ρ) + system1 = SimpleLiquid(dims, ρ, kBT, pot) + method = NgIteration(M=M, dr=dr, verbose=false) + sol1 = solve(system1, closure, method) + p1 = compute_virial_pressure(sol1, system1) + + dρ = sqrt(eps(ρ)) + system2 = SimpleLiquid(dims, ρ+dρ, kBT, pot) + sol2 = solve(system2, closure, method) + p2 = compute_virial_pressure(sol2, system2) + + dpdρ = (p2-p1)/dρ + + χ = compute_compressibility(sol1, system1) + inconsistency1 = dpdρ - 1/(ρ*χ) + + βu = OrnsteinZernike.evaluate_potential(system1.potential, sol1.r) + mayer_f = OrnsteinZernike.find_mayer_f_function.((system1,), βu) + + gammastar0 = sol1.gr[1] .- sol1.cr[1] .- 1 .+ ρ/2*mayer_f[1] + gamma0 = sol1.gr[1] .- sol1.cr[1] .- 1 + + η = ρ/6*π + grmax = (1-0.5η)/(1-η)^3 + + dBdgamma0 = -gammastar0*ζ - (α^2*gammastar0^3*ϕ*ζ)/(2*(1*+ α*gammastar0)^2) + (3*α*gammastar0^2*ϕ*ζ)/(2*(1 + α*gammastar0)) + inconsistency2 = dBdgamma0 - 1 + 1/grmax + + b = OrnsteinZernike.bridge_function(closure, sol1.r, mayer_f, sol1.gr .- sol1.cr .- 1.0) + mu = (8η-9η^2+3η^3)/(1-η)^3 + inconsistency3 = b[1] + gamma0 - mu + return [inconsistency1,inconsistency2, inconsistency3] + end + @show find_inconsistency(ρ, [1.2041, 0.9962, 1.0]) + func = α -> find_inconsistency(ρ, α) + s = nlsolve(func, [1.2041, 0.9962, 1.0], xtol=0.001, factor=0.1) + @show s + params = s.zero + ζ, ϕ, α = params + system = SimpleLiquid(dims, ρ, kBT, pot) + method = NgIteration(M=M, dr=dr, verbose=false) + closure = Lee(ζ, ϕ, α, ρ) + sol = solve(system, closure, method) + return system, sol, params +end + +ρ = 0.9 +sys, sol, params = find_self_consistent_solution_Lee(ρ, 1.0, 10^4, 0.001, 3, HardSpheres(1.0)) +ζ, ϕ, α= params +closure = Lee(ζ, ϕ, α, ρ) + +βu = OrnsteinZernike.evaluate_potential(sys.potential, sol.r) +mayer_f = OrnsteinZernike.find_mayer_f_function.((sys,), βu) + +br = OrnsteinZernike.bridge_function(closure, sol.r, mayer_f, sol.gr - sol.cr .- 1) +@show br[1] \ No newline at end of file diff --git a/checks/check_compressibility.jl b/checks/check_compressibility.jl new file mode 100644 index 0000000..c3abead --- /dev/null +++ b/checks/check_compressibility.jl @@ -0,0 +1,34 @@ +import Pkg; Pkg.activate(".") +using OrnsteinZernike, Plots + +# Make sure the discontinuity is a multiple of dr +Rmax = 20.0 +Ms = Rmax * round.(Int, 10 .^ (range(1,4,length=30))) + +χ1 = zeros(length(Ms)) +χ2 = zeros(length(Ms)) +ρ = 1.0 +kBT = 1.0 +dims = 3 + +pot = HardSpheres(1.0) +system = SimpleLiquid(dims, ρ, kBT, pot) +for (i,M) in enumerate(Ms) + @show i, M + dr = Rmax/M + method = NgIteration(M=M, dr=dr, verbose=false) + sol1 = solve(system, PercusYevick(), method) + χ1[i] = compute_compressibility(sol1, system) + @show χ1[i] + sol2 = solve(system, PercusYevick(), Exact(M=M, dr=dr)) + χ2[i] = compute_compressibility(sol2, system) +end +η = ρ/6*π +crint = ((-4 + η)* (2 + η^2))/(24(-1 + η)^4) +χ = 1/(1-ρ*4π*crint)/ρ/kBT + +@show χexact = ((1+2η)^-2 * (1-η)^4)/ρ/kBT +pl = scatter(Ms, abs.(χ1.-χexact)./χexact, label="Iterative") +plot!(Ms, abs.(χ1.-χexact)./χexact, label="from exact c(k)") +plot!(ylabel="log10(relative error)", xlabel="log10(M)", xscale=:log, yscale=:log) + diff --git a/checks/check_differentiation.jl b/checks/check_differentiation.jl new file mode 100644 index 0000000..5ab6300 --- /dev/null +++ b/checks/check_differentiation.jl @@ -0,0 +1,39 @@ +import Pkg; Pkg.activate(".") +using Revise +using OrnsteinZernike, Plots, Dierckx +import Roots, ForwardDiff + +function find_pressure_derivative(ρ, kBT, dims, pot, closure, method) + + function pressure(ρ) + system = SimpleLiquid(dims, ρ, kBT, pot) + sol = solve(system, closure, method) + p = compute_virial_pressure(sol, system) + return p + end + + dρ = sqrt(eps(ρ)) + + # @time dpdρ = ForwardDiff.derivative(pressure, ρ) + @time dpdρ2 = (pressure(ρ+dρ) - pressure(ρ))/(dρ) + @time dpdρ3 = (pressure(ρ+dρ) - pressure(ρ-dρ))/(2dρ) + @time dpdρ4 = (3kBT* (-72+π*ρ* (-60+π *ρ* (-18+π *ρ))))/(-6+π* ρ)^3 + @show dpdρ2, dpdρ3, dpdρ4 + + return +end + + + +println("Hard Spheres") +@profview for ρstar = [0.3] + ρ = ρstar + M = 100000 + dr = 10.0/M + kBT = 1.0 + dims = 3 + pot = HardSpheres(1.0) + closure = PercusYevick() + method = NgIteration(M=M, dr=dr, verbose=false) + find_pressure_derivative(ρ, kBT, dims, pot, closure, method) +end diff --git a/checks/check_self_consistency.jl b/checks/check_self_consistency.jl index 14543b2..acc9d04 100644 --- a/checks/check_self_consistency.jl +++ b/checks/check_self_consistency.jl @@ -3,101 +3,185 @@ using Revise using OrnsteinZernike, Plots import Roots -function find_self_consistent_solution(ρ, kBT, M, dr, dims, pot) +# function find_self_consistent_solution(ρ, kBT, method, dims, pot; lims=(0.1, 2.0)) - function RY_inconsistency(ρ, α) - system1 = SimpleLiquid(dims, ρ, kBT, pot) - method = NgIteration(M=M, dr=dr, verbose=false) - sol1 = solve(system1, RogersYoung(α), method) - p1 = compute_virial_pressure(sol1, system1) +# function RY_inconsistency(ρ, α) +# system1 = SimpleLiquid(dims, ρ, kBT, pot) +# sol1 = solve(system1, RogersYoung(α), method) +# p1 = compute_virial_pressure(sol1, system1) + +# dρ = sqrt(eps(ρ)) +# system2 = SimpleLiquid(dims, ρ+dρ, kBT, pot) +# sol2 = solve(system2, RogersYoung(α), method) +# p2 = compute_virial_pressure(sol2, system2) +# dpdρ = (p2-p1)/dρ + +# χ = compute_compressibility(sol1, system1) +# inconsistency = dpdρ/kBT - 1/(ρ*kBT*χ) +# return inconsistency +# end + +# func = α -> RY_inconsistency(ρ, α) +# α = Roots.find_zero(func, lims, Roots.Bisection(), atol=0.0001) +# system = SimpleLiquid(dims, ρ, kBT, pot) +# sol = solve(system, RogersYoung(α), method) +# return system, sol, α +# end +# println("Hard Spheres") +# for ρstar = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.654] +# ρ = ρstar*sqrt(2) +# M = 1000 +# dr = 10.0/M +# kBT = 1.0 +# dims = 3 +# pot = HardSpheres(1.0) +# method = NgIteration(dr=dr, M=M, verbose=false) +# system, sol, α = find_self_consistent_solution(ρ, kBT, method, dims, pot) +# P = compute_virial_pressure(sol, system)/ρ/kBT - 1 +# gmax = maximum(sol.gr) +# println("At ρ/√2 = $(ρstar), we find α = $(trunc(α,digits=2)), and βp/ρ - 1 = $(trunc(P,digits=4)).") +# end + +# ## 12 power fluid +# println("1/r^12 fluid") + +# for z = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.813] +# ρ = z*sqrt(2) +# M = 1000 +# dr = 20.0/M +# kBT = 1.0 +# dims = 3 +# Γ = (4π*sqrt(2)*z/3)^(4) +# σ = 1.0 +# ϵ = 1.0 + +# pot = PowerLaw(ϵ, σ, 12) +# method = NgIteration(dr=dr, M=M, verbose=false) +# system, sol, α = find_self_consistent_solution(ρ, kBT, method, dims, pot) +# P = compute_virial_pressure(sol, system)/ρ/kBT - 1 +# gmax = maximum(sol.gr) +# println("At z = $(z), we find Γ = $(trunc(Γ,digits=2)) α = $(trunc(α,digits=2)), and βp/ρ - 1 = $(trunc(P,digits=4)).") +# end + +# println("1/r^9 fluid") + +# for z = [0.1, 0.25, 0.5, 0.943] +# ρ = z*sqrt(2) +# M = 1000 +# dr = 20.0/M +# kBT = 1.0 +# dims = 3 +# Γ = (4π*sqrt(2)*z/3)^(3) +# σ = 1.0 +# ϵ = 1.0 + +# pot = PowerLaw(ϵ, σ, 9) +# method = NgIteration(dr=dr, M=M, verbose=false) +# system, sol, α = find_self_consistent_solution(ρ, kBT, method, dims, pot) +# P = compute_virial_pressure(sol, system)/ρ/kBT - 1 +# gmax = maximum(sol.gr) +# println("At z = $(z), we find Γ = $(trunc(Γ,digits=2)) α = $(trunc(α,digits=2)), and βp/ρ - 1 = $(trunc(P,digits=4)).") +# end + +# println("1/r^6 fluid") - dρ = sqrt(eps(ρ)) +# for z = [0.1, 0.25, 0.5, 1.0, 1.54] +# ρ = z*sqrt(2) +# M = 2000 +# dr = 20/M +# kBT = 1.0 +# dims = 3 +# Γ = (4π*sqrt(2)*z/3)^(2) +# σ = 1.0 +# ϵ = 1.0 + +# pot = PowerLaw(ϵ, σ, 6) +# method = NgIteration(dr=dr, M=M, verbose=false) +# system, sol, α = find_self_consistent_solution(ρ, kBT, method, dims, pot) +# P = compute_virial_pressure(sol, system)/ρ/kBT - 1 +# gmax = maximum(sol.gr) +# println("At z = $(z), we find Γ = $(trunc(Γ,digits=2)) α = $(trunc(α,digits=2)), and βp/ρ - 1 = $(trunc(P,digits=4)).") +# end + + +# println("Yukawa") +# #J. Chem. Phys. 128, 184507 共2008兲 +# for ϕ = [0.2, 0.3, 0.35, 0.4, 0.45, 0.5] +# ρ = ϕ*6/pi +# M = 10000 +# dr = 10.0/M +# kBT = 1.0 +# dims = 3 +# method = NgIteration(dr = dr, M = M, verbose=false, max_iterations=10^6, N_stages=2) +# pot = CustomPotential((r,p)->500*exp(-4r)/r, nothing) +# system, sol, α = find_self_consistent_solution(ρ, kBT, method, dims, pot; lims=(0.9, 2.0)) +# P = compute_virial_pressure(sol, system)/ρ/kBT - 1 +# println("At ϕ = $(ϕ), α = $(trunc(α,digits=4)), and βp/ρ = $(trunc(P,digits=4)).") +# end + + +using Optim +function find_self_consistent_solution_ERY(ρ, kBT, method, dims, pot; x0=[0.2, 0.2]) + function ERY_inconsistency(ρ, α2, a2) + α = sqrt(α2) + a = sqrt(a2) + closure = ExtendedRogersYoung(α, a) + dρ = ρ*0.001 + dkBT = sqrt(eps(kBT)) + + system0 = SimpleLiquid(dims, ρ-dρ, kBT, pot) + sol0 = solve(system0, closure, method) + system1 = SimpleLiquid(dims, ρ, kBT, pot) + sol1 = solve(system1, closure, method) system2 = SimpleLiquid(dims, ρ+dρ, kBT, pot) - sol2 = solve(system2, RogersYoung(α), method) + sol2 = solve(system2, closure, method) + system3 = SimpleLiquid(dims, ρ, kBT+dkBT, pot) + sol3 = solve(system3, closure, method) + + p0 = compute_virial_pressure(sol0, system0) p2 = compute_virial_pressure(sol2, system2) - dpdρ = (p2-p1)/dρ - χ = compute_compressibility(sol1, system1) - inconsistency = dpdρ/kBT - 1/(ρ*kBT*χ) + ρU0 = (ρ-dρ)*compute_excess_energy(sol0, system0) + ρU1 = ρ*compute_excess_energy(sol1, system1) + ρU2 = (ρ+dρ)*compute_excess_energy(sol2, system2) - return inconsistency - end - func = α -> RY_inconsistency(ρ, α) - α = Roots.find_zero(func, (0.1,5.0), Roots.Bisection(), atol=0.0001) - system = SimpleLiquid(dims, ρ, kBT, pot) - method = NgIteration(M=M, dr=dr, verbose=false) - sol = solve(system, RogersYoung(α), method) - return system, sol, α -end -println("Hard Spheres") -for ρstar = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.654] - ρ = ρstar*sqrt(2) - M = 1000 - dr = 10.0/M - kBT = 1.0 - dims = 3 - pot = HardSpheres(1.0) - system, sol, α = find_self_consistent_solution(ρ, kBT, M, dr, dims, pot) - P = compute_virial_pressure(sol, system)/ρ/kBT - 1 - gmax = maximum(sol.gr) - println("At ρ/√2 = $(ρstar), we find α = $(trunc(α,digits=2)), and βp/ρ - 1 = $(trunc(P,digits=4)).") -end + dpdρ = (p2-p0)/(2dρ) + d2ρUdρ2 = (ρU2 + ρU0 - 2ρU1)/(dρ^2) -## 12 power fluid -println("1/r^12 fluid") + ĉ0_1 = OrnsteinZernike.get_ĉ0(sol1, system1) + ĉ0_3 = OrnsteinZernike.get_ĉ0(sol3, system3) + dĉ0dkBT = (ĉ0_3-ĉ0_1)/dkBT + dĉ0dβ = - kBT^2 * dĉ0dkBT -for z = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.813] - ρ = z*sqrt(2) - M = 1000 - dr = 20.0/M - kBT = 1.0 - dims = 3 - Γ = (4π*sqrt(2)*z/3)^(4) - σ = 1.0 - ϵ = 1.0 - - pot = PowerLaw(ϵ, σ, 12) - system, sol, α = find_self_consistent_solution(ρ, kBT, M, dr, dims, pot) - P = compute_virial_pressure(sol, system)/ρ/kBT - 1 - gmax = maximum(sol.gr) - println("At z = $(z), we find Γ = $(trunc(Γ,digits=2)) α = $(trunc(α,digits=2)), and βp/ρ - 1 = $(trunc(P,digits=4)).") -end + χ = compute_compressibility(sol1, system1) + inconsistency1 = dpdρ/kBT - 1/(ρ*kBT*χ) -println("1/r^9 fluid") + inconsistency2 = d2ρUdρ2 + dĉ0dβ -for z = [0.1, 0.25, 0.5, 0.943] - ρ = z*sqrt(2) - M = 1000 - dr = 20.0/M - kBT = 1.0 - dims = 3 - Γ = (4π*sqrt(2)*z/3)^(3) - σ = 1.0 - ϵ = 1.0 + inconsistency = sum([inconsistency1, inconsistency2].^2) + @show inconsistency, α, a, d2ρUdρ2 + return inconsistency + end - pot = PowerLaw(ϵ, σ, 9) - system, sol, α = find_self_consistent_solution(ρ, kBT, M, dr, dims, pot) - P = compute_virial_pressure(sol, system)/ρ/kBT - 1 - gmax = maximum(sol.gr) - println("At z = $(z), we find Γ = $(trunc(Γ,digits=2)) α = $(trunc(α,digits=2)), and βp/ρ - 1 = $(trunc(P,digits=4)).") + func = x -> ERY_inconsistency(ρ, x[1]^2, x[2]^2) + α = optimize(func, [x0...]+rand(2)./10, NelderMead(), Optim.Options(x_tol=0.0001, f_tol=10^-6)).minimizer + system = SimpleLiquid(dims, ρ, kBT, pot) + sol = solve(system, ExtendedRogersYoung(α[1], α[2]), method) + return system, sol, α end -println("1/r^6 fluid") -for z = [0.1, 0.25, 0.5, 1.0, 1.54] - ρ = z*sqrt(2) - M = 2000 - dr = 20/M +x0s = [(0.323, 0.177), (0.375, 0.213), (0.398, 0.231), (0.423, 0.243), (0.446, 0.255), (0.467, 0.270)] +for (i,ϕ) = enumerate([0.2, 0.3, 0.35, 0.4, 0.45, 0.5]) + ρ = ϕ*6/pi + M = 10000 + dr = 10.0/M kBT = 1.0 dims = 3 - Γ = (4π*sqrt(2)*z/3)^(2) - σ = 1.0 - ϵ = 1.0 - - pot = PowerLaw(ϵ, σ, 6) - system, sol, α = find_self_consistent_solution(ρ, kBT, M, dr, dims, pot) + method = NgIteration(dr = dr, M = M, verbose=false, max_iterations=10^6, N_stages=2) + pot = CustomPotential((r,p)->500*exp(-4r)/r, nothing) + system, sol, α = find_self_consistent_solution_ERY(ρ, kBT, method, dims, pot; x0=x0s[i]) P = compute_virial_pressure(sol, system)/ρ/kBT - 1 - gmax = maximum(sol.gr) - println("At z = $(z), we find Γ = $(trunc(Γ,digits=2)) α = $(trunc(α,digits=2)), and βp/ρ - 1 = $(trunc(P,digits=4)).") + println("At ϕ = $(ϕ), (α,a) = $(trunc(α[1],digits=4)), $(trunc(α[2],digits=4)), and βp/ρ = $(trunc(P,digits=4)).") end \ No newline at end of file diff --git a/src/Closures.jl b/src/Closures.jl index 47a80e1..f6685aa 100644 --- a/src/Closures.jl +++ b/src/Closures.jl @@ -106,6 +106,32 @@ function bridge_function(closure::Verlet, _, _, γ) return @. -(A*γ^2/2)/(oneunit + B*γ/2) end +""" + ModifiedVerlet <: Closure + +Implements the modified Verlet Closure \$b(r) = -\\frac{\\gamma^2(r)/2}{1+\\alpha \\gamma(r)}\$. If \$\\gamma(r)<0\$, the closure reads \$-\\gamma^2/2\$. +Example: +```julia +closure = ModifiedVerlet(0.2) +``` +References: + +Verlet, Loup. "Integral equations for classical fluids: I. The hard sphere case." Molecular Physics 41.1 (1980): 183-190. + +""" +struct ModifiedVerlet{T} <: Closure + α::T +end + + +function bridge_function(closure::ModifiedVerlet, _, _, γ) + α = closure.α + oneunit = one.(γ) + B = @. ifelse(γ < 0, -(γ^2/2), -(γ^2/2)/(oneunit + α*γ/2)) + return B +end + + """ MartynovSarkisov <: Closure @@ -179,6 +205,51 @@ function bridge_function(closure::RogersYoung, r, _, γ) return b end +function closure_c_from_gamma(closure::RogersYoung, r, mayer_f, γ, _) + oneunit = one.(γ) + α = closure.α + @assert α > 0 + f = @. 1.0 - exp(-α*r) + term = @. (exp(f*γ)-oneunit)/f + + return @. (mayer_f + oneunit)*(oneunit + term) - γ - oneunit +end + + +""" + ExtendedRogersYoung <: Closure + +Implements the extended Rogers-Young closure \$b(r) = \\ln(a\\phi(r)^2 + \\phi(r) + 1) - γ(r) \$. +Here, \$\\phi(r) = \\frac{\\exp(f(r)\\gamma(r)) - 1}{f(r)}\$, and \$f(r)=1-\\exp(-\\alpha r)\$, in which \$\\alpha\$ is a free parameter, +that may be chosen such that thermodynamic consistency is achieved. +Example: +```julia +closure = ExtendedRogersYoung(0.5, 0.5) # order is α, a +``` + +References: +J. Chem. Phys. 128, 184507 (2008) + +""" +struct ExtendedRogersYoung{T} <: Closure + α::T + a::T +end + +function bridge_function(closure::ExtendedRogersYoung, r, _, γ) + oneunit = one.(γ) + α = closure.α + a = closure.a + @assert α > 0 + f = @. 1.0 - exp(-α*r) + ϕ = @. (exp(f*γ)-oneunit)/f + b = @. -γ + log1p(ϕ + a*ϕ^2) + return b +end + + + + """ ZerahHansen <: Closure @@ -225,7 +296,7 @@ References: struct DuhHaymet <: Closure end function bridge_function(::DuhHaymet, _, _, γstar) - oneunit = one.(γ) + oneunit = one.(γstar) b = @. ifelse(γstar>0, (-γstar^2) / (2 * (oneunit + (5γstar+11oneunit)/(7γstar+9oneunit) * γstar)), -γstar^2/2) return b end @@ -247,16 +318,17 @@ References: Lee, Lloyd L. "An accurate integral equation theory for hard spheres: Role of the zero‐separation theorems in the closure relation." The Journal of chemical physics 103.21 (1995): 9388-9396. """ -struct Lee{T} <: Closure +struct Lee{T, T2} <: Closure ζ::T ϕ::T α::T + ρ::T2 end function bridge_function(closure::Lee, _, mayerf, γ) oneunit = one.(γ) ρ = closure.ρ - γstar = @. γ - ρ*mayerf/2 + γstar = @. γ + ρ*mayerf/2 ζ = closure.ζ α = closure.α ϕ = closure.ϕ @@ -285,7 +357,7 @@ struct ChoudhuryGhosh{T} <: Closure end function bridge_function(closure::ChoudhuryGhosh, _, _, γstar) - oneunit = one.(γ) + oneunit = one.(γstar) α = closure.α return @. ifelse(γstar>0, (-γstar^2) / (2 * (oneunit + α * γstar)), -γstar^2/2) end @@ -358,7 +430,7 @@ struct CharpentierJackse{T} <: Closure end function bridge_function(closure::CharpentierJackse, _, _, γstar) - oneunit = one.(γ) + oneunit = one.(γstar) α = closure.α return @. (sqrt(oneunit+4α*γstar) - oneunit - 2α*γstar)/(2α) end @@ -383,7 +455,7 @@ struct BomontBretonnet{T} <: Closure end function bridge_function(closure::BomontBretonnet, _, _, γstar) - oneunit = one.(γ) + oneunit = one.(γstar) f = closure.f return @. sqrt(oneunit + 2γstar + f*γstar^2) - oneunit - γstar end @@ -414,9 +486,9 @@ function bridge_function(closure::Khanpour, _, _, γ) end """ - ReferenceHypernettedChain <: Closure +ModifiedHypernettedChain <: Closure -Implements the Reference Hypernetted Chain closure \$b(r) = b_{HS}(r) \$. Here \$b_{HS}(r)=\\left((a_1+a_2x)(x-a_3)(x-a_4)/(a_3 a_4)\\right)^2\$ for \$x 0, 3 + λ, 3exp(λ*r)) + + f = function (b) + ω = γ + b + if abs(ω) < 0.0001 + bfunc = e*(-(ω^2/6)+ω^4/360) + else + y = exp(ω) + bfunc = e*((2-ω)*y - 2 - ω)/(y - 1) + end + obj = b - bfunc + return obj + end + b = find_zero(f, γ) + + return b +end \ No newline at end of file diff --git a/src/OrnsteinZernike.jl b/src/OrnsteinZernike.jl index df92f05..8228a71 100644 --- a/src/OrnsteinZernike.jl +++ b/src/OrnsteinZernike.jl @@ -4,15 +4,16 @@ A generic solver package for Ornstein-Zernike equations from liquid state theory """ module OrnsteinZernike -using FFTW, StaticArrays, LinearAlgebra, Hankel, SpecialFunctions +using FFTW, StaticArrays, LinearAlgebra, Hankel, SpecialFunctions, Dierckx +using Roots: find_zero export solve export SimpleLiquid export OZSolution export Exact, FourierIteration, NgIteration, DensityRamp -export PercusYevick, HypernettedChain, MeanSpherical, ReferenceHypernettedChain, Verlet, MartynovSarkisov +export PercusYevick, HypernettedChain, MeanSpherical, ModifiedHypernettedChain, Verlet, MartynovSarkisov export SoftCoreMeanSpherical, RogersYoung, ZerahHansen, DuhHaymet, Lee, ChoudhuryGhosh, BallonePastoreGalliGazzillo -export VompeMartynov, CharpentierJackse, BomontBretonnet, Khanpour +export VompeMartynov, CharpentierJackse, BomontBretonnet, Khanpour, ModifiedVerlet, CarbajalTinoko, ExtendedRogersYoung export CustomPotential, PowerLaw, HardSpheres, LennardJones export compute_compressibility, compute_excess_energy, compute_virial_pressure export WCADivision diff --git a/src/Thermodynamics.jl b/src/Thermodynamics.jl index b1e2e8f..2ff98a9 100644 --- a/src/Thermodynamics.jl +++ b/src/Thermodynamics.jl @@ -191,6 +191,8 @@ function get_concentration_fraction(system) error("Unreachable code: file an issue!") end +_eachslice(a;dims=1) = eachslice(a,dims=dims) +_eachslice(a::Vector;dims=1) = a """ @@ -207,12 +209,17 @@ function compute_compressibility(sol::OZSolution, system::SimpleLiquid{dims, spe kBT = system.kBT x = get_concentration_fraction(system) ρ0 = sum(ρ) - ĉ = sol.ck T = typeof(ρ0) - k = sol.k - dcdk = (ĉ[2, :, :]-ĉ[1, :, :]) ./ (k[2]-k[1]) - ĉ0 = ĉ[1, :, :] - k[1] * dcdk + ĉ0 = get_ĉ0(sol, system) invρkBTχ = one(T) - ρ0 * sum((x*x') .* ĉ0) χ = (one(T)/invρkBTχ)/(kBT*ρ0) return χ -end \ No newline at end of file +end + +function get_ĉ0(sol::OZSolution, ::SimpleLiquid{dims, species, T1, T2, P}) where {dims, species, T1, T2, P} + spl = Spline1D(sol.r, _eachslice(sol.cr, dims=1).*sol.r.^(dims-1)) + ĉ0 = surface_N_sphere(dims)*integrate(spl, zero(eltype(sol.r)), maximum(sol.r)) + return ĉ0 +end + +# function pressure_inconsystency(system, closure)