# Set up MD method with two thermostats on tagged groups of atoms.
method = NVTBerendsen(
    time_step=1*femtoSecond,
    reservoir_temperature=[('region1', 300*Kelvin), ('region2', 600*Kelvin)],
    thermostat_timescale=100*femtoSecond,
    heating_rate=0*Kelvin/picoSecond,
    initial_velocity=None
)