InverseStatMech.jl

ornstein_zernike_v(dim::Int, ρ::Float64, g2::Function, s::Function, closure::String, r_vec::AbstractVector = 0.025:0.05:10,
k_vec::AbstractVector = 0.025:0.05:20) :: Interpolations.GriddedInterpolation

Compute the Ornstein-Zernike potential ($βv(r)$) using the Ornstein-Zernike equation.

Arguments

• dim::Int: The dimensionality of the system.
• ρ::Real: The number density (particles per unit volume) of the system.
• g2::Function: The pair correlation function $g_2(r)$ as a function of $r$.
• s::Function: The structure factor $S(k)$ as a function of $k$.
• closure::String: The closure relation to be used in the Ornstein-Zernike equation. Possible values are:
  • "PMF": Potential of mean force, basically $-\ln(g_2(r))$.
  • "MFA":Mean-Field Approximation (MFA).
  • "HNC": Hypernetted Chain closure.
  • "PY": Percus-Yevick closure.

Optional Arguments

• r_vec::AbstractVector: A vector of r values at which the g_2(r) function is valid. Default is 0.025:0.05:10.
• k_vec::AbstractVector: A vector of k values at which the S(k) function is valid. Default is 0.025:0.05:20.

Returns

• An instance of Interpolations.GriddedInterpolation representing the OZ potential $βv(r)$, which can be used as a regular function.

Note: The g2 and s functions should be provided as valid functions that return the pair correlation function and structure factor, respectively, as a function of r and k.

Example:

function pair_correlation_function(r)
    # Define the pair correlation function here

    # ...
end
function structure_factor(k)
    # Define the structure factor here

    # ...
end
closure_type = "PY"
oz_v = ornstein_zernike_v(3, 0.5, pair_correlation_function, structure_factor, closure_type)

optim_parametrized_pot(my_params, pot, dim, ρ, targ_g2, targ_s; 
    large_r_grid = missing, n::Int = 600, bin_size::Float64 = 0.05, r_range::Float64 = 10, k_range::Float64 = 10,
    g2_weight_range::Float64 = 2, s_weight_range::Float64 = 4, 
    n_threads::Int = 15, configs_per_thread::Int = 10, displace = 0.2, Ψ_tol::Float64 = 0.005, show_pb::Bool = true, test::Bool = false)

Using the Torquato-Wang algorithm to perform iterative optimization of potential parameters to match target pair correlation function and structure factor.

Arguments

• my_params::Vector{Float64}: Vector of initial potential parameters.
• pot::Function: Potential function that calculates the interaction potential between particles.
• dim::Int: Dimension of the system.
• ρ::Float64: Density of the system.
• targ_g2::Function: Function representing the target pair correlation function. Accepts a distance value r and returns the target g2 value at that distance.
• targ_s::Function: Function representing the target structure factor. Accepts a wave vector k and returns the target S value at that wave vector.

Keyword Arguments (all are optional)

• large_r_grid::Missing: Large-r grid for computation of long-ranged potentials. Default value is missing.
• n::Int: Number of boxes for simulation. Default value is 600.
• bin_size::Float64: Size of the bin for pair correlation function and structure factor calculations. Default value is 0.05.
• r_range::Float64: Range of r values for pair correlation function calculation. Default value is 10.
• k_range::Float64: Range of k values for structure factor calculation. Default value is 10.
• g2_weight_range::Float64: Weight range for pair correlation function in the objective function. Default value is 2.
• s_weight_range::Float64: Weight range for structure factor in the objective function. Default value is 4.
• n_threads::Int: Number of threads for parallel computation. Default value is 15.
• configs_per_thread::Int: Number of configurations to generate per thread for simulation. Default value is 10.
• displace::Float64: Kick size in the metropolis Monte Carlo simulation. Default value is 0.2.
• Ψ_tol::Float64: Tolerance for convergence of the objective function. Default value is 0.005.
• show_pb::Bool: Boolean indicating whether to display a progress bar during simulation. Default value is true.
• test::Bool: Boolean flag to indicate whether this is a test run and return a boolean indicating convergence. Default value is false.

Returns

• If test is true, returns true if convergence is achieved, false otherwise.
• If test is false, returns the optimized potential parameters.

Example

#form of the pair potential
pot(r, params) = params[1]*exp(-(r/params[2])^2)

#initial guess parameters
my_params = [1.0, 2.0] #write 1.0 instead of 1 to indicate that this is Float64

#target pair correlation function and structure factor
targ_g2(r) = 1 
targ_s(r) = 1

#optimize the parameters
InverseStatMech.optim_parametrized_pot(my_params, pot, 2, 1, targ_g2, targ_s)

iterative_boltzmann(pot::Function, dim::Int, ρ::Float64, targ_g2::Function, α = 1; n = 500, bin_size = 0.05, r_range = 10)::Function

Iteratively updates the pair potential pot using the Boltzmann inversion method to match the target pair correlation function targ_g2.

Arguments

• pot::Function: The initial pair potential function to be optimized.
• dim::Int: The dimensionality of the system.
• ρ::Float64: The number density of particles in the system.
• targ_g2::Function: The target pair correlation function $g_2(r)$.
• α::Float64: The update parameter for the potential. (Optional, default: 1)

Keyword Arguments (All are optional)

• n: Number of particles for each box in the simulation. (default: 500)
• bin_size: Bin size for the histograms of pair correlation functions. (default: 0.05)
• r_range: Maximum distance to compute the pair correlation function. (default: 10)
• n_threads: Number of threads to use for parallel computation. (default: 15)
• configs_per_thread: Number of configurations to generate for each thread. (default: 10)
• displace: Kick size in the metropolis Monte Carlo simulation. (default: 0.2).
• Ψ_tol: Tolerance of the distance metric between the target and simulated $g_2$ for stopping criterion. (default: 0.005)
• show_pb: Whether to show the progress bar during the simulation. (default: true)
• test: Whether to run the function in test mode. (default: false)

Returns

• If test=true, returns a boolean indicating whether the optimization is successful.
• Otherwise, returns the optimized pair potential function.

Example

optimized_potential = iterative_boltzmann(r -> 0, 2, 1.0, r -> 1 - exp(-π*r^2))

reverse_mc(dim::Int, n::Int, ρ::Float64, g2_targ::Function; initial_box = missing, bin_size = 0.05, range = 5, sweeps = 100, displace = 0.1, t_i = 1, t_f = 0.001, cooling_rate = 0.98)

Reverse Monte Carlo algorithm to generate equilibrium configurations that yield a target pair correlation function $g_2(r)$.

Arguments

• dim::Int: Dimensionality of the system.
• n::Int: Number of particles.
• ρ::Float64: Number density of the particles.
• g2_targ::Function: Target pair correlation function as a function g2_targ(r).

Keyword arguments

• initial_box::Box (optional): Initial configuration box. Default is missing which generates a random box.
• bin_size::Float64 (optional): Bin size for computing the pair correlation function. Default is 0.05.
• range::Float64 (optional): Range for the interaction potential. Default is 5.
• sweeps::Int (optional): Number of Monte Carlo sweeps at each temprature. Default is 100.
• displace::Float64 (optional): Maximum displacement for particle moves. Default is 0.1.
• t_i::Float64 (optional): Initial temperature. Default is 1.
• t_f::Float64 (optional): Final temperature. Default is 0.001.
• cooling_rate::Float64 (optional): Cooling rate for temperature reduction. Default is 0.98.

Returns

• b::Box: Generated equilibrium classical configuration box. Use b.particles' to get the particle positions.

Example

box = InverseStatMech.reverse_mc(2, 100, 0.5, r -> 1 - exp(-π*r^2))