ParallelTemperingSampler from OpenMM checkpoint file

[nope-sto]

I tried to run the example ParallelTemperingSampler script from a previously generated OpenMM checkpoint file. Unfortunately, I end up with following error: "Exception: All sampler states must have box_vectors defined if the system is periodic." How do I have to specify the positions to have box_vectors included? Thanks in advance!

simulation.loadCheckpoint('after_free_eq.chk')
state = simulation.context.getState(getPositions=True, getVelocities=True, enforcePeriodicBox=True)
positions = state.getPositions()

n_replicas = 3  # Number of temperature replicas.
T_min = 298.0 * unit.kelvin  # Minimum temperature.
T_max = 600.0 * unit.kelvin  # Maximum temperature.
reference_state = states.ThermodynamicState(system=simulation.system, temperature=T_min)
move = mcmc.GHMCMove(timestep=2.0*unit.femtoseconds, n_steps=50)
remd_simulation = mmt.multistate.ParallelTemperingSampler(mcmc_moves=move, number_of_iterations=2)
storage_path = tempfile.NamedTemporaryFile(delete=False).name + '.nc'
reporter = mmt.multistate.MultiStateReporter(storage_path, checkpoint_interval=10)

remd_simulation.create(reference_state,
                  states.SamplerState(positions),
                  reporter, min_temperature=T_min,
                  max_temperature=T_max, n_temperatures=n_replicas)
remd_simulation.run(n_iterations=1)

[amin-sagar]

I am also running into the same error. My script is as follows

modeller = Modeller(openmm_topology, openmm_positions)
modeller.addSolvent(forcefield, model='tip3p', padding=1*nanometer,boxShape='octahedron',ionicStrength=0.15*molar)
system = forcefield.createSystem(modeller.topology, nonbondedMethod=PME,nonbondedCutoff=1.0*nanometer, constraints=HBonds,hydrogenMass=4*amu)

temperature = 310 * unit.kelvin
collision_rate = 1.0 / unit.picoseconds
timestep = 4 * unit.femtoseconds
integrator = openmm.LangevinMiddleIntegrator(temperature, collision_rate, timestep)
simulation = Simulation(modeller.topology, system, integrator)
simulation.context.setPositions(modeller.positions)
simulation.minimizeEnergy(maxIterations=100)

NVTsteps=5000
NPTsteps=0
num_steps=NVTsteps+NPTsteps
positions = simulation.context.getState(getPositions=True).getPositions()
PDBFile.writeFile(simulation.topology, positions, open('Minimized_solvated.pdb', 'w'))
print('Minimization Done')
simulation.reporters.append(DCDReporter('output_solvated_min_NVT_SAMS.dcd', 1000))
simulation.reporters.append(StateDataReporter('output_solvated_min_NVT_SAMS.log', 1000, step=True,
        potentialEnergy=True, temperature=True, volume=True, progress=True, remainingTime=True, speed=True,totalSteps=num_steps))
print('Running NVT\n')
simulation.step(NVTsteps)

system.addForce(MonteCarloBarostat(1*bar, 300*kelvin))
simulation.context.reinitialize(preserveState=True)
positions = simulation.context.getState(getPositions=True).getPositions()

n_replicas = 3  # Number of temperature replicas.
T_min = 298.0 * unit.kelvin  # Minimum temperature.
T_max = 600.0 * unit.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)]
thermodynamic_states = [states.ThermodynamicState(system=system, temperature=T)
                        for T in temperatures]

move = mcmc.GHMCMove(timestep=2.0*unit.femtoseconds, n_steps=50)
simulation = SAMSSampler(mcmc_moves=move, number_of_iterations=100,
                         state_update_scheme='global-jump', locality=5,
                         update_stages='two-stage', flatness_criteria='logZ-flatness',
                         flatness_threshold=0.2, weight_update_method='rao-blackwellized',
                         adapt_target_probabilities=False)


storage_path = tempfile.NamedTemporaryFile(delete=False).name + '.nc'
reporter = multistate.MultiStateReporter(storage_path, checkpoint_interval=1)
simulation.create(thermodynamic_states=thermodynamic_states,
                  sampler_states=[states.SamplerState(positions)],
                  storage=reporter)

simulation.run()

The error that I get is

Warning: The openmmtools.multistate API is experimental and may change in future releases
Warning: The openmmtools.multistate API is experimental and may change in future releases
Traceback (most recent call last):
  File "/home/amin/Work/CD38/Bic/Parent/Parent_min_NVT_NPT_SAMS.py", line 84, in <module>
    simulation.create(thermodynamic_states=thermodynamic_states,
  File "/home/amin/anaconda3/envs/espaloma/lib/python3.9/site-packages/openmmtools/multistate/multistatesampler.py", line 535, in create
    self._pre_write_create(thermodynamic_states, sampler_states, storage,
  File "/home/amin/anaconda3/envs/espaloma/lib/python3.9/site-packages/openmmtools/multistate/sams.py", line 336, in _pre_write_create
    super()._pre_write_create(thermodynamic_states, sampler_states, storage=storage, **kwargs)
  File "/home/amin/anaconda3/envs/espaloma/lib/python3.9/site-packages/openmmtools/multistate/multistatesampler.py", line 767, in _pre_write_create
    raise Exception('All sampler states must have box_vectors defined if the system is periodic.')
Exception: All sampler states must have box_vectors defined if the system is periodic.
Exception ignored in: <function MultiStateSampler.__del__ at 0x7f1c2fedcb80>
Traceback (most recent call last):
  File "/home/amin/anaconda3/envs/espaloma/lib/python3.9/site-packages/openmmtools/multistate/multistatesampler.py", line 737, in __del__
TypeError: 'NoneType' object is not callable

I would be really grateful for any suggestions.

[amin-sagar]

I think I understand what was wrong. The test system is not periodic so the box vectors were not needed. For periodic systems,

simulation.create(thermodynamic_states=thermodynamic_states,
                  sampler_states=[states.SamplerState(positions,box_vectors=system.getDefaultPeriodicBoxVectors())],
                  storage=reporter)

this seems to work.