In a wide range of modern problems, it is required to estimate an influence of uncertainty of input parameters on uncertainty of an output value of a modeling function. In this contribution, we present algorithms for analyzing the sensitivity of a target function with respect to parameters in the problem of optimization of ReaxFF molecular-dynamic force field. In this particular case it allows one to effectively decrease the number of simultaneously optimized parameters. We compare the Sobol's global sensitivity indexes (SI) approach and the correlation analysis. Both methods are based on computations of the target function value on the set of pseudo- or quasi-randomly distributed points. The distribution derived is used for further computations of SI using Monte-Carlo technique and correlation coefficients. In the case of optimized ReaxFF force field one may spend up to several seconds to compute a value of the target function in a particular point. That is why it is important to perform calculations in parallel for multiple points. A parallel algorithm has been implemented in C++ using MPI. We compute Sobol's SI and coefficients of correlation of parameters variation and target function values variation while we optimize the force field for molecules and crystals of zinc hydroxide. We show that using of parameter set sorted by influence allows one to significantly increase convergence speed of the optimization algorithm and even completely exclude those parameters with relatively small influence.
In a wide range of modern problems, it is required to estimate the influence of the uncertainty of the input parameters on the uncertainty in the output value of the modeling function. In this paper there are presented algorithms of sensitivity analysis of an objective function to parameters in the problem of optimization of the molecular-dynamic force field ReaxFF. In this particular case, this effectively reduces the number of parameters simultaneously participating in the optimization. Two approaches are compared: one based on Sobol’s global sensitivity indexes (SI) by and the other one uses correlation analysis. A parallel program is presented, its scalability is studied, and the calculation results for a particular task are given.
A molecular-dynamic force field is characterized by a number of parameters: = ( 1, 2, … , ). These parameters are set before solving the Newton’s equations of motion and do not alter during the simulation. Their number ( ) may vary from 3 – 4 to several dozens. A procedure of a force field parameters search for a particular simulated system is referred to as a force field optimization. As a measure of optimality we use the parameter-dependent objective function (OF) ( 1, 2, … , ). The OF is determined by the deviation of any characteristics of the simulated system obtained using the methods of MD, from those obtained by more accurate methods. The process of optimization is to search the force field parameters that bring the OF the minimal value.
As the MD force field ( ⃗) we use ReaxFF (Reactive Force Field) that is able to simulate
chemical reactions [
The MGSA is able to simultaneously optimize relatively small number of parameters (the optimal variant is 4). However, optimizing ReaxFF for a particular system may need the number of parameters to be of several dozens. Therefore the parameters being searched are divided into groups, e. g. of 4 pieces, which participate in the optimization. The question arises of how to sort the parameters by groups, and also how to exclude parameters that almost do not affect the change in the OF. This task could be solved with a help of sensitivity analysis. In this paper we consider the correlation analysis and the method of Sobol’s sensitivity indices.
The objective function for searching MD force field parameters is a sum of terms depending on the parameters 1, 2, … , . The procedure of correlation analysis requires the following steps.
1. Generate a sample of (pseudo-)random sets of parameters 1, , 2, , … , , , 1 ≤ ≤ , that are uniformly distributed over the search domain of the OF.
2. For each set of parameters compute the OF value and record the values of particular terms. One will obtain a × matrix of OF terms’ values, where is the number of terms of the OF. In total, one has two matrices: the × matrix of parameters and the × matrix of OF terms . 3. Normalize matrices and : =
− ̅ √ ,
Here ̅ , ̅ are averages by columns, , ,
× matrix
, are dispersions by columns, is the sample size. = · of cross-correlations of parameters change and OF values change. The elements of the matrix, , are cross-correlation coefficients.
by the greatest absolute value of element in a row in descending order. Each row corresponds to a single parameter. monotonically, which is not always the case.
The correlation analysis will give good result only if the OF terms depend on the parameters =1 1<⋯< ( ) = 0 + ∑ ∑ 1,…, ( 1, … , ) =
in the unit cube . Then its ANOVA (Analysis of Variances) decomposition is
A more efficient technique for estimating the influence of variables on a function is the calculation of global sensitivity indices (SI) suggested by I. Sobol and A. Saltelli [3]. The method of sensitivity indices, in contrast to the correlation analysis, does not require the presence of a monotonic dependence of the function on the variables.
Let ( ) be a function of several variables
= ( 1, 2, … , ) defined and square integrable = 0 + ∑ ( ) + ∑ ( , ) + ⋯ + 1,2,…, ( 1, … , ), (1) 1
< The global sensitivity indices are the dispersion ratios
if 0 = ∫ ( ) , and ∫ 1,…, 0 = 0, when 1 ≤ ≤ . The inner sum in (1) is done by all the 1, … , , that satisfy inequalities 1 ≤ 1 < ⋯ < ≤ .
The decomposition (1) could be made with any complete orthonormal system of functions 0( ), 1( ), … , ( ), … , that includes the function 0( ) ≡ 1.
The quantities 1,…, = ∫ 21,…, ( 1, … , ) 1, … , are called the dispersions. Here and below the sign ∫ means integration from 0 to 1 by the corresponding variables. The quantity = ∫ 2( ) − 02 is called the total dispersion. It is also true that. = ∑ =1 ∑ 1<⋯< 1,…, .
1<⋯< = 1<⋯< ⁄ .
In application single-dimension sensitivity indices are most often used. With their help one is able to sort variables : the bigger is , the more influential is the variable .
Let us consider an arbitrary group of variables 1, … , , where 1 ≤ 1 < ⋯ < ≤ , 1 ≤ ≤ − 1. We shall denote them by a single letter = ( 1, … , ); and let be the aggregate of all the rest – variables. Thus, = ( , ). With let us denote the aggregate of indices ( 1, … , ). For the set let us introduce two types of global SI: 1 ≤ ≤ .
1 2 = ∑ 1,…, ,
tot = ∑ 1,…, .
In the summation is produced for all the groups of 1, … , , where all the ∈ . In tot the summation is produced for all the groups of 1, … , so that at least one index ∈ .
One is able to calculate the SI and tot by integrating using (pseudo/quasi) Monte Carlo technique. Let and ′ be the points of and let = ( , ), and ′ = ( ′, ′). In the paper [3] it is shown that = ∫ ( ) ( , ′)
′ − 02, tot = ∫[ ( ) − ( ′, )]2
′ = + 02 − ∫ ( ′) ( , ′) ′ .
In order to compute all the one-dimensional ( = ( )) SI and tot, in each sample of the Monte Carlo method one needs to use two n-dimensional random points and ′ and calculate the value of the function + 2 times: ( ), ( ′), and = ( 1′, … , ′ −1, , ′ +1, … , ′ ), when
One may select the points of parameters in a pseudo-random way, but it will not provide a
uniform cover of the whole unit cube . The alternative is to use LPτ sequences [
Let us compare a mesh of 16 points of this sequence with the simple cubic mesh for the square (Figure 1). The points of the simple mesh allow one to calculate function values only in 4 different values of each of the two variables. The mesh built by LPτ allows one to obtain function values in 16 different values of each variable for the same total number of points. Thus, LPτ provides more efficient distribution of points.
1 0
1 1 0
1
We implemented parallel programs for sensitivity analysis in C++ using MPI technology. The
use of sensitivity analysis in the problem of ReaxFF force field optimization has an important feature
that the time of calculation of a single OF value significantly exceeds the time of all other operations,
i. e. data exchange between processes of interaction with the file system. OF values are computed with
the help of LAMMPS molecular dynamic simulation library [
In Figure 2 the general scheme of the implemented programs is presented. In both cases (correlation analysis and sensitivity indices approach) the scheme is the same. At first one generates an array of N-dimensional points (ReaxFF parameters). The array is divided into equal parts between all the parallel processes and is sent between them to compute the OF values.
In the stage of OF values calculation the data exchange is not needed. Each process computes the values independently. The calculated OF values are gathered by the main process which computes correlation coefficients and Sobol’s SI.
Due to no data exchange in the main loop of the Figure 2. General scheme of programs programs good scalability is expected. This was verified on a cluster of the following configuration. 14 nodes: 2 × 4 core CPU Intel® Xeon® E5335 2.00 GHz, 16 GB RAM; 112 cores in total; OS: CentOS Linux 7 (Core); MPICH2 v1.4.1p1.
Hardware is provided by Resource Center Computer Center of St. Petersburg State University. We measured the dependence of acceleration of programs’ time of operation on the number of CPU cores in use for 1, 7, 14, 28, 56, and 112 pieces. The computations we done for 112 points of the search domain; when using all the 112 cores, we had one point per process. The average time of computation of a single OF value by a single process is 12.4 seconds. 7. Application
The programs presented have been applied to sensitivity analysis of the objective function to the parameters of ReaxFF force field in its optimization for Zn – O – H compounds.
The sensitivity indices (SI) and correlation coefficients (CC) calculation is performed for 58 parameters on 130000 points. The diagram on Figure 4 depicts relative values of SI and CC for 24 parameters of ReaxFF. For the unity we take the greatest SI in the first case, and the greatest CC in the second case. The parameters on the diagram are sorted by SI in descending order. On the chart for CC white hatching highlights those parameters that were not included in the first 24 in the correlation analysis. The number of such parameters is not large relative to the length of the whole list. One has to note that 5 parameters with the greatest CC are also located in the beginning of the list by SI. In general, the order of the parameters for the SI and CC is noticeably different.
AlphavdW_Zn-O
RvdW_Zn-O
Rsig_Zn-O EvdW_Zn-O
RvdW_Zn
ChiEEM_O AlphavdW_Zn
EtaEEM_O EtaEEM_Zn
Pbo1_Zn-O GammaEEM_O
Rsig_Zn
ChiEEM_Zn GammaEEM_Zn
Pbo2_Zn-O
RvdW_O
EvdW_Zn Pbe2_Zn-O
Pbe1_Zn-O AlphavdW_O
Rsig_O Pval2_O-Zn-O Pval1_Zn-O-Zn Pval4_O-Zn-O 1E-3 0.01 0.1 1 0.0 0.5 1.0
Using MGSA method we perform ReaxFF optimization by 24 parameters: first with the greatest CC, then with the greatest SI. The optimization is done in groups of 4 parameters. The sequence of the groups is defined by sorting the parameters by CC and SI, respectively. The loop of groups is repeated until the difference of OF values of the two consecutive iterations becomes less than a predefined tolerance. If the order of the groups in the optimization is defined by CC, the algorithm requires 2.5 loop passes to converge. When the order is defined by SIs, the algorithm needs
Sensitivity analysis is useful in the problem of optimizing of a large number of parameters of ReaxFF force field in groups. The approach allows one to exclude those parameters that have little influence on the objective function and also effectively sort them by influence.
The developed parallel program for Sobol sensitivity indices and correlation coefficients calculation is effectively scalable on a computational cluster.
To sort the parameters in ReaxFF optimization for zinc — oxygen — hydrogen compounds it is better to use Sobol sensitivity analysis than the correlation analysis, since the procedure converges faster in the first case. [3] Sobol’ I. M, Global’nye pokazateli chuvstvitelnosti dlya izucheniya nelineynykh matematicheskikh modeley [Global sensitivity indices for exploration of non-linear mathematical models] // Matematicheskoye modelirovanie [Mathematical modeling] — 2005. — Vol. 17 (9). — P. 43 – 52. (in Russian)