HTSTEvent

class HTSTEvent(neb, assumed_prefactor=None, minimum_displacement=PhysicalQuantity(0.05, Ang), finite_difference=PhysicalQuantity(0.01, Ang), finite_difference_method=<class 'NL.ComputerScienceUtilities.NLFlag._NLFlag.Central'>, constraints=None, saddle_image_index=None, spline_estimate=False, dynamical_matrix_log_filename_prefix='htst_')

Calculate harmonic transition state theory (HTST) reaction rates using configurations from a NudgedElasticBand object.

Parameters:
  • neb (NudgedElasticBand) – The nudged elastic band configuration to use.
  • assumed_prefactor (PhysicalQuantity of type per Second) – A fixed value for the prefactors (both forward and reverse). If assumed_prefactor is None, then the vibrational frequencies at the minima and saddle point will be calculated.
    Default: None
  • minimum_displacement (PhysicalQuantity of type length) – The minimum distance that an atom must move during the reaction in order to be included in the prefactor calculation. If an atom moves further than this distance, it and its neighbors will be included in the dynamical matrix calculation that is needed for calculating the prefactor. Setting this value to zero will include all degrees of freedom. Setting this value to a larger number will reduce the amount of worked needed to calculate the prefactor, but will decrease the accuracy. Note that this parameter is only used when assumed_prefactor is None.
    Default: 0.05 * Angstrom
  • finite_difference (PhysicalQuantity of type length) – The finite difference step size used when estimating the force constants.
    Default: 0.01*Angstrom
  • finite_difference_method (Central | Forward) – The finite difference method to use when estimating the force constants needed to calculate the harmonic prefactor.
    Default: Central
  • saddle_image_index (int) – The index of the image that is assumed to be the saddle point. If saddle_image_index is None, then the image with the largest energy will be chosen.
    Default: None
  • spline_estimate (bool) – Obtain a rough estimate of the prefactor using a spline to obtain the vibrational frequency at the minimum in the direction of the saddle point. This requires that the NEB configuration has either had its update method called or has already been optimized so that the energy and forces of each image are available. This option is incompatible with assumed_prefactor.
    Default: False
  • dynamical_matrix_log_filename_prefix (str | None) – Prefix for the filenames where the logging output for every displacement calculation is stored. The filenames are formed by appending a number and the file extension (".log"). A value of None means all logging will go to stdout.
    Default: "htst_"
forwardBarrier()
Returns:The forward (reactant to product) barrier.
Return type:PhysicalQuantity of type energy
forwardPrefactor()
Returns:The forward (reactant to product) vibrational prefactor.
Return type:PhysicalQuantity of type per second
forwardRate(temperature)

Calculate the forward (reactant to product) reaction rate using HTST.

Parameters:temperature (PhysicalQuantity of type temperature) – The reaction temperature.
Returns:The forward reaction rate.
Return type:PhysicalQuantity of type per second
classmethod fromConfigurations(reactant_configuration, saddle_configuration, product_configuration, assumed_prefactor=None, minimum_displacement=PhysicalQuantity(0.05, Ang), finite_difference=PhysicalQuantity(0.01, Ang), finite_difference_method=<class 'NL.ComputerScienceUtilities.NLFlag._NLFlag.Central'>, dynamical_matrix_log_filename_prefix='htst_')

This is an alternate constructor for this class that uses explicit configurations for the reactant minimum configuration, saddle point configuration, and product minimum configuration as opposed to extracting these configurations from a NudgedElasticBand object.

Parameters:
  • reactant_configuration (MoleculeConfiguration | BulkConfiguration | DeviceConfiguration | SurfaceConfiguration) – The configuration at the reactant/initial minimum.
  • saddle_configuration (MoleculeConfiguration | BulkConfiguration | DeviceConfiguration | SurfaceConfiguration) – The configuration at the saddle point.
  • product_configuration (MoleculeConfiguration | BulkConfiguration | DeviceConfiguration | SurfaceConfiguration) – The configuration at the product/final minimum.
  • assumed_prefactor (PhysicalQuantity of type per Second) – A fixed value for the prefactors (both forward and reverse). If assumed_prefactor is None, then the vibrational frequencies at the minima and saddle point will be calculated.
    Default: None
  • minimum_displacement (PhysicalQuantity of type length) – The minimum distance that an atom must move during the reaction in order to be included in the prefactor calculation. If an atom moves further than this distance, it and its neighbors will be included in the dynamical matrix calculation that is needed for calculating the prefactor. Setting this value to zero will include all degrees of freedom. Setting this value to a larger number will reduce the amount of worked needed to calculate the prefactor, but will decrease the accuracy. Note that this parameter is only used when assumed_prefactor is None.
    Default: 0.05*Angstrom
  • finite_difference (PhysicalQuantity of type length) – The finite difference step size used when estimating the force constants.
    Default: 0.01*Angstrom
  • finite_difference_method (Central | Forward) – The finite difference method to use when estimating the force constants needed to calculate the harmonic prefactor.
    Default: Central
  • dynamical_matrix_log_filename_prefix (str | None) – Prefix for the filenames where the logging output for every displacement calculation is stored. The filenames are formed by appending a number and the file extension (”.log”). A value of None means all logging will go to stdout.
    Default: "htst_"
static htstRate(barrier, temperature, prefactor)

Calculate the HTST reaction rate at the specified temperature.

Parameters:
  • barrier (PhysicalQuantity of type energy) – The energy barrier.
  • temperature (PhysicalQuantity of type temperature) – The reaction temperature.
  • prefactor (PhysicalQuantity of type per second) – The harmonic prefactor (attempt frequency).
Returns:

The HTST reaction rate.

Return type:

PhysicalQuantity of type per second

metatext()
Returns:The metatext of the object or None if no metatext is present.
Return type:str | unicode | None
nlprint(stream=None)

Print a string containing an ASCII table useful for plotting the AnalysisSpin object.

Parameters:stream (python stream) – The stream the table should be written to.
Default: NLPrintLogger()
productConfiguration()
Returns:The product configuration.
Return type:MoleculeConfiguration | BulkConfiguration | DeviceConfiguration | SurfaceConfiguration
productEnergy()
Returns:The energy of the product configuration.
Return type:PhysicalQuantity of type energy
productFrequencies()
Returns:The angular vibrational frequencies at the product configuration.
Return type:PhysicalQuantity of type radians per second
reactantConfiguration()
Returns:The reactant configuration.
Return type:MoleculeConfiguration | BulkConfiguration | DeviceConfiguration | SurfaceConfiguration
reactantEnergy()
Returns:The energy of the reactant configuration.
Return type:PhysicalQuantity of type energy
reactantFrequencies()
Returns:The angular vibrational frequencies at the reactant configuration.
Return type:PhysicalQuantity of type radians per second
reverseBarrier()
Returns:The reverse (product to reactant) barrier.
Return type:PhysicalQuantity of type energy
reversePrefactor()
Returns:The reverse (product to reactant) vibrational prefactor.
Return type:PhysicalQuantity of type per second
reverseRate(temperature)

Calculate the reverse (product to reactant) reaction rate using HTST.

Parameters:temperature (PhysicalQuantity of type temperature) – The reaction temperature.
Returns:The reverse reaction rate.
Return type:PhysicalQuantity of type per second
saddleConfiguration()
Returns:The saddle configuration.
Return type:MoleculeConfiguration | BulkConfiguration | DeviceConfiguration | SurfaceConfiguration
saddleEnergy()
Returns:The energy of the saddle configuration.
Return type:PhysicalQuantity of type energy
saddleFrequencies()
Returns:The angular vibrational frequencies at the saddle configuration.
Return type:PhysicalQuantity of type radians per second
setMetatext(metatext)

Set a given metatext string on the object.

Parameters:metatext (str | unicode | None) – The metatext string that should be set. A value of “None” can be given to remove the current metatext.

Logging

If the dynamical matrix is calculated at the minima and saddle points, the output from the calculator for each displacement is written to a separate file on the disk. By default, the displacements at the initial minimum will have the pattern htst_reactant_displacement_%i.log where the %i is an integer representing the displacement number. The displacements and the saddle are named htst_saddle_displacement_%i.log and displacements at the final minimum are named htst_product_displacement_%i.log. The prefix htst_ can be modified by changing the dynamical_matrix_log_filename_prefix parameter. Alternatively, all file-based logging can be disabled by setting dynamical_matrix_log_filename_prefix to None.

Notes

This class calculates reaction rates using harmonic transition state theory (HTST),

\[k_{\rm HTST} = \frac{\prod_i^{3N} \nu_i^{\rm min}}{\prod_i^{3N-1} \nu_i^\ddagger} \exp \left[ -\left(E^\ddagger - E^{\rm min}\right)/k_{\rm B} T \right],\]

where \(k_{\rm HTST}\) is the HTST reaction rate, \(N\) is the number of atoms, \(\nu^{\rm min}_i\) and \(\nu^\ddagger_i\) are the positive normal mode frequencies at the minimum and at the saddle point, \(E^{\rm min}\) \(E^\ddagger\) are the energies at the minimum and at the saddle point, \(k_{\rm B}\) is Boltzmann’s constant and \(T\) is the temperature. The ratio of the product of the normal mode frequencies (the term in front of the exponential) is called the prefactor or the attempt frequency.

Determining the prefactor is an expensive calculation. Formally, it requires calculating the eigenvalues of the dynamical matrix (without any repetitions) of the minima and of the saddle point. However, it can be approximated by reducing the dimensionality of the system to include only those atoms that move between the minimum and saddle point and their near neighbors. This approximation can be controlled using the minimum_displacement parameter. The default value of is 0.05 Angstrom and typically estimates the prefactor to within 5-10% of the true value. Note that the prefactor is temperature independent. This means that it only needs to be calculated once and then all calls to forwardRate(temperature) and reverseRate(temperature) can compute the rate inexpensively.

When calculating the eigenvalues of the dynamical matrix, the finite_difference_method parameter controls the order of the finite difference method used to compute the force constants. The finite_difference parameter controls the step length. The Forward method uses a first order finite difference scheme that only needs one displacement per degree of freedom, while the Central method uses a second order finite difference scheme that needs two displacements per degree of freedom, but is much more accurate. The central difference scheme is recommended when using DFT due to numerical noise that arises from having a discrete grid.

Alternatively, one can disable the prefactor calculation entirely by passing a frequency using the assume_prefactor parameter. In most systems, this value does not change significantly between different reactions and mostly depends upon the mass of the atoms involved. For lighter elements (about the first three rows of the periodic table) the prefactor is approximately \(10^{13}\,{\rm s}^{-1}\) while for heavier elements, it about \(10^{12}\,{\rm s}^{-1}\).

Lastly, there is a very crude way to estimate the prefactor by using a spline interpolation along the NEB reaction coordinate to determine the vibrational frequency in the direction of the saddle point. This option is enabled by setting spline_estimate to True. This very simple approach often gives a reasonable estimate of the prefactor, but is not guaranteed to give any particular level of accuracy and should not be relied upon.