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.

• constantforce $<$ Apply constant forces? $>$
  Acceptable Values: yes or no
  Default Value: no
  Description: Specifies whether or not constant forces are applied.

• consforcefile $<$ PDB file containing forces to be applied $>$
  Acceptable Values: UNIX filename
  Description: The X, Y, Z and occupancy (O) fields of this file are read to determine the constant force vector of each atom, which is (X,Y,Z)*O, in unit of Kcal/(mol*Å). The occupancy (O) serves as a scaling factor, which could expand the range of the force applied. (One may be unable to record very large or very small numbers in the data fields of a PDB file due to limited space). Zero forces are ignored.

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.

• eFieldOn $<$ apply electric field? $>$
  Acceptable Values: yes or no
  Default Value: no
  Description: Specifies whether or not an electric field is applied.

• eField $<$ electric field vector $>$
  Acceptable Values: vector of decimals (x y z)
  Description: Vector which describes the electric field to be applied. Units are ${\rm kcal} / ({\rm mol \; \AA} \; e)$, which is natural for simulations. This parameter may be changed between run commands, allowing a square wave or other approximate wave form to be applied.

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

$$\vec r(t) \; = \; \vec r_0 \, + \, \vec v t \,.$$ (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.

• movingConstraints $<$ Are moving constraints active $>$
  Acceptable Values: on or off
  Default Value: off
  Description: Should moving restraints be applied to the system. If set to on, then movingConsVel must be defined. May not be used with rotConstraints.

• movingConsVel $<$ Velocity of the reference position movement $>$
  Acceptable Values: vector in Å/timestep
  Description: The velocity of the reference position movement. Gives both absolute value and direction

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.

• rotConstraints $<$ Are rotating constraints active $>$
  Acceptable Values: on or off
  Default Value: off
  Description: Should rotating restraints be applied to the system. If set to on, then rotConsAxis, rotConsPivot and rotConsVel must be defined. May not be used with movingConstraints.

• rotConsAxis $<$ Axis of rotation $>$
  Acceptable Values: vector (may be unnormalized)
  Description: Axis of rotation. Can be any vector. It gets normalized before use. If the vector is 0, no rotation will be performed, but the calculations will still be done.

• rotConsPivot $<$ Pivot point of rotation $>$
  Acceptable Values: position in Å
  Description: Pivot point of rotation. The rotation axis vector only gives the direction of the axis. Pivot point places the axis in space, so that the axis goes through the pivot point.

• rotConsVel $<$ Angular velocity of rotation $>$
  Acceptable Values: rate in degrees per timestep
  Description: Angular velocity of rotation, degrees/timestep.

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:

• In harmonic constraints, each tagged atom is harmonically constrained to a reference point which moves with constant velocity. In SMD, it is the center of mass of the tagged atoms which is constrained to move with constant velocity.

• In harmonic constraints, each tagged atom is constrained in all three spatial dimensions. In SMD, tagged atoms are constrained only along the constraint direction.

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:

$$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.$$ (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.

• SMD $<$ Are SMD features active $>$
  Acceptable Values: on or off
  Default Value: off
  Description: Should SMD harmonic constraint be applied to the system. If set to on, then SMDk, SMDFile, SMDVel, and SMDDir must be defined. Specifying SMDOutputFreq is optional.

• SMDFile $<$ SMD constraint reference position $>$
  Acceptable Values: UNIX filename
  Description: File to use for the initial reference position for the SMD harmonic constraints. All atoms in this PDB file with a nonzero value in the occupancy column will be tagged as SMD atoms. The coordinates of the tagged SMD atoms will be used to calculate the initial center of mass. During the simulation, this center of mass will move with velocity SMDVel in the direction SMDDir.

• SMDk $<$ force constant to use in SMD simulation $>$
  Acceptable Values: positive real
  Description: SMD harmonic constraint force constant. Must be specified in kcal/mol/Å$^2$. The conversion factor is 1 kcal/mol = 69.479 pN Å.

• SMDVel $<$ Velocity of the SMD reference position movement $>$
  Acceptable Values: nonzero real, Å/timestep
  Description: The velocity of the SMD center of mass movement. Gives the absolute value.

• SMDDir $<$ Direction of the SMD center of mass movement $>$
  Acceptable Values: non-zero vector
  Description: The direction of the SMD reference position movement. The vector does not have to be normalized, it is normalized by NAMDbefore being used.

• SMDOutputFreq $<$ frequency of SMD output $>$
  Acceptable Values: positive integer
  Default Value: 1
  Description: The frequency in timesteps with which the current SMD data values are printed out.

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.

• IMDon $<$ is IMD active? $>$
  Acceptable Values: on or off
  Default Value: off
  Description: Specifies whether or not to listen for an IMD connection.

• IMDport $<$ port number to expect a connection on $>$
  Acceptable Values: positive integer
  Description: This is a free port number on the machine that node 0 is running on. This number will have to be entered into VMD.

• IMDfreq $<$ timesteps between sending coordinates $>$
  Acceptable Values: positive integer
  Description: This allows coordinates to be sent less often, which may increase NAMD performance or be necessary due to a slow network.

• IMDwait $<$ wait for an IMD connection? $>$
  Acceptable Values: yes or no
  Default Value: no
  Description: If no, NAMD will proceed with calculations whether a connection is present or not. If yes, NAMD will pause at startup until a connection is made, and pause when the connection is lost.

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:

• tclForces $<$ is Tcl interface 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 tclForcesScript parameters is executed.

• tclForcesScript $<$ input for Tcl interface $>$
  Acceptable Values: file or {script}
  Description: Must contain either the name of a Tcl script file or the script itself between { and } (may include multiple lines). This parameter may occur multiple times and scripts will be executed in order of appearance. The script(s) should perform any required initialization on the Tcl interpreter, including requesting data needed during the first timestep, and define a procedure calcforces { } to be called every timestep.

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

• print <anything>
  This command should be used instead of puts to display output. For example, ``SPMampprint Hello World&''.

• atomid <segname> <resid> <atomname>
  Determines atomid of an atom from its segment, residue, and name. For example, ``atomid br 2 N''.

• addatom <atomid>
  Request coordinates of this atom for next force evaluation. Request remains in effect until clearconfig is called. For example, ``addatom 4'' or ``addatom [atomid br 2 N]''.

• addgroup <atomid list>
  Request center of mass coordinates of this group for next force evaluation. Returns a group ID which is of the form gN where N is a small integer. This group ID may then be used to find coordinates and apply forces just like a regular atom ID. Aggregate forces may then be applied to the group as whole. Request remains in effect until clearconfig is called. For example, ``set groupid [addgroup { 14 10 12 }]''.

• clearconfig
  Clears the current list of requested atoms. After clearconfig, calls to addatom and addgroup can be used to build a new configuration.

• loadcoords <varname>
  Loads requested atom and group coordinates (in Å) into a local array. loadcoords should only be called from within the calcforces procedure. For example, ``loadcoords p'' and ``print p(4)''.

• loadforces <varname>
  Loads the forces applied in the previous timestep (in kcal mol$^{-1}$ Å$^{-1}$) into a local array. loadforces should only be called from within the calcforces procedure. For example, ``loadforces f'' and ``print f(4)''.

• loadmasses <varname>
  Loads requested atom and group masses (in amu) into a local array. loadmasses should only be called from within the calcforces procedure. For example, ``loadcoords m'' and ``print m(4)''.

• addforce <atomid|groupid> <force vector>
  Applies force (in kcal mol$^{-1}$ Å$^{-1}$) to atom or group. addforce should only be called from within the calcforces procedure. For example, ``addforce $groupid { 1. 0. 2. }''.

Several vector routines from the VMD Tcl interface are also defined.