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

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

# Set up lattice
lattice = SimpleOrthorhombic(5.546690740087016*Angstrom, 3.2093295699579034*Angstrom, 5.241822783047967*Angstrom)

# Define elements
elements = [Gallium, Gallium, Nitrogen, Nitrogen, Gallium, Gallium, Nitrogen,
            Nitrogen]

# Define coordinates
fractional_coordinates = [[ 0.338086773247, -0.            ,  0.005917824105],
                          [ 0.661913226753,  0.            ,  0.505917824105],
                          [ 0.329542526042,  0.            ,  0.379082175895],
                          [ 0.670457473958, -0.            ,  0.879082175895],
                          [ 0.838086773247,  0.5           ,  0.005917824105],
                          [ 0.161913226753,  0.5           ,  0.505917824105],
                          [ 0.829542526042,  0.5           ,  0.379082175895],
                          [ 0.170457473958,  0.5           ,  0.879082175895]]

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

# -------------------------------------------------------------
# Calculator
# -------------------------------------------------------------
#----------------------------------------
# Exchange-Correlation
#----------------------------------------
exchange_correlation = MGGA.SCAN

basis_set = [
    BasisGGAPseudoDojo.Nitrogen_Medium,
    BasisGGAPseudoDojo.Gallium_Medium,
    BasisGGAPseudoDojo.Carbon_Medium,
    ]

k_point_sampling = KpointDensity(
    density_a=4.0*Angstrom,
    )
numerical_accuracy_parameters = NumericalAccuracyParameters(
    density_mesh_cutoff=95.0*Hartree,
    k_point_sampling=k_point_sampling,
    )

calculator = LCAOCalculator(
    exchange_correlation=MGGA.SCAN,
    basis_set=basis_set,
    numerical_accuracy_parameters=numerical_accuracy_parameters,
    )

bulk_configuration.setCalculator(calculator)
nlprint(bulk_configuration)
bulk_configuration.update()
nlsave('charge_point_defect.hdf5', bulk_configuration)

# -------------------------------------------------------------
# Charged Point Defect
# -------------------------------------------------------------
# Defect
point_defect = Substitutional(
    element=Carbon,
    site_index=2,
    unit_cell_index=(0, 0, 0))

# Charge states
charge_states = [-1, 0]

# Atomic chemical potentials
atomic_chemical_potentials = [
    AtomicChemicalPotential(element=Carbon, internal_energy=None),
    AtomicChemicalPotential(element=Nitrogen, internal_energy=None)
]

# Supercell repetitions
supercell_repetitions_list = [
    [2, 3, 2],
]

# File name.
filename = 'charge_point_defect.hdf5'

band_gap_calculator = calculator(exchange_correlation=HybridGGA.HSE06)

charged_point_defect = ChargedPointDefect(
    bulk_configuration=bulk_configuration,
    formation_energy_calculator=calculator,
    band_gap_calculator=band_gap_calculator,
    filename=filename,
    object_id='charged_point_defect',
    point_defect=point_defect,
    charge_states=charge_states,
    atomic_chemical_potentials=atomic_chemical_potentials,
    supercell_repetitions_list=supercell_repetitions_list,
    relax_atomic_coordinates=True,
    dielectric_constant=8.9,
    log_filename_prefix='charged_point_defect_',
    number_of_processes_per_task=None,
    elastic_correction_method=None,
    optimize_geometry_parameters=OptimizeGeometryParameters(max_forces=0.01*eV/Angstrom),
    random_seed=1
)
charged_point_defect.update()