I am calculating the water permeation events through water channels of a protein aquaporinZ.Simulation of 100ns (10,000 frames) is performed for analysis which is based on a rectangular water box with water thickness of 22.5 and an aqpZ tetramer (having 4 water channels) embedded in lipid bilayer membrane. VMD tool is used for analysis. I got a relevant script in tcl but it does not give desired results. I want to find out only those water molecules which pass through the protein, while the script is calculating all the water molecules whose z coordinates are changing with time (due to periodic boundary conditions of MD simulation). As this code focus on z coordinates of water molecules with a defined upper (23.5) and lower end (-18.5) limit, I want to add a middle limit (z coordinate value "2.5") to calculate only those molecules that reach from upper end to lower end by passing through the middle limit and vice versa. Relevant code is:
# Specify the upper and lower boundaries of the nanotube layer
set upperEnd 6.7
set lowerEnd -6.7
# Start counting permeation events after this number of frames
set skipFrame 500
puts "Computing permeation events... (please wait)"
set wat [atomselect top "name OH2"]
set segList [$wat get segname]
set ridList [$wat get resid]
set labelList {}
foreach foo $segList {
lappend labelList 0
}
set num1 0
set num2 0
set numFrame [molinfo top get numframes]
# Each water molecule has a label, which has 5 possible values
# 2: Above the nanotube layer
# -2: Below the nanotube layer
# 1: Inside the nanotube layer, entering from upper surface
# -1: Inside the nanotube layer, entering from lower surface
# 0: Inside the nanotube layer from the beginning
# For every frame, the label of each water molecule is
# determined, and compared with its label in the previous frame.
# If the new label is +2 (or -2), while the old label is -1 (or +1),
# it means the water molecule has traversed the nanotube, thus a
# permeation event is reported and counted. If a water molecule
# is inside the nanotube layer in the current frame, its label
# will be determined by its old label.
for {set fr 0} {$fr < $numFrame} {incr fr} {
molinfo top set frame $fr
set oldList $labelList
set labelList {}
foreach z [$wat get z] oldLab $oldList segname $segList resid $ridList {
if {$z > $upperEnd} {
set newLab 2
if {$oldLab == -1} {
puts "$segname:$resid permeated through the nanotubes along +z direction at frame $fr"
if {$fr >= $skipFrame} {
incr num1
}
}
} elseif {$z < $lowerEnd} {
set newLab -2
if {$oldLab == 1} {
puts "$segname:$resid permeated through the nanotubes along -z direction at frame $fr"
if {$fr >= $skipFrame} {
incr num2
}
}
} elseif {abs($oldLab) > 1} {
set newLab [expr $oldLab / 2]
} else {
set newLab $oldLab
}
lappend labelList $newLab
}
}
puts ""
set nf [expr $numFrame-$skipFrame]
if {$nf >= 0} {
puts "The total number of permeation events during $nf frames in +z direction is: $num1"
puts "The total number of permeation events during $nf frames in -z direction is: $num2"
} else {
puts "The specified first frame ($skipFrame) is larger than the total number of frames ($numFrame)"
}
Is this how it works: oldlist starts populated with 0's, but I'm assuming water molecules start outside the layer, right? A molecules label will be 0 until it appears above or below the limit. then it's label is set to 2/-2. It is assumed in the next frame that if it's no longer above or below that limit, that it is in the tube (elseif {abs($oldLab) > 1}). If then it appears at the other side above the limit, the event is counted. So it seems to already count just molecules that have passed from upper limit to lower limit and vice versa, or am I missing something here?
If you want to do it differently, could you append all z coordinates for a molecule to a string, and then use regex to check for a pattern within the coordinates that includes $lowerend, $middlelimit, $upperend or the other way around?
Ok this makes sense, I will give it a try. Thanks
cross posted: https://stackoverflow.com/questions/73568397/