From: Olya Kravchenko (ovkrav_at_gmail.com)
Date: Mon Nov 28 2016 - 12:36:47 CST
Hi all,
I would appreciate if someone could comment on the script I made to push
ions away from the tagged ion during abf simulations inside the ion
channel. It seems to work the way it is now for simple equilibration, but I
am not sure whether this is correct, especially when combined with SMD
and/or ABF. If I apply force to the tagged ion, can it interfere somehow
with the tclBC that I apply to surrounding ions?
The simulation system is a protein channel solvated in NaCl and the
permeating (tagged) ion is sodium with a known atom id number (344742,
obtained by get serial command in VMD). I want to push all ions away from
the tagged ion and leave the protein unaffected. For that I distinguish
ions from the protein's atoms by mass. The sphere radius and the force
strength are set in the main configuration file. Here is a tclBC scrpt that
is called from within the main configuration file:
#Set a list of atom IDs of all the ions
set ionList {}
set atomID 1
set inStream [open $pdbSource r]
set sphereCenter {22.5 89.6 93.2} ; #this is the position of the tagged
ion
# go through pdb file and find entries for ions and save them to 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
# find the components of the sphere's center
foreach { x0 y0 z0 } $sphereCenter { break }
while {[nextatom]} {
set atomid [getid] ;# get the ID of the current atom
# check if this ID is listed in ionList
if { [lsearch $ionList $atomid] >= 0 } {
set rvec [getcoord] ;# get ion's coordinates
foreach { x y z } $rvec { break }
# get 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 { $atomid == 344742 || [getmass] < 22 || [getmass] > 36 } {
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"
}
}
}
Does it look correct? I would very much appreciate your comments! Thanks in
advance,
Olga
This archive was generated by hypermail 2.1.6 : Tue Dec 27 2016 - 23:22:38 CST