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

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

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

# Define coordinates
fractional_coordinates = [[ 0.            ,  0.            ,  0.179593327651],
                          [ 0.5           ,  0.5           ,  0.232355209221],
                          [ 0.            ,  0.            ,  0.28496297267 ],
                          [ 0.5           ,  0.5           ,  0.338624592794],
                          [ 0.            ,  0.            ,  0.391088997647],
                          [ 0.            ,  0.            ,  0.458574404379],
                          [ 0.5           ,  0.5           ,  0.45770745749 ],
                          [ 0.            ,  0.            ,  0.512774191059],
                          [ 0.5           ,  0.5           ,  0.512505467165]]

# Define ghost atoms
ghost_atoms = [7, 8]

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

# Add tags
bulk_configuration.addTags('Left Interface',  [0, 1, 2, 3, 4])
bulk_configuration.addTags('Right Interface', [5, 6, 7, 8])
bulk_configuration.addTags('Selection 0',     [0, 1, 7, 8])

# -------------------------------------------------------------
# Calculator
# -------------------------------------------------------------
#----------------------------------------
# Basis Set
#----------------------------------------
basis_set = [
    GGABasis.Oxygen_DoubleZetaPolarized,
    GGABasis.Magnesium_DoubleZetaPolarized,
    GGABasis.Silver_SingleZetaPolarized,
    ]

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

numerical_accuracy_parameters = NumericalAccuracyParameters(
    k_point_sampling=(11, 11, 1),
    )

iteration_control_parameters = IterationControlParameters(
    tolerance=1e-05,
    )

poisson_solver = FastFourier2DSolver(
    boundary_conditions=[[PeriodicBoundaryCondition,PeriodicBoundaryCondition],
                         [PeriodicBoundaryCondition,PeriodicBoundaryCondition],
                         [NeumannBoundaryCondition,DirichletBoundaryCondition]]
    )

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

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

# -------------------------------------------------------------
# Optimize Geometry
# -------------------------------------------------------------
constraints = [0, 1, 7, 8]

bulk_configuration = OptimizeGeometry(
        bulk_configuration,
        max_forces=0.01*eV/Ang,
        max_steps=200,
        max_step_length=0.2*Ang,
        constraints=constraints,
        trajectory_filename='MgO1LAg.nc',
        disable_stress=True,
        optimizer_method=LBFGS(),
        )
nlsave('MgO1LAg.nc', bulk_configuration)
nlprint(bulk_configuration)

# -------------------------------------------------------------
# Chemical Potential
# -------------------------------------------------------------
chemical_potential = ChemicalPotential(bulk_configuration)
nlsave('MgO1LAg.nc', chemical_potential)
nlprint(chemical_potential)

# -------------------------------------------------------------
# Effective Potential
# -------------------------------------------------------------
effective_potential = EffectivePotential(
    configuration=bulk_configuration,
    )
nlsave('MgO1LAg.nc', effective_potential)