Harmonic restraints

The harmonic biasing method may be used to enforce fixed or moving restraints, including variants of Steered and Targeted MD. Within energy minimization runs, it allows for restrained minimization, e.g. to calculate relaxed potential energy surfaces. In the context of the Colvars module, harmonic potentials are meant according to their textbook definition:

$$V(\xi) = \frac{1}{2} k \left(\frac{\xi - \xi_0}{w_{\xi}}\right)^2$$ (13.34)

There are two noteworthy aspects of this expression:

1. Because the standard coefficient of $1/2$ of the harmonic potential is included, this expression differs from harmonic bond and angle potentials historically used in common force fields, where the factor was typically omitted resulting in a non-standard definition of the force constant.
2. The variable $\xi$ is not only centered at $\xi_0$ , but is also scaled by its characteristic length scale $w_{\xi}$ (keyword width). The resulting dimensionless variable $z = (\xi - \xi_0)/w_{\xi}$ is typically easier to treat numerically: for example, when the forces typically experienced by $\xi$ are much smaller than $k/w_{\xi}$ and $k$ is chosen equal to $\kappa_{\mathrm{B}}T$ (thermal energy), the resulting probability distribution of $z$ is approximately a Gaussian with mean equal to 0 and standard deviation equal to 1.

   This property can be used for setting the force constant in umbrella-sampling ensemble runs: if the restraint centers are chosen in increments of $w_{\xi}$ , the resulting distributions of $\xi$ are most often optimally overlapped. In regions where the underlying free-energy landscape induces highly skewed distributions of $\xi$ , additional windows may be added as needed, with spacings finer than $w_{\xi}$ .

Beyond one dimension, the use of a scaled harmonic potential also allows a standard definition of a multi-dimensional restraint with a unified force constant:

$$V(\xi_{1}, \ldots, \xi_{M}) = \frac{1}{2} k \sum_{i=1}^{M} \left(\frac{\xi_{i} - \xi_0}{w_{\xi}}\right)^2$$ (13.35)

If one-dimensional or homogeneous multi-dimensional restraints are defined, and there are no other uses for the parameter $w_{\xi}$ , width can be left at its default value of $1$ .

A harmonic restraint is defined by a harmonic {...} block, which may contain the following keywords:

• name: see definition of name (biasing and analysis methods)
• colvars: see definition of colvars (biasing and analysis methods)
• outputEnergy: see definition of outputEnergy (biasing and analysis methods)
• writeTIPMF: see definition of writeTIPMF (biasing and analysis methods)
• writeTISamples: see definition of writeTISamples (biasing and analysis methods)
• stepZeroData: see definition of stepZeroData (biasing and analysis methods)

• forceConstant $\langle\,$ Scaled force constant (kcal/mol)$\,\rangle$
  Context: harmonic
  Acceptable values: positive decimal
  Default value: 1.0
  Description: This option defines a scaled force constant $k$ for the harmonic potential (eq. 13.36). To ensure consistency for multidimensional restraints, it is divided internally by the square of the specific width of each variable (which is 1 by default). This makes all values effectively dimensionless and of commensurate size. For instance, if this force constant is set to the thermal energy $\kappa_{\mathrm{B}}T$ (equal to $RT$ if molar units are used), then the amplitude of the thermal fluctuations of each variable $\xi$ will be on the order of its width, $w_{\xi}$ . This can be used to estimate the optimal spacing of umbrella-sampling windows (under the assumption that the force constant is larger than the curvature of the underlying free energy). The values of the actual force constants $k/w_{\xi}^2$ are always printed when the restraint is defined.

• centers $\langle\,$ Initial harmonic restraint centers$\,\rangle$
  Context: harmonic
  Acceptable values: space-separated list of colvar values
  Description: The centers (equilibrium values) of the restraint, $\xi_{0}$ , are entered here. The number of values must be the number of requested colvars. Each value is a decimal number if the corresponding colvar returns a scalar, a ``(x, y, z)'' triplet if it returns a unit vector or a vector, and a ``(q0, q1, q2, q3)'' quadruplet if it returns a rotational quaternion. If a colvar has periodicities or symmetries, its closest image to the restraint center is considered when calculating the harmonic potential.

Tip: A complex set of restraints can be applied to a system, by defining several colvars, and applying one or more harmonic restraints to different groups of colvars. In some cases, dozens of colvars can be defined, but their value may not be relevant: to limit the size of the colvars trajectory file, it may be wise to disable outputValue for such ``ancillary'' variables, and leave it enabled only for ``relevant'' ones.

Moving restraints: steered molecular dynamics

The following options allow to change gradually the centers of the harmonic restraints during a simulations. When the centers are changed continuously, a steered MD in a collective variable space is carried out.

• targetCenters $\langle\,$ Steer the restraint centers towards these targets$\,\rangle$
  Context: harmonic
  Acceptable values: space-separated list of colvar values
  Description: When defined, the current centers will be moved towards these values during the simulation. By default, the centers are moved over a total of targetNumSteps steps by a linear interpolation, in the spirit of Steered MD. If targetNumStages is set to a nonzero value, the change is performed in discrete stages, lasting targetNumSteps steps each. This second mode may be used to sample successive windows in the context of an Umbrella Sampling simulation. When continuing a simulation run, the centers specified in the configuration file $<$ colvarsConfig$>$ are overridden by those saved in the restart file $<$ colvarsInput$>$ . To perform Steered MD in an arbitrary space of colvars, it is sufficient to use this option and enable outputAccumulatedWork and/or outputAppliedForce within each of the colvars involved.

• targetNumSteps $\langle\,$ Number of steps for steering$\,\rangle$
  Context: harmonic
  Acceptable values: positive integer
  Description: In single-stage (continuous) transformations, defines the number of MD steps required to move the restraint centers (or force constant) towards the values specified with targetCenters or targetForceConstant. After the target values have been reached, the centers (resp. force constant) are kept fixed. In multi-stage transformations, this sets the number of MD steps per stage.

• outputCenters $\langle\,$ Write the current centers to the trajectory file$\,\rangle$
  Context: harmonic
  Acceptable values: boolean
  Default value: off
  Description: If this option is chosen and colvarsTrajFrequency is not zero, the positions of the restraint centers will be written to the trajectory file during the simulation. This option allows to conveniently extract the PMF from the colvars trajectory files in a steered MD calculation.

Note on restarting moving restraint simulations: Information about the current step and stage of a simulation with moving restraints is stored in the restart file (state file). Thus, such simulations can be run in several chunks, and restarted directly using the same colvars configuration file. In case of a restart, the values of parameters such as targetCenters, targetNumSteps, etc. should not be changed manually.

Moving restraints: umbrella sampling

The centers of the harmonic restraints can also be changed in discrete stages: in this cases a one-dimensional umbrella sampling simulation is performed. The sampling windows in simulation are calculated in sequence. The colvars trajectory file may then be used both to evaluate the correlation times between consecutive windows, and to calculate the frequency distribution of the colvar of interest in each window. Furthermore, frequency distributions on a predefined grid can be automatically obtained by using the histogram bias (see ).

To activate an umbrella sampling simulation, the same keywords as in the previous section can be used, with the addition of the following:

• targetNumStages $\langle\,$ Number of stages for steering$\,\rangle$
  Context: harmonic
  Acceptable values: non-negative integer
  Default value: 0
  Description: If non-zero, sets the number of stages in which the restraint centers or force constant are changed to their target values. If zero, the change is continuous. Each stage lasts targetNumSteps MD steps. To sample both ends of the transformation, the simulation should be run for targetNumSteps $\times$ (targetNumStages + 1).

Changing force constant

The force constant of the harmonic restraint may also be changed to equilibrate [76].

• targetForceConstant $\langle\,$ Change the force constant towards this value$\,\rangle$
  Context: harmonic
  Acceptable values: positive decimal
  Description: When defined, the current forceConstant will be moved towards this value during the simulation. Time evolution of the force constant is dictated by the targetForceExponent parameter (see below). By default, the force constant is changed smoothly over a total of targetNumSteps steps. This is useful to introduce or remove restraints in a progressive manner. If targetNumStages is set to a nonzero value, the change is performed in discrete stages, lasting targetNumSteps steps each. This second mode may be used to compute the conformational free energy change associated with the restraint, within the FEP or TI formalisms. For convenience, the code provides an estimate of the free energy derivative for use in TI. A more complete free energy calculation (particularly with regard to convergence analysis), while not handled by the Colvars module, can be performed by post-processing the colvars trajectory, if colvarsTrajFrequency is set to a suitably small value. It should be noted, however, that restraint free energy calculations may be handled more efficiently by an indirect route, through the determination of a PMF for the restrained coordinate.[76]

• targetForceExponent $\langle\,$ Exponent in the time-dependence of the force constant$\,\rangle$
  Context: harmonic
  Acceptable values: decimal equal to or greater than 1.0
  Default value: 1.0
  Description: Sets the exponent, $\alpha$ , in the function used to vary the force constant as a function of time. The force is varied according to a coupling parameter $\lambda$ , raised to the power $\alpha$ : $k_\lambda = k_0 + \lambda^\alpha (k_1 - k_0)$ , where $k_0$ , $k_\lambda$ , and $k_1$ are the initial, current, and final values of the force constant. The parameter $\lambda$ evolves linearly from 0 to 1, either smoothly, or in targetNumStages equally spaced discrete stages, or according to an arbitrary schedule set with lambdaSchedule. When the initial value of the force constant is zero, an exponent greater than 1.0 distributes the effects of introducing the restraint more smoothly over time than a linear dependence, and ensures that there is no singularity in the derivative of the restraint free energy with respect to lambda. A value of 4 has been found to give good results in some tests.

• targetEquilSteps $\langle\,$ Number of steps discarded from TI estimate$\,\rangle$
  Context: harmonic
  Acceptable values: positive integer
  Description: Defines the number of steps within each stage that are considered equilibration and discarded from the restraint free energy derivative estimate reported reported in the output.

• lambdaSchedule $\langle\,$ Schedule of lambda-points for changing force constant$\,\rangle$
  Context: harmonic
  Acceptable values: list of real numbers between 0 and 1
  Description: If specified together with targetForceConstant, sets the sequence of discrete $\lambda$ values that will be used for different stages.