CollectiveVariableHyperdynamics

class CollectiveVariableHyperdynamics(distortions, temperature=None, global_exponent=None, global_cut=None, gaussian_width=None, gaussian_height=None, gaussian_frequency=None, bias_damping_temperature=None, gaussian_limit=None, reaction_steps=None, optimize_new_state=None, measurement_frequency=None, tracked_atoms_tag=None, allow_no_variables=None)

Hook function implementing collective variable hyperdynamics.

Parameters:
  • distortions (Sequence of HyperdynamicsBaseDistortion) – The geometric distortions used to create the collective variable.

  • temperature (PhysicalQuantity of type temperature | None) – The temperature of the simulation, used for calculating the hypertime.
    Default: 300 Kelvin

  • global_exponent (float) – The exponent used to create the global distortion.
    Default: 6

  • global_cut – The value of the global distortion at the transition state. This can be used to set the position of the transition state where individual distortions do not range from 0 to 1. One example of this is when using a bond distortion without a maximum bond length.
    Default: 1

  • gaussian_width (float) – The width of the added Gaussian bias potential.
    Default: 0.025

  • gaussian_height (PhysicalQuantity of type energy) – The height of the added Gaussian bias potential.
    Default: 0.01 eV

  • gaussian_frequency (int) – The number of steps between adding new Gaussian bias potentials.
    Default: 1000

  • bias_damping_temperature (PhysicalQuantity of type temperature | None) – The temperature used to scale the bias potential using well-temperated metadynamics. None does not damp the added potential, resulting the in same size potential being added at each step.
    Default: None

  • gaussian_limit (float) – The upper limit in the collective variable where Gaussians are added to the hyperdynamic bias potential. Gaussians are not added at values above this limit in order to avoid biasing the transition state.
    Default: 1.0

  • reaction_steps (int) – The number of steps at the maximum collective variable before a reaction takes place.
    Default: 5000

  • optimize_new_state (bool) – Use an optimized configuration to reset the distortions after a reaction is detected.
    Default: True

  • measurement_frequency (int) – The step frequency with which measurements are added to the trajectory.
    Default: 10

  • tracked_atoms_tag (str) – Tag giving the indices of atoms for which the position is measured during the simulation. None results in no atoms being tracked.
    Default: None

  • allow_no_variables (bool) – Continue in cases where no variables are found for acceleration. The simulation continues the number of reaction steps without bias before trying again to find variables to accelerate.
    Default: True

allowNoVariables()
Returns:

True if the simulation continues in cases where no variables are found for acceleration.

Return type:

bool

biasDampingTemperature()
Returns:

The temperature used to scale the bias potential.

Return type:

PhysicalQuantity of type temperature

callInterval()
Returns:

The call interval of this hook function.

Return type:

int

distortions()
Returns:

The geometric distortions used to create the collective variable.

Return type:

Sequence of HyperdynamicsBaseDistortion

gaussianCenters()
Returns:

The center of the Gaussian potentials added to form the bias potential.

Return type:

ndarray

gaussianFrequency()
Returns:

The frequency at which new Gaussian bias potentials are added.

Return type:

int

gaussianHeight()
Returns:

The height of the Gaussian bias potential.

Return type:

PhysicalQuantity of type energy

gaussianHeights()
Returns:

The height of the Gaussian potentials added to form the bias potential.

Return type:

PhysicalQuantity of type energy

gaussianLimit()
Returns:

The limit in the collective variable above which Gaussians are not added.

Return type:

float

gaussianWidth()
Returns:

The width of the Gaussian bias potential.

Return type:

float

gaussianWidths()
Returns:

The widths of the Gaussian potentials added to form the bias potential.

Return type:

ndarray

globalCut()
Returns:

The value of the global distortion at the transition state.

Return type:

float

globalExponent()
Returns:

The global exponent used in the collective variable.

Return type:

float

measurementFrequency()
Returns:

The step frequency with which measurements are saved to the trajectory.

Return type:

int

numberOfGaussians()
Returns:

The number of Gaussian potentials added to form the hyperdynamic bias potential.

Return type:

int | None

numberOfVariables()
Returns:

The total number of geometric coordinates being accelerated. If the geometric coordinates are not set with a configuration None is returned.

Return type:

int | None

optimizeNewState()
Returns:

True if the configuration is optimized before setting new distortions.

Return type:

bool

reactionSteps()
Returns:

The number of steps at the maximum collective variable before a reaction takes place.

Return type:

int

temperature()
Returns:

The simulation temperature.

Return type:

PhysicalQuantity of type temperature

trackedAtomsTag()
Returns:

The tag for the atoms that are tracked during the simulation.

Return type:

str

uniqueString()

Return a unique string representing the state of the object.

Usage Examples

In this example collective variable hyperdynamics is used to study the diffusion of a single hydrogen atom in a copper crystal lattice. A position distortion is used to add a bias potential to push the atom through the lattice. The results of the simulation are recorded in the resulting MDTrajectory and can be viewed in the Movie Tool. The resulting diffusivity can also be calculated using the IrregularTimeMeanSquareDisplacement object.

../../../_images/cu_h_movie_tool.png

This simulation may take several hours to run. Note that the simulation performs a limited amount of sampling, using only a million steps or 1ns of simulation time. The hypertime of the simulation however is over 250ns, meaning that the simulation was able to effectively sample this time for the cost of only 1ns of simulation. At 300 Kelvin the diffusion of hydrogen in copper is around \(8.5 \times 10^{-14} m^2/s\) so very long trajectories are needed to measure the slow diffusion[1]. For accurate results it is necessary to use either more steps or to run multiple simulations. Running multiple simulations can enable a faster turn-around time, as the simulations can be run parametrically in parallel. The workflow for this example can be downloaded here with the corresponding Python script CVHD_CuH_Diffusion_Example.py.

Notes

Hyperdynamics is a method for sampling rare events using molecular dynamics. In hyperdynamics a bias potential is added to degrees of freedom separating different states. The bias potential aids in lowering the energy barrier between these states, allowing the system to sample different states more freely. In this way the biased degrees of freedom move at elevated temperature, while the rest of the system moves at the simulation temperature. The acceleration factor given to the system by the bias potential can be calculated as the ensemble average of the bias potential. When this boost factor is multiplied by the simulation time it gives the effective sampling time, known as the hypertime.

In using the hypertime to estimate reaction rates, it is important that the bias potential is not applied at the transition state. If bias is applied both at the reactant and transition state it obscures the relative free energy barrier between two states. This can lead to errors in estimations of the rate constant. This also applies to kinetic based processes such as diffusion.

In collective variable hyperdynamics the bias potential is built up during the simulation by adding Gaussian functions to the potential at the value of the collective variable during the simulation[2]. By building up the potential it avoids requiring knowledge of the the appropriate shape or height of the bias potential. This is distinct from StaticHyperdynamics, which requires some knowledge of the transition state barrier height. In this way collective variable hyperdynamics is similar to metadynamics. The main difference is that in metadynamics a single bias potential is built up during the simulation. This is converged to create an estimate of the total free energy of a process in the system. In hyperdynamics the bias potential is simply used to accelerate a defined reaction. Once the reaction occurs, the bias potential is reset to zero and is built up again at the new reactant. Rather than the free energy profile, the important results from hyperdynamics simulations are the trajectory and the effective hypertime of the simulation.

In order to create the bias potential, one or more important geometric features, known as distortions, are first identified in the configuration. These can be things such as the position of atoms or the bond lengths between specific atoms. These distortions are generally constructed so that they range from 0 to 1 when moving from reactant to transition state. The different distortions are then combined into a global distortion through the equation:

\[\chi_t = \left( \sum_{i=1}^{N} \chi_i^p \right)^{1/p}\]

Here \(\chi_t\) is the global distortion and \(\chi_i\) are the individual distortions. The power \(p\) increases the contribution of the larger distortions to the global distortion. In the CollectiveVariableHyperdynamics class the distortions are given as a list using the distortions argument. These can be either HyperdynamicsPositionDistortion objects or HyperdynamicsBondDistortion objects. The exponent \(p\) is set using the global_exponent keyword argument.

To ensure that the gradients at the boundaries are zero, the collective variable is then defined using a cosine function such that:

\[\begin{split}\eta = \begin{cases} \frac{1}{2} \left[1 - \cos\left(\pi \left(\frac{\chi_t}{\chi_c}\right)^2\right) \right] & \chi_t \le 1\\ 1 & \chi_t > 1 \\ \end{cases}\end{split}\]

Here \(\eta\) is the collective variable and \(\chi_c\) is a cutoff value of the global distortion. This is the value of the global distortion at the transition state. In the simplest case in which one distortion that ranges from 0 to 1 is used, this is simply one. If a combination of distortions are used, or if the distortions do not range from 0 to 1 the \(\chi_c\) is used to make the collective variable range from 0 to 1. This cutoff value is set using the global_cut keyword argument. The collective variable is also used to determine when a reaction has occurred in the configuration. When the collective variable is greater than 1 for a set number of steps the system is considered to have reacted, and the references for distortions and the bias potential is reset to zero. The number of steps used to determine if a reaction has occurred is set with the reaction_steps keyword argument. By default this is set to 5000 steps. In setting the new distortions the reference configuration can first be optimized using the current calculator. This ensures that the distortions are set from a local energy minimum. In the case of bonds it ensures that the correct bonds are selected for boosting. Likewise in the case of positions that the reference positions are at a potential minimum. This optimized reference configuration is not used in the dynamics simulation, as to do so would alter the sampling of different states. It is in general recommended to optimize the reference configuration. Enabling optimization after reactions is done with the optimize_new_state keyword argument, which is True by default.

After each reaction the distortion definitions are reset to the new reference configuration. In some cases, such as bonds being boosted, it means that a different number of distortions may be selected. This can also include cases where no available distortions are found that can be accelerated. In this case the simulation can simply continue the number of steps between reactions and try again to find some applicable distortions. Alternatively the simulation can simply end if no distortions indicates that something has gone wrong with the simulation. Stopping on no distortions being found is controlled with the allow_no_variables keyword argument. By default this is set to True, meaning that the simulation will continue.

The total bias potential \(\Delta V(\eta)\) is calculated by adding Gaussians at time intervals \(\tau_G\). This gives the expression for the total bias potential:

\[\Delta V(\eta) = \sum_{k=1}^{n_G} w_k \exp\left[-\frac{\left(\eta - \eta(k\tau_G)\right)^2}{2\delta^2}\right]\]

Here \(w_k\) is the height of the Gaussian and \(\delta\) is the Gaussian width. These can be set with the keyword arguments gaussian_height and gaussian_width respectively. The number of steps between adding Gaussians is also set with the gaussian_frequency keyword argument. Together with the Gaussian height and width this sets the rate of energy deposition into the bias during the simulation. The optimal energy deposition rate depends on the depth of the local potential well and the temperature of the simulation. The energy deposition rate should be low enough to allow for effective sampling of the reactant potential energy surface so that a reliable estimate of the hypertime is obtained.

The height of each added Gaussian can also be scaled down depending on the amount of bias already added at the collective variable.This is similar to the well-tempered metadynamics method, where the bias potential converges to a value based on a temperature used to set the rate of reducing the size of the added Gaussians. Using this method the temperature is proportional to the maximum effective temperature of the motion of the collective variable. Reducing the temperature reduces the energy deposition rate into the bias potential. The temperature used to set the height of the Gaussians is given with the bias_damping_temperature keyword argument. Here None means that no damping is done and that all added Gaussians have the same height. In order to avoid adding bias at the transition state, Gaussians can be added up to a limit in the collective variable. This limit depends on the width of the Gaussians and how close to the true transition state the collective variable extends to. This limit is set with the gaussian_limit keyword argument.

For each reaction the hypertime is calculated as the simulation time multiplied by the ensemble average of the bias potential. This is given as:

\[t_h = t_s \int_0^{t_s} \exp(\Delta V(\eta, t) / k_B T) dt\]

Here \(t_h\) is the hypertime, \(t_s\) is the simulation time and \(T\) is the simulation temperature. This is set with the temperature keyword argument. Note that as the bias potential is built up over the simulation, it is dependent on the time. When a new reaction is found the bias potential starts off as zero. At this point both the simulation time and hypertime move at the same rate. As the bias potential is built up the hypertime moves increasingly faster until a new reaction is found and the bias potential is reset. To get an accurate measure of the hypertime it is important to have sufficient sampling of areas of high bias. This can limit the bias energy deposition rate. If is is too high sufficient sampling of the bias potential may not occur.

When running a collective variable hyperdynamics simulation in MolecularDynamics, information about the simulation, such as the hypertime, current bias potential and collective variable value are added to the trajectory as measurements. These are shown in the Movie Tool when opening the trajectory. The number of steps between adding measurements is set using the measurement_frequency keyword argument. Additionally the positions of select atoms can be added as a measurement. This is useful in cases where the diffusion of single atom impurities is being studied. These are added with the measurement name tracked_atoms. As this is an array of atomic coordinates it cannot be displayed in the Movie Tool. The tag of the atoms to track is given with the tracked_atoms_tag keyword argument. If no tag is given atom tracking is not added to the trajectory.

To calculate diffusion from a hyperdynamics simulation the IrregularTimeMeanSquareDisplacement class can be used. This class can access both the hypertime and tracked atomic positions from the trajectory and calculated the resulting mean square displacement of the tracked atoms. As diffusion is a stochastic process, it is often useful to run several simulations to fully sample the motion of the atoms. This is especially true for amorphous systems or systems lacking long-range order.