# 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
unit_configuration = BulkConfiguration(
    bravais_lattice=lattice,
    elements=elements,
    fractional_coordinates=fractional_coordinates
    )
reps = 3
bulk_configuration = unit_configuration.repeat(reps, reps, reps)
bulk_configuration = bulk_configuration.copyAndDeleteAtoms(5)

# Setup the calculator.
calculator = LCAOCalculator(
    exchange_correlation=GGA.PBE,
    numerical_accuracy_parameters=NumericalAccuracyParameters(
        k_point_sampling=MonkhorstPackGrid(1, 1, 1))
    )

bulk_configuration.setCalculator(calculator)
bulk_configuration.update()

# Compute the electronic inverse participation ratio
electronic_ipr = ElectronicInverseParticipationRatio(
    configuration=bulk_configuration
    )

print('\nNumber of atoms N:', bulk_configuration.numberOfAtoms())
print(f'Expected IPR for delocalized states 1/N: {1 / bulk_configuration.numberOfAtoms():.3f}')

# Obtain results
energies = electronic_ipr.energies().inUnitsOf(eV)
ipr = electronic_ipr.inverseParticipationRatio()
idx = numpy.argmax(ipr)
print(f'State at energy {energies.ravel()[idx]:.2f} eV',
      f'has the highest IPR: {ipr.ravel()[idx]:.3f}')