setVerbosity(MinimalLog)

# -------------------------------------------------------------
# Bulk Configuration
# -------------------------------------------------------------

# Read in configuration
bulk_configuration = nlread('../sib-optimized.nc', BulkConfiguration)[-1]

# -------------------------------------------------------------
# Calculator
# -------------------------------------------------------------
#----------------------------------------
# Basis Set
#----------------------------------------
basis_set = [
    GGABasis.Boron_DoubleZetaPolarized,
    GGABasis.Silicon_SingleZetaPolarized,
    ]

#----------------------------------------
# Exchange-Correlation
#----------------------------------------
exchange_correlation = GGA.PBES

k_point_sampling = MonkhorstPackGrid(
    na=1,
    nb=1,
    nc=1,
    shift_to_gamma=[True, True, True],
    )
numerical_accuracy_parameters = NumericalAccuracyParameters(
    k_point_sampling=k_point_sampling,
    density_mesh_cutoff=50.0*Hartree
    )

parallel_parameters = ParallelParameters(
    processes_per_saddle_search=1,
    )

calculator = LCAOCalculator(
    basis_set=basis_set,
    exchange_correlation=exchange_correlation,
    numerical_accuracy_parameters=numerical_accuracy_parameters,
    parallel_parameters=parallel_parameters,
    )

bulk_configuration.setCalculator(calculator)
nlprint(bulk_configuration)
bulk_configuration.update()

# --------------------------------------------------------------
# AKMC
# --------------------------------------------------------------

# Reusing existing MarkovChain object if it already exists, otherwise a new one is created
if os.path.isfile('akmc_markov_chain.nc'):
    markov_chain = nlread('akmc_markov_chain.nc')[0]
else:
    markov_chain = MarkovChain(bulk_configuration, TotalEnergy(bulk_configuration).evaluate())

# Reusing existing KMC object if it already exists, otherwise a new one is created
if os.path.isfile('akmc_kmc.nc'):
    kmc = nlread('akmc_kmc.nc')[0]
else:
    kmc = None

# Modify the default maximum number of NEB images to 5.
saddle_search_parameters = SaddleSearchParameters(max_neb_images=5)

# Setup the AKMC simulation.
akmc = AdaptiveKineticMonteCarlo(markov_chain,
                                 calculator=calculator,
                                 kmc_temperature=300*Kelvin,
                                 md_temperature=3000*Kelvin,
                                 saddle_search_parameters=saddle_search_parameters,
                                 constraints=[15],
				 kmc=kmc,
                                 confidence=0.99,
				 assumed_prefactor=1e13/Second,
                                 write_searches=False,
                                 write_kmc=True,
                                 write_markov_chain=True,
                                 write_log=True)

# Run 200 saddle searches.
akmc.run(200)