# Solving Hartree-Fock systems with global optimization methods

Carlile Lavor, Leo Liberti, Nelson Maculan, Marco Antonio Chaer Nascimento

► To cite this version:

Carlile Lavor, Leo Liberti, Nelson Maculan, Marco Antonio Chaer Nascimento. Solving Hartree-Fock systems with global optimization methods. EPL - Europhysics Letters, 2007, 77 (5), 50006p1-50006p5. 10.1209/0295-5075/77/50006. [hal-00152406]

HAL Id: hal-00152406

https://hal.science/hal-00152406v1

Submitted on 4 Oct 2022

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Solving Hartree-Fock systems with global optimization methods

C. LAVOR¹(a), L. LIBERTI²(b), N. MACULAN³(c) and M. A. C. NASCIMENTO⁴(d)

¹ Department of Applied Mathematics (IMECC-UNICAMP), State University of Campinas - CP 6065, 13081-970, Campinas-SP, Brazil

² LIX, École Polytechnique - F-91128 Palaiseau, France

³ COPPE, Universidade Federal do Rio de Janeiro, UFRJ - CP 68511, Rio de Janeiro - RJ, 21945-970, Brazil

⁴ Departamento de Físico-Química, Instituto de Química, Universidade Federal do Rio de Janeiro, UFRJ, Rio de Janeiro - RJ, 21949-970, Brazil

received 21 September 2006; accepted in final form 11 January 2007
published online 23 February 2007

PACS 02.60.Pn – Numerical optimization

PACS 31.15.Ar – Ab initio calculations

PACS 71.15.Dx – Computational methodology (Brillouin zone sampling, iterative diagonalization, pseudopotential construction)

**Abstract** – The Hartree-Fock equations describe atomic and molecular electronic wave functions, based on the minimization of a functional of the energy. This can be formulated as a constrained global optimization problem involving nonconvex polynomials exhibiting many local minima. The traditional method of solving the Hartree-Fock problem does not provide a guarantee of global optimality and is very sensitive to the initial starting point. In this paper we show how to use a deterministic global optimization method to solve Hartree-Fock systems. The validity of the proposed approach was established by successfully computing the ground-state of the He and Be atoms.

Copyright © EPLA, 2007

**Introduction.** – The quantum behaviour of atoms and molecules, in the absence of relativistic effects, is determined by the time-independent Schrödinger equation:

$$H\Psi_n = E\Psi_n, \quad (1)$$

where *H*, the Hamiltonian operator of the system, represents the total energy (kinetic + potential) of all the particles of the system.

Analytical solutions for this equation are only possible for very simple systems. Hence, for the majority of the problems of interest, one has to rely on some approximate model. In the Hartree-Fock model, the electrons in atoms and molecules move independently of each other, the motion of each one of the electrons being determined by the attractive electrostatic potential of the nuclei and by a repulsive average field due to all the other electrons of the system. In this model, the approximate solutions $\Phi_n$ of eq. (1) are anti-symmetrized products of one-electron wave functions {$\varphi_i$} (also called orbitals), which

are solutions of the Hartree-Fock (HF) equations for the system under study [1]. This model gives rise to a set of coupled integro-differential equations which can only be solved numerically. Alternatively, each orbital $\varphi_i$ can be expanded in a complete basis set {$\chi_s$}$_{s=1}^\infty$. In order to transform the HF equations into a less cumbersome algebraic problem, we only consider a non-orthogonal finite subset {$\chi_s|s \le b$} of the basis, and we use it to approximate the orbitals. The larger we choose $b$, the better the approximation is likely to become.

The optimization problem considered in this paper arises because we need to find a set of coefficients $c_{si}$, for $s = 1, \dots, b$ and $i = 1, \dots, n$, such that for all $i \le n$ the function

$$\bar{\varphi}_i = \sum_{s=1}^{b} c_{si} \chi_s \quad (2)$$

is a good approximation of the $i$-th spatial orbital $\varphi_i$. A further requirement on the approximating set {$\bar{\varphi}_i$} is that it must be an orthogonal set. While not a necessary condition, orthogonality is always imposed in the Hartree-Fock method because the resulting equations are much easier to solve in an orthogonal basis of atomic or molecular orbitals [2]. The method most usually applied

<footnote>(a)E-mail: clavor@ime.unicamp.br</footnote>

<footnote>(b)E-mail: liberti@lix.polytechnique.fr</footnote>

<footnote>(c)E-mail: maculan@cos.ufrj.br</footnote>

<footnote>(d)E-mail: chaer@iq.ufrj.br</footnote>
to the Hartree-Fock equations iteratively solves a set of linear equations to find the coefficients $c_{si}$. This method, however, has two main limitations: (a) there is no guarantee that the set of coefficients $c_{si}$ found by the method is a globally optimal set; (b) it depends on an initial solution being available (starting guess); (c) the occupation number of all orbitals must be provided (electronic configuration).

We define the Hartree-Fock Problem (HFP) as the problem of finding a set of coefficients $c_{si}$ such that the $\bar{\phi}_i$ are the best possible approximations of the spatial orbitals. The objective function (quality of the approximation) is given by the energy function $E$ associated with the approximating set $\{\bar{\phi}_i\}$, which is guaranteed to be an upper bound to the energy function associated with the spatial orbitals [3]. The set $\{\bar{\phi}_i\}$ is required to be orthonormal for the reasons already mentioned. Furthermore, once each member of the set $\{\bar{\phi}_i\}$ is expanded in a complete orthonormal basis $\{\chi_s\}$, a normalization condition must be imposed on the expansion coefficients in order to preserve the probabilistic interpretation of the wave function. Thus, we need to minimize the energy function $E$ subject to orthonormality and normalization constraints. The decision variables of this mathematical programming problem are the coefficients $c_{si}$. For the orthonormal basis sets $\{\bar{\phi}_i\}$ and $\{\chi_s\}$, the problem can be expressed as follows:

$$ \left\{ \begin{array}{ll} \min_{c} & E(c) \\ \text{subject to (s.t.)} & \langle \bar{\phi}_i | \bar{\phi}_j \rangle = \delta_{ij} \quad \forall i \le j \le n \\ & \sum_{r=1}^{b} c_{ri}^2 = 1 \quad \forall i \le n \\ & c^L \le c \le c^U \end{array} \right. \qquad (3) $$

where $\delta_{ij}$ is the Kronecker delta-function, which is equal to 1, if $i=j$ and 0, otherwise. Problem (3) is a nonconvex, multi-extremal, polynomially constrained, polynomial programming problem, and falls therefore in the realm of Global Optimization (GO). We solve it by applying the spatial Branch-and-Bound (sBB) technique, a well-known deterministic GO method which yields an $\varepsilon$-guaranteed global optimum (for a given $\varepsilon > 0$) and does not need an initial solution, i.e. a set of starting guess coefficients. For validation purposes, this method has been applied to two instances of the HFP, namely to the helium and berillium atoms, with considerable success as regarding the CPU performance. In spite of the fact that He and Be are closed-shell systems, the sBB technique can be easily extended to open-shell systems (ROHF) and also to more sophisticate multiconfigurational wave functions, as will be discussed in future publications.

The rest of this paper is organized as follows. Next section summarizes the sBB algorithm for GO. The third section presents the mathematical programming formulation of the HFP. The fourth section describes the method used to compute the lower bound to the objective function value. The fifth section discusses the computational experience on the He and Be atoms. Last section concludes the paper.

**Global Optimization.** – A mathematical programming problem is formulated as follows:

$$ \left\{ \begin{array}{ll} \min_{x} & f(x) \\ \text{s.t.} & g(x) \le 0, \end{array} \right. \qquad (4) $$

where $x \in \mathbb{R}^n$ are the decision variables and $f: \mathbb{R}^n \to \mathbb{R}$ is the objective function to be minimized subject to a set of constraints $g: \mathbb{R}^n \to \mathbb{R}^m$ which may also include variable ranges or integrality constraints on the variables. Global Optimization is concerned with the solution of problems (4) where $f, g$ are non-convex nonlinear forms. A problem where $f, g$ are nonlinear is known as a Non Linear Programming problem (NLP); if some integrality constraints are present on the variable bounds the problem is known as a Mixed-Integer NLP (MINLP).

GO methods may be deterministic or stochastic according as to whether they provide a guarantee of optimality or they employ random search techniques. The most commonly employed deterministic GO method is an extension, called spatial Branch-and-Bound (sBB) [4–6], of the well-known Branch-and-Bound algorithm for implicit binary enumeration [8] to continuous spaces. In sBB, branching occurs by partitioning the continuous variable range in two or more sub-ranges. The recursive application of branching gives rise to a search tree. A node is fathomed (i.e. no further branching occurs on the node) either because the global minimum relative to the node has been found, or because the global minimum relative to the node cannot be better than the overall best solution found so far (the incumbent). In order to test these two conditions at each node, we compute a lower bound and an upper bound to the objective function value of the node’s problem restriction. The first condition is verified if these bounds differ by less than a pre-specified $\varepsilon > 0$ tolerance, and the second if the lower bound for the node is higher than the incumbent.

Because Branch-and-Bound algorithms have exponential worst-case complexity, the best way to run them efficiently is to provide a tight lower bound. In the specific case of the HFP, the lower bound is provided by solving a particular linear relaxation of the problem strengthened by several valid inequalities based on the quadratic structure of the problem and on the orthogonality constraints.

**Mathematical programming formulation of the problem.** – The expression for the Hartree-Fock electronic energy $E$ of a $2n$-electron molecule with closed shells is given by [3] as

$$ E = 2 \sum_{i=1}^{n} H_{ii}^{\text{core}} + \sum_{i=1}^{n} \sum_{j=1}^{n} (2J_{ij} - K_{ij}) + V_{NN}, \qquad (5) $$

where $H_{ii}^{\text{core}}$ contains the one-electron integrals, $J_{ij}$ is the Coulomb integral, $K_{ij}$ is the exchange integral, $V_{NN}$ is the nuclear repulsion term. The spatial orbitals $\{\phi_i | i \le n\}$ are expanded as linear combinations of a finite set of one-electron basis functions as per eq. (2). However,
the standard self-consistent procedure to solve the HF equations uses a non-orthogonal basis set {$\chi_s$}. Thus, in order to compare our results with the ones obtained by the standard procedure, the sBB technique will be also formulated in a non-orthogonal basis set {$\chi_s$}. The HF equations [3] are a set of b equations in the variables $c_{si}$:

$$ \sum_{s=1}^{b} c_{si} (F_{rs} - \varepsilon_i S_{rs}) = 0, \quad (6) $$

where $F_{rs} = \langle \chi_r | \hat{F} | \chi_s \rangle$, $S_{rs} = \langle \chi_r | \chi_s \rangle$, $\varepsilon_i$ is the orbital energy, and $\hat{F}$ is the Fock operator [3]. This is a nonlinear system, since the $\hat{F}$ operator depends on the orbitals {$\varphi_i$}, which in turn depend on the variables $c_{si}$.

It is possible to obtain an expression for $F_{rs}$ in terms of the coefficients {$c_{si}$} and a set of suitable integrals over the basis functions {$\chi_s$} [3]:

$$ F_{rs} = H_{rs}^{\text{core}} + \sum_{t=1}^{b} \sum_{u=1}^{b} \sum_{i=1}^{n} c_{ti}^{*} c_{ui} [2(rs|tu) - (ru|ts)], \quad (7) $$

where $(rs|tu)$ and $(ru|ts)$ stand for the Coulomb and exchange integrals between pairs of electrons. While the Coulomb integrals represent the quantum-mechanical equivalent of the classical Coulomb interaction between two charged particles, the exchange integrals are purely quantum entities, resulting from the fact that the total wave function for any multi-electronic system must be anti-symmetric (Pauli principle). Also, the wave functions representing the behaviour of atoms and molecules can be real or complex. However, since any observable (dynamic variable) must be real, the expectation value of the corresponding operator, $O$, must be taken as $\langle O \rangle = \int \phi^* O \phi dv$, which is the reason why the complex product $c_{ti}^* c_{ui}$ appears in eq. (7).

We write $H(r,s) = H_{rs}^{\text{core}}$ and $X(r,s,t,u) = (rs|tu)$, where $r=1, \dots, b$, $s=1, \dots, b$, $t=1, \dots, b$, and $u=1, \dots, b$. It is possible to rewrite eq. (5) in the form

$$ E = \sum_{r=1}^{b} \sum_{s=1}^{b} \sum_{i=1}^{n} (c_{ri}^{*} c_{si} (F_{rs} + H(r, s))) + V_{NN}. \quad (8) $$

Finally, using (7) and (8), we get

$$ E = \frac{1}{2} \sum_{r=1}^{b} \sum_{s=1}^{b} \sum_{t=1}^{b} \sum_{u=1}^{b} \left( P(r, s) P(t, u) \right. \\ \left. \times \left( X(r, s, t, u) - \frac{1}{2} X(r, u, t, s) \right) \right) \\ + \sum_{r=1}^{b} \sum_{s=1}^{b} (P(r, s) H(r, s)) + V_{NN}, \quad (9) $$

where

$$ P(j, k) = 2 \sum_{i=1}^{n} c_{ji}^{*} c_{ki} \quad (j = 1, \dots, b \text{ and } k = 1, \dots, b). $$

Note that expression (9) is actually a function of the coefficients $c_{si}$, since the integrals $H(r, s)$ and $X(r, s, t, u)$, as well as the value of $V_{NN}$, can be calculated once the basis {$\chi_s$} and the molecular geometry are defined.

As has been mentioned above, the decision variables of the HFP are the coefficients $c_{si}$ used in the basis expansion of the spatial orbitals. To further simplify the notation, we shall write the numerical problem parameters as

$$ \begin{aligned} \alpha_{rs}^{tu} &= X(r, s, t, u) - \frac{1}{2} X(r, u, t, s), \\ \beta_{rs} &= H(r, s), \\ \gamma &= V_{NN}. \end{aligned} $$

After simple term rearrangement, the objective function of the problem becomes

$$ E(c) = 2 \sum_{r,s=1}^{b} \left( \sum_{i=1}^{n} c_{ri} c_{si} \right) \\ \times \left( \sum_{t,u=1}^{b} \alpha_{rs}^{tu} \left( \sum_{i=1}^{n} c_{ti} c_{ui} \right) + \beta_{rs} \right) + \gamma. \quad (10) $$

The orthogonality constraints are $\langle \bar{\varphi}_i | \bar{\varphi}_j \rangle = \delta_{ij}$ for all $i \le j \le n$. Substituting $\bar{\varphi}_i = \sum_{r \le b} c_{ri} \chi_r$ for all $i \le n$ and $\langle \chi_r | \chi_s \rangle = S_{rs}$ for all $r, s \le b$ we obtain

$$ \sum_{r,s \le b} S_{rs} c_{ri} c_{sj} = \delta_{ij} \quad \forall i \le j \le n, \quad (11) $$

where normally $S_{rr} = 1$ for all $r \le b$. The variable bounds

$$ c^L \le c \le c^U \quad (12) $$

depend on the instance. The HFP can be succintly summarized as follows:

$$ \min\{E(c) \mid \text{s.t. } (11)-(12)\}, \quad (13) $$

as the minimization of a quartic objective function subject to quadratic constraints and variable bounds.

**Computing the lower bound.** – Recall that the most important step in the efficient running of the sBB algorithm is the computation of the lower bound at each node of the sBB search tree. In this section we shall discuss a method which produces reasonably tight bounds for the HFP. At any given sBB node, some of the decision variables are restricted to lie in subranges given by the recursive branching rules. We consider the restriction of the original problem (13) to the node’s variable ranges. We mention in passing that a reasonably good upper bound can be easily provided for by finding any local minimum of (13) restricted to the node (in practice, this is done by using a constrained nonlinear local solver such as SNOPT [9] on the problem formulation). The lower bound is calculated by finding a linear relaxation of the problem restricted to the node and solving it with Linear Programming (LP) tools (for example the
simplex method [10] implemented in CPLEX [11]). The linear relaxation is built in three steps: (a) reformulation, (b) relaxation, (c) tightening.

Step (a) reformulates the problem to a standard form [5] consisting of a linear objective function subject to linear constraints and a set of quadratic constraints. By introducing new variable sets $y, w$ with the following constraints:

$$y_{rs}^{ij} = c_{ri}c_{sj}, \quad \forall r, s \le b, \quad i, j \le n, \quad (14)$$

$$w_{rstu}^{ij} = y_{rs}^{ii} y_{tu}^{jj}, \quad \forall r, s, t, u \le b, \quad i, j \le n, \quad (15)$$

we can replace all the quartic and quadratic terms in the problem, obtaining

$$\left\{
\begin{array}{ll}
\min_y & 2 \sum_{r,s \le b} \left( \sum_{t,u \le b} \alpha_{rs}^{tu} \sum_{i,j \le n} w_{rstu}^{ij} + \beta_{rs} \sum_{i \le n} y_{rs}^{ii} \right) + \gamma \\
\text{s.t.} & \sum_{r,s \le b} S_{rs} y_{rs}^{ij} = \delta_{ij}, \quad \forall i \le j \le n \\
& \sum_{r \le b} y_{rr}^{ii} = 1, \quad \forall i \le n \\
& (14)-(15) \\
& c^L \le c \le c^U \\
& y^L \le y \le y^U \\
& w^L \le w \le w^U,
\end{array}
\right.
\quad (16)$$

where the bounds $y^L, y^U$ on $y$ and $w^L, w^U$ on $w$ are obtained through simple interval arithmetics using the bounds on $c$ and the quadratic relations (14)-(15).

Step (b) constructs a relaxation of (16) by replacing (14)-(15) by their convex envelopes. For simple constraints of the form $z = xt$ where $x^L \le x \le x^U$ and $t^L \le t \le t^U$ (such as (14)-(15)) the convex envelope is given by [12,13]:

$$z \le \min\{x^U t + t^L x - x^U t^L, x^L t + t^U x - x^L t^U\} \quad (17)$$

$$z \ge \max\{x^L t + t^L x - x^L t^L, x^U t + t^U x - x^U t^U\}, \quad (18)$$

which defines an enveloping tetrahedron around the points $(x^L, t^L)$, $(x^L, t^U)$, $(x^U, t^L)$, $(x^U, t^U)$. Purely quadratic constraints of the form $z = x^2$ are relaxed by the secant and the tangents of the parabola at $(x^L, (x^L)^2)$, $(x^U, (x^U)^2)$. This relaxation yields a linear problem whose optimal objective function value $\bar{f}$ is a lower bound to the globally optimal objective function value $f^*$ of (16), which is the same as that of (13). Since (17)-(18) vary as a function of the variable bounds, $\bar{f}$ depends on $c^L$, $c^U$, which means that the lower bound depends on the current sBB search tree node, as desired. In practice, it turns out that $\bar{f}$ is not a very tight bound. This is mostly due to the fact that the envelope (17)-(18) is generally not very close to the original surface (14)-(15). We try to improve this situation by adding some valid inequalities in the next step.

Step (c) generates valid constraints in order to tighten the linear relaxation obtained in steps (a) and (b).

<table><caption>Table I: Computational results for the He and Be atoms.</caption><thead><tr><th>Atom</th><th>sBB</th><th>VNS</th><th>SobolOpt</th></tr></thead><tbody><tr><td>He</td><td><b>0.26 s</b></td><td>0.116 s</td><td>0.14 s</td></tr><tr><td>Be</td><td><b>10 s</b></td><td>0.3 s</td><td>14 s</td></tr></tbody></table>

Recalling that $S_{rr} = 1$ for all $r \le b$, we first remark that when $i = j$, constraints (11) are

$$\sum_{r=1}^{b} c_{ri}^{2} + 2 \sum_{r<s} S_{rs} c_{ri} c_{si} = 1, \quad \forall i \le n, \quad (19)$$

which can be written in terms of the $y$ variables as

$$\sum_{r=1}^{b} y_{rr}^{ii} + 2 \sum_{r<s} S_{rs} y_{rs}^{ii} = 1, \quad \forall i \le n. \quad (20)$$

Secondly, we multiply (20) by problem variables $y_{tu}^{jj}$ for all $t, u \le b, j \le n$, obtaining

$$\sum_{r=1}^{b} y_{rr}^{ii} y_{tu}^{jj} + 2 \sum_{r<s} S_{rs} y_{rs}^{ii} y_{tu}^{jj} = y_{tu}^{jj}, \quad \forall t, u \le b, \quad i, j \le n.$$

We can now replace the quadratic products in the $y$ variables using the $w$ variables as per (15), obtaining a linear relation between the $y$ and the $w$ variables:

$$\sum_{r=1}^{b} w_{rrtu}^{ij} + 2 \sum_{r<s} S_{rs} w_{rstu}^{ij} = y_{tu}^{jj}, \quad \forall t, u \le b, \quad i, j \le n. \quad (21)$$

In the linear relaxation of (16) the $y$ and $w$ variables are only related through inequality constraints of type (17)-(18). Constraints (21) tighten the relaxation because in general a feasible region defined by equations has a smaller measure than one defined by inequalities ([14], theorem 4.1). Constraints (21) are called reduced RLT constraints [14–16] and form a subclass of the Reformulation-Linearization Technique (RLT) constraints described in [17].

**Computational results.** – We performed some preliminary computational tests on the ground state of the helium (He) and berillium (Be) atoms. For the He atom we considered an uncontracted basis of two s-type Gaussian functions and for the Be atom we used the STO-3G basis set, as detailed in [18]. The globally optimal solutions are $E_{\text{He}}^* = -2.7471h$ and $E_{\text{Be}}^* = -14.3519h$, and both were confirmed by the sBB algorithm. These results are in perfect agreement with the ones obtained through the standard self-consistent procedure to the HF problem. We performed some computational comparisons in terms of user CPU time (reported in table I) by solving the problem instances with two other global optimization solvers, namely Variable Neighbourhood Search (VNS) [19,20] and SobolOpt [21] (which are stochastic and do not provide an optimality guarantee).
All computational results have been obtained within the *ooOPS* optimization software framework [6,22] running Linux on an Intel Pentium III 850 MHz with 384 MB RAM.

**Conclusion.** – The usual way to solve Hartree-Fock equations has two limitations: (a) it provides a solution which is not guaranteed to be the optimal one, (b) it depends heavily on an initial solution being provided and (c) the occupation number of all orbitals must be provided (electronic configuration). We formulate this problem as a nonconvex optimization problem, which we solve using a spatial Branch-and-Bound algorithm for Global Optimization. The crucial step, *i.e.* the determination of the lower bound at each search tree node, relies on a convex relaxation which is further tightened by a suitable class of valid constraints. This approach overcomes both limitations of the usual method. The computational results presented in this paper refer to rather small test cases, but are nonetheless very promising. Computational work on larger cases is ongoing.

***

The authors are thankful to FAPERJ, FAPESP and CNPq for their support. One of the authors (LL) is thankful to Prof. D. KROB for financial support.

## REFERENCES

[1] ROOTHAN C. C. J., *Rev. Mod. Phy.*, **23** (1951) 69.

[2] LÖWDIN P. O., *J. Chem. Phys.*, **18** (1950) 365.

[3] LEVINE I. N., *Quantum Chemistry*, second edition (Prentice-Hall, Upper Saddle River, New Jersey) 2000.

[4] ADJIMAN C. S., DALLWIG S., FLOUDAS C. A. and NEUMAIER A., *Comput. Chem. Eng.*, **22** (1998) 1137.

[5] SMITH E. M. B. and PANTELIDES C. C., *Comput. Chem. Eng.*, **23** (1999) 457.

[6] LIBERTI L., *Writing global optimization software*, in [7], pp. 211-262.

[7] LIBERTI L. and MACULAN N. (Editors), *Global Optimization: From Theory to Implementation* (Springer, Berlin) 2006.

[8] AHO A. V., HOPCROFT J. E. and ULLMAN J. D., *Data Structures and Algorithms* (Addison-Wesley, Reading, MA) 1983.

[9] GILL P. E., *User's Guide for SNOPT 5.3*, Systems Optimization Laboratory, Department of EESOR, Stanford University, California, February 1999.

[10] DANTZIG G. B., *Linear Programming and Extensions* (Princeton University Press, Princeton, NJ) 1963.

[11] ILOG, *ILOG CPLEX 8.0 User's Manual*, ILOG S.A., Gentilly, France, 2002.

[12] MCCORMICK G. P., *Math. Program.*, **10** (1976) 146.

[13] AL-KHAYYAL F. A. and FALK J. E., *Math. Oper. Res.*, **8** (1983) 273.

[14] LIBERTI L. and PANTELIDES C. C., *J. Global Opt.*, **36** (2006) 161.

[15] LIBERTI L., *Int. Transact. Oper. Res.*, **11** (2004) 34.

[16] LIBERTI L., *J. Global Opt.*, **33** (2005) 157.

[17] SHERALI H. D. and ALAMEDDINE A., *J. Global Opt.*, **2** (1992) 379.

[18] LAVOR C., CARDOZO T. M. and NASCIMENTO M. A. C., *Inter. J. Quantum Chem.*, **103** (2005) 500.

[19] MLADENOVIĆ N., PETROVIĆ J., KOVAČEVIĆ-VUJČIĆ V. and ČANGALOVIĆ M., *Eur. J. Oper. Res.*, **151** (2003) 389.

[20] LIBERTI L. and DRAŽIC M., *Variable neighbourhood search for the global optimization of constrained NLPs*, in *Proceedings of GO Workshop, Almeria, Spain, 2005*.

[21] KUCHERENKO S. and SYTSKO YU., *Comput. Opt. Appl.*, **30** (2004) 297.

[22] LIBERTI L., TSIAKIS P., KEEPING B. and PANTELIDES C. C., *ooOPS*, Centre for Process Systems Engineering, Chemical Engineering Department, Imperial College, London, UK, 2001.