restrain orientation w/out affecting internal d.o.f.

[Daniel Konstantinovsky]

Good morning OpenMM forum!

I want to keep a DNA molecule pointed roughly along the z-axis throughout a long simulation. So far I've tried 1) restraining all DNA atoms with a very stiff spring constant and 2) introducing an energetic penalty to the DNA tilting and having x- and/or y- coordinates significantly bigger than zero. I want to ask if there is a better and cleaner way. What would you recommend I do to keep the molecule overall pointed in +z but allowing internal degrees of freedom to be roughly unperturbed so I can get dynamics within the DNA? Whatever the best way is, how can I use OpenMM to do it?

Thank you!

[Peter Eastman]

I can think of several ways of doing this.

As you suggested, you could use a CustomExternalForce to restrain all atoms to remain along the z axis with an energy expression like "k*(x^2+y^2)". That will tend to squash the molecule though. It would be better to use a flat bottomed potential that's zero as long as an atom remains inside a tube whose width is slightly larger than the size of the molecule: "k*(max(0, x-width)^2 + max(0, y-width)^2)".

Instead of applying this to every atom, you could just pick two atoms, one near each end, and only restrain them to stay near the axis. That pins the ends while leaving the rest of the molecule free.

You might also get less deformation of the molecule if you instead use a CustomCentroidBondForce. You would pick a cluster of atoms at each end of the molecule, and restrain their centroids to stay near the axis. That would spread out the restraining force over more atoms.

[Daniel Konstantinovsky]

Thank you for always answering so quickly. I will try the constraints first.

[Daniel Konstantinovsky]

Actually, I have one followup question. Say I wanted to penalize the molecule moving out of a "vertical tube" with walls defined by

force = CustomExternalForce("(step(x-x0) + step(-x0-x))*(k/2)*(x-x0)^2 + (step(y-y0) + step(-y0-y))*(k/2)*(y-y0)^2")
force.addGlobalParameter("k", 10 * kilocalories_per_mole / angstroms**2)
force.addPerParticleParameter("x0", 1.5 * nanometer)
force.addPerParticleParameter("y0", 1.5 * nanometer)

but in a simulation with PBC. periodicdistance() requires x, y, and z parameters. Is there a way I can make the above potential along x- and y- PBC-aware without including distances along z, since they don't matter?

force.addPerParticleParameter("z0", ???? * nanometer)

[Peter Eastman]

You can pass whatever arguments you want to periodicdistance. For example, periodicdistance(x1, y1, 0, x2, y2, 0) will give the distance between two points in the xy plane, ignoring z.

[Daniel Konstantinovsky]

Got it. Thank you!