11. Chain Of States Methods

A chain of states (COS) comprises a set of distinct states (images) and is usually spanned between two minima on a potential energy surface (PES).

When properly relaxed, a COS coincides with a minimum energy path (MEP), or is a good approximation to it. Tangents can be defined for every COS image and together they make up a discretized path that describes the reaction/chemical transformation. The tangents are also used to divide the COS gradient into perpendicular and parallel components, w.r.t the tangents. Relaxation (optimization) of the COS is achieved by minimizing the perpendicular component of the gradient, so only the parallel component remains.

pysisyphus offers different COS implementations, namely Nudged Elasic Band (NEB), including its Adaptive and Free-End (and Free-End-Adaptive) modificiations and different flavors of String methods (Growing String Method, GSM, and Simple Zero Temperature String, SZTS). The GSM implementation is also available for with internal coordinates.

11.1. String method

(When discussing String methods the COS 'images' are usually called 'nodes'.)

In the string method the whole COS is periodically (every n-th cycle) reparametrized by fitting a spline through all nodes and redistributing the nodes them along the spline, according to a predefined parametrization. By choosing between different parametrizations equal node spacing or higher resolution around the highest energy image (HEI) can be achieved. The tangents needed for the gradient projection are obtained as first derivatives of the spline.

Reparametrization every \(n\)-th cycle impedes efficient string optimization, and prevents the use of optimizers with some kind of history Conjugate Gradient (CG). The optimizer history is reset after each reparametrization and a simple Steepest Descent (SD) step is done after reparametrization.

11.2. Nudged Elastic Band

No reparametrization takes place in the NEB method. The parallel gradient component along the tangent is projected out and replaced by an artificial spring force. In principle, optimizing NEBs should be easier as there is no reparametrization and more sophisticated optimizers beyond SD can and should be employed.

11.3. General remarks

Converged COS produce good guesses for subsequent TS searches, when the (splined) HEI is determined. In pysisyphus subsequent TS searches are easily started by including tsopt: in the YAML input. Knowledge of the initial and final images of the COS is used to construct a more complete set of internal coordinates for the TS guess and it is less likely that important coordinates are missed. The initial imaginary mode to follow uphill is selected as the one featuring the highest overlap with HEI tangent.

11.4. YAML example(s)

Below you can find an example YAML-input including the most important options that the user may want to modify when running a GSM optimization.

precontr:                                # Preconditioning of translation & rotation
preopt:                                  # Preoptimize inital and final geometry
cos:
 type: gs                                # Do a growing string
 max_nodes: 9                            # Total string will have 9 + 2 == 11 images
 climb: False                            # Enable climbing image (CI), usually a good idea.
 climb_rms: 0.005                        # rms(forces) threshold for enabling CI
 climb_lanczos: False                    # Use tangent obtained from Lanczos algorithm for CI.
 climb_lanczos_rms: 0.005                # rms(forces) threshold for enabling Lanczos algorithm.
 reparam_check: rms                      # Criterian for growing new frontier nodes (rms/norm).
 perp_thresh: 0.05                       # Threshold for growing new frontier nodes.
 reparam_every: 2                        # Reparametrize every n-th cycle when NOT fully grown.
 reparam_every_full: 3                   # Reparametrize every n-th cycle when fully grown.
opt:
 type: string                            # Optimizer for GrowingString
 stop_in_when_full: -1                   # Stop string optimization N cycles after fully grown
                                         # Disabled by -1. Usually it's a better idea to
                                         # further converge the string, after is fully grown.
 align: False                            # Disable Kabsch algorithm. Should be True with
                                         # type: cart
 scale_step: global                      # Scale down step as whole (global) or per image
                                         # (per_image)
tsopt:
 type: rsprfo                            # Continue with TS-optimization of highest energy images
                                         # (HEI) using the RS-P-RFO algorithm
 do_hess: True                           # Calculate hessian at optimized TS geometry
 trust_max: 0.3
 thresh: gau_loose
calc:
 type: orca
 keywords: "b3lyp 6-31G* rijcosx"
 pal: 4
 charge: 0
 mult: 1
geom:
 type: dlc
 fn: [first_preopt.xyz, last_preopt.xyz] # Run GrowingString in delocalized internal coordinates
                                         # (preferred).

For NEB optimizations a different optimizer (not type: string) should be used, e.g., QuickMin type: qm in the beginning, and type: lbfgs later on. It is not yet possible to specify two different optimizers that are used in stages, so if is desired it must be done manually.

# Taken from examples/complex/06_diels_alder...
geom:
 type: cart
 fn: diels_alder.trj
calc:
 type: xtb
 charge: 0
 mult: 1
 pal: 4
preopt:
 max_cycles: 5
interpol:                       # In NEBs the whole path is interpolated beforehand.
 type: redund                   # Possible values: redund|idpp|lst|linear
 between: 10                    # Interpolate n-geometries between every pair of supplied
                                # geometries. For two this yields 1 + 10 + 1 == 12 images,
                                # for three geometries this yields 1 + 10 + 1 + 10 + 1 == 23
                                # geometries.
cos:
 type: neb
 climb: False
 align_fixed: True              # Align the fixed atoms of the initial and final images along the path.
opt:
 type: lbfgs
 align: True                   # Align the image of the current step with the image of the previous step
 align_factor: 0.9             # If full alignment is not desired, a factor between 0 and 1
                               # can be specified.
 rms_force: 0.01
 max_step: 0.04
tsopt:
 type: rsirfo
 do_hess: True
 max_cycles: 75
 thresh: gau_tight
 hessian_recalc: 7
irc:
 type: eulerpc
 rms_grad_thresh: 0.0005
endopt:

Further examples for COS optimizations from .yaml input can be found here.

11.5. General advice for COS optimizations

  • Start from optimized geometries or use the preopt: key in the YAML input.

  • Consider fixing the initial and final images fix_ends: True

  • Always use align: True when optimizing a COS in cartesian coordinates to remove translation and rotation. align: True must not be used when running a COS with DLC.

  • Don't over-converge a COS. It is usually a better idea to converge a COS loosely and use the highest energy image (HEI) as a guess for a subsequent transition state (TS) search.

  • If possible use a climbing image (CI) climb: True

  • When running a growing string calculation (type: gs) use stop_in_when_full: [n] in the opt: section with a small integer [n] to stop the COS relaxation after the string is fully grown

11.6. Chain Of States base class

Base class for chain of state methods

class pysisyphus.cos.ChainOfStates.ChainOfStates(images, fix_first=True, fix_last=True, align_fixed=True, climb=False, climb_rms=0.005, climb_lanczos=False, climb_lanczos_rms=0.005, climb_fixed=True, energy_min_mix=False, scheduler=None, progress=False)[source]
as_xyz(comments=None)[source]
property atoms
calculate_forces()[source]
property calculator
property cart_coords

Return a flat 1d array containing the cartesian coordinates of all images.

check_for_climbing_start(ref_rms)[source]
clear()[source]
concurrent_force_calcs(images_to_calculate, image_indices)[source]
property coords

Return a flat 1d array containing the coordinates of all images.

property coords3d
describe()[source]
property energy
property forces
get_climbing_forces(ind)[source]
get_climbing_indices()[source]
get_dask_client()[source]
get_fixed_indices()[source]
get_hei_index(energies=None)[source]

Return index of highest energy image.

get_image_calc_counter_sum()[source]
get_perpendicular_forces(i)[source]

[1] Eq. 12

get_splined_hei()[source]
get_tangent(i, kind='upwinding', lanczos_guess=None, disable_lanczos=False)[source]

[1] Equations (8) - (11)

get_tangents()[source]
property gradient
property image_coords
property image_inds
property last_index
log(message)[source]
logger = <Logger cos (DEBUG)>
property masses_rep
property max_image_num
property moving_images
property moving_indices

Returns the indices of the images that aren't fixed and can be optimized.

par_image_calc(image)[source]
property perpendicular_forces
prepare_opt_cycle(last_coords, last_energies, last_forces)[source]

Implements additional logic in preparation of the next optimization cycle.

Should be called by the optimizer at the beginning of a new optimization cycle. Can be used to implement additional logic as needed for AdaptiveNEB etc.

property results
rms(arr)[source]

Root mean square

Returns the root mean square of the given array.

Parameters:

arr (iterable of numbers) --

Returns:

rms -- Root mean square of the given array.

Return type:

float

set_climbing_forces(forces)[source]
set_coords_at(i, coords)[source]

Called from helpers.procrustes with cartesian coordinates. Then tries to set cartesian coordinate as self.images[i].coords which will raise an error when coord_type != "cart".

set_images(indices, images)[source]
set_vector(name, vector, clear=False)[source]
set_zero_forces_for_fixed_images()[source]

This is always done in cartesian coordinates, independent of the actual coord_type of the images as setting forces only work with cartesian forces.

valid_coord_types = ('cart', 'cartesian', 'dlc')
zero_fixed_vector(vector)[source]

11.7. Chain Of State Methods

11.7.1. Nudged Elastic Band (NEB)

class pysisyphus.cos.NEB.NEB(images, variable_springs=False, k_max=0.3, k_min=0.1, perp_spring_forces=None, bandwidth=None, **kwargs)[source]

Bases: ChainOfStates

fmt_k()[source]
property forces
get_parallel_forces(i)[source]
get_quenched_dneb_forces(i)[source]

See [3], Sec. VI and [4] Sec. D.

get_spring_forces(i)[source]
property parallel_forces
set_variable_springs()[source]
update_springs()[source]

11.7.2. Adaptive NEB

class pysisyphus.cos.AdaptiveNEB.AdaptiveNEB(images, adapt=True, adapt_fact=0.25, adapt_between=1, scale_fact=False, keep_hei=True, free_ends=True, **kwargs)[source]

Bases: NEB

__init__(images, adapt=True, adapt_fact=0.25, adapt_between=1, scale_fact=False, keep_hei=True, free_ends=True, **kwargs)[source]

(Free-End) Adaptive Nudged Elastic Band.

Parameters:
  • images (list of Geometry objects) -- Images of the band.

  • adapt (bool, default True) -- Whether to adapt the image number or not. This switch is included to support the FreeEndNEB class, that is just a thin wrapper around this class.

  • adapt_fact (positive float) -- Factor that is used to decide wether to adapt. The inital threshold is calculated by multiplying the RMS force of the band with this factor. When the RMS of the force falls below this threshold adaption takes place.

  • adapat_between (positive integer) -- Number of images to interpolate between the highest energy image and its neighbours. The number of starting images must be higher or equal then 2*adapt_between+3, as we reuse/transfer the calculators from the starting images onto the new ones.

  • scale_fact (bool, default False) -- Whether to increase adapt_fact in deeper levels. This may lead to earlier adapation.

  • keep_hei (bool, optional) -- Whether to keep the highest energy image (usually a very good idea) or to interpolate only between the neighbouring images.

  • free_ends (bool, default True) -- Whether to use modified forces on the end images.

adapt_this_cycle(forces)[source]

Decide wether to adapt.

Parameters:

forces (np.array) -- Forces of the previous optimization cycle.

Returns:

adapt -- Flag that indicates if adaption should take place in this cycle.

Return type:

bool

property forces

See Eq. (7) in [2].

prepare_opt_cycle(*args, **kwargs)[source]

Check for adaption and adapt if needed.

See ChainOfStates.prepare_opt_cycle for a complete docstring.

update_adapt_thresh(forces)[source]

Update the adaption threshold.

Parameters:

forces (np.array) -- Forces of the previous optimization cycle.

11.7.3. Free-End NEB

class pysisyphus.cos.FreeEndNEB.FreeEndNEB(*args, fix_first=False, fix_last=False, **kwargs)[source]

Bases: AdaptiveNEB

__init__(*args, fix_first=False, fix_last=False, **kwargs)[source]

Simple Free-End-NEB method.

Derived from AdaptiveNEB with disabled adaptation. Only implements Eq. (7) from [2]. For other implementations please see the commit 01bc8812ca6f1cd3645d43e0337d9e3c5fb0ba55. There the other variants are present but I think Eq. (7) in [2] is the simplest & best bet.

11.7.4. Simple Zero-Temperature String

class pysisyphus.cos.SimpleZTS.SimpleZTS(images, param='equal', **kwargs)[source]

Bases: ChainOfStates

reparametrize()[source]

11.8. Growing Chain Of States base class

Base class for growing chain of state methods

class pysisyphus.cos.GrowingChainOfStates.GrowingChainOfStates(images, calc_getter, max_nodes=10, **kwargs)[source]
property arc_dims
get_new_image_from_coords(coords, index)[source]
property max_image_num
new_node_coords(k)[source]
prepare_opt_cycle(*args, **kwargs)[source]

Implements additional logic in preparation of the next optimization cycle.

Should be called by the optimizer at the beginning of a new optimization cycle. Can be used to implement additional logic as needed for AdaptiveNEB etc.

set_new_node(k)[source]

11.9. Growing Chain Of State Methods

11.9.1. Growing String Method

class pysisyphus.cos.GrowingString.GrowingString(images, calc_getter, perp_thresh=0.05, param='equi', reparam_every=2, reparam_every_full=3, reparam_tol=None, reparam_check='rms', max_micro_cycles=5, reset_dlc=True, climb=False, **kwargs)[source]

Bases: GrowingChainOfStates

property forces
property full_string_image_inds
property fully_grown

Returns wether the string is fully grown. Don't count the first and last node.

get_additional_print()[source]
get_cur_param_density(kind=None)[source]
get_new_image(ref_index)[source]

Get new image by taking a step from self.images[ref_index] towards the center of the string.

get_tangent(i)[source]

[1] Equations (8) - (11)

property image_inds
property left_size
property lf_ind

Index of the left frontier node in self.images.

property nodes_missing

Returns the number of nodes to be grown.

reparam_cart(desired_param_density)[source]
reparam_dlc(desired_param_density, thresh=0.001)[source]
reparametrize()[source]
reset_geometries(ref_geometry)[source]
property rf_ind

Index of the right frontier node in self.images.

property right_size
set_coords(image, coords)[source]
spline(tangents=False)[source]
property string_size