Tcl Boundary Forces

While the tclForces interface described above is very flexible, it is only efficient for applying forces to a small number of pre-selected atoms. Applying forces individually to a potentially large number of atoms, such as applying boundary conditions, is much more efficient with the tclBC facility described below.

• tclBC $<$ are Tcl boundary forces active? $>$
  Acceptable Values: on or off
  Default Value: off
  Description: Specifies whether or not Tcl interface is active. If it is set to off, then no Tcl code is executed. If it is set to on, then Tcl code specified in the tclBCScript parameter is executed.

• tclBCScript $<$ input for Tcl interface $>$
  Acceptable Values: {script}
  Description: Must contain the script itself between { and } (may include multiple lines). This parameter may occur only once. The script(s) should perform any required initialization on the Tcl interpreter and define a procedure calcforces <step> <unique> [args...] to be called every timestep.

• tclBCArgs $<$ extra args for tclBC calcforces command $>$
  Acceptable Values: {args...}
  Description: The string (or Tcl list) provided by this option is appended to the tclBC calcforces command arguments. This parameter may appear multiple times during a run in order to alter the parameters of the boundary potential function.

The script provided in tclBCScript and the calcforces procedure it defines are executed in multiple Tcl interpreters, one for every processor that owns patches. These tclBC interpreters do not share state with the Tcl interpreter used for tclForces or config file parsing. The calcforces procedure is passed as arguments the current timestep, a ``unique'' flag which is non-zero for exactly one Tcl interpreter in the simulation (that on the processor of patch zero), and any arguments provided to the most recent tclBCArgs option. The ``unique'' flag is useful to limit printing of messages, since the command is invoked on multiple processors.

The print, vecadd, vecsub, vecscale, getbond, getangle, getdihedral, anglegrad, and dihedralgrad commands described under tclForces are available at all times.

The wrapmode <mode> command, available in the tclBCScript or the calcforces procedure, determines how coordinates obtained in the calcforces procedure are wrapped around periodic boundaries. The options are:

• patch, (default) the position in NAMD's internal patch data structure, requires no extra calculation and is almost the same as cell

• input, the position corresponding to the input files of the simulation

• cell, the equivalent position in the unit cell centered on the cellOrigin

• nearest, the equivalent position nearest to the cellOrigin

The following commands are available from within the calcforces procedure:

• nextatom
  Sets the internal counter to a new atom and return 1, or return 0 if all atoms have been processed (this may even happen the first call). This should be called as the condition of a while loop, i.e., while {[nextatom]} { ... } to iterate over all atoms. One one atom may be accessed at a time.

• dropatom
  Excludes the current atom from future iterations on this processor until cleardrops is called. Use this to eliminate extra work when an atom will not be needed for future force calculations. If the atom migrates to another processor it may reappear, so this call should be used only as an optimization.

• cleardrops
  All available atoms will be iterated over by nextatom as if dropatom had never been called.

• getcoord
  Returns a list {x y z} of the position of the current atom wrapped in the periodic cell (if there is one) in the current wrapping mode as specified by wrapmode.

• getcell
  Returns a list of 1-4 vectors containing the cell origin (center) and as many basis vectors as exist, i.e., {{ox oy oz} {ax ay az} {bx by bz} {cx cy cz}}. It is more efficient to set the wrapping mode than to do periodic image calculations in Tcl.

• getmass
  Returns the mass of the current atom.

• getcharge
  Returns the charge of the current atom.

• getid
  Returns the 1-based ID of the current atom.

• addforce {<fx> <fy> <fz>}
  Adds the specified force to the current atom for this step.

• addenergy <energy>
  Adds potential energy to the BOUNDARY column of NAMD output.

As an example, these spherical boundary condition forces:

sphericalBC        on
sphericalBCcenter  0.0,0.0,0.0
sphericalBCr1      48
sphericalBCk1      10
sphericalBCexp1    2

Are replicated in the following script:

tclBC on
tclBCScript {
  proc veclen2 {v1} {
    foreach {x1 y1 z1} $v1 { break }
    return [expr $x1*$x1 + $y1*$y1 + $z1*$z1]
  }

  # wrapmode input
  # wrapmode cell
  # wrapmode nearest
  # wrapmode patch ;# the default

  proc calcforces {step unique R K} {
    if { $step % 20 == 0 } {
      cleardrops
      # if $unique { print "clearing dropped atom list at step $step" }
    }
    set R [expr 1.*$R]
    set R2 [expr $R*$R]
    set tol 2.0
    set cut2 [expr ($R-$tol)*($R-$tol)]

    while {[nextatom]} {
      # addenergy 1 ; # monitor how many atoms are checked
      set rvec [getcoord]
      set r2 [veclen2 $rvec]
      if { $r2 < $cut2 } {
        dropatom
        continue
      }
      if { $r2 > $R2 } {
        # addenergy 1 ; # monitor how many atoms are affected
        set r [expr sqrt($r2)]
        addenergy [expr $K*($r - $R)*($r - $R)]
        addforce [vecscale $rvec [expr -2.*$K*($r-$R)/$r]]
      }
    }
  }
}

tclBCArgs {48.0 10.0}