Re: Colvar distance xy to restrain ions

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

particularly this section:
and I based my script on sample scripts that come along with it:

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

      if { [getid] == 344742 || [getid] < 344526 } {

      } 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.


On Thu, Jan 12, 2017 at 3:42 AM, Azadeh Alavizargar <> 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 : Sun Dec 31 2017 - 23:20:59 CST