# -*- coding: utf-8 -*-
setVerbosity(MinimalLog)

# -------------------------------------------------------------
# Left Electrode
# -------------------------------------------------------------

# Set up lattice
vector_a = [2.88902617589, 0.0, 0.0]*Angstrom
vector_b = [-1.44451308795, 2.50197006052, 0.0]*Angstrom
vector_c = [0.0, 0.0, 7.07663998448]*Angstrom
left_electrode_lattice = UnitCell(vector_a, vector_b, vector_c)

# Define elements
left_electrode_elements = [Silver, Silver, Silver]

# Define coordinates
left_electrode_coordinates = [[-0.            ,  1.667980040348,  1.179540023274],
                              [-0.            ,  0.            ,  3.538420018102],
                              [ 1.444513087947,  0.833990020174,  5.89730001293 ]]*Angstrom

# Set up configuration
left_electrode = BulkConfiguration(
    bravais_lattice=left_electrode_lattice,
    elements=left_electrode_elements,
    cartesian_coordinates=left_electrode_coordinates
    )

# -------------------------------------------------------------
# Central Region
# -------------------------------------------------------------

# Set up lattice
vector_a = [2.88902617589, 0.0, 0.0]*Angstrom
vector_b = [-1.44451308795, 2.50197006052, 0.0]*Angstrom
vector_c = [0.0, 0.0, 44.7683399716]*Angstrom
central_region_lattice = UnitCell(vector_a, vector_b, vector_c)

# Define elements
central_region_elements = [Silver, Silver, Silver, Silver, Silver, Silver, Silver, Silver,
                           Silver, Silver, Silver]

# Define coordinates
central_region_coordinates = [[ -0.            ,   1.667980040348,   1.179540023274],
                              [ -0.            ,   0.            ,   3.538420018102],
                              [  1.444513087947,   0.833990020174,   5.89730001293 ],
                              [ -0.            ,   1.667980040348,   8.256180007758],
                              [ -0.            ,   0.            ,  10.615060002586],
                              [  1.444513087947,   0.833990020174,  12.973939997414],
                              [ -0.            ,   1.667980040348,  15.332819992242],
                              [ -0.            ,   0.            ,  17.69169998707 ],
                              [  1.444513087947,   0.833990020174,  20.050579981898],
                              [ -0.            ,   1.667980040348,  22.409459976726],
                              [ -0.            ,   0.            ,  24.768339971554]]*Angstrom

# Set up configuration
central_region = BulkConfiguration(
    bravais_lattice=central_region_lattice,
    elements=central_region_elements,
    cartesian_coordinates=central_region_coordinates
    )

surface_configuration = SurfaceConfiguration(
    central_region,
    left_electrode,
    equivalent_electrode_length=7.07663998448*Angstrom,
    )

# -------------------------------------------------------------
# Calculator
# -------------------------------------------------------------
#----------------------------------------
# Basis Set
#----------------------------------------
basis_set = [
    LDABasis.Silver_SingleZetaPolarized,
    ]

#----------------------------------------
# Exchange-Correlation
#----------------------------------------
exchange_correlation = LDA.PZ

#----------------------------------------
# Numerical Accuracy Settings
#----------------------------------------
left_electrode_k_point_sampling = MonkhorstPackGrid(
    na=9,
    nb=9,
    nc=231,
    )
left_electrode_numerical_accuracy_parameters = NumericalAccuracyParameters(
    density_mesh_cutoff=45.0*Hartree,
    k_point_sampling=left_electrode_k_point_sampling,
    )

device_k_point_sampling = MonkhorstPackGrid(
    na=9,
    nb=9,
    nc=231,
    )
device_numerical_accuracy_parameters = NumericalAccuracyParameters(
    density_mesh_cutoff=45.0*Hartree,
    k_point_sampling=device_k_point_sampling,
    )

#----------------------------------------
# Poisson Solver Settings
#----------------------------------------
left_electrode_poisson_solver = FastFourier2DSolver(
    boundary_conditions=[[PeriodicBoundaryCondition(),PeriodicBoundaryCondition()],
                         [PeriodicBoundaryCondition(),PeriodicBoundaryCondition()],
                         [PeriodicBoundaryCondition(),PeriodicBoundaryCondition()]]
    )

#----------------------------------------
# Contour Integral Settings
#----------------------------------------
equilibrium_contour = SemiCircleContour(
    integral_lower_bound=1.17597836274*Hartree,
    circle_points=30,
    )
contour_parameters = ContourParameters(
    equilibrium_contour=equilibrium_contour,
    )

#----------------------------------------
# Electrode Calculators
#----------------------------------------
left_electrode_calculator = LCAOCalculator(
    basis_set=basis_set,
    exchange_correlation=exchange_correlation,
    numerical_accuracy_parameters=left_electrode_numerical_accuracy_parameters,
    poisson_solver=left_electrode_poisson_solver,
    )

#----------------------------------------
# Device Calculator
#----------------------------------------
calculator = DeviceLCAOCalculator(
    basis_set=basis_set,
    exchange_correlation=exchange_correlation,
    numerical_accuracy_parameters=device_numerical_accuracy_parameters,
    contour_parameters=contour_parameters,
    electrode_calculators=
        [left_electrode_calculator, left_electrode_calculator],
    )

surface_configuration.setCalculator(calculator)
nlprint(surface_configuration)
surface_configuration.update()
nlsave('Surface Silver 111 LDA.hdf5', surface_configuration)

# -------------------------------------------------------------
# Surface Bandstructure
# -------------------------------------------------------------
surface_bandstructure = SurfaceBandstructure(
    configuration=surface_configuration,
    route=['X', 'G', 'Y'],
    points_per_segment=150,
    energies=numpy.linspace(-1.5, 1.5, 300)*eV,
    projection_list=ProjectionList(atoms=All),
    contributions=All,
    self_energy_calculator=RecursionSelfEnergy(storage_strategy=NoStorage()),
    energy_zero_parameter=AverageFermiLevel,
    infinitesimal=1e-03*eV,
    )
nlsave('Surface Silver 111 LDA.hdf5', surface_bandstructure)
nlprint(surface_bandstructure)