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

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

# Set up lattice
lattice = SimpleTetragonal(2.866*Angstrom, 45.34172044506112*Angstrom)

# Define elements
elements = [Iron, Iron, Iron, Iron, Iron, Iron, Iron, Iron, Iron, Iron, Iron,
            Iron, Iron, Iron, Iron, Iron, Iron, Iron, Iron, Iron, Iron]

# Define coordinates
fractional_coordinates = [[ 0.248796793653,  0.248796793651,  0.015236052172],
                          [ 0.749240954235,  0.749240954234,  0.045902009605],
                          [ 0.248800401505,  0.248800401504,  0.076577038988],
                          [ 0.749235217307,  0.749235217306,  0.107154778896],
                          [ 0.248894116126,  0.248894116127,  0.137846827345],
                          [ 0.749635120547,  0.749635120548,  0.168439264182],
                          [ 0.249677932107,  0.249677932108,  0.199221654535],
                          [ 0.750847797286,  0.750847797287,  0.229845723164],
                          [ 0.251186736882,  0.251186736884,  0.260914588726],
                          [ 0.752082597185,  0.752082597186,  0.292712488503],
                          [ 0.25222413526 ,  0.252224135259,  0.322027504254],
                          [ 0.248434621128,  0.248434621125,  0.708899029361],
                          [ 0.74905639095 ,  0.749056390949,  0.738051844042],
                          [ 0.248850456433,  0.24885045643 ,  0.769766600765],
                          [ 0.749370078138,  0.749370078136,  0.800746712211],
                          [ 0.248874627157,  0.248874627152,  0.831322119996],
                          [ 0.749284193202,  0.7492841932  ,  0.862064593914],
                          [ 0.248816407684,  0.248816407683,  0.892637429232],
                          [ 0.749242343327,  0.749242343324,  0.923315868123],
                          [ 0.248804076995,  0.248804076993,  0.953895871831],
                          [ 0.749225236779,  0.749225236777,  0.984566314048]]

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

# Add tags
bulk_configuration.addTags('Selection 0', [0, 1, 2, 3, 4, 5, 16, 17, 18, 19, 20])

# -------------------------------------------------------------
# Calculator
# -------------------------------------------------------------
#----------------------------------------
# Basis Set
#----------------------------------------
basis_set = [
    BasisGGAPseudoDojoSO.Iron_Medium,
    ]

#----------------------------------------
# Exchange-Correlation
#----------------------------------------
exchange_correlation = SOGGA.PBE

k_point_sampling = KpointDensity(
    density_a=7.0*Angstrom,
    force_timereversal=False,
    )
numerical_accuracy_parameters = NumericalAccuracyParameters(
    density_mesh_cutoff=140.0*Hartree,
    k_point_sampling=k_point_sampling,
    occupation_method=FermiDirac(300.0*Kelvin*boltzmann_constant),
    )

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

bulk_configuration.setCalculator(calculator)
nlprint(bulk_configuration)
nlsave('Fe21-vacuum.hdf5', bulk_configuration)

# -------------------------------------------------------------
# Magnetic Anisotropy Energy
# -------------------------------------------------------------
# Kpoint sampling
kpoint_grid = KpointDensity(
    density_a=17.0*Angstrom,
    )

# File name.
filename = 'Fe21-vacuum.hdf5'

magnetic_anisotropy_energy = MagneticAnisotropyEnergy(
    configuration=bulk_configuration,
    filename=filename,
    object_id='magnetic_anisotropy_energy',
    theta_angles=numpy.linspace(0, 90, 2)*Degrees,
    phi_angles=numpy.linspace(0, 0, 1)*Degrees,
    kpoints=kpoint_grid,
    projections=ProjectOnSites,
    energy_minimum=-1000*eV,
    energy_maximum=1*eV,
    bands_per_electron=1.1,
    log_filename_prefix='magnetic_anisotropy_energy_',
    number_of_processes_per_task=None,
)

magnetic_anisotropy_energy.update()