Restraints and their derivatives

The chain rule is used to find the partial derivatives of the feature pdf with respect to the atomic coordinates. Thus, only the derivatives of the pdf with respect to the features are listed here.

Single Gaussian restraint

The pdf for a geometric feature $f$ (e.g., distance, angle, dihedral angle) is

$$p = \frac{1}{\sigma \sqrt{2 \pi}} \exp \left[-\frac{1}{2} \left(\frac{f-\bar{f}}{\sigma}\right)^2\right] \; .$$ (A.62)

A corresponding restraint $c$ in the sum that defines the objective function $F$ is

$$c = -\ln p = \frac{1}{2} \left(\frac{f-\bar{f}}{\sigma}\right)^2 - \ln \frac{1}{\sigma \sqrt{2 \pi}}$$ (A.63)

(Note that since the second term is constant for a given restraint, it is ignored. $c$ is also scaled by $RT$ in kcal/mol with $T = 297.15K$ to allow these scores to be summed with CHARMM energies.)

The first derivatives with respect to feature $f$ are:

$$\frac{ \; {d}c}{ \; {d}f} = \frac{f-\bar{f}}{\sigma} \; \frac{1}{\sigma} \; .$$ (A.64)

The relative heavy violation with respect to $f$ is given as:

$$\frac{f - \bar{f}}{\sigma}$$ (A.65)

Multiple Gaussian restraint

The polymodal pdf for a geometric feature $f$ (e.g., distance, angle, dihedral angle) is

$$p = \sum_{i=1}^n \omega_i p_i = \sum_{i=1}^n \omega_i \frac{1}{\s... ...\exp \left[-\frac{1}{2} \left(\frac{f-\bar{f}_i}{\sigma_i}\right)^2\right] \; .$$ (A.66)

A corresponding restraint $c$ in the sum that defines the objective function $F$ is (as before, this is scaled by $RT$):

$$c = -\ln p = -\ln \sum_{i=1}^n \omega_i p_i$$ (A.67)

The first derivatives with respect to feature $f$ are:

$$\frac{ \; {d}c}{ \; {d}f} = \frac{1}{p} \sum_{i=1}^n \omega_i p_i \cdot \left[\frac{f-\bar{f}_i}{\sigma_i}\frac{1}{\sigma_i}\right] \; .$$ (A.68)

When any of the normalized deviations $v_i = (f-\bar{f}_i) / \sigma_i$ is large, there are numerical instabilities in calculating the derivatives because $v_i$ are arguments to the exp function. Robustness is ensured as follows. The ‘effective’ normalized deviation is used in all the equations above when the magnitude of normalized violation $v$ is larger than cutoff rgauss1 (10 for double precision). This scheme works up to rgauss2 (200 for double precision); violations larger than that are ignored. This trick is equivalent to increasing the standard deviation $\sigma_i$. A slight disadvantage is that there is a discontinuity in the first derivatives at rgauss1. However, if continuity were imposed, the range would not be extended (this is equivalent to linearizing the Gaussian, but since it is already linear for large deviations, a linearization with derivatives smoothness would not introduce much change at all).

$$M = 37 \ \ \ ; \ \ \ \mbox{$M^2/2$\ has to be smaller than the largest argument to $\exp$}$$ (A.69)

$$A = \frac{1}{M} \frac{\mbox{\tt rgauss2}-M} {\mbox{\tt rgauss2} -\mbox{\tt rgauss1}}$$ (A.70)

$$B = \frac{\mbox{\tt rgauss2}}{M} \frac{M-\mbox{\tt rgauss1}} {\mbox{\tt rgauss2}-\mbox{\tt rgauss1}}$$ (A.71)

$$v = \frac{f - \bar{f}_i}{\sigma_i}$$ (A.72)

$$F = A \left\vert v \right\vert + B$$ (A.73)

$$v' = v / F$$ (A.74)

Now, Eqs. A.66-A.68 are used with $v'$ instead of $v$. For single precision, $M = 12$, rgauss1 = 4, rgauss2 = 100.

The relative heavy violation with respect to $f$ is given as:

$$\max_{\omega_i} (f-\bar{f}_i)/\sigma_i$$ (A.75)

Multiple binormal restraint

The polymodal pdf for a geometric feature $(f_1, f_2)$ (e.g., a pair of dihedral angles) is

$$p = \sum_{i=1}^n \omega_i p_i \; = \; \sum_{i=1}^n \omega_i \frac{1}{2 \pi \sigma_{1i} \sigma_{2i} \sqrt{(1-\rho_i^2)}} \; \cdot$$

$$\exp \left\{-\frac{1}{2 (1-\rho_i^2)} \left[ \left(\frac{f_1-\bar... ...i}} + \left(\frac{f_2-\bar{f}_{2i}}{\sigma_{2i}}\right)^2 \right] \right\} \; .$$ (A.76)

where $\rho < 1$. $\rho$ is the correlation coefficient between $f_1$ and $f_2$. MODELLER actually uses the following series expansion to calculate $p$:

$$p = \sum_{i=1}^n \omega_i \frac{1}{2 \pi \sigma_{1i} \sigma_{2i} \sqrt{(1-\rho_i^2)}} \; \cdot$$

$$\exp \left\{-\frac{1}{1-\rho_i^2} \left[ \frac{1-\cos(f_1-\bar{f}... ...ma_{2i}} + \frac{1-\cos(f_2-\bar{f}_{2i})}{\sigma_{2i}^2} \right] \right\} \; .$$ (A.77)

A corresponding restraint $c$ in the sum that defines the objective function $F$ is (as before, this is scaled by $RT$):

$$c = -\ln p = -\ln \sum_{i=1}^n \omega_i p_i$$ (A.78)

The first derivatives with respect to features $f_1$ and $f_2$ are:

$$\frac{\partial c}{\partial f_1} = \frac{1}{p} \sum_{i=1}^n \left[ \omega_i p_i \cdot \frac{1}{\sigm... ...frac{\cos(f_1-\bar{f}_{1i})\sin(f_2-\bar{f}_{2i})}{\sigma_{2i}} \right) \right]$$ (A.79)

$$\frac{\partial c}{\partial f_2} = \frac{1}{p} \sum_{i=1}^n \left[ \omega_i p_i \cdot \frac{1}{\sigm... ...\cos(f_2-\bar{f}_{2i})\sin(f_1-\bar{f}_{1i})}{\sigma_{1i}} \right) \right] \; .$$ (A.80)

The relative heavy violation with respect to $f$ is given as:

$$\max_{\omega_i} \sqrt{ -\frac{1}{2 (1-\rho_i^2)} \left[ \left(\fr... ...} {\sigma_{2i}} + \left(\frac{f_2-\bar{f}_{2i}}{\sigma_{2i}}\right)^2 \right] }$$ (A.81)

Lower bound

This is like the left half of a single Gaussian restraint:

$$p = \left\{ \begin{array}{ll} p_{gauss} \; ; & f < \bar{f} \\ 0 \; ; & f \geq \bar{f} \\ \end{array} \right.$$ (A.82)

where $\bar{f}$ is a lower bound and $p_{gauss}$ is given in Eq. A.62. A similar equation relying on the first derivatives of a Gaussian $p$ holds for the first derivatives of a lower bound.

Upper bound

This is like the right half of a single Gaussian restraint:

$$p = \left\{ \begin{array}{ll} p_{gauss} \; ; & f > \bar{f} \\ 0 \; ; & f \leq \bar{f} \\ \end{array} \right.$$ (A.83)

where $\bar{f}$ is an upper bound and $p_{gauss}$ is given in Eq. A.62. A similar equation relying on the first derivatives of a Gaussian $p$ holds for the first derivatives of an upper bound.

Cosine restraint

This is usually used for dihedral angles $f$ (improper dihedrals generally use a Gaussian restraint instead):

$$c = \vert b \vert - b \cos(n f + a)$$ (A.84)

where $b$ is CHARMM force constant, $a$ is phase shift (tested for 0 and 180°), and $n$ is periodicity (tested for 1, 2, 3, 4, 5, and 6). The CHARMM phase value from the CHARMM parameter library corresponds to $a-180$°. The force constant $b$ can be negative, in effect offsetting the phase $a$ for 180° compared to the same but positive force constant.

$$\frac{ \; {d}c}{ \; {d}f} = b n \sin(n f + a)$$ (A.85)

Coulomb restraint

$$c = \frac{1}{\epsilon_r}\frac{q_i q_j}{f} s(f,f_1,f_2)$$ (A.86)

$$s(f,f_1,f_2) = \left\{ \begin{array}{ll} 1 \; ; & f \leq f_1 \\ \frac{(f_2 - f)... ...f_2-f_1)^3}\; ; & f_o < f \leq f_2 \\ 0 \; ; & f > f_2 \\ \end{array} \right.$$ (A.87)

where $q_i$ and $q_j$ are the atomic charges of atoms $i$ and $j$, obtained from the CHARMM topology file, that are at a distance $f$. $\epsilon_r$ is the relative dielectric, controlled by the EnergyData.relative_dielectric variable. Function $s(f,f_1,f_2)$ is a switching function that smoothes the potential down to zero in the interval from $f_1$ to $f_2$ ($f_2> f_1$). The total Coulomb energy of a molecule is a sum over all pairs of atoms that are not in the same bonds or bond angles. 1-4 energy for the 1-4 atom pairs in the same dihedral angle corresponds to the ELEC14 MODELLER term; the remaining longer-range contribution corresponds to the ELEC term.

The first derivatives are:

$$\frac{ \; {d}c}{ \; {d}f} = -\frac{c}{f} + \frac{c}{s} \frac{ \; {d}s}{ \; {d}f}$$ (A.88)

$$\frac{ \; {d}s}{ \; {d}f} = \left\{ \begin{array}{ll} 0 \; ; & f \leq f_1 \\ \frac{6 (f_2-f)... ..._2-f_1)^3} \; ; & f_1 < f \leq f_2 \\ 0 \; ; & f > f_2 \\ \end{array} \right.$$ (A.89)

The violations of this restraint are always reported as zero.

Lennard-Jones restraint

Usually used for non-bonded distances:

$$c = \left[\left(\frac{A}{f}\right)^{12} - \left(\frac{B}{f}\right)^6 \right] s(f,f_1,f_2)$$ (A.90)

The parameters $f_1$ and $f_2$ of the switching function can be different from those in Eq. A.87. The parameters $A$ and $B$ are obtained from the CHARMM parameter file (NONBOND section) where they are given as $E_i$ and $r_j$ such that $E_{ij}(f) = -4 \sqrt{E_i E_j} [(\rho_{ij}/f)^{12} - (\rho_{ij}/f)^6]$ in kcal/mole for $f$ in angstroms and $\rho = (r_i + r_j)/2^{1/6}$; the minimum of $E$ is $-\sqrt{E_i E_j}$ at $f=(r_i + r_j)$, and its zero is at $f=\rho$. The total Lennard-Jones energy should be evaluated over all pairs of atoms that are not in the same bonds or bond angles. The parameters $A$ and $B$ for 1-4 pairs in dihedral angles can be different from those for the other pairs; they are obtained from the second set of $E_i$ and $r_i$ in the CHARMM parameter file, if it exists. 1-4 energy corresponds to the LJ14 MODELLER term; the remaining longer-range contribution corresponds to the LJ term.

The first derivatives are:

$$\frac{ \; {d}c}{ \; {d}f} = \frac{C s}{f} - C \frac{ \; {d}s}{ \; {d}f}$$ (A.91)

$$C = -12 \left(\frac{A}{f}\right)^{12} + 6 \left(\frac{B}{f}\right)^6$$ (A.92)

As $f$ tends toward zero, the repulsive part of the energy dominates, and approaches infinity. Near-infinite forces result in unstable trajectories during optimization. This is particularly a problem in the first few steps of optimization starting from randomized, interpolated, or otherwise non-physical atomic coordinates. To avoid this, the potential is simply artificially truncated: if $A/f$ exceeds 6, $f$ is treated as being equal to $A/6$.

The violations of this restraint are always reported as zero.

Spline restraint

Any restraint form can be represented by a cubic spline [Press et al., 1992]:

$$c = A c_j + B c_{j+1} + C c''_j + D c''_{j+1}$$ (A.93)

$$A = \frac{f_{j+1} - f}{f_{j+1} - f_j}$$ (A.94)

$$B = 1 - A$$ (A.95)

$$C = \frac{1}{6}(A^3-A)(f_{j+1} - f_j)^2$$ (A.96)

$$D = \frac{1}{6}(B^3-B)(f_{j+1} - f_j)^2$$ (A.97)

where $f_j \le f \le f_{j+1}$.

The first derivatives are:

$$\frac{ \; {d}c}{ \; {d}f} = \frac{c_{j+1} - c_j}{f_{j+1} - f_j} -... ...c{3A^2 - 1}{6} (f_{j+1}-f_j) c''_j + \frac{3B^2 - 1}{6} (f_{j+1}-f_j) c''_{j+1}$$ (A.98)

The values of $c$ and $c'$ beyond $f_1$ and $f_n$ are obtained by linear interpolation from the termini. A violation of the restraint is calculated by finding the global minimum. A relative violation is estimated by using a standard deviation (e.g., force constant) obtained by fitting a parabola to the global minimum.

Variable spacing of spline points could be used to save on memory. However, this would increase the execution time, so it is not used.

To calculate the relative heavy violation, the feature value $\bar{f}$ that results in the smallest value of the restraint is obtained by interpolation, and a Gaussian function is fitted locally around this value to obtain the standard deviation $\sigma$. These are then used in Eq. A.65.

Symmetry restraint

The asymmetry penalty added to the objective function is defined as

$$F_{symm} = \sum_{i<j} \omega_i \omega_j (d_{ij} - d'_{ij})^2$$ (A.99)

where the sum runs over all pairs of equivalent atoms $ij$, $\omega_i$ is an atom weight for atom $i$, $d_{ij}$ is an intra-molecular distance between atoms $ij$ in the first segment, and $d'_{ij}$ is the equivalent distance in the second segment.

For each $i<j$, the first derivatives are:

$$\frac{\partial c}{\partial \vec{d}_{ij}} = 2 \omega_i \omega_j (d_{ij} - d'_{ij}) \frac{\vec{d}_{ij}}{d_{ij}}$$ (A.100)

$$\frac{\partial c}{\partial \vec{d}'_{ij}} = - 2 \omega_i \omega_j (d_{ij} - d'_{ij}) \frac{\vec{d}'_{ij}}{d'_{ij}}$$ (A.101)

Thus, the total first derivatives are obtained by summing the two expressions above for all $i$ and $j>i$ distances.