Package potentials_from_particle_insertion
Description
Python package with a set of functions for calculating the g(r) using the Widom test-particle insertion method, and/or for fitting the g(r) to solve for the pair-potential using particle coordinates from e.g. confocal microscopy or cryo-TEM.
The codes for this package were developed as part of the work in reference [1], where the iterative use of test-particle insertion to solve for the pair potential was based on [2] with correction for periodic or finite boundary conditions based on equations described in [3,4]
References
[1] Bransen, M. (2024). Measuring interactions between colloidal (nano)particles. PhD thesis, Utrecht University.
[2] Stones, A. E., Dullens, R. P. A., & Aarts, D. G. A. L. (2019). Model-Free Measurement of the Pair Potential in Colloidal Fluids Using Optical Microscopy. Physical Review Letters, 123(9), 098002. https://doi.org/10.1103/PhysRevLett.123.098002
[3] Markus Seserno (2014). How to calculate a three-dimensional g(r) under periodic boundary conditions. https://www.cmu.edu/biolphys/deserno/pdf/gr_periodic.pdf
[4] Kopera, B. A. F., & Retsch, M. (2018). Computing the 3D Radial Distribution Function from Particle Positions: An Advanced Analytic Approach. Analytical Chemistry, 90(23), 13909–13914. https://doi.org/10.1021/acs.analchem.8b03157
License
MIT license
Copyright (c) 2024 Maarten Bransen
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.
Sub-modules
potentials_from_particle_insertion.multicomponent
-
This file contains codes for multicomponent compatible versions of functions, i.e. for more than one 'type' of particle …
Functions
def load_TPI_results(filename)
-
loads text file containing TPI results as written by
save_tpi_results()
Parameters
filename
:str
- name of the TPI result file to load
Returns
rmin
:float
- lower bound on interparticle distance
rmax
:float
- upper bound on interparticle distance
dr
:float
- bin width of interparticle distance
bincentres
:list
offloat
- centre position of each interparticle distance bin
dhrdf
:list
offloat
- g(r) values from DH calculation for each bin
rdfs
:list
oflist
offloat
- g(r) values from TPI for each bin for each iteration
potentials
:list
oflist
offloat
- u(r) values for each bin for each iteration
counts
:list
oflist
offloat
- bin counts for each bin for each iteration
errors
:list
offloat
- chi squared error for each iteration
def rdf_dist_hist_2d(coordinates, rmin=0, rmax=10, dr=None, handle_edge='rectangle', boundary=None, density=None, quiet=False, neighbors_upper_bound=None, workers=1)
-
calculates g(r) via a 'conventional' distance histogram method for a set of 2D coordinate sets. Provided for convenience.
Parameters
coordinates
:numpy.array
orlist-like
ofnumpy.array
- list of sets of coordinates, where each item along the 0th dimension is
a N*2 numpy.array of particle coordinates, where each array is an
independent set of coordinates (e.g. one TEM image, a time step from a
video, etc.), with each element of the array of form
[y,x]
. Each set of coordinates is not required to have the same number of particles but must have the same particle density. When no boundaries or a single boundary is given, it is assumed all coordinate sets share these same boundaries. rmin
:float
, optional- lower bound for the pairwise distance, left edge of 0th bin. The default is 0.
rmax
:float
, optional- upper bound for the pairwise distance, right edge of last bin. The default is 10.
dr
:float
, optional- bin width for the pairwise distance bins. The default is (rmax-rmin)/20
handle_edge
:one
of['none','rectangle','periodic rectangle','circle', (list of) callable]
, optional-
specifies how to correct for edge effects in the radial distribution function. These edge effects occur due to particles closer than
rmax
to any of the boundaries missing part of their neighbour shells in the dataset. The following options for correcting for this are available:'none'
orNone
orFalse
(all equivalent): do not correct for edge effects at all. Note that in order to calculate the particle density, cuboidal boundaries are assumed even whenboundary
is not specified. This can be overridden by explicitely giving the particle density using thedensity
parameter.'rectangle'
: correct for the missing volume in a square or rectangular boundary.'periodic rectangle'
: like'rectangle'
, except in periodic boundary conditions (i.e. one side wraps around to the other), based on ref [1].'circle'
: correct for missing volume in a spherical boundary.- a custom callable function (or list thereof) can be given to
correct for arbitrary and/or mixed boundary conditions. This
function must take three arguments: a numpy array of bin edges, an
N×2 numpy array of coordinates and a boundary as specified in
boundary
, and return an N ×len(bin edges)-1
numpy array with a value between 0 and 1 specifying the fraction of the volume of each circular shell that is within the boundary.
The default is 'rectangle'.
boundary
:array-like
, optional-
positions of the walls that define the bounding box of the coordinates. The format must be given as
((ymin,ymax),(xmin,xmax))
whenhandle_edge
is'none'
,'rectangle'
or'periodic rectangle'
. Whenhandle_edge='circle'
the coordinates of the origin and radius of the bouding circle must be given as(y,x,radius)
. For custom boundary handling,boundary
can be any format as required by the edge edge correction function(s) passed tohandle_edge
.If all coordinate sets share the same boundary, a single such boundary can be given, otherwise a list of such array-likes of the same length as coordinates can be given for specifying boundaries of each set in
coordinates
separately. The defaultNone
which determines the smallest set of boundaries encompassing a set of coordinates. density
:float
orcallable
orlist
ofcallable
, optional- number density of particles in the box to use for normalizing the
values. The default is the average density based on
coordinates
andboundary
. When using custom edge correction inhandle_edge
it is possible to instead pass a (list of) callable that takes the number of particles and boundary and returns the number density. quiet
:bool
, optional- if True, no output is printed to the terminal by this function call. The default is False.
neighbors_upper_bound
:int
, optional- upper limit on the number of neighbors expected within rmax from a particle. Useful for reducing memory consumption in datasets with dimensions much larger than rmax. The default is None, which takes the total number of particles in the coordinate set as upper limit.
workers
:int
, optional- number of workers to use for parallel processing during the neighbour
detection step. If -1 is given all CPU threads are used. Note: this is
ignored when
handle_edge='periodic rectangle'
. The default is 1.
Returns
rvals
:numpy.array
- bin-edges of the radial distribution function.
bincounts
:numpy.array
- values for the bins of the radial distribution function
References
[1] Markus Seserno (2014). How to calculate a three-dimensional g(r) under periodic boundary conditions. https://www.cmu.edu/biolphys/deserno/pdf/gr_periodic.pdf
def rdf_dist_hist_3d(coordinates, rmin=0, rmax=10, dr=None, handle_edge='cuboid', boundary=None, density=None, quiet=False, neighbors_upper_bound=None, workers=1)
-
calculates g(r) via a 'conventional' distance histogram method for a set of 3D coordinate sets. Provided for convenience.
Parameters
coordinates
:numpy.array
orlist-like
ofnumpy.array
- list of sets of coordinates, where each item along the 0th dimension is
a N*3 numpy.array of particle coordinates, where each array is an
independent set of coordinates (e.g. one z-stack, a time step from a
video, etc.), with each element of the array of form
[z,y,x]
. Each set of coordinates is not required to have the same number of particles but must have the same particle density. When no boundaries or a single boundary is given, it is assumed all coordinate sets share these same boundaries. rmin
:float
, optional- lower bound for the pairwise distance, left edge of 0th bin. The default is 0.
rmax
:float
, optional- upper bound for the pairwise distance, right edge of last bin. The default is 10.
dr
:float
, optional- bin width for the pairwise distance bins. The default is (rmax-rmin)/20
handle_edge
:one
of['none','cuboid','periodic cuboid','sphere', (list of) callable]
, optional-
specifies how to correct for edge effects in the radial distribution function. These edge effects occur due to particles closer than
rmax
to any of the boundaries missing part of their neighbour shells in the dataset. The following options for correcting for this are available:'none'
orNone
orFalse
(all equivalent): do not correct for edge effects at all. Note that in order to calculate the particle density, cuboidal boundaries are assumed even whenboundary
is not specified. This can be overridden by explicitely giving the particle density usingdensity
.'cuboid'
: correct for the missing volume in cuboidal boundary conditions, e.g. a 3D rectangular box with right angles. Based on ref. [1]'periodic cuboid'
: like'cuboid'
, except in periodic boundary conditions (i.e. one side wraps around to the other). Based on ref. [2].'sphere'
: correct for missing volume in spherical boundary conditions- a custom callable function (or list thereof) can be given to
correct for arbitrary and/or mixed boundary conditions. This
function must take three arguments: a numpy array of bin edges, an
N×3 numpy array of coordinates and a boundary as specified in
boundary
, and return an N ×len(bin edges)-1
numpy array with a value between 0 and 1 specifying the fraction of the volume of each spherical shell that is within the boundary.
The default is 'cuboid'.
boundary
:array-like
, optional-
positions of the walls that define the bounding box of the coordinates. The format must be given as
((zmin,zmax),(ymin,ymax),(xmin,xmax))
whenhandle_edge
is'none'
,'cuboid'
or'periodic cuboid'
. Whenhandle_edge='sphere'
the coordinates of the origin and radius of the bouding sphere must be given as(z,y,x,radius)
.If all coordinate sets share the same boundary a single such boundary can be given, otherwise a list of such array-likes of the same length as coordinates can be given for specifying boundaries of each set in
coordinates
separately. The default isNone
which determines the smallest set of boundaries encompassing a set of coordinates. density
:float
orcallable
orlist
ofcallable
, optional- number density of particles in the box to use for normalizing the
values. The default is the average density based on
coordinates
andboundary
. When using custom edge correction inhandle_edge
it is possible to instead pass a (list of) callable that takes the number of particles and boundary and returns the number density. quiet
:bool
, optional- if True, no output is printed to the terminal by this function call. The default is False.
neighbors_upper_bound
:int
, optional- upper limit on the number of neighbors expected within rmax from a particle. Useful for reducing memory consumption in datasets with dimensions much larger than rmax. The default is None, which takes the total number of particles in the coordinate set as upper limit.
workers
:int
, optional- number of workers to use for parallel processing during the neighbour
detection step. If -1 is given all CPU threads are used. Note: this is
ignored when
handle_edge='periodic cuboid'
. The default is 1.
Returns
rvals
:numpy.array
- bin-edges of the radial distribution function.
bincounts
:numpy.array
- values for the bins of the radial distribution function
References
[1] Kopera, B. A. F., & Retsch, M. (2018). Computing the 3D Radial Distribution Function from Particle Positions: An Advanced Analytic Approach. Analytical Chemistry, 90(23), 13909–13914. https://doi.org/10.1021/acs.analchem.8b03157
[2] Markus Seserno (2014). How to calculate a three-dimensional g(r) under periodic boundary conditions. https://www.cmu.edu/biolphys/deserno/pdf/gr_periodic.pdf
def rdf_insertion_binned_2d(coordinates, pairpotential, rmin=0, rmax=10, dr=None, handle_edge='rectangle', boundary=None, pairpotential_binedges=None, n_ins=1000, interpolate=True, avoid_boundary=False, avoid_coordinates=False, neighbors_upper_bound=None, workers=1, testparticle_func=None, **kwargs)
-
Calculate g(r) from insertion of test-particles into sets of existing 2D coordinates, averaged over bins of width dr, and based on the pairwise interaction potential u(r) (in units of kT).
Parameters
coordinates
:numpy.array
orlist-like
ofnumpy.array
- list of sets of coordinates, where each item along the 0th dimension is
a N*2 numpy.array of particle coordinates, where each array is an
independent set of coordinates (e.g. one TEM image, a time step from a
video, etc.), with each element of the array of form
[y,x]
. Each set of coordinates is not required to have the same number of particles but must have the same particle density. When no boundaries or a single boundary is given, it is assumed all coordinate sets share these same boundaries. pairpotential
:iterable
- list of values for the pairwise interaction potential. Must have length
of
len(pairpotential_binedges)-1
and be in units of thermal energy kT rmin
:float
, optional- lower bound for the pairwise distance, left edge of 0th bin. The default is 0.
rmax
:float
, optional- upper bound for the pairwise distance, right edge of last bin. The default is 10.
dr
:float
, optional- bin width for the pairwise distance bins. The default is (rmax-rmin)/20
handle_edge
:one
of['none','rectangle','periodic rectangle','circle', (list of) callable]
, optional-
specifies how to correct for edge effects in the radial distribution function. These edge effects occur due to particles closer than
rmax
to any of the boundaries missing part of their neighbour shells in the dataset. The following options for correcting for this are available:'rectangle'
: correct for the missing volume in a square or rectangular boundary.'periodic' rectangle'
: like'rectangle'
, except in periodic boundary conditions (i.e. one side wraps around to the other), based on ref [1].'circle'
: correct for missing volume in a spherical boundary.- a custom callable function (or list thereof) can be given to
correct for arbitrary and/or mixed boundary conditions. This
function must take three arguments: a numpy array of bin edges, an
N×2 numpy array of coordinates and a boundary as specified in
boundary
, and return an N ×len(bin edges)-1
numpy array with a value between 0 and 1 specifying the fraction of the volume of each spherical shell that is within the boundary.
The default is 'rectangle'.
boundary
:array-like
orlist
ofsuch
, optional-
positions of the walls that define the bounding box of the coordinates. The format must be given as
((ymin,ymax),(xmin,xmax))
whenhandle_edge
is`'rectangle'
or'periodic rectangle'
. Whenhandle_edge='circle'
the coordinates of the origin and radius of the bouding circle must be given as(y,x,radius)
. For custom boundary handling,boundary
can be any format as required by the edge edge correction function(s) passed tohandle_edge
.If all coordinate sets share the same boundary, a single such boundary can be given, otherwise a list of such array-likes of the same length as coordinates can be given for specifying boundaries of each set in
coordinates
separately. The default isNone
which determines the smallest set of boundaries encompassing a set of coordinates. pairpotential_binedges
:iterable
, optional- bin edges corresponding to the values in
pairpotential
. The default is None, which uses the bins defined byrmin
,rmax
anddr
. n_ins
:int
, optional- the number of test-particles to insert into each item in
coordinates
. The default is 1000. interpolate
:bool
, optional- whether to use linear interpolation for calculating the interaction of
two particles using the values in
pairpotential
. The default is True. If False, the nearest bin value is used. avoid_boundary
:bool
, optional- if True, all test-particles are inserted at least
rmax
away from any of the surfaces defined inboundary
to avoid effects of the finite volume of the bounding box. The default is False, which uses an analytical correction factor for missing volume of test-particles near the boundaries. This parameter is ignored when using custom edge handling. avoid_coordinates
:bool
, optional- whether to insert test-particles at least
rmin
away from the center of any of the 'real' coordinates incoordinates
. The default is False. neighbors_upper_bound
:int
, optional- upper limit on the number of neighbors expected within rmax from a particle. Useful for reducing memory consumption in datasets with dimensions much larger than rmax. The default is None, which takes the total number of particles in the coordinate set as upper limit.
workers
:int
, optional- number of workers to use for parallel processing during the neighbour
detection step. If -1 is given all CPU threads are used. Note: this is
ignored when
handle_edge='periodic rectangle'
. The default is 1. testparticle_func
:callable
orlist thereof
- function that generates test-particle coordinates in the bounding box,
required when using custom edge handling (and ignored otherwise). When
avoid_coordinates=True
this function must take 4 arguments: a boundary as given inboundary
, a N×2 numpy.array of coordinates to avoid, a float specifing the radius around the coordinates to avoid, and an integer giving the number of coordinates to generate. Otherwise, it must only take the boundary and the number of coordinates to generate as arguments. It must return ann_ins
× 2 numpy array of coordinates.
Returns
rvals
:numpy.array
- bin-edges of the radial distribution function.
pair_correlation
:numpy.array
- values for the rdf / pair correlation function in each bin
counter
:numpy.array
- number of pair counts that contributed to the (mean) values in each bin
References
[1] Markus Seserno (2014). How to calculate a three-dimensional g(r) under periodic boundary conditions. https://www.cmu.edu/biolphys/deserno/pdf/gr_periodic.pdf
def rdf_insertion_binned_3d(coordinates, pairpotential, rmin=0, rmax=10, dr=None, handle_edge='cuboid', boundary=None, pairpotential_binedges=None, n_ins=1000, interpolate=True, avoid_boundary=False, avoid_coordinates=False, neighbors_upper_bound=None, workers=1, testparticle_func=None, **kwargs)
-
Calculate g(r) from insertion of test-particles into sets of existing 2D coordinates, averaged over bins of width dr, and based on the pairwise interaction potential u(r) (in units of kT).
Parameters
coordinates
:numpy.array
orlist-like
ofnumpy.array
- list of sets of coordinates, where each item along the 0th dimension is
a N*3 numpy.array of particle coordinates, where each array is an
independent set of coordinates (e.g. one z-stack, a time step from a
video, etc.), with each element of the array of form
[z,y,x]
. Each set of coordinates is not required to have the same number of particles but must have the same particle density. When no boundaries or a single boundary is given, it is assumed all coordinate sets share these same boundaries. pairpotential
:iterable
- list of values for the pairwise interaction potential. Must have length
of
len(pairpotential_binedges)-1
and be in units of thermal energy kT rmin
:float
, optional- lower bound for the pairwise distance, left edge of 0th bin. The default is 0.
rmax
:float
, optional- upper bound for the pairwise distance, right edge of last bin. The default is 10.
dr
:float
, optional- bin width for the pairwise distance bins. The default is (rmax-rmin)/20
handle_edge
:one
of['none','cuboid','periodic cuboid','sphere', (list of) callable]
, optional-
specifies how to correct for edge effects in the radial distribution function. These edge effects occur due to particles closer than
rmax
to any of the boundaries missing part of their neighbour shells in the dataset. The following options for correcting for this are available:'none'
orNone
orFalse
(all equivalent): do not correct for edge effects at all. Note that in order to calculate the particle density, cuboidal boundaries are assumed even whenboundary
is not specified. This can be overridden by explicitely giving the particle density usingdensity
.'cuboid'
: correct for the missing volume in cuboidal boundary conditions, e.g. a 3D rectangular box with right angles. Based on ref. [1]'periodic' cuboid'
: like'cuboid'
, except in periodic boundary conditions (i.e. one side wraps around to the other). Based on ref. [2].'sphere'
: correct for missing volume in spherical boundary conditions- a custom callable function (or list thereof) can be given to
correct for arbitrary and/or mixed boundary conditions. This
function must take three arguments: a numpy array of bin edges, an
N×3 numpy array of coordinates and a boundary as specified in
boundary
, and return an N ×len(bin edges)-1
numpy array with a value between 0 and 1 specifying the fraction of the volume of each spherical shell that is within the boundary.
The default is 'cuboid'.
boundary
:array-like
, optional-
positions of the walls that define the bounding box of the coordinates. The format must be given as
((zmin,zmax),(ymin,ymax),(xmin,xmax))
whenhandle_edge
is'none'
,'cuboid'
or'periodic cuboid'
. Whenhandle_edge='sphere'
the coordinates of the origin and radius of the bouding sphere must be given as(z,y,x,radius)
.If all coordinate sets share the same boundary a single such boundary can be given, otherwise a list of such array-likes of the same length as coordinates can be given for specifying boundaries of each set in
coordinates
separately. The default isNone
which determines the smallest set of boundaries encompassing a set of coordinates. pairpotential_binedges
:iterable
, optional- bin edges corresponding to the values in
pairpotential
. The default is None, which uses the bins defined byrmin
,rmax
anddr
. n_ins
:int
, optional- the number of test-particles to insert into each item in
coordinates
. The default is 1000. interpolate
:bool
, optional- whether to use linear interpolation for calculating the interaction of
two particles using the values in
pairpotential
. The default is True. If False, the nearest bin value is used. avoid_boundary
:bool
, optional- if True, all test-particles are inserted at least
rmax
away from any of the surfaces defined inboundary
to avoid effects of the finite volume of the bounding box. The default is False, which uses an analytical correction factor for missing volume of test-particles near the boundaries. This parameter is ignored when using custom edge handling. avoid_coordinates
:bool
, optional- whether to insert test-particles at least
rmin
away from the center of any of the 'real' coordinates incoordinates
. The default is False. neighbors_upper_bound
:int
, optional- upper limit on the number of neighbors expected within rmax from a particle. Useful for reducing memory consumption in datasets with dimensions much larger than rmax. The default is None, which takes the total number of particles in the coordinate set as upper limit.
workers
:int
, optional- number of workers to use for parallel processing during the neighbour
detection step. If -1 is given all CPU threads are used. Note: this is
ignored when
handle_edge='periodic rectangle'
. The default is 1. testparticle_func
:callable
orlist thereof
- function that generates test-particle coordinates in the bounding box,
required when using custom edge handling (and ignored otherwise). When
avoid_coordinates=True
this function must take 4 arguments: a boundary as given inboundary
, a N×3 numpy.array of coordinates to avoid, a float specifing the radius around the coordinates to avoid, and an integer giving the number of coordinates to generate. Otherwise, it must only take the boundary and the number of coordinates to generate as arguments. It must return ann_ins
×3 numpy array of coordinates.
Returns
rvals
:numpy.array
- bin-edges of the radial distribution function.
pair_correlation
:numpy.array
- values for the rdf / pair correlation function in each bin
counter
:numpy.array
- number of pair counts that contributed to the (mean) values in each bin
References
[1] Kopera, B. A. F., & Retsch, M. (2018). Computing the 3D Radial Distribution Function from Particle Positions: An Advanced Analytic Approach. Analytical Chemistry, 90(23), 13909–13914. https://doi.org/10.1021/acs.analchem.8b03157
[2] Markus Seserno (2014). How to calculate a three-dimensional g(r) under periodic boundary conditions. https://www.cmu.edu/biolphys/deserno/pdf/gr_periodic.pdf
def rdf_insertion_exact_3d(coordinates, pairpotential, rmax, dr, boundary, pairpotential_binedges=None, gen_prob_reps=1000, shell_prob_reps=10, interpolate=True, use_numba=True, rmin=0)
-
Warning
This function underwent very minimal testing and is computationally inefficient. Use at your own risk.
calculate g(r) from particle insertion method using particle coordinates and pairwise interaction potential u(r) (in units of kT). Inserts test- particles at a specific r for every real particle. Implementation based on ref [1]
Parameters
coordinates
:list
ofnumpy.array
offloats
ofshape nt*n*3
- The coordinates of particles in the system
pairpotential
:list
offloats
- the value of the interparticle potential in each bin of 0 to rmax+dr in steps of dr. The value is assumed to correspond to the centrepoint of the bins.
rmax
:float
- cut-off value for the radius. Interactions beyond this range are considered negligable. Right edge of last bin in the pair potential
dr
:float
- bin width/step size of the interparticle distance used for g(r) and the pair potential.
boundary
:tuple
offloats
ofform ((zmin,zmax),(ymin,ymax),(xmin,xmax))
- defines the boundaries of the box in which coordinates may exist.
gen_prob_reps
:int
, optional- number of trial particles used to evaluate the general probability of
placing a particle at random coordinates in the box defined by
boundary
. The default is 1000. shell_prob_reps
:int
, optional- number of trial particles per reference coordinate in
coordinates
used to evaluate the probability of placing a particle in each distance bin. The default is 10. interpolate
:bool
, optional- If true, the pair potential is linearly interpolated between bin centres to calculate the interactions between all pairs of particles. If false, the value from the nearest bin is taken which is slower to compute. The default is True.
Returns
pair_correlation
:list
offloat
- the pair correlation / radial distribution functions evaluated in the bins whose edges are defined by numpy.arange(0,rmax+dr,dr)
counters
:list
ofint
- The number of trialparticles evaluated for each distance bin
References
[1] Stones, A. E., Dullens, R. P. A., & Aarts, D. G. A. L. (2018). Contact values of pair distribution functions in colloidal hard disks by test-particle insertion. The Journal of Chemical Physics, 148(24), 241102. https://doi.org/10.1063/1.5038668
def run_iteration(coordinates, pair_correlation_func, initial_guess=None, rmin=0, rmax=10, dr=None, convergence_tol=1e-05, max_iterations=100, zero_clip=1e-20, regulate=False, **kwargs)
-
Run the algorithm to solve for the pairwise potential that most accurately reproduces the radial distribution function using test-particle insertion, based on the ideas in ref. [1].
Parameters
coordinates
:list-like
ofnumpy.array
- List of sets of coordinates, where each item along the 0th dimension is
a n*3 numpy.array of particle coordinates, where each array is an
independent set of coordinates (e.g. one z-stack, a time step from a
video, etc.), with each element of the array of form
[z,y,x]
or[y,x]
in case of 2D data. Each set of coordinates is not required to have the same number of particles but all stacks must share the same bounding box as given byboundary
, and all coordinates must be within this bounding box. pair_correlation_func
:list
offloat
- bin values for the true pair correlation function that the algorithm will try to match iteratively.
initial_guess
:list
offloat
, optional- Initial guess for the particle potential on the 0th iteration. The default is None which gives 0 in each bin.
rmin
:float
, optional- left edge of the smallest bin in interparticle distance r to consider. The default is 0.
rmax
:float
, optional- Right edge of the largest bin in interparticle distance r to consider. The default is 20.
dr
:float
, optional- Stepsize or bin width in interparticle distance r. The default is 0.5.
convergence_tol
:float
, optional- target value for χ², if it dips below this value the iteration is
considered to be converged and ended. The default is
1e-5
. max_iterations
:int
, optional- Maximum number of iterations after which the algorithm is ended. The default is 100.
zero_clip
:float
, optional- values below the value of zero-clip are set to this value to avoid
devision by zero errors. The default is
1e-20
. regulate
:bool
, optional- if True, use regularization to more gently nudge towards the input g(r)
at the cost of slower convergence. Experimental option. The default is
False
. **kwargs
:key=value
- Additional keyword arguments are passed on to
rdf_insertion_binned_2d()
orrdf_insertion_binned_3d()
Returns
χ²
:list
offloat
- summed squared error in the pair correlation function for each iteration
pairpotential
:list
oflist
offloat
- the values for the pair potential in each bin for each iteration
paircorrelation
:list
oflist
offloat
- the values for the pair correlation function from test-particle insertion for each iteration
counts
:list
oflist
ofint
- number of pair counts contributing to each bin in each iteration
References
[1] Stones, A. E., Dullens, R. P. A., & Aarts, D. G. A. L. (2019). Model- Free Measurement of the Pair Potential in Colloidal Fluids Using Optical Microscopy. Physical Review Letters, 123(9), 098002. https://doi.org/10.1103/PhysRevLett.123.098002
See Also
rdf_insertion_binned_2d()
- 2D routine for g(r) from test-particle insertion
rdf_insertion_binned_3d()
- 3D routine for g(r) from test-particle insertion
def run_iterator_fitfunction(coordinates, pair_correlation_func, potential_func, initial_guess=None, fit_bounds=None, rmin=0, rmax=20, dr=None, convergence_tol=1e-05, max_iterations=100, max_func_evals=100, zero_clip=1e-20, **kwargs)
-
Run the algorithm to solve for the pairwise potential that most accurately reproduces the radial distribution function using test-particle insertion, using a fit function of known form instead of a generalized binned potential based on scipy.optimize.curve_fit.
Parameters
coordinates
:list-like
ofnumpy.array
- List of sets of coordinates, where each item along the 0th dimension is
a n*3 numpy.array of particle coordinates, where each array is an
independent set of coordinates (e.g. one z-stack, a time step from a
video, etc.), with each element of the array of form
[z,y,x]
or[y,x]
in case of 2D data. Each set of coordinates is not required to have the same number of particles but all stacks must share the same bounding box as given byboundary
, and all coordinates must be within this bounding box. pair_correlation_func
:list
offloat
- bin values for the true pair correlation function that the algorithm will try to match iteratively.
potential_func
:callable
- fit funtion to use, where the first argument is the interparticle distance r, and any subsequent arguments are optimized for.
initial_guess
:list
, optional- values for the function parameters to use for the first iteration. The default is None.
fit_bounds
:2-tuple
oflist_like
, optional- Each element of the tuple must be either an array with the length equal to the number of parameters, or a scalar (in which case the bound is taken to be the same for all parameters). Use np.inf with an appropriate sign to disable bounds on all or some parameters.
rmin
:float
, optional- left edge of the smallest bin in interparticle distance r to consider. The default is 0.
rmax
:float
, optional- Right edge of the largest bin in interparticle distance r to consider. The default is 20.
dr
:float
, optional- Stepsize or bin width in interparticle distance r. The default is 0.5.
convergence_tol
:float
, optional- target value for χ², if it dips below this value the iteration is
considered to be converged and ended. The default is
1e-5
. max_iterations
:int
, optional- Maximum number of iterations after which the algorithm is ended. The default is 100.
max_func_evals
:int
, optional- Maximum number of times the TPI function is evaluated. The default is 100.
zero_clip
:float
, optional- values below the value of zero-clip are set to this value to avoid
devision by zero errors. The default is
1e-20
. **kwargs
:key=value
- Additional keyword arguments are passed on to
rdf_insertion_binned_2d()
orrdf_insertion_binned_3d()
Returns
χ²
:list
offloat
- summed squared error in the pair correlation function for each iteration
pairpotential
:list
ofcallable
- potential functions giving the potential for each iteraction
fit_params
:list
oflist
offloat
- the values for the function parameters in each iteration
paircorrelation
:list
oflist
offloat
- the values for the pair correlation function from test-particle insertion for each iteration
counts
:list
oflist
ofint
- number of pair counts contributing to each bin in each iteration
See Also
rdf_insertion_binned_2d()
- 2D routine for g(r) from test-particle insertion
rdf_insertion_binned_3d()
- 3D routine for g(r) from test-particle insertion
def save_TPI_results(rmin, rmax, dr, bincentres, dhrdf, rdfs, potentials, counts, errors, filename='particle_insertion_result.txt')
-
utility function for saving results from TPI as txt file
Parameters
rmin
:float
- lower bound on interparticle distance
rmax
:float
- upper bound on interparticle distance
dr
:float
- bin width of interparticle distance
bincentres
:list
offloat
- centre position of each interparticle distance bin
dhrdf
:list
offloat
- g(r) values from DH calculation for each bin
rdfs
:list
oflist
offloat
- g(r) values from TPI for each bin for each iteration
potentials
:list
oflist
offloat
- u(r) values for each bin for each iteration
counts
:list
oflist
offloat
- bin counts for each bin for each iteration
errors
:list
offloat
- chi squared error for each iteration
filename
:str
, optional- filename to use. The default is 'particle_insertion_result.txt'.