Advice implementing this restraint

[lewis martin]

Hi all,
I'm hoping for some advice on how to implement some restraints.

I have a protein with five-fold axial symmetry (ion channel with five subunits). I have a ligand that binds at the interface of two of any of the subunits. I'd like to have restraining walls that go through each subunit, restraining the ligand to a wedge-shape of the cell that includes just one interface. I would apply these walls such that five ligands are restrained to five of the interfaces in a single simulation.

Is there a way, in OpenMM, to restrain the angle around the z axis, for instance? This would be more elegant than xyz restraints in that it reduces the total number of restraints and thus has better performance. If not then I will look at PLUMED!
cheers for your thoughts
Lewis

[Peter Eastman]

You should be able to do that with a CustomExternalForce or CustomCompoundBondForce (if you want to apply it to individual atoms), or with a CustomCentroidBondForce (if you want to apply it to the centroid of the ligand). Suppose you want to restrain it to the range of angles +/- pi/5. You can compute the angle as theta=atan(y/x), then apply a flat bottomed harmonic force along the lines of k*(max(0, abs(theta)-pi/5))^2. You'll need to be careful of numerical stability. That expression will blow up if it gets too close to x=0. A more conventional restraint that just keeps it close to a particular region of space may seem less elegant, but it's more numerically stable.

[lewis martin]

cheers Peter!
Didn't know about the customcentroidbondforce. Ive implemented this and will watch and wait for explosions. Another possible restraint Ill try is a cylinder encompassing the interface, like k*(max(0, (x-x0)^2+(y-y0)^2-150)). Thanks!