Function

MODELLER minimizes the objective function $F$ with respect to Cartesian coordinates of $\sim 10,000$ atoms (3D points) that form a system (one or more molecules):

$$F = F({\bf R}) = F_{symm} + \sum_i c_i ({\bf f}_i, {\bf p}_i)$$ (5.5)

where $F_{symm}$ is an optional symmetry term defined in Eq. 5.72, $\bf R$ are Cartesian coordinates of all atoms, $c$ is a restraint $i$, ${\bf f}$ is a geometric feature of a molecule, and ${\bf p}$ are parameters. For a 10,000 atom system there can be on the order of 200,000 restraints. The form of $c$ is simple; it includes a quadratic function, cosine, a weighted sum of a few Gaussian functions, Coulomb law, Lennard-Jones potential, cubic splines, and some other simple functions. The geometric features presently include a distance, an angle, a dihedral angle, a pair of dihedral angles between two, three, four atoms and eight atoms, respectively, the shortest distance in the set of distances (not documented further), solvent accessibility in $\mbox{\AA}^2$, and atom density expressed as the number of atoms around the central atom. A pair of dihedral angles can be used to restrain such strongly correlated features as the mainchain dihedral angles $\Phi$ and $\Psi$. Each of the restraints also depends on a few parameters ${\bf p}_i$ that generally vary from a restraint to a restraint. Some restraints can restrain pseudo-atoms such as a gravity center of several atoms.

MODELLER allows some atoms to be fixed during optimization; i.e., only selected atoms are allowed to be moved. Similarly, MODELLER also allows only a subset of all restraints to be actually used in the calculation of the objective function. Each subset is indicated by a list of indices specifying the selected atoms or restraints.

There are two kinds of restraints, static and dynamic, that both contribute to the objective function as indicated in Eq. 5.5:

$$F = F_{symm} + F_s + F_d \; .$$ (5.6)

The static restraints and their parameters are pre-defined; i.e., they are given before the call to the optimizer and are not changed during optimization. The dynamic restraints are re-generated repeatedly during optimization. Usually, the CPU time is spent evenly between the two kinds of restraints, although the dynamic restraints become more important as the size of the system increases. All dynamic restraints are always selected and they can restrain only pairs of atoms. In all other respects, the two kinds of restraints are the same.

The dynamic restraints are obtained from a dynamic pairs list (the non-bonded pairs list). Each dynamic pair corresponds to at least one restraint, which may or may not be violated. The dynamic pairs list includes only the pairs of atoms that satisfy the following three conditions: (1) One or both atoms in a pair are allowed to move. (2) The two atoms are not connected through one, two, or three chemical bonds. (3) The two atoms are closer than a preset cutoff distance (e.g., 4 Å). There are on the order of 5000 atom pairs in the dynamic pairs list when only soft-sphere overlap restraints are used. Currently, the restraint types on the dynamic atom pairs that can be selected include the soft-sphere overlap, Lennard-Jones, Coulombic interactions, and MODELLER non-bonded spline restraints. xx atom density?

The existence of the dynamic pairs list is justified by the fact that dynamic pairs are usually a small fraction of all possible atom-atom pairs ( $N \cdot (N-1)/2$, where $N$ is the number of atoms in a system). The use of the dynamic pairs list becomes especially beneficent as the size of the system increases.

The actual algorithm for creating the dynamic pairs list varies with the size of the system, whether or not all atoms are allowed to move, or whether or not the user wants to include the fixed environment in the calculation of non-bonded restraints involving the selected atoms. See Section 3.6.35 for more information.

The hash-function algorithm is used to determine whether or not two atoms are a dynamic atom pair. This algorithm is about 20 times slower than a lookup table but it requires much less memory and still spends a negligible fraction of the total CPU time. A hash-function table is prepared only once before the start of the optimization and any other operation involving an evaluation of the objective function (e.g., model.optimize(), model.energy(), and model.pick_hot_atoms()).

The dynamic pairs list is not necessarily re-generated each time the objective function is evaluated, although the contribution of the restraint to the objective function is calculated in each call to the objective function routine with the current values of the Cartesian coordinates. The dynamic pairs list is re-generated only when maximal atomic shifts accumulate to a value larger than a preset cutoff. This cutoff is chosen such that there cannot be a violation of a restraint without having its atom pair on the dynamic pairs list. The dynamic pairs list is recalculated in $\sim 20\%$ and $\sim 2\%$ of the objective function calls at the beginning and the end of optimization, respectively.

Each evaluation of the objective function or of its first derivatives with respect to the Cartesian coordinates involves the following steps:

1. Calculate non-fixed pseudo-atoms from the current atomic positions
   (routine objfunc:pseudo).

2. Update the dynamic pairs list, if necessary (routine objfunc:upddyn).

3. Calculate the violations of selected restraints and all other quantities that are shared between the calculations of the objective function and its derivatives (routine objfunc:getviol).

4. Sum the contributions of all violated restraints to the objective function and the derivatives (routine objfunc:getviol).