Package `trajectories_to_forces`

Description

Python package with a set of functions for extracting pairwise interaction forces from particle trajectories under overdamped (Brownian) or inertial (molecular) dynamics.

The codes for this package were developed as part of the work in reference [1], based on the concepts in reference [2].

References

[1] Bransen, M. (2024). Measuring interactions between colloidal (nano)particles. PhD thesis, Utrecht University.

[2] Jenkins, I. C., Crocker, J. C., & Sinno, T. (2015). Interaction potentials from arbitrary multi-particle trajectory data. Soft Matter, 11(35), 6948–6956. https://doi.org/10.1039/C5SM01233C

License

MIT license

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.

Functions

def convert_cylindrical_force_axes(binmean_z, binmean_rho, G_z, G_rho)

convert the force coefficients along rho and z axes to the absolute value and force as well as the coefficient values parallel and perpendicular to the particle-particle vector.

Parameters

binmean_z : np.ndarray of float: count weighted mean z distance for all pairs in each bin.
binmean_rho : np.ndarray of float: count weighted mean rho distance for all pairs in each bin.
G_z : np.ndarray of float: force coefficients along z.
G_rho : np.ndarray of float: force coefficients along rho.

Returns

G_abs : np.ndarray of float: absolute value of the force coefficients, i.e. net force scalar.
G_ang : np.ndarray of float: (clockwise) angle of the force vector w.r.t. the z axis in radians.
G_para : np.ndarray of float: force coefficients for the force parallel to the particle-particle vector.
G_perp : np.ndarray of float: force coefficients for the force perpendicular to the particle-particle vector.

def filter_msd(coordinates, times=None, pos_cols=('z', 'y', 'x'), msd_min=0.01, msd_max=1, interval=1)

calculates the squared displacement between subsequent frames for each particle individually, meaned over all intervals, and filters the indices by some threshold displacement

Parameters

coordinates : list of list of pandas.DataFrame: A pandas dataframe containing coordinates for each timestep. Must be indexed by particle (with each particle having a unique identifyer that matches between different time steps) and contain coordinates along each dimension in a separate column, with column names matching those given in pos_cols.
times : list, optional: timestamps for the sets of coordinates
pos_cols : tuple of str, optional: names of the columns of the DataFrames in coordinates containing the particle coordinates along each dimension. The length (i.e. number of dimensions) must match len(boundary) and the number of columns in coordinates. The default is ('z','y','x').
msd_min : float, optional: Lower limit for the mean squared displacement between two time steps below which a particle is considered stationary. The default is 0.01.
msd_max : float, optional: Upper limit on the MSD between two timesteps. The default is 1.
interval : int, optional: time interval in integer index units used for the calculations, e.g. averaging over neigbouring steps or over larger intervals than the sampling rate in the data. The default is 1.

Returns

indices : set: particle identifiers for the particles which are NOT considered stationary, i.e. those whose MSD is LARGER than the threshold.
msds : list of float: full list of the MSD values for all particles in all timesteps, useful e.g. for plotting a histogram to estimate a value for threshold.

def load_forceprofile(filename)

load results file as stored by save_forceprofile()

Parameters

filename : str: filename with extension to store data in

Returns

m : np.ndarray of int: bin indices
r_cent : np.ndarray of float: center r position of each bin.
G : np.ndarray of float: force coefficients.
G_err : np.ndarray of float: residual error in G.
binmean_r : np.ndarray of float: count weighted mean distance for all pairs in each bin.
counts : np.ndarray of int or float: number of pair counts in each bin weighted for the basis functions.
M : int: number of discretisation steps
rmax : float: cut-off distance for the force.
basis_function : str: type of basis functions used.
gamma : float: friction factor used for the calculations.
periodic_boundary : bool: whether periodic boundary conditions were used.

def load_forceprofile_cylindrical(filename)

load results file as stored by save_forceprofile_cylindrical()

Parameters

filename : str: filename with extension to store data in

Returns

m_z : np.ndarray of int: z bin indices
m_rho : np.ndarray of int: rho bin indices.
z_cent : np.ndarray of float: z position of center of each bin.
rho_cent : np.ndarray of float: rho position of center of each bin.
G_z : np.ndarray of float: force coefficients along z.
G_rho : np.ndarray of float: force coefficients along rho.
G_z_err : np.ndarray of float: residual error in G_z.
G_rho_err : np.ndarray of float: residual error in G_rho.
binmean_z : np.ndarray of float: count weighted mean z distance for all pairs in each bin.
binmean_rho : np.ndarray of float: count weighted mean rho distance for all pairs in each bin.
counts : np.ndarray of int or float: number of pair counts in each bin weighted for the basis functions.
M_z : int: number of discretisation steps along the z axis
M_rho : int: number of discretisation steps along the rho axis.
rmax : float: cut-off distance for the force.
basis_function : str: type of basis functions used.
gamma : float: friction factor used for the calculations.
periodic_boundary : bool: whether periodic boundary conditions were used.

def run_inertial(coordinates, times, boundary=None, mass=1, rmax=1, M=20, pos_cols=['z', 'y', 'x'], periodic_boundary=False, bruteforce=False, remove_near_boundary=False, solve_per_dim=False, return_data=False, neighbour_upper_bound=None)

Run the analysis for inertial dynamics (molecular dynamics like), iterates over all subsequent sets of three timesteps and obtains forces from the accellerations of the particles. Based on [1].

Parameters

coordinates : list of pandas.DataFrame: A pandas dataframe containing coordinates for each timestep. Must be indexed by particle (with each particle having a unique identifyer that matches between different time steps) and contain coordinates along each dimension in a separate column, with column names matching those given in pos_cols.
times : list of float: list timestamps corresponding to the coordinates with length matching len(coordinates).
boundary : list or tuple, optional: boundaries of the box in which the coordinates are defined in the form ((d0_min,d0_max),(d1_min,d1_max),…) with a length (number) and order of dimensions matching pos_cols. The default is None, which uses the min and max value found in the entire set of coordinates along each axis.
mass : float, optional: particle mass for calculation of F=m*a. The default is 1.
pos_cols : list of str, optional: names of the columns of the DataFrames in coordinates containing the particle coordinates along each dimension. The length (i.e. number of dimensions) must match len(boundary) and the number of columns in coordinates. The default is `['z','y','x'].
rmax : float, optional: cut-off radius for calculation of the pairwise forces. The default is 1.
M : int, optional: The number of discretization steps for the force profile, i.e. the number of bins from 0 to rmax into which the data will be sorted. The default is 20.
periodic_boundary : bool, optional: Whether the box has periodic boundary conditions. If True, the boundary must be given. The default is False.
bruteforce : bool, optional: If True, the coefficients are calculated in a naive brute-force approach with a nested loop over all particles. The default is False, which uses a scipy.spatial.cKDTree based approach to only evaluate pairs which are <rmax apart.
remove_near_boundary : bool, optional: If true, particles which are closer than rmax from any of the boundaries are not analyzed, but still accounted for when analyzing particles closer to the center of the box in order to only analyze particles for which the full spherical shell up to rmax is within the box of coordinates, and to prevent erroneous handling of particles which interact with other particles outside the measurement volume. Only possible when periodic_boundary=False. The default is True.
solve_per_dim : bool, optional: if True, the matrix is solved for each dimension separately, and a force vector and error are returned for each dimension.
return_data : bool, optional: If True, the full coefficient matrix and force vector are returned together with the force vector, error and the list of bincounts. The default is False.
neighbour_upper_bound : int: upper bound on the number of neighbours within rmax a particle may have to limit memory use and computing time in the pair finding step. The default is the total number of particles in each time step.

Returns

G : numpy.array of length M: discretized force vector, the result of the computation.
G_err : numpy.array of length M: errors in G based on the least_squares solution of the matrix equation
counts : numpy.array of length M: number of individual force evaluations contributing to the result in each bin.
coefficients : numpy.array of M by 3n*(len(times)-1): coefficient matrix of the full dataset as specified in [1]. This is only returned when return_data=True
forces : numpy.array of length 3n*(len(times)-1): vector of particle forces of form [t0p0z,t0p0y,t0p0x,t0p1z,t0p1y, …, tn-1pnx]. This is only returned when return_data=True

References

[1] Jenkins, I. C., Crocker, J. C., & Sinno, T. (2015). Interaction potentials from arbitrary multi-particle trajectory data. Soft Matter, 11(35), 6948–6956. https://doi.org/10.1039/C5SM01233C

def run_overdamped(coordinates, times, boundary=None, gamma=1, rmax=1, M=20, pos_cols=('z', 'y', 'x'), eval_particles=None, periodic_boundary=False, basis_function='constant', bruteforce=False, remove_near_boundary=True, constant_particles=False, solve_per_dim=False, neighbour_upper_bound=None, newline=False, use_gpu=False)

Run the analysis for overdamped dynamics (brownian dynamics like), iterates over all subsequent sets of two timesteps and obtains forces from the velocity of the particles as a function of the distribution of the particles around eachother. Based on ref. [1].

Parameters

coordinates : (list of) list of pandas.DataFrame: A pandas dataframe containing coordinates for each timestep as a series of consecutive timesteps of at least 2 items (i.e. 1 time interval). Multiple nonconsecutive series may be given as list of lists of DataFrames. DataFrames must be indexed by particle (with each particle having a unique identifyer that matches between different time steps) and contain coordinates along each dimension in a separate column, with column names matching those given in pos_cols.
times : (list of) list of float: list(s) of timestamps corresponding to the coordinate sets
boundary : tuple, list of tuple or list of list of tuple, optional: boundaries of the box in which the coordinates are defined in the form ((d0_min,d0_max),(d1_min,d1_max),…) with a length (number) and order of dimensions matching pos_cols. A single set of boundaries may be given for all timesteps, or a (list of) list of boundaries for each timestep may be specified. The default is None, which uses the min and max value found in the entire set of coordinates along each axis.
gamma : float, optional: damping/friction coefficient (kT/D) for calculation of F=V*kT/D. The default is 1.
rmax : float, optional: cut-off radius for calculation of the pairwise forces. The default is 1.
M : int, optional: The number of discretization steps for the force profile, i.e. the number of bins from 0 to rmax into which the data will be sorted. The default is 20.
pos_cols : tuple of str, optional: names of the columns of the DataFrames in coordinates containing the particle coordinates along each dimension. The length (i.e. number of dimensions) must match len(boundary) and the number of columns in coordinates. The default is ('z','y','x').
eval_particles : set: set of particle id's (matching the indices in coordinates) to use in the force evaluation, such that forces are not calculated for any particle not in the set of eval_particles. Note that all particles are always used to calculate the coefficients.
periodic_boundary : bool, optional: Whether the box has periodic boundary conditions. If True, the boundary must be given. The default is False.
basis_function : ['constant', 'linear']: the type of basis functions to use, where 'constant' uses square wave basis functions which assume the force is constant over each bin, and 'linear' uses linear wave basis functions where each pair contributes to the nearby bins with linearly interpolated weights. The default is 'constant'.
bruteforce : bool, optional: If True, the coefficients are calculated in a naive brute-force approach with a nested loop over all particles. The default is False, which uses a scipy.spatial.cKDTree based approach to only evaluate pairs which are <rmax apart.
remove_near_boundary : bool, optional: If true, particles which are closer than rmax from any of the boundaries are not analyzed, but still accounted for when analyzing particles closer to the center of the box in order to only analyze particles for which the full spherical shell up to rmax is within the box of coordinates, and to prevent erroneous handling of particles which interact with other particles outside the measurement volume. Only possible when periodic_boundary=False. The default is True.
constant_particles : bool, optional: when the same set of particles is present in each timestep, i.e. the indices of coordinates are identical for all time steps after selecting eval_particles, more efficient (indexing) algorithms can be used
solve_per_dim : bool, optional: if True, the matrix is solved for each dimension separately, and a force vector and error are returned for each dimension.
neighbour_upper_bound : int, optional: upper bound on the number of neighbours within rmax a particle may have to limit memory use and computing time in the pair finding step. The default is the total number of particles in each time step.
newline : bool, optional: whether to print output for each series on a new line. The default is True.
use_gpu : bool, optional: if True, matrix operations for the least squares solver are offloaded to the gpu via CuPy, which requires a cuda-compatible gpu. For large numbers of bins (>>100) this may result in significantly better performance but there is considerable overhead in moving data between cpu and gpu, so for most use cases this is not faster.

Returns

G : numpy.array of length M: discretized force vector, the result of the computation.
G_err : numpy.array of length M: errors in G based on the least_squares solution of the matrix equation
counts : numpy.array of length M: number of individual force evaluations contributing to the result in each bin.
mean_rho : numpy.array of length M: the average interparticle distance of all particle pairs contributing to each bin, i.e. the mean of the distances for which the forces were evaluated.

References

[1] Jenkins, I. C., Crocker, J. C., & Sinno, T. (2015). Interaction potentials from arbitrary multi-particle trajectory data. Soft Matter, 11 (35), 6948–6956. https://doi.org/10.1039/C5SM01233C

def run_overdamped_cylindrical(coordinates, times, boundary=None, gamma=1, rmax=1, M_z=20, M_rho=20, pos_cols=('z', 'y', 'x'), eval_particles=None, periodic_boundary=False, basis_function='constant', remove_near_boundary=True, constant_particles=False, neighbour_upper_bound=None, use_gpu=False, check_coordinates=True, newline=True)

Run the analysis for overdamped dynamics (brownian dynamics like) in a cylindrical coordinate system where interactions are radially averaged in the xy plane (dimensions 2 and 1) but not along the cylinder (z) axis.

Parameters

coordinates : (list of) list of pandas.DataFrame: A pandas dataframe containing coordinates for each timestep as a series of consecutive timesteps of at least 2 items (i.e. 1 time interval). Multiple nonconsecutive series may be given as list of lists of DataFrames. DataFrames must be indexed by particle (with each particle having a unique identifyer that matches between different time steps) and contain coordinates along each dimension in a separate column, with column names matching those given in pos_cols.
times : (list of) list of float: list(s) of timestamps corresponding to the coordinate sets
boundary : tuple, list of tuple or list of list of tuple, optional: boundaries of the box in which the coordinates are defined in the form ((d0_min,d0_max),(d1_min,d1_max),…) with a length (number) and order of dimensions matching pos_cols. A single set of boundaries may be given for all timesteps, or a (list of) list of boundaries for each timestep may be specified. The default is None, which uses the min and max value found in the entire set of coordinates along each axis.
gamma : float, optional: damping/friction coefficient (kT/D) for calculation of F=V*kT/D. The default is 1.
rmax : float, optional: cut-off radius for calculation of the pairwise forces. The default is 1.
M_z : int, optional: The number of discretization steps for the force profile along the z (cylinder) direction, i.e. the number of bins from 0 to rmax into which the data will be sorted. The default is 20.
M_rho : int, optional: The number of discretization steps along the axial (rho) direction. The default is 20.
pos_cols : tuple of str, optional: names of the columns of the DataFrames in coordinates containing the particle coordinates along each dimension. The first is assumed to be the axial (z) direction and the second and third define the plane for the rho component. The default is ('z','y','x').
eval_particles : set, optional: set of particle id's (matching the indices in coordinates) to use in the force evaluation, such that forces are not calculated for any particle not in the set of eval_particles. Note that all particles are always used to calculate the coefficients. The default is to evaluate all particles in coordinates.
periodic_boundary : bool, optional: Whether the box has periodic boundary conditions. If True, the boundary must be given. The default is False.
basis_function : ['constant', 'linear']: the type of basis functions to use, where 'constant' uses square wave basis functions which assume the force is constant over each bin, and 'linear' uses linear wave basis functions where each pair contributes to the nearby bins with linearly interpolated weights. The default is 'constant'.
remove_near_boundary : bool, optional: If true, particles which are closer than rmax from any of the boundaries are not analyzed, but still accounted for when analyzing particles closer to the center of the box in order to only analyze particles for which the full spherical shell up to rmax is within the box of coordinates, and to prevent erroneous handling of particles which interact with other particles outside the measurement volume. Only possible when periodic_boundary=False. The default is True.
constant_particles : bool, optional: when the same set of particles is present in each timestep, i.e. the indices of coordinates are identical for all time steps after selecting eval_particles, more efficient (indexing) algorithms can be used
neighbour_upper_bound : int: upper bound on the number of neighbours within rmax a particle may have to limit memory use and computing time in the pair finding step. The default is the total number of particles in each time step.
use_gpu : bool, optional: if True, matrix operations for the least squares solver are offloaded to the gpu via CuPy, which requires a cuda-compatible gpu. For large numbers of bins (M_z×M_rho>>100) this may result in significantly better performance. There is considerable overhead in moving data between cpu and gpu, so for small numbers of bins this is not faster.
check_coordinates : bool, optional: by default, coordinates are checked for being within boundaries. This can be bypassed by setting this to False.
newline : bool, optional: whether to print output for each series on a new line. The default is True.

Returns

m_z : numpy.array of ints of shape (M_z,M_rho): z axis bin indices corresponding to the data arrays.
m_rho : numpy.array of ints of shape (M_z,M_rho): rho axis bin indices corresponding to the data arrays.
G_z : numpy.array of floats of shape (M_z,M_rho): values for the axial (z) component of the force in each bin, i.e. the magnitude associated with each basis function.
G_rho : numpy.array of floats of shape (M_z,M_rho): values for the in-plane (rho) component of the force in each bin, i.e. the magnitude associated with each basis function.
G_z_err : numpy.array of floats of shape (M_z,M_rho): residual least-squares error in each bin in G_z
G_rho_err : numpy.array of floats of shape (M_z,M_rho): residual least-squares error in each bin in G_rho
mean_z : numpy.array of floats of shape (M_z,M_rho): the average absolute z distance of all particle pairs contributing to each bin, i.e. the mean of the distances for which the forces were evaluated.
mean_rho : numpy.array of floats of shape (M_z,M_rho): the average absolute planar (rho) distance of all particle pairs contributing to each bin, i.e. the mean of the distances for which the forces were evaluated.
counts : numpy.array of floats of shape (M_z,M_rho): number of individual force evaluations contributing to the result in each bin. When basis_function='linear' this is a total weight rather than an integer count, since each particle pair has its unity weight divided over the 4 surrounding bins.

def save_forceprofile(filename, rmax, M, basis_function, gamma, periodic_boundary, binmean_r, G, G_err, counts)

Saves the results from trajectory analysis to a text file

Parameters

filename : string: filename to use for results file
binmean_r : list of float: mean value of all pairs in each bin, i.e. a 'count weighted' bin center
forces : list of float: list of force values as obtained from the trajectory analysis
counts : list of int: the number of evaluations used for each bin
M : int: number of bins
rmax : float: cut-off radius for force
gamma : float: the friction coefficient used in the calculation

def save_forceprofile_cylindrical(filename, rmax, M_z, M_rho, basis_function, gamma, periodic_boundary, m_z, m_rho, G_z, G_rho, G_z_err, G_rho_err, binmean_z, binmean_rho, counts)

save results from cylindrical trajectory analysis to a file

Parameters

filename : str: filename with extension to store data in.
rmax : float: cut-off distance for the force.
M_z : int: number of discretisation steps along the z axis
M_rho : int: number of discretisation steps along the rho axis.
basis_function : str: type of basis functions used.
gamma : float: friction factor used for the calculations.
periodic_boundary : bool: whether periodic boundary conditions were used.
m_z : np.ndarray of int: z bin indices
m_rho : np.ndarray of int: rho bin indices.
G_z : np.ndarray of float: force coefficients along z.
G_rho : np.ndarray of float: force coefficients along rho.
G_z_err : np.ndarray of float: residual error in G_z.
G_rho_err : np.ndarray of float: residual error in G_rho.
binmean_z : np.ndarray of float: count weighted mean z distance for all pairs in each bin.
binmean_rho : np.ndarray of float: count weighted mean rho distance for all pairs in each bin.
counts : np.ndarray of int or float: number of pair counts in each bin weighted for the basis functions.