InverseStatMech.jl
Efficient inverse statistical mechanical algorithms to generate effective potentials or configurations from pair correlation functions or structure factors.
Package Features
- Torquato-Wang algorithm (Torquato and Wang, PRE 2020) to find optimized parametrized potentials given target $g_2(r)$ and $S(k)$
- Iterative Boltzmann Inversion (Soper, Chem. Phys. 1996) to find optimimzed binned potentials given target $g_2(r)$
- Reverse Monte-Carlo algorithm (McGreevy, J. Phys. Cond. Matter 2001) to find configuration that match target $g_2(r)$
- More to come!
Function Documentation
InverseStatMech.optim_parametrized_pot
— Functionoptim_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_box::Int = 10, Ψ_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 valuer
and returns the target g2 value at that distance.targ_s::Function
: Function representing the target structure factor. Accepts a wave vectork
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 ismissing
.n::Int
: Number of boxes for simulation. Default value is600
.bin_size::Float64
: Size of the bin for pair correlation function and structure factor calculations. Default value is0.05
.r_range::Float64
: Range of r values for pair correlation function calculation. Default value is10
.k_range::Float64
: Range of k values for structure factor calculation. Default value is10
.g2_weight_range::Float64
: Weight range for pair correlation function in the objective function. Default value is2
.s_weight_range::Float64
: Weight range for structure factor in the objective function. Default value is4
.n_threads::Int
: Number of threads for parallel computation. Default value is15
.configs_per_box::Int
: Number of configurations per box for simulation. Default value is10
.Ψ_tol::Float64
: Tolerance for convergence of the objective function. Default value is0.005
.show_pb::Bool
: Boolean indicating whether to display a progress bar during simulation. Default value istrue
.test::Bool
: Boolean flag to indicate whether this is a test run and return a boolean indicating convergence. Default value isfalse
.
Returns
- If
test
is true, returnstrue
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)
InverseStatMech.iterative_boltzmann
— Functioniterative_boltzmann(pot, dim, ρ, targ_g2, α = 1; n = 500, bin_size = 0.05, r_range = 10)
Iteratively updates the pair potential pot
using the Boltzmann inversion method to match the target pair correlation function targ_g2
.
Arguments
pot
: The initial pair potential function to be optimized.dim
: The dimensionality of the system.ρ
: The number density of particles in the system.targ_g2
: The target pair correlation functiong_2(r)
.α
: The update parameter for the potential. (Optional, default: 1)
Keyword Arguments (All are optional)
n
: Number of configurations 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_box
: Number of configurations to generate for each thread. (default: 10)Ψ_tol
: Tolerance 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 -> exp(-π*r^2))
InverseStatMech.reverse_mc
— Functionreverse_mc(dim, n, ρ, g2_targ; 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 functiong2_targ(r)
.
Keyword arguments
initial_box::Box
(optional): Initial configuration box. Default ismissing
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. Useb.particles'
to get the particle positions.
Example
box = InverseStatMech.reverse_mc(2, 100, 0.5, r -> 1 - exp(π*r^2))