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.

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} (2)

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 interface

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

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]
    
    # 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]
  }
}

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