DirectMaxFlux

class dmf.dmf.DirectMaxFlux(ref_images, coefs=None, nsegs=4, dspl=3, remove_rotation_and_translation=True, mass_weighted=False, parallel=False, world=None, t_eval=None, w_eval=None, n_vel=None, n_trans=None, n_rot=None, eps_vel=0.01, eps_rot=0.01, beta=10.0, nmove=5, update_teval=False, params_t_update={'de': 0.15, 'dia': 1.0, 'dib': 0.2, 'epsb': 0.02, 'max_alpha0': 0.1, 'mua': 5.0, 'mub': 5.0})

Bases: VariationalPathOpt

Variational reaction path/transition states optimization based on the direct MaxFlux method.

Ref. 1.

S.-i. Koda and S. Saito, Locating Transition States by Variational Reaction Path Optimization with an Energy-Derivative-Free Objective Function J. Chem. Theory Comput. 20, 2798–2811 (2024).

This class implements the MaxFlux variational principle in the large-β (low-temperature) regime for locating transition states (TSs) and approximating minimum-energy paths (MEPs), following the formulation of Ref. 1.

The reaction path ( x(t) ) is represented by a B-spline expansion, and the following functional is minimized:

\[\tilde{I}[x] = \beta^{-1} \log I[x],\]

where

\[I[x] = \int_0^1 dt\, \vert \dot{x}(t) \vert \, e^{\beta E(x(t))}.\]

With large ( beta ), the highest-energy point along the optimized path approximates the TS geometry.

The method requires only first-order atomic forces, because the objective contains no derivatives of the potential energy.

Additional features include: - construction of initial B-spline coefficients from ref_images - optional removal of translational and rotational redundancy - parallel energy/force evaluation (threads or MPI) - optional adaptive refinement of t_eval near the high-energy region

Parameters:

• ref_images (list of ase.Atoms) – List of atomic structures representing an initial guess for the path. If coefs is not provided, a piecewise linear interpolation through ref_images is constructed, and B-spline coefficients are obtained by fitting this interpolated path. If coefs is provided, no interpolation is performed: ref_images[0] is used only to extract atomic numbers, masses, cell, and PBC settings.

• coefs (ndarray of shape (nbasis, natoms, 3), optional) – Initial B-spline coefficients. If provided, interpolation from ref_images is skipped and these coefficients define the initial path. Default: None.

• nsegs (int, optional) – Number of B-spline segments. The number of basis functions per Cartesian degree of freedom is nbasis = nsegs + dspl. See Ref. 1 for details. Default: 4.

• dspl (int, optional) – Polynomial degree of the B-spline basis. Default: 3.

• remove_rotation_and_translation (bool, optional) – If True, remove global translational and rotational motion using nonlinear constraints. Default: True.

• mass_weighted (bool, optional) – If True, the velocity norm ( vert dot{x}(t) vert ) uses mass-weighted coordinates. Default: False.

• parallel (bool, optional) – Evaluate energies and forces in parallel using threads or MPI. Default: False.

• world (MPI communicator, optional) – Communicator used when parallel=True. Defaults to ase.parallel.world.

• t_eval (ndarray of shape (nmove+2,), optional) – Initial evaluation points in ( t in [0,1] ). If omitted or update_teval is True, an even distribution np.linspace(0.0,1.0,nmove+2) is generated.

• w_eval (ndarray of shape (nmove+2,), optional) –

  Initial quadrature weights for evaluating the integral

  \[I[x] \approx \sum_i w_i\, \vert \dot{x}(t) \vert \, e^{\beta E(x(t_i))}.\]

  If omitted, trapezoidal weights are used.

• n_vel (int, optional) – Number of discretized velocity constraints. Default: 4 * nsegs.

• n_trans (int, optional) – Number of translational constraints. Default: 2 * nsegs.

• n_rot (int, optional) – Number of rotational constraints. Default: 2 * nsegs.

• eps_vel (float, optional) – Tolerance for velocity constraints. Default: 0.01.

• eps_rot (float, optional) – Tolerance for rotational constraints. Default: 0.01.

• beta (float, optional) – Reciprocal temperature ( beta ) (in 1/eV) used in the MaxFlux functional. Default: 10.0.

• nmove (int, optional) – Number of movable interior evaluation points. Total number of images = nmove + 2 (including both endpoints). Default: 5.

• update_teval (bool, optional) – If True, t_eval is adaptively updated toward the high-energy region during optimization. Default: False.

• params_t_update (dict, optional) – Parameters controlling the update of t_eval. Includes keys such as max_alpha0, de, dia, mua, dib, mub, epsb (Defaults are the same as in Ref. 1).

# ---- Path representation ----

images

Atomic structures at the current t_eval. The length of this list is nmove + 2 (including both endpoints).

Type:

list of ase.Atoms

coefs

Current B-spline coefficients defining the variational path.

Type:

ndarray of shape (nbasis, natoms, 3)

angs

Euler angles used when removing rotational redundancy.

Type:

ndarray of shape (3,)

# ---- Energies and forces ----

energies

Energies evaluated at the current t_eval (shifted by e0).

Type:

ndarray of shape (nmove+2,)

forces

Forces evaluated at the current t_eval.

Type:

ndarray of shape (nmove+2, natoms, 3)

# ---- Evaluation grid ----

t_eval

Current energy evaluation points along the path. If update_teval=True, these differ from the initial values.

Type:

ndarray of shape (nmove+2,)

w_eval

Current quadrature weights associated with t_eval.

Type:

ndarray of shape (nmove+2,)

# ---- Constraint configuration ----

n_vel

Number of velocity constraints.

Type:

int

n_trans

Number of translational constraints.

Type:

int

n_rot

Number of rotational constraints.

Type:

int

eps_vel

Tolerance for velocity constraints.

Type:

float

eps_rot

Tolerance for rotational constraints.

Type:

float

remove_rotation_and_translation

Whether translational/rotational redundancy is removed.

Type:

bool

# ---- B-spline representation ----

nsegs

Number of B-spline segments.

Type:

int

dspl

Degree of the B-spline basis.

Type:

int

nbasis

Number of B-spline basis functions per Cartesian degree of freedom. nbasis = nsegs + dspl.

Type:

int

# ---- MaxFlux functional parameters ----

beta

Reciprocal temperature ( beta ).

Type:

float

nmove

Number of movable interior evaluation points.

Type:

int

update_teval

Whether t_eval is adaptively updated.

Type:

bool

params_t_update

Parameters controlling the update of t_eval.

Type:

dict

# ---- Optimization ----

ipopt_options

IPOPT options used for the optimization.

Type:

dict

history

Container storing iteration-by-iteration quantities.

Type:

HistoryDMF

add_ipopt_options(dict_options)

Add or update IPOPT options.

This method updates self.ipopt_options with the key–value pairs given in dict_options and forwards them to IPOPT via self.add_option.

Parameters:

dict_options (dict) – Dictionary of IPOPT options (e.g., {"tol": 1e-3}).

add_option(keyword, val)

Add a keyword/value option pair to the problem.

See the Ipopt documentaion for details on available options.

Parameters:

• keyword (str) – Option name.

• val (str, int or float) – Value of the option. The type of val should match the option definition as described in the Ipopt documentation.

close()

Deallocate memory resources used by the Ipopt package.

Called implicitly by the Problem class destructor.

constraints(x)

IPOPT callback: nonlinear constraint values.

This method is not intended to be called directly by users. It is invoked internally by IPOPT during the optimization process. The argument x is the flattened optimization variable vector, and the return value is the array of constraint values at that state.

property e0

float: Minimum endpoint energy used to shift the energy scale.

get_current_iterate(scaled=False)

Return the current iterate vectors during an Ipopt solve

The iterate contains vectors for primal variables, bound multipliers, constraint function values, and constraint multipliers. Here, the constraints are treated as a single function rather than separating equality and inequality constraints. This method can only be called during an intermediate callback.

Only supports Ipopt >=3.14.0

Parameters:

scaled (bool) – Whether the scaled iterate vectors should be returned

Returns:

:py:class:`dict` or `None` – A dict containing the iterate vector with keys "x", "mult_x_L", "mult_x_U", "g", and "mult_g". If iterate vectors cannot be obtained, None is returned.

get_current_violations(scaled=False)

Return the current violation vectors during an Ipopt solve

Violations returned are primal variable bound violations, bound complementarities, the gradient of the Lagrangian, constraint violation, and constraint complementarity. Here, the constraints are treated as a single function rather than separating equality and inequality constraints. This method can only be called during an intermediate callback.

Only supports Ipopt >=3.14.0

Parameters:

scaled (bool) – Whether to scale the returned violations

Returns:

:py:class:`dict` or `None` – A dict containing the violation vector with keys "x_L_violation", "x_U_violation", "compl_x_L", "compl_x_U", "grad_lag_x", "g_violation", and "compl_g". If violation vectors cannot be obtained, None is returned.

get_positions(t=None, P=None, nu=0)

Evaluate the positions (or their derivatives) along the path.

Normally, users provide only t; however, advanced users may supply precomputed basis values P (from _get_basis_values()) to avoid repeated evaluations.

If both t and P are provided, P takes priority.

Parameters:

• t (ndarray, optional) – 1D array of parameter values in [0, 1] at which positions (or derivatives) are evaluated. If omitted, t_eval is used.

• P (ndarray, optional) – Precomputed B-spline basis values from _get_basis_values(). Default: None.

• nu (int, optional) – Derivative order with respect to t (0, 1, or 2). Default: 0.

Returns:

ndarray – Array of shape (len(t), natoms, 3) containing the positions (nu = 0) or the nu-th derivatives of the path.

get_x()

Return the flattened optimization variable vector used by IPOPT.

Although mainly intended for internal use, this method is exposed because solve() returns the optimized variable vector x. The returned array contains all internal degrees of freedom:

• the flattened interior B-spline coefficients coefs[1:-1] (endpoint coefficients are fixed), and

• the rotation angles angs if remove_rotation_and_translation=True.

Returns:

x (ndarray of shape (nvar,)) – Flattened optimization variable vector.

gradient(x)

IPOPT callback: gradient of the objective.

This method is not intended to be called directly by users. It is invoked internally by IPOPT during the optimization process. The argument x is the flattened optimization variable vector, and the return value is the gradient of the objective with respect to x.

interpolate_energies(t_eval=None, energies=None, forces=None, coefs=None, delta_e=None)

Construct a piecewise-cubic interpolation of the energy along the path.

This method reconstructs a smooth interpolation \(\tilde{E}(t)\) of the discrete energy values evaluated at t_eval. The interpolation is C^1-continuous and uses both energies and their first derivatives.

Optionally, the method can also locate the values of t satisfying

\[\tilde{E}(t) = E_{\max} - \Delta E,\]

for user-specified delta_e.

See Ref. 1 for details.

Parameters:

• t_eval (ndarray, optional) – 1D array of parameter values at which energies/forces were evaluated. If omitted, self.t_eval is used. Only the region t_eval <= 1 is used internally.

• energies (ndarray, optional) – Energy values at t_eval. If omitted, self.energies is used.

• forces (ndarray, optional) – Forces at t_eval with shape (len(t_eval), natoms, 3). If omitted, self.forces is used.

• coefs (ndarray of shape (nbasis, natoms, 3), optional) – B-spline control-point coefficients. If omitted, self.coefs is used.

• delta_e (list of float, optional) –

  Energy offsets \(\Delta E\). If provided, this method also returns the corresponding parameter values t satisfying

  \[\tilde{E}(t) = E_{\max} - \Delta E.\]

Returns:

• polys (ndarray of shape (len(t_eval) - 1, 4)) – Polynomial coefficients defining the piecewise cubic interpolation. Each segment corresponds to:

  \[\tilde{E}(t) = c_0 + c_1 t + c_2 t^2 + c_3 t^3.\]

• t_max (float) – The parameter value t at which the interpolated energy \(\tilde{E}(t)\) attains its maximum.

• e_max (float) – The maximum interpolated energy \(\tilde{E}(t_{\max})\).

• t_de (list of ndarray, optional) – Returned only when delta_e is provided. t_de[j] contains all roots satisfying \(\tilde{E}(t) = E_{\max} - \Delta E_j\).

jacobian(x)

IPOPT callback: Jacobian of the constraints.

This method is not intended to be called directly by users. It is invoked internally by IPOPT during the optimization process. The argument x is the flattened optimization variable vector, and the return value is the Jacobian matrix of the constraint functions.

objective(x)

IPOPT callback: objective function.

This method is not intended to be called directly by users. It is invoked internally by IPOPT during the optimization process. The argument x is the flattened optimization variable vector, and the return value is the scalar objective evaluated at that state.

set_coefs_angs(coefs=None, angs=None)

Update the B-spline coefficients and/or rotation angles.

This method updates coefs and angs if the corresponding arguments are provided. After updating the angles, the final B-spline control point (coefs[-1]) is recomputed as

\[\mathrm{coefs}[-1] = \mathrm{coefs}_0[-1] \, R_x R_y R_z,\]

where R_x, R_y, R_z are the rotation matrices generated from self.angs. This ensures that the endpoint geometry is kept consistent under rotational constraints.

Parameters:

• coefs (ndarray of shape (nbasis, natoms, 3), optional) – New B-spline coefficients. If omitted, the current coefficients are preserved.

• angs (ndarray of shape (3,), optional) – Rotation angles used to the final endpoint alignment. If omitted, the current angles are preserved.

set_positions(coefs=None, angs=None)

Update the positions of all images along the path.

This method first updates the B-spline coefficients and/or rotation angles by calling set_coefs_angs(). It then recomputes the atomic positions along the path using get_positions(), and writes these positions into the existing self.images objects.

Note that this method does not change the number of images; it only updates their positions according to the current path parameters.

Parameters:

• coefs (ndarray of shape (nbasis, natoms, 3), optional) – New B-spline coefficients. If omitted, the existing coefficients are preserved.

• angs (ndarray of shape (3,), optional) – Rotation angles used for endpoint alignment. If omitted, the existing angles are preserved.

set_problem_scaling(obj_scaling=1.0, x_scaling=None, g_scaling=None)

Optional function for setting scaling parameters for the problem.

To use the scaling parameters set the option nlp_scaling_method to user-scaling.

Parameters:

• obj_scaling (float) – Determines, how Ipopt should internally scale the objective function. For example, if this number is chosen to be 10, then Ipopt solves internally an optimization problem that has 10 times the value of the original objective. In particular, if this value is negative, then Ipopt will maximize the objective function instead of minimizing it.

• x_scaling (array-like, shape(n, )) – The scaling factors for the variables. If None, no scaling is done.

• g_scaling (array-like, shape(m, )) – The scaling factors for the constrains. If None, no scaling is done.

set_t_eval(t_eval)

Set the energy evaluation points t_eval.

This also updates the cached B-spline basis values used for evaluating positions and derivatives.

Parameters:

t_eval (ndarray) – 1D array of parameter values in the interval [0, 1]. Its length must match the length of the initial t_eval used at initialization, because the number of images is fixed.

set_w_eval(w_eval=None)

Set the quadrature weights w_eval used in the action integral.

If w_eval is not provided, trapezoidal weights are generated from the current values of t_eval. The number of weights must match the number of energy evaluation points, which is fixed after initialization.

Parameters:

w_eval (ndarray, optional) – 1D array of quadrature weights corresponding to t_eval. Its length must match that of t_eval. If omitted, trapezoidal-rule weights are constructed automatically.

set_x(x)

Update coefs and angs from the flattened optimization vector.

This method is normally not called directly by end-users; it is invoked internally during IPOPT callbacks such as objective(), gradient(), constraints(), and jacobian().

The input vector x must be exactly the one produced by get_x(). The method reconstructs:

• coefs[1:-1] (interior B-spline control points), and

• angs (rotation angles, if enabled),

and then updates all image positions via set_positions().

Parameters:

x (ndarray of shape (nvar,)) – Flattened optimization variable vector.

solve(tol='tight')

Solve the variational optimization problem using IPOPT.

The current path parameters are flattened into a 1D variable vector x and passed to IPOPT. After optimization, the updated vector is written back via set_x. The return values are those provided by cyipopt.Problem.solve.

The argument tol provides a convenient shortcut for adjusting the IPOPT option dual_inf_tol using the presets from Ref. 1:

• 'tight' → dual_inf_tol = 0.04

• 'middle' → dual_inf_tol = 0.10

• 'loose' → dual_inf_tol = 0.20

• a float value directly sets dual_inf_tol to that number.

Parameters:

tol ({'tight', 'middle', 'loose'} or float, optional) – Desired dual infeasibility tolerance. Default is 'tight'.

Returns:

• x_opt (ndarray) – Optimized 1D variable array.

• info (dict) – IPOPT information dictionary.

get_forces()

Evaluate forces and energies for all images along the path.

For each image at t = t_eval[i], this method computes the atomic forces and potential energy using the calculator attached to the corresponding ase.Atoms object. Forces and energies at the endpoints are obtained from cached values (_f_ends and _e_ends) with a small tolerance region near t = 0 and t = 1. Interior images are evaluated serially, using threads, or using MPI depending on the settings of parallel and world.

After calling this method, the arrays self.forces and self.energies are updated in place.

Returns:

ndarray – Array of shape (len(t_eval), natoms, 3) containing the forces for all images along the current path.

Notes

• Endpoint forces are stored in _f_ends and are reused for t < eps_t and t > 1 - eps_t.

• If remove_rotation_and_translation is enabled, forces at the final endpoint are rotated to match the aligned coordinate frame.

• When MPI is used, each rank evaluates a subset of interior images; results are broadcast so that all ranks obtain full arrays.

intermediate(alg_mod, iter_count, obj_value, inf_pr, inf_du, mu, d_norm, regularization_size, alpha_du, alpha_pr, ls_trials)

IPOPT callback: per-iteration monitor.

This method is not intended to be called directly by users. It is invoked internally by IPOPT at the end of each iteration. The arguments are provided by IPOPT and follow its callback interface specification.

In addition to the default IPOPT behavior, this method records iteration-by-iteration quantities into self.history. See HistoryDMF for details.

If update_teval=True, this method also adaptively updates the energy evaluation points t_eval based on the current energy profile and dual infeasibility. See Ref. 1 for details.

Parameters:

• alg_mod – Values supplied directly by IPOPT at each iteration. These are passed through unchanged and are not meant to be modified by the user.

• iter_count – Values supplied directly by IPOPT at each iteration. These are passed through unchanged and are not meant to be modified by the user.

• obj_value – Values supplied directly by IPOPT at each iteration. These are passed through unchanged and are not meant to be modified by the user.

• inf_pr – Values supplied directly by IPOPT at each iteration. These are passed through unchanged and are not meant to be modified by the user.

• inf_du – Values supplied directly by IPOPT at each iteration. These are passed through unchanged and are not meant to be modified by the user.

• mu – Values supplied directly by IPOPT at each iteration. These are passed through unchanged and are not meant to be modified by the user.

• d_norm – Values supplied directly by IPOPT at each iteration. These are passed through unchanged and are not meant to be modified by the user.

• regularization_size – Values supplied directly by IPOPT at each iteration. These are passed through unchanged and are not meant to be modified by the user.

• alpha_du – Values supplied directly by IPOPT at each iteration. These are passed through unchanged and are not meant to be modified by the user.

• alpha_pr – Values supplied directly by IPOPT at each iteration. These are passed through unchanged and are not meant to be modified by the user.

• ls_trials – Values supplied directly by IPOPT at each iteration. These are passed through unchanged and are not meant to be modified by the user.

HistoryDMF

class dmf.dmf.HistoryDMF

Bases: object

Container storing the optimization history of the DirectMaxFlux method.

This object collects various physical and numerical quantities evaluated along the reaction path during the optimization. At each IPOPT iteration, the DirectMaxFlux.intermediate method appends the current values of these quantities to the corresponding lists below.

forces

History of DirectMaxFlux.forces.

Type:

list of ndarray

energies

History of DirectMaxFlux.energies.

Type:

list of ndarray

coefs

History of DirectMaxFlux.coefs.

Type:

list of ndarray

angs

History of DirectMaxFlux.angs.

Type:

list of ndarray

t_eval

History of DirectMaxFlux.t_eval.

Type:

list of ndarray

tmax

History of the location t_max corresponding to the maximum interpolated energy along the path. See Ref. 1 for details.

Type:

list of float

images_tmax

History of the atomic structure at t = t_max, providing an approximate transition-state geometry at each iteration.

Type:

list of ase.Atoms

duals

History of the scaled dual infeasibility (IPOPT diagnostic).

Type:

list of float