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

Subsections

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 USMD and a corresponding force $\vec f_{SMD}$ are added for the specified atom. The potential is computed according to

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

and the force acting on atom i is (by differentiating USMD)
\begin{displaymath} \vec f_{SMD}(t) \; = \; \, k \, [\vec r(t) \, - \, \vec r_i(t)] \, , \end{displaymath} (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
\begin{displaymath} \vec r(t) \; = \; \vec r_0 \, + \, \vec v t \,. \end{displaymath} (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-x0)d and the force with F = d k (x-x0), 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.

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 R0 (at t=0), R = M (R0 - 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 x 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. 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 t0 which can be specified through SMDTStamp parameter, defaulting to firstTimestep. The reference position is calculated then by

(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 t0 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 vmin (SMDVmin) and the averaging time (SMDVminTave) have to be specified. Then every timesteps NAMD will compute the average velocity in the current direction over the past timesteps, and compare it to vmin. The average velocity is computed simply by

(5)

In the case that is less than vmin, a new direction is chosen, and the simulation proceeds.

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

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

where t is the current timestep, and Fmin (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 Fmin and the direction of $\vec n_{new}$. The time stamp t0 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
\begin{displaymath} p(\theta) \, = \, \frac{1}{\sigma \sqrt{2\pi}} \exp(- \frac{\theta^2}{\sigma^2}) \, . \end{displaymath} (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

\begin{displaymath} \vec n_{new} \; = \; \vec a \cos\varphi \sin\theta \, + \, \vec b \sin\varphi \sin\theta \, + \, \vec n_0 \cos\theta \, . \end{displaymath} (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 Fmin (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 vmax (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 vmax, the reference position is reset according to Eq. 6. The time stamp t0 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.

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:

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

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