From: Olya Kravchenko (ovkrav_at_gmail.com)
Date: Thu Jan 12 2017 - 10:39:00 CST
Hi Azadeh,
I recently asked similar question here about restricting the flow of ions
through the pore in a convex surface of a spherical protein, you can try to
find this discussion. I ended up setting up an "exclusion sphere" with a
centre near the middle of the channel. The ions that enter the sphere are
pushed out except for the one I want to stay inside.
I used the "User-Defined Forces in NAMD" tutorial:
http://www.ks.uiuc.edu/Training/Tutorials/science/forces/forces-tutorial-html/forces-tutorial.html
particularly this section:
http://www.ks.uiuc.edu/Training/Tutorials/science/forces/forces-tutorial-html/node3.html
and I based my script on sample scripts that come along with it:
http://www.ks.uiuc.edu/Training/Tutorials/science/forces/forces-tutorial-files/tclBCfiles/
I noticed that it significantly slows down my calculations. If you come up
with a more efficient way of removing the ions, please share.
Here is the script that works in my system:
set ionList {} ;# a list of atom IDs of the ions
set atomID 1
set inStream [open $pdbSource r] ;# open PDB file
#set sphereCenter {34.0 98.4 92.3} ;# this sets can be set in the
configuration file
# find all sodium and chloride ions in the pdb file and add their IDs into
the list
foreach line [split [read $inStream] \n] {
set string1 [string range $line 0 3]
set string2 [string range $line 17 20]
set string2 [string trim $string2]
if { [string equal $string1 {ATOM}] || \
[string equal $string1 {HETA}] } {
if { [string equal $string2 {CLA}] || \
[string equal $string2 {SOD}] } {
lappend ionList $atomID
}
incr atomID
}
}
close $inStream
wrapmode cell
proc calcforces {step unique} {
global atomID sphereCenter sphereRadius maxForce ionList
foreach { x0 y0 z0 } $sphereCenter { break }
while {[nextatom]} {
set atomid [getid]
# check whether atom ID is on the list, then check it's position
if { [lsearch $ionList $atomid] >= 0 } {
set rvec [getcoord]
foreach { x y z } $rvec { break }
# find the distance between the ion and the sphere's center
set rho [expr sqrt(($x-$x0)*($x-$x0) + ($y-$y0)*($y-$y0) + \
($z-$z0)*($z-$z0))]
if { [getid] == 344742 || [getid] < 344526 } {
dropatom
continue
} else {
set forceX [expr $maxForce * ($x-$x0) / $rho]
set forceY [expr $maxForce * ($y-$y0) / $rho]
set forceZ [expr $maxForce * ($z-$z0) / $rho]
}
addforce "$forceX $forceY $forceZ"
}
}
}
In my case all ions are listed in the end of the pdb file, hence all IDs
that are less than 344526 are protein atoms that I want to be left
untouched. The force is applied to the rest of the ions if they enter the
sphere of the radius and centre specified in the main configuration file.
Hope that helps.
Olga
On Thu, Jan 12, 2017 at 3:42 AM, Azadeh Alavizargar <
azadeh.alavizargar_at_gmail.com> wrote:
>
> Hello everyone
>
> I have a protein-membrane system with lots of ions. I want to prevent the
> ions to enter a hypothetical cylinder aligned with the pore of the protein
> using the upperBoundary of the distance xy of the colvar module. Do I have
> to have a colvar for each ion or there is a way to do it for all the ions
> with a single colvar?
> If not, what is the best way to do it?
>
> Thank you very much for your help.
>
> Kind regards,
> Azadeh
>
This archive was generated by hypermail 2.1.6 : Mon Dec 31 2018 - 23:20:01 CST