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

# Set up lattice
lattice = Hexagonal(3.187*Angstrom, 7.72599*Angstrom)

# Define elements
elements = [Lithium, Lithium, Lithium, Lithium, Oxygen, Oxygen, Oxygen, Oxygen]

# Define coordinates
fractional_coordinates = [[ 0.            ,  0.            ,  0.            ],
                          [ 0.333333333333,  0.666666666667,  0.5           ],
                          [ 0.666666666667,  0.333333333333,  0.25          ],
                          [ 0.666666666667,  0.333333333333,  0.75          ],
                          [ 0.            ,  0.            ,  0.401         ],
                          [ 0.            ,  0.            ,  0.599         ],
                          [ 0.333333333333,  0.666666666667,  0.1           ],
                          [ 0.333333333333,  0.666666666667,  0.9           ]]

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

# -------------------------------------------------------------
# Calculator
# -------------------------------------------------------------
#----------------------------------------
# Basis Set
#----------------------------------------
basis_set = [
    GGABasis.Lithium_DoubleZetaPolarized,
    GGABasis.Oxygen_DoubleZetaPolarized,
    ]

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

numerical_accuracy_parameters = NumericalAccuracyParameters(
    k_point_sampling=(9, 9, 5),
    )

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

bulk_configuration.setCalculator(calculator)
nlprint(bulk_configuration)
bulk_configuration.update()
nlsave('Li2O2_opt.nc', bulk_configuration)

# -------------------------------------------------------------
# Optimize Geometry
# -------------------------------------------------------------
bulk_configuration = OptimizeGeometry(
        bulk_configuration,
        max_forces=0.02*eV/Ang,
        max_stress=0.002*eV/Ang**3,
        max_steps=200,
        max_step_length=0.2*Ang,
        trajectory_filename='Li2O2_opt.nc',
        optimizer_method=LBFGS(),
        )
nlsave('Li2O2_opt.nc', bulk_configuration)
nlprint(bulk_configuration)

# -------------------------------------------------------------
# Bandstructure
# -------------------------------------------------------------
bandstructure = Bandstructure(
    configuration=bulk_configuration,
    route=['G', 'M', 'L', 'A', 'G', 'K', 'H', 'A'],
    points_per_segment=101,
    bands_above_fermi_level=All
    )
nlsave('Li2O2_opt.nc', bandstructure)

# -------------------------------------------------------------
# Density Of States
# -------------------------------------------------------------
density_of_states = DensityOfStates(
    configuration=bulk_configuration,
    kpoints=MonkhorstPackGrid(9,9,5),
    energy_zero_parameter=FermiLevel,
    bands_above_fermi_level=40,
    )
nlsave('Li2O2_opt.nc', density_of_states)
nlprint(density_of_states)

# -------------------------------------------------------------
# Transmission Spectrum
# -------------------------------------------------------------
transmission_spectrum = TransmissionSpectrum(
    configuration=bulk_configuration,
    energies=numpy.linspace(-3,3,201)*eV,
    kpoints=MonkhorstPackGrid(25,25),
    energy_zero_parameter=AverageFermiLevel,
    infinitesimal=1e-06*eV,
    self_energy_calculator=RecursionSelfEnergy(),
    )
nlsave('Li2O2_opt.nc', transmission_spectrum)
nlprint(transmission_spectrum)