Leo Liberti
email: liberti@lix.polytechnique.fr

Carlile Lavor
email: clavor@ime.unicamp.br

Nelson Maculan May

N Maculan
email: maculan@cos.ufrj.br

A mathematical programming formulation for the Hartree-Fock Problem on open-shell systems

Keywords: quantum chemistry, mathematical programming, linear independence, rank constraint

 we proposed a new Mathematical Programming (MP) approach which enhanced the likelihood of attaining globally optimal approximations, limited to closed-shell atomic systems. In this paper, we discuss an extension to open-shell systems: this is nontrivial as it requires the expression of a rank constraint within an MP formulation. We achieve this by explicitly modelling eigenvalues and requiring them to be nonzero. Although our approach might not necessarily scale well, we show it works on two open-shell systems (lithium and boron).

Introduction

The wave function of a quantum many-body system of electrons in atoms and molecules satisfies the time-independent Schrödinger equation:

HΦ = EΦ. ( 1
)
This equation is the projection on the position/momentum variables of the time-dependent Schrödinger equation [START_REF] Schrödinger | Quantisierung als Eigenwertproblem[END_REF], which describes the changes over time of a quantum physical system. The time-dependent Schrödinger equation describes the wave nature of particles and explains phenomena such as quantum jumps, quantum state superposition, and probabilistic interpretations of particle position and momentum [START_REF] Penrose | The road to reality[END_REF]. The time-independent version given in Eq. ( 1) describes the standing waves corresponding to systems in some form of equilibrium. H is a differential operator which represents the total energy of the system. Solutions to Eq. ( 1) are usually impossible to derive analytically (aside from the very simplest cases), so approximate numerical solutions are usually sought instead. The Hartree-Fock method is used to approximate Φ by expressing it as a linear combination of finitely many basis functions. We propose a Mathematical Programming (MP) formulation to decide the value of these coefficients. More specifically, the approximations of Φ are anti-symmetrized products of a set { φi | i ≤ n} of one-electron wave functions called spatial orbitals, each of which is approximated by a linear combination of a finite subset B ⊆ N of a complete basis set {χ r } ∞ r=1 of the Hilbert space. The larger |B| is, the better the approximation is likely to become. See [18, §1] for more information.

The Hartree-Fock problem (HFP) was defined in [START_REF] Lavor | Solving Hartree-Fock systems with global optimization methods[END_REF]. To summarize, we want to find a set of coefficients c ri , for r ∈ B and i ∈ N = {1, . . . , n} ⊆ B = {1, . . . , b}, such that:

∀i ∈ N ϕ i = ∑ r∈B c ri χ r (2)
are good approximations of the spatial orbitals φi . The quality of the approximation is inversely proportional to the energy function E associated with the approximating set {ϕ i }, which is guaranteed to be an upper bound to the energy function associated with the spatial orbitals [START_REF] Levine | Quantum Chemistry[END_REF]. The most common procedure used to solve the HFP is the Self-Consistent Field (SCF) methodology [START_REF] Szabo | Modern Quantum Chemistry[END_REF]. This is essentially an iterative optimization algorithm that takes an initial solution as input and improves it until a local optimum is reached. As such, its main limitations are that: (a) it needs a good starting point, and (b) it may fail to find the global optimum.

The expression for the energy function E varies according to whether the considered system is open-or closed-shell. Unlike a closed-shell system, an open-shell system has a molecular orbital with a valence shell which is not completely filled with electrons, or that has not given all of its valence electrons through chemical bonds (see en.wikipedia.org/ wiki/Open_shell).

The HFP was first cast as a nonconvex Nonlinear Program (NLP) in [START_REF] Lavor | Using an interval branch-and-bound algorithm in the Hartree-Fock method[END_REF]: in that paper, an interval Branch-and-Bound (iBB) method was proposed to find {ϕ i } for the He (helium) and Be (beryllium) closed-shell atomic systems. In [START_REF] Lavor | Solving Hartree-Fock systems with global optimization methods[END_REF][START_REF] Liberti | Reformulation in mathematical programming: an application to quantum chemistry[END_REF], a different type of Branch-and-Bound (BB) algorithm, called spatial BB (sBB), was tested on the same atomic systems He and Be. Both iBB and sBB provide an optimality guarantee to within a given ε > 0 tolerance, but sBB turns out to be more computationally efficient than the iBB on the HFP. More guaranteed results (H 2 , LiH and H + 3 ) were obtained in [START_REF] Janes | Rigorous numerical approaches in electronic structure theory[END_REF][START_REF] Janes | Deterministic global optimization in ab-initio quantum chemistry[END_REF] using a specially adjusted version of the new, open-source sBB solver COUENNE [START_REF] Belotti | Branching and bounds tightening techniques for non-convex MINLP[END_REF]. Some work in the same direction was also carried out in [START_REF] Zorn | Hartree fock self-consistent calculations: Global optimization of electronic structure[END_REF]. All the MP-based methods reviewed above refer to closedshell atomic systems: the main properties of such systems are (a) they can be approximated by expressing them as linear combinations of a finite-dimensional function basis (b) they minimize a certain energy function and (c) the set {ϕ i } turns out to be an orthonormal set (see [START_REF] Lavor | Solving Hartree-Fock systems with global optimization methods[END_REF] and [19, §3.5]). These properties are exploited in MP formulations as summarized in Sect. 2 below.

The present paper proposes an MP-based method to address the extension of the HFP to open-shell systems. In such systems, the orthonormality of {ϕ i } is not a physical requirement: in fact, all that is required is linear independence [START_REF] Martin | Electronic Structure: Basic Theory and Practical Methods[END_REF]. Orthonormality is usually imposed to simplify the SCF calculations [START_REF] Hollauer | A generalized multistructural wave function[END_REF]. A consequence of this fact is that the SCF methodology is even more likely to find locally optimal solutions which fail to be global. Linear independence restrictions amount to a rank constraint on a matrix of decision variables, and is very difficult to express in MP formulations. A general modelling methodology was proposed in [START_REF] Liberti | Modelling rank constraints in mathematical programming[END_REF]. In this paper we derive a new MP formulation for the HFP on openshell systems (Sect. 3), and solve it with a combination of reformulations [START_REF] Liberti | Reformulations in mathematical programming: Definitions and systematics[END_REF] and Global Optimization (GO) methods (Sect. 4) on the B (boron) and Li (lithium) atomic systems (Sect. 5).

Formulation for closed-shell systems

In this section we recall from [START_REF] Lavor | Solving Hartree-Fock systems with global optimization methods[END_REF][START_REF] Liberti | Reformulation in mathematical programming: an application to quantum chemistry[END_REF] the MP formulation for the HFP on closed-shell systems, which is based on orthogonality constraints on the orbitals. We employ an approximating basis set X = {χ r | r ∈ B}. We consider three sets of parameters α, β , γ related to the energy of the considered system: α depends on four indices r, s,t, u ∈ B, β depends on two indices r, s ∈ B and γ is a scalar. Since we have to decide how the approximated orbitals ϕ i depend on the approximating basis X, by (2) we have to find an assignment of feasible values for the coefficients c ri (for r ∈ B, i ∈ N). In other words, c ri are the decision variables of our MP formulation. The molecular energy is given in [START_REF] Lavor | Using an interval branch-and-bound algorithm in the Hartree-Fock method[END_REF][START_REF] Lavor | Solving Hartree-Fock systems with global optimization methods[END_REF]:

E closed (c) = 2 ∑ r,s∈B ∑ i∈N c ri c si ∑ t,u∈B α tu rs ∑ i∈N c ti c ui + β rs + γ,
so that the objective function of the closed-shell formulation is:

min c∈R bn E closed (c),
where b = |B|. We remark that E closed (c) refers to the energy of closed-shell systems assuming that {ϕ i | i ∈ N} is an orthonormal set. We therefore impose orthonormality constraints:

∀i ≤ j ∈ N ϕ i |ϕ j = δ i j ,
where δ i j is the function with value 1 if and only if i = j. Replacing ϕ i = ∑ r∈B c ri χ r for all i ∈ N and χ r |χ s = S rs for all r, s ∈ B, we obtain a constraint expressed in function of the decision variables c:

∀i ≤ j ∈ N ∑ r,s≤b S rs c ri c s j = δ i j .
We remark that S rr = 1 for all r ∈ B. Further constraints to the problem are given by the decision variable bounds:

c L ≤ c ≤ c U ,
which depend on the instance. Common values for these bounds are ±1.5 or ±2. The HFP can thus be formulated as follows [START_REF] Lavor | Solving Hartree-Fock systems with global optimization methods[END_REF][START_REF] Liberti | Reformulation in mathematical programming: an application to quantum chemistry[END_REF]:

min c∈R bn E closed (c) ∀i ≤ j ∈ N ∑ r,s≤b S rs c ri c s j = δ i j c L ≤ c ≤ c U ,        (3)
i.e. the minimization of a quartic objective function subject to quadratic constraints and variable bounds.

3 Formulation for open-shell systems Orthonormality of the {ϕ i } set derives from physical considerations in closed-shell systems, but is only imposed to simplify SCF computations in open-shell systems [START_REF] Martin | Electronic Structure: Basic Theory and Practical Methods[END_REF][START_REF] Hollauer | A generalized multistructural wave function[END_REF]. In such systems, {ϕ i } is simply required to be a linearly independent set. The molecular energy expression changes accordingly [START_REF] Cardozo | Private communication[END_REF]:

E open (c) = ∑ i∈N r,s∈B W i β rs c ri c si + ∑ i, j∈N r,s,t,u∈B A i j α tu rs c ri c si c t j c u j -∑ i, j∈N r,s,t,u∈B D i j α tu rs c ri c s j c t j c ui , (4)
where W is a given constant n-vector and A, D are given constant n × n matrices, where n = |N|. Moreover, since we no longer impose the orthonormality of {ϕ i }, we need to explicitly require that each ϕ i be normalized:

∀i ∈ N ∑ r∈B S rr c 2 ri = 1.
Modelling linear independence of a general set of vectors

Y = {ϕ i ∈ R b | i ∈ N}
, where b = |B|, within a MP formulation is far from easy. Specifically, the direct formulation:

∀λ ∈ R n ∑ i∈N λ i ϕ i = 0 (5)
is an uncountable set of constraints, each of which defines an open set in R b , i.e. the entire space minus a proper affine subspace. The uncountable intersection of open sets is certainly not guaranteed to be closed, which means that not only the corresponding MP formulation would have a practically unmanageable number of constraints, but also that no minimum might exist within the feasible set described by it. It is relatively easy to find examples that show that replacing [START_REF] Fletcher | User manual for FILTER[END_REF] with any finite subset of its constraints fails to guarantee linear independence. An altogether different idea stems from the following equivalent condition to the linear independence of Y : the Cayley-Menger determinant [START_REF] Blumenthal | Theory and Applications of Distance Geometry[END_REF] relative to Y in R n has nonzero value if and only if Y is linearly independent. This is equivalent to requiring that the (n -1)-dimensional simplex defined by the projection of Y onto R n-1 should have positive volume. Since explicit formulae for computing an n×n determinant are only valid for a given constant n, or else are recursive, it is not clear how to express this condition in practice in a MP formulation in such a way as to be valid for any n. Moreover, even if a valid expression were found, in general it would yield a polynomial MP of degree n, which is considered to be very hard to solve for even moderately high values of n.

In general, the linear independence constraint is a special case of the more general "rank constraint", i.e. requiring that a matrix of decision variables should have a given rank r. Although a famous relaxation of this constraint yields the well-known Semidefinite Programming (SDP) problem [START_REF] Rendl | Semidefinite relaxations for partitioning, assignment and ordering problems[END_REF], the paper [START_REF] Liberti | Modelling rank constraints in mathematical programming[END_REF] is, to the best of our knowledge, the only existing reference treating this problem explicitly. The formulation given in [START_REF] Liberti | Modelling rank constraints in mathematical programming[END_REF] is a nonconvex quadratic MP formulation whose constraints require the projection of Y on R n to yield a square symmetric matrix with n nonzero eigenvalues. In the present setting, after replacing ϕ i with the corresponding expression in terms of the c (see Eq. ( 2)), we obtain the following eigensystem:

∀k, i ∈ N ∑ j∈N r,s∈B S rs c ri c s j v k j = λ k v ki , (6)
where {λ k ∈ R | k ∈ N} is the set of eigenvalues, and {v k ∈ R n | k ∈ N} is the set of corresponding eigenvectors. We remark that since the matrix is expressed in terms of decision variables c, both λ and v are also decision variables of the new MP formulation. Without loss of generality, we require the eigenvectors to be an orthonormal set:

∀k ≤ h ∈ N ∑ j∈N v k j v h j = δ kh .
Finally, we introduce a constant ε > 0 and require that all squared eigenvalues should have greater value: ∀k ∈ N λ 2 k ≥ ε. This ensures all eigenvalues are nonzero, which is equivalent to Y being a linearly independent set.

In summary, the following:

min c,λ ,v E open (c) ∀i ∈ N ∑ r∈B S rr c 2 ri = 1 ∀k, i ∈ N ∑ j∈N r,s∈B S rs c ri c s j v k j = λ k v ki ∀k ≤ h ∈ N ∑ j∈N v k j v h j = δ kh ∀k ∈ N λ 2 k ≥ ε c L ≤ c ≤ c U                          (7)
is a valid formulation for the HFP on open-shell systems for any ε > 0. It consists of the minimization of a quartic polynomial subject to a set of cubic and quadratic polynomials. We remark that since orthonormal vectors are also linearly independent normalized vectors, the feasible set of (7) contains the feasible set of (3).

Reformulations and solution methods

Both formulations (3) and ( 7) belong to the class of nonconvex NLPs: they are multiextremal, involve nonlinear terms in both the objective function and the constraints, and are defined over continuous variables only. This calls for GO solution methods. We had previously [START_REF] Lavor | Using an interval branch-and-bound algorithm in the Hartree-Fock method[END_REF][START_REF] Lavor | Solving Hartree-Fock systems with global optimization methods[END_REF][START_REF] Liberti | Reformulation in mathematical programming: an application to quantum chemistry[END_REF] employed spatial Branch-and-Bound (sBB) algorithms in order to solve the small closed-shell systems He and Be: the advantage of the sBB is that it provides a guarantee of ε-approximation for any given ε > 0. For the larger atomic systems tested in this paper, however, the sBB is too slow. Some computational experience on HFP instances showed that the closed-shell MP formulation (3) is numerically rather well-behaved notwithstanding its nonconvexities: a few well-known local NLP solvers, such as SNOPT [START_REF] Gill | User's Guide for SNOPT 5.3. Systems Optimization Laboratory[END_REF] and IPOPT [START_REF]COIN-OR: Introduction to IPOPT: A tutorial for downloading, installing, and using IPOPT[END_REF], usually identify the global optima if restarted from a few different starting points. This suggests that a simple MS metaheuristic [START_REF] Schoen | Two-phase methods for global optimization[END_REF] is appropriate to solve [START_REF] Cardozo | Private communication[END_REF].

The open-shell MP formulation [START_REF] Gill | User's Guide for SNOPT 5.3. Systems Optimization Laboratory[END_REF], however, poses severe numerical problems to all local NLP solvers we tried (SNOPT, IPOPT, Filter [START_REF] Fletcher | User manual for FILTER[END_REF], MINOS [START_REF] Murray | GAMS/MINOS: A solver for large-scale nonlinear optimization problems[END_REF] and DONLP2 [START_REF] Spellucci | DONLP2 Users Guide[END_REF]). A little testing showed that solvers had a hard time trying to satisfy [START_REF] Fourer | The AMPL Book[END_REF]. We therefore reformulated [START_REF] Gill | User's Guide for SNOPT 5.3. Systems Optimization Laboratory[END_REF] by adding slack variables s ki (for k, i ∈ N) to (6) whilst minimizing their squares added to the objective (4). This results in the following formulation:

min c,λ ,v,s E open (c) + ∑ k,i∈N s 2 ki ∀i ∈ N ∑ r∈B S rr c 2 ri = 1 ∀k, i ∈ N ∑ j∈N r,s∈B S rs c ri c s j v k j = λ k v ki + s ki ∀k ≤ h ∈ N ∑ j∈N v k j v h j = δ kh ∀k ∈ N λ 2 k ≥ ε c L ≤ c ≤ c U ,                          (8)
whose numerical behaviour with respect to local NLP solvers is improved. Naturally, unless the computed optimum (c * , λ * , v * , s * ) of ( 8) has ∑ k,i∈N s * ki = 0, we cannot be sure that linear independence is attained. In practice, however, this is not a problem, as linear independence of the solution can be verified a posteriori. Should the solution fail to be linearly independent, a new run of the (stochastic) MS algorithm can be launched. Because the MS converges to the global optimum with probability 1 in infinite time [START_REF] Solis | Minimization of random search techniques[END_REF], the chance that all computed optima turn out to be linearly dependent is low.

It is interesting to remark that (8) exhibits some symmetry: with the computational methods discussed in [START_REF] Liberti | Reformulations in mathematical programming: Automatic symmetry detection and exploitation[END_REF] and the generalization from instance to problem proposed in [START_REF] Liberti | Symmetry in mathematical programming[END_REF], the formulation group of ( 8) is found to be the full symmetric permutation group S n (where n = |N|), and corresponds to the fact that no constraint in (8) enforces the eigenvalues λ k to be ordered. Adjoining constraints

λ 1 ≤ . . . ≤ λ n (9)
to [START_REF] Hollauer | A generalized multistructural wave function[END_REF] to the formulation yields a reformulation of the narrowing type (see [START_REF] Liberti | Reformulations in mathematical programming: Definitions and systematics[END_REF]) with trivial group. This is generally useful when the solution approach is of the sBB type. For metaheuristics, as remarked in [START_REF] Liberti | Symmetry in mathematical programming[END_REF], it is better to keep formulation symmetry.

Computational results

We implemented the formulations (3) and ( 8) using the MP language AMPL [START_REF] Fourer | The AMPL Book[END_REF] and solved four closed-shell and two open-shell instances using the MS approach with both SNOPT and IPOPT. It is worth remarking that there is nothing complicated about our MS: it simply computes random starting points, then runs the SNOPT local solver (which is, itself, a very sophisticated piece of mathematical software) a finite number I of times (in all our experiments, I = 10). We want to convey the message that the model-and-solve paradigm we are using is advantageous: a sophisticated model may well drive a general-purpose solver to an early and successful termination, thereby yielding an efficient overall method. The accent is on modelling the problem rather than devising and implementing a solution algorithm.

In Table 1 we recorded (in column order): the atomic symbol, the number V of variables in the formulation, the number C of constraints in the formulation, the value of the optimal energy E * expressed in atomic units (au), the seconds of user CPU time computed using a single core of an Intel Core i7-7500 CPU running at 2.70GHz with 16GB RAM, and the MS iteration (over I = 10) at which the putative global optimum was found. The energy values E * are not as precise as those given in the literature. This is a consequence of two factors. First, the data we used (in particular the values for the S rs symbols) was obtained through a simplified model. Second, the MP formulations we employ allow a 1 × 10 -3 tolerance in order to improve the numerical stability of the associated NLP solvers. Because the energy functions E closed and E open depend quartically on the decision variables, the imprecision in the objectives generally corresponds to smaller imprecisions in the values of the c variables. This need not be a limitation: the solutions we obtain are intended to provide very high quality starting points for further refinement using the established SCF methodology. In order to validate our numerical results, the MP formulations were solved by MS using very different local solvers, namely SNOPT (our main solver, to which the CPU times refer), IPOPT, KNitro, Filter, MINOS, Loqo. The (putative) global optima were confirmed by all solvers (except for some nonconvergence failures of MINOS on Ne and Mg). In Table 1, however, we only reported results relative to the SNOPT solver, which we found to be the most efficient (with Knitro and Loqo) for this data set.

Conclusion

In this paper, we extended the Mathematical Programming approach for solving the Hartree-Fock equations to open-shell systems. This methodology was first proposed in [START_REF] Lavor | Solving Hartree-Fock systems with global optimization methods[END_REF] limited to closed-shell systems. The present extension is innovative as far as quantum chemistry is concerned, because it dispenses with a limitation imposed by numerical requirements of the SCF methodology, i.e. orthogonality of the ϕ i s. As far as mathematical programming is concerned, the approach is innovative because insofar as it discusses the "linear independence constraint", which is sometimes useful but not easy to encode in MP terms.

Table 1

 1 Instances at www.lix.polytechnique.fr/~liberti/hartreefock-instances.zip.

	Atom	Name	V	C	E * (au) CPU (s) Opt. It.
			Closed-shell systems
	He	helium	2	1 -2.74706	0.01	8
	Be	beryllium	4	3 -14.3519	0.02	8
	Ne	neon	25 15 -126.603	0.11	6
	Mg	magnesium 54 21 -196.687	3.95	7
			Open-shell systems
	Li	lithium	20 11 -12.2584	0.08	9
	B	boron	36 23	-57.334	0.17	4

Acknowledgments

The authors are thankful to Brazilian research agencies FAPERJ, FAPESP, CAPES and CNPq for financial support, and to M.A. Chaer Nascimento and T.M. Cardozo for useful discussions leading to the formulation of Eq. ( 4). The first author (LL) has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement n. 764759 "MINOA".