#!/usr/bin/env python
import math
import sys
import configparser
import shutil
from time import sleep
import time as unix_time
import optparse
import logging
import logging.handlers
logger = logging.getLogger('akmc')
import numpy as np
from pathlib import Path
np.seterr(divide="raise", over="raise", under="print", invalid="raise")
from eon import version
from eon.config import ConfigClass
from eon import communicator
from eon import locking
from eon import akmcstatelist
from eon import explorer
from eon import fileio as io
from eon import atoms
from eon import superbasinscheme
from eon import askmc
from eon import movie
[docs]
def akmc(config: ConfigClass = None, steps=0):
if config is None:
raise TypeError("akmc requires a ConfigClass instance")
"""Poll for status of AKMC clients and possibly make KMC steps.
Returns the number of KMC steps in the current run at the end of
this function. The argument "steps" gives the number of KMC steps
done before calling this function, defaulting to zero.
"""
#log version information
logger.info('Eon version %s', version)
# Here's what this does:
# 1) Read in the state of our calculation from last time
# 2) Initialize necessary data structures (statelist, communicator, displace)
# 3) Get any results that have come in
# 4) Possibly take a KMC step
# 5) Make new work units
# 6) Write out the state of the simulation
# First of all, does the root directory even exist?
if not Path(config.path_root).is_dir():
raise FileNotFoundError(
f"root directory does not exist: {config.path_root}"
)
# If we are saving debug results, create the directory if it does not exist.
if config.debug_keep_all_results:
rp = Path(config.path_root) / config.debug_results_path
rp.mkdir(exist_ok=True)
# Define constants.
kT = config.main_temperature/11604.5 #in eV
# Load metadata, the state list, and the current state.
start_state_num, time, previous_state_num, first_run, previous_temperature = get_akmc_metadata(config)
states = get_statelist(kT, config=config)
current_state = states.get_state(start_state_num)
# --- START: DYNAMIC ATOM LIST SCRIPT EXECUTION ---
if config.displace_atom_kmc_state_script:
current_state.get_displacement_atom_list(config)
# --- END: DYNAMIC ATOM LIST SCRIPT EXECUTION ---
if first_run:
previous_temperature = config.main_temperature
# Did the temperature change? Then the existing superbasins are invalid now.
if abs(previous_temperature - config.main_temperature) > 1e-6:
# Remove superbasin data.
if Path(config.sb_path).is_dir():
shutil.rmtree(config.sb_path)
for state_dir in Path(config.path_states).iterdir():
if state_dir.name == 'state_table':
continue
superbasin_file = state_dir / config.sb_state_file
if superbasin_file.is_file():
superbasin_file.unlink()
# Keep the new temperature.
previous_temperature = config.main_temperature
if previous_state_num == -1:
previous_state = current_state
else:
previous_state = states.get_state(previous_state_num)
# If the Novotny-based superbasining scheme is being used, initialize it.
if config.sb_on:
superbasining = get_superbasin_scheme(states, config)
sb = superbasining.get_containing_superbasin(current_state)
# If we are exploring states in a superbasin, we will always
# explore the state with the lowest confidence.
if sb:
explore_state = sb.get_lowest_confidence_state()
else:
explore_state = current_state
else:
superbasining = None
sb = None
explore_state = current_state
state_explorer = explorer.get_minmodexplorer(config)(
states, previous_state, explore_state, superbasin=sb, config=config)
state_explorer.explore()
# Take a KMC step, if it's time.
current_state, previous_state, time, steps = kmc_step(current_state, states, time, kT, superbasining, steps, config=config)
parser = configparser.ConfigParser()
if previous_state.number != current_state.number:
previous_state_num = previous_state.number
write_akmc_metadata(parser, current_state.number, time, previous_state_num, previous_temperature)
io.write_info_txt(config, parser)
io.save_prng_state(io.prng_state_path(config))
return steps
[docs]
def get_statelist(kT, config: ConfigClass = None):
if config is None:
raise TypeError("get_statelist requires a ConfigClass instance")
initial_state_path = str(Path(config.path_root) / "pos.con")
return akmcstatelist.AKMCStateList(
kT,
config.akmc_thermal_window,
config.akmc_max_thermal_window,
initial_state_path,
filter_hole=config.disp_moved_only,
config=config,
)
def _resolve_akmc_state(states, number):
"""Return a state by number, or None when the list has no such state."""
getter = getattr(states, "get_state", None)
if getter is None:
return None
try:
return getter(int(number))
except Exception:
return None
[docs]
def get_superbasin_scheme(states, config):
kT = config.main_temperature / 11604.5
if config.sb_scheme == 'transition_counting':
superbasining = superbasinscheme.TransitionCounting(
config.sb_path, states, kT, config.sb_tc_ntrans, config=config)
elif config.sb_scheme == 'energy_level':
superbasining = superbasinscheme.EnergyLevel(
config.sb_path, states, kT, config.sb_el_energy_increment, config=config)
elif config.sb_scheme == 'rate':
superbasining = superbasinscheme.RateThreshold(
config.sb_path, states, kT, config.sb_rt_rate_threshold, config=config)
else:
raise ValueError(f"Unknown superbasin scheme: {config.sb_scheme}")
return superbasining
[docs]
def kmc_step(current_state, states, time, kT, superbasining, steps=0, config: ConfigClass = None):
if config is None:
raise TypeError("kmc_step requires a ConfigClass instance")
t1 = unix_time.time()
previous_state = current_state
# If the Chatterjee & Voter superbasin acceleration method is being used
if config.askmc_on:
pass_rec_path = None
asKMC = askmc.ASKMC(kT, states, config.askmc_confidence, config.askmc_alpha,
config.askmc_gamma, config.askmc_barrier_test_on,
config.askmc_connections_test_on, config.sb_recycling_on,
config.path_root, config.akmc_thermal_window,
recycle_path = pass_rec_path)
# The system might be in a superbasin
if config.sb_on:
sb = superbasining.get_containing_superbasin(current_state)
else:
sb = None
# discover_decide does not wait on the repeat-count confidence.
# One exit is taken and the product is left for the outer search.
# A state that already meets the threshold keeps the multi-step loop.
discover_decide = bool(getattr(config, "amsel_discover_decide", False))
unconfident_exit_taken = False
while steps < config.akmc_max_kmc_steps or config.akmc_max_kmc_steps == 0:
if unconfident_exit_taken:
break
confident = (
(not sb and current_state.get_confidence() >= config.akmc_confidence)
or (sb and sb.get_confidence() >= config.akmc_confidence)
)
if not confident and not discover_decide:
break
used_superbasin = False
amsel_exit = None
# discover_decide runs from [amsel] alone. use_mcamc is not required.
if discover_decide:
from eon.amsel_superbasin_gate import amsel_discover_exit
amsel_exit = amsel_discover_exit(
current_state,
lambda number: _resolve_akmc_state(states, number),
config,
uniform=np.random.random_sample,
)
if not confident:
if amsel_exit is None:
break
unconfident_exit_taken = True
# Do a KMC step.
steps += 1
if amsel_exit is not None:
mean_time = amsel_exit.mean_time
exit_state = _resolve_akmc_state(states, amsel_exit.exit_number)
if exit_state is None:
exit_state = current_state
current_state = exit_state
next_state = states.get_product_state(
amsel_exit.exit_number, amsel_exit.proc_id
)
elif config.sb_on and sb:
from eon.amsel_superbasin_gate import AmselSuperbasinReject
try:
mean_time, current_state, next_state, sb_proc_id_out, sb_id = sb.step(
current_state, states.get_product_state
)
used_superbasin = True
except AmselSuperbasinReject as amsel_exc:
# Fall back to ordinary single-state KMC for this step.
logger.info(
"amsel rejected superbasin; falling back to single-state KMC: %s",
amsel_exc,
)
sb = None
if amsel_exit is None and not used_superbasin:
if config.askmc_on:
rate_table = asKMC.get_ratetable(current_state)
else:
rate_table = current_state.get_ratetable()
if len(rate_table) == 0:
logger.error("No processes in rate table, but confidence "
"has been reached")
ratesum = sum((row[1] for row in rate_table), 0.0)
if len(rate_table) == 0 or ratesum == 0.0:
# The confidence gate already counted this attempt.
steps -= 1
break
u = np.random.random_sample()
p = 0.0
nsid = 1.1 # Next state process id, will throw exception if remains unchanged.
# A boolean false and the historical string False both mean off.
# The string compare alone treats the boolean default as on.
target = config.debug_target_trajectory
if target and str(target) != "False":
# Get the Dynamics objects.
owndynamics = io.Dynamics(str(Path(config.path_results) / "dynamics.txt")).get()
targetdynamics = io.Dynamics(str(Path(config.debug_target_trajectory) / "dynamics.txt")).get()
# Get the current_step.
try:
current_step = len(owndynamics)
except:
current_step = 0
# Get the target step process id.
if current_step > 0:
stateid = targetdynamics[current_step]['reactant']
else:
stateid = 0
try:
procid = targetdynamics[current_step]['process']
except (IndexError, KeyError, TypeError) as exc:
raise RuntimeError(
"Can no longer follow target trajectory"
) from exc
# Load the con file for that process saddle.
procdata = Path(config.debug_target_trajectory) / "states" / str(stateid) / "procdata"
targetSaddleCon = io.loadcon(str(procdata / ("saddle_%d.con" % procid)))
targetProductCon = io.loadcon(str(procdata / ("product_%d.con" % procid)))
ibox = np.linalg.inv(targetSaddleCon.box)
# See if we have this process
for i in range(len(rate_table)):
p1 = current_state.get_process_saddle(rate_table[i][0])
for dist in atoms.per_atom_norm_gen(p1.free_r() - targetSaddleCon.free_r(), targetSaddleCon.box, ibox):
if dist > config.comp_eps_r:
break
else:
p1 = current_state.get_process_product(rate_table[i][0])
for dist in atoms.per_atom_norm_gen(p1.free_r() - targetProductCon.free_r(), targetProductCon.box, ibox):
if dist > config.comp_eps_r:
break
else:
nsid = i
break
else:
raise RuntimeError("Can no longer follow target trajectory")
# We are not following another trajectory:
else:
for i, row in enumerate(rate_table):
p += row[1]/ratesum
if p>u:
nsid = i
break
else:
logger.warning("Warning: Failed to select rate; p = " + str(p))
break
next_state = states.get_product_state(current_state.number, rate_table[nsid][0])
mean_time = 1.0/ratesum
# Accounting for time. An FPTA exit is already one sample.
if amsel_exit is not None and amsel_exit.time_is_sample:
step_time = mean_time
elif config.debug_use_mean_time:
step_time = mean_time
else:
#np.random.random_sample() uses [0,1)
#which could produce issues with math.log()
step_time = -mean_time*math.log(1 - np.random.random_sample())
time += step_time
# Pass transition information to extension schemes
if config.askmc_on:
asKMC.register_transition(current_state, next_state)
if config.sb_on:
superbasining.register_transition(current_state, next_state)
if amsel_exit is not None:
proc_id_out = amsel_exit.proc_id
elif config.sb_on and sb:
proc_id_out = -1
else:
proc_id_out = rate_table[nsid][0]
# Write data to disk
dynamics = io.Dynamics(str(Path(config.path_results) / "dynamics.txt"))
if proc_id_out != -1:
proc = current_state.get_process(proc_id_out)
dynamics.append(current_state.number, proc_id_out, next_state.number, step_time, time, proc['barrier'], proc['rate'], current_state.get_energy())
logger.info("KMC step from state %i through process %i to state %i ", current_state.number, proc_id_out, next_state.number)
else:
# Superbasin hop: process id is not a single-state table row.
dynamics.append_sb(current_state.number, sb_proc_id_out, next_state.number, step_time, time, sb_id, 1.0/mean_time, current_state.get_energy())
logger.info("SB step from state %i through process %i to state %i ", current_state.number, sb_proc_id_out, next_state.number)
#criterion used to stop the job: currently an energy limit is used as criterion
if current_state.get_energy() > config.debug_stop_criterion:
raise RuntimeError(
f"energy {current_state.get_energy()} exceeds the stop criterion"
)
previous_state = current_state
current_state = next_state
# The system might be in a superbasin
if config.sb_on:
sb = superbasining.get_containing_superbasin(current_state)
else:
sb = None
if config.sb_on:
superbasining.write_data()
if not sb:
logger.info("Currently in state %i with confidence %.6f", current_state.number,
current_state.get_confidence())
else:
logger.info("Currently in state %i (superbasin %i) with confidence %.6f",
current_state.number,
sb.id,
sb.get_confidence())
t2 = unix_time.time()
logger.debug("KMC finished in %.4f seconds", (t2-t1))
if t2 - t1 > 0:
logger.debug("%.2f KMC steps per second", float(steps)/(t2-t1))
return current_state, previous_state, time, steps
[docs]
def main(config: ConfigClass = None):
if config is None:
config = ConfigClass()
optpar = optparse.OptionParser(usage = "usage: %prog [options] config.ini")
optpar.add_option("-C", "--continuous", action="store_true", dest="continuous", default=False, help="don't quit")
optpar.add_option("-R", "--reset", action="store_true", dest="reset", default = False, help="reset the aKMC simulation, discarding all data")
optpar.add_option("-f", "--force", action="store_true", dest="force", default = False, help="force a reset, no questions asked")
optpar.add_option("-r", "--restart", action="store_true", dest="restart", default = False, help="restart the aKMC simulations from a clean dynamics.txt file")
optpar.add_option("-s", "--status", action="store_true", dest="print_status", default = False, help = "print the status of the simulation and currently running jobs")
optpar.add_option("-q", "--quiet", action="store_true", dest="quiet", default=False,help="only write to the log file")
optpar.add_option("-m", "--movie", action="store", dest="movie_type", default = "", help="Specify the type of movie to make [dynamics, states, fastestpath, fastestfullpath, graph, processes]. Process movies are specified like so: --movie processes,statenumber,processlimit. Where processes is the string processes, statenumber is the number of the state that you want to view, and process limit is the maximum number of processes you would like in the movie. The returned processes are reverse sorted by rate such that the fastest processes is the first in the movie.")
optpar.add_option("-M", "--separate-movie-files", action="store_true", dest="separate_movie_files", default=False, help="Do not write the movie into a single file but use a separate POSCAR for every state. These are created in the \"movies\" directory. Only useful in conjunction with --movie.")
optpar.add_option("-n", "--no-submit", action="store_true", dest="no_submit", default=False,help="don't submit searches; only register finished results")
(options, args) = optpar.parse_args()
if len(args) > 1:
print("akmc.py takes only one positional argument")
config.init_from_cli(args)
if options.no_submit:
config.comm_job_buffer_size = 0
# The log and the lock both live under path_results. Metadata
# creates that directory only after the lock is held. The lock
# path is resolved from the working directory, so create that
# parent as well as path_results.
results = Path(config.path_results)
results.mkdir(parents=True, exist_ok=True)
lock_path = (Path.cwd() / results / "lockfile").absolute()
lock_path.parent.mkdir(parents=True, exist_ok=True)
#setup logging
logging.basicConfig(level=logging.DEBUG,
filename=str(results / "akmc.log"),
format="%(asctime)s %(levelname)s:%(name)s: %(message)s",
datefmt="%F %T")
logging.raiseExceptions = False
if not options.quiet:
rootlogger = logging.getLogger('')
console = logging.StreamHandler()
console.setLevel(logging.INFO)
formatter = logging.Formatter("%(message)s")
console.setFormatter(formatter)
rootlogger.addHandler(console)
lock = locking.LockFile(str(lock_path))
# Some options are mutually exclusive. Let's check them now.
exclusive_options = {}
if len(options.movie_type) > 0:
exclusive_options['movie_type'] = True
else:
exclusive_options['movie_type'] = False
exclusive_options['print_status'] = options.print_status
exclusive_options['reset'] = options.reset
if sum(exclusive_options.values()) > 1:
offending_options = [ k for k,v in list(exclusive_options.items()) if v ]
optpar.error("Options %s are mutually exclusive" % ", ".join(offending_options))
if len(options.movie_type) > 0:
states = get_statelist(config.main_temperature / 11604.5, config=config)
movie.make_movie(options.movie_type, config.path_root, states,
options.separate_movie_files)
sys.exit(0)
# From the config file: The Novotny and C&V (ASKMC) methods should not be used together.
if config.sb_on and config.askmc_on:
logger.error("Both superbasin methods should not be used at the same time")
sys.exit(1)
# From the config file: dynamics confidence can not be used with MinMode searches
if config.saddle_method=="min_mode" and config.akmc_confidence_scheme=="dynamics":
logger.error("Dynamics confidence scheme can not be used with MinMode saddle searches")
sys.exit(1)
if options.print_status:
states = get_statelist(config.main_temperature / 11604.5, config=config)
start_state_num, time, previous_state_num, first_run, previous_temperature =\
get_akmc_metadata(config)
current_state = states.get_state(start_state_num)
if config.sb_on:
sb_scheme = get_superbasin_scheme(states, config)
sb = sb_scheme.get_containing_superbasin(current_state)
else:
sb = None
print()
print("General")
print("-------")
if not sb:
print("Current state:", start_state_num)
else:
print("Current state:", start_state_num, "in superbasin", sb.id)
print("Number of states:",states.get_num_states())
print("Time simulated: %.3e seconds" % time)
print()
print("Current State")
print("-------------")
if not sb:
print("Confidence: %.4f" % current_state.get_confidence())
else:
print("Superbasin Confidence: %.4f" % sb.get_confidence())
non_ignored_states = set(sb._get_filtered_states())
for s in sb.states:
if s in non_ignored_states:
ignore_string = ""
else:
ignore_string = " (no exit from superbasin found)"
print(" %4i: %.4f%s" % (s.number, s.get_confidence(sb), ignore_string))
print("Unique Saddles:", current_state.get_unique_saddle_count())
print("Good Saddles:", current_state.get_good_saddle_count())
print("Bad Saddles:", current_state.get_bad_saddle_count())
print("Percentage bad saddles: %.1f" % (float(current_state.get_bad_saddle_count())/float(max(current_state.get_bad_saddle_count() + current_state.get_good_saddle_count(), 1)) * 100))
print()
comm = communicator.get_communicator(config)
print("Saddle Searches")
print("---------------")
print("Searches in queue:", comm.get_queue_size())
print()
if config.sb_on:
print("Superbasins")
print("-----------")
for i in sb_scheme.superbasins:
print("%s: %s" % (i.id, i.state_numbers))
sys.exit(0)
elif options.reset:
if options.force:
res = 'y'
else:
res = input("Are you sure you want to reset (all data files will be lost)? (y/N) ").lower()
if len(res)>0 and res[0] == 'y':
def attempt_removal(thing):
if thing is None:
return
p = Path(thing)
if p.is_dir():
io.remove_tree_and_empty_parents(p)
elif p.is_file():
p.unlink()
root = Path(config.path_root)
results = Path(config.path_results)
rmthings = [config.path_jobs_out,
config.path_jobs_in,
config.path_incomplete,
config.path_states,
config.path_scratch,
config.kdb_name,
config.kdb_path,
config.kdb_scratch_path,
config.sb_path,
config.sb_recycling_path,
config.debug_results_path,
root / "searchdata",
results / "askmc_data.txt",
results / "searches.log",
results / "dynamics.txt",
io.info_txt_path(config),
results / "akmc.log",
results / "jobs.tbl",
root / "results",
io.prng_state_path(config),
root / "explorer.pickle",
root / "temperatures.dat",
root / "client.log",
root / "lockfile",
]
for thing in rmthings:
attempt_removal(thing)
if not options.quiet:
print("Reset")
sys.exit(0)
else:
print("Not resetting")
sys.exit(1)
elif options.restart:
string_sb_clear = ""
if options.force:
res = 'y'
else:
res = input("Are you sure you want to restart (remove dynamics.txt, info.txt and akmc.log)? (y/N) ").lower()
if len(res)>0 and res[0] == 'y':
results = Path(config.path_results)
for i in [
results / "dynamics.txt",
io.info_txt_path(config),
results / "akmc.log",
results / "jobs.tbl",
]:
p = Path(i)
if p.is_file():
p.unlink()
if config.sb_on:
if options.force:
res = 'y'
else:
res = input("Should the superbasins be removed? (y/N) ").lower()
# remove superbasin data (specific for a trajectory)
if len(res)>0 and res[0] == 'y':
if Path(config.sb_path).is_dir():
io.remove_tree_and_empty_parents(config.sb_path)
states_root = Path(config.path_states)
if states_root.is_dir():
for state_dir in states_root.iterdir():
if state_dir.name == 'state_table':
continue
superbasin_file = state_dir / config.sb_state_file
if superbasin_file.is_file():
superbasin_file.unlink()
string_sb_clear = " with directory 'superbasins' and files named '"
string_sb_clear += str(config.sb_state_file) + "' removed"
if not options.quiet:
print("Restart"+string_sb_clear+".")
sys.exit(0)
else:
print("Not restarting")
sys.exit(1)
if lock.acquirelock():
if options.continuous or config.comm_type == 'mpi':
# define a wait method.
if config.comm_type == 'mpi':
from eon.mpiwait import mpiwait
wait = lambda: mpiwait(config.mpi_poll_period)
elif options.continuous:
if config.comm_type == "local":
# In local, everything is synchronous, so no need to wait here.
wait = lambda: None
else:
wait = lambda: sleep(10.0)
else:
raise RuntimeError("You have found a bug in eOn!")
# Run a specified number of steps or forever.
steps = 0
while True:
steps = akmc(config, steps)
if (config.akmc_max_kmc_steps > 0 and
steps >= config.akmc_max_kmc_steps):
break
wait()
if config.comm_type == 'mpi':
comm = communicator.get_communicator(config)
comm.stop_clients()
else:
akmc(config)
else:
logger.info("Server is locked by pid %i" % lock.pid)
sys.exit(1)
if __name__ == '__main__':
main()
eOn