choderalab / openmmtools

A batteries-included toolkit for the GPU-accelerated OpenMM molecular simulation engine.
http://openmmtools.readthedocs.io
MIT License
244 stars 76 forks source link

[Questions] How to report to stdout while running ReplicaExchangeSampler #737

Open 4doom4 opened 1 month ago

4doom4 commented 1 month ago

I am using MultiStateReporter ("file.nc", checkpoint_interval=interval) as the storage in a simulation from ReplicaExchangeSampler. This works to have the states as netcdf once the simulation is done. Is there a way to attach multiple reports to the simulation like in native openmm to report to stdout information on the progress of the simulation? I saw https://openmmtools.readthedocs.io/en/stable/devtutorial.html#id18 but I am not sure how to use this.

Thanks.

ijpulidos commented 3 weeks ago

Some progress of the simulation should already be given in stdout by default. Can you share the code that you are using that's not reporting anything on stdout?

4doom4 commented 3 weeks ago

I basically followed the REST2 protocol from here: https://openmm.github.io/openmm-cookbook/latest/notebooks/tutorials/Running_a_REST_simulation.html

The relevant part of the script:

print("Configuring REST2 protocol...")
# Set temperatures for each thermodynamic state
n_replicas = 16  # Number of temperature replicas
T_min = 300 * kelvin  # Minimum temperature (i.e., temperature of desired distribution)
T_max = 600 * kelvin  # Maximum temperature
temperatures = [T_min + (T_max - T_min) * (math.exp(float(i) / float(n_replicas-1)) - 1.0) / (math.e - 1.0)
                for i in range(n_replicas)]

# Create reference thermodynamic state
rest_state = RESTState.from_system(rest_system)
thermostate = ThermodynamicState(rest_system, temperature=T_min)
compound_thermodynamic_state = CompoundThermodynamicState(thermostate, composable_states=[rest_state])

# Create thermodynamic states
sampler_state =  SamplerState(coord.positions, box_vectors=rest_system.getDefaultPeriodicBoxVectors())
beta_0 = 1/(kB*T_min)
thermodynamic_state_list = []
sampler_state_list = []
for temperature in temperatures:
    # Create a thermodynamic state with REST interactions scaled to the given temperature
    beta_m = 1/(kB*temperature)
    compound_thermodynamic_state_copy = copy.deepcopy(compound_thermodynamic_state)
    compound_thermodynamic_state_copy.set_rest_parameters(beta_m, beta_0)
    thermodynamic_state_list.append(compound_thermodynamic_state_copy)

    # Generate a sampler_state with minimized positions
    context, integrator = cache.global_context_cache.get_context(compound_thermodynamic_state_copy)
    sampler_state.apply_to_context(context, ignore_velocities=True)
    openmm.LocalEnergyMinimizer.minimize(context)
    sampler_state.update_from_context(context)
    sampler_state_list.append(copy.deepcopy(sampler_state))

# Set up sampler
swap_attempts = 500
move = mcmc.LangevinDynamicsMove(timestep=0.002*picoseconds, n_steps=1000) # each move is 2 ps
simulation = ReplicaExchangeSampler(mcmc_moves=move, number_of_iterations=swap_attempts)

# Run repex
start, start_time = timer()
print(f"Start REST2 protocol at {start_time}...")
reporter = multistate.MultiStateReporter(f"{prefix}_rest2.nc", checkpoint_interval=5)

simulation.create(thermodynamic_states=thermodynamic_state_list,
                  sampler_states=sampler_state_list,
                  storage=reporter)

simulation.run()

end, end_time = timer()
print(f"REST2 protocol completed at {end_time}. Elapsed time: {end - start:.2f} seconds")