Merge branch 'main' of https://github.com/IlianPihlajamaa/OrnsteinZer…

…nike.jl

Project.toml

@@ -1,24 +1,29 @@
 name = "OrnsteinZernike"
 uuid = "51399321-c693-48c0-994d-05e3df0250f2"
 authors = ["Ilian Pihlajamaa <[email protected]>"]
-version = "0.1.3"
+version = "0.1.4"
 
 [deps]
+Bessels = "0e736298-9ec6-45e8-9647-e4fc86a2fe38"
 Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+FunctionZeros = "b21f74c0-b399-568f-9643-d20f4fa2c814"
 Hankel = "74863788-d124-456e-a676-9b76578dd39e"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"
 
 [compat]
+Bessels = "0.2"
+Dierckx = "0.5"
 FFTW = "1"
+ForwardDiff = "0.10"
+FunctionZeros = "0.3"
 Hankel = "0.5"
+Roots = "2.0"
 SpecialFunctions = "2"
 StaticArrays = "1"
 julia = "1.7"
-Roots = "2.0"
-ForwardDiff = "0.10"
-Dierckx = "0.5"

README.md

@@ -2,7 +2,7 @@
 A generic solver for Ornstein-Zernike equations from liquid state theory
 
 [![Build status (Github Actions)](https://github.com/IlianPihlajamaa/OrnsteinZernike.jl/workflows/CI/badge.svg)](https://github.com/IlianPihlajamaa/OrnsteinZernike.jl/actions)
-[![codecov.io](http://codecov.io/github/IlianPihlajamaa/OrnsteinZernike.jl/coverage.svg?branch=main)](http://codecov.io/github/IlianPihlajamaa/OrnsteinZernike.jl?branch=main)
+[![codecov](https://codecov.io/github/IlianPihlajamaa/OrnsteinZernike.jl/graph/badge.svg?token=ZFM4LKORIS)](https://codecov.io/github/IlianPihlajamaa/OrnsteinZernike.jl)
 [![](https://img.shields.io/badge/docs-stable-blue.svg)](https://IlianPihlajamaa.github.io/OrnsteinZernike.jl/stable)
 [![](https://img.shields.io/badge/docs-dev-blue.svg)](https://IlianPihlajamaa.github.io/OrnsteinZernike.jl/dev)
 

src/Closures.jl

@@ -30,9 +30,13 @@ closure = PercusYevick()
 struct PercusYevick <: Closure end
 
 function closure_cmulr_from_gammamulr(::PercusYevick, r::Number, mayer_f::T, Γmulr::T, βuLR) where T
-    return  @. r*(mayer_f+1)*exp(βuLR)*(1+Γmulr/r-βuLR) - Γmulr - r#mayer_f*(r + Γmulr)
+    return  @. (mayer_f+1)*exp(βuLR)*(r+Γmulr-βuLR*r) - Γmulr - r#mayer_f*(r + Γmulr)
 end
 
+# function closure_c_from_gamma(::PercusYevick, r::Number, mayer_f::T, Γmulr::T, βuLR) where T
+#     return  @. mayer_f*(r + Γmulr)
+# end
+
 function bridge_function(::PercusYevick, _, _, γ)
     B  = @. log1p(γ) - γ
     return B

src/FourierTransforms.jl

@@ -15,6 +15,16 @@ struct My1DPlan{T2,T3,T4, T5, T6}
     M::T4
 end
 
+struct MyNDPlan{dims, T, T2}
+    r::Vector{T}
+    k::Vector{T}
+    M::Int
+    plan_forward::Matrix{T2}
+    plan_backward::Matrix{T2}
+end
+
+
+
 function get_fourier_plan(::SimpleLiquid{1, species, T1, T2, P}, method, F) where {species, T1, T2, P}
     M = length(F)
     dr = method.dr
@@ -65,6 +75,27 @@ function find_fourier_plan_nd(::SimpleLiquid{ndims, species, T1, T2, P}, F::Vect
     M = length(F)
     plan = QDHT(ndims/2-1, 1, M*dr, M)
     return plan
+    # M = length(F)
+    # plan_forward = zeros(eltype(eltype(F)), M, M)
+    # plan_backward = zeros(eltype(eltype(F)), M, M)
+    # println("hi")
+    # p = ndims/2 - 1
+    # λ_i = besselj_zero.(p, 1:(M+1))
+    # Rmax = dr * M
+    # Kmax = λ_i[end]/Rmax
+    # k = λ_i[1:end-1]/Rmax
+    # r = λ_i[1:end-1]/Kmax
+
+    # for i = 1:M
+    #     for j = 1:M
+    #         plan_forward[j, i] = 2 * (2π)^(p+1) / Kmax^2 *  besselj(p, k[j]*r[i])/besselj(p+1, Kmax*r[i])^2
+    #         plan_backward[j, i] = 2 / ( Rmax^2 * (2π)^(p+1) ) * besselj(p, k[i]*r[j])/besselj(p+1, Rmax*k[i])^2
+    #     end
+    # end
+    # Tr = eltype(r)
+    # Tq = eltype(plan_backward)
+    # return MyNDPlan{ndims, Tr, Tq}(r, k, M, plan_forward, plan_backward)
+
 end
 
 function inverse_radial_fourier_transform_3d(F̂, r, k)
@@ -181,6 +212,17 @@ function fourier!(F̂::AbstractVector{T}, F::AbstractVector{T}, Q::Hankel.QDHT{p
     @. F = F / r^(halfdims-2)
 
     @. F̂ = F̂ * (2π)^halfdims / k^(halfdims-2)
+
+    # p = dims/2 - 1
+    # k = Q.k
+    # r = Q.r
+
+    # @. F = F * r^(p-1)
+    # mul!(F̂, Q.plan_forward, F)
+    # @. F = F / r^(p-1)
+
+    # @. F̂ = F̂ / k^(p-1)
+
 end
 
 """
@@ -199,6 +241,15 @@ function inverse_fourier!(F::AbstractVector{T}, F̂::AbstractVector{T}, Q::Hanke
     @. F̂ = F̂ / k^(halfdims-2)
     @. F ./= Q.scaleRK^2
     @. F = F * (2π)^(-halfdims) / r^(halfdims-2)
+
+    # p = dims/2 - 1
+    # k = Q.k
+    # r = Q.r
+
+    # @. F̂ = F̂ * k^(p-1)
+    # mul!(F, Q.plan_backward, F̂)
+    # @. F̂ = F̂ / k^(p-1)
+    # @. F = F / r^(p-1)
 end
 
 

src/OrnsteinZernike.jl

@@ -5,6 +5,8 @@ A generic solver package for Ornstein-Zernike equations from liquid state theory
 """
 module OrnsteinZernike
 using FFTW, StaticArrays, LinearAlgebra, Hankel, SpecialFunctions, Dierckx 
+using Bessels: besselj
+using FunctionZeros
 using Roots: find_zero
 
 export solve

src/Solvers.jl

@@ -19,7 +19,7 @@ Right now, the implemented exact methods are
 
 Construct using 
 
-```Exact(;  M=2^10, dr = sqrt(π/(M+1))/(2π))```
+```Exact(;  M=1000, dr = 10.0/M)```
 
 Here, `M` is the number of points that the exact solution is evaluated on, and `dr` is the grid spacing.
 These are used to perform Fourier transformations.
@@ -42,7 +42,7 @@ struct Exact <: Method
     dr::Float64
 end
 
-function Exact(;  M=2^10, dr = sqrt(π/(M+1))/(2π))
+function Exact(;  M=1000, dr = 10.0/M)
     return Exact(M, dr)
 end
 
@@ -68,7 +68,7 @@ Arguments:
 - `tolerance::Float64`: tolerance to be reached
 - `verbose::Bool`: whether or not to print convergence information
 
-Default: `FourierIteration(; mixing_parameter=0.5, max_iterations=10^5, tolerance=10^-6, verbose=true, M=2^10, dr=sqrt(π/(M+1))/(2π))`
+Default: `FourierIteration(; mixing_parameter=0.5, max_iterations=10^5, tolerance=10^-6, verbose=true, M=1000, dr=10.0/M)`
 """
 struct FourierIteration <: Method 
     M::Int
@@ -79,7 +79,7 @@ struct FourierIteration <: Method
     verbose::Bool
 end
 
-function FourierIteration(; mixing_parameter=0.5, max_iterations=10^5, tolerance=10^-10, verbose=true, M=2^10, dr=sqrt(π/(M+1))/(2π))
+function FourierIteration(; mixing_parameter=0.5, max_iterations=10^5, tolerance=10^-10, verbose=true, M=1000, dr=10.0/M)
     @assert max_iterations > 0 
     @assert tolerance > 0 
     @assert 0 <= mixing_parameter <= 1
@@ -108,7 +108,7 @@ Arguments:
 - `tolerance::Float64`: tolerance to be reached
 - `verbose::Bool`: whether or not to print convergence information
 
-Default: `NgIteration(; N_stages=3, max_iterations=10^3, tolerance=10^-6, verbose=true, M=2^10, dr=sqrt(π/(M+1))/(2π))`
+Default: `NgIteration(; N_stages=3, max_iterations=10^3, tolerance=10^-6, verbose=true, M=1000, dr=10.0/M)`
 
 References:
 Ng, K. C. (1974). Hypernetted chain solutions for the classical one‐component plasma up to Γ= 7000. The Journal of Chemical Physics, 61(7), 2680-2689.
@@ -122,7 +122,7 @@ struct NgIteration <: Method
     verbose::Bool
 end
 
-function NgIteration(; N_stages=3, max_iterations=10^3, tolerance=10^-10, verbose=true, M=2^10, dr=sqrt(π/(M+1))/(2π))
+function NgIteration(; N_stages=3, max_iterations=10^3, tolerance=10^-10, verbose=true, M=1000, dr=10.0/M)
     @assert max_iterations > 0 
     @assert tolerance > 0 
     @assert N_stages > 0
@@ -156,13 +156,15 @@ function DensityRamp(method, densities::AbstractVector{T}; verbose=true) where T
         return DensityRamp(method, densities, verbose)
     elseif T <: AbstractVector
         @assert issorted(sum.(densities))
+        Ns = length(densities[1])
+        @assert allequal(length.(densities))
         densities = Diagonal.(SVector{Ns}.(densities))
         return DensityRamp(method, densities, verbose)
-    elseif T <: AbstractMatrix
+    elseif T <: Diagonal
         @assert issorted(sum.(densities))
-        @assert all(isdiag.(densities))
         return DensityRamp(method, densities, verbose)
     end
+    error("Invalid type for densities.")
 end
 
 """

src/Solvers/DensityRamp.jl

@@ -1,3 +1,16 @@
+function recast_γ(γ::Array{T, 3}) where {T}
+    Ns = size(γ, 2)
+    Nk = size(γ, 1)
+    γ2 = zeros(SMatrix{Ns, Ns, T, Ns^2}, Nk)
+    for i = 1:Nk
+        γ2[i] = γ[i, :, :]
+    end
+    return γ2
+end
+
+recast_γ(γ::Array{T, 1}) where {T} = γ
+
+
 function solve(system::SimpleLiquid, closure::Closure, method::DensityRamp)
     densities = method.densities
     ρtarget = system.ρ
@@ -14,15 +27,17 @@ function solve(system::SimpleLiquid, closure::Closure, method::DensityRamp)
             println("\nSolving the system at ρ = $(densities[i]).\n")
         end
         system.ρ = densities[i]
-        sol = solve(system, closure, method.method, init=γ_old)
+        γ_old2 = recast_γ(γ_old)
+        sol = solve(system, closure, method.method, init=γ_old2)
         push!(sols, sol)
         @. γ_old =  sol.gr - one(eltype(sol.gr)) - sol.cr
     end
     if method.verbose
         println("\nSolving the system at ρ = $(ρtarget).\n")
     end
     system.ρ = ρtarget
-    sol = solve(system, closure, method.method, init=γ_old)
+    γ_old2 = recast_γ(γ_old)
+    sol = solve(system, closure, method.method, init=γ_old2)
     push!(sols, sol)
     return sols
 end 

src/Solvers/FourierIteration.jl

@@ -57,7 +57,7 @@ function solve(system::SimpleLiquid{dims, species, T1, T2, P}, closure::Closure,
     end
     if method.verbose 
         print("Converged after $iteration iterations, ")
-        println("the error is $(round(err, digits=ceil(Int, 1-log10(tolerance)))).")
+        println("the error is $(round(err, digits=ceil(Int, 1-log10(err)))).")
     end
     c = C ./ r
     g = find_g_from_c_and_Γ(c, Γ_new, r)

src/Solvers/NgIteration.jl

@@ -81,7 +81,7 @@ function solve(system::SimpleLiquid{dims, 1, T1, T2, P}, closure::Closure, metho
 
     if method.verbose 
         print("Converged after $iteration iterations, ")
-        println("the error is $(round(err, digits=ceil(Int, 1-log10(tolerance)))).")
+        println("the error is $(round(err, digits=ceil(Int, 1-log10(err)))).")
     end
     c = C ./ r
     g = find_g_from_c_and_Γ(c, Γ_new, r)
@@ -122,21 +122,28 @@ function solve(system::SimpleLiquid{dims, species, T1, T2, P}, closure::Closure,
     Ĉ = copy(mayer_f)
     Γ_new = copy(mayer_f)
 
+
     gn = initialize_vector_of_vectors(TT, N_stages+1, Ns*Ns*Nr) # first element is g_n, second is g_{n-1} etc
     fn = initialize_vector_of_vectors(TT, N_stages+1, Ns*Ns*Nr)
     dn = initialize_vector_of_vectors(TT, N_stages+1, Ns*Ns*Nr)
     d0n = initialize_vector_of_vectors(TT, N_stages, Ns*Ns*Nr) # first element is d01 second is d02 etc
-    if !(isnothing(init))
-        fn[end] .= init.*r
-    else
-        fn[end] .= zero(eltype(eltype(fn)))
-    end
+    # the elements of fn etc are vectors of length Ns*Ns*Nr
+
+
+
     A = zeros(TT, N_stages, N_stages)
     b = zeros(TT, N_stages)
     Γ_new_full = reshape(reinterpret(reshape, TT, Γ_new), Ns*Ns*Nr) # for going back and forth between vec{float} and vec{Smat}
     gn_red = reinterpret_vector_of_vectors(gn, T, Ns*Ns, Nr)#[reinterpret(reshape, T, reshape(gn[i], (Ns*Ns, Nr))) for i in eachindex(gn)]
     fn_red = reinterpret_vector_of_vectors(fn, T, Ns*Ns, Nr)#[reinterpret(reshape, T, reshape(fn[i], (Ns*Ns, Nr))) for i in eachindex(fn)]
 
+    if !(isnothing(init))
+        fn_red[end] .= init.*r
+    else
+        @show zero(eltype(eltype(fn_red)))
+        fn_red[end] .= (zero(eltype(eltype(fn_red))), )
+    end
+
     max_iterations = method.max_iterations
     tolerance = method.tolerance
 
@@ -189,7 +196,7 @@ function solve(system::SimpleLiquid{dims, species, T1, T2, P}, closure::Closure,
 
     if method.verbose 
         print("Converged after $iteration iterations, ")
-        println("the error is $(round(err, digits=ceil(Int, 1-log10(tolerance)))).")
+        println("the error is $(round(err, digits=ceil(Int, 1-log10(err)))).")
     end
     c = C ./ r
     g = find_g_from_c_and_Γ(c, Γ_new, r)

test/test_NgIteration_HS.jl

@@ -78,9 +78,9 @@ kBT = 1.0
 pot = HardSpheres(1.0)
 system = SimpleLiquid(dims, ρ, kBT, pot)
 closure = PercusYevick()
-method = NgIteration(tolerance=10^-10, N_stages=5, M=M, verbose=false, max_iterations=10^3)
+method = NgIteration(tolerance=10^-10, N_stages=5, M=M, verbose=false, max_iterations=10^3, dr=sqrt(π/(M+1))/(2π))
 sol = solve(system, closure, method)
-sol2 = solve(system, closure, Exact(M=M))
+sol2 = solve(system, closure, Exact(M=M, dr=sqrt(π/(M+1))/(2π)))
 
 atol = 0.1
 
@@ -92,4 +92,4 @@ atol = 0.1
 
 end 
 
-test2()
+test2()