# -*- coding: utf-8 -*-
# -------------------------------------------------------------
# Bulk Configuration
# -------------------------------------------------------------

# Set up lattice
vector_a = [4.74148287577, -0.000127269878148, 0.00783573590836]*Angstrom
vector_b = [-0.000164005204931, 6.08030371318, -0.00256586481387]*Angstrom
vector_c = [0.0172414754673, -0.00445585886655, 10.3479508435]*Angstrom
lattice = UnitCell(vector_a, vector_b, vector_c)

# Define elements
elements = [Iron, Iron, Iron, Iron, Lithium, Lithium, Lithium, Phosphorus,
            Phosphorus, Phosphorus, Phosphorus, Oxygen, Oxygen, Oxygen, Oxygen,
            Oxygen, Oxygen, Oxygen, Oxygen, Oxygen, Oxygen, Oxygen, Oxygen,
            Oxygen, Oxygen, Oxygen, Oxygen]

# Define coordinates
fractional_coordinates = [[ 0.974025963402,  0.249803584653,  0.280873622179],
                          [ 0.474211659358,  0.750288731051,  0.220126866786],
                          [ 0.025975867482,  0.75019410548 ,  0.719126139844],
                          [ 0.525787461722,  0.249708336893,  0.779873152613],
                          [-0.000000330507, -0.000000548041, -0.000000353828],
                          [-0.00000055133 ,  0.500000328126,  0.000000490067],
                          [ 0.499999881749,  0.499999951689,  0.500000126016],
                          [ 0.418846558454,  0.248986560763,  0.092260345574],
                          [ 0.913583304052,  0.750820097461,  0.409195811152],
                          [ 0.581159529367,  0.751010865365,  0.907739715577],
                          [ 0.086411561798,  0.249176441361,  0.590800007724],
                          [ 0.744833815076,  0.248286207754,  0.091792522207],
                          [ 0.23907370572 ,  0.751598759587,  0.407281696103],
                          [ 0.255167963397,  0.751715186538,  0.908207940704],
                          [ 0.760924380509,  0.248402201525,  0.592720597394],
                          [ 0.221254366546,  0.248646397503,  0.451790852888],
                          [ 0.715394564227,  0.75169850186 ,  0.046471669566],
                          [ 0.778740044253,  0.751355302298,  0.548207293946],
                          [ 0.284608330972,  0.24830259825 ,  0.953523862775],
                          [ 0.282266525522,  0.04494155692 ,  0.165733669804],
                          [ 0.775242881473,  0.955583835088,  0.337294744727],
                          [ 0.717126972383,  0.546413954203,  0.834427511096],
                          [ 0.224205383627,  0.454437406911,  0.662603384446],
                          [ 0.717727886822,  0.955065433049,  0.834266883231],
                          [ 0.224763387198,  0.04441993238 ,  0.66270686361 ],
                          [ 0.282866738823,  0.453583101291,  0.165574151583],
                          [ 0.775802983568,  0.545560195365,  0.337400548719]]

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

# -------------------------------------------------------------
# Calculator
# -------------------------------------------------------------
#----------------------------------------
# Basis Set
#----------------------------------------
basis_set = [
    GGABasis.Lithium_DoubleZetaPolarized,
    GGABasis.Oxygen_DoubleZetaPolarized,
    GGABasis.Phosphorus_DoubleZetaPolarized,
    GGABasis.Iron_DoubleZetaPolarized,
    ]

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

k_point_sampling = MonkhorstPackGrid(
    na=7,
    nb=5,
    nc=3,
    )
numerical_accuracy_parameters = NumericalAccuracyParameters(
    k_point_sampling=k_point_sampling,
    )

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

bulk_configuration.setCalculator(calculator)
nlprint(bulk_configuration)
bulk_configuration.update()
nlsave('initial_C.nc', bulk_configuration)

# -------------------------------------------------------------
# Optimize Geometry
# -------------------------------------------------------------
bulk_configuration = OptimizeGeometry(
    bulk_configuration,
    max_forces=0.05*eV/Ang,
    max_steps=200,
    max_step_length=0.2*Ang,
    trajectory_filename='initial_C.nc',
    disable_stress=True,
    optimizer_method=LBFGS(),
)
nlsave('initial_C.nc', bulk_configuration)
nlprint(bulk_configuration)

# -------------------------------------------------------------
# Total Energy
# -------------------------------------------------------------
total_energy = TotalEnergy(bulk_configuration)
nlsave('initial_C.nc', total_energy)
nlprint(total_energy)