User guide
This section offers simple explanations for the most important
parameters defined in the KMC.sh file.
Launching k-ART simulation
Before launching a k-ART simulation, it is necessary to determine the
length of the simulation, either in number of steps (NBRE_KMC_STEPS )
or in total time spent (TOTAL_TIME). The simulation will stop when the
first of these two occurrences is met. Default values for these two
quantities are large and it is recommended to fix at least one of the
two.
It is also necessary to set the TEMPERATURE. It is defined in Kelvin
(considering that internal energy units are in eV). There is no default
value and it must be set.
setenv NBRE_KMC_STEPS 100 # Max number of KMC steps to be executed (def: 50 000)
setenv TOTAL_TIME 1.0 # Maximum simulation time in seconds (def: 20 s)
setenv TEMPERATURE 600.0 # Simulated temperature in kelvin (no def)
setenv ELAPSED_TIME 0.0 # Elapsed time before first KMC step
# (will be added to the sim. time)
Repeating the simulations, specially to find a bug, might be necessary. For this, it is preferable to start with the same seed for the random number generator. It can be fixed with the following command. When not defined, it is extracted from the computer clock.
setenv Random_seed -1103 # Random number generator see (use only when debugging)
Second level of parallelization
k-ART is now able to use two level of parallelization. the first one is about distributing ARTnouveau search throughout a specific group of submaster, and the second level is about computing force throughtout each group of workers associated to a specific submaster. The implementation is done in such a way that one has to define the following variable.
setenv NTRAVAILLEUR k
With k is the number of workers that will compute forces. Be careful of the value that you might want to use because the number of submaster is calculated base on the number of worker and the number of cores asked. If this division is not a integer the code will stop.
Describing the simulated system
The description of the simulated system includes, of course, its size, number of species, box type and the force field used to describe it.
Details of the system: atomic types and active/inactive species
setenv NUMBER_ATOMS 1018 # Total number of atoms
setenv NSPECIES 2 # NUmber of atomic species (default: 2)
setenv ATOMIC_SYMBOLS "Fe C:active Ni:inactive" # Atomic symbols of atoms used in the
# simulation, they must be given in
# the same order than in the potential.
setenv ATOMIC_TYPE_FILE 'Atomic_types' # List of atomic types (read at start of sim)
# (def: Atomic_types)
We can set the atomic symbols using the the string ATOMIC_SYMBOLS. It
is set as a list with the names of each atomic symbol. Also we can say
if each atom type is active or inactive, that is, if to launch for
researching of events on it or not. We can type for example: “Fe” (or
“Fe:active”, although by default all elements are active, so this is
not necessary) and if inactive: “Fe:inactive”. As another example, in
Si-H, the H could be declared inactive such topologies are never
centered on it as the relevant moves involve Si rebonding.
ATOMIC_SYMBOLS” can be set to a string of a Maximum of 30 atomic
symbols.
If the list becomes large we can use the file Atomic_types which also
contains the name of the atomic species as well as whether they are
active or inactive. If neither of these variables are set, default
values are used. If both are used, the Atomic_types is used. Comments
can also added to this file by starting the line with #. The format of
the Atomic_types file is the following (see the Si-H directory in
EXAMPLES):
# Atomic name, state (launch for search of events on it if active)
Si active
H inactive
Maintaining different subgroups of the same element as active or inactive
In some cases, it is useful to select as active only a subset of all
atoms of the same species. To do so, simply give different names to the
various subsets of the same element. For example, with three regions for
a block of silicon. In KMC.sh write:
setenv NSPECIES 3 # The number of different atom types (default: 2)
setenv ATOMIC_SYMBOLS "Si Sit:inactive Sib:inactive"
For lammps, you will need to modify the potential. For example, for
Silicon:
pair_style sw
pair_coeff * * Si.sw Si Si Si
and define the masses
1 28.0855
2 28.0855
3 28.0855
Soft or strict activity
By default, the code uses so-called “soft activity”, i.e. no event are launched from inactive species but events can be recentered around these atoms.
Setting SOFT_ACTIVITY to .false., no event can be centered on
inactive atomic species, even if the initial deformation was made on
an active atom. This, however, can create problems when, for example,
the displacement of the interesting species involves the motion of
inactive species on the same scale. This setting can be made by
inserting the following into KMC.sh:
Note: This option biases the kinetics and must be used carefully.
A system might use both soft and strict activity for different atom
types. For example, if a big systems has many atoms that never will be
moved by events thoughout the whole simulation, these atoms can be
declared strictly inactive to not analyze their topologies and therefore
speed up the simulation’s initialization time. If also certain atoms
won’t yield relevant events, they should be declared softly inactive.
This can be achieved by not using the SOFT_ACTIVITY option and having
an Atomic_types file as follows:
# Atomic name, state (launch for search of events on it if active)
Si active
H inactive
C strictlyinactive
Now, C behaves as if under the option SOFT_ACTIVITY .false., H as if
SOFT_ACTIVITY was .true., and Si is active.
The keyword “strictlyinactive” does not work in ATOMIC_SYMBOLS yet.
This should be implemented in the future and, for better consistency,
the option SOFT_ACTIVITY should probably be removed.
Box size and boundary conditions
Currently, k-ART allows orthorombic as well as triclinic boxes with periodic conditions in all directions. The box size is defined in length units and should be the full box size.
For a cubic box, simply use SIMULATION_BOX. For other boxes, it is
necessary to specify the box length for each direction separately. When
simulating a surface, for example, simply use a very large Z_BOX
setenv SIMULATION_BOX 29.0638147 # Size of the simulation box (x, y and z)
or
setenv X_BOX 28.02 # Simulation box size in the x dimension
setenv Y_BOX 14.17 # Simulation box size in the y dimension
setenv Z_BOX 68.45 # Simulation box size in the z dimension
These parameters may also be defined in the k-ART input configuration file. For a parallelipedic box, only one line is needed after the energy line. The terms represent, respectively, the size of the box in the x, y and z directions.
total energy: -2213.34826631256
21.7200000000000 21.72000000000 21.7200000000000
However, when simulating a triclinic box, three lines are needed. The triclinic parameters matrix is defined following the LAMMPS formalism. The documentation can be found on the
total energy : -2213.34826631256
21.7200000000000 1.000000000000 2.00000000000000
0.00000000000000 21.72000000000 3.00000000000000
0.00000000000000 0.000000000000 21.7200000000000
When using lammps (ENERGY_CALC set to LAM), these parameters are not needed as the box size is already defined in the input of Lammps; so you can safely comment these lines. If they are defined they will be ignored. You may use the following format to input the boxe size to LAMMPS in the LAMMPS input configuration file. For an orthorombic box, the three lines after the number of atomic types line are needed. For a triclinic box, the four lines, including the xy, xz and yz term, are needed.
511 atoms
1 atom types
0.0 21.72 xlo xhi
0.0 21.72 ylo yhi
0.0 21.72 zlo zhi
1.0 0.0 0.0 xy xz yz
Frozen boundary conditions
It is also possible to fix a subset of atoms, i.e., to prevent them from moving and keep their original coordinates. This can be used to impose an external deformation, to look at one surface, fixing the other surface, for example. To allow for most flexibility, the code does not analyse the geometry of the system (surfaces or other) and simply requires a list of atoms that are prevented to move.
To imposed fixed boundaries first set the variable FIXED_BC to TRUE.
By default, the list of fixed atoms must be placed in the file Fixed_atoms. This name can be changed with the envrionnment variable ``FIXED_ATOMS_FILE`:
setenv FIXED_BC .true.
setenv FIXED_ATOMS_FILE new_name.dat # Name of file that contains the list of atoms fixed during the simulation.
# This is optional, one can used "Fixed_atoms' which is the default.
The file must contain the atomic number, one atom by row, followed by free or fixed. More specifically: the file must respect the following format
atomic id (a number) in the first 8 colums
followed by
fixedorfreestarting in column 9
For example, to fix atoms 1 and 3, the file will contain (omit the first line, that shows colum position).
12345678901234567890 # Column number, DO NOT PUT IN THE FILE
1 fixed
2 free
3 fixed
4 free
5 free
Atoms that are not in the list are, by defaults, set to free, so it is sufficient, in general, to only list atoms that are fixed.
Open boundary conditions
The kART code requires periodic boundary conditions. However, to set open boundary condition, it is sufficient to define a simulation box that exceeded the region occupied by the atoms. For example if you want to simulate a cluster with a 10 Å-radius, you can define your box with a size of 50 or 100 Å, to ensure that there will be no interaction between images.
Choice of forcefield/neighbour lists
A few potentials are directly implemented in k-ART
Stillinger-Weber with the original parameters (SWP) or the modified ones by Vink et al. for simulating amorphous silicon (SWA);
The EDIP potential, still for silicon (EDP);
The MEAM potential, implemented by Noël Jackse for Si/Au (MEA)
The EAM potential used for simple metals (EAM).
A mixture of modified Stillinger-Weber potential coupled with a Slater-Buckingham for H-H interactions. This is not great for H diffusion (SIH) (see Ref. [@mousseau90]).
Of those, only Stillinger-Weber is program both with a global and local calculation of forces.
When those forcefields are not sufficient, it is possible to turn to
LAMMPS to access its vast set of forcefields. To do so, one has to
select the LAM option for ENERGY_CALC and to provide the appropriate
LAMMPS input file.
LAMMPS works in a number of different energy and length units. Even though k-ART does not formally care, all default parameters are defined for eV/Å set of units. If the potential is in kcal, it is possible to change it directly at the k-ART level. Other options can be implemented but are not at the moment.
System Potential Command Neighbours
silicon Stillinger-Weber SWP TTL or PAR
Modified SW SWA TTL or PAR
EDIP EDP TTL
metals EAM EAM TTL
Si/metals MEAM MEA TTL
silicon/hydrogen Stillinger-Weber and Slater-Buckingham SIH TTL generic Lammps LAM TTL
pt
setenv ENERGY_CALC SIH # choose between EDP, SWP, SWA, SIH or LAM (Lammps)
setenv FORCE_CALC TTL # choose TTL (total force) or PAR (partial force)
setenv INPUT_LAMMPS_FILE 'in.lammps' # LAMMPS input file when ENERGY_TYPE = LAM
setenv UNITS_CONVERSION metal # Converts the energy and distance units
# to ensure eV and Angstrom in kART
When using an internal potential, it is also necessary to parametrize
the neighboring list. These options are dependent on the type of
forcefield (global or local). Except for Stillinger-Weber, it is
therefore necessary to select : TTL or VER.
setenv UPDATE_VER_NEIB TTL # TTL for total force or PAR for partial force
setenv NEIB_CALC ALL # ALL or VER
setenv UPDATE_TOPO TTL # TTL or PAR
setenv PAR_DISP_THRESH2 0.00001 #max displacement squared which triggers an update of
# neighbor list with VER (default: 0.00001)
Managing output: options and file names
A number of options are possible for outputs, depending both on the amount of information you want and how you want to manage it.
Here is the list along with brief description.
The following allow to follow the run itself. PRINT_DETAILS, for exemple, shows the details of activation and minimization for each attempted event. It is particularly useful when adjusting initial parameters.
If PRINT_DETAILS = .true., then adding MINSAD_DETAILS, force the code to write in a director MINDSAD_DIR all saddle and mininimum files generated. While this later option generates a large number of files, it might be useful for understanding problems with some events.
setenv PRINT_DETAILS .true. # Prints the details of activation and minimization steps in sortieproc (def: false)
setenv MINSAD_DETAILS .false. # Prints the details of activation and minimization (def: false)
The STATISTICS option generates a number of files describing the simulation, including : Energies.dat, that presents information about accepted events, including energies, barriers, rates and displacements; Gen_ev.dat, that gives informations about generic events searches and more at each step, Spec_ev.dat, that focuses on specific events, selec_ev.dat and topo_stat.dat, that provide information about selected events and topologies at each step of the KMC simulation.
setenv STATISTICS .true. # Write statistics about force and event calculation (def: true)
Event catalog - detailed events
The next one present informations about events per se. Generic and specific events, i.e. all details, are saved in a non-formatted file to save space (generic and specific are saved in a different file). It is possible, however, to also write these files in a text form.
setenv USE_TXT_EVENTFILE .true. # If true, writes each generic event in a different file (def: false)
setenv OUTPUT_SPECIFIC .true. # If true, writes each specific configuration in a different file (def: false)
setenv USE_TXT_SPECFILE .true. # If true, writes each specific events in a different file (def: false)
NOTE: THe option USE_TXT_SPECFILE is not yet operational. At the moment, specific files are only write to text format.
It is possible to impose a name for the event files and directories (if USE_TXT_EVENTFILE is TRUE):
setenv UFT_EVENTFILENAME ‘events.uft’ # The name of the uft event file (default: events.uft)
setenv EVENTS_DIR ‘EVENT_DIR’ # The name of the directory where the event txt files are stored (default: EVENT_DIR)
To convert the catalog from one format (txt or utf) to the other, one may use the following commands in the KMC.sh file. When this is done, the code converts the catalog and stops immediately, without any event search.
setenv CONVERT_TO_UFT .false. # IF true, the program will read the txt event files and create a new uft file
setenv CONVERT_TO_TXT .false. # IF true, the program will read the uft file and copy the catalogue to txt files (use for vizualization)
Event catalog - overall properties
The succinct event catalog describing the general properties of events available at each step provide an information unique to k-ART. With it, it is possible to reconstruct the complete local landscape. By default, the output is written in a multiple text files. If a single file is wanted (to reduce the number of files, for example), set SPLIT_EVENT_FILE to FALSE.
setenv OUTPUT_CONFIG_EVENTS .true. # If true, creates an event catalog (all events available at each step - def: true)
setenv SPLIT_EVENT_FILE .false. # If true, puts the event catalog in separate files for each event (def: true)
Note that OUTPUT_CONFIG_EVENTS must be TRUEfor the second option to mean something.
The next option defines the completeness of the output configuraitons. By default, init-saddle-accepted configuration are written along the trajectory. It is also possible to save only minima along the pathway (set ALLCONF_WITH_SADDLEto FALSE):
setenv ALLCONF_WITH_SADDLE .true. # If true, write full configuration events including init-saddle-conf. (def: true)
Other outputs
setenv OUTPUT_NEB_GEN_EVENT .true. # Can be useful
setenv SAVE_FULL_EVENTS .true. # save full events during catalog and after one file per node
Restarting
It is possible to restart and extend a simulation. This simulation can start from the last point, using the catalogue already computed or reconstruct one if there are problems with the one at hand.
It is also possible to start a new simulation using a previously
calculated catalog. In this case, one has only to start with
RESTART_IMPORT set to .true.
By default, all these options are set to .false.
In particularly large or complex systems building the initial catalog can be expensive. Moreover, if the simulation is interrupted during the construction of the catalog kART will have to restart it from nothing. To have the possibility to build the catalog in multiple runs the following parameter in the environment file must be set to true (in the initial run).
setenv SAVE_SEARCH_LIST .true.
This will write the formatted text file search_list.txt. To resume
building the catalog in a subsequent run the following parameter must be
set to true in the environment file:
setenv RESUME_BUILD_CATALOG .true.
Basin parameters
KMC methods are limited in time by the lowest-energy barrier available at any time. When this low-energy leads to a state with a significantly lower energy, it belongs to the kinetics of the systems. However, these low-energy barriers often connect non-diffusive states of similar energy, i.e., states that do not lead to the structural evolution of the system. Such states are called flickers and represent the bane of KMC simulations.
In kinetic ART, we handle both types of events specifically.
Low-energy barrier events
Very low-energy barriers leading to lower-energy states could be handled directly, without going through a KMC choice. Indeed, these barriers typically involve time scales of 10$^{-10}$ to 10$^{-13}$ s and so do not influence the time scale for the problem typically studied with k-ART. There is therefore no impact in applying them directly. The advantage would be that doing so decreases considerably the cost of handling these low barriers, allowing us to focus on the kinetics of interest.
Further improvements: very low-energy barriers.
At the moment, the whole catalog of generic events is constructed before applying a move to the lowest energy barrier (i.e., below the threshold for real KMC). In some disordered systems, it could be preferable to perform the moves as these low barriers are met and to update the catalog only when no new very low-energy barrier are found.
This is yet another modification that we are planning for k-ART in the near future.
Flickers
There are many ways to handle these flickers. In k-ART, we use the autoconstructing basin mean-rate method (bac-MRM), a modified version of the mean-rate method proposed by Puchala et al.. [@puchala]. It is described in some details in Ref. [@beland2011].
setenv OSCILL_TREAT BMRM # choose between BMRM, TABU or NONE
setenv MIN_SIG_BARRIER 0.1 # Max height of dir. and inv. barrier for an event
# to be considered inside a basin
setenv BASIN_LOCAL .false. # If true, local basins are used (default: false)
setenv BASIN_RADIUS 10.0 # If using local basin, radius of a local move
We have also implemented a TABU algorithm. However, for many systems, this algorithm biases significantly the kinetics and we do not really recommend it. A description of the method can be found in Ref. [@vocks2005; @chubynsky2006].
The following environment parameter define its use.
setenv OSCILL_TREAT TABU # choose between BMRM, TABU or NON
setenv FLICKER_DIR 'FLICKERS' # Directory where the TABU event files are stored
setenv EVENTS_MEMORY 100 # Number of stored event TABU (default: 100)
setenv BASIN_MEMORY 10 # Number of stored basins for TABU (default: 10)
Manually adding basin events
In some cases, it is helpful to add events with high energy barriers to
a basin. For example, it can happen that a system flickers between two
states connected by a large $\Delta E$ because any other event has
significantly higher barriers. Treating this transition as a basin
solves this. When using BMRM, it is possible to load a file with the
default name Event_basin and the following structure:
# initial saddle final
1247195 80279 1247195
...
The file name can be changed with
setenv EVENT_BASIN_FILE Event_basin # File that is read for additional basin events
Using this file means that events will be basin events that have energy
barriers below MIN_SIG_BARRIER or topologies from the list or both.
The reverse event to any list entry will automatically be considered.
Basin auto analysis
Because it can be tedious to manually add events to the event basin file, there is the basin auto analysis. Its parameters are
setenv BAA_STEPS 3 # Number of steps analysed (default 0)
setenv BAA_RADIUS 0.3 # Radius for atom to move in, in angstrom (default 0.3)
setenv BAA_MONITORED 100 # Particle to be monitored for the analysis
# (default NUMBER_ATOMS, i.e. the last atom)
The method remembers the last $N=$ BAA_STEPS positions of the atom
with identifier BAA_MONITORED and if all of these positions lie within
a sphere with radius BAA_RADIUS, the event last executed will be added
to the EVENT_BASIN_FILE file. If the parameter BAA_STEPS is not
given or identical to 0, the basin auto analysis is disabled.
Topological identification of local geometries
K-ART relies on the assumption that we can identify the local geometry around each atom in the system by looking at the topology. K-ART works by generating the connectivity graph of a cluster of atoms centered on the target atom. While k-ART has functions to deal when this one-to-one correspondence between geometry and topology breaks down, it is important to choose well the parameters used to construct the connectivity graph. A topology label is associated with each connectivity graph to identify it.
Two file names can be defined for the outfile statistics regarding topologies. These files are not used for restarting.
setenv TOPOLOGY_FILE 'Topologies' # Store info about topologies
setenv TOPO_STAT_FILE 'topos.list' # Store statistics about topologies
Defining graph size: MAX_TOPO_CUTOFF
The three parameters used defining the graph size are:
MAX_TOPO_CUTOFF, TOPO_RADIUS and MAX_NODES_GRAPH.
MAX_TOPO_CUTOFF defines the maximum distance for which 2 atoms are
considered connected. This value is the same for all atom types. A good
value of MAX_TOPO_CUTOFF usually can be found by looking at the radial
distribution function (RDF) of the system. The value should generally be
between the 1st and 2nd peaks. It is important that minute variations in
bond length do not affect the connectedness of two neighboring atoms. In
fact, this is good because the program will be able to identify the same
structures even if the actual bond lengths vary slightly. It is easier
to find this value in crystalline systems since the RDF has well defined
peaks. It is more difficult in amorphous systems, where there is no
clear separation between the peaks. Even so, choosing a value where the
RDF is a minimum between the 1st and 2nd peaks was found to be optimal.
In the case where k-ART cannot recover the geometry from a topology
label, it does what we call splitting the topology. In this case, the
original topology label is called the primary label and we associated
with it secondary labels by changing the cut-off used to connect atoms.
The parameter MIN_TOPO_CUTOFF determines the smallest cut-off value
used. Each time the topology is split, the cut-off will be reduced from
its original MAX_TOPO_CUTOFF by a distance of
(MAX_TOPO_CUTOFF-MIN_TOPO_CUTOFF)/10
and a new topology label is calculated. This is done until a different
topology label is found. The value of MIN_TOPO_CUTOFF should not
ideally be smaller that the minimum bonding distance observed but small
enough so that the connectivity distance changes enough to find new
topology labels.
Improvement for alloys: Multiple MAX_TOPO_CUTOFF
You may also define multiple cutoff distances between atoms for multi-atomic crystals.
In this case, these may be defined directly with both the
MAX_TOPO_CUTOFF and MIN_TOPO_CUTOFF in the environment file:
setenv MAX_TOPO_CUTOFF "2.0 3.2 2.5 2.7 2.9 3.1"
setenv MIN_TOPO_CUTOFF "1.8 2.9 2.1 2.5 2.6 2.7"
The cutoffs can also be placed in a separate file defined by the
environment variable TOPO_CUTOFF_FILE. If this variable exist, it
supersedes the MAX_TOPO_CUTOFF and MIN_TOPO_CUTOFF variable. In the
file, the first line represents the MAX_TOPO_CUTOFF values and the
second the MIN_TOPO_CUTOFF ones:
"2.0 3.2 2.5 2.7 2.9 3.1"
"1.8 2.9 2.1 2.5 2.6 2.7"
For a crystal with two atom types, the format is 1-1, 1-2, 2-2 cutoff lengths. For a crystal with three atom types, as shown in the above examples, the format is 1-1, 1-2, 1-3, 2-2, 2-3, 3-3 cutoff length. The same logic applies for crystals with more atomic types. Note that with both methods the cutoffs lengths must be defined between quotation marks.
Defining graph size: TOPO_RADIUS
The second important parameter is TOPO_RADIUS, which defines the
radius of the sphere centered on the target atom that is used to include
neighboring atoms in the graph. The larger the sphere, the more atoms
are included and the more complex the graph is. A too small value will
result in a graph that does not properly correlates with the local
geometry. A too large value will result in clusters that are too
specific and include atoms situated far away that do not affect the
transition states centered on the target atom. We find that in c-Si
and a-Si, a value of 6 Angstroms works well. Note that as the volume
of this sphere increases, so does the probability of having atoms near
the edge of the sphere. These atoms’ positions can fluctuate in and out
of the surface and change the number of atoms in the graph, resulting a
multiple changes in the topology. Again, the program will work but will
be less efficient because for each new topology, a new series of
(redundant in this case) event searches are launched.
setenv TOPO_RADIUS 6.0 # by default, set to 6 Angst.
setenv MAX_TOPO_CUTOFF 2.8 # maximum length cutoff
setenv MIN_TOPO_CUTOFF 2.2 # Also sets a minimum distance
Further improvements for alloys: MAX_NODES_GRAPH
The cut-off defining the limits of a graph used in the topological classification is not ideal for an alloy. In this case, particularly when the atomic size varies greatly, it would be preferable to use a fixed number of atoms around the center.
This is possible with the command MAX_NODES_GRAPH, which imposes a
maximum number of nodes (atoms) to the graphical analysis. To use this
command, you first must set a cut-off size for TOPO_RADIUS that
corresponds to the correct size of a cluster containing only your
biggest atom (which would have the lowest number of atoms for a given
radius).
You can then use MAX_NODES_GRAPH to impose a maximum number of nodes
(atoms) for the graph.
setenv MAX_NODES_GRAPH 42
To optimize these parameters, it is best to first run kART on a system
with only the largest species and to identify the size of the graph
corresponding to TOPO_RADIUS. Look, for this, into the sortieproc.0,
which gives the cluster size (number of nodes) for the largest atoms.
Fetching topology label for atom: 1( cluster size: 45)
Topology has changed : 0--> 617416
Assigning topoId 617416 to atom 1
Here, for example, setting MAX_NODES_GRAPH to 45 would ensure that
most clusters have a similar size, facilitating classification and
cataloguing
The special case of random alloys: UNIQUE_CRYST_STRUCT
In the case of an alloy whose atoms are randomly distributed on the lattice sites, without any short or long range order (solid solutions for example), the crystallinity should not depend on the atom types but only on the positions of the atoms. The default calculations in kART are not performed this way but one can decide to make the crystallinity check on atomic positions only by setting:
setenv UNIQUE_CRYST_STRUCT = .true.
By doing so, even if the topologies id found by kART do not correspond
to the CRYST_TOPOID given as a parameter, if they are structurally
equivalent they will be considered as crystalline and will be ignored.
Leaving some topologies aside.
K-ART’s performance is proportional to the number of different topologies in the system and not necessarily with the number of atoms. However, to avoid spending time on environments where no event is expected, such as perfectly crystalline regions, it is possible to remove these events from the catalog.
One way to do this is to include the file Topo_ignore in the
simulation directory. This file should contain the topology labels (one
per line) associated with region to be ignored.
The previous solution is applicable to neglect any type of topologies. To avoid spending times on barely non-crystalline environments (i.e. those with a defect at the edge), it is also possible to set, in the environment file, up to five crystalline topologies along with a small cutoff for ignoring them. Either or both of them can be left out if not needed.
The use of many crystalline topologies is useful in the case of an ordered compound, for example.
To get this label, it is necessary to run k-ART on a perfect crystal,
for example, with TOPO_RADIUS set to the same value of
CRYST_TOPO_RADIUS. Run for some seconds, and from sortiproc.0 go to
line “Listing topologies in system”, and get the number (if system is
not perfect crystal, the most repeated number is the CRYST_TOPOID).
IMPORTANT: After this, reset TOPO_RADIUS to its previous value and
check that the crystalline topology is already removed by rerunning.
Then, in sortiproc.0 at the end of “Listing topologies in system”,
look for “topo has no events x”. Verify that the number $x$
corresponds indeed to the crystalline topology (realize that $x$ is
different to the number set in CRYST_TOPOID as TOPO_RADIUS is
reset).
Note that all new topology labels are written in the topos.list and Topologies files. The file topo_stat.dat keeps track of the statistics associated with topology labels per KMC step.
Leaving specific events aside
In some cases it is good idea to ignore some problematic events. There are two ways to do this.
If an undesired events (sometimes an error) has popped in the catalog,
then it might be useful to remove it. This is achieved by creating a
filed called Event_remove that contains the list of all undesired
events. It has the following structure.
# Event id (one per line)
330236
576443
...
576443
where the numbers can be found, for example, in the event_list_xx
files.
Because it removes the events from a previously constructed catalog, it needs the file .utf@ in order to work. If there is no previous catalog, the code will crash.
Now, it is also possible prevent given events to be stored in the catalog. Because the events ids are not uniquely defined, it is then necessary to define the events in terms of their topological identifier at the initial, the saddle and the final states. To allow for further distinction (for example, a rare, but unrealistic events that involves a jump to the second neighbour site but through the same saddle point), it is also necessary to define a minimum displacement threshold. Any event involving the same topologies but with a smaller displacement will not be ignored.
The default file name is Event_ignore and has the following structure
# initial saddle final min_dr
1247195 80279 1247195 0.0
...
For both files, comments can be added to that file using the character
# at the beginning of line.
Activated Event generation
For each topology label found, k-ART launches a series of ART nouveau searches to find all the activated events. The events associated with a topology label are called generic events. K-ART assumes that all the atoms sharing the same topology label have the same generic events. While the actual values of the activated energy and displacements can change slightly atom to atom, they all have them.
The important non-ART nouveau parameters for the generic event
generation are SEARCH_FREQUENCY and USE_LOG_SEARCH. These control
how many event searches are launched per topology. If USE_LOG_SEARCH
is false, then this number is simply SEARCH_FREQUENCY. If
USE_LOG_SEARCH is set to true (default value), then the number of
searches launched are a function of how many times we have seen this
topology before. The number of searches is calculated every turn and is
equal to
$$\label{log_searches} max(SEARCH_FREQ \times (1 + \log_{10} (COUNTS)) - ATTEMPTS, 0)
where $COUNTS$ is the number of times the topology has been seen and $ATTEMPTS$ is the previous number of event searches launched. This logarithmic search function is used so that topologies seen more often are searched more. Moreover, topologies are also searched multiple times during a simulation (and not just when first found), albeit less and less often.
All generic events found are written to a binary file called
events.uft by default. This file can be reused to restart a simulation
or to start a new one. The topology file topos.list also needs to be
provided. If the parameter USE_TXT_EVENTFILE is set to true, a folder
called EVENT_DIR will also be created and a text version of the event
files will be created during the simulation. Note that in the TOOLS
directory, the program kART_uft2txt.f90 is provided. This stand-alone
program can convert an events.uft file to text format for
visualization purposes (see the READ_ME.txt file).
Generic events
Generic events are stored in the catalog. They are used to construct, at every step, the specific events that are put in the KMC tree.
Every time a new topology is found, a search is launched to associate
generic events with this topology. The number of searches is defined by
SEARCH_FREQUENCY.
As a topology is found more often, new searches are made to minimize the
probability of having missed an important barrier. This is done on a
logarithmic scale, i.e., the number of addition searches is given by the
logarithm of a topology’s occurrence multiplied by SEARCH_FREQUENCY.
Specific events
When arriving in a new minimum, k-ART first updates its list of topologies and completes the list of generic events.
But default all generic events contributing to 99.99 % (1 in 10 000) of the rate are reconstructed for specific configurations, with the saddle points fully relaxed. This ensures that elastic and geometric effects are fully taken into account
It is possible to modify this proportion with the SPECIFIC_THRESHOLD paramter.
This one is set, by default, as :
setenv SPECIFIC_THRESHOLD 99.99
To reconstruct as specific only events having a propability of 1 in 1000, for example, enter:
setenv SPECIFIC_THRESHOLD 99,9
A number of environment variables relate to the reconstruction and relaxation of specific events. In general, the default values can be set.
However, for topologies with high symmetries (e.g. relaxed interstitial
in a crystalline environment), there can be problems with the mapping of
the events onto the system. For $n$-fold symmetry, the chances of
correct mapping are $1/n$ . If in the sortieproc.0 file there is “Mapping
… failed”, it can help to increase the parameter REFINE_ATTEMPTS to
for example $N = 50$, reducing the overall probability of mapping
failure $(1-1/n)^N \ll 1$. Increasing REFINE_ATTEMPTS comes at almost
no cost.
setenv REFINE_ATTEMPTS 2 # Default value for this parameter
Event analysis
Using local forces for accelerating simulations of large systems
Because events are intrinsically local, it is not necessary to compute the force on all atoms to generate generic events when the box is sufficiently large, since these will be further refined, taking full account of the complete system, when relevant.
In this case, kART can therefore compute forces only on a limited subset of atoms, independently of system size, which can considerably accelerated the generation of a catalog and even the relaxation of specific events as described below.
In the next subsection, we describe the parameters associated with this option. In the following one, its implementation.
Selecting parameters
The option LOCAL_FORCE determines the use of local force. By default,
this option is turned down.
When this option is , the global simulation box is cut into cells of
width defined by CELL_LENGTH. The same length is used in all
directions, irrespective of the box size. In general, as explained below
this length should be of the order or larger than the potential cut-off.
The second option, REGION_WIDTH, defines the number of cells included
in the local region on each side of the central cell. For example, if
REGION_WIDTH is set to 2, the total length of the local region for
force calculations is 5 cells (2 on each side plus the central cell).
This allow the selection of a spherical region around the central atom defined with a first inner region with a radius of ( 2*REGION_WIDTH -1)/2 followed by an outer shell of width REGION_WIDTH to (2*REGION_WIDTH+1)/2.
Since, the local region is set-up with open boundary conditions, as
shown in Fig. [[fig:local]] fig:local{reference-type=”ref”
reference=”fig:local”}, forces on the atoms in the outer shell cannot be
used. This is why, only atoms in the inner sphere are considered for
activation and relaxation. This explains why it is recommended that the
CELL_LENGTH be at least as long as the potential cut-off (or, at
least, sufficiently large to ensure that the local-region’s boundary
does not affect the results.
If REGION_WIDTH is set to 2, as per default, this means that the inner
core, where atoms can move with local forces, is 3 cells large. When
elastic deformation is more long range, it might be necessary to
increase the region width.
Adapt these parameters if the saddle points found with local force are systematically lost during the final global force convergence.
To check that what these parameters do, the outputs show the number of
atoms in the local region : nat_inner, nat_outer and
nat_local = nat_inner + nat_outer are printed every time a local
region is defined.
Similarly, the convergence to saddle point for specific events indicates whether local or global forces are used, allowing a check on the validity of the local region.
If there is little change between saddle points found after local and
global convergence, or if you do not need high accuracy on the saddle
point, you might disable global convergence at the saddle (global
convergence is still used at the final minimum). In this case, set
GLOBAL_CONVERGENCE to .false.@
[[fig:local]]{#fig:local label=”fig:local”}
For large box, it might be difficult to converge the eigenvalue and
eigenvector surrouding a local deformation because of the way the
initial random vector in the Lanczos calculation is set up. To
circumvent this problem it might be useful to set LOCAL_LANCOS to
true. In this case, global forces are systematically used, but the
initial deformation for construction the Lanczos matrix is done on a
local region defined by the same parameters ass for LOCAL_FORCE.
Implementation
When set-up, local forces are used for finding all generic events. In LAMMPS, the local region is positionned in the same box as for the full system and atoms are not recentered. Periodic boundary conditions determined by the large box are used to ensure the continuity of the local region, allowing the algorithm to work even when there are surfaces and interfaces in the original box.
Convergence of specific events is also first done with local forces, once the local force is below threshold, convergence is relaunched with global forces, using the eigenvector obtained with local forces as seed.
Breaking the kinetics : approach to get specific results
Kinetic ART implements a number of approach that break the correct kinetics but allow to generate specific desired events. While, these cases, the kinetics and the sequence of events may be wrong, it can provide useful information difficult to obtain otherwise.
The following options are possible:
Active/inactive species
A number of options are available to set some active and inactive species as described in detail here: https://kart-doc.readthedocs.io/en/latest/user_guide.html#details-of-the-system-atomic-types-and-active-inactive-species
Since this option is a very useful way to focus on the right events without too much bias, it is the most used one of this section.
Adding additional conditions for event selection
It is possible to program in specific conditions on events to be added to the catalog. This is useful, for example, when one is interested in specific type of mechanisms. For example, if one want to impose diffusion from the surface to the bulk, it is possible to assess, after each event is generated, whether this event is compatible with insertion, something that might not be possible to impose by the initial displacement or even looking at the saddle.
Setting ADDITIONAL_CONDITIONS to TRUE will force the code to call the additional_conditions subroutine (to be found in src/lib) after an event generated. It is then possible to perform additional tests to decide whether or not the event will be added to the catalog or not.
The additional_conditions.f90 file shows an example for an atom diffusing in the bulk. The first lines of the routines define positions and energies at the initial, saddle and final states for use on the various conditions. To prevent an event from being added to the catalog, set:
success = .false.
Otherwise, if the event should be kepts, then set success = .true.
Remove detailed balance
By default, all events are checked for connectivity. From the saddle point, the configuration is pushed forward and backward to converge both on the new (final) energy minimum and the original (initial) one.
Two parameters are important here:
setenv CHECK_INI_SAD_CONNECTIVITY .true. # Impose a check on the connectivity (detailed balance; def: .true.)
setenv SADDLE_PUSH_PARAM 0.25 # Fraction of displacement for push back to initial point
CHECK_INI_SAD_CONNECTIVITY is a logical parameter than enforce the construction of a connected path between initial-saddle-final states. If this parameter is TRUE, then any event that is not fully reversible is not added to the catalogue.
To reconstruct the pathway, the system is first pushed away from the saddle point along the direction of negative curvature. SADDLE_PUSH_PARAM determines the size of this push (typically betwen a fraction between 0.1 and 0.25 (10 to 25 %) of the displacement from the initial minimum to the saddle). If the reconstruction fails, this fraction is decreased automatically and a new attemp is made by the code to converge to the initial minimum.
When not to impose detailed balance?
For some very rough landscape, it might be appropriate to set the CHECK_INI_SAD_CONNECTIVITY parameter to FALSE. In general, though, this is not recommended.
Bias the initial displacement
By default, the initial random displacement is selected from a complete hypersphere. It is also possible to bias the displacement by forcing the initial displacement to be in the negative z direction to favor motion into a surface which is normal to +z. This option is only applied when the central atom belows to an active species. This is down with :
setenv INITIAL_PUSH_DOWNWARDS .true. # Push in the negative z direction
# (xy anything) (default: .false.)
This should be used with caution as it will then apply to all atom types in the box.
Directed pathway
Finally, kART provides an option for constructed a pathway step by step. Setting the DIRECTED_PATHWAY to TRUE allows the use to choose the next step among events available in the catalog. This means that once the event catalog for a configuration is completed, the code will stop and wait for the user to indicate what event to select.
To do that first set
setenv DIRECTED_PATHWAY .true.
in KMC.sh
Then launch the simulation.
One the event catalog for the current state is completed, The code stops, providing the following instructions:
KMC step no: 1 ( 1 for this run)
**************************************************************************
Waiting for instructions regarding the next event, event number : 1
The file, named, 'imposed_event', must have the following format:
event_number AtomId IniTopoId eventId Spec_id
The code then waits for a few seconds, attempting to read the file. If now file is created, then the code stops.
To move to the next step, a file imposed_event much be created. With the information taken from either Event_catalogue.dator from an event catalog in EVLIST_DIR.
For example, to select the 4th event from
#TypeId AtomId IniTopoId SadTopoId FinTopoId eventId Spec_id barrier inv_bar asy_ener true_rate inisad_dr inifin_dr refined
#****** ******* ******* ******* ******* ******* ******* ******* ********* ********* ********* ********* ********* *******
1 331 236619 1253516 1108930 421298 1 2.351033 0.084875 2.266158 2.007E-11 1.610530 1.972890 F
1 331 236619 624389 296339 1223290 2 2.854080 0.495831 2.358248 1.707E-16 1.627572 2.170290 F
1 357 236619 1253516 1108930 421298 1 2.351033 0.084875 2.266158 2.007E-11 1.610530 1.972890 F
1 357 236619 624389 296339 1223290 2 2.854080 0.495831 2.358248 1.707E-16 1.627572 2.170290 F
1 363 236619 1253516 1108930 421298 1 2.351033 0.084875 2.266158 2.007E-11 1.610530 1.972890 F
1 363 236619 624389 296339 1223290 2 2.854080 0.495831 2.358248 1.707E-16 1.627572 2.170290 F
1 365 236619 1253516 1108930 421298 1 2.351033 0.084875 2.266158 2.007E-11 1.610530 1.972890 F
We create a file imposed_event with the line:
1 357 236619 1223290 2
Then relaunch the code, making sure to set : setenv RESTART_KMC .true.