In a NAMD configuration file, one sets up the simulation parameters, such as the temperature, pressure, or a uniform electric field, for the molecular system as a whole. In many cases, though, you will need to impose boundary conditions, apply a non-uniform electric field, or selectively apply forces to a group of atoms of the given type or located in a particular area. For this purpose, you can use a tclBC (which stands for Tcl Boundary Conditions) script. A tclBC script is written in the Tcl programming language, a rudimentary knowledges of which would be enough for most tasks. The script will be interpreted line-by-line, therefore slowing down the simulation, sometimes considerably. For this reason you will want to keep your scripts simple and efficient.
![\framebox[\textwidth]{
\begin{minipage}{.2\textwidth}
\includegraphics[width=2...
...ation about the whole system, you will need {\tt tclForces}.
}
\end{minipage} }](img21.gif)
The way the coordinates of atoms are wrapped around periodic boundaries in a tclBC script is controlled using the command wrapmode, which that may have one of the following arguments:
For most purposes, mode cell is a good choice.
Imagine that you have a water box, and you want to create a spherical bubble of vacuum. You can do that by applying forces to every atom found inside this sphere to push it out. To avoid instability, you would want to start with a small bubble and then increase it until the desired size is reached. The rate of increasing the size of the bubble is up to you. If this rate is too large, some atoms will be pushed too hard and therefore move too fast, causing the simulation to crush. Fast-moving atoms may also break the structure of a double-stranded DNA or another molecular complex held together by relatively weak non-bonded interactions (van der Waals and electrostatic forces). Even if the simulation remains stable, you should watch the changes in temperature and pressure to make sure that the energy influx, which is equal to the work done by the force you are applying, has enough time to dissipate.
![\framebox[\textwidth]{
\begin{minipage}{.2\textwidth}
\includegraphics[width=2...
...ion or increasing the temperature or pressure beyond reason.
}
\end{minipage} }](img22.gif)
You should find the following code:
tclBC on
tclBCScript {
set bubbleCenter "0.0 0.0 0.0"
set tclBCScript bubble.tcl
source $tclBCScript
}
tclBCArgs {0. 15. 0.01 5.}
Since a tclBC script is called from NAMD, it is referenced NAMD configuration file, which is also a good place to set the script parameters.
The NAMD command tclBC on turns the TclBC interface on.
In our example, tclBCScript {...} contains the initialization of a key variable and a reference to the file that contains the script itself: source $tclBCScript.
If the script is short, it can be placed entirely within the body of the command tclBCScript in the NAMD configuration file.
Finally, the command tclBCArgs is used to pass a list of variables to the main TclBC procedure calcforces (we will talk about it very soon).
In this case, the four arguments found in curly brackets have the meaning: ``make a bubble starting from radius 0 Å and increase it to 15 Å at a rate of 0.01 Å per simulation step, by applying forces of 5 kcal/mol
Å.''
bubble.tcl contains the following code:
# Two first agruments of calcforces are automatically forwarded
# to it by NAMD. The other 4 arguments match the list of 4 values
# from command tclBCArgs.
proc calcforces {step unique Rstart Rtarget Rrate K} {
global bubbleCenter ;# defined in tclBCScript{ ... }
# increase R, starting from $Rstart, by $Rrate at each step,
# until it reaches $Rtarget; then keep it constant
set R [expr $Rstart + $Rrate * $step]
if { $R > $Rtarget } { set R $Rtarget }
# let only the main processor print the output
if { $unique } {
print "step $step, bubble radius = $R"
}
# get the components of the bubble center vector
foreach { x0 y0 z0 } $bubbleCenter { break }
# pick atoms of the given patch one by one
while {[nextatom]} {
set rvec [getcoord] ;# get the atom's coordinates
# get the components of the vector
foreach { x y z } $rvec { break }
# find the distance between the atom and the bubble center
# (long lines can be broken by a backslash and continued
# on the next line)
set rho [expr sqrt(($x-$x0)*($x-$x0) + ($y-$y0)*($y-$y0) + \
($z-$z0)*($z-$z0))]
# if the atom is inside the sphere, push it away radially
if { $rho < $R } {
set forceX [expr $K * ($x-$x0) / $rho]
set forceY [expr $K * ($y-$y0) / $rho]
set forceZ [expr $K * ($z-$z0) / $rho]
addforce "$forceX $forceY $forceZ"
}
}
}
![\framebox[\textwidth]{
\begin{minipage}{.2\textwidth}
\includegraphics[width=2...
...rgs} line in the NAMD configuration file specifies the rest.
}
\end{minipage} }](img23.gif)
In the code bubble.tcl, we use the unique flag to print something from within calcforces once at each step. Without if { $unique } clause, the print command would be executed by every processor.
The list of variables Rstart Rtarget Rrate K forwarded to calcforces matches the list of values tclBCArgs {0.0 10.0 0.01 5.0} of the NAMD configuration file. You can use more or fewer variables, or none at all, as long as the lengths of the two lists are the same. The automatic variable step and unique don't count, so calcforces will always have two more arguments than to tclBCArgs.
The key command in calcforces is nextatom, which tells NAMD to select the next atom of the given patch. nextatom returns a value which is nonzero if the next atom has been successfully selected, or zero if no more atoms are left to process. Thus, while {[nextatom]} {...} is a loop that will be repeated until all atoms of the current patch are processed.
Once an atom is selected by nextatom, a number of commands can be used to get its properties:
set atomid [getid] ;# 1-based atom ID number set coords [getcoord] ;# 3-component coordinate vector (A) set mass [getmass] ;# mass (atomic units) set charge [getcharge] ;# charge (atomic units)
Using these parameters, you can decide what to do to that atom.
In the script bubble.tcl, we find the coordinates of the atom using the command foreach { x y z } $rvec { break }; this is a standard trick to get the components of a vector or a list in Tcl. Using the same method, we find the coordinates of the bubble center. Then we check whether the atom is inside the bubble of radius $R and center $bubbleCenter. If so ($rho < $R), we calculate the components of the force vector such that it is directed radially from the bubble's center. In this simple example, the absolute value of the force is always $K whenever the force is applied. Finally, addforce tells NAMD to actually push the given atom by adding the force vector ``$forceX $forceY $forceZ'' to the force that would act on that atom otherwise.
namd2 eq04.bubble.namd > eq04.bubble.log
(for standard installations of NAMD.)
Sample snapshots of the simulation trajectory are shown in Fig. 5. The directory tclBCmovies contains MPEG movies bubble*.mpg that show the bubble formation trajectory animation at 0.1 ps per frame.
![\begin{figure}\begin{center}
\par
\par
\latex{
\includegraphics[scale=0.7]{pictures/bubble1all}
}
\end{center}
\end{figure}](img24.gif) |
The above script bubble.tcl is far from being perfect. The same job can be done much more efficiently by decreasing the number of atoms processed at each step, as well as the number of arithmetic operations. Note that Tcl arithmetic is much slower than hard-coded C++ arithmetic of the NAMD. A more efficient script is given in file BUBBLE/bubbleFast.tcl.
![\framebox[\textwidth]{
\begin{minipage}{.2\textwidth}
\includegraphics[width=2...
...ardrops} occasionally and then drop unnecessary atoms again.
}
\end{minipage} }](img25.gif)
wrapmode cell
# this line is moved above calcforces to avoid doing this
# transformation more than once
foreach { x0 y0 z0 } $bubbleCenter { break }
set TOL 3. ;# distance tolerance parameter
proc calcforces {step unique Rstart Rtarget Rrate K} {
global x0 y0 z0 TOL
set R [expr {$Rstart + $Rrate * $step}]
if { $R > $Rtarget } { set R $Rtarget }
if { $unique } {
print "step $step, bubble radius = $R"
}
# Atoms found at a distance larger than $TOL from the surface
# of the bubble, or $RTOL from the center of the bubble, will
# be ignored (dropped) for the rest of a 100-step cycle.
set RTOL [expr {$R + $TOL}]
# restore the complete list of atoms to consider
if { $step % 100 == 0 } { cleardrops }
while {[nextatom]} {
set rvec [getcoord]
foreach { x y z } $rvec { break }
set rho [expr {sqrt(($x-$x0)*($x-$x0) + ($y-$y0)*($y-$y0) + \
($z-$z0)*($z-$z0))}]
# Atoms at distances 0 to $R from the bubble center are pushed,
# atoms father than $RTOL are dropped. Atoms between $R to $RTOL,
# that is with a layer of thickness $TOL, are neither pushed nor
# dropped at this step, so that they will be considered again at
# the next step(s). They may come closer to the bubble and then
# they will have to be pushed.
if { $rho < $R } {
set forceX [expr {$K * ($x-$x0) / $rho}]
set forceY [expr {$K * ($y-$y0) / $rho}]
set forceZ [expr {$K * ($z-$z0) / $rho}]
addforce "$forceX $forceY $forceZ"
} elseif { $rho > $RTOL } {
dropatom ;# no longer consider this atom until "cleardrop"
}
}
}
Since most water molecules are far from the bubble, we don't need to check their coordinates every step. It's enough to do that once in a while (every 100 steps in this example). At all other times, only a small number of atoms that are still inside the bubble or within distance $TOL from it, and only these are selected by nextatom each step.
namd2 eq04.bubbleFast.namd > eq04.bubbleFast.log.
![\framebox[\textwidth]{
\begin{minipage}{.2\textwidth}
\includegraphics[width=2...
...and, is essential to keeping up high simulation performance.
}
\end{minipage} }](img26.gif)
A more complex initialization of the TclBC script is necessary when it is supposed to process a selected group of atoms.
Imagine that, for whatever reason, you want to gather all K
and Cl
ions within a sphere with the same center and radius as in the above example.
For this purpose, we will push all ions toward the center of the sphere.
The force applied will be constant for all ions beyond the sphere, and proportional to the distance between the ion and the sphere center for atoms inside the sphere.
As before, the tclBC script is set up in the NAMD configuration file:
tclBC on
tclBCScript {
set sphereCenter "0.0 0.0 0.0"
set sphereRadius 10.0
set maxForce 5.0
set pdbSource waterbox40-0.2M.pdb
set tclBCScript concentrateIons.tcl
source $tclBCScript
}
tclBCArgs { }
In this case, we don't explicitly pass any extra arguments to calcforces (the variables step and unique will still be passed implicitly). Instead, we set up global variables sphereCenter, sphereRadius, maxForce, and pdbSource. They will be made visible in calcforces using command global.
To identify the ions, we will read the system's PDB file, locate the ions' entries, and make a list of their IDs.
The ion entries in the PDB file for a large system may look like
ATOM ***** CLA CLA 170 -36.100 -48.734 7.045 1.00 0.00 ION CL ... ATOM ***** POT POT 1 -0.574 -49.643 34.378 1.00 0.00 POT K
In the above example, atom IDs for the ions are above 99999, and therefore not shown in the PDB file by a decimal. Thus, to find out the atom IDs of the ions, we will need to read the PDB file line by line while counting the atom entries. An entry for each atom must begin with four letters ATOM or HETA. In the PDB file above, each ion is identified by the residue name CLA or POT, which occupies positions 17-19 of the line. Note that a 4-letter residue name would occupy positions 17-20 of the line. To make the script more general, we will read symbols 17-20 and then trim them (remove spaces).
The file concentrateIons.tcl begins with the following code:
set ionList {} ;# a list of atom IDs of the ions
set atomID 1 ;# in NAMD, atom IDs are 1-based numbers
set inStream [open $pdbSource r] ;# opening the PDB file
# reading the PDB file line after line:
foreach line [split [read $inStream] \n] {
# symbols 0-3 should be ATOM or HETA
set string1 [string range $line 0 3]
# symbols 17-20 contain the residue name
set string2 [string range $line 17 20]
set string2 [string trim $string2] ;# trimming the residue name
if { [string equal $string1 {ATOM}] || \
[string equal $string1 {HETA}] } {
# so it's a valid atom entry; let's see if this is an ion
if { [string equal $string2 {CLA}] || \
[string equal $string2 {POT}] } {
# yes it's an ion; append its index to the list
lappend ionList $atomID
}
# increase the atomID counter by 1 (default increment)
incr atomID
}
}
close $inStream
This part of the code, preceding proc calcforces, will be executed once by each processor. Identical instances of the list ionList of atom IDs for the ions will also be kept in memory at each processor.
Now we are ready to write the main part of the code:
wrapmode cell
# in this case the command tclBCArgs {} located in the NAMD config
# file doesn't pass any arguments to calcforces, therefore only
# step and unique will be listed below:
proc calcforces {step unique} {
# list all global variables set outside this procedure
# (except those we won't need, like pdbSource)
global 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 it's found,
# lsearch will return the position of the search pattern
# in the list (a number >= 0), otherwise -1 is returned
if { [lsearch $ionList $atomid] >= 0 } {
set rvec [getcoord] ;# get the ion's coordinates
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))]
# Apply same force $maxForce to each ion if it's outside the
# sphere. The components of the force vector are chosen so
# that the vector is directed toward the sphere's center,
# and the vector's norm is $maxForce.
if { $rho > $sphereRadius } {
set forceX [expr -$maxForce * ($x-$x0) / $rho]
set forceY [expr -$maxForce * ($y-$y0) / $rho]
set forceZ [expr -$maxForce * ($z-$z0) / $rho]
} else {
# If the ion in already inside the sphere, scale the force
# by a factor of $rho/$sphereRadius, so that the force
# decreases from $maxForce to 0 as the ion approaches the
# sphere's center:
set forceX [expr -$maxForce * ($x-$x0) / $sphereRadius]
set forceY [expr -$maxForce * ($y-$y0) / $sphereRadius]
set forceZ [expr -$maxForce * ($z-$z0) / $sphereRadius]
}
# Finally, applying the force calculated above
# to the current atom
addforce "$forceX $forceY $forceZ"
}
}
}
and 12 (lime) for K ions.
Sample snapshots of the ion concentration trajectory are shown in Fig. 6. The directory tclBCmovies contains MPEG movie ionConcentration.mpg that shows the process at 0.1 ps per frame.
![\begin{figure}\begin{center}
\par
\par
\latex{
\includegraphics[scale=0.7]{pictures/ionsConc1all}
}
\end{center}
\end{figure}](img29.gif) |
The above code can be made both more computationally efficient, as well as more concise and elegant. In our case, we are pushing only the ions, therefore all other atoms can be safely dropped using the statement
if { [lsearch $ionList $atomid] == -1 } {
dropatom
continue
}
This way all non-ion atoms will be dropped at the first step and never selected again by nextatom until cleardrops is called, thus making the code more efficient. Some extra time could be saved if we didn't need to continue pushing the atoms that are already inside the sphere:
if { $rho < $sphereRadius } {
dropatom ;# don't process this atom in the future
continue ;# don't process this atom now, either
}
# time to time (every 100 steps in this example),
# let's restore the complete list of atoms
# to catch ions escaping from the sphere
if { $step % 100 == 0 } { cleardrops }
Using the vector routines, instead of robust but lengthy (and therefore error-prone) code
foreach { x0 y0 z0 } $sphereCenter { break }
foreach { x y z } $rvec { break }
set rho [expr sqrt(($x-$x0)*($x-$x0) + ($y-$y0)*($y-$y0) + \
($z-$z0)*($z-$z0))] ;# find the distance between two points
# calculate the force components
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" ;# apply the force
we could write briefly and clearly:
set relativePosition [vecsub $rvec $sphereCenter]
set rho [veclength $relativePosition]
set forceVector [vecscale $relativePosition \
[expr -$maxForce / $rho]]
addforce $forceVector
or even more briefly (though maybe less clearly):
set relativePosition [vecsub $rvec $sphereCenter]
addforce [vecscale [vecscale $relativePosition \
[expr -$maxForce / [veclength $relativePosition]]]]
Another trick can simplify a tclBC script if you need to process all atoms of a certain type, as K and Cl in the above example.
Since such atoms have unique masses, you could avoid reading the PDB file, and make the list of atoms by writing, e.g.:
while { [nextatom] } {
if { [getmass] > 35 && [getmass] < 40 } { ... }
}
Likewise, you could use if { [getcharge] } as a selection criterion. Just make sure that no other atoms could match the conditions of the selection.
# we will use the following variable to calculate and print the
# average number of ions found outside the sphere at each step
set avgNumIons 0
set forceCoef [expr -$maxForce/$sphereRadius]
wrapmode cell
##################################################
proc calcforces {step unique} {
global sphereCenter sphereRadius maxForce avgNumIons
if { $step > 0 && $step % 100 == 0 } {
# Calculate and print the average number ions found outside
# the sphere
set avgNumIons [expr $avgNumIons / 100.]
print "Step $step, average number of ions outside \
the sphere: $avgNumIons"
set avgNumIons 0
cleardrops
}
while {[nextatom]} {
if { [getmass] < 35 || [getmass] > 40 } {
dropatom ;# not a K+ or Cl-, forget about it
continue
}
# vector between the ion and the sphere's center
set relativePosition [vecsub [getcoord] $sphereCenter]
set rho [veclength $relativePosition]
if { $rho > $sphereRadius } {
addforce [vecscale $relativePosition [expr -$maxForce/$rho]]
incr avgNumIons
} else {
dropatom ;# this ion is already inside the sphere
}
}
}
As previously mentioned, the main feature of TclBC scripts is their ability to apply geometry-based forces.
Using a TclBC script, it is possible to impose arbitrary fields of forces, thus imitating complex geometries and interactions while keeping the molecular system relatively small and simple.
For example, a solution containing DNA in a microfluidic channel undergoes a shear flow, so that parallel layers of the liquid move with different velocities.
We can simulate the shear flow by applying force to oxygen atoms of water molecules found within thin layers in XY planes at and around 
= zLo and = zHi, which variables will be set in the NAMD configuration file.
The following initialization is made in the NAMD configuration file:
tclBC on
tclBCScript {
set zLo -15. ;# lower plane where forces are applied
set zHi 15. ;# top plane where forces are applied
set dz 3. ;# half-width of the layer
set TOL 6. ;# drop atoms that far from either layer
set force 5.
set pdbSource dsDNA6_solv.pdb
set tclBCScript shear.tcl
source $tclBCScript
}
tclBCArgs { }
set atomList {} ;# a list of atom IDs of water oxygens
set atomID 1 ;# in NAMD, atom IDs are 1-based numbers
set inStream [open $pdbSource r] ;# opening the PDB file
# reading the PDB file line after line:
foreach line [split [read $inStream] \n] {
# symbols 0-3 should be ATOM or HETA
set string1 [string range $line 0 3]
# symbols 13-16 contain the atom's name
set string2 [string range $line 13 16]
# trimming the atom's name
set string2 [string trim $string2]
if { [string equal $string1 {ATOM}] || \
[string equal $string1 {HETA}] } {
# so it's a valid atom entry; let's see if
# this is a water's oxygen
if { [string equal $string2 {OH2}] } {
# yes it is; append its index to the list
lappend atomList $atomID
}
incr atomID
}
}
close $inStream
puts "[llength $atomList] water oxygens atoms found"
set totalTorque 0.
wrapmode cell
####################################################
proc calcforces {step unique} {
global atomList shearStress zLo zHi dz TOL force
# the value of a global variable will be stored
# between calls to this procedure, while a local
# variable would be lost
global totalTorque
if { $step % 100 == 0 } {
if { $totalTorque != 0. } {
print "Step $step, total torque applied: $totalTorque"
set totalTorque 0.
}
cleardrops
}
while {[nextatom]} {
# check if this atom should be considered at all
if { [lsearch $atomList [getid]] == -1 } {
dropatom
continue
}
# now check if it's within bounds zLo to zHi
foreach { x y z } [getcoord] { break } ;# get coordinates
if { $z >= [expr {$zLo-$dz}] && $z <= [expr {$zLo+$dz}] } {
addforce "[expr {-$force}] 0.0 0.0"
set totalTorque [expr {$totalTorque - $force * $z}]
} elseif { $z >= [expr {$zHi-$dz}] && \
$z <= [expr {$zHi+$dz}] } {
addforce "$force 0.0 0.0"
set totalTorque [expr {$totalTorque + $force * $z}]
} elseif { ( $z >= [expr {$zLo-$TOL}] && \
$z <= [expr {$zLo+$TOL}] ) || \
( $z >= [expr {$zHi-$TOL}] && \
$z <= [expr {$zHi+$TOL}] ) } {
continue ;# keep an eye on it, it may come closer next time
} else {
dropatom ;# this atom is too far, forget about it for now
}
}
}
Sample snapshots of the simulation trajectory are shown in Fig. 7.
MPEG movies (found in directory tclBCmovies) shearSlow.small.mpg and
shearSlow.big.mpg show the trajectory animation at 0.1 ps per frame.
At this slow framerate, the movement of water layers is clearly visible.
In movies shearFast.small.mpg and shearFast.big.mpg, the frame rate is 1 ps per frame, which allows one to see the slow rotation and dissociation of DNA in the shear flow.
![\begin{figure}\begin{center}
\par
\par
\latex{
\includegraphics[scale=0.85]{pictures/shear1all}
}
\end{center}
\end{figure}](img30.gif) |