# Define a benzene ring
cartesian_coordinates = [[ 0.0     ,  1.399795, 0.0],
                         [ 1.212253,  0.69976 , 0.0],
                         [ 1.212253, -0.69976 , 0.0],
                         [ 0.0     , -1.399795, 0.0],
                         [-1.212253, -0.69976 , 0.0],
                         [-1.212253,  0.69976 , 0.0],
                         [ 0.0     ,  2.500678, 0.0],
                         [ 2.165671,  1.250363, 0.0],
                         [ 2.165671, -1.250363, 0.0],
                         [ 0.0     , -2.500678, 0.0],
                         [-2.165671, -1.250363, 0.0],
                         [-2.165671,  1.250363, 0.0]]*Angstrom

molecule_configuration = MoleculeConfiguration(
    elements=[Carbon,]*6+[Hydrogen,]*6,
    cartesian_coordinates=cartesian_coordinates
    )

# Attach a calculator
molecule_configuration.setCalculator(BrennerCalculator())

initial_velocity = MaxwellBoltzmannDistribution(
        temperature=300.0*Kelvin
    )

length_of_simulation = 20*femtoSecond

# Generate trajectories for different time steps
for time_step in [0.1, 0.5, 1, 2]*femtoSecond:

    method = NVEVelocityVerlet(
        time_step=time_step,
        initial_velocity=initial_velocity
        )

    # Generate trajectory
    molecular_dynamics = MolecularDynamics(
        molecule_configuration,
        trajectory_filename='trajectory.nc',
        steps=int(length_of_simulation/time_step),
        log_interval=max(1,int(2*femtoSecond/time_step)),
        method=method
        )