Optimization of Geometry at Hartree-Fock Level Using the Generalized Simulated Annealing

This work presents a procedure to optimize the molecular geometry at the Hartree-Fock level, based on a global optimization method—the Generalized Simulated Annealing. The main characteristic of this methodology is that, at least in principle, it enables the mapping of the energy hypersurface as to guarantee the achievement of the absolute minimum. This method does not use expansions of the energy, nor of its derivates, in terms of the conformation variables. Distinctly, it performs a direct optimization of the total Hartree-Fock energy through a stochastic strategy. The algorithm was tested by determining the Hartree-Fock ground state and optimum geometries of the H2, LiH, BH, Li2, CH, OH, FH, CO, CH, NH, OH and O2 systems. The convergence of our algorithm is totally independent of the initial point and do not require any previous specification of the orbital occupancies.


Introduction
The global optimization problem is a subject of intense current interest.Stochastic optimization methods have been utilized to solve this kind of problem.Essentially these methods consist of performing a direct optimization of a given function, denominated cost function, E, within a determinated stochastic strategy.The Monte Carlo method (MC) is a well-known example of this kind of method.It has been proposed by Metropolis and Ulam [1] whose presented it as a general purpose tool 1 .In the stochastic strategy applied by Metropolis and Ulam, it is used a function, g, to calculate the visiting probability of the hypersurface definited by the cost function, In the original concept of the MC method, the system configurations randomly chosen and used in the calculations of E, are supposed equiprobable.Metropolis [2] proposed a modification in the MC algorithm, i.e. he gave distinct weights to distinct configurations.This method is known in the literature concerned as the Metropolis method.Kirkpatrick et al. [3,4] proposed a new procedure denominated the Simulated Annealing (SA) method, which is a modification of the Metropolis method [2].In the SA method, g is a gaussian function and the T parameter is no longer considered a constant and changes according to     0 log 1 T t T t   , where t enumerates the cycles of the process.In the literature this method is referred to as the Classical Simulated Annealing (CSA) method or Boltzmann Machine.Szu and Hartley [5] proposed a modification in the CSA method where the g function is a Cauchy-Lorentz function and T varies according to  .The Szu and Hartley procedure became known in the literature as the Fast Simulated Annealing (FSA) method or Cauchy Machine.These SA methods have been applied in distinct situations such as restoration of degraded images [6] and microprocessor circuitry design [4].
The Generalized Simulated Annealing method (GSA) [7], has been developed and includes both procedures, the FSA and CSA, as special cases.The GSA approach uses the Tsallis statistics [8,9] to define the visiting distribution function g and has been applied successfully in the description of a variety of global extremization problems.In the domain of the atomic and molecular aggregates, for example, the discovery of the lowest-energy conformations for biological macromolecules or crystal structures for systems with known composition is a frequent goal.In particular, the GSA approach has been used with success in the prediction of new three-dimen-sional protein structure and protein folding [10,11], fitting the potential energy surface for path reaction and chemical reaction dynamics [12,13], gravimetric problem [14], mechanical properties in alloys [15][16][17], in electronic structure problems [18][19][20], among others.
It is important to point out that those typical methodologies used to treat optimization problems based on solving nonlinear necessary-condition equations do not guarantee the achievement of the absolute minimum.This is the case with some variational electronic structure methods, for instance, the Hartree-Fock (HF), multiconfiguration selfconsistent field (MCSCF), molecular geometry determination problems and the correspondent methods in the scope of the Nuclear-Electronic Orbital theory (NEO) [21][22][23][24][25][26].Moreover, it should be observed that the absolute minimum of the functional energy in a given class of functions is the best description of the ground state, as the energy is concerned, within that given class.These observations suggest the importance of developing direct optimization methods for studying these classes of extremal problems.
In previous works [18,19], the GSA algorithm was used to study the problem of determining the absolute minimum of the restricted Hartree-Fock-Roothaan (RHF) [27] and of the unrestricted Hartree-Fock-Pople-Nesbet (UHF) [28] functionals.In another work the GSA algorithm was applied to construct atomic bases [20].The method presented in this work is also based on the GSA algorithm, and it is used to determine the absolute minimum and optimum geometry at the Hartree-Fock (HF) level.This geometry optimization method (hereafter referred to as HF g , RHF g or UHF g ) was tested by determining the HF ground state and optimum geometries of the H 2 , LiH, BH, Li 2 , CH + , OH − , FH, CO, CH, NH, OH and O 2 molecules, using minimal, double-zeta and double-zeta with polarization basis functions (d functions for Li, B, F and p functions for H).The main characteristic of this methodology is that it enables the mapping of the potential energy hypersurface in order to guarantee, at least in principle, that the absolute minimum of the functional in focus is achieved.This methodology does not use expansions of the energy, or of its derivatives, in terms of the conformation variables [29,30].Distinctly, a direct optimization is performed of the total Hartree-Fock energy function through a stochastic strategy, the GSA method.A detailed discussion about the multiple HF extrema, the HF absolute minimum and the GSA algorithm can be found in [18,[31][32][33][34][35][36][37].

The Real HF g Functional and the Constraint Equations
Since the Roothaan and Pople-Nesbet problems are very well known and documented in the literature, only the more general UHF g functional and constraint conditions will be presented.Consider a molecular system with nu- and n n n , , , , , , , , , , , . The number m of atomic functions must satisfy the condition m n   and m n   .The electronic energy functional (in atomic units), in the real UHF approximation, is given by, with the constraint conditions given by, 0 ; 1, , ; , .
In the above equations are the usual overlap, kinetic energy plus nuclear attraction, and electronic repulsion integrals, respectively, that depend of the nuclear coordinates     X R .For the HF g functional, the total energy E, is given by,

The HF g Algorithm
As in references [18][19][20], the GSA algorithm used here includes two additional modifications relative to the original version described in reference [7].The first one is the introduction of constraint conditions in the structure of the algorithm (steps 4) and 9)) to treat our variational problem of constrained extrema.The second modification was the introduction of a new independent parameter, T , to construct the temperature function defined in the step 6). q The procedure used to search for the global and local minima or to map the cost function hypersurface consists in comparing the total energy for two consecutive values of the and X obtained with the GSA and X, for two consecutive GSA steps, are related to the previous ones via random perturbations on the LCAO-MOs coefficients and on the molecular conformation, respectively.In each cycle, and X are simultaneous and independently generated.

 
In summary, the whole UHF g algorithm for mapping and searching for the global minimum of the total energy surface is: 1) Fix the q A , q V and q T parameters relative to acceptance and visitation probability-distribution functions and temperature function, respectively; 2) Start at t = 1, the first step in the iterative process, with an arbitrary initial matrix guess t , an arbitrary molecular conformation and a high enough value for the "temperature" ; 3) Calculate the integrals , , ,  according to Equations 1; 5) Calculate the total energy using Equation 2; Calculate a new temperature as follows [7]: , where and , are randomly generated by using the visiting probability distribution [7] to obtain the new matrix t and the new molecular conformation according to Equations 1; 10) Calculate the total energy using Equation 2. The new energy value will be accepted or not according to the rule: and by , , if , the acceptance probability [7] defined by retain t and , otherwise, replace and by and ;

11) Take
and return to 6) until the convergence of After convergence is achieved, the orbital energies 1 , , , , is the Fock's matrix constructed with the converged UHF g occupied matrix .Also, it is always possible to obtain the virtual canonical vectors,  and the respective virtual orbital energies Note that, while for the standard RHF/UHF-SCF [27,28] calculations one needs to specify, a priory, the orbital occupancy, no ad hoc orbital occupation rule is needed for the RHF g and UHF g calculations The following stopping criterion was adopted for the HF g iterative process convergence was assumed if the difference between the current total energy value and the lowest total energy previously obtained during the process was less than a pre-established value  E   for a certain number of consecutive steps (nstop).The HF g calculations were performed in atomic units and was used and .The algorithm HF g described above is illustrated in the flowchart of Figure 1.
Calculate the new temperature ( ) Calculate the overlap, core and 2e integrals at geometry t X Ortho-normalize T A V T q q q q T , nstop, ΔE and setup initial matrix   Figure 1.HF g -GSA Flowchart.

Discussion
To test the HF g method, we calculated the HF ground state energy and optimum geometry for the H 2 , LiH, BH, Li 2 , CH + , OH -, FH, CO, CH, NH, OH and O 2 molecules.The calculations were carried out using minimal (STO-6G), double-zeta, and double-zeta with polarization functions (d functions for Li, B, F and p functions for H) basis sets.Tables 1-3 show the point and spin symmetry classes of the ground state, the kind of calculations performed, the geometry and correspondent energy obtained and the experimental geometry extracted from the Herzberg book [38].We performed several RHF g and UHF g calculations combining different initial values with distinct sets of the parameters q A , q V , q T and T 0 .We found that the narrow ranges of values of the parameters q V and q T , leading to a better convergence (smallest number of HF g cycles), namely, T q e similar to those obtained in previous works [18,19].In particular, for the minimal basis set, the best convergence is achieved for 2.9 V q  ar  and q 1.9 T  , which are quite close to the best values obtained for the RHF-GSA and UHF-GSA problems [18,19].In addition, we performed several calculations using different values of q A , including acc 0 A  , wh the convergence, therefore, been adopted in step 10) of the HF g algorithm, an acceptance probability, acc ich led to A , equal zero for all calculations.The genera to onvergence be l g algorit c havior of the HF hm is similar to that of the RHF-GSA and UHF-GSA methods [18,19].For all systems and bases sets employed, it was always possible to obtain the global minimum, with several distinct combinations of these parameters, each set of parameters requiring a different number of HF g cycles.Also in all calculations, the HF g energies initially show a strong oscillatory behavior but soon afterwards the energy starts to smoothly converge towards the absolute minimum.Figures 2 and 3 present the RHF g and UHF g convergence profiles for the CH + and OH molecules, indicating the values of the parameters acc A , q V , q T , T 0 , the atomic basis sets, the type of guess r the initial values of    .Similar convergence profiles were obtained for thers molecules.In order to verify the accuracy of the calculation all the o s, the RH F g and UHF g results were compared with those obtained by the standard gradient RHF/UHF geometry calculation method, for all molecules considered, using the  program GAMESS [39].Three choices for the initial ma- The Monte Carlo trix  were considered when performing the RHF/ UHF AMESS calculations: the eigenvectors of core hamiltonian (H core ) 4 , the eigenvectors of an extended Huckel calculations (Huckel) 5 , and the eigenvectors of a previous RHF g or UHF g calculation.Furthermore, whenever necessary, the Direct Inversion in the Iterative Subspace (DIIS) convergence acceleration technique [40,41] was also used for the RHF/UHF-GAMESS calculations.In all cases we examined, the RHF g and UHF g calcula ns converge to the global minimum with any randomly generated initial  and X values, what is not observed for the RHF/UH GAMESS calculations.Besides, the HF g method do not need any previous specification of the orbital occupancies.

Figure 3 .
Figure 3. OH double-zeta basis HF g convergence process.