MetropolisMonteCarloIntegrator does not work when I apply steps over time.

[BlackPianoCat]

I am trying to create a Monte Carlo simulation that I have some structure and I apply Gaussian displacements over time. I thought that MetropolisMonteCarloIntegrator is appropriate class to do something like that. I have a code like the following one,

pdb = PDBxFile(self.save_path+'MultiEM_init.cif') if init_struct_path==None or build_init_mmcif else PDBxFile(init_struct_path)
self.mass_center = np.average(get_coordinates_mm(pdb.positions),axis=0)
forcefield = ForceField('forcefields/ff.xml')
self.system = forcefield.createSystem(pdb.topology)
integrator = MetropolisMonteCarloIntegrator(Temperature,0.2,5)
# integrator = mm.LangevinIntegrator(Temperature, 0.05, 10 * mm.unit.femtosecond)

 # Import forcefield
print('\nImporting forcefield...')
self.add_forcefield()
print('---Done!---')

# Run simulation / Energy minimization
print('\nEnergy minimization...')
platform = mm.Platform.getPlatformByName(pltf)
simulation = Simulation(pdb.topology, self.system, integrator, platform)
simulation.reporters.append(StateDataReporter(stdout, (MD_steps*sim_step)//20, step=True, totalEnergy=True, potentialEnergy=True, temperature=True))
simulation.reporters.append(DCDReporter(self.save_path+'MultiEM_annealing.dcd', 5))
simulation.context.setPositions(pdb.positions)
simulation.context.setVelocitiesToTemperature(Temperature, 0)
current_platform = simulation.context.getPlatform()
print(f"Simulation will run on platform: {current_platform.getName()}.")
start_time = time.time()
simulation.minimizeEnergy()
state = simulation.context.getState(getPositions=True)
PDBxFile.writeFile(pdb.topology, state.getPositions(), open(self.save_path+'MultiEM_minimized.cif', 'w'))
print(f"--- Energy minimization done!! Executed in {(time.time() - start_time)/60:.2f} minutes. :D ---\n")

Until here everything is ok. Energy is minimized. And now I do simulated annealing, and I use it as I use any other integrator in OpenMM,

print('Running Simulated Annealing...')
start = time.time()
for i in range(MD_steps):
    T_i  = Temperature-i/MD_steps*(Temperature-Temp_f)
    simulation.integrator.setTemperature(T_i)
    simulation.step(sim_step)
    if (i*sim_step)%10==0:
        self.state = simulation.context.getState(getPositions=True)
        PDBxFile.writeFile(pdb.topology, self.state.getPositions(), open(self.save_path+f'ensembles/ens_{i//100+1}.cif', 'w'))
end = time.time()
elapsed = end - start
state = simulation.context.getState(getPositions=True)
PDBxFile.writeFile(pdb.topology, state.getPositions(), open(self.save_path+'MultiEM_afterMD.cif', 'w'))
print(f'Everything is done! Simulation finished succesfully!\nMD finished in {elapsed/60:.2f} minutes.\n')

And here the structure remains unchanged and the potential energy remains the same. what is wrong? Is there any simple way to integrate like Langevin integrator? Or should I use these samplers that you have in documentation, instead of Simulation class? I would like to have both .dcd trajectory and an ensemble of structures if it is possible.