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

# Set up lattice
lattice = FaceCenteredOrthorhombic(10.467*Angstrom, 12.87*Angstrom, 24.493*Angstrom)

# Define elements
elements = [Sulfur, Sulfur, Sulfur, Sulfur, Sulfur, Sulfur, Sulfur, Sulfur,
            Sulfur, Sulfur, Sulfur, Sulfur, Sulfur, Sulfur, Sulfur, Sulfur,
            Sulfur, Sulfur, Sulfur, Sulfur, Sulfur, Sulfur, Sulfur, Sulfur,
            Sulfur, Sulfur, Sulfur, Sulfur, Sulfur, Sulfur, Sulfur, Sulfur]

# Define coordinates
fractional_coordinates = [[ 0.0488,  0.8544,  0.8564],
                          [ 0.1436,  0.2596,  0.9512],
                          [ 0.8564,  0.7404,  0.0488],
                          [ 0.2596,  0.1436,  0.1456],
                          [ 0.7404,  0.8564,  0.8544],
                          [ 0.1456,  0.9512,  0.2596],
                          [ 0.8544,  0.0488,  0.7404],
                          [ 0.9512,  0.1456,  0.1436],
                          [ 0.322 ,  0.8306,  0.7382],
                          [ 0.2618,  0.3908,  0.678 ],
                          [ 0.7382,  0.6092,  0.322 ],
                          [ 0.3908,  0.2618,  0.1694],
                          [ 0.6092,  0.7382,  0.8306],
                          [ 0.1694,  0.678 ,  0.3908],
                          [ 0.8306,  0.322 ,  0.6092],
                          [ 0.678 ,  0.1694,  0.2618],
                          [ 0.2766,  0.7314,  0.6824],
                          [ 0.3176,  0.1904,  0.7234],
                          [ 0.6824,  0.8096,  0.2766],
                          [ 0.1904,  0.3176,  0.2686],
                          [ 0.8096,  0.6824,  0.7314],
                          [ 0.2686,  0.7234,  0.1904],
                          [ 0.7314,  0.2766,  0.8096],
                          [ 0.7234,  0.2686,  0.3176],
                          [ 0.2501,  0.0079,  0.5645],
                          [ 0.4355,  0.3225,  0.7499],
                          [ 0.5645,  0.6775,  0.2501],
                          [ 0.3225,  0.4355,  0.9921],
                          [ 0.6775,  0.5645,  0.0079],
                          [ 0.9921,  0.7499,  0.3225],
                          [ 0.0079,  0.2501,  0.6775],
                          [ 0.7499,  0.9921,  0.4355]]

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

# -------------------------------------------------------------
# Calculator
# -------------------------------------------------------------

potentialSet = ReaxFF_LiS_2015(strict_bondpairs = None)
calculator = TremoloXCalculator(parameters=potentialSet)
calculator.setVerletListsDelta(0.25*Angstrom)

bulk_configuration.setCalculator(calculator)
nlprint(bulk_configuration)
bulk_configuration.update()
nlsave('a-Sulfur.hdf5', bulk_configuration)

# -------------------------------------------------------------
# Optimize Geometry
# -------------------------------------------------------------
bulk_configuration = OptimizeGeometry(
    bulk_configuration,
    max_forces=0.001*eV/Ang,
    max_stress=0.01*GPa,
    max_steps=200,
    max_step_length=0.2*Ang,
    trajectory_filename='a-Sulfur.hdf5',
    optimizer_method=LBFGS(),
    constrain_bravais_lattice=True,
)
nlsave('a-Sulfur.hdf5', bulk_configuration)
nlprint(bulk_configuration)

# -------------------------------------------------------------
# Total Energy
# -------------------------------------------------------------
total_energy = TotalEnergy(bulk_configuration)
nlsave('a-Sulfur.hdf5', total_energy)
nlprint(total_energy)