PBCs and translation of molecules

[Charles Brooks]

I am seeing an error that I want to attribute to atoms/molecules being translated "into" a periodic box as part of the PBCs but am unsure exactly what is implemented in OpenMM (looked through the code and it appears it doesn't translate atoms back into box). Can someone let me know how this works in OpenMM?

My observation is the following: I have a layer of atoms fixed, i.e., their masses are set to 0 so that the the eoms for these atoms are not computed, I have long-chain alkane-like molecules with harmonic tethers adjacent to these atoms and I am occasionally getting both very large values for the harmonic restraint energy and NaN values for the bond energy, which I am assuming arise from the long-chain hydrocarbon having one of it's atoms "translated" across the periodic box and the restraint therefore being "stretched" as well as the bonds adjacent to the atom that is translated being stretched.

Ideas and suggestions welcome.

Charles Brooks

[Lee-Ping Wang]

Hi Charles,

I had issues like this when running simulations in TINKER 6.1 a couple of years ago. My solution was to apply the minimum-image convention to the displacement vector between the atom and its tethered site. It was a small change to the TINKER code and it's in the latest version (6.3).

It would be nice if OpenMM also supported this, since we never intend for distance vectors in our CustomForces to change by a periodic box vector. I'm not sure which internal variables OpenMM is using to actually compute this.

Thanks,

- Lee-Ping

[Charles Brooks]

Thanks for the comments. We already do that. I believe it may be the bond that is leading the problem here and the restraint energy is just a consequence. I am just not sure how OpenMM deals with the bonds for more extended systems, e.g., long hydrocarbon chains, when there are periodic boundary conditions.

Charlie

[Jason Swails]

As far as I can tell, OpenMM does not do any internal "imaging" on the internal positions of its Context. When you call "getState" on the context, you can specify "enforcePeriodicBox" which has the effect of wrapping all particles back into the central box, but does so on a copy of the positions rather than the actual position array.

It does this by molecule, which are defined in OMM as a closed topological set with atoms connected by bonds (only the CustomBondForce and HarmonicBondForce classes actually define Bonds) and/or constraints. I also don't see anywhere that the Integrators wrap molecules either -- all of the "imaging" will be done in the pairlist building, but that routine doesn't seem to actually change any of the positions in the Context either (it just uses the MIC to build the pairlist). But this pairlist is only used for Forces that implement a periodic potential (e.g., NonbondedForce, CustomNonbondedForce, GBSAOBCForce, and CustomGBForce, and AmoebaMultipoleForce).

I believe if you want to see exactly what positions are being used internally by OpenMM, you should set "enforcePeriodicBox=false" in the call to Context.getState, I believe. But at the end of the day, as long as you've actually defined a "bond" between two atoms, those atoms will be considered part of the same molecule and will be both scaled together during the MC Barostat moves and imaged together when reporting "wrapped" coordinates via calls to getState.

One thing I thought of, though, is if your restraint "bond" length is greater than the length of one of your periodic cell dimensions, it's possible that the nonbonded forces that enforce periodicity could use a closer image than the one used by the CustomBondForce or HarmonicBondForce -- assuming, of course, that the two atoms connected by that "restraint" are not in each others' exclusion lists. In this case, the solution is probably to just use a bigger box so this can't happen.

All the best,
Jason

[Peter Eastman]

If you have a barostat on the system, then each time it attempts a step it will first translate all the molecules into a single periodic box. So that's probably what's happening. As Jason described, this is done by molecules, so a molecule should never get split into pieces by it.

How are your restraints implemented? Is this with a CustomExternalForce or a plugin?

Peter

[Charles Brooks]

Hi Peter,

This info has been helpful regarding how PBCs are implemented in OpenMM. No, I don't have a barostat on the system. However, as I noted earlier, I do have a layer of "fixed" atoms, i.e., atoms whose mass I set to zero. This represents a layer of gold to which the hydrocarbon chains are anchored, as represented by absolute harmonic restraints. this is implemented as a custom external force of form:

<Force energy="k * (px^2 + py^2 + pz^2); px=min(dx, bx-dx); py=min(dy, by-dy); pz=min(dz, bz-dz); dx=abs(x-x0); dy=abs(y-y0); dz=abs(z-z0)" forceGroup="6" type="CustomExternalForce" version="1">

where bx, by, bz are the box dimensions.

Also, I tested and this large change (NaN occurrence in the bond energy and and very large value for the energy of the restraint) occurs even if I "turn-off" the pbc treatment in the above force, and this puzzles me significantly.

Any further insights would be great.

Thanks,

Charlie

[Peter Eastman]

Ok, so this (probably) doesn't actually relate to periodic boundary conditions. Can you determine what the cause of the large energies is? For example, the positions of the massless particles should never change. So query their positions and verify that they're still where you put them. You also can loop over the atoms you've applied the restraint to and calculate the distance of each one from its anchor. Is there one whose position has somehow gotten messed up?

Another thought: how quickly does this happen once you start the simulation? I'm wondering whether it would be reasonable to record the state at every time step so you can identify the first point at which something goes wrong. Does the problem appear suddenly in a single time step, or do the forces and energies grow over multiple time steps? If the latter, which Force does the problem first appear in?

Peter