# Set up configuration
bulk_configuration = BulkConfiguration(
    bravais_lattice=FaceCenteredCubic(5.4306*Angstrom),
    elements=[Silicon, Silicon],
    fractional_coordinates=[[0.0, 0.0, 0.0], [0.25, 0.25, 0.25]]
    )

# Set calculator
calculator = TremoloXCalculator(parameters=Tersoff_Si_1988b())
bulk_configuration.setCalculator(calculator)

# Set up MD method
method = NPTBernettiBussi(
    time_step=1*femtoSecond,
    initial_velocity=MaxwellBoltzmannDistribution(temperature=300*Kelvin),
    # Thermostat settings
    reservoir_temperature=300*Kelvin,
    thermostat_timescale=100*femtoSecond,
    heating_rate=0*Kelvin/picoSecond,
    # Barostat settings
    reservoir_pressure=1.0*bar,
    barostat_timescale=1000.0*femtoSecond,
    compressibility=1.0e-4*bar**-1,
)

# Run MD simulation
md_trajectory = MolecularDynamics(
    bulk_configuration,
    constraints=[],
    trajectory_filename='trajectory.hdf5',
    steps=50,
    log_interval=10,
    method=method
)