add closures and checks

Project.toml

@@ -4,9 +4,14 @@ authors = ["Ilian Pihlajamaa <[email protected]>"]
 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"
 

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
+
+
+

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))

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]

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)
+

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