Writing out coordinates that respect PBCs with MCMC Sampler

[aclyde11]

Hello,

I have a naive question...

I want to be able to pull a PDB File out every few steps or so (for some clustering/checking forces on ligand/etc).

For instance, when using a standard OpenMM simulation:

coords =  self.simulation.context.getState(getPositions=True,enforcePeriodicBox=True).getPositions(asNumpy=True)
app.PDBFile.writeFile(self.topology, coords, file=f)

but when I try it with MCMCSampler I can't get the periodic boundary conditions enforced in the positions (the PDB is correctly writing out the crystal header, but protein/ligand go over the edge).

I know I can get the coordinates out by

sampler = MCMCSampler(...)
sampler.run(10)
coords =  sampler.sampler_state.positions
# or 
coords = cache.global_context_cache.get_context(self.thermodynamic_state[0].getState(getPositions=True, enforcePeriodicBox=True).getPositions()

but in either case it doesn't work... except after minimization the second method works, but not the first. Once run/step is called, it doesn't seem to work anymore. I've checked all the objects, everyone knows about the box vector and states it periodic...

I would guess I'm just doing something I shouldn't be or not grabbing the context correctly.

Thanks!

[aclyde11]

It seems from a bit of digging there is an immediate solution with mdtraj #297 , is this still the best way to handle this use case?

[jchodera]

Tagging @andrrizzi who may have some insight as well!

There are two different operations that may be of interest:

1. Wrapping the centers of mass of the molecules into the box (which may still cause the protein to protrude from the box)
2. Imaging the molecules so that the macromolecule is in the central box and any larger non-solvent ligands are in the same periodic image, so that the protein and ligand are not split into separate images

Do you just need (1), or do you need (2) as well?

MDTraj has a method for (2) with Trajectory.image_molecules(), and as you point out, (1) should simply require getting the state---though I'm curious why that does not work.

Do you have a test case you can upload?
And would sampler_state.to_mdtraj() (feature request #297) work from you if we quickly implemented this?

[aclyde11]

@jchodera just 2 is useful. I was able to get going with my work using a link to the code from #297. So no need to rush, but in terms of features/documentation I found it pretty difficult to figure out using MCMC samplers what is the best way to pull out that information.

I think adding #297 would make this package instantly more useable--but again this isn't blocking on my end right now since manually doing #297 was pretty easy.

Thanks!