Collision Detection

[Yan Zhang]

When the distance between two particles is smaller than certain value, could openmm detect this event and allow me to run some codes before the next simulation step?

[Peter Eastman]

What sort of code do you want to run when that happens?

[Yan Zhang]

The codes could be removing and adding some external forces.

[Peter Eastman]

If you mean you want to actually add and remove Force objects in the System, that requires reinitializing the Context. There's no way to do that from inside of a running simulation. So what you would have to do is take one step by calling step(1), then check the positions yourself and take appropriate action before taking the next step.

But do you really need to do that? Adding and removing forces is rarely the best solution to a problem. What's the actual behavior you want to implement?

[Yan Zhang]

I want to simulate the effect of topoisomerase which allows a DNA segment to pass through another one. Because the enzyme has an orientation and is fixed on the DNA, so the passing-through effect is also directional. In the last post I wanted to use a forcefield to simulate this but I find it too difficult for me to implement within a week. So now I came up with this idea that once a collision between a DNA bead and an enzyme bead is found I will remove the LJ of the latter and apply a force on the DNA to drag it through the crack.

Does that make sense? It's an ongoing project and our experimental collaborators don't want us to leak too much information. I hope that was enough information for you to understand my dilemma

[Yan Zhang]

Picture1.png (93.29 KiB) Viewed 497 times

This is a diagram of what I wanted to implement. I hope it is informative enough

[Peter Eastman]

Thanks, that gives me a much better idea of what you're trying to accomplish. I assume this is a coarse grained model, not an atomic detail one? What level of detail does it use? One particle per base? What's the total number of particles in the system?

Another important question is how you want to model the behavior of the enzyme. For example, can any base on one strand pass through any base on the other? Or is there one particular base on one strand that can pass through one particular base on the other? Are you treating the behavior up until that point as random, so nothing happens until the correct bases just happen to bump into each other? Or are you specifically pulling them together?

In any case, rather than adding and removing force objects (expensive, requires reinitializing the context), you can just modify the parameters of the particles in question. Modify them in the Force object, then call updateParametersInContext() on it. That's a lot more efficient.

[Yan Zhang]

Thanks, that gives me a much better idea of what you're trying to accomplish. I assume this is a coarse grained model, not an atomic detail one? What level of detail does it use? One particle per base? What's the total number of particles in the system?

Another important question is how you want to model the behavior of the enzyme. For example, can any base on one strand pass through any base on the other? Or is there one particular base on one strand that can pass through one particular base on the other? Are you treating the behavior up until that point as random, so nothing happens until the correct bases just happen to bump into each other? Or are you specifically pulling them together?

In any case, rather than adding and removing force objects (expensive, requires reinitializing the context), you can just modify the parameters of the particles in question. Modify them in the Force object, then call updateParametersInContext() on it. That's a lot more efficient.

Yes, this is a coarse grained model, where one particle represents a 10nm segment on a chromosome (a few kilo-basepairs as I remember). The total number of particles will be around a thousand if not more.

There are two regions on this chromosome fiber, as shown in the diagram above, the helix region (blue) and the loop region (green). The helix region is where the enzymes are located, and thus only a bead in the loop region can pass through a bead in the helix region. So far I don't have any attraction in the system, so yes, nothing happens untill they bump into each other. I checked that around a few thousand steps, such collision event will likely happen.

The easy part about my simulation is that, the helix region is treated as a static object, so the coordinates of all the helix beads are constant. Now the expression of my force F (in the diagram) is:

Code: Select all

force_enz = CustomExternalForce('kt*1*r; r=sqrt((x-x0)*(x-x0) + (y-y0)*(y-y0) + (z-z0)*(z-z0))')
ikt = force_enz .addPerParticleParameter("kt")
ix = force_enz .addPerParticleParameter("x0")
iy = force_enz .addPerParticleParameter("y0")
iz = force_enz .addPerParticleParameter("z0")
system.addForce(force_enz )

As you can see the force is a constant force dragging the particles towards (x0, y0, z0). Now after reading the API, I'm still not sure how to set the values for per-particle parameters. E.g. if I want to add this force to particle i, should I do:

Code: Select all

force_enz.addParticle(i, [kt, x0, y0, z0])

Also if I want to set "kt" for atom j to be 10.0 during the simulation, should I do:

Code: Select all

force_enz.setParticleParameters(ikt, i, 10.0)

?
The API says:

Code: Select all

setParticleParameters(self, index, particle, parameters)

But I'm not sure what each parameter exactly means here.

Thanks!

[Peter Eastman]

CustomExternalForce can be applied to an arbitrary subset of particles, so that's why there's two different indices. For example, if you apply it to particles 10, 50, and 321, you would call addParticle() three times with those indices to create three "particle terms". When you call setParticleParameters(), the first argument says which of those particle terms you're modifying, so it will be 0, 1, or 2. The second argument says which particle that term is applied to, so it will be 10, 50, or 321. The third argument is the full list of per-particle parameter values, just as you passed it to addParticle().

[Yan Zhang]

CustomExternalForce can be applied to an arbitrary subset of particles, so that's why there's two different indices. For example, if you apply it to particles 10, 50, and 321, you would call addParticle() three times with those indices to create three "particle terms". When you call setParticleParameters(), the first argument says which of those particle terms you're modifying, so it will be 0, 1, or 2. The second argument says which particle that term is applied to, so it will be 10, 50, or 321. The third argument is the full list of per-particle parameter values, just as you passed it to addParticle().

I'm not sure I understand what "particle terms" exactly means. I understand that it is a separate indexing for particle 10, 50, and 321, but then in setPerParticleParameter why do we need two indices instead of just one? E.g.:

Code: Select all

force.setPerParticleParameter(50, 1, some_value)