Tables for
Volume B
Reciprocal space
Edited by U. Shmueli

International Tables for Crystallography (2006). Vol. B, ch. 3.3, pp. 377-380   | 1 | 2 |

Section Modelling methods

R. Diamonda*

aMRC Laboratory of Molecular Biology, Hills Road, Cambridge CB2 2QH, England
Correspondence e-mail: Modelling methods

| top | pdf |

Fundamental to the design of any software for molecular modelling are the choices of modelling criteria, and of parameterization. Criteria which may be adopted might include the fitting of electron density, the minimization of an energy estimate or the matching of complementary surfaces between a pair of molecules. Parameterizations which may be adopted include the use of Cartesian coordinates of atoms as independent variables, or of internal coordinates, such as dihedral angles, as independent variables with atomic positions being dependent on these. Systems designed to suit energetic criteria usually use Cartesian coordinates since all aspects of the structure, including bond lengths, must be treated as variables and be allowed to contribute to the energy estimate. Systems designed to fit a model to observed electron density, however, may adequately meet the stereochemical requirements of modelling on either parameterization, and examples of both types appear below.

Inputs to modelling systems vary widely. Systems intended for use mainly with proteins or other polymeric structures usually work with a library of monomers which the software may develop into a polymer. Systems intended for smaller molecules usually develop the molecular structure atom-by-atom rather than a residue at a time, and systems of this kind require a very general form of input. They may accept a list of atom types and coordinates if measurement and display of a known molecule is the objective, or they may accept `sketch-pad' input in the form of a hand-drawn two-dimensional sketch of the type conventional in chemistry, if the objective is the design of a molecule. Sketch-pad input is a feature of some systems with quantum-mechanical capabilities. Methods based on conformational variables

| top | pdf |

Suppose that t represents a vector from the current position of an atom in the model to a target position then (see Section[link], to first order, the observational equations are [t_{IA} = D_{IpA} \theta_{p} + v_{IA}] in which [\boldtheta] represents changes to conformational variables which may include dihedral angles, bond angles, bond lengths, and parameters determining overall position and orientation of the molecule as a whole. If every such parameter is included the model acquires 3n degrees of freedom for n atoms, in which case the methods of the next section[link] are more appropriate, but if bond lengths and some or all bond angles are being treated as constants then the above equation becomes the basis of the treatment. [D_{IPA} = {\partial r_{IA}\over \partial \theta_{P}} = \varepsilon_{Ijk} n_{jP} (r_{kA} - r_{kP})] in which [{\bf n}_{P}] is a unit vector defining the axis of rotation for an angular variable [\theta_{P}], [{\bf r}_{A}] and [{\bf r}_{P}] are position vectors of the atom A and the site of the parameter P, and [{\bf v}_{A}] represents a residual vector.

[\sigma = v_{ia} v_{ia}] is minimized by [{\boldtheta} = {\bi M}^{-1} {\bf V}] in which [\eqalign{ M_{PQ} &= D_{iPa}D_{iQa}\cr V_{P} &= D_{iPa}t_{ia}.}] More generally, if σ represents any scalar quantity which is to be minimized, e.g. an energy, then [\eqalign{ V_{P} &= -{1\over 2} {\partial \sigma\over \partial \theta_{P}}\cr M_{PQ} &= {1\over 2} {\partial^{2} \sigma\over \partial \theta_{P} \partial \theta_{Q}}.}]

It is beyond the scope of this chapter to review the methods available for evaluating [\boldtheta] from these equations. Difficulties may arise from two sources:

  • (i) Inversion of M may be difficult if M is large or ill conditioned and impossible if M is singular.

  • (ii) Successful evaluation of [{\bi M}^{-1}{\bf V}] will not minimize σ in one step if t is not linearly dependent on [\boldtheta] or, equivalently, [\partial^{2}\sigma / \partial \theta_{P}\partial \theta_{Q}] is not constant, and substantial changes [\boldtheta] are involved. Iteration is then necessary.

Difficulties of the first kind may be overcome by gradient methods, for example the conjugate gradient method without searches if M is available or with searches if it is not available, or they may be overcome by methods based on eigenvalue decompositions. If non-linearity is serious less dependence should be placed on M and gradient methods using searches are more valuable. In this connection Diamond (1966[link]) introduced a sliding filter technique which produced rapid convergence in extreme conditions of non-linearity. These topics have been reviewed elsewhere (Diamond, 1981[link], 1984b[link]) and are the subject of many textbooks (Walsh, 1975[link]; Gill et al., 1981[link]; Luenberger, 1984[link]).

Warme et al. (1972[link]) have developed a similar system using dihedral angles as variables and a variety of alternative optimization algorithms.

The modelling of flexible rings or lengths of chain with two or more fixed parts is sometimes held to be a difficulty in methods using conformational variables, although a simple two-stage solution does exist. The principle involved is the sectioning of the space of the variables into two orthogonal subspaces of which the first is used to satisfy the constraints and the second is used to perform the optimization subject to those constraints.

The algebra of the method may be outlined as follows, and is given in more detail by Diamond (1971[link], 1980a[link],b[link]). Parametric shifts [{\boldtheta}_{1}] which satisfy the constraints are solutions of [{\bf V}_{1} = {\bi M}_{1} {\boldtheta}_{1}] in which [{\bf V}_{1}] and [{\bi M}_{1}] depend only on the target vectors, [{\bf t}_{1}], of the atoms with constrained positions and on the corresponding derivatives. We then find a partitioned orthogonal matrix [({\bi A}_{1}{\bi B}_{1})] satisfying [({\bi A}_{1}{\bi B}_{1}) = ({\bi E}_{\lambda}{\bi E}_{0}) \pmatrix{{\bi R}_{\lambda} &{\bi 0}\cr {\bi 0} &{\bi R}_{0}\cr}] in which [{\bi E}_{\lambda}] are the eigenvectors of [{\bi M}_{1}] having positive eigenvalues, [{\bi E}_{0}] are those having zero eigenvalues, and [{\bi R}_{\lambda}] and [{\bi R}_{0}] are arbitrary orthogonal matrices. Then [\openup-3pt\displaylines{\eqalign{\pmatrix{{\bi A}_{1}^{T}\cr {\bi B}_{1}^{T}\cr} {\bf V}_{1} &= \pmatrix{{\bi A}_{1}^{T}\cr {\bi B}_{1}^{T}\cr} {\bi M}_{1}({\bi A}_{1}{\bi B}_{1}) \pmatrix{{\bi A}_{1}^{T}\cr {\bi B}_{1}^{T}\cr} {\boldtheta}_{1}\cr \noalign{\vskip3pt} &= \pmatrix{{\bi A}_{1}^{T}{\bi M}_{1}{\bi A}_{1} &{\bi 0}\cr {\bi 0} &{\bi 0}\cr} \pmatrix{{\bi A}_{1}^{T}\cr {\bi B}_{1}^{T}\cr} {\boldtheta}_{1}\cr} \cr  {\bi A}_{1}^{T}{\boldtheta}_{1} = ({\bi A}_{1}^{T}{\bi M}_{1}{\bi A}_{1})^{-1} {\bi A}_{1}^{T}{\bf V}_{1}}] in which the matrix to be inverted is positive definite. [{\bi A}_{1}], however, is rectangular so that multiplying on the left by [{\bi A}_{1}] does not necessarily serve to determine [{\boldtheta}_{1}], but we may write [{\boldtheta}_{1} = ({\bi A}_{1}{\bi B}_{1}) \pmatrix{{\boldvarphi}_{1}\cr {\boldpsi}_{1}\cr}] giving [\displaylines{\pmatrix{{\bi A}_{1}^{T}\cr {\bi B}_{1}^{T}\cr} {\bf V}_{1} = \pmatrix{{\bi A}_{1}^{T}{\bi M}_{1}{\bi A}_{1} &{\bi 0}\cr {\bi 0} &{\bi 0}\cr} \pmatrix{{\boldvarphi}_{1}\cr {\boldpsi}_{1}\cr}\cr \noalign{\vskip3pt} {\boldvarphi}_{1} = ({\bi A}_{1}^{T}{\bi M}_{1}{\bi A}_{1})^{-1} {\bi A}_{1}^{T}{\bf V}_{1}}] and [{\boldpsi}_{1}] is indeterminate and free to adopt any value. We therefore adopt [{\boldtheta}_{1} = {\bi A}_{1}{\boldvarphi}_{1} = {\bi A}_{1}({\bi A}_{1}^{T}{\bi M}_{1}{\bi A}_{1})^{-1} {\bi A}_{1}^{T}{\bf V},] which is the smallest vector of parametric shifts which will satisfy the constraints, and allow [{\boldpsi}_{1}] to be determined by the remaining observational equations in which the target vectors, t, are now modified to [{\bf t}_{2}] according to [{\bf t}_{2} = {\bf t} - {\bi D}_{2}{\boldtheta}_{1},] [{\bi D}_{2}] and [{\bf t}_{2}] being the derivatives and effective target vectors for the unconstrained atoms. We then solve [{\bf V}_{2} = {\bi M}_{2}{\boldtheta}_{2}] in which [{\boldtheta}_{2}] is required to be of the form [{\boldtheta}_{2} = {\bi B}_{1}{\boldpsi}_{1}] giving [{\boldtheta}_{2} = {\bi B}_{1} ({\bi B}_{1}^{T}{\bi M}_{2}{\bi B}_{1})^{-1} {\bi B}_{1}^{T}{\bf V}_{2}] and apply the total shifts [{\boldTheta} = {\boldtheta}_{1} + {\boldtheta}_{2}] to obtain a structure which is optimized within the restrictions imposed by the constraints.

It may happen that [{\bi B}_{1}^{T} {\bi M}_{2} {\bi B}_{1}] is itself singular because there are insufficient data in the vector [{\bf t}_{2}] to control the structure and the parametric shifts contained in [{\boldtheta}_{2}] fully. In this event the same process may be applied again, basing the solution for [{\boldtheta}_{2}] on [\openup3pt\pmatrix{{\bi A}_{2}^{T}\cr {\bi B}_{2}^{T}\cr} {\bi B}_{1}^{T}{\bi M}_{2}{\bi B}_{1} ({\bi A}_{2}{\bi B}_{2})] so that the vectors in [{\bi B}_{2}] represent the degrees of freedom which remain uncommitted. This method of application of constraints by subspace sectioning may be nested to any depth and is completely general.

A valid matrix [{\bi A}_{1}] may be found from [{\bi M}_{1}] by using the fact that the columns of [{\bi M}_{1}] are all linear combinations of the columns of [{\bi E}_{\lambda}] and are void of any contribution from [{\bi E}_{0}]. It follows that [{\bi A}_{1}] may be found by using the columns of [{\bi M}_{1}] as priming vectors in the Gram–Schmidt process [Section[link] (i[link])] until the normalizing step involves division by zero. [{\bi A}_{1}] is then complete if all the columns of [{\bi M}_{1}] have been tried. [({\bi A}_{1}{\bi B}_{1})] may then be completed by using arbitrary vectors as primers.

Manipulation of a ring of n atoms may be achieved by treating it as a chain of [(n + 2)] atoms [having [(n + 1)] bond lengths, n bond angles and [(n - 1)] dihedral angles] in which atom 1 is required to coincide with atom [(n + 1)] and atom 2 with [(n + 2)]. [{\bf t}_{1}] then contains two vectors, namely the lack-of-closure vectors at these points, and is typically zero. [{\bi A}_{1}] is then found to have five columns corresponding to the five degrees of freedom of two points of fixed separation; [\boldtheta_{1}] contains only zeros if the ring is initially closed, and contains ring-closure corrections if, through non-linearity or otherwise, the ring has opened. [{\bi B}_{1}] contains [(p - 5)] columns if the chain of [(n + 2)] points has p variable parameters. It follows, if bond lengths and bond angles are treated as constants, that the seven-membered ring is the smallest ring which is flexible, that the six-membered ring (if it can be closed with the given bond angles) has no flexibility (though it may have discrete alternatives) and that it may be impossible to close a five-membered ring. Therefore some variation of bond angles and/or bond lengths is essential for the modelling of flexible five- and six-membered rings. Treating the ring as a chain of [(n + 1)] atoms is less satisfactory as there is then no control over the bond angle at the point of ring closure.

A useful concept for the modelling of flexible five-membered rings with near-constant bond angles is the concept of the pseudorotation angle P, and amplitude [\theta_{m}], for which the jth dihedral angle is given by [\theta_{j} = \theta_{m} \cos \left(P + {4\pi j\over 5}\right).] This formulation has the property [\textstyle{\sum}_{j = 0}^{4} \theta_{j} = 0], which is not exactly true; nevertheless, [\theta_{j}] values measured from observed conformations comply with this formulation within a degree or so (Altona & Sundaralingam, 1972[link]).

Software specialized to the handling of condensed ring systems has been developed by van der Lieth et al. (1984[link]) (Section[link]) and by Cohen et al. (1981[link]) (Section[link]). Methods based on positional coordinates

| top | pdf |

Modelling methods in which atomic coordinates are the independent variables are mathematically simpler than those using angular variables especially if the function to be minimized is a quadratic function of interatomic distances or of distances between atoms and fixed points. The method of Dodson et al. (1976[link]) is representative of this class and it may be outlined as follows. If d is a column vector containing ideal values of the scalar distances from atoms to fixed target points or to other atoms, and if l is a column vector containing the prevailing values of these quantities obtained from the model, then [{\bf d} = {\bf l} + {\bi D}\bolddelta {\bf x} + \boldvarepsilon ] in which the column matrix [{\bolddelta}{\bf x}] contains alterations to the atomic coordinates, [\boldvarepsilon] contains residual discrepancies and D is a large sparse rectangular matrix containing values of [\partial l/\partial x], of which there are not more than six non-zero values on any row, consisting of direction cosines of the line of which l is the length. [\boldvarepsilon^{T}{\bi W} \boldvarepsilon] is then minimized by setting [{\bi D}^{T}{\bi W} ({\bf d} - {\bf l}) = {\bi D}^{T}{\bi WD}\bolddelta {\bf x},] which they solve by the method of conjugate gradients without searches. This places reliance on the linearity of the observational equations (Diamond, 1984b[link]). It also works entirely with the sparse matrix [{\bi W}^{1/2}{\bi D}], the dense matrix [{\bi D}^{T}{\bi WD}], and its inverse, being never calculated.

The method is extremely efficient in annealing a model structure for which an initial position for every atom is available, especially if the required shifts are within the quasi-linear region, but is less effective when large dihedral-angle changes are involved or when many atoms are to be placed purely by interpolation between a few others for which target positions are available. Interbond angles are controlled by assigning d values to second-nearest-neighbour distances; this is effective except for bond angles near 180° so that, in particular, planar groups require an out-of-plane dummy atom to be included which has no target position of its own but does have target values of distances between itself and atoms in the planar group. The method requires a d value to be supplied for every type of nearest- and next-nearest-neighbour distance in the structure, of which there are many, together with W values which are the inverse variances of the distances concerned as assessed by surveys of the corresponding distances in small-molecule structures, or from estimates of their accuracy, or from estimates of accuracy of the target positions.

Hermans & McQueen (1974[link]) published a similar method which differs in that it moves only one atom at a time, in the environment of its neighbours, these being considered fixed while the central atom is under consideration. This is inefficient in the sense that in any one cycle one atom moves only a small fraction (∼3%) of the distance it will ultimately be required to move, but individual atom cycles are so cheap and simple that many cycles can be afforded. The method was selected for inclusion in Frodo by Jones (1978[link]) (Section[link]) and is an integral part of the GRIP system (Tsernoglou et al., 1977[link]; Girling et al., 1980[link]) (Section[link]) for which it was designed. Approaches to the problem of multiple minima

| top | pdf |

Modelling methods which operate by minimizing an objective function of the coordinates (whether conformational or positional) suffer from the fact that any realistic objective function representing the potential energy of the structure is likely to have many minima in the space of the variables for any but the simplest problems. No general system has yet been devised that can ensure that the global minimum is always found in such cases, but we indicate here two approaches to this problem.

The first approach uses dynamics to escape from potential-energy minima. Molecular-mechanics simulations allow each atom to possess momentum as well as position and integrate the equations of motion, conserving the total energy. By progressively removing energy from the simulation by scaling down the momentum vectors some potential-energy minimum may be found. Conversely, a minimization of potential energy which has led to a minimum thought not to be the global minimum may be continued by introducing atomic momenta sufficient to overcome potential-energy barriers between minima, and subsequently attenuate the momenta again. In this way a number of minima may be found (Levitt & Warshel, 1975[link]). It is equivalent to initializing a potential-energy minimization from a number of different conformations but it has the property that the minima so found are separated by energy barriers for which an upper limit is known so that the possibility exists of exploring transition pathways.

A second approach (Purisima & Scheraga, 1986[link]) is relatively new. If the objective function to be minimized can be expressed in terms of interatomic distances, and if each atom is given coordinates in a space of [n - 1] dimensions for n atoms, then a starting structure may be postulated for which the interatomic distances all take their ideal values and the objective function is then necessarily at an absolute minimum. This multidimensional structure is then projected into a space of fewer dimensions, within which it is again optimized with respect to the objective function. The dimensionality of the model is thus progressively reduced until a three-dimensional model is attained at a low energy. This means that the minimum so attained in three dimensions is approached from beneath, having previously possessed a lower value in a higher-dimensional space. This, in itself, does not guarantee that the three-dimensional minimum-energy structure so found is at the global minimum, but it is not affected by energy barriers between minima in the same way, and it does appear to reach very low minima, and frequently the global one. Because it is formulated entirely in terms of interatomic distances it offers great promise for modelling molecules on the basis of data from nuclear magnetic resonance.


Altona, C. & Sundaralingam, M. (1972). Conformational analysis of the sugar ring in nucleosides and nucleotides. A new description using the concept of pseudorotation. J. Am. Chem. Soc. 94(23), 8205–8212.Google Scholar
Cohen, N. C., Colin, P. & Lemoine, G. (1981). Script: interactive molecular geometrical treatments on the basis of computer-drawn chemical formula. Tetrahedron, 37, 1711–1721.Google Scholar
Diamond, R. (1966). A mathematical model-building procedure for proteins. Acta Cryst. 21, 253–266.Google Scholar
Diamond, R. (1971). A real-space refinement procedure for proteins. Acta Cryst. A27, 436–452.Google Scholar
Diamond, R. (1980a). BILDER: a computer graphics program for biopolymers and its application to the interpretation of the structure of tobacco mosaic virus protein discs at 2.8 Å resolution. In Biomolecular structure, conformation, function and evolution, Vol. 1, edited by R. Srinivasan, pp. 567–588. Oxford: Pergamon Press.Google Scholar
Diamond, R. (1980b). Some problems in macromolecular map interpretation. In Computing in crystallography, edited by R. Diamond, S. Ramaseshan & K. Venkatesan, pp. 21.01–21.19. Bangalore: Indian Academy of Sciences for the International Union of Crystallography.Google Scholar
Diamond, R. (1981). A review of the principles and properties of the method of least squares. In Structural aspects of biomolecules, edited by R. Srinivasan & V. Pattabhi, pp. 81–122. Delhi: Macmillan India Ltd.Google Scholar
Diamond, R. (1984b). Least squares and related optimisation techniques. In Methods and applications in crystallographic computing, edited by S. R. Hall & T. Ashida, pp. 174–192. Oxford University Press.Google Scholar
Dodson, E. J., Isaacs, N. W. & Rollett, J. S. (1976). A method for fitting satisfactory models to sets of atomic positions in protein structure refinements. Acta Cryst. A32, 311–315.Google Scholar
Gill, P. E., Murray, W. & Wright, M. H. (1981). Practical optimization. Orlando, Florida: Academic Press.Google Scholar
Girling, R. L., Houston, T. E., Schmidt, W. C. Jr & Amma, E. L. (1980). Macromolecular structure refinement by restrained least-squares and interactive graphics as applied to sickling deer type III hemoglobin. Acta Cryst. A36, 43–50.Google Scholar
Hermans, J. & McQueen, J. E. (1974). Computer manipulation of (macro) molecules with the method of local change. Acta Cryst. A30, 730–739.Google Scholar
Jones, T. A. (1978). A graphics model building and refinement system for macromolecules. J. Appl. Cryst. 11, 268–272.Google Scholar
Levitt, M. & Warshel, A. (1975). Computer simulation of protein folding. Nature (London), 253, 694–698.Google Scholar
Lieth, C. W. van der, Carter, R. E., Dolata, D. P. & Liljefors, T. (1984). RINGS – a general program to build ring systems. J. Mol. Graphics, 2, 117–123.Google Scholar
Luenberger, D. G. (1984). Linear and nonlinear programming. Reading, Massachusetts: Addison Wesley.Google Scholar
Purisima, E. O. & Scheraga, H. A. (1986). An approach to the multiple-minima problem by relaxing dimensionality. Proc. Natl Acad. Sci. USA, 83, 2782–2786.Google Scholar
Tsernoglou, D., Petsko, G. A., McQueen, J. E. & Hermans, J. (1977). Molecular graphics: application to the structure determination of a snake venom neurotoxin. Science, 197, 1378–1381.Google Scholar
Walsh, G. R. (1975). Methods of optimization. London: John Wiley.Google Scholar
Warme, P. K., Go, N. & Scheraga, H. A. (1972). Refinement of X-ray data of proteins. 1. Adjustment of atomic coordinates to conform to a specified geometry. J. Comput. Phys. 9, 303–317.Google Scholar

to end of page
to top of page