Initial commit

Author: IlianPihlajamaa

.gitattributes

@@ -0,0 +1,2 @@
+# Auto detect text files and perform LF normalization
+* text=auto

.gitignore

@@ -0,0 +1,24 @@
+# Files generated by invoking Julia with --code-coverage
+*.jl.cov
+*.jl.*.cov
+
+# Files generated by invoking Julia with --track-allocation
+*.jl.mem
+
+# System-specific files and directories generated by the BinaryProvider and BinDeps packages
+# They contain absolute paths specific to the host computer, and so should not be committed
+deps/deps.jl
+deps/build.log
+deps/downloads/
+deps/usr/
+deps/src/
+
+# Build artifacts for creating documentation generated by the Documenter package
+docs/build/
+docs/site/
+
+# File generated by Pkg, the package manager, based on a corresponding Project.toml
+# It records a fixed state of all packages used by the project. As such, it should not be
+# committed for packages, but should be committed for applications that require a static
+# environment.
+Manifest.toml

LICENSE

@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2023 Ilian Pihlajamaa
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

Project.toml

@@ -0,0 +1,12 @@
+name = "OrnsteinZernike"
+uuid = "51399321-c693-48c0-994d-05e3df0250f2"
+authors = ["Ilian Pihlajamaa <ilian.pihlajamaa@gmail.com>"]
+version = "0.1.0"
+
+[deps]
+BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
+FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"

README.md

@@ -0,0 +1,2 @@
+# OrnsteinZernike.jl
+ A generic solver for Ornstein-Zernike equations from liquid state theory

src/Closures.jl

@@ -0,0 +1,46 @@
+abstract type Closure end
+
+struct PercusYevick <: Closure end
+struct HypernettedChain <: Closure end
+struct MeanSpherical <: Closure end
+
+
+
+# function bridge_function(::PercusYevick, γ, _, _)
+#     oneunit = one.(γ)
+#     B = @. log(Complex(oneunit + γ)) - γ
+#     return B
+# end
+
+function bridge_function(::HypernettedChain, γ, _, _)
+    zerounit = zero.(γ)
+    B = zerounit
+    return B
+end
+
+function bridge_function(::MeanSpherical, γ, u_long_range, _) 
+    oneunit = one.(γ)
+    s = @. γ - u_long_range # temp s = γ*
+    B = @. log(oneunit + s) - s
+    return B
+end
+
+function c_closure_from_γ(closure, r, mayer_f, γ, u_long_range)
+    B = bridge_function(closure, γ, u_long_range, r)
+    myone = one.(B)
+    c = @. -myone - γ + (mayer_f + myone)*exp(γ)*real(exp(B))
+    return c
+end
+
+
+function cmulr_closure_from_Γmulr(closure::Closure, r, mayer_f, Γmulr, u_long_range)
+    γ = Γmulr/r
+    return r*c_closure_from_γ(closure, r, mayer_f, γ, u_long_range) 
+end
+
+
+function cmulr_closure_from_Γmulr(::PercusYevick, r, mayer_f::T, Γmulr::T, u_long_range::T) where T
+    return  @. mayer_f*(r + Γmulr)
+end
+
+

src/FourierTransforms.jl

@@ -0,0 +1,82 @@
+function get_forward_and_backward_plan(::SimpleLiquid{3, species, T1, T2, P}, F) where {species, T1, T2, P}
+    return find_fourier_plans_3d(F)
+end
+
+function get_forward_and_backward_plan(::SimpleLiquid{2, species, T1, T2, P}, F) where {species, T1, T2, P}
+    return find_fourier_plans_2d(F)
+end
+
+function find_fourier_plans_3d(F::Vector{T}) where T
+    if T <: Number
+        forward = FFTW.plan_r2r!(copy(F), FFTW.RODFT00; flags=FFTW.ESTIMATE)
+        backward = FFTW.plan_r2r!(copy(F), FFTW.RODFT00; flags=FFTW.ESTIMATE)
+        return forward, backward
+    elseif T <: AbstractArray
+        F2 = copy(F)
+        Nspecies = size(F2[1], 1)
+        elT = eltype(T)
+        F2 = reinterpret(reshape, elT, F2)
+        if Nspecies == 1
+            forward = FFTW.plan_r2r!(F2, FFTW.RODFT00; flags=FFTW.ESTIMATE)
+            backward = FFTW.plan_r2r!(F2, FFTW.RODFT00; flags=FFTW.ESTIMATE)
+            return forward, backward
+        else
+            forward = FFTW.plan_r2r!(F2, FFTW.RODFT00, 2, flags=FFTW.ESTIMATE)
+            backward = FFTW.plan_r2r!(F2, FFTW.RODFT00, 2, flags=FFTW.ESTIMATE)
+            return forward, backward
+        end
+    end
+end
+
+function inverse_radial_fourier_transform_3d(F̂, r, k)
+    M = length(r)
+    @assert length(k) == length(F̂) ==  M
+    dk = k[2] - k[1]
+    dr = r[2] - r[1]
+    # @assert dk*dr ≈ π/(M+1)
+    F = FFTW.r2r(F̂, FFTW.RODFT00)*dk/(4π^2)
+    return F
+end
+
+function radial_fourier_transform_3d(F, r, k)
+    M = length(r)
+    @assert length(k) == length(F) ==  M
+    dk = k[2] - k[1]
+    dr = r[2] - r[1]
+    # @assert dk*dr ≈ π/(M+1)
+    F̂ = FFTW.r2r(F, FFTW.RODFT00)*2π*dr
+    return F̂ 
+end
+
+
+"""
+    fourier!(F̂, F, plan, equation::OZEquation{3, T1, T2, T3}) where {T1, T2, T3}
+
+computes the three dimensional radial fourier transform of F = r*f(r), returning F̂ = k*f̂(k), where f̂ is the fourier transform of f.
+it uses the discrete sine tranform provided by the r2r function of FFTW internally.
+"""
+function fourier!(F̂::Vector{T}, F::Vector{T}, plan::FFTW.r2rFFTWPlan, dr) where T
+    @. F̂ = F*2π*dr
+    if T <: Number
+        plan*F̂
+    elseif T<:AbstractMatrix
+        F̂2 = reinterpret(reshape, eltype(T), F̂)
+        plan*F̂2
+    end
+end
+
+"""
+inverse_fourier!(F, F̂, plan, equation::OZEquation{3, T1, T2, T3}) where {T1, T2, T3}
+
+computes the three dimensional radial fourier transform of F̂ = k*f̂(k), returning F = r*f(r), where f̂ is the fourier transform of f.
+it uses the discrete sine tranform provided by the r2r function of FFTW internally.
+"""
+function inverse_fourier!(F::AbstractVector{T}, F̂::AbstractVector{T}, plan::FFTW.r2rFFTWPlan, dk) where T
+    @. F = F̂ * dk/(4π^2)
+    if T <: Number
+        plan*F
+    elseif T<:AbstractMatrix
+        F2 = reinterpret(reshape, eltype(T), F)
+        plan*F2
+    end
+end

src/OrnsteinZernike.jl

@@ -0,0 +1,12 @@
+module OrnsteinZernike
+using FFTW, QuadGK, StaticArrays, LinearAlgebra, Plots, BlockArrays
+
+include("Systems.jl")
+include("Potentials.jl")
+include("Closures.jl")
+include("FourierTransforms.jl")
+include("Solutions.jl")
+include("Solvers.jl")
+include("Thermodynamics.jl")
+
+end # module OrnsteinZernike

src/Potentials.jl

@@ -0,0 +1,81 @@
+abstract type Potential end
+
+struct SingleComponentHardSpheres <: Potential end
+
+function evaluate_potential(::SingleComponentHardSpheres, r::Number)
+    return ifelse(r < oneunit(r), Inf64, 0.0)*oneunit(r)
+end
+
+struct MultiComponentHardSpheres{N_components, T} <: Potential 
+    D::SVector{N_components, T}
+end
+
+function MultiComponentHardSpheres(D::AbstractVector{T}) where T<:Number
+    N = length(D)
+    return MultiComponentHardSpheres{N,T}(SVector{N,T}(D))
+end
+
+function evaluate_potential(potential::MultiComponentHardSpheres{N,T}, r::Number) where {N,T}
+    pot = MMatrix{N, N, T, N*N}(undef)
+    D = potential.D
+    for species1 = 1:N
+        for species2 = 1:N
+            pot[species2, species1] = ifelse(r < (D[species1]+D[species2])/2, Inf64, 0.0)*oneunit(r)
+        end
+    end
+    return SMatrix(pot)
+end
+
+"""
+exp(- beta * u) - 1.
+"""
+function find_mayer_f_function(system::SimpleLiquid{dims, species, T1, T2, P}, r::Number, β::Number) where {dims, species, T1, T2, P}
+    U = evaluate_potential(system.potential, r)
+    return @. exp(-β*U) - 1.0
+end
+
+function find_mayer_f_function(system::SimpleLiquid{dims, species, T1, T2, P}, r::AbstractArray, β::Number) where {dims, species, T1, T2, P}
+    return find_mayer_f_function.((system, ), r, β)
+end
+
+evaluate_potential_derivative(potential::SingleComponentHardSpheres, r) = 0.0
+
+function evaluate_potential_derivative(potential::Potential, r)
+    #check for discontinuities
+    ϵ = sqrt(eps(r))
+
+
+    discs = discontinuities(potential)
+    for discontinuity in discs
+        if any(abs.(discontinuity .- r) .< ϵ)
+            error("Trying to evaluate the derivative of the potential at the discontinuity. To fix, define a specialized method for `evaluate_potential_derivative(potential::MyPotential, r)`")
+        end
+    end 
+
+    u2 = evaluate_potential(potential, r+ϵ)
+    u1 = evaluate_potential(potential, r-ϵ)
+    if isinf(u2) && isinf(u1)
+        return zero(u2)
+    end
+    return (u2-u1)/(2ϵ)
+end
+
+function discontinuities(::Potential)
+    return Float64[]
+end
+
+function discontinuities(::SingleComponentHardSpheres)
+    return [1.0]
+end
+
+function discontinuities(potential::MultiComponentHardSpheres)
+    D = potential.D
+    T = eltype(D)
+    disc = MMatrix{N, N, T, N*N}(undef)
+    for species1 = 1:N
+        for species2 = 1:N
+            disc[species2, species1] =  (D[species1]+D[species2])/2
+        end
+    end
+    return [SMatrix(disc)]
+end

src/Solutions.jl

@@ -0,0 +1,27 @@
+struct OZSolution{T1, T2}
+    r::T1
+    k::T1
+    gr::T2
+    Sk::T2
+    Ck::T2
+    Cr::T2
+end
+
+function convert_vecofmat_to_3darr(a)
+    elT = eltype(eltype(a))
+    Ns1, Ns2 = size(a[1])
+    Nr = length(a)
+    a = reinterpret(reshape, elT, a)
+    a = reshape(a, Ns1, Ns2, Nr)
+    a = permutedims(a, (3,1,2))
+    return a
+end
+
+
+function OZSolution(r::T1, k::T1, gr::T, Sk::T, Ck::T, Cr::T) where {T1,T<:Vector{<:AbstractMatrix}}
+    gr = convert_vecofmat_to_3darr(gr)
+    Sk = convert_vecofmat_to_3darr(Sk)
+    Ck = convert_vecofmat_to_3darr(Ck)
+    Cr = convert_vecofmat_to_3darr(Cr)
+    OZSolution(r, k, gr, Sk, Ck, Cr)
+end