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

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

# Set up lattice
vector_a = [2.866, 0.0, 0.0]*Angstrom
vector_b = [0.0, 2.866, 0.0]*Angstrom
vector_c = [0.0, 0.0, 46.23210000000001]*Angstrom
lattice = UnitCell(vector_a, vector_b, vector_c)

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

# Define coordinates
fractional_coordinates = [[ 0.25          ,  0.25          ,  0.015497889994],
                          [ 0.75          ,  0.75          ,  0.046493669983],
                          [ 0.25          ,  0.25          ,  0.077489449971],
                          [ 0.75          ,  0.75          ,  0.108485229959],
                          [ 0.25          ,  0.25          ,  0.139481009948],
                          [ 0.75          ,  0.75          ,  0.170476789936],
                          [ 0.25          ,  0.25          ,  0.201472569924],
                          [ 0.75          ,  0.75          ,  0.232468349913],
                          [ 0.25          ,  0.25          ,  0.263464129901],
                          [ 0.75          ,  0.75          ,  0.294459909889],
                          [ 0.25          ,  0.25          ,  0.325455689878],
                          [ 0.25          ,  0.25          ,  0.373041674508],
                          [ 0.75          ,  0.75          ,  0.373041674508],
                          [ 0.25          ,  0.25          ,  0.42052707967 ],
                          [ 0.75          ,  0.75          ,  0.42052707967 ],
                          [ 0.25          ,  0.25          ,  0.468012484832],
                          [ 0.75          ,  0.75          ,  0.468012484832],
                          [ 0.25          ,  0.25          ,  0.515497889994],
                          [ 0.75          ,  0.75          ,  0.515497889994],
                          [ 0.25          ,  0.25          ,  0.562983295156],
                          [ 0.75          ,  0.75          ,  0.562983295156],
                          [ 0.25          ,  0.25          ,  0.610468700319],
                          [ 0.75          ,  0.75          ,  0.610468700319],
                          [ 0.25          ,  0.25          ,  0.657954105481],
                          [ 0.75          ,  0.75          ,  0.657954105481],
                          [ 0.25          ,  0.25          ,  0.705540090111],
                          [ 0.75          ,  0.75          ,  0.736535870099],
                          [ 0.25          ,  0.25          ,  0.767531650087],
                          [ 0.75          ,  0.75          ,  0.798527430076],
                          [ 0.25          ,  0.25          ,  0.829523210064],
                          [ 0.75          ,  0.75          ,  0.860518990052],
                          [ 0.25          ,  0.25          ,  0.891514770041],
                          [ 0.75          ,  0.75          ,  0.922510550029],
                          [ 0.25          ,  0.25          ,  0.953506330017],
                          [ 0.75          ,  0.75          ,  0.984502110006]]

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

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

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

k_point_sampling = KpointDensity(
    density_a=4.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-MgO7-unrelaxed-300K-140Ha.hdf5', bulk_configuration)

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

# File name.
filename = 'Fe21-MgO7-unrelaxed-300K-140Ha.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()