\subsection{Features and their derivatives}

\subsubsection{Distance}

Distance is defined by points $i$ and $j$:
\begin{equation}
d = \sqrt{\vec{r}_{ij} \cdot \vec{r}_{ij}} = |\vec{r}_{ij}| = r_{ij}
\end{equation}
where
\begin{equation}
  \vec{r}_{ij} = \vec{r}_i - \vec{r}_j \; .
\end{equation}

The first derivatives of $d$ with respect to Cartesian coordinates are:
\begin{eqnarray}
  \frac{\partial d} {\partial \vec{r}_i} & = & 
      \frac{\vec{r}_{ij}}{|\vec{r}_{ij}|} \\[3mm]
  \frac{\partial d} {\partial \vec{r}_j} & = &  
     -\frac{\partial d} {\partial \vec{r}_i}
\end{eqnarray}


\subsubsection{Angle}

Angle is defined by points $i$, $j$, and $k$, and spanned by vectors 
$ij$ and $kj$:
\begin{equation}
\alpha = \arccos \frac{\vec{r}_{ij} \cdot \vec{r}_{kj}}
                      {|\vec{r}_{ij}| |\vec{r}_{kj}|} \; .
\end{equation}
It lies in the interval from 0 to 180\degr. Internal 
\MODELLER\ units are radians.

The first derivatives of $\alpha$ with respect to Cartesian coordinates are:
\begin{eqnarray}
  \frac{\partial \alpha}{\partial \vec{r}_i} & = & 
      \frac{\partial \alpha}{\partial \cos \alpha} \; 
      \frac{\partial \cos \alpha}{\partial \vec{r}_i} = 
      \frac{1}{\sqrt{1-\cos^2 \alpha}} \;
      \frac{1}{r_{ij}} \;
      \left( \frac{\vec{r}_{ij}}{r_{ij}} \cos \alpha -
             \frac{\vec{r}_{kj}}{r_{kj}} \right) \\[4mm]
  \frac{\partial \alpha}{\partial \vec{r}_k} & = & 
      \frac{\partial \alpha}{\partial \cos \alpha} \; 
      \frac{\partial \cos \alpha}{\partial \vec{r}_k} = 
      \frac{1}{\sqrt{1-\cos^2 \alpha}} \;
      \frac{1}{r_{kj}} \;
      \left( \frac{\vec{r}_{kj}}{r_{kj}} \cos \alpha -
             \frac{\vec{r}_{ij}}{r_{ij}} \right) \\[4mm]
  \frac{\partial \alpha} {\partial \vec{r}_j} & = &  
     -\frac{\partial d} {\partial \vec{r}_i}
     -\frac{\partial d} {\partial \vec{r}_k}
\end{eqnarray}

These equations for the  derivatives have a numerical instability when 
the angle goes to 0 or to 180\degr.
Presently, the problem is `solved' by testing for the size
of the angle; if it is too small, the derivatives are set to 0
in the hope that other restraints will eventually pull the angle
towards well behaved regions. Thus, angle restraints of 0 or 
180\degr\ should not be used in the conjugate gradients or molecular dynamics
optimizations.

\subsubsection{Dihedral angle}

Dihedral angle is defined by points $i$, $j$, $k$, and $l$ ($ijkl$):
\begin{equation}
  \chi = \mbox{sign}(\chi) \arccos\frac{
                           (\vec{r}_{ij} \times \vec{r}_{kj}) \cdot 
                           (\vec{r}_{kj} \times \vec{r}_{kl})}
                          {|\vec{r}_{ij} \times \vec{r}_{kj}|
                           |\vec{r}_{kj} \times \vec{r}_{kl}|}
\end{equation}
where
\begin{equation}
  \mbox{sign}(\chi) = \mbox{sign}[\vec{r}_{kj} \cdot 
                        (\vec{r}_{ij} \times \vec{r}_{kj}) \times
                        (\vec{r}_{kj} \times \vec{r}_{kl})] \; .
\end{equation}

The first derivatives of $\chi$ with respect to Cartesian coordinates are:
\begin{equation}
  \frac{\dd \chi} {\dd \vec{r}} = 
    \frac{\dd \chi}{\dd \cos \chi} \frac{\dd \cos \chi}{\dd \vec{r}}
\end{equation}
where
\begin{equation}
  \frac{\dd \chi}{\dd \cos \chi} = \left(\frac{\dd \cos \chi}{\dd \chi}\right)^{-1} = 
                                   -\frac{1}{\sin \chi}
\end{equation}
and
\begin{eqnarray}
  \frac{\partial \cos \chi}{\partial \vec{r}_i} & = & \vec{r}_{kj} \times \vec{a} \\[3mm]
  \frac{\partial \cos \chi}{\partial \vec{r}_j} & = & \vec{r}_{ik} \times \vec{a} -
                                        \vec{r}_{kl} \times \vec{b} \\[3mm]
  \frac{\partial \cos \chi}{\partial \vec{r}_k} & = & \vec{r}_{jl} \times \vec{b} -
                                        \vec{r}_{ij} \times \vec{a} \\[3mm]
  \frac{\partial \cos \chi}{\partial \vec{r}_l} & = & \vec{r}_{ij} \times \vec{b} \\[3mm]
  \vec{a} & = & \frac{1}{|\vec{r}_{ij} \times \vec{r}_{kj}|} \;
                     \left(\frac{\vec{r}_{kj} \times \vec{r}_{kl}}
                     {|\vec{r}_{kj} \times \vec{r}_{kl}|} - 
                     \cos \chi  \;
                     \frac{\vec{r}_{ij} \times \vec{r}_{kj}}
                          {|\vec{r}_{ij} \times \vec{r}_{kj}|}\right) \\[6mm]
  \vec{b} & = & \frac{1}{|\vec{r}_{kj} \times \vec{r}_{kl}|} \;
                     \left(\frac{\vec{r}_{ij} \times \vec{r}_{kj}}
                     {|\vec{r}_{ij} \times \vec{r}_{kj}|} - 
                     \cos \chi  \;
                     \frac{\vec{r}_{kj} \times \vec{r}_{kl}}
                          {|\vec{r}_{kj} \times \vec{r}_{kl}|}\right) \; .
\end{eqnarray}

These equations for the derivatives have a numerical instability when 
the angle goes to 0. Thus, the following set of equations is used
instead \cite{Sch93}:

\begin{eqnarray}
\vec{r}_{mj} & = & \vec{r}_{ij} \times \vec{r}_{kj} \\
\vec{r}_{nk} & = & \vec{r}_{kj} \times \vec{r}_{kl} \\
\frac{\partial \chi}{\partial \vec{r}_i} & = &
      \frac{r_{kj}}{r^2_{mj}} \vec{r}_{mj} \\
\frac{\partial \chi}{\partial \vec{r}_l} & = &
      -\frac{r_{kj}}{r^2_{nk}} \vec{r}_{nk} \\
\frac{\partial \chi}{\partial \vec{r}_j} & = &
      \left(\frac{\vec{r}_{ij} \cdot \vec{r}_{kj}}{r^2_{kj}} - 1 \right) \;
      \frac{\partial \chi}{\partial \vec{r}_i} - 
      \frac{\vec{r}_{kl} \cdot \vec{r}_{kj}}{r^2_{kj}}
      \frac{\partial \chi}{\partial \vec{r}_l} \\
\frac{\partial \chi}{\partial \vec{r}_k} & = &
      \left(\frac{\vec{r}_{kl} \cdot \vec{r}_{kj}}{r^2_{kj}} - 1 \right) \;
      \frac{\partial \chi}{\partial \vec{r}_l} -
      \frac{\vec{r}_{ij} \cdot \vec{r}_{kj}}{r^2_{kj}}
      \frac{\partial \chi}{\partial \vec{r}_i}
\end{eqnarray}

The only possible instability in these equations is when the length of
the central bond of the dihedral, $r_{kj}$, goes to 0. In such a case,
which should not happen, the derivatives are set to 0. The expressions for
an improper dihedral angle, as opposed to a dihedral or dihedral angle,
are the same, except that indices $ijkl$ are permuted to $ikjl$.
In both cases, covalent bonds $ij$, $jk$, and $kl$ are defining
the angle.

\subsubsection{Atomic solvent accessibility}

xx

\subsubsection{Atomic density}

xx

\subsection{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.

\subsubsection{Single Gaussian restraint}

The pdf for a geometric feature $f$ (\eg, distance, angle, 
dihedral angle) is
\begin{equation}
  p = \frac{1}{\sigma \sqrt{2 \pi}} \exp \left[-\frac{1}{2}
               \left(\frac{f-\bar{f}}{\sigma}\right)^2\right] \; .
  \label{eq:gauss}
\end{equation}
A corresponding restraint $c$ in the sum that defines the objective 
function $F$ is
\begin{equation}
  c = -\ln p = \frac{1}{2} \left(\frac{f-\bar{f}}{\sigma}\right)^2 -
               \ln \frac{1}{\sigma \sqrt{2 \pi}}
  \label{eq:gaussian}
\end{equation}

The first derivatives with respect to feature $f$ are:
\begin{equation}
  \frac{\dd c}{\dd f} = \frac{f-\bar{f}}{\sigma} \; \frac{1}{\sigma} \; .
\end{equation}



\subsubsection{Multiple Gaussian restraint}

The polymodal pdf for a geometric feature $f$ (\eg, distance, angle, 
dihedral angle) is
\begin{equation}
  p = \sum_{i=1}^n \omega_i p_i = \sum_{i=1}^n \omega_i 
      \frac{1}{\sigma_i \sqrt{2 \pi}} \exp \left[-\frac{1}{2}
               \left(\frac{f-\bar{f}_i}{\sigma_i}\right)^2\right] \; .
  \label{eq:polygauss}
\end{equation}
A corresponding restraint $c$ in the sum that defines the objective 
function $F$ is
\begin{equation}
  c = -\ln p = -\ln \sum_{i=1}^n \omega_i p_i
\end{equation}

The first derivatives with respect to feature $f$ are:
\begin{eqnarray}
  \frac{\dd c}{\dd 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]
%      = \nonumber \\ & &
%      \frac{1}{p} \sum_{i=1}^n \left\{ \omega_i
%      \frac{1}{\sigma_i \sqrt{2 \pi}} \exp \left[-\frac{1}{2}
%               \left(\frac{f-\bar{f}_i}{\sigma_i}\right)^2\right] \cdot
%      \left[\frac{f-\bar{f}_i}{\sigma_i}\frac{1}{\sigma_i}\right] \right\} 
       \; .
  \label{eq:dpolygauss}
\end{eqnarray}

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 {\it 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 {\tt rgauss1} (10 for double precision). This scheme works up 
to {\tt 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 {\tt 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).

\begin{eqnarray}
 M  & = & 37 \ \ \ ; \ \ \ \ 
    \mbox{$M^2/2$ has to be smaller than the largest argument to $\exp$} \\
 A  & = &  \frac{1}{M} \frac{\mbox{\tt rgauss2}-M}
                            {\mbox{\tt rgauss2} -\mbox{\tt rgauss1}}\\
 B  & = &  \frac{\mbox{\tt rgauss2}}{M} \frac{M-\mbox{\tt rgauss1}}
                                       {\mbox{\tt rgauss2}-\mbox{\tt rgauss1}}\\
 v  & = & \frac{f - \bar{f}_i}{\sigma_i} \\
 F  & = & A \left| v \right| + B \\
 v' & = & v / F
\end{eqnarray}

Now, Eqs.~\ref{eq:polygauss}--\ref{eq:dpolygauss} are used with $v'$ instead 
of $v$. For single precision, $M = 12$, {\tt rgauss1} = 4, {\tt rgauss2} = 100.


\subsubsection{Multiple binormal restraint}

The polymodal pdf for a geometric feature $(f_1, f_2)$ (\eg, a pair of
dihedral angles) is

\begin{eqnarray}
  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 
      \nonumber \\
      & & \exp \left\{-\frac{1}{2 (1-\rho_i^2)}
               \left[
               \left(\frac{f_1-\bar{f}_{1i}}{\sigma_{1i}}\right)^2 - 
               2 \rho_i \frac{f_1-\bar{f}_{1i}} {\sigma_{1i}}
                        \frac{f_2-\bar{f}_{2i}} {\sigma_{2i}} +
               \left(\frac{f_2-\bar{f}_{2i}}{\sigma_{2i}}\right)^2
              \right] 
              \right\} \; .
  \label{eq:poly2gauss}
\end{eqnarray}
where $\rho < 1$. $\rho$ is the correlation coefficient between $f_1$ and $f_2$.
A corresponding restraint $c$ in the sum that defines the objective 
function $F$ is
\begin{equation}
  c = -\ln p = -\ln \sum_{i=1}^n \omega_i p_i
\end{equation}

The first derivatives with respect to features $f_1$ and $f_2$ are:
\begin{eqnarray}
  \frac{\partial c}{\partial f_1} & = & \frac{1}{p}
      \sum_{i=1}^n \left[ \omega_i p_i \cdot
      \frac{1}{\sigma_{1i} (1-\rho_i^2)}
      \left(\frac{f_1-\bar{f}_{1i}}{\sigma_{1i}} -
            \rho_i\frac{f_2-\bar{f}_{2i}}{\sigma_{2i}}
      \right)
      \right] \\[5mm]
  \frac{\partial c}{\partial f_2} & = & \frac{1}{p}
      \sum_{i=1}^n \left[ \omega_i p_i \cdot
      \frac{1}{\sigma_{2i} (1-\rho_i^2)}
      \left(\frac{f_2-\bar{f}_{2i}}{\sigma_{2i}} -
            \rho_i\frac{f_1-\bar{f}_{1i}}{\sigma_{1i}}
      \right)
      \right] \; . 
  \label{eq:dpoly2gauss}
\end{eqnarray}


\subsubsection{Lower bound}

This is like the left half of a single Gaussian restraint:
\begin{equation}
  p = \left\{ \begin{array}{ll} p_{gauss} \; ; & f < \bar{f} \\ 
                                0 \; ;         & f \geq \bar{f} \\ 
                \end{array} \right.
  \label{eq:lgaussian}
\end{equation}
where $\bar{f}$ is a lower bound and $p_{gauss}$ is given in 
Eq.~\ref{eq:gauss}. 
A similar equation relying on the first derivatives of a Gaussian $p$ 
holds for the first derivatives of a lower bound.


\subsubsection{Upper bound}

This is like the right half of a single Gaussian restraint:
\begin{equation}
  p = \left\{ \begin{array}{ll} p_{gauss} \; ; & f > \bar{f} \\ 
                                0 \; ;         & f \leq \bar{f} \\ 
               \end{array} \right.
  \label{eq:rgaussian}
\end{equation}
where $\bar{f}$ is an upper bound and $p_{gauss}$ is given 
in Eq.~\ref{eq:gauss}. 
A similar equation relying on the first derivatives of a Gaussian $p$ 
holds for the first derivatives of an upper bound.


\subsubsection{Cosine restraint}

This is usually used for dihedral angles $f$:

\begin{equation}
c = | b | - b \cos(n f + a)
\label{eq:cosine}
\end{equation}
where $b$ is \CHARMM\ force constant, $a$ is phase shift (tested for
0 and 180\degr), 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$\degr. The force constant $b$ can be negative, in effect offsetting
the phase $a$ for 180\degr\ compared to the same but positive force constant.

\begin{equation}
\frac{\dd c}{\dd f} = b n \sin(n f + a)
\end{equation}

\subsubsection{Coulomb restraint}

\begin{eqnarray}
c & = & \frac{q_i q_j}{f} s(f,f_1,f_2) \\
s(f,f_1,f_2) & = & \left\{ \begin{array}{ll} 1 \; ;
                                          & f \leq f_1 \\ 
                                  \frac{(f_2 - f)^2 (f_2+2f-3f_1)}
                                       {(f_2-f_1)^3}\; ;
                                          & f_o  < f \leq f_2  \\
                                  0 \; ;
                                          & f > f_2 \\ 
                \end{array} \right.
\label{eq:coulomb}
\end{eqnarray}
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$. 
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 \V{ELEC14} \MODELLER\ term; the remaining 
longer-range contribution corresponds to the \V{ELEC} term.

The first derivatives are:
\begin{eqnarray}
\frac{\dd c}{\dd f} & = & -\frac{c}{f} + \frac{c}{s} \frac{\dd s}{\dd f} \\
\frac{\dd s}{\dd f} & = & \left\{ \begin{array}{ll} 0 \; ;
                                          & f \leq f_1 \\ 
                                  \frac{6 (f_2-f)(f_1-f)}{(f_2-f_1)^3} \; ;
                                          & f_1  < f \leq f_2  \\
                                  0 \; ;
                                          & f > f_2 \\ 
                \end{array} \right.
\end{eqnarray}


\subsubsection{Lennard-Jones restraint}

Usually used for non-bonded distances:
\begin{equation}
  c = \left[(\frac{A}{f})^{12} - (\frac{B}{f})^6 \right] s(f,f_1,f_2)
 \label{eq:lennard}
\end{equation}
The parameters $f_1$ and $f_2$ of the switching function can be
different from those in Eq.~\ref{eq:coulomb}. 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 \V{LJ14} \MODELLER\ term; the
remaining longer-range contribution corresponds to the \V{LJ} term.

The first derivatives are:
\begin{eqnarray}
\frac{\dd c}{\dd f} & = & \frac{C s}{f} - C \frac{\dd s}{\dd f} \\
C & = & -12 (\frac{A}{f})^{12} + 6 (\frac{B}{f})^6
\end{eqnarray}


\subsubsection{Spline restraint}

Any restraint form can be represented by a cubic spline \cite{PreTeuVet92}:

\begin{eqnarray}
  c & = & A c_j + B c_{j+1} + C c''_j + D c''_{j+1} \\
  A & = & \frac{f_{j+1} - f}{f_{j+1} - f_j} \\
  B & = & 1 - A \\
  C & = & \frac{1}{6}(A^3-A)(f_{j+1} - f_j)^2 \\
  D & = & \frac{1}{6}(B^3-B)(f_{j+1} - f_j)^2
  \label{eq:spline}
\end{eqnarray}
where $f_j \le f \le f_{j+1}$.

The first derivatives are:
\begin{equation}
\frac{\dd c}{\dd f} = \frac{c_{j+1} - c_j}{f_{j+1} - f_j} -
                      \frac{3A^2 - 1}{6} (f_{j+1}-f_j) c''_j +
                      \frac{3B^2 - 1}{6} (f_{j+1}-f_j) c''_{j+1}
\end{equation}

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 (\eg, 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.


\subsubsection{Symmetry restraint}

The asymmetry penalty added to the objective function is defined as

\begin{equation}
F_{symm} = \sum_{i<j} \omega_i \omega_j (d_{ij} - d'_{ij})^2
\label{eq:symmetry}
\end{equation}
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:
\begin{eqnarray}
\frac{\partial c}{\partial \vec{d}_{ij}} & = &
   2 \omega_i \omega_j (d_{ij} - d'_{ij})
     \frac{\vec{d}_{ij}}{d_{ij}} \\
\frac{\partial c}{\partial \vec{d}'_{ij}} & = & -
   2 \omega_i \omega_j (d_{ij} - d'_{ij})
     \frac{\vec{d}'_{ij}}{d'_{ij}}
\end{eqnarray}
Thus, the total first derivatives are obtained by summing the two expressions
above for all $i$ and $j>i$ distances.
