AdaptiveKineticMonteCarlo

class AdaptiveKineticMonteCarlo(markov_chain, kmc_temperature, calculator, kmc=None, superbasin_threshold=None, saddle_offset_magnitude=None, saddle_search_method=None, optimize_geometry_parameters=None, htst_parameters=None, confidence=None, filename_prefix=None, write_log=None, write_searches=None, write_markov_chain=None, write_kmc=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.

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

• 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_offset_magnitude (SaddleSearchParameters) – The parameters for the saddle search.

• saddle_search_method (LanczosSaddleSearch | ARTnSaddleSearch | MDSaddleSearch) –

• optimize_geometry_parameters (OptimizeGeometryParameters) –

• 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 the log that records the outcome for each search should be written.
  Default: True.

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

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

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

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

Below are two examples of modelling diffusion using AKMC. The first example models the diffusion of Hydrogen in a Copper lattice using the Lanczos saddle search method to locate new states. In the second example, two platinum atoms are diffusing on a platinum surface This example uses the MD saddle search method to locate new states.

# Set up the direction generator for the saddle search.
direction_generator = HypersphereDirection(
    order=Furthest,
    magnitude=(0.5*Angstrom, 1.0*Angstrom),
)

# Set up the LanczosSaddleSearch object. We are interested in the hydrogen atom, so we set the
# index_selector to the index of the hydrogen atom.
hydrogen_index = copper.indicesFromTags(['hydrogen'])
saddle_search = LanczosSaddleSearch(
    initial_direction_generator=direction_generator,
    constraints=optimized_configuration.indicesFromTags(['constrain']),
    index_selector=hydrogen_index,
)

akmc = AdaptiveKineticMonteCarlo(
    markov_chain=markov_chain,
    kmc_temperature=300.0 * Kelvin,
    calculator=optimized_configuration.calculator(),
    kmc=kmc,
    saddle_search_method=saddle_search
)

akmc.run(max_searches=10, max_kmc_steps=1000)

hydrogen_diffusion_example.py

saddle_search = MDSaddleSearch(
    temperature=1000.0*Kelvin,
    constraints=constraints,
)

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

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

There are three different saddle search methods available. The LanczosSaddleSearch and ARTnSaddleSearch methods are based on the Lanczos and ARTn algorithms respectively. The MDSaddleSearch method is based on molecular dynamics simulations. The saddle search methods are used to locate saddle points on the potential energy surface (PES) between two given geometries. If no saddle point is found, a number of different exceptions can be raised, depending on why the search failed.

We recommend using the LanczosSaddleSearch method for most systems (which is also the default). As an alternative, the ARTnSaddleSearch method can be used.

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.

Note that there is a slight difference in the meaning of the max_searches argument depending on which saddle search method used. When using MDSaddleSearch, saddle searches are performed each KMC step until either the specified confidence has been reached, or the maximum number of searches has been performed. When using the other two saddle search methods the saddle searches are terminated when max_searches have been performed, since in this case there is no concept of confidence.

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.

AKMC With MDSaddleSearch and Active Learning

When running AKMC calculation with the MDSaddleSearch method, one can also utilize active learning (see ActiveLearningSimulation) to simulate the dynamics of solid-state systems while training MTPs. While AKMC with active learning can be utilized for different purposes, the main purpose is to populate a training set for an ActiveLearningSimulation object with diverse configurations encountered during AKMC molecular dynamics or optimization steps. AKMC with active learning can also be used for production runs but, for this use case, the quality of the results depends on the system at hand as well as the variance and degree of coverage of the initial active learning training set in respect to the actual system (as described by the reference calculator).

As with normal AKMC usage, the initial configuration needs to be preoptimized and have its TotalEnergy evaluated before starting the AKMC state searches. When passing an ActiveLearningSimulation object to the AKMC object, the AKMC simulation will automatically use the active learning object’s methods for MD and optimization. Hence, the calculator input is ignored and will not be used during the AKMC simulation. Configurations and their associated energies are therefore evaluated via the MTP calculator as defined from the current active learning object/training set state. Here it should be noted that the AKMC simulation does not restart in the event that new configurations are added to the training data as is the default for other active learning applications. Instead, the AKMC simulation continues from the latest transition state search but with an updated MTP calculator.

Since the fitted MTP can change over the course of the AKMC simulation, certain pitfalls should be noted - especially if used for production runs. If the approach is mainly used for diversifying MTP training data, the implications on the quality of the AKMC output are not as severe. States found at previous cycles in the AKMC simulation and their corresponding energies might no longer be accurate according to the newer MTP version. Assuming a simulation is started with an ActiveLearningSimulation object that was created from a small TrainingSet or a training set that does not represent the initially encountered solid state dynamics well, the credibility of the early AKMC results for energies and states is limited. Nevertheless, as the simulation carries on, more and more structures will be added to the training set from other phase-space regions of the solid assuming that the given reference calculator can describe the structures with sufficient precision and level of detail. If the conditions for a such scenario is met in a production run, it should be considered to either use a larger and more diverse initial training set, or starting a new, separate AKMC simulation using the updated MTP from the beginning. Another significant caveat is that the saddle searches cannot utilize the parallelization scheme described above. In order to ensure consistency between subsequent state transition searches, the AKMC simulation is only allowed to run saddle searches in serial mode when using active learning. In spite of these caveats, AKMC with active learning is a promising approach for speeding up simulations of solid state dynamics in complex systems.

A usage example showing how the method could be utilized is sketched in the following code snippet, where a silicon carbide supercell with an interstitial magnesium atom is modelled. A training set for the active learning object has been created beforehand. The code snippet shows how the setup deviates from a normal AKMC simulation in order to have the initial configuration have a calculator and total energy based on the MTP created from the ActiveLearningSimulation object.

# First AKMC simulation: Initial configuration setup for an AKMC simulation
# with active learning that should be started from scratch

# Set up configuration
configuration_0 = BulkConfiguration(
    bravais_lattice=lattice,
    elements=elements,
    fractional_coordinates=fractional_coordinates
    )

# Add tags - constraints needed for AKMC
configuration_0.addTags('fixed', [36, 37])

configuration_name_0 = "optimized_configuration_sic_mginterstitial"

configuration_table = Table()
configuration_table.addInstanceColumn(key='configuration', types=BulkConfiguration)
configuration_table.append(configuration_0)

configuration_names_table = Table()
configuration_names_table.addStringColumn(key='configuration_name')
configuration_names_table.append(configuration_name_0)


# %% ActiveLearningOptimizeGeometry

optimized_configuration_table = runActiveLearningOptimizeGeometry(
    active_learning_simulation=active_learning_simulation,
    configurations=configuration_table,
    max_forces=0.01 * eV / Angstrom,
    constraints=[BravaisLatticeConstraint()],
    trajectory_filename='AKMC-AL_SiC_Mg-interstitial.hdf5',
    trajectory_object_id='alopt',
    log_filename_prefix='AKMC-AL_pre-optimize_geometry',
)

nlsave(
    'AKMC-AL_SiC_Mg-interstitial.hdf5', optimized_configuration_table, object_id='alopt_1'
)

training_set_table = active_learning_simulation.trainingSetTable()


# Custom block for extraction the configuration from the table


def custom(table):
    # This custom script supports atkpython syntax
    # and can perform almost any procedure.

    configuration = table[0, 0]
    configuration.update()
    return configuration


configuration = custom(optimized_configuration_table)

nlsave('AKMC-AL_SiC_Mg-interstitial.hdf5', configuration)


# %% TotalEnergy

total_energy = TotalEnergy(configuration=configuration)
nlsave('AKMC-AL_SiC_Mg-interstitial.hdf5', total_energy)

After the initial setup, the AKMC simulation is scripted as usual but with an extra parameter for the active learning simulation object.

saddle_search_method = MDSaddleSearch(
    temperature=2900.0 * Kelvin,
    saddle_search_parameters=saddle_search_parameters,
    active_learning=active_learning_simulation,
    constraints=configuration.indicesFromTags(['fixed']),
)

akmc = AdaptiveKineticMonteCarlo(
    markov_chain=markov_chain,
    kmc_temperature=300.0 * Kelvin,
    calculator=configuration.calculator(),
    kmc=kmc,
    filename_prefix='akmc-al',
    saddle_search_method=saddle_search_method,
)

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

The full script is available for download: AKMC-AL_example.py

The restarting with AKMC-AL needs to be handled with more care than for regular AKMC where a static calculator is used. Due to the (potential) updates of the MTP, there is no guarantee that one AKMC simulation script run is going to produce results that are consistent with the next run. For production runs, if many new candidates are added during a script run, it may be favorable to use the most recent ActiveLearningSimulation object to create a new AKMC simulation (with new Markov Chain and Kinetic Monte Carlo filenames). If the MTP on the other hand is already so mature that it is rarely updated or if training set updates do not alter predictions significantly, the AKMC simulation can be resumed by loading the saved initial configuration and total energy objects and use them as AKMC inputs instead of preoptimizing the initial configuration and evaluating its total energy again. The restarting syntax is also outlined in the script, but is commented out and should remain so until the first submission of the script has finished running.

[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.