If your output frequency is sufficiently small, you should be able to
accomplish this by running a calculation over the z coordinate of each
molecule (assuming the membrane normal is parallel to the z axis.
Consider a single water molecule.
There are 3 regions it can be located in.
(region label), location
(-1), Z > pore upper bound
(0), pore lower bound < Z < pore upper bound
(1), Z < pore lower bound
In order for this to work your output frequency must be sufficiently short
such that a water molecule cannot
1) move from A up through C (on the next periodic cell above) and into B
from below, or vice versa (go from C down trhough A and into B) in a single
frame
2) travel across the channel (B) in a single frame
Or to put it another way, no molecule may not traverse an entire region in
a single frame.
So first, use an array to track what region a molecule is currently in at
each frame (i.e. -1, 0, or 1), call it X0
Then use a second array to track what region the molecule was in prior to
the region it is currently in, call it XP
Then create a third array to track what region the molecule will be in on
the next frame, call it X1
next for each frame (i), compare XP[i], X0[i], and X1[i].
If XP[i]==1, X0[i]==0, and X1[i]==-1, then a negative flux occurred
(molecule traversed channel going down)
If XP[i]==-1, X0[i]==0, and X1[i]==1, then a positive flux occurred
(molecule traversed channel going up)
note that XP may be tricky to construct since it wants the region that the
molecule was in before it got to the region it is currently in at each
frame.
So while X1 is just an offset of X0, XP must be constructed with more care
Below is a python function that can do this.
Using pytraj, you would have to take the z coordinates of each water
molecule's trajectory, then
use that to create the appropriate 'X0' and feed them to the function. This
will return a pair of values
corresponding to the number of times that molecule traversed the channel in
the positive and negative directions respectively
The function can also return the array structures it uses to count the
number of channel traversals.
def getFluxCounts(X0,giveFluxVector=False,giveTrackingArray=False):
X1=X0[1:] #next location
#XP - Tracks previous location
#X0 - Tracks current location
#X1 - Tracks next location
#Need to construct xp specially since it only changes
#when X0 changes to a new location
xp=X0[0]
XP=[]
#construct XP array
for ix,x0 in enumerate(X0[1:]):
#check if we are in a new region and update xp if so
if not (X0[ix+1]==X0[ix]):
xp=X0[ix]
XP.append(xp)
#Since we defined XP specially to always track the location that
#the molecule was in prior to being in the region it is currently
#in, we can define flux in terms of 2 cases
# +1: positive flux across channel when XP[i]=-1, X0[i]=0, X1[i]=1
# -1: negative flux across channel when XP[i]=1, X0[i]=0, X1[i]=-1
# and 0 otherwise
# in the interest of speed this can be accomplished using boolean
binary operators
# and applied using the 'map' command. This should be a bit faster than
using
# using a for loop (which has more overhead).
XX=np.array([XP[:-1],X0[1:-1],X1[1:]]).T
Xflux=map(lambda x: (x[2]-x[0])/2*((x[0]+x[2])==0)*(x[1]==0),XX)
Xflux=np.concatenate([[0],Xflux,[0]])
if not (giveFluxVector or giveTrackingArray):
return( (np.sum(Xflux>0)*1.,np.sum(Xflux<0)*1.) )
else:
outDict={}
if giveFluxVector:
outDict['FluxVector']=Xflux
if giveTrackingArray:
outDict['TrackingArray']=XX
return( (np.sum(Xflux>0)*1.,np.sum(Xflux<0)*1.,outDict) )
Hope that helps
On Wed, Jun 12, 2019 at 9:31 AM Chetna Tyagi <cheta231.gmail.com> wrote:
> Dear Amber users,
>
> I would like to calculate the number of water molecules that crossed
> through the pore as a function of simulation time from one end to another.
>
> The hexameric pore is set in a DOPC:DOPG bilayer membrane. I have tried
> 'diffusion' calculation but it doesn't calculate the number of water
> molecules passing through the pore.
>
> Any ideas will help me a lot. Thanks.
>
> --
> Best wishes
> Chetna
> _______________________________________________
> AMBER mailing list
> AMBER.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber
>
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Wed Jun 12 2019 - 13:30:02 PDT
