diff --git a/dev/.documenter-siteinfo.json b/dev/.documenter-siteinfo.json index e97a28a..adf7b57 100644 --- a/dev/.documenter-siteinfo.json +++ b/dev/.documenter-siteinfo.json @@ -1 +1 @@ -{"documenter":{"julia_version":"1.9.3","generation_timestamp":"2023-10-23T07:41:52","documenter_version":"1.1.1"}} \ No newline at end of file +{"documenter":{"julia_version":"1.9.3","generation_timestamp":"2023-10-23T07:43:54","documenter_version":"1.1.1"}} \ No newline at end of file diff --git a/dev/API.html b/dev/API.html index 3e6a923..f9a835a 100644 --- a/dev/API.html +++ b/dev/API.html @@ -1,6 +1,6 @@ -API · OrnsteinZernike.jl Documentation
GitHub

The OrnsteinZernike Module

OrnsteinZernike — Module

A generic solver package for Ornstein-Zernike equations from liquid state theory

source

Module Index

Detailed API

OrnsteinZernike.OZSolution — Type
OZSolution

Holds the solution of an Ornstein Zernike problem.

Fields:

  • r: vector of distances
  • k: vector of wave numbers
  • gr: radial distribution function
  • Sk: static structure factor
  • ck: direct correlation function in k space
  • cr: direct correlation function in real space

if the system was a single-component system, gr, Sk, ck and cr are vectors. If instead the system was a multicomponent one, they are three dimensional vectors, where the first dimension contains the values along r, and the second and third dimension contain the data for the species.

source
OrnsteinZernike.SimpleLiquid — Type
SimpleLiquid{dims, ...} <: System

Holds information about a homogeneous, isotropic system with radially symmetric pair interactions. dims is the dimensionality.

Construct using

SimpleLiquid(dims, ρ, kBT, potential)

Fields:

  • ρ: number density, must be either a Number in case of a single component system, or a Vector in case of a mixture. In the latter case, each element contains the number density of the respective component.
  • kBT: thermal energy
  • potential::Potential: the interaction potential.

Examples:

ρ = 0.5; kBT = 1.1; dims = 3
+API · OrnsteinZernike.jl Documentation
GitHub
The OrnsteinZernike Module
OrnsteinZernike — Module
A generic solver package for Ornstein-Zernike equations from liquid state theory
source
Module Index
Detailed API
OrnsteinZernike.OZSolution — Type
OZSolution
Holds the solution of an Ornstein Zernike problem. 
Fields:
  • r: vector of distances
  • k: vector of wave numbers
  • gr: radial distribution function    
  • Sk: static structure factor
  • ck: direct correlation function in k space
  • cr: direct correlation function in real space
if the system was a single-component system, gr, Sk, ck and cr are vectors.  If instead the system was a multicomponent one, they are three dimensional vectors,  where the first dimension contains the values along r, and the second and third dimension contain the data for the species.
source
OrnsteinZernike.SimpleLiquid — Type
SimpleLiquid{dims, ...} <: System
Holds information about a homogeneous, isotropic system with radially symmetric pair interactions. dims is the dimensionality.
Construct using
SimpleLiquid(dims, ρ, kBT, potential)
Fields:
  • ρ: number density, must be either a Number in case of a single component system, or a Vector in case of a mixture. In the latter case, each element contains the number density of the respective component.
  • kBT: thermal energy
  • potential::Potential: the interaction potential.  
Examples:
ρ = 0.5; kBT = 1.1; dims = 3
 pot = SingleComponentHardSpheres()
 system = SimpleLiquid(dims, ρ, kBT, pot)
ρ = [0.5, 0.1]; kBT = 5.2; dims = 3
 pot = MultiComponentHardSpheres([1.0, 0.8])
-system = SimpleLiquid(dims, ρ, kBT, pot)
source
OrnsteinZernike.compute_compressibility — Method
compute_compressibility(sol::OZSolution, system::SimpleLiquid)
Computes the isothermal compressibility χ of the system
uses the formula 1/ρkBTχ = 1 - ρ ĉ(k=0) for single component systems and 1/ρkBTχ = 1 - ρ Σᵢⱼ ĉᵢⱼ(k=0) for mixtures.  Eq. (3.6.16) in Hansen and McDonald
source
OrnsteinZernike.compute_excess_energy — Method
compute_excess_energy(sol::OZSolution, system::SimpleLiquid)
Computes the excess energy per particle Eₓ, such that E = (dims/2kBT + Eₓ)N.
uses the formula Eₓ = 1/2 ρ ∫dr g(r) u(r) for single component systems and Eₓ = 1/2 ρ Σᵢⱼ xᵢxⱼ ∫dr gᵢⱼ(r) uᵢⱼ(r) for mixtures. Here x is the concentration fraction xᵢ=ρᵢ/sum(ρ).
source
OrnsteinZernike.compute_virial_pressure — Method
compute_virial_pressure(sol::OZSolution, system::SimpleLiquid)
Computes the pressure via the virial route
uses the formula p = kBTρ - 1/6 ρ^2 ∫dr r g(r) u'(r) for single component systems and p =  kBT Σᵢρᵢ - 1/6 Σᵢⱼ ρᵢρⱼ ∫dr r gᵢⱼ(r) u'ᵢⱼ(r) for mixtures.
It handles discontinuities in the interaction potential analytically if discontinuities(potential) is defined. For additional speed/accuracy define a method of evaluate_potential_derivative(potential, r::Number) that analytically computes du/dr.  By default this is done using finite differences.
source
OrnsteinZernike.solve — Function
solve(system::SimpleLiquid, closure::Closure, method::Method)
Solves the system system using the closure closure with method method.
solve(system::SimpleLiquid, closure::Closure)
Solves the system system using the closure closure with the default method NgIteration().
source

+system = SimpleLiquid(dims, ρ, kBT, pot)
source
OrnsteinZernike.compute_compressibility — Method
compute_compressibility(sol::OZSolution, system::SimpleLiquid)

Computes the isothermal compressibility χ of the system

uses the formula 1/ρkBTχ = 1 - ρ ĉ(k=0) for single component systems and 1/ρkBTχ = 1 - ρ Σᵢⱼ ĉᵢⱼ(k=0) for mixtures. Eq. (3.6.16) in Hansen and McDonald

source
OrnsteinZernike.compute_excess_energy — Method
compute_excess_energy(sol::OZSolution, system::SimpleLiquid)

Computes the excess energy per particle Eₓ, such that E = (dims/2kBT + Eₓ)N.

uses the formula Eₓ = 1/2 ρ ∫dr g(r) u(r) for single component systems and Eₓ = 1/2 ρ Σᵢⱼ xᵢxⱼ ∫dr gᵢⱼ(r) uᵢⱼ(r) for mixtures. Here x is the concentration fraction xᵢ=ρᵢ/sum(ρ).

source
OrnsteinZernike.compute_virial_pressure — Method
compute_virial_pressure(sol::OZSolution, system::SimpleLiquid)

Computes the pressure via the virial route

uses the formula p = kBTρ - 1/6 ρ^2 ∫dr r g(r) u'(r) for single component systems and p = kBT Σᵢρᵢ - 1/6 Σᵢⱼ ρᵢρⱼ ∫dr r gᵢⱼ(r) u'ᵢⱼ(r) for mixtures.

It handles discontinuities in the interaction potential analytically if discontinuities(potential) is defined. For additional speed/accuracy define a method of evaluate_potential_derivative(potential, r::Number) that analytically computes du/dr. By default this is done using finite differences.

source
OrnsteinZernike.solve — Function
solve(system::SimpleLiquid, closure::Closure, method::Method)

Solves the system system using the closure closure with method method.

solve(system::SimpleLiquid, closure::Closure)

Solves the system system using the closure closure with the default method NgIteration().

source
diff --git a/dev/Accuracy-22557a30.svg b/dev/Accuracy-ba25c58c.svg similarity index 80% rename from dev/Accuracy-22557a30.svg rename to dev/Accuracy-ba25c58c.svg index 4e4702e..b83d983 100644 --- a/dev/Accuracy-22557a30.svg +++ b/dev/Accuracy-ba25c58c.svg @@ -1,53 +1,53 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/Accuracy.html b/dev/Accuracy.html index 97960fd..da24b5a 100644 --- a/dev/Accuracy.html +++ b/dev/Accuracy.html @@ -40,4 +40,4 @@ The pressure = 0.58613053 with M = 100000 gridpoints.

We can see that the method has well-behaved second order convergence, and that with $M=10^4$, we get almost 6 digits of relative accuracy.

η = ρ/6*π
 p_exact = (1+2η+3η^2)/(1-η)^2-1.0
 scatter(M_array, abs.(p.-p_exact)./p_exact)
-plot!(ylabel="relative error", xlabel="M", xscale=:log, yscale=:log)
Example block output

In principle, these results can be extrapolated to improve the accuracy further.

+plot!(ylabel="relative error", xlabel="M", xscale=:log, yscale=:log)Example block output

In principle, these results can be extrapolated to improve the accuracy further.

diff --git a/dev/Closures.html b/dev/Closures.html index 2e9c49b..f0c648f 100644 --- a/dev/Closures.html +++ b/dev/Closures.html @@ -1,4 +1,4 @@ -Closures · OrnsteinZernike.jl Documentation
GitHub

Closures

This package defines several closure relations that can be used out of the box. It is straightforward to implement your own closure, see Defining your own closure.

Some closures, for example the Rogers-Young closure, include free parameters that may be fixed by the requirement of thermodynamic consistency. See the Thermodynamic Consistency page for an example.

Some closures use a renormalized indirect correlation function $\gamma^*(r) = \gamma(r) - u_{LR}(r)$ instead of the standard one. Here, $u_{LR}(r)$ is the long range tail of the interaction potential. There are several ways in which the interaction potential can be split into a short-range and long range part, the most common one is the Weeks-Chandler-Andersen construction. In order to use these ...

Implemented Closures

Below is an alphabetical list of implemented closures. We use the notation shown in the Theory section.

OrnsteinZernike.BallonePastoreGalliGazzillo — Type
BallonePastoreGalliGazzillo <: Closure

Implements the Ballone-Pastore-Galli-Gazzillo closure $b(r) = (1 + s \gamma(r))^{1/s} - \gamma(r) - 1 $. Here, $s$ is a free parameter that can be determined with thermodynamic consistency.

Example:

closure = BallonePastoreGalliGazzillo(1.5)

References:

Ballone, P., et al. "Additive and non-additive hard sphere mixtures: Monte Carlo simulation and integral equation results." Molecular Physics 59.2 (1986): 275-290.

source
OrnsteinZernike.BomontBretonnet — Type
BomontBretonnet <: Closure

Implements the Bomont-Bretonnet closure $b(r) = \sqrt{1+2\gamma^*(r) + f \gamma^*(r)^2} - \gamma^*(r) - 1 $. Here $\gamma^* = \gamma - u_{LR}$ , in which $ u_{LR}$ is the long range tail of the potential, and $f$ is a free parameter that can be determined with thermodynamic consistency.

Example:

closure = BomontBretonnet(0.5)

References:

Bomont, J. M., and J. L. Bretonnet. "A self-consistent integral equation: Bridge function and thermodynamic properties for the Lennard-Jones fluid." The Journal of chemical physics 119.4 (2003): 2188-2191.

source
OrnsteinZernike.CarbajalTinoko — Type

CarbajalTinoko <: Closure

Implements the Carbajal-Tinoko closure $e(r;\alpha)\left[(2-y(r))e^{y(r)}-2-y(r)\right]/\left(e^{y(r)}-1\right)$ where $e(r;\alpha)=3\exp(\alpha r)$ for $\alpha <0$ and $e(r;\alpha) = 3+\alpha$ otherwise.

Example:

closure = CarbajalTinoko(0.4)

References:

Carbajal-Tinoco, Mauricio D. "Thermodynamically consistent integral equation for soft repulsive spheres." The Journal of chemical physics 128.18 (2008).

source
OrnsteinZernike.CharpentierJackse — Type
CharpentierJackse <: Closure

Implements the Charpentier-Jackse closure $b(r) = \frac{1}{2\alpha}\left(\sqrt{1+4\alpha\gamma^*(r) } - 2\alpha\gamma^*(r) - 1\right) $. Here $\gamma^* = \gamma - u_{LR}$ , in which $ u_{LR}$ is the long range tail of the potential, and $\alpha$ is a free parameter that can be determined with thermodynamic consistency.

Example:

closure = CharpentierJackse(0.5)

References:

Charpentier, I., and N. Jakse. "Exact numerical derivatives of the pair-correlation function of simple liquids using the tangent linear method." The Journal of Chemical Physics 114.5 (2001): 2284-2292.

source
OrnsteinZernike.ChoudhuryGhosh — Type
ChoudhuryGhosh <: Closure

Implements the Choudhury-Ghosh closure $b(r) = -\frac{-\gamma^*(r)^2}{2(1+\alpha \gamma^*(r))} $ for $\gamma^*(r) >0$, and $b(r)=-\gamma^*(r)^2/2$ otherwise. Here $\gamma^* = \gamma - u_{LR}$ , in which $ u_{LR}$ is the long range tail of the potential, and $\alpha$ is a free parameter that is determined by an empirical relation.

Example:

α(ρ) = 1.01752 - 0.275ρ # see the reference for this empirical relation
-closure = ChoudhuryGhosh(α(0.4))

References:

Choudhury, Niharendu, and Swapan K. Ghosh. "Integral equation theory of Lennard-Jones fluids: A modified Verlet bridge function approach." The Journal of chemical physics 116.19 (2002): 8517-8522.

source
OrnsteinZernike.DuhHaymet — Type
DuhHaymet <: Closure

Implements the Duh-Haymet closure $b(r) = \frac{-\gamma^*(r)^2}{2\left[1+\left(\frac{5\gamma^*(r)+11}{7\gamma^*(r)+9}\right)\gamma^*(r)\right]} $ for $\gamma^*(r) >0$, and $b(r)=-\gamma^*(r)^2/2$ otherwise. Here $\gamma^* = \gamma - u_{LR}$ , in which $ u_{LR}$ is the long range tail of the potential,

Example:

closure = DuhHaymet()

References:

source
OrnsteinZernike.ExtendedRogersYoung — Type
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:

closure = ExtendedRogersYoung(0.5, 0.5) # order is α, a

References: J. Chem. Phys. 128, 184507 (2008)

source
OrnsteinZernike.HypernettedChain — Type
HypernettedChain

Implements the Hypernetted Chain closure $c(r) = (f(r)+1)\exp(\gamma(r)) - \gamma(r) - 1$, or equivalently $b(r) = 0$.

Example:

closure = HypernettedChain()
source
OrnsteinZernike.Khanpour — Type
Khanpour <: Closure

Implements the Khanpour closure $b(r) = \frac{1}{\alpha}\ln(1+\alpha\gamma(r)) - \gamma $. Here $\alpha$ is a free parameter that can be determined with thermodynamic consistency.

Example:

closure = Khanpour(0.5)

References:

Khanpour, Mehrdad. "A unified derivation of Percus–Yevick and hyper-netted chain integral equations in liquid state theory." Molecular Physics 120.5 (2022): e2001065.

source
OrnsteinZernike.Lee — Type
Lee <: Closure

Implements the Lee closure $b(r) = -\frac{\zeta\gamma^*(r)^2}{2} \left( 1- \frac{\phi \alpha \gamma^*(r)}{1 + \alpha\gamma^*(r)} \right) $. Here $\gamma^* = \gamma + ρ f(r)/2$ , in which $ f(r)$ is the Mayer-f function and $\rho$ the density. Additionally, $\zeta$,$\phi$, and, $\alpha$ are free parameters that can be determined with thermodynamic consistency or zero-separation theorems.

Example:

closure = Lee(1.073, 1.816, 1.0, 0.4) # ζ, ϕ, α, ρ

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.

source
OrnsteinZernike.MartynovSarkisov — Type
MartynovSarkisov <: Closure

Implements the Martynov-Sarkisov Closure b(r) = -\sqrt{1+2\gamma}-1-\gamma, which is constructed for the 3d hard-sphere liquid.

Example:

closure = MartynovSarkisov()

References:

Martynov, G. A., and G. N. Sarkisov. "Exact equations and the theory of liquids. V." Molecular Physics 49.6 (1983): 1495-1504.

source
OrnsteinZernike.MeanSpherical — Type
MeanSpherical

Implements the MSA closure $c(r) = -\beta u(r)$, or equivalently $b(r) = \ln(\gamma(r) - \beta u(r) + 1) - γ(r) + \beta u(r)$.

Example:

closure = MeanSpherical()
source
OrnsteinZernike.ModifiedHypernettedChain — Type

ModifiedHypernettedChain <: Closure

Implements the Modified 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<a_4$ and $b_{HS}(r)=\left(A_1 \exp(-a_5(x-a_4))\sin(A_2(x-a_4))/r\right)^2$ is the hard sphere bridge function found in Malijevský & Labík. The parameters are defined as

\[x = r-1\]

\[A_1 = (a_1+a_2 a_4)(a_4-a_3)(a_4+1)/(A_2 a_3 a_4)\]

\[A_2 = \pi / (a_6 - a_4 - 1)\]

\[a_1 = \eta (1.55707 - 1.85633\eta) / (1-\eta)^2\]

\[a_2 = \eta (1.28127 - 1.82134\eta) / (1-\eta)\]

\[a_3 = (0.74480 - 0.93453\eta)\]

\[a_4 = (1.17102 - 0.68230\eta)\]

\[a_5 = 0.15975/\eta^2\]

\[a_6 = (2.69757 - 0.86987\eta)\]

and $\eta$ is the volume fraction of the hard sphere reference system. This closure only works for single component systems in three dimensions.

Example:

closure = ModifiedHypernettedChain(0.4)

References:

Lado, F. "Perturbation correction for the free energy and structure of simple fluids." Physical Review A 8.5 (1973): 2548.

Malijevský, Anatol, and Stanislav Labík. "The bridge function for hard spheres." Molecular Physics 60.3 (1987): 663-669.

source
OrnsteinZernike.ModifiedVerlet — Type
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:

closure = ModifiedVerlet(0.2)

References:

Verlet, Loup. "Integral equations for classical fluids: I. The hard sphere case." Molecular Physics 41.1 (1980): 183-190.

source
OrnsteinZernike.PercusYevick — Type
PercusYevick

Implements the Percus-Yevick closure $c(r) = f(r)(1+\gamma(r))$, or equivalently $b(r) = \ln(1 + \gamma(r)) - γ(r)$.

Example:

closure = PercusYevick()
source
OrnsteinZernike.RogersYoung — Type
RogersYoung <: Closure

Implements the Rogers-Young closure $b(r) = \ln(\frac{\exp(f(r)\gamma(r)) - 1}{f(r)} + 1) - γ(r) $. Here $f(r)=1-\exp(-\alpha r)$, in which $\alpha$ is a free parameter, that may be chosen such that thermodynamic consistency is achieved. Example:

closure = RogersYoung(0.5)

References:

source
OrnsteinZernike.SoftCoreMeanSpherical — Type
SoftCoreMeanSpherical <: Closure

Implements the soft core mean spherical closure $b(r) = \ln(\gamma^*(r) + 1) - \gamma^*(r)$. Here $\gamma^* = \gamma - u_{LR}$ , in which $u_{LR}$ is the long range tail of the potential.

Example:

closure = SoftCoreMeanSpherical()

References:

source
OrnsteinZernike.Verlet — Type
Verlet <: Closure

Implements the Verlet Closure $b(r) = -\frac{A\gamma^2(r)/2}{1+B \gamma(r)/2}$, where by default $A=1$ and $B=8/5$. These values are tuned by the virial coefficients of the 3d hard sphere liquid.

Example:

closure = Verlet()
-closure = Verlet(A=3.0, B=4.0)

References:

Verlet, Loup. "Integral equations for classical fluids: I. The hard sphere case." Molecular Physics 41.1 (1980): 183-190.

source
OrnsteinZernike.VompeMartynov — Type
VompeMartynov <: Closure

Implements the Vompe-Martynov closure $b(r) = \sqrt{1+2\gamma^*(r) } - \gamma^*(r) - 1 $. Here $\gamma^* = \gamma - u_{LR}$ , in which $ u_{LR}$ is the long range tail of the potential, and $\alpha$ is a free parameter that can be determined with thermodynamic consistency.

Example:

closure = VompeMartynov()

References:

Vompe, A. G., and G. A. Martynov. "The bridge function expansion and the self‐consistency problem of the Ornstein–Zernike equation solution." The Journal of chemical physics 100.7 (1994): 5249-5258.

source
OrnsteinZernike.ZerahHansen — Type

ZerahHansen <: Closure

Implements the Zerah-Hansen (HMSA) (HNC-SMSA) closure $b(r) = \ln(\frac{\exp(f(r)\gamma^*(r)) - 1}{f(r)} + 1) - γ^*(r) $. Here $\gamma^* = \gamma - u_{LR}$ , in which $ u_{LR}$ is the long range tail of the potential, 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:

closure = ZerahHansen(0.5)

References:

Zerah, Gilles, and Jean‐Pierre Hansen. "Self‐consistent integral equations for fluid pair distribution functions: Another attempt." The Journal of chemical physics 84.4 (1986): 2336-2343.

source
+Closures · OrnsteinZernike.jl Documentation
GitHub

Closures

This package defines several closure relations that can be used out of the box. It is straightforward to implement your own closure, see Defining your own closure.

Some closures, for example the Rogers-Young closure, include free parameters that may be fixed by the requirement of thermodynamic consistency. See the Thermodynamic Consistency page for an example.

Some closures use a renormalized indirect correlation function $\gamma^*(r) = \gamma(r) - u_{LR}(r)$ instead of the standard one. Here, $u_{LR}(r)$ is the long range tail of the interaction potential. There are several ways in which the interaction potential can be split into a short-range and long range part, the most common one is the Weeks-Chandler-Andersen construction. In order to use these ...

Implemented Closures

Below is an alphabetical list of implemented closures. We use the notation shown in the Theory section.

OrnsteinZernike.BallonePastoreGalliGazzillo — Type
BallonePastoreGalliGazzillo <: Closure

Implements the Ballone-Pastore-Galli-Gazzillo closure $b(r) = (1 + s \gamma(r))^{1/s} - \gamma(r) - 1 $. Here, $s$ is a free parameter that can be determined with thermodynamic consistency.

Example:

closure = BallonePastoreGalliGazzillo(1.5)

References:

Ballone, P., et al. "Additive and non-additive hard sphere mixtures: Monte Carlo simulation and integral equation results." Molecular Physics 59.2 (1986): 275-290.

source
OrnsteinZernike.BomontBretonnet — Type
BomontBretonnet <: Closure

Implements the Bomont-Bretonnet closure $b(r) = \sqrt{1+2\gamma^*(r) + f \gamma^*(r)^2} - \gamma^*(r) - 1 $. Here $\gamma^* = \gamma - u_{LR}$ , in which $ u_{LR}$ is the long range tail of the potential, and $f$ is a free parameter that can be determined with thermodynamic consistency.

Example:

closure = BomontBretonnet(0.5)

References:

Bomont, J. M., and J. L. Bretonnet. "A self-consistent integral equation: Bridge function and thermodynamic properties for the Lennard-Jones fluid." The Journal of chemical physics 119.4 (2003): 2188-2191.

source
OrnsteinZernike.CarbajalTinoko — Type

CarbajalTinoko <: Closure

Implements the Carbajal-Tinoko closure $e(r;\alpha)\left[(2-y(r))e^{y(r)}-2-y(r)\right]/\left(e^{y(r)}-1\right)$ where $e(r;\alpha)=3\exp(\alpha r)$ for $\alpha <0$ and $e(r;\alpha) = 3+\alpha$ otherwise.

Example:

closure = CarbajalTinoko(0.4)

References:

Carbajal-Tinoco, Mauricio D. "Thermodynamically consistent integral equation for soft repulsive spheres." The Journal of chemical physics 128.18 (2008).

source
OrnsteinZernike.CharpentierJackse — Type
CharpentierJackse <: Closure

Implements the Charpentier-Jackse closure $b(r) = \frac{1}{2\alpha}\left(\sqrt{1+4\alpha\gamma^*(r) } - 2\alpha\gamma^*(r) - 1\right) $. Here $\gamma^* = \gamma - u_{LR}$ , in which $ u_{LR}$ is the long range tail of the potential, and $\alpha$ is a free parameter that can be determined with thermodynamic consistency.

Example:

closure = CharpentierJackse(0.5)

References:

Charpentier, I., and N. Jakse. "Exact numerical derivatives of the pair-correlation function of simple liquids using the tangent linear method." The Journal of Chemical Physics 114.5 (2001): 2284-2292.

source
OrnsteinZernike.ChoudhuryGhosh — Type
ChoudhuryGhosh <: Closure

Implements the Choudhury-Ghosh closure $b(r) = -\frac{-\gamma^*(r)^2}{2(1+\alpha \gamma^*(r))} $ for $\gamma^*(r) >0$, and $b(r)=-\gamma^*(r)^2/2$ otherwise. Here $\gamma^* = \gamma - u_{LR}$ , in which $ u_{LR}$ is the long range tail of the potential, and $\alpha$ is a free parameter that is determined by an empirical relation.

Example:

α(ρ) = 1.01752 - 0.275ρ # see the reference for this empirical relation
+closure = ChoudhuryGhosh(α(0.4))

References:

Choudhury, Niharendu, and Swapan K. Ghosh. "Integral equation theory of Lennard-Jones fluids: A modified Verlet bridge function approach." The Journal of chemical physics 116.19 (2002): 8517-8522.

source
OrnsteinZernike.DuhHaymet — Type
DuhHaymet <: Closure

Implements the Duh-Haymet closure $b(r) = \frac{-\gamma^*(r)^2}{2\left[1+\left(\frac{5\gamma^*(r)+11}{7\gamma^*(r)+9}\right)\gamma^*(r)\right]} $ for $\gamma^*(r) >0$, and $b(r)=-\gamma^*(r)^2/2$ otherwise. Here $\gamma^* = \gamma - u_{LR}$ , in which $ u_{LR}$ is the long range tail of the potential,

Example:

closure = DuhHaymet()

References:

source
OrnsteinZernike.ExtendedRogersYoung — Type
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:

closure = ExtendedRogersYoung(0.5, 0.5) # order is α, a

References: J. Chem. Phys. 128, 184507 (2008)

source
OrnsteinZernike.HypernettedChain — Type
HypernettedChain

Implements the Hypernetted Chain closure $c(r) = (f(r)+1)\exp(\gamma(r)) - \gamma(r) - 1$, or equivalently $b(r) = 0$.

Example:

closure = HypernettedChain()
source
OrnsteinZernike.Khanpour — Type
Khanpour <: Closure

Implements the Khanpour closure $b(r) = \frac{1}{\alpha}\ln(1+\alpha\gamma(r)) - \gamma $. Here $\alpha$ is a free parameter that can be determined with thermodynamic consistency.

Example:

closure = Khanpour(0.5)

References:

Khanpour, Mehrdad. "A unified derivation of Percus–Yevick and hyper-netted chain integral equations in liquid state theory." Molecular Physics 120.5 (2022): e2001065.

source
OrnsteinZernike.Lee — Type
Lee <: Closure

Implements the Lee closure $b(r) = -\frac{\zeta\gamma^*(r)^2}{2} \left( 1- \frac{\phi \alpha \gamma^*(r)}{1 + \alpha\gamma^*(r)} \right) $. Here $\gamma^* = \gamma + ρ f(r)/2$ , in which $ f(r)$ is the Mayer-f function and $\rho$ the density. Additionally, $\zeta$,$\phi$, and, $\alpha$ are free parameters that can be determined with thermodynamic consistency or zero-separation theorems.

Example:

closure = Lee(1.073, 1.816, 1.0, 0.4) # ζ, ϕ, α, ρ

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.

source
OrnsteinZernike.MartynovSarkisov — Type
MartynovSarkisov <: Closure

Implements the Martynov-Sarkisov Closure b(r) = -\sqrt{1+2\gamma}-1-\gamma, which is constructed for the 3d hard-sphere liquid.

Example:

closure = MartynovSarkisov()

References:

Martynov, G. A., and G. N. Sarkisov. "Exact equations and the theory of liquids. V." Molecular Physics 49.6 (1983): 1495-1504.

source
OrnsteinZernike.MeanSpherical — Type
MeanSpherical

Implements the MSA closure $c(r) = -\beta u(r)$, or equivalently $b(r) = \ln(\gamma(r) - \beta u(r) + 1) - γ(r) + \beta u(r)$.

Example:

closure = MeanSpherical()
source
OrnsteinZernike.ModifiedHypernettedChain — Type

ModifiedHypernettedChain <: Closure

Implements the Modified 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<a_4$ and $b_{HS}(r)=\left(A_1 \exp(-a_5(x-a_4))\sin(A_2(x-a_4))/r\right)^2$ is the hard sphere bridge function found in Malijevský & Labík. The parameters are defined as

\[x = r-1\]

\[A_1 = (a_1+a_2 a_4)(a_4-a_3)(a_4+1)/(A_2 a_3 a_4)\]

\[A_2 = \pi / (a_6 - a_4 - 1)\]

\[a_1 = \eta (1.55707 - 1.85633\eta) / (1-\eta)^2\]

\[a_2 = \eta (1.28127 - 1.82134\eta) / (1-\eta)\]

\[a_3 = (0.74480 - 0.93453\eta)\]

\[a_4 = (1.17102 - 0.68230\eta)\]

\[a_5 = 0.15975/\eta^2\]

\[a_6 = (2.69757 - 0.86987\eta)\]

and $\eta$ is the volume fraction of the hard sphere reference system. This closure only works for single component systems in three dimensions.

Example:

closure = ModifiedHypernettedChain(0.4)

References:

Lado, F. "Perturbation correction for the free energy and structure of simple fluids." Physical Review A 8.5 (1973): 2548.

Malijevský, Anatol, and Stanislav Labík. "The bridge function for hard spheres." Molecular Physics 60.3 (1987): 663-669.

source
OrnsteinZernike.ModifiedVerlet — Type
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:

closure = ModifiedVerlet(0.2)

References:

Verlet, Loup. "Integral equations for classical fluids: I. The hard sphere case." Molecular Physics 41.1 (1980): 183-190.

source
OrnsteinZernike.PercusYevick — Type
PercusYevick

Implements the Percus-Yevick closure $c(r) = f(r)(1+\gamma(r))$, or equivalently $b(r) = \ln(1 + \gamma(r)) - γ(r)$.

Example:

closure = PercusYevick()
source
OrnsteinZernike.RogersYoung — Type
RogersYoung <: Closure

Implements the Rogers-Young closure $b(r) = \ln(\frac{\exp(f(r)\gamma(r)) - 1}{f(r)} + 1) - γ(r) $. Here $f(r)=1-\exp(-\alpha r)$, in which $\alpha$ is a free parameter, that may be chosen such that thermodynamic consistency is achieved. Example:

closure = RogersYoung(0.5)

References:

source
OrnsteinZernike.SoftCoreMeanSpherical — Type
SoftCoreMeanSpherical <: Closure

Implements the soft core mean spherical closure $b(r) = \ln(\gamma^*(r) + 1) - \gamma^*(r)$. Here $\gamma^* = \gamma - u_{LR}$ , in which $u_{LR}$ is the long range tail of the potential.

Example:

closure = SoftCoreMeanSpherical()

References:

source
OrnsteinZernike.Verlet — Type
Verlet <: Closure

Implements the Verlet Closure $b(r) = -\frac{A\gamma^2(r)/2}{1+B \gamma(r)/2}$, where by default $A=1$ and $B=8/5$. These values are tuned by the virial coefficients of the 3d hard sphere liquid.

Example:

closure = Verlet()
+closure = Verlet(A=3.0, B=4.0)

References:

Verlet, Loup. "Integral equations for classical fluids: I. The hard sphere case." Molecular Physics 41.1 (1980): 183-190.

source
OrnsteinZernike.VompeMartynov — Type
VompeMartynov <: Closure

Implements the Vompe-Martynov closure $b(r) = \sqrt{1+2\gamma^*(r) } - \gamma^*(r) - 1 $. Here $\gamma^* = \gamma - u_{LR}$ , in which $ u_{LR}$ is the long range tail of the potential, and $\alpha$ is a free parameter that can be determined with thermodynamic consistency.

Example:

closure = VompeMartynov()

References:

Vompe, A. G., and G. A. Martynov. "The bridge function expansion and the self‐consistency problem of the Ornstein–Zernike equation solution." The Journal of chemical physics 100.7 (1994): 5249-5258.

source
OrnsteinZernike.ZerahHansen — Type

ZerahHansen <: Closure

Implements the Zerah-Hansen (HMSA) (HNC-SMSA) closure $b(r) = \ln(\frac{\exp(f(r)\gamma^*(r)) - 1}{f(r)} + 1) - γ^*(r) $. Here $\gamma^* = \gamma - u_{LR}$ , in which $ u_{LR}$ is the long range tail of the potential, 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:

closure = ZerahHansen(0.5)

References:

Zerah, Gilles, and Jean‐Pierre Hansen. "Self‐consistent integral equations for fluid pair distribution functions: Another attempt." The Journal of chemical physics 84.4 (1986): 2336-2343.

source
diff --git a/dev/ExtendingClosures-4819516b.svg b/dev/ExtendingClosures-40c93170.svg similarity index 89% rename from dev/ExtendingClosures-4819516b.svg rename to dev/ExtendingClosures-40c93170.svg index c1a91f0..65f6821 100644 --- a/dev/ExtendingClosures-4819516b.svg +++ b/dev/ExtendingClosures-40c93170.svg @@ -1,46 +1,46 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/ExtendingClosures.html b/dev/ExtendingClosures.html index df1441b..058c9b1 100644 --- a/dev/ExtendingClosures.html +++ b/dev/ExtendingClosures.html @@ -19,7 +19,7 @@ closure = MyHNC() sol = solve(system, closure) using Plots -plot(sol.r, sol.gr, xlims=(0,5), xlabel="r", ylabel="g(r)")Example block output

which can be compared to that of First steps.

Mixtures

In the case of multicomponent systems, instead of a number the function that is overloaded should return a StaticMatrix containing either values for $c_{ij}$ or $b_{ij}$. Since, in that case also the inputs are StaticMatrix, we can make use of Julia's broadcasting syntax to perform an closure elementwise. 

import OrnsteinZernike.closure_c_from_gamma
+plot(sol.r, sol.gr, xlims=(0,5), xlabel="r", ylabel="g(r)")
Example block output

which can be compared to that of First steps.

Mixtures

In the case of multicomponent systems, instead of a number the function that is overloaded should return a StaticMatrix containing either values for $c_{ij}$ or $b_{ij}$. Since, in that case also the inputs are StaticMatrix, we can make use of Julia's broadcasting syntax to perform an closure elementwise. 

import OrnsteinZernike.closure_c_from_gamma
 function OrnsteinZernike.closure_c_from_gamma(::MyHNC, _, mayer_f, γ, _)
     return @. (mayer_f + 1) * exp(γ) - γ - 1 # note the @.
-end
+end diff --git a/dev/ExtendingPotentials-d5f3d268.svg b/dev/ExtendingPotentials-1d0940e2.svg similarity index 93% rename from dev/ExtendingPotentials-d5f3d268.svg rename to dev/ExtendingPotentials-1d0940e2.svg index 23ea281..4eedeb0 100644 --- a/dev/ExtendingPotentials-d5f3d268.svg +++ b/dev/ExtendingPotentials-1d0940e2.svg @@ -1,52 +1,52 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/ExtendingPotentials-1490660e.svg b/dev/ExtendingPotentials-ef6e94fb.svg similarity index 88% rename from dev/ExtendingPotentials-1490660e.svg rename to dev/ExtendingPotentials-ef6e94fb.svg index 40ec07d..2c97a71 100644 --- a/dev/ExtendingPotentials-1490660e.svg +++ b/dev/ExtendingPotentials-ef6e94fb.svg @@ -1,52 +1,52 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/ExtendingPotentials.html b/dev/ExtendingPotentials.html index a3800e1..c50152a 100644 --- a/dev/ExtendingPotentials.html +++ b/dev/ExtendingPotentials.html @@ -10,7 +10,7 @@ closure = HypernettedChain() sol = solve(system, closure) using Plots -plot(sol.r, sol.gr, xlims=(0,5), xlabel="r", ylabel="g(r)")Example block output

Mixtures

In the case of multicomponent systems, instead of a number the function should return a StaticMatrix from the StaticArrays package containing values for $u_{ij}$.

Example

Suppose we want to implement the same potential for the multicomponent case: $u_{ij}(r) = \epsilon_{ij} (\frac{\sigma_{ij}}{r_{ij}})^6.$

While one could implement this in one line with broadcasting, here the function is written out fully for clarity:

using OrnsteinZernike, StaticArrays
+plot(sol.r, sol.gr, xlims=(0,5), xlabel="r", ylabel="g(r)")
Example block output

Mixtures

In the case of multicomponent systems, instead of a number the function should return a StaticMatrix from the StaticArrays package containing values for $u_{ij}$.

Example

Suppose we want to implement the same potential for the multicomponent case: $u_{ij}(r) = \epsilon_{ij} (\frac{\sigma_{ij}}{r_{ij}})^6.$

While one could implement this in one line with broadcasting, here the function is written out fully for clarity:

using OrnsteinZernike, StaticArrays
 
 function mypotential(r, p)
     # we can construct a mutable sized matrix first
@@ -36,4 +36,4 @@
 using Plots
 plot(sol.r, sol.gr[:, 1, 1], xlims=(0,5), xlabel="r", ylabel="g(r)", label="g11(r)")
 plot!(sol.r, sol.gr[:, 1, 2], xlims=(0,5), xlabel="r", ylabel="g(r)", label="g12(r)")
-plot!(sol.r, sol.gr[:, 2, 2], xlims=(0,5), xlabel="r", ylabel="g(r)", label="g22(r)")
Example block output +plot!(sol.r, sol.gr[:, 2, 2], xlims=(0,5), xlabel="r", ylabel="g(r)", label="g22(r)")Example block output diff --git a/dev/FromPython.html b/dev/FromPython.html index b4588ae..0ce0af9 100644 --- a/dev/FromPython.html +++ b/dev/FromPython.html @@ -13,4 +13,4 @@ plt.plot(sol.r, sol.gr) plt.xlim(0, 5) -plt.show()

image

See the documentation of juliacall for more information on how to call Julia from Python.

+plt.show()

image

See the documentation of juliacall for more information on how to call Julia from Python.

diff --git a/dev/GeneralWorkflow.html b/dev/GeneralWorkflow.html index b5ee340..8ac1ce6 100644 --- a/dev/GeneralWorkflow.html +++ b/dev/GeneralWorkflow.html @@ -1,2 +1,2 @@ -General Workflow · OrnsteinZernike.jl Documentation
GitHub

General Workflow

In general, solving an Ornstein-Zernike equation requires the definition of following items. The linked pages describe these in more details.

  1. The interaction potential: Interaction Potentials
  2. The system parameters: density, temperature and dimensionality: Systems
  3. The closure relation: Closures
  4. The solver: Solvers

Once these have been defined, solve can be called to solve the system at hand.

+General Workflow · OrnsteinZernike.jl Documentation
GitHub

General Workflow

In general, solving an Ornstein-Zernike equation requires the definition of following items. The linked pages describe these in more details.

  1. The interaction potential: Interaction Potentials
  2. The system parameters: density, temperature and dimensionality: Systems
  3. The closure relation: Closures
  4. The solver: Solvers

Once these have been defined, solve can be called to solve the system at hand.

diff --git a/dev/HardSphereMixture-ae3948c1.svg b/dev/HardSphereMixture-7daf2594.svg similarity index 94% rename from dev/HardSphereMixture-ae3948c1.svg rename to dev/HardSphereMixture-7daf2594.svg index 6ed2dee..7891a8d 100644 --- a/dev/HardSphereMixture-ae3948c1.svg +++ b/dev/HardSphereMixture-7daf2594.svg @@ -1,58 +1,58 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/HardSphereMixture.html b/dev/HardSphereMixture.html index 90c08d7..8f6828d 100644 --- a/dev/HardSphereMixture.html +++ b/dev/HardSphereMixture.html @@ -30,4 +30,4 @@ 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) -pExample block output +pExample block output diff --git a/dev/HighDensities-6afdfb14.svg b/dev/HighDensities-39eb1792.svg similarity index 98% rename from dev/HighDensities-6afdfb14.svg rename to dev/HighDensities-39eb1792.svg index c247aab..5f4462a 100644 --- a/dev/HighDensities-6afdfb14.svg +++ b/dev/HighDensities-39eb1792.svg @@ -1,52 +1,52 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/HighDensities.html b/dev/HighDensities.html index 5076c04..cf7028e 100644 --- a/dev/HighDensities.html +++ b/dev/HighDensities.html @@ -35,7 +35,7 @@ Sk = [0.00031758964535433254 0.00031762625283637824 … 1.0091523150484354 1.0091605655037994]
using Plots
 plot(sol.r, sol.gr, xlims=(0,5), xlabel="r", ylabel="g(r)", lw=4, label="iterative")
 sol2 = solve(system, closure, Exact(M=M; dr=dr));
-plot!(sol2.r, sol2.gr, lw=2, color=:black, label="exact")
Example block output

Which converges even though the density is clearly nonphysically high.

Sometimes convergence is accelerated if the solution of the same equations at a slightly lower density is used as an initial condition. This is especially useful if the solution is needed at many different densities. While this can be done by hand using the init keyword argument of the solve function, a convenient method DensityRamp is implemented to do this automatically. For example, consider a hard sphere system at ρ = 1.2.

using OrnsteinZernike
+plot!(sol2.r, sol2.gr, lw=2, color=:black, label="exact")
Example block output

Which converges even though the density is clearly nonphysically high.

Sometimes convergence is accelerated if the solution of the same equations at a slightly lower density is used as an initial condition. This is especially useful if the solution is needed at many different densities. While this can be done by hand using the init keyword argument of the solve function, a convenient method DensityRamp is implemented to do this automatically. For example, consider a hard sphere system at ρ = 1.2.

using OrnsteinZernike
 ρ = 1.2
 kBT = 1.0
 dims = 3
@@ -105,4 +105,4 @@
 After iteration 60, the error is 4.2188e-7.
 After iteration 70, the error is 3.1e-9.
 Converged after 80 iterations, the error is 9.0e-11.
-  0.071976 seconds (17.20 k allocations: 4.204 MiB)

While, in total, this has increased the number of iterations done, we have obtained the solution at two different densities as well.

+ 0.071976 seconds (17.20 k allocations: 4.204 MiB)

While, in total, this has increased the number of iterations done, we have obtained the solution at two different densities as well.

diff --git a/dev/OtherDimensions.html b/dev/OtherDimensions.html index 48770a4..bb0a3f0 100644 --- a/dev/OtherDimensions.html +++ b/dev/OtherDimensions.html @@ -14,4 +14,4 @@ sol = solve(system, closure, method) plot!(p, sol.r, sol.gr .+ i, label="d = $dimension") end -p

image

Currently, however, noninteger dimensions below 2.0 are not supported.

+p

image

Currently, however, noninteger dimensions below 2.0 are not supported.

diff --git a/dev/Potentials-38629949.svg b/dev/Potentials-ccf07e85.svg similarity index 85% rename from dev/Potentials-38629949.svg rename to dev/Potentials-ccf07e85.svg index 1fc4188..4c1b136 100644 --- a/dev/Potentials-38629949.svg +++ b/dev/Potentials-ccf07e85.svg @@ -1,41 +1,41 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/Potentials.html b/dev/Potentials.html index 0f76e02..77e54eb 100644 --- a/dev/Potentials.html +++ b/dev/Potentials.html @@ -3,6 +3,6 @@ r = 0.9:0.01:4.0 potential = LennardJones(1.0, 1.0) u = OrnsteinZernike.evaluate_potential(potential, r) -plot(r, u, xlabel="r", ylabel="u(r)", ylims=(-1,1), label=nothing)Example block output

Implemented interaction potentials

Below is a list of implemented closures. We use the notation shown in the Theory section.

OrnsteinZernike.CustomPotential — Type
CustomPotential

Implements a potential that evaluates a user defined function.

Expects values f, and p, which respecively are a callable and a list of parameters. The function should be called f(r::Number, p) and it should produce either a Number, in the case of a single-component system, or an SMatrix, in the case of a multicomponent system.

Example:

f = (r, p) -> 4*p[1]*((p[2]/r)^12 -  (p[2]/r)^6)
-potential = CustomPotential(f, (1.0, 1.0))
source
OrnsteinZernike.HardSpheres — Type
HardSpheres

Implements the hard-sphere pair interaction for single component systems $ u(r) = \infty$ for $r < 1$ and $u(r) = 0$ for $r > 1$, or $u_{ij}(r) = \infty$ for $r < D_{ij}$ and $u_{ij}(r) = 0$ for $r > D_{ij}$ for mixtures.

For mixtures expects either a vector $D_i$ of diameters for each of the species in which case an additive mixing rule is used $\left(D_{ij} = (D_{i}+D_{j})/2\right)$ or a matrix $D_ij$ of pair diameters.

Example:

potential = HardSpheres(1.0)

Example:

potential = HardSpheres([0.8, 0.9, 1.0])
Dij = rand(3,3)
-potential = HardSpheres(Dij)
source
OrnsteinZernike.LennardJones — Type
LennardJones

Implements the Lennard-Jones pair interaction $u(r) = 4\epsilon [(\sigma/r)^{12} - (\sigma/r)^6]$.

Expects values ϵ and σ, which respecively are the strength of the potential and particle size.

Example:

potential = LennardJones(1.0, 2.0)
source
OrnsteinZernike.PowerLaw — Type
PowerLaw

Implements the power law pair interaction $u(r) = \epsilon (\sigma/r)^{n}$.

Expects values ϵ, σ, and n, which respecively are the strength of the potential and particle size.

Example:

potential = PowerLaw(1.0, 2.0, 8)
source
OrnsteinZernike.WCADivision — Type
WCADivision{P<:Potential, T, UC}

fields

  • potential::Potential
  • cutoff : ($r_{c}$)
  • U_c : $U(r=r_c)$

Splits the potential at the cutoff point using the WCA splitting rule. This means

\[u(r) = u_{SR}(r) + U_{LR}(r),\]

where $U_{LR}(r) = u(r)$ for $r > r_{c}$, and $U(r_{c})$ for $r < r_{c}$, and $USR(r) = 0$ for $r > r_{c}$, and $U(r) - U(r_{c})$ for $r < r_{c}$.

source
+plot(r, u, xlabel="r", ylabel="u(r)", ylims=(-1,1), label=nothing)Example block output

Implemented interaction potentials

Below is a list of implemented closures. We use the notation shown in the Theory section.

OrnsteinZernike.CustomPotential — Type
CustomPotential

Implements a potential that evaluates a user defined function.

Expects values f, and p, which respecively are a callable and a list of parameters. The function should be called f(r::Number, p) and it should produce either a Number, in the case of a single-component system, or an SMatrix, in the case of a multicomponent system.

Example:

f = (r, p) -> 4*p[1]*((p[2]/r)^12 -  (p[2]/r)^6)
+potential = CustomPotential(f, (1.0, 1.0))
source
OrnsteinZernike.HardSpheres — Type
HardSpheres

Implements the hard-sphere pair interaction for single component systems $ u(r) = \infty$ for $r < 1$ and $u(r) = 0$ for $r > 1$, or $u_{ij}(r) = \infty$ for $r < D_{ij}$ and $u_{ij}(r) = 0$ for $r > D_{ij}$ for mixtures.

For mixtures expects either a vector $D_i$ of diameters for each of the species in which case an additive mixing rule is used $\left(D_{ij} = (D_{i}+D_{j})/2\right)$ or a matrix $D_ij$ of pair diameters.

Example:

potential = HardSpheres(1.0)

Example:

potential = HardSpheres([0.8, 0.9, 1.0])
Dij = rand(3,3)
+potential = HardSpheres(Dij)
source
OrnsteinZernike.LennardJones — Type
LennardJones

Implements the Lennard-Jones pair interaction $u(r) = 4\epsilon [(\sigma/r)^{12} - (\sigma/r)^6]$.

Expects values ϵ and σ, which respecively are the strength of the potential and particle size.

Example:

potential = LennardJones(1.0, 2.0)
source
OrnsteinZernike.PowerLaw — Type
PowerLaw

Implements the power law pair interaction $u(r) = \epsilon (\sigma/r)^{n}$.

Expects values ϵ, σ, and n, which respecively are the strength of the potential and particle size.

Example:

potential = PowerLaw(1.0, 2.0, 8)
source
OrnsteinZernike.WCADivision — Type
WCADivision{P<:Potential, T, UC}

fields

  • potential::Potential
  • cutoff : ($r_{c}$)
  • U_c : $U(r=r_c)$

Splits the potential at the cutoff point using the WCA splitting rule. This means

\[u(r) = u_{SR}(r) + U_{LR}(r),\]

where $U_{LR}(r) = u(r)$ for $r > r_{c}$, and $U(r_{c})$ for $r < r_{c}$, and $USR(r) = 0$ for $r > r_{c}$, and $U(r) - U(r_{c})$ for $r < r_{c}$.

source
diff --git a/dev/SingleCompLJ-299cf6a8.svg b/dev/SingleCompLJ-aab386ed.svg similarity index 89% rename from dev/SingleCompLJ-299cf6a8.svg rename to dev/SingleCompLJ-aab386ed.svg index 2e44812..1f69371 100644 --- a/dev/SingleCompLJ-299cf6a8.svg +++ b/dev/SingleCompLJ-aab386ed.svg @@ -1,46 +1,46 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/SingleCompLJ.html b/dev/SingleCompLJ.html index 2652536..e8fc652 100644 --- a/dev/SingleCompLJ.html +++ b/dev/SingleCompLJ.html @@ -17,7 +17,7 @@ gr = [-1.0242739989507754e-10 -1.0231282487893623e-10 … 0.9999999263517393 0.9999999263010866]' ck = [-16.004861141732388 -15.596565462363577 … 2.469275944499878e-16 -8.487600743862082e-17]' Sk = [0.09431367138921176 0.09654430029148264 … 1.0000000000000002 1.0]

Note that we have not specified the method by which we do so. In this case, a simple default method is chosen that works well in most cases.

The sol object contains fields for the radial distribution function gr, direct correlation function (in real and Fourier space) cr and ck, the static structure factor Sk, and arrays for r and k. For example, we can plot the radial distribution function as follows:

using Plots
-plot(sol.r, sol.gr, xlims=(0,5), xlabel="r", ylabel="g(r)")
Example block output

Full code:

using OrnsteinZernike
+plot(sol.r, sol.gr, xlims=(0,5), xlabel="r", ylabel="g(r)")
Example block output

Full code:

using OrnsteinZernike
 ϵ = 1.0
 σ = 1.0
 n = 8
@@ -29,4 +29,4 @@
 closure = HypernettedChain()
 sol = solve(system, closure)
 using Plots
-plot(sol.r, sol.gr, xlims=(0,5), xlabel="r", ylabel="g(r)")
+plot(sol.r, sol.gr, xlims=(0,5), xlabel="r", ylabel="g(r)") diff --git a/dev/Solvers.html b/dev/Solvers.html index 75108be..7e07ad8 100644 --- a/dev/Solvers.html +++ b/dev/Solvers.html @@ -1,2 +1,2 @@ -Solvers · OrnsteinZernike.jl Documentation
GitHub

Solvers

Having defined a SimpleLiquid and a Closure, one may choose a method by which to solve the equations. The implemented methods are Exact, FourierIteration, NgIteration. If no method is given the solve function, it will use the default NgIteration.

All solvers in some way, need to define a grid on which to solve the equations. This is done using the keyword arguments for M and dr, which represent respectively the number of grid points, and the spacing between them. Some solvers have additional settings, such as a tolerance. Using powers of 2 for M typically give the best performance (for 3D systems). For $N$ dimensional systems, the grid is constructed such that all M points lie below M*dr, but dr may not be the exact grid spacing between all points.

For the implemented cases, Exact solves the system exactly, or throws an error if the method is not implemented.

OrnsteinZernike.Exact — Type
Exact <: Method

Solves the system exactly. This is only implemented for specific systems.

Right now, the implemented exact methods are

  • three-dimensional single-component system of hard spheres with the Percus Yevick closure [1]
  • three-dimensional multi-component system of additive hard spheres with the Percus Yevick closure [2]
  • one-dimensional single-component system of hard spheres with the Percus Yevick closure [3]
  • five-dimensional single-component system of hard spheres with the Percus Yevick closure [3]

Construct using

Exact(; M=2^10, dr = sqrt(π/(M+1))/(2π))

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.

Examples

method = Exact()

method = Exact(M=1000)

method = Exact(M=1000, dr=0.01)

References:

  1. Wertheim, M. S. "Exact solution of the Percus-Yevick integral equation for hard spheres." Physical Review Letters 10.8 (1963): 321.
  2. Baxter, R. J. "Ornstein–Zernike relation and Percus–Yevick approximation for fluid mixtures." The Journal of Chemical Physics 52.9 (1970): 4559-4562.
  3. Leutheusser, E. "Exact solution of the Percus-Yevick equation for a hard-core fluid in odd dimensions." Physica A: Statistical Mechanics and its Applications 127.3 (1984): 667-676.
source

The methods FourierIteration and NgIteration both use recursive iteration to find improved estimates of the solution using Fourier Transforms. NgIteration uses a scheme to accelerate convergence, see Ref. [1].

OrnsteinZernike.FourierIteration — Type
FourierIteration <: Method

Solves the system by recursive iteration in Fourier Space. Essentially, the algorithm is:

  1. guess an initial γ(r)
  2. find c(r) using the closure relation
  3. fourier transform to get ĉ(k)
  4. find γ(k) using the OZ-eq in k-space
  5. compute γ(r) with a inverse fourier transform
  6. compare with previous value, if not converged go to 2.

Optionally, a mixing rule is used to mix the new and previous iteration of c(r) in step 2.

Arguments:

  • M::Int: number of points discretize the solution on
  • dr::Float64: grid spacing in real space
  • mixing_parameter::Float64: mixing parameter for iteration mixing. A value of 1 is no mixing. Must be between 0 and 1.
  • max_iterations::Int64: maximal number of iterations
  • 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π))

source
OrnsteinZernike.NgIteration — Type
NgIteration <: Method

Solves the system by recursive iteration in Fourier Space, and uses the Ng acceleration method. Essentially, the algorithm is:

  1. guess an initial γ(r)
  2. find c(r) using the closure relation
  3. fourier transform to get ĉ(k)
  4. find γ(k) using the OZ-eq in k-space
  5. compute γ(r) with a inverse fourier transform
  6. use Ng's method to provide a next guess for γ
  7. compare with previous value, if not converged go to 2.

Arguments:

  • M::Int: number of points discretize the solution on
  • dr::Float64: grid spacing in real space
  • N_stages::Int: Number of previous values to take into account for step 6. A higher number should lead to faster convergence, yet more computation time per iteration.
  • max_iterations::Int64: maximal number of iterations
  • 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π))

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.

source

Meta-solvers

OrnsteinZernike.DensityRamp — Type
DensityRamp <: Method

Solves the system by iteratively solving systems of increasing density, using the previous solution as initial guess at a higher density. This may help deal with convergence issues.

Arguments

  • method: method by which to solve the system for individual densities.
  • densities: densities to consider. Must be a vector of increasing values.
  • verbose: whether to print information.

Example: DensityRamp(NgIteration(), [0.1, 0.3, 0.4]; verbose=false)

source
+Solvers · OrnsteinZernike.jl Documentation
GitHub

Solvers

Having defined a SimpleLiquid and a Closure, one may choose a method by which to solve the equations. The implemented methods are Exact, FourierIteration, NgIteration. If no method is given the solve function, it will use the default NgIteration.

All solvers in some way, need to define a grid on which to solve the equations. This is done using the keyword arguments for M and dr, which represent respectively the number of grid points, and the spacing between them. Some solvers have additional settings, such as a tolerance. Using powers of 2 for M typically give the best performance (for 3D systems). For $N$ dimensional systems, the grid is constructed such that all M points lie below M*dr, but dr may not be the exact grid spacing between all points.

For the implemented cases, Exact solves the system exactly, or throws an error if the method is not implemented.

OrnsteinZernike.Exact — Type
Exact <: Method

Solves the system exactly. This is only implemented for specific systems.

Right now, the implemented exact methods are

  • three-dimensional single-component system of hard spheres with the Percus Yevick closure [1]
  • three-dimensional multi-component system of additive hard spheres with the Percus Yevick closure [2]
  • one-dimensional single-component system of hard spheres with the Percus Yevick closure [3]
  • five-dimensional single-component system of hard spheres with the Percus Yevick closure [3]

Construct using

Exact(; M=2^10, dr = sqrt(π/(M+1))/(2π))

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.

Examples

method = Exact()

method = Exact(M=1000)

method = Exact(M=1000, dr=0.01)

References:

  1. Wertheim, M. S. "Exact solution of the Percus-Yevick integral equation for hard spheres." Physical Review Letters 10.8 (1963): 321.
  2. Baxter, R. J. "Ornstein–Zernike relation and Percus–Yevick approximation for fluid mixtures." The Journal of Chemical Physics 52.9 (1970): 4559-4562.
  3. Leutheusser, E. "Exact solution of the Percus-Yevick equation for a hard-core fluid in odd dimensions." Physica A: Statistical Mechanics and its Applications 127.3 (1984): 667-676.
source

The methods FourierIteration and NgIteration both use recursive iteration to find improved estimates of the solution using Fourier Transforms. NgIteration uses a scheme to accelerate convergence, see Ref. [1].

OrnsteinZernike.FourierIteration — Type
FourierIteration <: Method

Solves the system by recursive iteration in Fourier Space. Essentially, the algorithm is:

  1. guess an initial γ(r)
  2. find c(r) using the closure relation
  3. fourier transform to get ĉ(k)
  4. find γ(k) using the OZ-eq in k-space
  5. compute γ(r) with a inverse fourier transform
  6. compare with previous value, if not converged go to 2.

Optionally, a mixing rule is used to mix the new and previous iteration of c(r) in step 2.

Arguments:

  • M::Int: number of points discretize the solution on
  • dr::Float64: grid spacing in real space
  • mixing_parameter::Float64: mixing parameter for iteration mixing. A value of 1 is no mixing. Must be between 0 and 1.
  • max_iterations::Int64: maximal number of iterations
  • 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π))

source
OrnsteinZernike.NgIteration — Type
NgIteration <: Method

Solves the system by recursive iteration in Fourier Space, and uses the Ng acceleration method. Essentially, the algorithm is:

  1. guess an initial γ(r)
  2. find c(r) using the closure relation
  3. fourier transform to get ĉ(k)
  4. find γ(k) using the OZ-eq in k-space
  5. compute γ(r) with a inverse fourier transform
  6. use Ng's method to provide a next guess for γ
  7. compare with previous value, if not converged go to 2.

Arguments:

  • M::Int: number of points discretize the solution on
  • dr::Float64: grid spacing in real space
  • N_stages::Int: Number of previous values to take into account for step 6. A higher number should lead to faster convergence, yet more computation time per iteration.
  • max_iterations::Int64: maximal number of iterations
  • 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π))

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.

source

Meta-solvers

OrnsteinZernike.DensityRamp — Type
DensityRamp <: Method

Solves the system by iteratively solving systems of increasing density, using the previous solution as initial guess at a higher density. This may help deal with convergence issues.

Arguments

  • method: method by which to solve the system for individual densities.
  • densities: densities to consider. Must be a vector of increasing values.
  • verbose: whether to print information.

Example: DensityRamp(NgIteration(), [0.1, 0.3, 0.4]; verbose=false)

source
diff --git a/dev/Systems.html b/dev/Systems.html index 4e47d7e..109a23f 100644 --- a/dev/Systems.html +++ b/dev/Systems.html @@ -18,4 +18,4 @@ ρ = [0.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.1 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.1 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.1 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.1 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.1 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.1 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.1] kBT = 1.0 potential = HardSpheres{StaticArraysCore.SMatrix{10, 10, Float64, 100}}([0.1 0.15000000000000002 0.2 0.25 0.3 0.35 0.39999999999999997 0.45 0.5 0.55; 0.15000000000000002 0.2 0.25 0.30000000000000004 0.35 0.4 0.44999999999999996 0.5 0.55 0.6; 0.2 0.25 0.3 0.35 0.4 0.44999999999999996 0.5 0.55 0.6 0.65; 0.25 0.30000000000000004 0.35 0.4 0.45 0.5 0.55 0.6000000000000001 0.65 0.7; 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75; 0.35 0.4 0.44999999999999996 0.5 0.55 0.6 0.6499999999999999 0.7 0.75 0.8; 0.39999999999999997 0.44999999999999996 0.5 0.55 0.6 0.6499999999999999 0.7 0.75 0.8 0.85; 0.45 0.5 0.55 0.6000000000000001 0.65 0.7 0.75 0.8 0.8500000000000001 0.9; 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.8500000000000001 0.9 0.95; 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1.0]) - + diff --git a/dev/Theory.html b/dev/Theory.html index 5702686..63fbcaf 100644 --- a/dev/Theory.html +++ b/dev/Theory.html @@ -1,2 +1,2 @@ -Theory · OrnsteinZernike.jl Documentation
GitHub

Theory

In this section, we describe our conventions and notation.

Single component systems

The Ornstein zernike equation $h(r) = c(r) + \rho \int d\textbf{r}' c(\textbf{r}')h(|\textbf{r} - \textbf{r}'|) $ links the pair correlation function $g(r) = h(r)+1$ to the number density $\rho$ and the direct correlation function $c(r)$. In order to solve it, a second relation between $h(r)$ and $c(r)$ must be specified. This is called the closure relation.

In practise, the closure relation typically takes the form $c(r) = f(\gamma(r), r, u(r))$, where $u(r)$ is the interaciton potential and $\gamma(r) = h(r) - c(r)$ is the indirect correlation function. Sometimes, instead the closure is defined for the bridge function $b(r)$, and the closure relation is then given by $h(r) - 1 = \exp\left(-\beta u(r) + h(r) - c(r) + b(r) \right)$, where $\beta = 1/k_BT$.

Mixtures

Everything above generalizes to the mixture case:

The Ornstein zernike equation

\[h_{ij}(r) = c_{ij}(r) + \sum_l \rho_l \int d\textbf{r}' c_{il}(\textbf{r}')h_{lj}(|\textbf{r} - \textbf{r}'|)\]

links the pair correlation function $g_{ij}(r) = h_{ij}(r)+1$ to the species specific number density $\rho_{i}$ and the direct correlation function $c_{ij}(r)$. In order to solve it, a second relation between $h_{ij}(r)$ and $c_{ij}(r)$ must be specified. This is called the closure relation.

In practise, the closure relation typically takes the form $c_{ij}(r) = f(\gamma_{ij}(r), r, u_{ij}(r))$, where $u_{ij}(r)$ is the interaciton potential and $\gamma_{ij}(r) = h_{ij}(r) - c_{ij}(r)$ is the indirect correlation function. Sometimes, instead the closure is defined for the bridge function $b(r)$, and the closure relation is then given by $h_{ij}(r) - 1 = \exp\left(-\beta u_{ij}(r) + h_{ij}(r) - c_{ij}(r) + b_{ij}(r) \right)$, where $\beta = 1/k_BT$.

Fourier Transforms

The ability to numerically solve the Ornstein-Zernike equation relies heavily on doing repeated Fourier Transforms. In arbitrary dimensions, these Fourier transforms can be written as Hankel transforms in the case that the argument is a radially symmetric function. In particular, in $d$ dimensions, the radial Fourier transform and its inverse are given by

\[\hat{F}(k) = (2\pi)^{d/2} k ^{1-d/2}\int_0^\infty dr r^{d/2}J_{d/2-1}(kr)F(r)\]

\[F(r) = (2\pi)^{-d/2} r ^{1-d/2}\int_0^\infty dk k^{d/2}J_{d/2-1}(kr)\hat{F}(k),\]

in which $J_{d/2-1}(x)$ is the bessel function of order $d/2-1$. In the special cases of 1 and 3 dimensions, the transform simplifies into:

\[\hat{F}(k) = 2\int_0^\infty dr \cos(kr)F(r)\]

\[F(r) = \frac{1}{\pi}\int_0^\infty dk \cos(kr)\hat{F}(k),\]

in 1$d$, and

\[\hat{F}(k) = \frac{4\pi}{k}\int_0^\infty dr r \sin(kr)F(r)\]

\[F(r) = \frac{1}{2\pi^2r}\int_0^\infty dk k\sin(kr)\hat{F}(k),\]

in 3$d$. This package uses discrete versions of all of the above.

Thermodynamic properties

Using the structure as determined by this package, several thermodynamic properies can be computed. In particular, this package contains methods to compute the (virial) pressure $p$, the isothermal compressibility $\chi$, and the excess internal energy per particle $E_x$.

For mixtures, they are computed respectively from the following definitions

\[p = k_BT \rho_0\sum_i x_i - 1/6 \rho_0^2 \sum_{ij} x_i x_j \int d\textbf{r} r g_{ij}(r) u'_{ij}(r)\]

\[1/(\rho_0 k_BT \chi) = 1 - ρ_0 \sum_{ij} x_i x_j \hat{c}_{ij}(k\to0)\]

,

\[E_x = 1/2 \rho_0 \sum_{ij} x_i x_j \int d\textbf{r} g_{ij}(r) u_{ij}(r)\]

in which $\rho_0=N/V$. The functions to use are compute_virial_pressure, compute_compressibility, and, compute_excess_energy.

+Theory · OrnsteinZernike.jl Documentation
GitHub

Theory

In this section, we describe our conventions and notation.

Single component systems

The Ornstein zernike equation $h(r) = c(r) + \rho \int d\textbf{r}' c(\textbf{r}')h(|\textbf{r} - \textbf{r}'|) $ links the pair correlation function $g(r) = h(r)+1$ to the number density $\rho$ and the direct correlation function $c(r)$. In order to solve it, a second relation between $h(r)$ and $c(r)$ must be specified. This is called the closure relation.

In practise, the closure relation typically takes the form $c(r) = f(\gamma(r), r, u(r))$, where $u(r)$ is the interaciton potential and $\gamma(r) = h(r) - c(r)$ is the indirect correlation function. Sometimes, instead the closure is defined for the bridge function $b(r)$, and the closure relation is then given by $h(r) - 1 = \exp\left(-\beta u(r) + h(r) - c(r) + b(r) \right)$, where $\beta = 1/k_BT$.

Mixtures

Everything above generalizes to the mixture case:

The Ornstein zernike equation

\[h_{ij}(r) = c_{ij}(r) + \sum_l \rho_l \int d\textbf{r}' c_{il}(\textbf{r}')h_{lj}(|\textbf{r} - \textbf{r}'|)\]

links the pair correlation function $g_{ij}(r) = h_{ij}(r)+1$ to the species specific number density $\rho_{i}$ and the direct correlation function $c_{ij}(r)$. In order to solve it, a second relation between $h_{ij}(r)$ and $c_{ij}(r)$ must be specified. This is called the closure relation.

In practise, the closure relation typically takes the form $c_{ij}(r) = f(\gamma_{ij}(r), r, u_{ij}(r))$, where $u_{ij}(r)$ is the interaciton potential and $\gamma_{ij}(r) = h_{ij}(r) - c_{ij}(r)$ is the indirect correlation function. Sometimes, instead the closure is defined for the bridge function $b(r)$, and the closure relation is then given by $h_{ij}(r) - 1 = \exp\left(-\beta u_{ij}(r) + h_{ij}(r) - c_{ij}(r) + b_{ij}(r) \right)$, where $\beta = 1/k_BT$.

Fourier Transforms

The ability to numerically solve the Ornstein-Zernike equation relies heavily on doing repeated Fourier Transforms. In arbitrary dimensions, these Fourier transforms can be written as Hankel transforms in the case that the argument is a radially symmetric function. In particular, in $d$ dimensions, the radial Fourier transform and its inverse are given by

\[\hat{F}(k) = (2\pi)^{d/2} k ^{1-d/2}\int_0^\infty dr r^{d/2}J_{d/2-1}(kr)F(r)\]

\[F(r) = (2\pi)^{-d/2} r ^{1-d/2}\int_0^\infty dk k^{d/2}J_{d/2-1}(kr)\hat{F}(k),\]

in which $J_{d/2-1}(x)$ is the bessel function of order $d/2-1$. In the special cases of 1 and 3 dimensions, the transform simplifies into:

\[\hat{F}(k) = 2\int_0^\infty dr \cos(kr)F(r)\]

\[F(r) = \frac{1}{\pi}\int_0^\infty dk \cos(kr)\hat{F}(k),\]

in 1$d$, and

\[\hat{F}(k) = \frac{4\pi}{k}\int_0^\infty dr r \sin(kr)F(r)\]

\[F(r) = \frac{1}{2\pi^2r}\int_0^\infty dk k\sin(kr)\hat{F}(k),\]

in 3$d$. This package uses discrete versions of all of the above.

Thermodynamic properties

Using the structure as determined by this package, several thermodynamic properies can be computed. In particular, this package contains methods to compute the (virial) pressure $p$, the isothermal compressibility $\chi$, and the excess internal energy per particle $E_x$.

For mixtures, they are computed respectively from the following definitions

\[p = k_BT \rho_0\sum_i x_i - 1/6 \rho_0^2 \sum_{ij} x_i x_j \int d\textbf{r} r g_{ij}(r) u'_{ij}(r)\]

\[1/(\rho_0 k_BT \chi) = 1 - ρ_0 \sum_{ij} x_i x_j \hat{c}_{ij}(k\to0)\]

,

\[E_x = 1/2 \rho_0 \sum_{ij} x_i x_j \int d\textbf{r} g_{ij}(r) u_{ij}(r)\]

in which $\rho_0=N/V$. The functions to use are compute_virial_pressure, compute_compressibility, and, compute_excess_energy.

diff --git a/dev/ThermodynamicConsistency-d3bf8645.svg b/dev/ThermodynamicConsistency-357dec04.svg similarity index 97% rename from dev/ThermodynamicConsistency-d3bf8645.svg rename to dev/ThermodynamicConsistency-357dec04.svg index 3162216..5cc4743 100644 --- a/dev/ThermodynamicConsistency-d3bf8645.svg +++ b/dev/ThermodynamicConsistency-357dec04.svg @@ -1,56 +1,56 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/ThermodynamicConsistency.html b/dev/ThermodynamicConsistency.html index 569bac6..73a5084 100644 --- a/dev/ThermodynamicConsistency.html +++ b/dev/ThermodynamicConsistency.html @@ -16,7 +16,7 @@ plot(solRY.r, solRY.gr, label="Rogers-Young") plot!(solPY.r, solPY.gr, label="Percus-Yevick") plot!(solHNC.r, solHNC.gr, label="Hypernetted Chain") -plot!(xlims=(0.9,2.3), xlabel="r", ylabel="g(r)")Example block output

Here, we can see, that indeed, the Rogers-Young closure finds a middle ground between the two others.

Using thermodynamics

There are a number of thermodynamic relations that allow a consistent choice of $\alpha$. Here we opt for the requirement $\frac{1}{\rho\chi_T}=\left(\frac{\partial p}{\partial \rho}\right)_T$, which physically means that we require the virial and compressibility route of obtaining the pressure to give the same results. For simplicity, we use finite differences to obtain the derivative of the pressure. See the Theory section for more details. Subsequently, we use the bisection method from Roots.jl to find the optimal parameter. The results can be compared to Table 1 in Ref. 1. 

using OrnsteinZernike, Plots
+plot!(xlims=(0.9,2.3), xlabel="r", ylabel="g(r)")
Example block output

Here, we can see, that indeed, the Rogers-Young closure finds a middle ground between the two others.

Using thermodynamics

There are a number of thermodynamic relations that allow a consistent choice of $\alpha$. Here we opt for the requirement $\frac{1}{\rho\chi_T}=\left(\frac{\partial p}{\partial \rho}\right)_T$, which physically means that we require the virial and compressibility route of obtaining the pressure to give the same results. For simplicity, we use finite differences to obtain the derivative of the pressure. See the Theory section for more details. Subsequently, we use the bisection method from Roots.jl to find the optimal parameter. The results can be compared to Table 1 in Ref. 1. 

using OrnsteinZernike, Plots
 import Roots
 function find_self_consistent_solution(ρ, kBT, M, dr, dims, pot)
 
@@ -63,4 +63,4 @@
 At ρ/√2 = 0.4, we find α = 0.26, and βp/ρ - 1 = 2.8665.
 At ρ/√2 = 0.5, we find α = 0.26, and βp/ρ - 1 = 4.7362.
 At ρ/√2 = 0.6, we find α = 0.26, and βp/ρ - 1 = 7.7674.
-At ρ/√2 = 0.654, we find α = 0.25, and βp/ρ - 1 = 10.1957.

References:

  1. Rogers, Forrest J., and David A. Young. "New, thermodynamically consistent, integral equation for simple fluids." Physical Review A 30.2 (1984): 999.
+At ρ/√2 = 0.654, we find α = 0.25, and βp/ρ - 1 = 10.1957.

References:

  1. Rogers, Forrest J., and David A. Young. "New, thermodynamically consistent, integral equation for simple fluids." Physical Review A 30.2 (1984): 999.
diff --git a/dev/index.html b/dev/index.html index e1ed16c..914dfb4 100644 --- a/dev/index.html +++ b/dev/index.html @@ -1,3 +1,3 @@ Index · OrnsteinZernike.jl Documentation
GitHub

OrnsteinZernike.jl

OrnsteinZernike.jl provides generic solvers for Ornstein-Zernike equations from liquid state theory.

Installation

To install the package run 

import Pkg
-Pkg.add("OrnsteinZernike")
+Pkg.add("OrnsteinZernike")