tclBC during SMD and ABF simulations

From: Olya Kravchenko (
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

# 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) + \

      if { $atomid == 344742 || [getmass] < 22 || [getmass] > 36 } {

      } 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


This archive was generated by hypermail 2.1.6 : Tue Dec 27 2016 - 23:22:38 CST