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

# Set up lattice
lattice = Hexagonal(4.138*Angstrom, 28.64*Angstrom)

# Define elements
elements = [Selenium, Selenium, Selenium, Bismuth, Bismuth, Bismuth, Bismuth,
            Bismuth, Bismuth, Selenium, Selenium, Selenium, Selenium, Selenium,
            Selenium]

# Define coordinates
fractional_coordinates = [[ 0.        ,  0.        ,  0.        ],
                          [ 0.33333333,  0.66666667,  0.66666667],
                          [ 0.66666667,  0.33333333,  0.33333333],
                          [ 0.66666667,  0.33333333,  0.73233333],
                          [ 0.33333333,  0.66666667,  0.26766667],
                          [ 0.        ,  0.        ,  0.601     ],
                          [ 0.66666667,  0.33333333,  0.93433333],
                          [ 0.        ,  0.        ,  0.399     ],
                          [ 0.33333333,  0.66666667,  0.06566667],
                          [ 0.33333333,  0.66666667,  0.87266667],
                          [ 0.66666667,  0.33333333,  0.12733333],
                          [ 0.33333333,  0.66666667,  0.46066667],
                          [ 0.        ,  0.        ,  0.794     ],
                          [ 0.        ,  0.        ,  0.206     ],
                          [ 0.66666667,  0.33333333,  0.53933333]]

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

# -------------------------------------------------------------
# Calculator
# -------------------------------------------------------------
#----------------------------------------
# Basis Set
#----------------------------------------

SeleniumBasis = OpenMXBasisSet(
    element=PeriodicTable.Selenium,
    filename="openmx/pao/Se7.0.pao.zip",
    atomic_species="s2p2d1",
    pseudopotential=NormConservingPseudoPotential("normconserving/upf2/Se_PBE13.upf.zip"),
    )


BismuthBasis = OpenMXBasisSet(
    element=PeriodicTable.Bismuth,
    filename="openmx/pao/Bi8.0.pao.zip",
    atomic_species="s2p2d2f1",
    pseudopotential=NormConservingPseudoPotential("normconserving/upf2/Bi_PBE13.upf.zip"),
    )

basis_set = [
    SeleniumBasis,
    BismuthBasis,
    ]

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

numerical_accuracy_parameters = NumericalAccuracyParameters(
    k_point_sampling=(9, 9, 3),
    density_mesh_cutoff=150.0*Hartree,
    )

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

bulk_configuration.setCalculator(calculator)

# -------------------------------------------------------------
# Initial State
# -------------------------------------------------------------

bulk_configuration.update()
nlsave('Bi2Se3_bulk.nc', bulk_configuration)
nlprint(bulk_configuration)

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

# -------------------------------------------------------------
# Calculator
# -------------------------------------------------------------
#----------------------------------------
# Basis Set
#----------------------------------------

SeleniumBasis = OpenMXBasisSet(
    element=PeriodicTable.Selenium,
    filename="openmx/pao/Se7.0.pao.zip",
    atomic_species="s2p2d1",
    pseudopotential=NormConservingPseudoPotential("normconserving/upf2/Se_PBE13.upf.zip"),
    )


BismuthBasis = OpenMXBasisSet(
    element=PeriodicTable.Bismuth,
    filename="openmx/pao/Bi8.0.pao.zip",
    atomic_species="s2p2d2f1",
    pseudopotential=NormConservingPseudoPotential("normconserving/upf2/Bi_PBE13.upf.zip"),
    )

basis_set = [
    SeleniumBasis,
    BismuthBasis,
    ]

#----------------------------------------
# Exchange-Correlation
#----------------------------------------
exchange_correlation = SOGGA.PBE

numerical_accuracy_parameters = NumericalAccuracyParameters(
    k_point_sampling=(9, 9, 3),
    density_mesh_cutoff=150.0*Hartree,
    )

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

bulk_configuration.setCalculator(calculator)

# -------------------------------------------------------------
# Initial State
# -------------------------------------------------------------
initial_spin = InitialSpin(scaled_spins=[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
old_calculation = nlread('Bi2Se3_bulk.nc', BulkConfiguration)[0]
bulk_configuration.setCalculator(
    calculator,
    initial_spin=initial_spin,
    initial_state=old_calculation,
)
bulk_configuration.update()
nlsave('Bi2Se3_bulk.nc', bulk_configuration)
nlprint(bulk_configuration)

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