# -------------------------------------------------------------
# Bulk configuration
# -------------------------------------------------------------

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

# Define elements
elements = [Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon,
            Silicon]

# Define coordinates
fractional_coordinates = [[0.0 ,  0.0 ,  0.0 ],
                          [0.25,  0.25,  0.25],
                          [0.5 ,  0.5 ,  0.0 ],
                          [0.75,  0.75,  0.25],
                          [0.5 ,  0.0 ,  0.5 ],
                          [0.75,  0.25,  0.75],
                          [0.0 ,  0.5 ,  0.5 ],
                          [0.25,  0.75,  0.75]]

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

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

potentialSet = StillingerWeber_Si_1985()
calculator = TremoloXCalculator(parameters=potentialSet)

bulk_configuration.setCalculator(calculator)
bulk_configuration.update()

# -------------------------------------------------------------
# Molecular Dynamics
# -------------------------------------------------------------

initial_velocity = MaxwellBoltzmannDistribution(
    temperature=300.0*Kelvin,
    remove_center_of_mass_momentum=True
)

method = NVEVelocityVerlet(
    time_step=1*femtoSecond,
    initial_velocity=initial_velocity
)

md_trajectory = MolecularDynamics(
    bulk_configuration,
    trajectory_filename='trajectory.nc',
    steps=500,
    log_interval=50,
    method=method
)

bulk_configuration = md_trajectory.lastImage()

# Save the final configuration
nlsave('final_configuration.nc', bulk_configuration)

# Get the list of the kinetic energies of the snapshots
kinetic_energies = md_trajectory.kineticEnergies()

# Get the list of the potential energies of the snapshots
potential_energies = md_trajectory.potentialEnergies()