Periodicity of Bond Forces

[Gideon Simpson]

I've been playing with a CustomBondForce, and I've noticed that the pairs of atoms can drift outside of the periodic box, despite calling setUsesPeriodicBoundaryConditions(True). On the the other hand, the center of the bonds remains in the periodic domain. Is this the expected behavior? Is there a way to use bond forces that keeps everything in the periodic box?

[Peter Eastman]

Periodic boundary conditions affect how forces are computed. You can still place particles at any position in space, but for purposes of computing the displacements between particles, it uses the shortest possible distance between all periodic pairs of the two particles.

Independently of that, when you call getState() you can use the enforcePeriodicBox argument to tell it to return the positions of particles wrapped into a single periodic box. However, it will never break molecules into pieces. It places the centroid of each molecule inside a single periodic box, which means individual particles may sometimes still be outside it.

[Gideon Simpson]

I'm still a bit concerned about this. Consider the following two figures for my problem with a dimer joined by a CustomBondForce and the rest interacting through a nonbonded force. In both figures, there is the one blue atom outside of the periodic domain [0,10]X[0,10]. I would expect this atom to behave, through periodicity, as though it were in the upper right, within the simulation cell. I obtained these figures by dumping the data to Numpy with and without enforcePeriodicBox=True. This flag is definitely getting the nonbonded atoms, but the one blue atom is still outside the box.

To put this another way, I have two concerns:

• First, I'd just like to confirm that the periodicity is actually respected by the CustomBondForce. The only place I have indicated this is periodic is with a system.setDefaultPeriodicBoxVectors command and a force.setUsesPeriodicBoundaryConditions(True) command. Is that sufficient.
  
  Second, it would be nice to see some confirmation of this when the data is dumped for visualization. enforcePeriodicBox=True does not appear to work as expected.

positiosn1.pdf

(12.5 KiB) Downloaded 28 times

positiosn2.pdf

(12.44 KiB) Downloaded 26 times

[Peter Eastman]

Once again, there are two completely different issues here. It's important to understand what periodic boundary conditions mean and what effect they have.

When you are simulating a periodic system, every "particle" really represents an infinite grid of particles repeating throughout space. Suppose you have a cubic box of size 2. If a "particle" is at (1, 0, 0), that means there also are particles at (3, 0, 0), (-1, 0, 0), (3, 2, -2), etc. Likewise, if you specify that a particle has any one of those positions, that means there also is a particle at (1, 0, 0). It makes no difference which of those positions you specify. They all represent the same grid of particles.

Normally you apply periodic boundary conditions to nonbonded forces but not bonded ones. For the nonbonded interaction between two "particles" that means all possible periodic copies of the first one interact with all possible periodic copies of the second one. That's physically correct if, for example, you are computing the Coulomb interaction between charged particles. But it isn't appropriate for most bonded interactions. If two particles are bonded, they interact only with each other. One particle doesn't interact with the infinite periodic copies of the other, because it isn't bonded to them. This is why, by default, bonded forces don't apply periodic boundary conditions. There are some situations where you do want them to, but they're unusual.

Now let's talk about a different issue: writing coordinates to files. As far as nonbonded forces are concerned, it doesn't matter whether it writes (1, 0, 0) or (3, 0, 0). They both represent the same grid of repeating particles. But it makes visualization easier if it consistently reports coordinates that are in a single periodic box. So that's what the enforcePeriodicBox option is for. It tells it to move everything into a single box before writing them.

But remember, for most bonded forces it does matter which position it writes. If you split a molecule down the middle and reported one copy of one half, and a different copy of the other half, that would be both confusing for visualization and incorrect for computing forces. So it never splits a molecule in half. It always outputs a single periodic copy of the entire molecule. If the molecule happens to straddle an edge of the box, that means part of it will stick outside the box. That's a little confusing, but much better than breaking it in half.

What is a "molecule"? Simply a set of particles that are connected by bonds. In this case you have a CustomBondForce that connects the two atoms. Since there is a bond between them, they are treated as a single molecule and they will never get split apart. The midpoint between them is always in the periodic box, but that means one of them may stick outside it.

[Istvan Kolossvary]

I have a follow-up question regarding to call getState() with the enforcePeriodicBox argument. I have a simple reporter that writes out coordinates in (x,y,z) format (for some post processing) but I have a triclinic simulation cell and want to write out the coordinates enforcing PBC. The reporter looks like this:

Code: Select all

class PositionReporter(object):
    def __init__(self, file, reportInterval):
        self._out = open(file, 'w')
        self._reportInterval = reportInterval

    def __del__(self):
        self._out.close()

    def describeNextReport(self, simulation):
        steps = self._reportInterval - simulation.currentStep%self._reportInterval
        return (steps, True, False, False, False, None)

    def report(self, simulation, state):
        positions = state.getPositions().value_in_unit(nanometer)
        for p in positions:
            self._out.write('%g %g %g\n' % (p[0], p[1], p[2]))

What is not clear, how to set enforcePeriodicBox since I am not calling getState() directly. I tried

positions = simulation.context.getState(getPositions=True, enforcePeriodicBox=True).value_in_unit(nanometer)

but it gave me this error message.

Traceback (most recent call last):
File "run-openmm_diagnostics.py", line 181, in <module>
simulation.step(steps)
File ".../simtk/openmm/app/simulation.py", line 132, in step
self._simulate(endStep=self.currentStep+steps)
File ".../simtk/openmm/app/simulation.py", line 233, in _simulate
self._generate_reports(wrapped, True)
File ".../simtk/openmm/app/simulation.py", line 254, in _generate_reports
reporter.report(self, state)
File "run-openmm_diagnostics.py", line 55, in report
positions = simulation.context.getState(getPositions=True, enforcePeriodicBox=True).value_in_unit(nanometer)
File ".../simtk/openmm/openmm.py", line 14300, in <lambda>
__getattr__ = lambda self, name: _swig_getattr(self, State, name)
File ".../simtk/openmm/openmm.py", line 74, in _swig_getattr
return _swig_getattr_nondynamic(self, class_type, name, 0)
File ".../simtk/openmm/openmm.py", line 69, in _swig_getattr_nondynamic
return object.__getattr__(self, name)
AttributeError: type object 'object' has no attribute '__getattr__'

I probably overlooked something trivial. Could you tell me how to do this correctly?

Thanks,

Istvan

[Peter Eastman]

That should be

Code: Select all

positions = simulation.context.getState(getPositions=True, enforcePeriodicBox=True).getPositions().value_in_unit(nanometer)

You were missing the call to getPositions() on the State.

What is not clear, how to set enforcePeriodicBox since I am not calling getState() directly.

The last value returned by describeNextReport() specifies whether to use enforcePeriodicBox. You're returning None, which gets interpreted as True for a periodic system or False for a nonperiodic one. See http://docs.openmm.org/latest/userguide ... other-data.

[Istvan Kolossvary]

Thank you, Peter. I was wondering what the 5th return value of describeNextReport() was. It looks like, unknowingly, I was already writing out PBC coordinates.

[Xinyu Gu]

Dear Peter,

This clarification is very helpful to me. Many thanks!

By the way, I have another question about the Periodicity of Bond Forces. I have a CustomBondForce which is used for calculating a collective variable between two diffenent molecules. Then I used that collective variable as a reaction coordinate to run umbrella sampling simulations with PBC.

I'm going to use "setUsesPeriodicBoundaryConditions(true)" to turn on PBC for above CustomBondForce, so that we only consider the shortest "bond" given all periodic copies. However, I'm worried about the definition of "molecule" in my case. I'm wondering that whether the two molecules are considered as one "molecule" if there is a CustomBondForce(PBC) between them?

If they are considered as one "molecule", we can not split that big "molecule" and my CustomBondForce would consider the "bond" inside the big "molecule" instead of the shortest "bond" among all periodic copies. That could be a problem for me. Would you please help me figure out the definition of "molecule" in my case?

Thank you in advance for your time and helps, I sincerely appreciate it!


Best,
Xinyu

Once again, there are two completely different issues here. It's important to understand what periodic boundary conditions mean and what effect they have.

When you are simulating a periodic system, every "particle" really represents an infinite grid of particles repeating throughout space. Suppose you have a cubic box of size 2. If a "particle" is at (1, 0, 0), that means there also are particles at (3, 0, 0), (-1, 0, 0), (3, 2, -2), etc. Likewise, if you specify that a particle has any one of those positions, that means there also is a particle at (1, 0, 0). It makes no difference which of those positions you specify. They all represent the same grid of particles.

Normally you apply periodic boundary conditions to nonbonded forces but not bonded ones. For the nonbonded interaction between two "particles" that means all possible periodic copies of the first one interact with all possible periodic copies of the second one. That's physically correct if, for example, you are computing the Coulomb interaction between charged particles. But it isn't appropriate for most bonded interactions. If two particles are bonded, they interact only with each other. One particle doesn't interact with the infinite periodic copies of the other, because it isn't bonded to them. This is why, by default, bonded forces don't apply periodic boundary conditions. There are some situations where you do want them to, but they're unusual.

Now let's talk about a different issue: writing coordinates to files. As far as nonbonded forces are concerned, it doesn't matter whether it writes (1, 0, 0) or (3, 0, 0). They both represent the same grid of repeating particles. But it makes visualization easier if it consistently reports coordinates that are in a single periodic box. So that's what the enforcePeriodicBox option is for. It tells it to move everything into a single box before writing them.

But remember, for most bonded forces it does matter which position it writes. If you split a molecule down the middle and reported one copy of one half, and a different copy of the other half, that would be both confusing for visualization and incorrect for computing forces. So it never splits a molecule in half. It always outputs a single periodic copy of the entire molecule. If the molecule happens to straddle an edge of the box, that means part of it will stick outside the box. That's a little confusing, but much better than breaking it in half.

What is a "molecule"? Simply a set of particles that are connected by bonds. In this case you have a CustomBondForce that connects the two atoms. Since there is a bond between them, they are treated as a single molecule and they will never get split apart. The midpoint between them is always in the periodic box, but that means one of them may stick outside it.

[Peter Eastman]

If you add a bond then yes, it will consider them a single molecule. But that might not matter. If you call setUsesPeriodicBoundaryConditions(true) on a force, that particular force will always apply periodic boundary conditions. It will compute the interaction based on the shortest distance between periodic copies of the atoms, whether or not they're in the same "molecule". I think that's what you want?

Your output might look a little strange, since it will always keep those two molecules together as if they were a single larger molecule. But that's just for output. It doesn't affect how forces are computed.

[Xinyu Gu]

Dear Peter,

That makes sense to me! Thank you so much for your help!!

Best,
Xinyu