# -*- coding: utf-8 -*-
# -------------------------------------------------------------
# Two-probe Configuration
# -------------------------------------------------------------

filename = '../dyn-mat/ff-device-si.hdf5'
device_configuration = nlread(filename, DeviceConfiguration)[-1]

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

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

#----------------------------------------
# Numerical Accuracy Settings
#----------------------------------------
left_electrode_k_point_sampling = MonkhorstPackGrid(
    na=9,
    nb=9,
    nc=101,
    )
left_electrode_numerical_accuracy_parameters = NumericalAccuracyParameters(
    occupation_method=FermiDirac(300.0*Kelvin*boltzmann_constant),
    k_point_sampling=left_electrode_k_point_sampling,
    density_mesh_cutoff=75.0*Hartree,
    )

right_electrode_k_point_sampling = MonkhorstPackGrid(
    na=9,
    nb=9,
    nc=101,
    )
right_electrode_numerical_accuracy_parameters = NumericalAccuracyParameters(
    occupation_method=FermiDirac(300.0*Kelvin*boltzmann_constant),
    k_point_sampling=right_electrode_k_point_sampling,
    density_mesh_cutoff=75.0*Hartree,
    )

device_k_point_sampling = MonkhorstPackGrid(
    na=9,
    nb=9,
    nc=101,
    )
device_numerical_accuracy_parameters = NumericalAccuracyParameters(
    occupation_method=FermiDirac(300.0*Kelvin*boltzmann_constant),
    k_point_sampling=device_k_point_sampling,
    density_mesh_cutoff=75.0*Hartree,
    )

#----------------------------------------
# Iteration Control Settings
#----------------------------------------
left_electrode_iteration_control_parameters = IterationControlParameters(
    number_of_history_steps=10,
    max_steps = 1000,
    tolerance=1e-07,
    )

right_electrode_iteration_control_parameters = IterationControlParameters(
    number_of_history_steps=10,
    max_steps = 1000,
    tolerance=1e-07,
    )

device_iteration_control_parameters = IterationControlParameters(
    number_of_history_steps=10,
    max_steps = 1000,
    damping_factor = 0.01,
    )

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

right_electrode_poisson_solver = ParallelConjugateGradientSolver(
    boundary_conditions=[[PeriodicBoundaryCondition(),PeriodicBoundaryCondition()],
                         [PeriodicBoundaryCondition(),PeriodicBoundaryCondition()],
                         [PeriodicBoundaryCondition(),PeriodicBoundaryCondition()]]
    )

device_poisson_solver = ParallelConjugateGradientSolver(
    boundary_conditions=[[PeriodicBoundaryCondition(),PeriodicBoundaryCondition()],
                         [PeriodicBoundaryCondition(),PeriodicBoundaryCondition()],
                         [DirichletBoundaryCondition(),DirichletBoundaryCondition()]]
    )

#----------------------------------------
# Contour Integral Settings
#----------------------------------------
equilibrium_contour = SemiCircleContour(
    integral_lower_bound=1.48099775058*Hartree,
    )
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,
    iteration_control_parameters=left_electrode_iteration_control_parameters,
    poisson_solver=left_electrode_poisson_solver,
    )

right_electrode_calculator = LCAOCalculator(
    basis_set=basis_set,
    exchange_correlation=exchange_correlation,
    numerical_accuracy_parameters=right_electrode_numerical_accuracy_parameters,
    iteration_control_parameters=right_electrode_iteration_control_parameters,
    poisson_solver=right_electrode_poisson_solver,
    )

#----------------------------------------
# Device Calculator
#----------------------------------------
calculator = DeviceLCAOCalculator(
    basis_set=basis_set,
    exchange_correlation=exchange_correlation,
    numerical_accuracy_parameters=device_numerical_accuracy_parameters,
    iteration_control_parameters=device_iteration_control_parameters,
    poisson_solver=device_poisson_solver,
    contour_parameters=contour_parameters,
    electrode_calculators=
        [left_electrode_calculator, right_electrode_calculator],
    )

device_configuration.setCalculator(calculator)
nlprint(device_configuration)
device_configuration.update()
nlsave('0K-zero-bias.hdf5', device_configuration)

# -------------------------------------------------------------
# Projected Local Density Of States
# -------------------------------------------------------------
kpoint_grid = MonkhorstPackGrid(
    na=21,
    nb=21,
    )

projected_local_density_of_states = ProjectedLocalDensityOfStates(
    configuration=device_configuration,
    method=DeviceDensityOfStates,
    energies=numpy.linspace(-1, 1, 101)*eV,
    kpoints=kpoint_grid,
    contributions=All,
    self_energy_calculator=RecursionSelfEnergy(),
    energy_zero_parameter=AverageFermiLevel,
    infinitesimal=1e-06*eV,
    )
nlsave('0K-zero-bias.hdf5', projected_local_density_of_states)