next up previous contents index
Next: Free Energy of Conformational Up: Additional Simulation Parameters Previous: Pressure Control   Contents   Index

Subsections

Applied Forces and Analysis

There are several ways to apply external forces to simulations with NAMD. These are described below.

Constant Forces

NAMD provides the ability to apply constant forces to some atoms. There are two parameters that control this feature.

External Electric Field

NAMD provides the ability to apply a constant electric field to the molecular system being simulated. Energy due to the external field will be reported in the MISC column and may be discontinuous in simulations using periodic boundary conditions if, for example, a charged hydrogen group moves outside of the central cell. There are two parameters that control this feature.

Moving Constraints

Moving constraints feature works in conjunction with the Harmonic Constraints (see an appropriate section of the User's guide). The reference positions of all constraints will move according to

\begin{displaymath} \vec r(t) \; = \; \vec r_0 \, + \, \vec v t \,. \end{displaymath} (1)

A velocity vector $\vec v$ (movingConsVel) needs to be specified.

The way the moving constraints work is that the moving reference position is calculated every integration time step using Eq. 1, where $\vec v$ is in Å/timestep, and $t$ is the current timestep (i.e., firstTimestep plus however many timesteps have passed since the beginning of NAMD run). Therefore, one should be careful when restarting simulations to appropriately update the firstTimestep parameter in the NAMD configuration file or the reference position specified in the reference PDB file.

NOTE: NAMD actually calculates the constraints potential with $U = k (x-x_0)^d$ and the force with $F = d k (x-x_0)$, where $d$ is the exponent consexp. The result is that if one specifies some value for the force constant $k$ in the PDB file, effectively, the force constant is $2 k$ in calculations. This caveat was removed in SMD feature.

The following parameters describe the parameters for the moving harmonic constraint feature of NAMD.

Rotating Constraints

The constraints parameters are specified in the same manner as for usual (static) harmonic constraints. The reference positions of all constrained atoms are then rotated with a given angular velocity about a given axis. If the force constant of the constraints is sufficiently large, the constrained atoms will follow their reference positions.

A rotation matrix $M$ about the axis unit vector $v$ is calculated every timestep for the angle of rotation corresponding to the current timestep. angle = $\Omega t$, where $\Omega$ is the angular velocity of rotation.

From now on, all quantities are 3D vectors, except the matrix $M$ and the force constant $K$.

The current reference position $R$ is calculated from the initial reference position $R_0$ (at $t=0$), $R = M (R_0 - P) + P$, where $P$ is the pivot point.

Coordinates of point N can be found as $N = P + ( (R - P) \cdot v ) v$. Normal from the atom pos to the axis is, similarly, normal $= ( P + ( (X - P) \cdot v ) v ) - X$ The force is, as usual, $F = K (R - X)$; This is the force applied to the atom in NAMD (see below). NAMD does not know anything about the torque applied. However, the torque applied to the atom can be calculated as a vector product torque $= F \times normal$ Finally, the torque applied to the atom with respect to the axis is the projection of the torque on the axis, i.e., $torque_{proj} = torque \cdot v$

If there are atoms that have to be constrained, but not moved, this implementation is not suitable, because it will move all reference positions.

Only one of the moving and rotating constraints can be used at a time.

Using very soft springs for rotating constraints leads to the system lagging behind the reference positions, and then the force is applied along a direction different from the "ideal" direction along the circular path.

Pulling on N atoms at the same time with a spring of stiffness K amounts to pulling on the whole system by a spring of stiffness NK, so the overall behavior of the system is as if you are pulling with a very stiff spring if N is large.

In both moving and rotating constraints the force constant that you specify in the constraints pdb file is multiplied by 2 for the force calculation, i.e., if you specified $K = 0.5 \; {\rm kcal}/{\rm mol}/{\rm\AA}^2$ in the pdb file, the force actually calculated is $F = 2 K (R-X) = 1 \; {\rm kcal}/{\rm mol}/{\rm\AA}^2 \; (R-X)$. SMD feature of namd2 does the calculation without multiplication of the force constant specified in the config file by 2.

Targeted Molecular Dynamics (TMD)

In TMD, subset of atoms in the simulation is guided towards a final 'target' structure by means of steering forces. At each timestep, the RMS distance between the current coordinates and the target structure is computed (after first aligning the target structure to the current coordinates). The force on each atom is given by the gradient of the potential

\begin{displaymath} U_{TMD} = \frac{1}{2} \frac{k}{N} \left[ RMS(t) - RMS^*(t) \right]^2 \end{displaymath} (2)

where $RMS(t)$ is the instantaneous best-fit RMS distance of the current coordinates from the target coordinates, and $RMS^*(t)$ evolves linearly from the initial RMSD at the first TMD step to the final RMSD at the last TMD step. The spring constant $k$ is scaled down by the number $N$ of targeted atoms.

Steered Molecular Dynamics (SMD)

The SMD feature is independent from the harmonic constraints, although it follows the same ideas. In both SMD and harmonic constraints, one specifies a PDB file which indicates which atoms are 'tagged' as constrained. The PDB file also gives initial coordinates for the constraint positions. One also specifies such parameters as the force constant(s) for the constraints, and the velocity with which the constraints move.

There are two major differences between SMD and harmonic constraints:

The center of mass of the SMD atoms will be harmonically constrained with force constant $k$ (SMDk) to move with velocity $v$ (SMDVel) in the direction $\vec n$ (SMDDir). SMD thus results in the following potential being applied to the system:

\begin{displaymath} U(\vec r_1, \vec r_2, ..., t) \; = \; \frac{1}{2} k\left[vt - (\vec R(t) - \vec R_0)\cdot \vec n \right]^2. \end{displaymath} (3)

Here, $t \equiv N_{ts} dt$ where $N_{ts}$ is the number of elapsed timesteps in the simulation and $dt$ is the size of the timestep in femtoseconds. Also, $\vec R(t)$ is the current center of mass of the SMD atoms and $R_0$ is the initial center of mass as defined by the coordinates in SMDFile. Vector $\vec n$ is normalized by NAMD before being used.

Output

NAMD provides output of the current SMD data. The frequency of output is specified by the SMDOutputFreq parameter in the configuration file. Every SMDOutputFreq timesteps NAMD will print the current timestep, current position of the center of mass of the restrained atoms, and the current force applied to the center of mass (in piconewtons, pN). The output line starts with word SMD

Parameters

The following parameters describe the parameters for the SMD feature of NAMD.

Interactive Molecular Dynamics (IMD)

NAMD now works directly with VMD to allow you to view and interactively steer your simulation. With IMD enabled, you can connect to NAMD at any time during the simulation to view the current state of the system or perform interactive steering.

Tcl Forces and Analysis

NAMD provides a limited Tcl scripting interface designed for applying forces and performing on-the-fly analysis. This interface is efficient if only a few coordinates, either of individual atoms or centers of mass of groups of atoms, are needed. In addition, information must be requested one timestep in advance. To apply forces individually to a potentially large number of atoms, use tclBC instead as described in Sec. 6.6.9. The following configuration parameters are used to enable the Tcl interface:

At this point only low-level commands are defined. In the future this list will be expanded. Current commands are:

With the commands above and the functionality of the Tcl language, one should be able to perform any on-the-fly analysis and manipulation. To make it easier to perform certain tasks, some Tcl routines are provided below.

Several vector routines (vecadd, vecsub, vecscale) from the VMD Tcl interface are defined. Please refer to VMD manual for their usage.

The following routines take atom coordinates as input, and return some geometry parameters (bond, angle, dihedral).

The following routines calculate the derivatives (gradients) of some geometry parameters (angle, dihedral).

As an example, here's a script which applies a harmonic constraint (reference position being 0) to a dihedral. Note that the ``addenergy'' line is not really necessary - it simply adds the calculated constraining energy to the MISC column, which is displayed in the energy output. 

tclForcesScript {
  # The IDs of the four atoms defining the dihedral
  set aid1 112
  set aid2 123
  set aid3 117
  set aid4 115
  
  # The "spring constant" for the harmonic constraint
  set k 3.0
  
  addatom $aid1
  addatom $aid2
  addatom $aid3
  addatom $aid4
  
  set PI 3.1416
  
  proc calcforces {} {
  
    global aid1 aid2 aid3 aid4 k PI
    
    loadcoords p
    
    # Calculate the current dihedral
    set phi [getdihedral $p($aid1) $p($aid2) $p($aid3) $p($aid4)]
    # Change to radian
    set phi [expr $phi*$PI/180]
    
    # (optional) Add this constraining energy to "MISC" in the energy output
    addenergy [expr $k*$phi*$phi/2.0]
    
    # Calculate the "force" along the dihedral according to the harmonic constraint
    set force [expr -$k*$phi]
    
    # Calculate the gradients
    foreach {g1 g2 g3 g4} [dihedralgrad $p($aid1) $p($aid2) $p($aid3) $p($aid4)] {}
    
    # The force to be applied on each atom is proportional to its
    # corresponding gradient
    addforce $aid1 [vecscale $g1 $force]
    addforce $aid2 [vecscale $g2 $force]
    addforce $aid3 [vecscale $g3 $force]
    addforce $aid4 [vecscale $g4 $force]
  }
}


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.

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:

The following commands are available from within the calcforces procedure:

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}

next up previous contents index
Next: Free Energy of Conformational Up: Additional Simulation Parameters Previous: Pressure Control   Contents   Index
namd@ks.uiuc.edu