softspot

The Python library softspot analyzes low-energy quasilocalized excitations in computer glasses, based solely on the harmonic approximation of the potential energy. It accompanies the research paper A simple and broadly-applicable definition of shear transformation zones by David Richard, Geert Kapteijns, Julia Giannini, Lisa Manning, and Edan Lerner.

See the GitLab repository for more information, and the iPython tutorial to get started with a simple example.

class softspot.softspot.SoftSpot(ndim, npart, hessian, pair_lookup=None)
Parameters:
  • ndim (int) – The dimensionality of the system under consideration.
  • npart (int) – The number of particles.
  • hessian (scipy.sparse.csr_matrix) – Hessian matrix in scipy’s Compressed Sparse Row representation.
  • pair_lookup (array_like) – Numpy array of shape (num_contacts, 2) with dtype=int containing particle pairs for calculating sum_zij_squared (and its gradient) — the denominator of Eq. 2 of the paper. For example, pair_lookup[0] == [4, 5] indicates that particles 4 and 5 will contribute to sum_zij_squared. Default is None: in most use cases — i.e., short-range potentials — the relevant particle pairs can be automatically derived from the Hessian matrix, and will coincide with the interacting pairs. In this case, the user does not have to specify this field. In the case of long-range interactions, the user might want to manually specify only nearby pairs to calculate sum_zij_squared, in order to effectively penalize modulating (phonon) fields in the cost function (see discussion below Eq. 3 in the paper).

Example

import numpy as np

# Example for a high-parent-temperature 2D SWAP system
# consisting of 4096 particles.
# See iPython tutorial for more details.
np.random.seed(1)
ndim = 2
npart = 4096
z_init = np.random.rand((ndim*npart))
soft_spot = SoftSpot(ndim=ndim, npart=npart, hessian=hessian)
result = soft_spot.find(z_init, mode='cg')
# => {
#    'pi': array([-0.00445838,  0.00169628, -0.00328873, ...,  0.00055183, -0.00468846,  0.00158765]),
#    'kappa': 3.7748011979346128,
#    'cost_function': 36.194330960051445
#    'converged': True,
#    'tol': 9.87820288575708e-11,
#    'iterations': 1085,
# }
find(z_init, mode='cg', options={})

Finds a pseudo-harmonic mode of the system. This is the main interface for regular users.

Parameters:
  • z_init (array_like) – Initial direction on the energy landscape from which to start the optimization. Must be a 1-D numpy array of size npart * ndim. Before the optimization starts, z_init will be made translation-free (since by assumption the potential energy is translationally invariant) and will be normalized.
  • mode (str) – Either 'cg' (default) or 'mapping'. If 'cg', will minimize the cost function with a conjugate gradient minimizer. If 'mapping', will repeatedly apply the mapping (Eq. 4 of the paper) until convergence.
  • options (dict, optional) –

    Dictionary of solver options.

    When mode is 'cg', the most pertinent options are:

    tol : float
    Tolerance criterion for the dimensionless norm of the gradient. Default: 1e-10.
    max_iterations : int
    Maximum number of iterations. Default: 100000.

    In general, all possible options of the Macopt conjugate gradient minimizer may be supplied, except for the normalize_function, which is overridden.

    When mode is 'mapping', the possible options are

    tol : float
    Tolerance criterion for the convergence of the mapping. Default: 1e-10.
    tol_cg: float
    Tolerance criterion for the conjugate gradient iterative linear solver of the equation H x = zeta (see Eqs. 4 and 5 of the paper). Default: 1e-10.
    max_iterations : int
    Maximum times to iteratively apply the mapping. Default: 20000.
    max_iterations_cg : int
    Maximum iterations for the conjugate gradient iterative linear solver of the equation H x = zeta (see Eqs. 4 and 5 of the paper). Default: 10000.
Returns:

result –

Dictionary with the result of the optimization.

pi : array_like

Pseudo-harmonic mode resulting from the optimization.

kappa : float

Stiffness of pi.

cost_function : float

Value of the cost function at the minimum pi.

iterations : int

Number of iterations.

tol : float

Final tolerance.

converged : bool

False if the optimization did not converge, True otherwise.
Return type:

dict
cost_function(z)

Returns the value of the cost function C(z) = (H:zz)^2 / sum_zij_squared (see Eq. 2 of the paper). This function is not called by the user in the regular use case of this class, but may be used by to those who wish to implement their own optimization schemes.

Parameters:z (array_like) – Direction on the energy landscape. Must be a 1-D numpy array of size npart * ndim.
gradient_cost_function(z)

Returns the gradient of the cost function C(z) = (H:zz)^2 / sum_zij_squared (Eq. 2 of the paper). This function is not called by the user in the regular use case of this class, but may be used by to those who wish to implement their own optimization schemes.

Parameters:z (array_like) – Direction on the energy landscape. Must be a 1-D numpy array of size npart * ndim.
Returns:
  • gradient (array_like) – Gradient of the cost function.
  • convergence (float) – Convergence criterion; this is the dimensionless norm of the gradient:

    np.sqrt( np.dot(gradient, gradient)  / np.dot(hesOnZ, hesOnZ)) / kappa, where hesOnZ is the Hessian applied to z, and kappa is the stiffness of z.
mapping(z, z0_cg, options)

Returns the result of applying the mapping once to the vector z (see Eq. 4 of the paper). This function is not called by the user in the regular use case of this class, but may be used by to those who wish to implement their own optimization schemes.

Parameters:
  • z (array_like) – Direction on the energy landscape. Must be a 1-D numpy array of size npart * ndim.
  • z0_cg (array_like) – Initial vector for the conjugate gradient iterative linear solver of the equation H x = zeta (See Eqs. 4 and 5 of the paper).
  • options (dict) – Options for the conjugate gradient iterative linear solver of the equation H x = zeta (See Eqs. 4 and 5 of the paper). Must contain the options tol_cg and max_iterations_cg described in the find method.
Returns:

  • mapped_z (array_like) – Result of applying the mapping to z.
  • convergence (float) – Convergence criterion np.linalg.norm(mapped_z - z) / npart.