Reflective wall to prevent water from penetrating the lipid bilayer

[Ramon Mendoza]

I am seeking to implement a

Code: Select all

 CustomExternalForce 

applied to water particles to prevent them from penetrating the membrane in the early stage of the simulation.

This is simple with a flat-bottom harmonic potential, however the issue that I find is with PBC. This excerpt from the

Code: Select all

 CustomExternalForce 

documentation has caused me a lot of nightmares,

Special care is needed in systems that use periodic boundary conditions. In that case, each particle really represents an infinite set of particles repeating through space. The variables x, y, and z contain the coordinates of one of those periodic copies, but there is no guarantee about which.

and equally horrific from the FAQ,

Periodic boundary conditions are applied while computing forces that depend on the displacements between particles.

For the case of the harmonic potential the particles position affects the force. For example, for a linear force solely on z with a general expression of

Code: Select all

 " f * z " 

the force will be constant and will be independent of the coordinates of whichever particle OpenMM decides to use, although this is not true with the harmonic energy expression. Therefore, I think the usual implementation with the global coordinate is dangerous. For example, the flat-bottom harmonic restraint can be implemented with the expression,

Code: Select all

 " select(a, a, b) * (k/2) * ( z - select(a, zUpper, zLower) )^2; a=step(z - zUpper); b=step(zLower - z) " 

where

Code: Select all

 zUpper 

and

Code: Select all

 zLower 

are the locations of the reflective walls. The difficulty here is with PBC. I have considered the use of

Code: Select all

 periodicdistance(0, 0, z, 0, 0, zUpper) 

and

Code: Select all

 periodicdistance(0, 0, z, 0, 0, zLower) 

, however this does not lead to a simple solution. There is also an increase of complexity here because

Code: Select all

 periodicdistance(x1, y1, z1, x2, y2, z2) 

seems to give the absolute distance, which for the case of an arbitrary plane at

Code: Select all

 z = zUpper 

will give the minimum distance of particles in both the above and below the plane with the loss of information on which side of the plane the particle was on.

Does anyone have suggestions to implement this restraint? It is also a possible that I have misunderstood the

Code: Select all

 CustomExternalForce 

documentation. In any case, any help will be appreciated.

[Peter Eastman]

Can you give the exact equation for what you want the potential to be as a function of z? If you can define it mathematically, it should be easy to translate into a custom force.

[Ramon Mendoza]

Sure, here is the expression of the flat-bottom potential

equation_edit.png (10.52 KiB) Viewed 244 times

The z-upper and z-lower bounds will be defined by the head groups of the lipids.

**Edit: Fixed typo in energy equation

[Peter Eastman]

That expression is not periodic. It just depends on the absolute z coordinate, and has no dependence on the size of the box. Are you sure that's what you want?

[Ramon Mendoza]

That seems to be what I need to be able to prevent water from entering the membrane region.

I can define the boundaries of the membrane region and then I just want to reflect all waters. But doesn't the customExternalForce use periodic coordinates? So, the direct implementation would not be correct.

[Peter Eastman]

There aren't really "periodic coordinates", just coordinates. But if you intend your system to be periodic, you need to consider the periodicity correctly. Otherwise you'll get unexpected results if the coordinate isn't in the range you expect.

Periodic boundary conditions are complicated and confusing! They seem like they should be simple until you try to write down an equation. The equation you gave is a simple sum of two harmonic terms. There's one potential if z is above zupper, and a different potential if z is below zlower. But with periodic boundary conditions, every point is both above and below both edges!

I think the following may be something like what you want. Let zcenter=(zupper+zlower)/2 be the center of the potential, and let zwidth=(zupper-zlower)/2 be half the width of the flat region. You could create a force as

Code: Select all

force = CustomExternalForce('0.5*k*r^2; r=max(0, zperiodic-zwidth); zperiodic=periodicdistance(0, 0, z, 0, 0, zcenter)')

That will work as long as the particles don't get too far outside the flat region. But be aware that there's a discontinuity in the forces where the two harmonic terms meet, half a box width away from zcenter. Bad things could happen if anything crosses it (e.g. different atoms in a molecule being strongly pulled in opposite directions).

[Ramon Mendoza]

Thanks! I also came to that realization. Your point really sums up well

every point is both above and below both edges

This is what I was trying to say in my first post. Well really thank you, I will give this a shot.