Monte Carlo and subsets of energy terms

[Justin Spiriti]

As part of my efforts to implement mixed-resolution modeling in OpenMM, I'm working on an efficient and general Monte Carlo simulator. So far I have coded the Metropolis algorithm and backbone dihedral moves in Python. My code is attached in "mc-test4.txt" (the forum wouldn't allow me to attach files with a .py extension) and profiling output obtained using cProfile is attached in mc-prof-output4b.txt.

I am assuming that the calls to "_openmm.Context__getStateAsLists" which appear to be taking up about 2/3 of the total time for the simulation represent energy calls. My code calculates the energy in full on every trial move. Is it possible to use OpenMM to calculate the energy for only a subset of the atoms, or for an arbitrary subset of energy terms?

Thank you very much for your help.

Justin Spiriti

[Peter Eastman]

There are a few things you can do, but I don't know if any of them would be useful. You can restrict the energy evaluation to only a subset of the Force objects in the System. To do that, put different Forces into different force groups by calling setForceGroup() on them. Then when you call getState(), use the "groups" argument to specify which force groups to include. Since the bond and angle energies aren't affected by a dihedral rotation, you could skip calculating them. But I doubt it would make a difference. Unless you're dealing with a tiny system, bonds and angles take a tiny amount of time compared to the nonbonded interactions.

Another possibility is that if you use a CustomNonbondedForce for the nonbonded interactions, you can use interaction groups to compute only a subset of the nonbonded interactions. But in your case, I don't think that's really practical. It's mainly useful when you have a specific small subset of interactions you want to compute repeatedly. But a dihedral rotation will change a large fraction of all nonbonded interactions, and that subset will be different for every move. Furthermore, it looks like you're using implicit solvent, which is a many particle (not just pairwise) interaction. So after every move, you really do need to recompute all GB interactions for all atoms.

Peter

[Mikael Lund]

I'm also interested in MC simulations where only a subset of energies are calculated. My idea is to use setInteractionGroupParameters() for each MC move, but so far I haven't found a way to do this dynamically during simulation (the context is not updated). Is there a reasonably inexpensive way to do so?

UPDATE: After changing the interaction group(s), the following pushes the changes back to the context:

Code: Select all

context.reinitialize()
context.setPositions(positions)

I'm not quite sure why the positions need updating, though(?)

Mikael

[Peter Eastman]

When you call reinitialize(), it throws out all the internal data structures that represent the context and rebuilds them from scratch. That means all the state information (such as positions) associated with the context is lost, and you need to set it again. Other things that get lost include velocities, periodic box vectors (if you've changed them from the default values set in the System), and context parameters.

As you'd guess, reinitialize() is a somewhat expensive operation. It's fine to do it every now and then, but you wouldn't want to do it for every MC move. It would be faster to just evaluate the entire energy.

Peter