apply force along an internally defined direction

From: Vlad Cojocaru (
Date: Mon Oct 29 2007 - 07:24:02 CDT

Dear NAMD users,

I would like to use Tcl forces to move an atom with constant velocity
along a directon defined by two other atoms in the protein (an internal
direction which would change between time steps in function of how the
protein translates and/or rotates). Below I attach an attempt of script
to do that. However, I am not sure whether this is a correct script and
NAMD complains about the "vecnorm" command.

First I would like to ask whether somebody has performed such a task in
NAMD and if the script below seems reasonablefor this? If yes, could
somebody indicate me how to calculate the direction in each time step
without using "vecnorm" (using only routines implemented in namd)?

Thanks in advance
Best wishes

----tcl script ---------

set id1 7614
set id2 1193
set id3 3101

set atoms {}
lappend atoms $id1 $id2 $id3
foreach atom $atoms {
 addatom $atom
# set the SMD output frequency, initialize the time counter,
set Tclfreq 50
set t 0
#set force constant in kcal/mol AA^2 and pulling vel. in AA/ts;
set k 7.2
#set velocity to be 35 A/100 ps
set v 0.00035

# initial coordinates of pulled atom (t = 0)
set x0 52.842
set y0 50.161
set z0 51.822

set outfilename smd_tcl.out
open $outfilename w

proc calcforces {} {

 global atoms id1 id2 id3 Tclfreq t k v x0 y0 z0 outfilename
 # form local array coords

 loadcoords coords

# calculate direction of pulling
 set n [vecnorm [expr 0.5*[vecsub $coords($id2) $coords($id3)]]]

# get coords of pulled atom
 set r $coords($id1)

 # get the new x,y,z components and calculate the force

 set x1 [lindex $r 0]
 set y1 [lindex $r 1]
 set z1 [lindex $r 2]
 set nx [lindex $n 0]
 set ny [lindex $n 1]
 set nz [lindex $n 2]
 set vx [expr $v*$nx]
 set vy [expr $v*$ny]
 set vz [expr $v*$nz]
 # calculate forces
 set fx [expr $k*($x0+$vx*$t-$x1)]
 set fy [expr $k*($y0+$vy*$t-$y1)]
 set fz [expr $k*($z0+$vz*$t-$z1)]
 lappend f $fx $fy $fz
 # add forces

 addforce $id1 $f
 # set the output

 set foo [expr $t % $Tclfreq]
 if { $foo == 0 } {
  set outfile [open $outfilename a]
  set time [expr $t/1000.0]
  puts $outfile "$time $r $f"
  close $outfile
 incr t

Dr. Vlad Cojocaru
EML Research gGmbH
Schloss-Wolfsbrunnenweg 33
69118 Heidelberg
Tel: ++49-6221-533266
Fax: ++49-6221-533298
EML Research gGmbH
Amtgericht Mannheim / HRB 337446
Managing Partner: Dr. h.c. Klaus Tschira
Scientific and Managing Director: Prof. Dr.-Ing. Andreas Reuter

This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:45:26 CST