Strongin's multifactorial global search algorithm (MGSA) allows one to find an absolute minimum of a function of multiple variables on a mesh. In this contribution a parallel program is presented that implements the algorithm above applied to ReaxFF MD force field parameters search. In case of ReaxFF optimization, computation time of an objective function value significantly exceeds time of data exchange between parallel processes. One is able to speed up computation by obtaining not only one but several function values in various points simultaneously. Our software implements two levels of parallelism. To deal with function of multiple variables, one uses a scan for mapping a multi dimensional domain of definition of a function into a one-dimensional segment. To decrease the effect of losing information of multi-dimensional points proximity, scans are used. Function values in different mesh points are computed in parallel. This is the first level of parallelism. To define a mesh point of a next iteration, MGSA finds a subinterval with the most probable location of the minimum and computes an objective function value in a certain point of this subinterval. Function values are also calculated in parallel in ( − 1) subintervals with less probability. This is the second level of parallelism. Thus the two levels allow one to compute ∙ function values in parallel each iteration. In this contribution we research scalability of our MGSA implementation, namely, the dependence of the number of algorithm iterations and time it needs to converge on the number of CPU cores used, separately for each level of parallelism.
hardware.
Strongin’s multifactorial global search algorithm (MGSA) allows one to obtain an absolute minimum of a function of multiple variables on a mesh. This paper investigates scalability of our parallel MGSA implementation, developed to search for parameters of ReaxFF molecular-dynamic force field. The feature of the implementation is two-level parallelism. The aim of the work is to justify the choice of the optimal parameters of the algorithm for a specific calculation on the available Classical molecular dynamics (MD) method is to numerically solve Newton’s equations for every atom in a simulated system: ⃗ = −∇ ( ⃗ ), where ⃗ are impulses of atoms in the system, is time, ⃗ are coordinates of atoms in the system, and ∇ ( ⃗ ) is a gradient of a force field. The force field where
is the domain of parameters search having dimension of . depends on a set of parameters:
= ( 1, 2, … , ). These parameters are defined before the solution of Newton’s equations and do not alter during the simulation. The procedure of searching the force field parameters for a particular simulated system is referred to as optimization. The measure of optimality is the objective function (OF) that depend on the parameters: ( 1, 2, … , ). The OF is determined by the deviation of any characteristics of the simulated system obtained methods of MD, from those obtained by other, more accurate methods. The process of optimization is to search the force field parameters that bring the OF the minimal value: 1∗, 2∗, … , ∗ : ( 1∗, 2∗, … , ∗ ) = mi n ,
As the MD force field ( ⃗) we use ReaxFF (Reactive Force Field). It is the sum of several terms, each responsible for their own type of interaction, and has the general formula [1]: ReaxFF({ }, { }, { }, { }, {
}) = bond + lp + over + under + val + pen +
+ coa + tors + conj+ hbond + vdWaals + Coulomb.
(
As the measure of optimality we use the objective function (
are components of forces acting on every atom of the models, is the
Here are potential energies of models of the training-set (simple chemical compounds used
number of models in the training-set, is the number of atoms in the model , and are the weight
factors. Indices QC and ReaxFF mean that the corresponding characteristics are obtained with
quantum chemistry methods and MD methods, respectively. kReaxFF and kRαeiaxFF depend on the
parameters 1, 2, … , of ReaxFF force field. In our multifactorial algorithm in addition to the main
objective function (
To optimize the objective function we use Strongin global search algorithm [2] (GSA). It allows one to find an absolute minimum of a function on a segment and is bases on probability approach. Initially GSA is formulated for functions of one variable defined on a segment. The OF ( ) must satisfy the generalized Lipschitz condition on the entire search domain: | ( 1) − ( 2)| ≤ ( 1, 2), where 1 и 2 are any numbers from the search interval, is a constant, and is a metric in the space 1, ( 1, 2) = (| 1 − 2|), and the function has its own inverse −1. Initial step of GSA is to compute values of the OF at the boundaries of the search domain and, if necessary, in different arbitrary points inside it.
The general scheme of a single iteration is below. 1. Sort the points of previous iterations in ascending order of their coordinates: = 0 < ⋯ < < ⋯ < = . 2. Compute for every interval ( −1, ), 1 ≤ ≤ the characteristic ( ) that defines probability of finding the minimum on the given interval. 3. Sort ( ) in descending order: ( 1) > ⋯ > ( ). 4. For all the intervals having numbers 1, … , , < , compute OF in points + = ( ) ∈ ( −1, ), = ̅1̅̅,̅̅̅. Computations are done in parallel. Here ( ) is the position of the mathematical expectation of a minimum. 5. Check the stop condition: | 1 − 1−1| < .
The presented scheme allows one to execute GSA in parallel using processes. In the multivariate version of the algorithm a multivariable function is reduced to a function of one variable by means of scans of the Peano curve type. From here one can obtain another level of parallelism, applying the method rotating scans [3]. Each parallel process works with its own scan rotated relative to the main one by the angles ± in some pair of dimensions. In total one is able to do · ( − 1) of 2 such rotations for -dimensional domain of definition of a function. Thus, the maximal number of processes at this level of parallelism is max = · ( − 1) + 1. Each of them executed GSA and at every iteration communicates its result to all the rest of the processes. Working with several scans simultaneously compensates consequences of information loss about points proximity when a multidimensional domain is mapped to a segment.
In the current work we use the non-injective Peano type curve [2] when every point of a search domain of dimension may have up to 2 prototypes on a segment. In this variant GSA the minimum is being searched on a mesh. The number of mesh points on a segment is 2 ( +1) − 2 . Here m is the scan parameter that defines its complexity and the number of points of multidimensional mesh for a single dimension (2 + 1). More about the MGSA in the problem of ReaxFF optimization see in [4].
The program is written in C++ with use of MPI and has two levels of parallelism: the first one is parallel work with scans when mapping the search domain to a segment; the second one is parallel calculation of the OF value for greatest probability characteristics ( ) (see GSA general scheme). Thus, the program operates with a rectangular matrix of processes of size · .
Let us note that in case of ReaxFF optimization the computation time of the OF value significantly exceeds both the data exchange time and the time of all the rest procedures inside a single MGSA iteration.
On the first level of parallelism (Figure 1) each block corresponds to a single scan. All blocks are equivalent and perform the same sequence of routines, exchanging data between each other three times every iteration. ones mean exchange between 2nd-level processes inside a 1st-level block.
Notations:
— search domain points (type double[D]); — objective function values (type double);
— MGSA indices values (type int); — mesh points in segment [0, 1], to which the serach domain is mapped (type double).
The triplets { , , } are stored in a special structure which simultaneously is a linked list for all , common red-black tree for , and separate red-black trees for with the same .
This structure significantly reduces the time of a single iteration of the algorithm.
limited cases
The first level of parallelism
has maximal number of processes. The = 11 and
= 4 with averaging over 100 iterations were investigated. In Figure 3 the plots of average time of a single MGSA iteration dependencies on the number of used cores = ·
are shown for 7 different cases. is the number of cores on the 1st level of parallelism, and —
on the 2nd. The number of OF values calculated in a single iteration is equal to . So although the time of a single iteration increases with , the 17 n16 o ita15 r e itA14 S G13 M computation of an OF value accelerates. An OF value is computed for several dozens of molecules with the help of LAMMPS [5] MD simulation library. It takes from 10 to 13 seconds. Each process consumes less than 100 MB of RAM.
The chart shows two stages of growth. In the first stage, when 1 ≤ ≤ 4 (within 1 CPU), the increase of the duration of a single iteration is due to an increase in the calculation time of the value of the OF. This is a feature of the hardware and LAMMPS library. In the second stage, when > 4, the increase of duration is caused by overhead costs for data exchange between processes and interaction with the file system. Let us note that the duration increase does not depend on the distribution of cores by the levels of parallelism. For
< 11 the charts’ slopes are smaller, since overheads are declining with .
and the number of mesh nodes of the order of 106 the mesh on Let us plot the dependencies of MGSA acceleration by the number of its iterations needed for convergence on the number of parallel processes level and the number of scans of the 2nd = 1 for 1 ≤ ( ) (8), where ( ) is the number of MGSA
≤ 5. The acceleration ( ) = iterations needed for convergence.
The number of points of the segment is chosen so that for different it is of the same order — 106, and the algorithm converged in less than
For the case of = 4 and three fixed numbers of scans
= 1, 7, 13 let us compare the dependency of MGSA acceleration by the number of iterations and time on the number of cores 5). ( , ) =
on the 2nd level (Figure
(
MGSA
iterations and its operation time until convergence. In Figure 5 the blue lines indicate the dependencies by iterations ( ( )), the red lines — by time ( ( )). Linear approximation using least square method tells us of the constant scalability efficiency (∆ /∆ ) no matter how many cores are used. By the difference between the plots of 16 14 12 1 =10 N ,on 8 it rea 6 l
AI, Group 1 AI, Group 2 AT, Group 1
AT, Group 2
Analytical formulas of individual
terms of ReaxFF (
It is shown that the implementation of the multifactorial global search algorithm (MGSA) with meshes more than 106 nodes and executed on a cluster of 128 cores is scalable with constant efficiency up to 71 % when the number of optimized parameters > 2.
For the given hardware it is optimal to operate with groups of 4 ReaxFF force field parameters, since MGSA converges in less than an hour.
It is more efficient to distribute resources to the 2nd level of parallelism with only one scan.
The introduction of the 2nd level of parallelism (by probability characteristics of MGSA) and the increase in the number of cores to 128 reduces the calculation time by 30 – 50 times (compared to the previously used single-level MGSA implementation [6]).