markov_chain = nlread('akmc_markov_chain.nc')[0]

initial_conf = nlread('initial_conf.nc',BulkConfiguration)[-1]
final_conf = nlread('final_conf.nc',BulkConfiguration)[-1]

initial_state_id = 0
final_state_id = 17

# Make the NEB configuration

saddle_conf = markov_chain.getSaddleConfiguration(initial_state_id, final_state_id)

configuration_list = [initial_conf, saddle_conf, final_conf]

neb_configuration = NudgedElasticBand(configuration_list, image_distance=1.0*Angstrom, generate_images=True, wrap_pbc=True)

nlsave('neb-configuration.nc',neb_configuration)

# -------------------------------------------------------------
# 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
    )

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

# Optimize the NEB

neb_configuration.setCalculator(calculator)

optimized_neb_configuration = OptimizeNudgedElasticBand(
        neb_configuration,
        max_forces=0.01*eV/Ang)

nlsave('optimized-neb-configuration.nc',optimized_neb_configuration)

ci_optimized_neb_configuration = OptimizeNudgedElasticBand(
        optimized_neb_configuration,
	climbing_image=True,
        max_forces=0.01*eV/Ang)

nlsave('ci-optimized-neb-configuration.nc',ci_optimized_neb_configuration)