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 , force constant k, velocity of the reference position movement and the number of the constrained atom i. Then during the simulation a potential USMD and a corresponding force are added for the specified atom. The potential is computed according to

 (1)

and the force acting on atom i is (by differentiating USMD)
 (2)

where t is time, is the position of the atom i and is the current reference position (position of the restraint point), given by
 (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 (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 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.

• 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 = , where 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 . Normal from the atom pos to the axis is, similarly, normal 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.,

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 in the pdb file, the force actually calculated is . 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 (SMDRefPos), the absolute value of the velocity of movement v (SMDVel), and the direction of movement (SMDDir). The velocity is then given by . Vector 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

 (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 the force applied to the atom has an absolute value of Fmin and the direction of . 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 , the angle of the cone which has 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 and defining the new direction. Angle is a random number between 0 and 360 degrees, and is a random number between 0 and . Gaussian method uses uniformly distributed random numbers to find and normally distributed random numbers with variance (SMDGaussW) to find . All values that exceed are ignored. The distribution of values is given by
 (7)

After the angles and are calculated, the new direction is found by constructing a right orthonormal system of coordinates with , the initial direction specified by SMDDir in the configuration file, being the z direction. The other two vectors and , orthogonal to and to each other are found by first producing , or in the case when is collinear to , normalizing , and, finally, producing . Then the new direction is given by

 (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 over a given period of time (SMDVmaxTave) exceeds the specified maximum allowed velocity vmax (SMDVmax) in the direction of reference position movement . The average velocity is calculated similarly to Eq. 5, and in the case when 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.

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

### 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''.

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]''.

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.

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

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