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

# Set up lattice
lattice = Rhombohedral(4.131*Angstrom, 54.167*Degrees)

# Define elements
elements = [Arsenic, Arsenic]

# Define coordinates
fractional_coordinates = [[ 0.226,  0.226,  0.226],
                          [ 0.774,  0.774,  0.774]]

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

# -------------------------------------------------------------
# Calculator
# -------------------------------------------------------------
#----------------------------------------
# Basis Set
#----------------------------------------
basis_set = [
    GGABasis.Arsenic_DoubleZetaPolarized,
    ]

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

k_point_sampling = MonkhorstPackGrid(
    na=15,
    nb=15,
    nc=15,
    )
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('bulk-As.nc', bulk_configuration)

# -------------------------------------------------------------
# Optimize Geometry
# -------------------------------------------------------------
bulk_configuration = OptimizeGeometry(
        bulk_configuration,
        max_forces=0.01*eV/Ang,
        max_stress=0.01*GPa,
        max_steps=200,
        max_step_length=0.2*Ang,
        trajectory_filename=None,
        optimizer_method=LBFGS(),
        constrain_bravais_lattice=True,
        )
nlsave('bulk-As.nc', bulk_configuration)
nlprint(bulk_configuration)

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