# -*- coding: utf-8 -*-
from bulks import si, ge, sige
setVerbosity(MinimalLog)

configurations = [si,ge,sige]
labels         = ['Si','Ge','SiGe']

# Basis set for Silicon
SiliconBasis_projector_shift = PseudoPotentialProjectorShift(
    s_orbital_shift=21.33*eV,
    p_orbital_shift=-1.43*eV,
    d_orbital_shift=0.0*eV,
    f_orbital_shift=0.0*eV,
    g_orbital_shift=0.0*eV
    )
SiliconBasis = BasisGGASG15.Silicon_Medium(projector_shift=SiliconBasis_projector_shift)

# Basis set for Germanium
GermaniumBasis_projector_shift = PseudoPotentialProjectorShift(
    s_orbital_shift=14.55*eV,
    p_orbital_shift=0.20*eV,
    d_orbital_shift=-2.17*eV,
    f_orbital_shift=0.0*eV,
    g_orbital_shift=0.0*eV
    )
GermaniumBasis = BasisGGASG15.Germanium_Medium(projector_shift=GermaniumBasis_projector_shift)

# Combined basis sets
basis_sets = [[SiliconBasis], [GermaniumBasis], [SiliconBasis, GermaniumBasis]]

for bulk_configuration,label,basis_set in zip(configurations,labels,basis_sets):
    outfile = "%s_PPS.hdf5" % label

    # -------------------------------------------------------------
    # Calculator
    # -------------------------------------------------------------
    k_point_sampling = MonkhorstPackGrid(
        na=9,
        nb=9,
        nc=9,
        )
    numerical_accuracy_parameters = NumericalAccuracyParameters(
        k_point_sampling=k_point_sampling,
        density_mesh_cutoff=90.0*Hartree,
        )

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

    bulk_configuration.setCalculator(calculator)
    nlprint(bulk_configuration)
    bulk_configuration.update()
    nlsave(outfile, bulk_configuration)

    # -------------------------------------------------------------
    # Optimize Geometry
    # -------------------------------------------------------------
    fix_atom_indices_0 = [0, 1]
    constraints = [FixAtomConstraints(fix_atom_indices_0)]

    bulk_configuration = OptimizeGeometry(
        bulk_configuration,
        max_forces=0.05*eV/Ang,
        max_stress=0.1*GPa,
        max_steps=200,
        max_step_length=0.2*Ang,
        constraints=constraints,
        trajectory_filename=None,
        optimizer_method=LBFGS(),
        constrain_bravais_lattice=True,
    )
    nlsave(outfile, bulk_configuration)
    nlprint(bulk_configuration)

    # -------------------------------------------------------------
    # Bandstructure
    # -------------------------------------------------------------
    bandstructure = Bandstructure(
        configuration=bulk_configuration,
        route=['L', 'G', 'X'],
        points_per_segment=50
        )
    nlsave(outfile, bandstructure)