# -*- coding: utf-8 -*-

setVerbosity(MinimalLog)

# -------------------------------------------------------------
# Two-probe Configuration
# -------------------------------------------------------------

filename = '../0K-zero-bias/0K-zero-bias.hdf5'
filesave = 'ham-der.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=1,
    nb=1,
    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=1,
    nb=1,
    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=1,
    nb=1,
    nc=1,
    )
    
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)

# -------------------------------------------------------------
# Hamiltonian Derivatives
# -------------------------------------------------------------
hamiltonian_derivatives = HamiltonianDerivatives(
    configuration=device_configuration,
    repetitions=(3, 3, 1),
    atomic_displacement=0.01*Angstrom,
    constraints=[0,1,2,3,44,45,46,47],
    use_equivalent_bulk=False,
    )
nlsave(filesave, hamiltonian_derivatives)