Applied Forces and Analysis

Currently, there are two ways to steer simulations with NAMD. One is called ``moving constraints'' and is based on the harmonic constraints feature. The other is a stand-alone feature called ``SMD''. In both cases the user specifies a reference position $\vec r_0$, force constant $k$, velocity of the reference position movement $\vec v$ and the number of the constrained atom $i$. Then during the simulation a potential $U_{SMD}$ and a corresponding force $\vec f_{SMD}$ are added for the specified atom. The potential is computed according to

$$U_{SMD}(t) \; = \; \frac{k \, [\vec r(t) \, - \, \vec r_i(t)]^2}{2} \, ,$$ (1)

and the force acting on atom $i$ is (by differentiating $U_{SMD}$)

$$\vec f_{SMD}(t) \; = \; \, k \, [\vec r(t) \, - \, \vec r_i(t)] \, ,$$ (2)

where t is time, $\vec r_i$ is the position of the atom $i$ and $\vec r$ is the current reference position (position of the restraint point), given by

$$\vec r(t) \; = \; \vec r_0 \, + \, \vec v t \,.$$ (3)

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 Eq. 3. 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. 3, 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. This feature allows the restraint reference positions for all restrained atoms to be moved. Moving the restraint reference positions for a single atom is provided by the SMD feature.

• 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. One has to specify the force constant $k$ (SMDk), the number (1-based indexing) of the atom (SMDAtom), which is restrained to the moving reference position, the initial reference position $\vec r_0$ (SMDRefPos), the absolute value of the velocity of movement $v$ (SMDVel), and the direction of movement $\vec n$ (SMDDir). The velocity $\vec v$ is then given by $\vec v \, = \, v \vec n$. Vector $\vec n$ is normalized by NAMD before being used. Optionally, the frequency of SMD data output can be specified.

The time used in the reference position calculation starts at time $t_0$ which can be specified through SMDTStamp parameter, defaulting to firstTimestep. The reference position is calculated then by

$$\vec r \; = \; \vec r_0 \, + \, \vec n v \, (t - t_0) \, ,$$ (4)

where $t$ is the current timestep (i.e., firstTimestep plus however many timesteps have passed since the beginning of NAMD run). When restarting the simulation it may be useful to set $t_0$ to whatever time the reference position specified in the configuration file refers to.

Changing direction of reference position movement

One may want to change the direction of the reference position movement, e.g., if the restrained atom is moving too slow in the given direction. To check for such slow movement, the minimum allowed average velocity $v_{min}$ (SMDVmin) and the averaging time $\tau_{min}$ (SMDVminTave) have to be specified. Then every $\tau_{min}$ timesteps NAMD will compute the average velocity in the current direction over the past $\tau_{min}$ timesteps, and compare it to $v_{min}$. The average velocity is computed simply by

$$\langle v\rangle_{min} \; = \; \frac{\vec n \cdot [\vec r_i(t) \, - \, \vec r_i(t - \tau_{min}) ]}{\tau_{min}} \, .$$ (5)

In the case that $\langle v\rangle_{min}$ is less than $v_{min}$, a new direction is chosen, and the simulation proceeds.

Every time the direction is changed, the reference point is also changed to

$$\vec r_{0,new} \; = \; \vec r_i(t) + \vec n_{new} F_{min}/k \, ,$$ (6)

where $t$ is the current timestep, and $F_{min}$ (SMDFmin) is the minimum force to which the force imposed by the restraint is reset. Thus, after the direction is changed to $\vec n_{new}$ the force applied to the atom has an absolute value of $F_{min}$ and the direction of $\vec n_{new}$. The time stamp $t_0$ is set to current time.

Choice of new direction

The choice of the new (random) direction is based on the specified initial direction $\vec n$, the angle $\theta_c$ of the cone which has $\vec n$ as its axis, and the method by which the random numbers are generated. The new direction is guaranteed to lie within the specified cone. The method can be ``uniform'' (default) or ``gaussian''. Uniform method uses uniformly distributed random numbers to find angles $\varphi$ and $\theta$ defining the new direction. Angle $\varphi$ is a random number between 0 and 360 degrees, and $\theta$ is a random number between 0 and $\theta_c/2$. Gaussian method uses uniformly distributed random numbers to find $\varphi$ and normally distributed random numbers with variance $\sigma$ (SMDGaussW) to find $\theta$. All $\theta$ values that exceed $\theta_c/2$ are ignored. The distribution of $\theta$ values is given by

$$p(\theta) \, = \, \frac{1}{\sigma \sqrt{2\pi}} \exp(- \frac{\theta^2}{\sigma^2}) \, .$$ (7)

After the angles $\varphi$ and $\theta$ are calculated, the new direction is found by constructing a right orthonormal system of coordinates with $\vec n_0$, the initial direction specified by SMDDir in the configuration file, being the $z$ direction. The other two vectors $\vec a$ and $\vec b$, orthogonal to $\vec n_0$ and to each other are found by first producing $\vec a \, = \, \vec n_0 \times \vec{(0, 1, 1)}$, or $\vec a \, = \, \vec n_0 \times \vec{(1, 1, 1)}$ in the case when $\vec n_0$ is collinear to $\vec{(0, 1, 1)}$, normalizing $\vec a$, and, finally, producing $\vec b \, = \, \vec n_0 \times \vec a$. Then the new direction is given by

$$\vec n_{new} \; = \; \vec a \cos\varphi \sin\theta \, + \, \vec b \sin\varphi \sin\theta \, + \, \vec n_0 \cos\theta \, .$$ (8)

Resetting the applied force

Sometimes, especially when a small enough $k$ is employed, the movement of the restrained atom proceeds in steps, resulting in a fast movement of the atom over a short period of time. When this happens, it may be desirable to reset (reduce) the applied force. Resetting of the force to a specified value $F_{min}$ (SMDFmin) should happen when the average velocity of atom movement $\langle v\rangle_{max}$ over a given period of time $\tau_{max}$ (SMDVmaxTave) exceeds the specified maximum allowed velocity $v_{max}$ (SMDVmax) in the direction of reference position movement $\vec n$. The average velocity is calculated similarly to Eq. 5, and in the case when $\langle v\rangle_{max}$ is larger than $v_{max}$, the reference position is reset according to Eq. 6. The time stamp $t_0$ is set to current time.

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 restrained atom, the current force applied to the restrained atom (in piconewtons, pN), the reference position at the time of the last change of direction or force resetting (initially, the reference position specified in the configuration file), current direction of the reference position movement, and the time stamp, indicating the timestep when the direction was changed or force reset (initially, SMDTStep specified in the configuration file, or its default value). The output line starts with word SMD

The energy of the SMD restraint is printed out along with other energies in the column titled SMD. The frequency of this output is defined by outputEnergies parameter in the configuration file.

When the direction of the reference position movement changes or the force is reset, this change is indicated by an output line starting with SMDChange and showing the current timestep, the current atom position, the new reference position after the change, and the new direction of the reference position movement. The current timestep becomes the time stamp for the following simulation until the next change occurs.

Parameters

The following parameters describe the parameters for the SMD feature of NAMD. This feature allows a harmonic restraining force to be applied to any ONE atom in the simulation and the restraining reference position to be moved. It also allows to randomly change the direction of the reference position movement if the restrained atom progress is too slow, and to reset the applied force to a given value if the progress of the restrained atom is too fast.

• 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, SMDRefPos, SMDVel, SMDDir and SMDAtom must be defined. Also, SMDOutputFreq, SMDChDir and associated parameters, SMDChForce and associated parameters, and SMDTStamp can be optionally defined.

• SMDexp $<$ exponent to use in SMD harmonic constraint energy function $>$
  Acceptable Values: positive, even integer
  Default Value: 2
  Description: Exponent to be used in the SMD harmonic constraint energy function. This value must be a positive integer, and only even values really make sense. This parameter is only used if SMD is set to on.

• SMDRefPos $<$ SMD constraint reference position $>$
  Acceptable Values: vector
  Description: Vector to use for the initial reference position for the SMD harmonic constraints. The atom that is specified by SMDAtom will be initially constrained to this reference position. During the simulation, this reference position 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: positive real, Å/timestep
  Description: The velocity of the SMD reference position movement. Gives the absolute value.

• SMDDir $<$ Direction of the SMD reference position 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.

• SMDAtom $<$ Index (1-based) of the atom constrained to the moving reference position $>$
  Acceptable Values: positive integer
  Description: The index of the atom constrained to the moving reference position. If the atom is numbered $N$ in the PDB file, $N$ must be specified.

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

• SMDTStamp $<$ time of the last change in reference position $>$
  Acceptable Values: non-negative integer
  Default Value: firstTimestep
  Description: The timestep at which the reference position $\vec r_0$ was last reset. Normally, it would be the firstTimestep, but if the force or direction was reset during the simulation, this parameter may be useful for restarts.

• SMDFmin $<$ value of the force to reset to $>$
  Acceptable Values: non-negative real
  Default Value: 0
  Description: The value in pN to which the force gets reset when the direction changes or the average velocity exceeds the given SMDVmax value. This parameter is only used if SMDChDir or SMDChForce are on.

• SMDChDir $<$ Is changing direction option of SMD active $>$
  Acceptable Values: on, off
  Default Value: off
  Description: If turned on, this option allows to change the direction of the reference position movement when the average velocity of the restrained atom in this direction is smaller than the given SMDVmin value. Uses settings of SMDVmin, SMDVminTave, SMDConeAngle, SMDChDirMethod.

• SMDVmin $<$ minimum allowed average velocity of the restrained atom $>$
  Acceptable Values: non-negative real
  Description: The minimum allowed average velocity in Å/timestep of the restrained atom in the direction of the reference position movement. The averaging time is given by SMDVminTave.

• SMDVminTave $<$ averaging time for velocity to compare to Vmin $>$
  Acceptable Values: positive integer
  Description: The averaging time in timesteps for calculation of the average velocity of the restrained atom and for comparison of this average velocity to SMDVmin.

• SMDConeAngle $<$ angle of the cone to choose the new direction from $>$
  Acceptable Values: non-negative
  Default Value: 360
  Description: The angle (in degrees) of the cone with the axis along the initial direction of the reference position movement SMDDir. The new direction will come from the vertex of the cone and will lie within the cone. This angle is twice the angle between the cone's axis and any line going through the cone's vertex on the cone's surface. Omitting this parameter will lead to choosing a direction without any restrictions (all directions are possible).

• SMDChDirMethod $<$ method to use when choosing a new direction $>$
  Acceptable Values: gaussian, uniform
  Default Value: uniform
  Description: The method to choose the angle between the new direction and the initial direction of the reference position movement. When gaussian is specified, it is recommended to set the variance of the distribution with SMDGaussW.

• SMDGaussW $<$ variance of the gaussian distribution used to choose new directions $>$
  Acceptable Values: positive real
  Default Value: 360
  Description: Variance of the gaussian distribution (in degrees) used for choosing the angle between the new direction and the initial direction of the reference position movement.

• SMDChForce $<$ is resetting force SMD option active $>$
  Acceptable Values: on, off
  Default Value: off
  Description: If turned on, this option allows to reset the reference position (so that the applied force becomes SMDFmin) when the average velocity of the restrained atom in the direction of the reference position movement is larger than a given SMDVmax value. Uses settings of SMDVmax and SMDVmaxTave.

• SMDVmax $<$ maximum allowed average velocity of the restrained atom $>$
  Acceptable Values: non-negative real
  Description: The maximum allowed average velocity in Å/timestep of the restrained atom in the direction of the reference position movement. The averaging time is given by SMDVmaxTave.

• SMDVmaxTave $<$ averaging time for velocity to compare to Vmax $>$
  Acceptable Values: positive integer
  Description: The averaging time in timesteps for calculation of the average velocity of the restrained atom and for comparison of this average velocity to SMDVmax.

Interactive Molecular Dynamics (IMD)

NAMD now works directly with VMD to allow you to view and interactively steer your simulation. NAMD will wait for VMD to connect on startup.

• IMDon $<$ is IMD active? $>$
  Acceptable Values: on or off
  Default Value: off
  Description: Specifies whether or not to expect and wait for a 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. NAMD will wait for a connection on this port. 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.

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 <molname> <resid> <atomname>
  Determines atomid of an atom from its molecule, 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 reconfig 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 reconfig is called. For example, ``set groupid [addgroup { 14 10 12 }]''.

• reconfig
  Signals that new atoms are being requested. addatom and addgroup calls during calcforces will be ignored unless reconfig is called. Old configuration is replaced by new configuration. reconfig should only be called from within the calcforces procedure.

• 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)''.

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