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] \; .$$ (5.38)

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}}$$ (5.39)

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

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

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 \fr... ...{2} \left(\frac{f-\bar{f}_i}{\sigma_i}\right)^2\right] \; .$$ (5.41)

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

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

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] \; .$$ (5.43)

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$}$$ (5.44)

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

$$B = \frac{\mbox{\tt rgauss2}}{M} \frac{M-\mbox{\tt rgauss1}} {\mbox{\tt rgauss2}-\mbox{\tt rgauss1}}$$ (5.46)

$$v = \frac{f - \bar{f}_i}{\sigma_i}$$ (5.47)

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

$$v' = v / F$$ (5.49)

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

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\} \; .$$ (5.50)

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 \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\} \; .$$ (5.51)

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

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

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]$$ (5.53)

$$\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] \; .$$ (5.54)

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.$$ (5.55)

where $\bar{f}$ is a lower bound and $p_{gauss}$ is given in Eq. 5.38. 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.$$ (5.56)

where $\bar{f}$ is an upper bound and $p_{gauss}$ is given in Eq. 5.38. 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$:

$$c = \vert b \vert - b \cos(n f + a)$$ (5.57)

where $b$ is CHARMM force constant, $a$ is phase shift (tested for 0 and 180${}^{o}$), 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$${}^{o}$. The force constant $b$ can be negative, in effect offsetting the phase $a$ for 180${}^{o}$ compared to the same but positive force constant.

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

Coulomb restraint

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

$$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.$$ (5.60)

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 RELATIVE_DIELECTRIC TOP 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}$$ (5.61)

$$\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.$$ (5.62)

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)$$ (5.63)

The parameters $f_1$ and $f_2$ of the switching function can be different from those in Eq. 5.60. 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}$$ (5.64)

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

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}$$ (5.66)

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

$$B = 1 - A$$ (5.68)

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

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

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+1}-f_j) c''_j + \frac{3B^2 - 1}{6} (f_{j+1}-f_j) c''_{j+1}$$ (5.71)

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.

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$$ (5.72)

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}}$$ (5.73)

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

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