New Real-Coded Genetic Algorithm Operators for Minimization of Molecular Potential Energy Function

Abstract

The global minimum of the potential energy of a molecule corresponds to its most stable conformation and it dictates most of its properties. Due to the extensive search space and the massive number of local minima that propagate exponentially with molecular size, determining the global minimum of a potential energy function could prove to be significantly challenging. This study demonstrates the application of newly designed real-coded genetic algorithm (RCGA) called RX-STPM, which incorporates the use of Rayleigh crossover (RX) and scale-truncated Pareto mutator (STPM) as defined earlier for minimizing molecular potential energy functions. Computational results for problems with up to 100 degrees of freedom are compared with five other existing methods from the literature. The numerical results indicate the underlying reliability (robustness) and efficiency of the proposed approach compared to other existing algorithms with low computational costs.

INTRODUCTION

Molecular mechanics (MM) is a mathematical model used to predict the accurate structures and energies of molecules. Its fundamental notion is aimed at optimizing the geometry of a molecule through computed bonds and bond angles to assume its natural meanings and values. The MM energy expression explains the correlation between chemical structure and its energy. Having said that, one of the setbacks in the field of computational chemistry and applied mathematics is to locate the most stable conformation of a molecule with the minimum potential energy, as derived by a mathematical function. Determining the most stable conformation of a molecule can be formulated as a global optimization and belongs to a class of NP problems. Consequently, this problem has recently gained recognition among researchers in the field of optimization and computer science. The potential energy function is a nonconvex expression that has many local minima. The core challenge of this problem is that the number of local minima grows exponentially with molecular size (Wales and Scheraga 1999).

Based on a review of past studies, both deterministic and stochastic approaches have been applied to this problem. Deterministic approaches ensure a theoretical guarantee of finding the global minimum if not the local minimum whose objective function value differs by at worst ε from the global one for a given ε > 0, whereas stochastic techniques offer only a guarantee in probability as it varies in the search procedure. On one hand, the deterministic techniques used include diffusion equation methods, distance scaling methods, and packet annealing techniques (Piela, Kostrowicki, and Scheraga 1989; Kostrowicki and Scheraga 1992; Pillardy et al. 2001). To streamline the problem, the approaches mentioned attempt to change the shape of the potential energy surface to reduce the number of local minima. Examples of other deterministic techniques used are αBB (branch and bound; Maranas and Floudas 1994a, 1994b) and interval analysis (Lavor 2003; Lin and Stadtherr 2005). Stochastic methods used to solve this problem are simulated annealing (Kawai, Kikuchi, and Okamoto 1989; Moret et al. 1998), Tabu search, genetic algorithms (GAs), real-coded GAs (RCGAs), and hybrid GAs (Ali and Hassanien 2013; Brodmeier and Pretsch 1994; Barbosa et al. 2000; Barbosa, Lavor, and Raupp 2005; Deep et al. 2011; Hedar, Ali, and Hamid 2011), swarm intelligence algorithm (Agrawal and Silakari 2014; Bansal, Shashi, and Katiyar 2010), Monte Carlo approaches (Li and Scheraga 1988), and continuous neighborhood search (Dražić et al. 2008).

This study is aimed at investigating the effect of a newly developed RCGA called Rayleigh crossover scale-truncated Pareto mutator (RX-STPM) in solving a scalable, highly complex molecular potential energy function for up to 100 degrees of freedom with known global minimum. This function is minimized based on an independent set of internal coordinates that involve only torsion angles. The performance of RX-STPM was tested based on the performance criteria of reliability and efficiency. The experimental results are extensively compared with those from GA.

This article is presented as follows: section “The Problem: Optimization of Molecular Structures” presents the mathematical formulation of the molecular potential energy function; section “Proposed Real-coded Genetic Algorithms” describes the proposed RCGA for solving the problem; section “Computational Steps and Experimental Setup” discusses the computational steps, experimental setup, and performance evaluation criteria used by the proposed algorithms; section “Results and Discussion” includes discussions and explanation of computational results; and conclusions are presented in section “Conclusion.”

THE PROBLEM: OPTIMIZATION OF MOLECULAR STRUCTURES

The pro and contra of many different force fields used for modeling molecular structures are reported in Allinger, Yuh, and Lii (1989), Hamida and Petrowski (2000), Mühlenbein, Schomisch, and Born (1991), and Neumaier (1997). Lavor’s potential energy function (Lavor and Maculan 2004) is used as a simplified molecular model to represent the force field in this study. Changes in the internal motions such as the bending of angles, stretching of bonds, and rotation about a single bond help generate the potential energy of a molecule. The global minimum of the molecular potential energy function is linked to the lowest energy structure of a molecule. The significance of determining this structure is to discern most of the physical and chemical properties of the molecule. This model can be illustrated through a linear chain of n beads centered at X₁, …, X_n in a three-dimensional space. The total potential energy is given as

E₁ is the contribution of bond stretching energy for every pair of consecutive beads, X_i and X_i+1. r_i,_i+1 refers to the bond length, which is the Euclidean distance between the beads. is the bond force constant, and is the equilibrium bond length.

E₂ is the contribution of angle bending energy over every three consecutive beads, X_i, X_i+1, X_i+2. Let θ_i, _i+2 be the bond angle corresponding to the relative position of the third bead with respect to the line containing the previous two. is the angle bending force constant, and is the equilibrium bond angle.

E₃ is the contribution due to the dihedral angle over every four consecutive beads, X_i, X_i+1, X_i+2, X_i+3. Let ω_i,_i+3 be the angle, called the torsion angle, between the normals through the planes determined by the beads X_i, X_i+1, X_i+2 and X_i+1, X_i+2, X_i+3. is the torsion force constant and is the equilibrium angle. M_i, i = 1, 2, 3, represents the set of pair of atoms separated by i covalent bonds.

The two-body interaction between every pair of beads separated by more than two covalent bonds along the chain is characterized by E₄. r_ij is the Euclidean distance between the beads X_i and X_j. The challenge lies in minimizing the total molecular potential energy, E = E₁ + E₂ + E₃ + E₄, correspond to finding the optimal spatial configuration of the beads. In this problem, all of the bond lengths and bond angles followed the characteristic of a hydrocarbon chain and are fixed at the equilibrium values = 1.526 Å and = 1.91 and = 1 and = 0. All of the parameters are clearly defined in Lavor and Maculan (2004) and the energy function takes the following form:

The trick is in finding the global minimum of this potential function, which depends only on the dihedral angles ω_ij, (i, j) ∈ M₃. For the N-atom case, the global minimizer can be traced analytically by alternating the sequence of dihedral angles 1.039195303 and π, i = 1, 3, 5, …, N. For any size of n, the difference between the global minimum value, E_x, and the second lowest value, E_y, will always satisfy the following relation:

This simplified molecular model is a feasible yet challenging test function. Given the constrain of ω_ij; 0 ≤ ω_ij ≤ 5, there exists only one global minimum and a huge number (2ⁿ⁻³) of local minima, where n = number of beads in a molecule. For example, a chain with 50 beads has 2⁴⁷ = 1.407374884 × 10¹⁴ local minima.

PROPOSED REAL-CODED GENETIC ALGORITHMS

GA refers to a group of biologically driven optimization techniques that allow a population of individuals who would thrive in the survival of the fittest to evolve while going into the next generation. In regards to this problem, the variables of the optimization problem are the dihedral angles and the objective function value is the potential energy. The central operations of GA are reproduction, crossover, and mutation on a population of chromosomes. Each individual of a population contributes their genotype to their suitability of their expressed phenotype in the form of offspring. Through a process of mating whereby the crossover operator takes two genotypes and combines them (either by merging or by exchanging the values of the genes), the next generation is hence produced. The mutation operator then modifies one or multiple genes. This evolutionary iterative process continues until termination criteria are achieved. GAs that incorporate real number vector representations of chromosomes are called RCGAs. A review of RCGAs is reported in Herrera, Lozano, and Verdegay (1998). RCGA is used in mnay applications and is recommended for optimization problems (Michalewicz 1996; Ono, Satoh, and Kobayashi 1999). There are many advantages of using RCGAs in numerical function optimization (Wright 1990). Over the years, RCGA research has focused on developing new crossover and mutating operators to improve the function optimization and GA performance (Deep and Thakur 2007a; Tang and Tseng 2012).

Rayleigh Crossover

We defined RX based on a Rayleigh distribution (Lim et al. 2014a). This is a continuous probability distribution that randomly populates offspring for a real number GA. The Rayleigh density function is given as

(1)

where s > 0 is the scale parameter of the distribution. By experiment, s is best kept at 3.0. The distribution function is

(2)

To use a Rayleigh distribution, two parents and are taken to produce two offspring and in the following equations:

(3)

(4)

From Eq. (3), offspring is set closer to parent 1, , yet at the same time inherits elements from parent 2, . This is the same for ; offspring is set closer to parent 2, , yet inherits from as well. Log is introduced to set the boundary of x, for |x|, {0 > x > 1}. The Rayleigh-distributed number x is generated by inverting the distribution function of the Rayleigh distribution as follows:

x is suggested to take only the positive values. Hence:

Scale-Truncated Pareto Mutation

We defined STPM based on the formula to generate Pareto random variables (Lim et al. 2014b). The truncated Pareto distribution has three parameters, α, L, and H. α determines the shape, L denotes the lower bound, and H denotes the upper bound of the evaluation function to be optimized. The probability density function is

Applying inverse transformation, the equation for U as a uniformly distributed function is defined as

In addition, is truncated Pareto distributed as

A modulus, | |, is used to eliminate the possible imaginary number produce by the algorithm. To apply the truncated Pareto distribution as the mutation operator, an adjustable scale, , is introduced. is added to make sure that is not overweighted. If is overweighted, there will be a chance of good chromosomes being altered in the process. Thus, the newly mutated offspring m is defined as

where P is the parent, and , to determine the direction of the mutation. By experiment, and α are best kept at and 3.0, respectively.

COMPUTATIONAL STEPS AND EXPERIMENTAL SETUP

Computational Steps

Computational steps of the proposed RX-STPM used are as follows:

• Step 1 (Initialization). Generate an initial population of chromosomes randomly. Set Generation = 0.

• Step 2 (Evaluation). Evaluate the fitness of each individual of the current population.

• Step 3 (Termination check). If any member of the current population meets the termination criteria then stop. Else, go to the next step.

• Step 4 (Selection). Use Tournament selection operator to extract members from the current population to generate a mating pool.

• Step 5 (Crossover and mutation). Crossover and mutate the population in the mating pool with optimum probability to make new population.

• Step 6. Increase generation; replace old population with new population; repeat the loop by going back to Step 2.

Experimental Setup

The performance analysis of RX-STPM and all of the GAs used in the comparative studies on the molecular potential energy function, E, is deliberated in section “The Problem: Optimization of Molecular Structures.” E is minimized in the specified search space [0, 5]ⁿ, where n refers to the total number of beads in a system or the number of dimensions. Table 1 shows the global minimum values for the function E corresponding to 20, 40, 60, 80, and 100 dimensions (Barbosa, Lavor, and Raupp 2005).

TABLE 1 Global minimum values, E, for potential molecular energy function for different dimensions, n = 20 to 100

n	20	40	60	80	100
E	−0.822366	−1.644732	−2.467098	−3.289464	−4.111830

Proper parameter setting in a GA will significantly affect the performance. GA parameters are tuned using decision support system (Pakath and Zaveri 1995), full-factorial experimental design (Gupta, Gupta, and Kumar 1993), forming trial chromosomes at each GA iteration (Leung and Wang 2001; Liu, Tsai, and Chou 2006), etc.

Schaffer et al. (1989) conducted one of the most extensive studies on this and concluded that optimal parameter setting is problem dependent. Thus far, no general consensus has been reached in regards of the optimum parameterization of GA operators. In addition, there are no GA tools that could generalize adjustment of control parameters. Researchers are to design their ideal operators, parameters, and mechanisms to tackle the problems (Digalakis and Margaritis 2001). Ideally, the range for crossover probability (P_x) and mutation rate (P_m) are 0.2 to 1 and 0.0001 to 0.1, respectively. If P_x is too low, the GA may encounter difficulty in searching forward, hence decreasing the performance or destroying the high fitness of the individual. A low P_m, however, makes it hard for the GA to produce new individuals, or the algorithm may become a general random search and produce too many new individuals (Cao and Zhang 2010).

In GA, selection is a process of choosing parents for the next generation based on their scaled values from the fitness scaling function. The proposed RX-STPM and all of the RCGAs used in the comparative studies are based on tournament selection. Typical range for tournament selection (T_s) is 2 to 6. If T_s is too low, the convergence rate will diminish and, hence, the GA will consume more time seeking the optimal solution. On the contrary, if T_s is too high, the GA is more likely to converge prematurely, thus leading to an inaccurate (suboptimal) solution (Miller and Goldberg 1995). A heavy experimental study was conducted to determine the optimal parameter setting for RX-STPM. We identified the best combination values of P_c, P_m, and T_s for RX-STPM as 0.60, 0.003, and 2, respectively. The population size for RX-STPM is assumed to be 10 times the number of variables. The termination conditions and number of runs are identical for RX-STPM and the four RCGAs. Termination conditions rely on either the predetermined maximum number of generations (5,000) reached or known global minimum achieved, whichever precedes the other. The number of runs is 100 times with the same initial populations; each run is initiated using a different set of initial population. Table 2 shows the final parameter settings for all the algorithms. The proposed algorithm is implemented in MATLAB 2012 and the experiments are performed on a Core i7 Processor with 2.00-GHz speed and 8.00-GB RAM on a Windows 7platform.

TABLE 2 Parameter setting for all the algorithms

	GA name
GA parameter	WX-PM	LX-PM	WX-LLM	LX-LLM	GNMA	RX-STPM
Pop size	5N	10N	5N	5N	25	10N
Pc	0.7	0.65	0.75	0.7	0.6	0.6
Pm	0.08	0.05	0.04	0.06	0.01	0.003

N denotes number of variables.

RESULTS AND DISCUSSIONS

In this study, the performance is gauged based on the proposed algorithm for its reliability and efficiency. Reliability is measured by the number of successes in finding the global minimum or at least approaching it sufficiently, and efficiency is measured by the number of function evaluations required. Both of the following performance evaluation criteria are calculated based on successful runs only. A run is considered successful if a threshold of 99% of the global minimum function value is achieved.

1. Success rate (SR) = × 100.

2. Average number of function evaluations (AFE).

The proposed algorithm was tested against a hybrid GA namely Genetic Nelder-Mead Algorithm (GNMA; Ali and Hassanien 2013) and four RCGAs (WX-PM, LX-PM, WX-LLM, LX-LLM) from the literature. GNMA was based on a global search GA with a local search Nelder-Mead algorithm. The four RCGAs used Weibull crossover (WX), Laplace crossover (LX; Deep and Thakur 2007a), log logistic mutator (LLM; Deep and Katiyar 2012), and power mutator (PM; Deep and Thakur 2007b).

Table 3 depicts the computational results of the success rate and the average number of function evaluations for RX-STPM along with earlier published results applying GA methods. RX-STPM recorded a 100% success rate for the molecular potential energy function, E, for 20 and 40 degrees of freedom; that is, n = 20 and 40. As the chain size or dimension of the molecular model becomes larger, RX-STPM still manages to achieve the global minimum by having success rate of 98, 85, and 91% for n = 60, 80, and 100, respectively. As can be seen, the overall success rate for RX-STPM is significantly higher than the other five GA methods. This proves RX-STPM to be a more reliable (robust) algorithm. Through a comparative study on RX-STPM with the other four RCGAs (WX-PM, LX-PM, WX-LLM, LX-LLM), RX-STPM reportedly requires the least average number of function evaluation in all cases. What is even more astonishing is that RX-STPM also outperformed GNMA for the average number of function evaluation for n = 20 and 40. The results in Table 3 indicate the supremacy of RX-STPM over most of the algorithms in every facet. RX-STPM is hence proven to be a far more reliable and efficient RCGA.

TABLE 3 Computational results for potential molecular energy function for dimension n = 20 to 100

		GA name
Performance evaluation criteria	n	WX-PM^a	LX-PM^a	WX-LLM^a	LX-LLM^a	GNMA^b	RX-STPM
SR	20	100	100	100	100	90	100
AFE		15,574	23,257	28,969	14,586	8,982	8,654
SR	40	100	100	98	100	90	100
AFE		59,999	71,336	89,478	39,366	21,107	19,781
SR	60	93	91	89	95	100	98
AFE		175,865	280,131	225,008	105,892	35,391	97,832
SR	80	71	79	78	82	100	85
AFE		302,011	326,287	372,836	237,621	48,789	198,964
SR	100	89	88	82	86	90	91
AFE		369,376	379,998	443,786	320,146	56,399	298,460

^aDeep et al. (2011).

^bAli and Hassanien (2013).

CONCLUSION

To conclude, our adaptation of RCGA operators, RX and STPM, were successfully applied on a molecular potential energy function with known global minima. This energy function is a simplified molecular model. One of the limitations present is that the number of local minima grows exponentially with molecular size, hence making it difficult to find the global minima. We compared the computational results of RX-STPM with earlier published results using arithmetical crossover, Weibull crossover, Laplace crossover, log logistic mutation, power mutation, and the Nelder-Mead algorithm. Computational tests are executed with degrees of freedom or dimensions, n, ranging from 20 to 100. A heavy experimental study was conducted to set up the optimal parameter for the proposed algorithms to run the experiments. Numerical results presented in Table 3 shows that RX-STPM outperformed the majority of the algorithms, especially WX-PM, LX-PM, WX-LLM, LX-LLM, in all cases and hence redeems itself as a reliable (in terms of success rate) and efficient (in terms of average number of function evaluation) algorithm. The success of RX-STPM comes from detailed testing of the probability distribution parameters and relies on the proper design of RX as the key search operator. The offspring produced by RX inherit the favorable elements from both parents, leading to a steadfast convergence. The unique design of STPM, on the other hand, maintains the good chromosomes and diversity in the population. It is of utmost importance to have the correct combination of GA operators to ensure effective resolution to the problem. In this study, the GA was tested on a few operators and it has been found that the RX crossover is the most appropriate complement to the STPM mutator to reliably and efficiently resolve the problem. Therefore, future studies could look into the improvement of RX-STPM and subsequently applying it in other potential energy functions.