Mathematical Challenges from Theoretical/Computational Chemistry

Introduction

Mathematical optimization is a branch of mathematics that seeks to answer the question "What is best?" for problems in which the quality of the answer can be expressed as a numerical value (see, e.g., Gill et al., 1983b; Fletcher, 1987; Ciarlet, 1989). This question might refer to the "best" approximation in some local sense (i.e., a local solution) or to the global solution over the entire feasible space (i.e., the global minimum) (see, e.g., Nemhauser et al., 1989; Floudas and Pardalos, 1991). A com mon problem arises when a complex physical system is described by a collection of particles, or combinations of states, in a multidimensional phase space. An energy or cost function is associated with each different configuration, and the challenge is to find sets of points that minimize (or maximize) the objective function. Such applications arise frequently in molecular modeling, rational drug design, quantum mechanical calculations, mathematical biology models, neural networks, combinatorial problems, financial investment planning, engineering, electronics, meteorology, and computational geometry. In applications that arise in computational chemistry (Scheraga, 1992; Schlick, 1992), the feasible space is often very high in dimensionality and complexi ty, so both local and global minima are of interest.

There are many optimization techniques available for the computational scientist. Nonetheless, implementation of the more sophisticated techniques requires considerable computing experience, algorithm familiarity, and intuition. While software vendors offer a variety of "black- box" codes, serious practitioners frequently discover that a good deal of understanding and modification is required for successful applicati ons. Such modifications involve tailoring the algorithm to features of the problem at hand--such as function separability--or exploiting available experimental information that might guide the optimization path--such as nuclear magnetic resonance (NMR) d istance restraints in molecular mechanics. Moreover, successful new optimization schemes may be not be known or available to nonspecialist mathematicians, let alone to scientists in allied fields.

Thus, the transfer of knowledge, its application to real problems, and its further developments will greatly benefit from increased interdisciplinary interactions. In particular, it may be useful to stimulate algorithmic developments in optimization to ward important scientific problems, such as arise in chemistry, that would involve synergistic efforts on both parts: the application- oriented scientist and the algori thm developer. Such collaborations are likely to be fruitful to both parties, since testing of new methods will be possible on real- life problems and might generate a n evolving body of solutions that take into account the available physical data. There is recent evidence (e.g., in special sessions of meetings of the Society for Industrial and Applied Mathematics) that mathematicians have discovered the challenges in "mathematical chemistry problems" and protein folding, but many frontiers lie ahead.

Problem Classification

The available optimization algorithms are classified according to the features of the target problem. The objective function to be minimized (or maximized) may be formulated in terms of integer variables (discrete optimization), integers in p ermutations (combinatorial optimization), continuous real numbers (continuous optimization), or both continuous real numbers and integers (mixed integer optimization). Examples from these four classes involve, respectively, order planning for organizations (the integers may denote, for example, the number of units of each item to be purchased monthly for a restaurant); the traveling salesman problem (the ordered list of N integers represents a cyclical itinerary for visiting N ci ties); molecular structure prediction (the real numbers may denote nuclear or electronic positions of the particles, or a set of internal variables describing the molecular system); and airline crew scheduling (the integers may identify particular flight routes and the real numbers may refer to the hours of shift for the flight crew). For computational chemistry, continuous optimization is the most important type of problem.

In addition to the nature of the control (or independent) variables, the objective function may be linear, quadratic, or nonlinear (the latter in varying extent). The problem may be formulated as unconstrained or constrained, with constraints involving equality or inequality conditions, which may be linear, quadratic, or nonlinear. Thus, for the above examples, constrained formulations may introduce upper and lower bounds for the Cartesian positions or specified values for certain internal variables t hat should remain fixed; the airline crew scheduling problem will incorporate into the optimization formulation the total number of scheduled flights, lower and upper bounds for the lengths of shifts, enforced limits on gaps between transatlantic flights, and so on. In addition to functional form and constraints, other important considerations involve the cost of evaluating the objective function and the availability (or lack) and associated cost of derivatives. In some cases, the derivatives may be dis continuous, and special techniques may be required. Derivative information can often be exploited significantly for the optimization algorithm, but the benefits must be balanced with the additional costs involved.

The Complexity of Computational Chemistry Problems

Optimization problems frequently arise in molecular and quantum mechanical calculations in chemistry. These problems are typical of optimization applications seeking favorable configurational states of a physical system. The large-scale nature of these problems together with the lack of convexity rules out exhaustive sampling in the feasible space except for very small systems. Therefore, clever optimization methods are a necessity, and their improvement translates into the ability to model larger phys ical systems and generate important structural predictions.

The expense of calculating the function and the associated derivatives also introduces difficulties that limit the type of algorithm that may be utilized. In many molecular mechanics applications, it may be tedious but possible to calculate the derivati ves; often, the additional computational cost involved in computing the gradient is only a small factor more (e.g., 4 to 5) than computing the function (and guaranteed by automatic differentiation, which also saves coding efforts; see Box 4.5).

Computational intensity often stems from the long- range interactions among the N particles in the system (e.g., Coulombic forces). In molecular mechanics, the di rect evaluation requires on the order of N operations, and even if a cutoff radius is introduced, computation of the nonbonded terms dominates computation time. Implementation of fast particle methods (Green gard, 1994) in molecular mechanics and dynamics calculations (Grubmüller et al., 1991; Board et al., 1992) is clearly important for reducing the severity of this problem and allowing more accurate representation of the long-range interactions; the ad vantage of such an approach has already been demonstrated in other scientific applications (Greengard, 1994), for example in the context of integral equations in engineering problems (Nabors et al., 1994).

The multiple-minimum problem is a severe hurdle in many large-scale optimization applications. The state of the art today is such that only for small and reasonable problems do suitable algorithms exist for finding all local minima for linear and nonli near functions. For larger problems, however, many trials are generally required to find local minima, and finding the global minimum cannot be ensured. These features have prompted research in conformational-search techniques independent of, or in comb ination with, minimization (Leach, 1991). To illustrate, consider a simple model for an alkane chain of m units (residues). From combinations or rough partitions in favorable structures of the individual building blocks, the number of possible starting p oints produces 3m starting configurations. For polypeptides and polynucleotides, the flexibility of the monomer (building block) configurations increases, producing a rough range of 10 to 25 reasonable starting points by coarse subdomain partition (e.g., combinations of typical side chain, main chain, backbone, or sugar dihedral angles). Exhaustive searches are clearly not feasible.

The buildup technique is a related configurational search strategy, used in studies of proteins (Pincus et al., 1982) and nucleic acids (Hingerty et al., 1989). Reasonable starting points are constructed by combining minima of conformatio nal building blocks. This rational strategy has performed rather well in practice, but there is no guarantee that all biologically important local minima, much less the global minimum, are revealed. One of the problems is the nonlocal nature of the inte ractions in the folded macromolecule. That is, segments far apart in the linear sequence will make close contact upon folding; thus, the collective minimum may not correspond to any minima of the constituent building blocks. Furthermore, the number of s tarting points is still exponential in the number of building blocks. This buildup technique might be an interesting mathematical area to explore further, perhaps through techniques of interval analysis (see below under global optimization methods).

Molecular dynamics, discussed on pages 54-58, can also be viewed as a technique for obtaining structural information (e.g., mean atomic fluctuations, dynamical pathways, isomerization rates) that is complementary to potential energy minimization. While in theory information on all thermally accessible states should be observable, the restriction of the integration time step to a very small value with respect to time scales of collective biomolecular motions limits the scope of molecular dynamics in pra ctice.

Local Optimization Methods

Local methods are defined by an iterative procedure that generates iterates {x, x,...,x,...} intended to converge to a loc al minimum x from a given x. Their performance is clearly sensitive to the choice of starting point in addition to search direction and algorithmic details. In the line- search subclass, a search vector p is computed at each step by a given strategy, and the objective funct ion f is minimized approximately along that direction so that "sufficient decrease" is obtained (see, e.g., Dennis and Schnabel, 1983; Luenberger, 1984). In trust-region approaches, a local quadratic model of the function is minimized at every step using current Hessian information, and an optimal step is chosen to lie within the "trust region" of the quadratic model (Dennis and Schnabel, 1983).

Local deterministic optimization methods have experienced extensive development in the last decade (e.g., Nocedal, 1991; Wright, 1991). Studies have produced a range of robust and reliable techniques tailored to problem size, smoothness, complexity, an d memory considerations. Many variants of Newton's method have been produced that extend applicability far beyond small or sparse problems. Nonderivative methods are generally not competitive, but significant developments have been made in nonlinear con jugate gradient (CG) methods (generally recommended for very large problems whose function is very expensive to evaluate) and Newton methods.

The classes and extensions of Newton's method, the prototype of second-derivative algorithms, include discrete Newton, quasi- Newton (QN) (also termed variable metric), and truncated Newton (TN) (e.g., Dennis and Schnabel, 19 83; Gill et al., 1983b). Historically, because of the O(n)memory requirements, where n is the number of variables in the objective function, and the O(n) computation associated with solving a linear system directly, Newton methods have been most widely used (1) for small problems, (2) for problems with special sparsity patterns, or (3) when near a solution, after a gradient method has been applied. Fo rtunately, advances in computing technology and algorithmic developments have made the Newton approach feasible for a wide range of problems. Indeed, effective strategies have been tailored to available storage and computation, exhibiting good performanc e in theory and practice, and this trend will undoubtedly intensify.

Two specific classes are emerging as the most powerful techniques for large- scale applications: limited- memory quasi-Newton (LMQN) and truncated Newton methods. LMQN methods attempt to retain the modest storage and computational requirements of CG methods while approaching the superlinear convergence properties of standard (i.e., full memory) QN methods (Gilbert and Lemaréchal, 1989; Liu and Nocedal, 1989; Nash and Nocedal, 1991; Zou et al., 1993). Similarly, TN algorithms attempt to retain the rapid quadratic convergence rate of classic Newton methods whil e making computational requirements feasible for large-scale functions (Dembo and Steihaug, 1983; Nash, 1985; Schlick and Overton, 1987). With advances in automatic differentiation (see Box 4.5), the appeal of these methods will undoubtedly increase even further (Dixon, 1991).

Both limited- memory QN and TN methods are promising for computational chemistry problems. Moreover, they can be adapted to both constrained and unconstrained formula tions and can exploit the special composition (distinct components) of the potential energy function to accelerate convergence (Derreumaux et al., 1994). This issue involves a natural separation of the objective function into components of differing comp lexity (e.g., local and nonlocal interactions). This special composition can be exploited to construct banded or other sparse preconditioners in the context of CG and TN. Such problem tailoring requires some familiarity with the algorithmic modules and a lso demands knowledge of the theoretical and practical strengths and weaknesses of the different minimization methods. With rapidly growing improvements in high- perfo rmance vector and massively parallel machines, application-tailored software may be even more important in combination with parallel architectures whose design is motivated by specific applications.

Global Optimization Methods

In their attempt to find a global rather than local minimum, global optimization methods tend to explore larger regions of function space (see, e.g., Dixon and Szegö, 1975; Floudas and Pardalos, 1991). The global minimum of a function can be sough t through two classes of approaches: deterministic and stochastic. Deterministic methods usually require the objective function to satisfy certain smoothness properties; they construct a sequence of points converging to lower and lower loca l minima. Ideally, they attempt to "tunnel" through local barriers. Local minimization methods are often required repeatedly in the framework; hence, developments in local methods are likely to have an important impact on global techniques as well. Comp utational effort tends to be very large, and a guarantee of success can be obtained only under specific assumptions.

Stochastic global methods, on the other hand, involve systematic manipulation of randomly selected points (Nemhauser et al., 1989; Rinnooy Kan and Timmer, 1989; Schnabel, 1989; Törn and Zilinskas, 1989; Byrd et al., 1990). Success can be guaranteed only in an asymptotic, stochastic sense, although in practice many applications are very promising.

In the early days of global optimization (mid-1970s), most efforts focused on stochastic or heuristic approaches (Dixon and Szegö, 1975). In chemical applications, simulated annealing (Metropolis et al., 1953; Kirkpatrick et al., 1983; Dekkers and Aarts, 1991) is an appealing method of this class and is effective for small to medium molecular systems. It is also very easy to implement and generally requires no derivative computations. Indeed, there has been a wide application of this method to c hemical systems.

More recent efforts have focused also on deterministic global optimization methods. Interesting examples include the tunneling method (Levy and Gomez, 1985; Levy and Montalvo, 1985) and several innovative deterministic approaches in chemical application s (Purisima and Scheraga, 1986; Piela et al., 1989; Scheraga, 1992; Shalloway, 1992). In particular, in the mathematical community, two recent powerful methods have been identified that might be useful to chemical applications. One exploits convex prope rties and is based on differences of convex functions (Pardalos and Rosen, 1987; Horst and Tuy, 1993); the other is based on interval analysis (Hansen, 1980, 1992; Neumaier, 1990; Schnepper and Stadtherr, 1993). The convexity approach has been successful for global quadratic problems of up to approximately 300 variables and 50 constraints (Pardalos and Rosen, 1987; Horst and Tuy, 1993). Interval analysis, a field little known even to mathematicians, was pioneered by Hansen, among others. It involves co mputation of strict bounds to bracket the global minimum of a function. The algorithms involve various branch and bound techniques that recursively split the configuration space, aiming at bracketing the minimum as tightly as possible. Other information , such as bounds on derivatives, may also be generated. This class of methods can be applied to the solution of nonlinear systems, as well as global constrained and unconstrained optimization. However, these methods require second-derivative information (Hessians for optimization problems) and, moreover, the inverse of a preconditioning matrix to produce realistic bounds. For these reasons, interval analysis has been applied only to relatively small problems thus far. However, future research may be pr omising with preconditioning techniques that are now well developed for local optimization.

Perspective

In sum, the optimization applications that arise in computational chemistry offer challenging and rewarding problems to mathematicians. There is a need for the development of both local and global methods (the latter stochastic as well as deterministic ) and for transferring the technology rapidly from one discipline to another. In particular, optimization schemes will be more effective when all available chemical information (e.g., function separability, availability of derivatives, additional experim ental data) is taken into account in design of the algorithm, as is possible by preconditioning in both limited-memory quasi-Newton and truncated-Newton algorithms. Multigrid approaches (Kuruvila et al., 1994) and functional transformations (e.g., Piela et al., 1989; Wu, 1994) appear promising to global optimization problems in computational chemistry, and further developments might be fruitful.

Areas of mathematics that may have an important impact on the field are interval analysis and automatic differentiation. While the field of deterministic global optimization is still in its infancy in terms of general large-scale applicability, it is a nticipated that the exploitation of vector and massively parallel computing environments for algorithm design will lead to significant progress in the coming years. Technological advances will clearly improve the range of global optimization strategies t hat can be considered, but greater efforts in parallel programming skills will be essential so that these high-performance platforms will have a true impact on these important scientific problems.

References

Bischof, C.H., A. Carle, G.F. Corliss, A. Griewank, and P. Hovland, 1992, ADIFOR: Generating derivative codes from FORTRAN programs, Scientific Programming 1:1-29.

Board, J.A., Jr., J.W. Causey, T.F. Leathrum, Jr., A. Windemuth, and K. Schulten, 1992, Accelerated molecular dynamics simulations with the parallel fast multiple algorithm, Chem. Phys. Lett. 198:89-94.

Byrd, R.H., C.L. Dert, A.H.G. Rinnooy Kan, and R.B. Schnabel, 1990, Concurrent stochastic methods for global optimization, Math. Program. 46:1-29.

Ciarlet, P.G., 1989, Introduction to Numerical Linear Algebra and Optimization, Cambridge University Press, Cambridge.

Dekkers, A., and E. Aarts, 1991, Global optimization and simulated annealing, Math. Program. 50:367-393.

Dembo, R.S., and T. Steihaug, 1983, Truncated- Newton algorithms for large- scal e unconstrained optimization, Math. Prog. 26:190-212.

Dennis, Jr., J.E., and R.B. Schnabel, 1983, Numerical Methods for Unconstrained Optimization and Nonlinear Equations, Prentice- Hall, Englewood Cliffs, N.J.

Derreumaux, P., G. Zhang, B. Brooks, and T. Schlick, 1994, A truncated Newton minimizer adapted for CHARMM and biomolecular applications, J. Comput. Chem. 15:532-552.

Dixon, L.C.W., 1991, On the impact of automatic differentiation on the relative performance of parallel truncated Newton and variable metric algorithms, SIAM J. Opt. 1:475-486.

Dixon, L.C.W., and G.P. Szegö, 1975, Towards Global Optimization, Elsevier, New York.

Fletcher, R., 1987, Practical Methods of Optimization, Second Edition, John Wiley & Sons, New York.

Floudas, C.A., and P.M. Pardalos, eds., 1991, Recent Advances in Global Optimization, Princeton Series in Computer Science, Princeton University Press, Princeton, N.J.

Gilbert, J.C., and C. Lemaréchal, 1989, Some numerical experiments with variable- storage quasi- Newton algorithms, Math. Prog. 45:407-435.

Gill, P.E., W. Murray, M.A. Saunders, and M.H. Wright, 1983a, Computing forward- difference intervals for numerical optimization, SIAM J. Sci. Stat. Comput. 4:3 10-321.

Gill, P.E., W. Murray, and M.H. Wright, 1983b, Practical Optimization, Academic Press, N.Y.

Greengard, L., 1994, Fast algorithms for classical physics, Science 265:909-914.

Griewank, A., 1988, On automatic differentiation, in Mathematical Programming 1988, Kluwer Academic Publishers, Norwell, Mass. pp. 83-107.

Griewank, A., and G.F. Corliss, eds., 1991, Automatic Differentiation of Algorithms: Theory, Implementation, and Application, Society for Industrial and Applied Mathematics, Philadelphia, Pa.

Grubmüller, H., H. Heller, A. Windemuth, and K. Schulten, 1991, Generalized Verlet algorithm for efficient molecular dynamics simulations with long- range interac tions, Mol. Simul. 6:121-142.

Hansen, E.R., 1980, Global optimization using interval analysis--The multidimensional case, Numer. Math. 34:247-270.

Hansen, E.R., 1992, Global Optimization Using Interval Analysis, Dekker, New York.

Hingerty, B.E., S. Figueroa, T.L. Hayden, and S. Broyde, 1989, Prediction of DNA structure from sequence: A buildup technique, Biopolymers 28:1195-1222.

Horst, R., and H. Tuy, 1993, Global Optimization: Deterministic Approaches, second edition, Springer, Berlin.

Kirkpatrick, S., C.D. Gelatt, Jr., and M.P. Vecchi, 1983, Optimization by simulated annealing, Science 220:671-680.

Kuruvila, G., S. Ta'asan, and M.D. Salas, 1994, Airfoil optimization by the one shot method, Lecture Notes on Optimum Design Methods for Aerodynamics, AGARD FDP/VKI Special Course, von Karman Institute for Fluid Dynamics, Rhode- Saint- Genese, Belgium, April, 25-29.

Leach, A.R., 1991, A survey of methods for searching the conformational space of small and medium- sized molecules, in Reviews in Computational Chemistry, Vol. II, K.B. Lipkowitz and D.B. Boyd, eds., VCH Publishers, New York.

Levy, A.V., and S. Gomez, 1985, The tunneling method applied to global optimization, in Numerical Optimization 1984, P.T. Boggs, R.H. Byrd, and R.B. Schnabel, eds., SIAM, Philadelphia, pp. 213-244.

Levy, A.V., and A. Montalvo, 1985, The tunneling algorithm for the global minimization of functions, SIAM J. Sci. Stat. Comput. 6:15-29.

Liu, D.C., and J. Nocedal, 1989, On the limited memory BFGS method for large scale optimization, Math. Program. 45:503-528.

Luenberger, D.G., 1984, Linear and Nonlinear Programming, Second Edition, Addison- Wesley, Reading, Mass.

Metropolis, N., A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller, and E. Teller, 1953, Equation of state calculations by fast computing machines, J. Chem. Phys. 21:1087-1092.

Nabors, K., F.T. Korsmeyer, F.T. Leighton, and J. White, 1994, Multipole accelerated preconditioned iterative methods for three- dimensional potential integral equatio ns of the first kind, SIAM J. Sci. Comput. 15:713-735.

Nash, S.G., 1985, Solving nonlinear programming problems using truncated- Newton techniques in SIAM Numerical Optimization 1984, P.T. Boggs, R.H. Byrd, and R.B. Schnabel, eds., Philadelphia, pp. 119-136.

Nash, S.G., and J. Nocedal, 1991, A numerical study of the limited memory BFGS method and the truncated- Newton method for large- scale optimization, SIAM J. Opt. 1:358-372.

Nemhauser, G.L., A.H.G. Rinnooy Kan, and M.J. Todd, eds., 1989, Handbook in Operations Research Management Science, Vol. 1, Elsevier Science Publishers (North- Holland), Amsterdam.

Neumaier, A., 1990, Interval Methods for Systems of Equations, Cambridge University Press, Cambridge.

Nocedal, J., 1991, Theory of algorithms for unconstrained optimization, Acta Numerica 1:199-242.

Pardalos, P.M., and J.B. Rosen, 1987, Constrained Global Optimization: Algorithms and Applications, Lecture Notes in Computer Science 268, Springer, Berlin.

Piela, L., J. Kostrowicki, and H.A. Scheraga, 1989, The multiple- minima problem in conformational analysis of molecules. Deformation of the potential energy hypersurf ace by the diffusion equation method, J. Phys. Chem. 93:3339-3346.

Pincus, M., R. Klausner, and H.A. Scheraga, 1982, Calculation of the three- dimensional structure of the membrane bound portion of melittin from its amino acids, Pr oc. Natl. Acad. Sci. USA 79:5107-5110.

Purisima, E.O., and H.A. Scheraga, 1986, An approach to the multiple- minima problem by relaxing dimensionality, Proc. Natl. Acad. Sci. USA 83:2782-2786.

Rall, L.B., 1981, Automatic Differentiation--Techniques and Applications, Lecture Notes in Computer Science 120, Springer- Verlag, Berlin.

Rinnooy Kan, A.H.G., and G.T. Timmer, 1989, Global optimization, in Handbooks in Operations Research and Management Science, Vol. 1, G.L. Nemhauser, A.H.G. Rinnooy Kan, and M.J. Todd, eds., Elsevier Science Publishers (North- Holland), Amsterdam.

Scheraga, H.A., 1992, Predicting three- dimensional structures of oligopeptides, in Reviews in Computational Chemistry, Vol. III, K.B. Lipkowitz and D.B. Boyd, eds., VCH Publishers, New York, pp. 73-142.

Schlick, T., 1992, Optimization methods in computational chemistry, in Reviews in Computational Chemistry, Vol. III, K.B. Lipkowitz and D.B. Boyd, eds., VCH Publishers, New York, pp. 1-71.

Schlick, T., and A. Fogelson, 1992, TNPACK--A truncated Newton minimization package for large-scale problems: I. Algorithm and usage, and II. Implementation examples, ACM Trans. Math. Softw. 14:46-111.

Schlick, T., and M. Overton, 1987, A powerful truncated Newton method for potential energy minimization, J. Comput. Chem. 8:1025-1039.

Schnabel, R.B., 1989, Sequential and parallel methods for unconstrained optimization, in Mathematical Programming, M. Iri and K. Tanabe, eds., Kluwer Academic Publishers, Norwell, Mass., pp. 227-261.

Schnepper, C.A., and M.A. Stadtherr, 1993, Application of a parallel interval Newton/generalized bisection algorithm to equation- based chemical process flowsheeting, Interval Computations 4:40-64.

Shalloway, D., 1992, Application of the renormalization group to deterministic global minimization of molecular conformation energy functions, J. Global Opt. 2:281-311.

Törn, A., and A. Zilinskas, 1989, Global Optimization, Lecture Notes in Computer Science 350, Springer-Verlag, Berlin.

van Laarhoven, P.J.M., and E.H.L. Aarts, 1987, Simulated Annealing: Theory and Applications, D. Reidel, Dordrecht.

Wright, M.H., 1991, Interior methods for constrained optimization, Acta Numerica 1:341-407.

Wu, Z., 1994, The Effective Energy Transformation Scheme as a Special Continuation Approach to Global Optimization with Application to Molecular Conformation, Argonne National Laboratory report MCS-P442-0694, July.

Zou, X., I.M. Navon, F.X. Le Dimet, A. Nouailler, and T. Schlick, 1993, A comparison of efficient large- scale minimization algorithms for optimal control applications in meteorology, SIAM J. Opt. 3:582-608.

NAS Home Page | NAP Home Page | Reading Room | Report Home Page