AdaptiveKineticMonteCarlo

class AdaptiveKineticMonteCarlo(markov_chain, kmc_temperature, md_temperature, calculator, constraints, kmc=None, superbasin_threshold=None, saddle_search_parameters=None, htst_parameters=None, confidence=None, filename_prefix=None, write_log=None, write_searches=None, write_markov_chain=None, write_kmc=None, lkmc_parameters=None)

Class for representing AdaptiveKineticMonteCarlo simulations.

Parameters:

• markov_chain (MarkovChain) – The Markov chain.

• kmc_temperature (PhysicalQuantity of type temperature) – The temperature of the kinetic Monte Carlo simulation.

• md_temperature (PhysicalQuantity of type temperature) – The reservoir temperature for the molecular dynamics search.

• calculator (Calculator) – The calculator object that should be attached to the configuration.

• constraints (list of type int) – The atoms that should be constrained during optimization. At least one atom must be constrained.

• kmc (KineticMonteCarlo) – Optional, existing kinetic Monte Carlo simulation, useful for resuming a calculation.
  Default: None.

• superbasin_threshold (int) – The number of steps between two states before they are merged together and treated with the coarse graining algorithm. A value of None disables coarse graining.
  Default: None.

• saddle_search_parameters (SaddleSearchParameters) – The parameters for the saddle search.

• htst_parameters (HTSTParameters) – The parameters used when calculating the reaction rate based on harmonic transition-state-theory (HTST). By default the prefactors will be approximated by only including atoms that move more than 0.05 Angstrom and the force constants will be calculated using a central finite difference scheme. A fixed prefactor can instead be assumed for all reactions: HTSTParameters(assumed_prefactor=1e12/Second)
  Default: None.

• confidence (float) – The confidence in the kinetic Monte Carlo rate table for each state in the simulation. This confidence is the probability that the “correct” kinetic Monte Carlo step was taken. The confidence should be in the range (0, 1).
  Default: 0.99.

• filename_prefix (str) – The filename prefix for saved files.
  Default: ‘akmc’.

• write_log (bool) – Whether or not the log that records the outcome for each search should be written.
  Default: True.

• write_searches (bool) – Whether or not the objects for each search (e.g. saddle point, NEB, or MD) should be written.
  Default: True.

• write_markov_chain (bool) – Whether or not the MarkovChain object should be written.
  Default: True.

• write_kmc (bool) – Whether or not the KineticMonteCarlo object should be written.
  Default: True.

• lkmc_parameters (LKMCParameters | None) – The localized AKMC parameters or None if the global AKMC should be used.
  Default: None.

kmc()

Returns:

The KineticMonteCarlo object.

Return type:

KineticMonteCarlo

log()

Returns:

The AKMC log object.

Return type:

AKMCLog

markovChain()

Returns:

The Markov chain object

Return type:

MarkovChain

run(max_searches, max_kmc_steps=None)

Run the adaptive kinetic Monte Carlo simulation.

Parameters:

• max_searches (int) – The AKMC simulation will terminate after this many saddle searches.

• max_kmc_steps (int) – The AKMC simulation will terminate after this many kinetic Monte Carlo steps.
  Default: 10000

state()

Returns:

The current state id.

Return type:

int

Usage Examples

Model the dynamics of two adatoms on a Pt(100) surface.

# -------------------------------------------------------------
# Adaptive Kinetic Monte Carlo
# -------------------------------------------------------------

htst_parameters = HTSTParameters(
    assumed_prefactor=1e14/Second,
    )

# If the the markov chain already exists, read it from file,
# otherwise create a new object.
if os.path.isfile('akmc_markov_chain.nc'):
    markov_chain = nlread('akmc_markov_chain.nc')[0]
else:
    markov_chain = MarkovChain(
        configuration=bulk_configuration,
        configuration_energy=TotalEnergy(bulk_configuration).evaluate(),
        )

# If the the kinetic Monte Carlo already exists, read it from file to resume it,
# otherwise create a new kinetic Monte Carlo object.
if os.path.isfile('akmc_kmc.nc'):
    kmc = nlread('akmc_kmc.nc')[0]
else:
    kmc = None

constraints = [ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12,
               13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25,
               26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
               39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51,
               52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64,
               65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77,
               78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90,
               91, 92, 93, 94, 95, 96, 97]

akmc = AdaptiveKineticMonteCarlo(
    markov_chain=markov_chain,
    kmc_temperature=300.0*Kelvin,
    md_temperature=1000.0*Kelvin,
    calculator=bulk_configuration.calculator(),
    kmc=kmc,
    constraints=constraints,
    confidence=0.99,
    htst_parameters=htst_parameters,
    filename_prefix='akmc',
    )

akmc.run(max_searches=50, max_kmc_steps=10000)

akmc.py

Notes

In solid state systems, atoms reside near their equilibrium positions and reactions are rare events. MD simulations are typically limited to nanosecond timescale while room temperature solid-state reactions of interest may occur on the millisecond timescale or longer.

Adaptive kinetic Monte Carlo (AKMC) is a tool for simulating the dynamics of solid-state systems. [1][2] The goal to is automatically construct a kinetic Monte Carlo (KMC) simulation without prior knowledge of the reaction rates. Unlike lattice KMC, the atomic positions are not limited to lattice sites. Instead the geometries are the stable minima of the potential energy surface obtained from geometry optimization.

The reaction rates are calculated using harmonic transition state theory (HTST). See the HTSTEvent notes for more details. HTST calculates the reaction rate between stable states (minima) if the saddle point between them is known. Thus in AKMC the main goal (and computational expense) is to locate new stable states and the saddle points between then.

Algorithm Details

Saddle points are located using high temperature (controlled by the md_temperature parameter) molecular dynamics (MD) simulations, starting from a geometry-optimized configuration. These trajectories are periodically minimized to determine if a reaction has occurred. A reaction has occurred if the minimized structure is different from the initial configuration.

If a previously discovered minimum is found the saddle search can terminate early, otherwise, the saddle point between the two geometries must be located. This is done in two steps. The first is a NudgedElasticBand (NEB) optimization is run. The converged NEB calculation is then used to make an initial guess at the location of the saddle point. This initial guess is then optimized using minimum-mode following saddle optimization algorithm.

Once a saddle point has been found a check is made to ensure that the saddle point is connected to the initial state. It is possible that the saddle point that was found does not describe a reaction for the initial minima. Two geometry optimization starting from the saddle point (displaced in the position and negative directions of the unstable direction) are performed. If one of the minima is the initial state then the saddle is connected. If neither minima is the starting minima then the saddle is disconnected and the saddle search ends.

Once a connected saddle point is located for the first time, the HTST prefactor must be calculated for the forward and reverse reactions. The prefactor can either be set to a fixed value for all reactions or it can be estimated. The prefactor calculation is controlled by the htst_parameters argument.

The prefactor involves calculating the dynamical matrix at both minima and the saddle point. This is a computationally expensive calculation and for many systems it is unnecessary, because there is a typical value for the prefactor that does not vary much between reactions. This value can be determined by calculating the prefactor for some reactions in the system.

When numerically calculating the prefactor, not all atoms need to be included in the dynamical matrix to get a good result. Instead, only the atoms that move significantly along the reaction pathway (and their nearest neighbors) need to be considered. The minimium_displacement parameter stored in the HTSTParameters object allows you to specify how much an atom needs to move by to be included in the dynamical matrix.

Parallel Saddle Searches

AKMC simulations are highly parallelizable since each saddle search is an independent calculation. Thus, many saddle searches can be run in parallel. When multiple MPI processes are used to run an AKMC calculation, QuantumATK automatically determines how many saddle searches can be run in parallel. It is important to note that one process is always used as a “master” that coordinates the work of the parallel saddle searches. This makes calculating the total number of MPI processes to use a little trickier than normal.

The number of saddle searches that will be run simultaneously is given by:

\[\frac{N_{\rm p} - 1}{N_{\rm s}}\]

where \(N_{\rm p}\) is the total number of MPI Processes and \(N_{\rm s}\) is the number of processes to use per saddle search (this is defined by the ParallelParameters object attached to the calculator that is being used).

Resuming Previous Calculations

The Python script that is generated by the QuantumATK Scripter supports resuming the AKMC simulation. The typical way to run the calculation would be to set the max_searches parameter to a value small enough that the script will finishing running in a reasonable period of time (perhaps a day or so for DFT calculators) and then re-running the script as needed to model the dynamics of interest.

When each saddle search is completed all of the files on disk are updated impenitently. This means that it is possible to monitor the progress of the simulation as it progresses.

[1]

Graeme Henkelman and Hannes Jónsson. Long time scale kinetic monte carlo simulations without lattice approximation and predefined event table. The Journal of Chemical Physics, 115(21):9657–9666, 2001. doi:10.1063/1.1415500.

[2]

Samuel T. Chill and Graeme Henkelman. Molecular dynamics saddle search adaptive kinetic monte carlo. The Journal of Chemical Physics, 140:214110, 2014. doi:10.1063/1.4880721.