## Periodicity of Bond Forces

The functionality of OpenMM will (eventually) include everything that one would need to run modern molecular simulation.
Gideon Simpson
Posts: 21
Joined: Tue Sep 26, 2017 11:15 pm

### Periodicity of Bond Forces

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
Posts: 1912
Joined: Thu Aug 09, 2007 1:25 pm

### Re: Periodicity of Bond Forces

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
Posts: 21
Joined: Tue Sep 26, 2017 11:15 pm

### Re: Periodicity of Bond Forces

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
positiosn2.pdf
Last edited by Gideon Simpson on Wed Jul 17, 2019 8:02 am, edited 1 time in total.

Peter Eastman
Posts: 1912
Joined: Thu Aug 09, 2007 1:25 pm

### Re: Periodicity of Bond Forces

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
Posts: 28
Joined: Fri Jul 20, 2018 1:48 pm

### Re: Periodicity of Bond Forces

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
Posts: 1912
Joined: Thu Aug 09, 2007 1:25 pm

### Re: Periodicity of Bond Forces

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
Posts: 28
Joined: Fri Jul 20, 2018 1:48 pm

### Re: Periodicity of Bond Forces

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.