The MMFF potential expression is:

\[V_{MMFF} = \sum_{b}V_{b} + \sum_{a}V_{a} + \sum_{s}V_{s} + \sum_{o}V_{o} + \sum_{t}V_{t} + \sum_{v}V_{v} + \sum_{c}V_{c}\]

where the seven terms respectively describe bond stretching (b), angle bending (a), stretch-bending (s), out-of-plane bending (o), torsion (t), Van der Waals (v) and electrostatic (c) interactions. Their functional forms are given below.

For a bond b between atoms i and j the stretching potential is:

\[V_b = 143.9325 \frac{k_{ij}}{2} \Delta r_{ij}^2(1 + c_s\Delta r_{ij} + 7/12c_s^2 \Delta r_{ij}^2)\]

where \(k_{ij}\) is the force constant, \(\Delta r_{ij}\) is the difference between actual and reference bond lengths, and \(c_s\) is so called “cubic-stretch” constant for which the value is -2 Å\(^{-1}\).

The bending potential of a bond angle a made by the bonds between atoms i, j and atoms j, k is given by:

\[V_a = 0.043844\frac{k_{ijk}}{2} \Delta \vartheta_{ijk}^{2}(1+c_b \Delta \vartheta_{ijk})\]

where \(k_{ijk}\) is the force constant, \(\Delta \vartheta_{ijk}\) is the difference between actual and reference bond angles, and \(c_b\) is the so called “cubic-bend” constant for which the value is \(-0.4 rad^{-1}\).

The coupling between the stretching potential of two bonds forming a bond angle and bending that angle is described by:

\[V_s = 2.5121(k_{ijk} \Delta r_{ij} + k_{kji} \Delta r_{kj})\Delta \vartheta_{ijk}\]

where \(k_{ijk}\) and \(k_{kji}\) are force constants which couple stretches of i-j and k-j to the i-j-k bend respectively. \(\Delta r\) and \(\Delta \vartheta\) are defined above.

For a trigonal center j, the potential of displacement for an atom l bonded to atom j out of plane i-j-k is:

\[V_o = 0.043844\frac{k_{ijkl}}{2} \chi_{ijkl}^2\]

where \(k_{ijkl}\) is the force constant and \(\chi_{ijkl}\) is an angle formed by the bond j-l and the plane i-j-k.

For every four bonded atoms i-j-k-l the torsion interaction is described by the term:

\[V_t = 0.5(V_1(1 + \cos\Phi) + V_2(1 - \cos2\Phi) +V_3(1+ \cos3\Phi)\]

where \(V_1\), \(V_2\) and \(V_3\) are constants depending on atoms i, j, k, l, and \(\Phi\) is the dihedral angle formed by bonds i-j and k-l.

For a pair of atoms i and j separated by three or more bonds, where the distance between them is \(r_{ij}\), MMFF adopts the following Van der Waals potential:

\[V_v = \epsilon_{ij} \left(\frac{1.07R_{ij}}{r_{ij}+0.07R_{ij}} \right) ^7 \left( \frac{1.12R_{ij}^7}{r_{ij}^7+0.12R_{ij}} -2 \right)\]

where \(R_{ij}\) and \(\epsilon_{ij}\) are defined as follows:

\[R_{i} = A_i\alpha_i^{0.25}\]

\[R_{ij} = 0.5(R_i + R_j)(1+B(1-\exp(-12\gamma_{ij}^2)))\]

\[\gamma_{ij} = (R_i - R_j)/(R_i + R_j)\]

\[\epsilon_{ij} = \frac{181.16G_iG_j \alpha_i \alpha_j}{(\alpha_i/N_i)^{0.5} + (\alpha_j/N_j)^{0.5}} \frac{1}{R_{ij}^6}\]

where \(\alpha_i\) is atomic polarizability of atom i, B is 0.2 or 0.0 if one of the atoms is polar hydrogen, \(N_i\) and \(N_j\) are the Slater-Kirkwood effective numbers of valence electrons, \(G_i\), \(G_j\) and \(A_i\) are scale factors.

The electrostatic interaction between two charged atoms i and j separated by at least three bonds is calculated from the standard Coulombic expression:

\[V_c = f\frac{332.0716q_iq_j}{D(r_{ij} + \delta)}\]

where D is the dielectric constant for which the default value is 1, \(q_i\) and \(q_j\) are the MMFF partial charges on atoms i and j, \(r_{ij}\) is the interatomic distance, \(\delta\) is the “electrostatic buffering” constant of 0.05 Å. Scaling factor f is 0.75 for 1,4 interactions, and 1.0 otherwise.

Protein-ligand interaction can be described by the Amber force field instead of MMFF potential:

\[V_{Amber} = \sum_{i,j} \left\{ \epsilon_{ij}\left[\left(\frac{R_{ij}}{r_{ij}}\right)^{12}
-2\left(\frac{R_{ij}}{r_{ij}}\right)^6\right] +
\frac{q_iq_j}{4\pi\epsilon_0r_{ij}}\right\}\]

where the summation is over protein-ligand atom pairs. \(r_{ij}\) is the interatomic distance, \(R_{ij}\) is the Vdw distance for a pair of atoms, \(q_i\) and \(q_j\) are the Amber partial charges on atoms i and j, and \(\epsilon_0\) is the vacuum permittivity. Vdw parameters \(R_{ij}\) and \(\epsilon_{ij}\) and partial charges are taken from ([Wang-2000]).

Ligand entropy is evaluated as a sum of configurational entropy (\(S_c\)) and solvation entropy (\(\Delta S_s\)):

\[S = S_c + \Delta S_s\]

Configurational entropy is calculated as:

\[S_c = kN \Bigg[ 1 + \ln \Bigg(\frac{q}{N} \Bigg) + \frac{T}{q} \frac{\partial q}{\partial T} \Bigg]\]

where q is the conformation dependent partition function:

\[q = q_t \sum_{i=1}^{n_c} e^{-\frac{\epsilon_i}{kT}} q_{iv} q_{ir}\]

Here \(q_t\), \(q_ir\) and \(q_iv\) are the translational, rotational and vibrational partition functions respectively, \(n_c\) is the number of unique conformations in the ensemble. All 3 partition functions are calculated from the classical statistical mechanics expressions which could be found in [McQuarrie-1976]. Vibrational frequencies for each conformation, needed for evaluation of \(q_{ir}\) are derived from diagonalization of a Hessian matrix obtained from Quasi-Newton optimization when convergence is achieved. Eigenvalues \(\lambda_i\) of the mass-weighted Hessian:

\[\mathbf{H}^{m} = \mathbf{M^{-1/2}HM^{-1/2}}\]

are converted into wavenumbers \(\tilde{\nu}_i\) according to:

\[\tilde{\nu}_i = \frac{1}{2\pi c} \sqrt{\lambda_i}\]

Solvation entropy is split into electrostatic and hydrophobic parts:

\[\Delta S_s = \Delta S_{s,elec} + \Delta S_{s,hyd}\]

The electrostatic part of solvation entropy is divided in to the bulk component and tight electrostatic polar solute - water interactions (hydrogen bonds). The bulk contribution is estimated from the temperature dependence of the solvent dielectric constant as:

\[\Delta S_{s,elec\_bulk} = -\bigg(\frac{\partial \Delta G_s}{\partial \epsilon_{solv}}\bigg) \bigg(\frac{\partial \epsilon_{solv}}{\partial T}\bigg)\]

The second term of the electrostatic solvation entropy is estimated as a constant of 28 J/(mol K).

The hydrophobic term, \(\Delta S_{s,hyd}\), is evaluated as:

\[\Delta S_{s,hyd} = -\bigg(\frac{\partial \Delta G_{s,hyd}}{\partial T}\bigg)\]

where \(\Delta G_{s,hyd}\) consists of 3 components:

\[\Delta G_{s,hyd} = \Delta G_{cav} + \Delta G_{VdW} + \Delta G_{Ind}\]

describing the free energy of cavitation, solute-solvent van der Waals and inductive terms respectively. The cavity formation term is calculated from Scaled Particle Theory [Pierotti-1976]. Analytical expressions for \(\Delta G_{VdW}\) and \(\Delta G_{Ind}\) terms are also taken from the 1976 Pierotti review.

Configurational entropy of a protein bound ligand is calculated totally as vibrational entropy for 3N modes, assuming that 3 rotational and 3 translational degrees of freedom of a solution ligand are transformed into low-vibrational degrees of freedom for the bound ligand.

Solvation entropy for a ligand in the active site is assumed to be a sum of its fractional value in solution determined by the percentage of the ligand surface exposed to the solvent, \(f\Delta S_s\), and a partial desolvation entropy of the protein active site, \(\Delta S_{des}\)

\[S_{protein} = S_v + f\Delta S_s + \Delta S_{des}\]

where f is the fraction of ligand surface exposed to the solvent. It is important to notice that \(S_{protein}\) is not an experimentally measurable value, and that only the difference between \(S_{protein}\) and \(S = S_c + \Delta S_s\) might be compared with experimental binding entropy.