# %% Silicon (alpha)

# Set up lattice
lattice = FaceCenteredCubic(5.4306*Angstrom)

# Define elements
elements = [Silicon, Silicon]

# Define coordinates
fractional_coordinates = [[ 0.  ,  0.  ,  0.  ],
                          [ 0.25,  0.25,  0.25]]

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


# %% Set SemiEmpiricalCalculator

# %% SemiEmpiricalCalculator

hamiltonian_parametrization = SlaterKosterHamiltonianParametrization(
    basis_set=Bassani.Si_Basis)

k_point_sampling = KpointDensity(
    density_a=4.0*Angstrom,
    density_b=4.0*Angstrom,
    density_c=4.0*Angstrom
)

numerical_accuracy_parameters = NumericalAccuracyParameters(
    k_point_sampling=k_point_sampling
)

calculator = SemiEmpiricalCalculator(
    hamiltonian_parametrization=hamiltonian_parametrization,
    numerical_accuracy_parameters=numerical_accuracy_parameters,
    checkpoint_handler=NoCheckpointHandler
)


# %% Set Calculator

silicon_alpha.setCalculator(calculator)

silicon_alpha.update()

nlsave('results.hdf5', silicon_alpha)


# %% BlochStatesGenerator

bloch_states_generator = calculateBlochStates(
    configuration=silicon_alpha,
    quantum_numbers=(0, 2, 3)
)

for bloch_state in bloch_states_generator:
    quantum_number = bloch_state.quantumNumber()
    nlsave('results.hdf5', bloch_state, object_id=f'bloch_state_{quantum_number}')