Proceedings of the VIII International Conference "Distributed Computing and Grid-technologies in Science and
             Education" (GRID 2018), Dubna, Moscow region, Russia, September 10 - 14, 2018




   SENSITIVITY ANALYSIS IN A PROBLEM OF REAXFF
  MOLECULAR-DYNAMIC FORCE FIELD OPTIMIZATION
                              K. S. Shefov a, M. M. Stepanova b
         St. Petersburg State University, 7-9 University emb, St. Petersburg, 199034, Russia

                         E-mail: a k.s.shefov@gmail.com, b mstep@mms.nw.ru


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.

Keywords: sensitivity analysis, reactive force field, molecular dynamics, parameter optimization,
parallel algorithm, scalability

                                                       © 2018 Konstantin S. Shefov, Margarita M. Stepanova




                                                                                                        595
Proceedings of the VIII International Conference "Distributed Computing and Grid-technologies in Science and
             Education" (GRID 2018), Dubna, Moscow region, Russia, September 10 - 14, 2018




1. Introduction
         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.


2. The Problem
         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 [1]. To optimize the OF we use the multifactorial global search algorithm (MGSA)
suggested by Strongin [2].
         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.


3. Correlation analysis
         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 𝐹:
                                            𝑝𝑖𝑗 − 𝑝̅𝑖                  ̅𝑖
                                                                 𝑓𝑖𝑗 − 𝑓
                                      𝑝𝑖𝑗 =           ,   𝑓𝑖𝑗 =           .
                                             √𝑅𝜎𝑝,𝑖              √𝑅𝜎𝑓,𝑖

                  ̅𝑖 are averages by columns, 𝜎𝑝,𝑖 , 𝜎𝑓,𝑖 are dispersions by columns, 𝑅 is the sample size.
       Here 𝑝̅𝑖 , 𝑓
       4. Calculate the 𝑀 × 𝑁 matrix 𝐶 = 𝑃𝑇 · 𝐹 of cross-correlations of parameters change and
OF values change. The elements of the matrix, 𝑐𝑖𝑗 , are cross-correlation coefficients.
       5. Sort the rows of the matrix 𝐶 by the greatest absolute value of element in a row in
descending order. Each row corresponds to a single parameter.
       The correlation analysis will give good result only if the OF terms depend on the parameters
monotonically, which is not always the case.

                                                                                                        596
Proceedings of the VIII International Conference "Distributed Computing and Grid-technologies in Science and
             Education" (GRID 2018), Dubna, Moscow region, Russia, September 10 - 14, 2018




4. Sobol sensitivity analysis
         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
in the unit cube 𝐾𝑛. Then its ANOVA (Analysis of Variances) decomposition is
                          𝑛

         𝑓(𝑥) = 𝑓0 + ∑ ∑ 𝑓𝑖1 ,…,𝑖𝑠 (𝑥𝑖1 , … , 𝑥𝑖𝑠 ) =
                         𝑠=1 𝑖1 <⋯<𝑖𝑠

                              = 𝑓0 + ∑ 𝑓𝑖 (𝑥𝑖 ) + ∑ 𝑓𝑖𝑗 (𝑥𝑖 , 𝑥𝑗 ) + ⋯ + 𝑓1,2,…,𝑛(𝑥1 , … , 𝑥𝑛),        (1)
                                        𝑖      𝑖<𝑗
                                 1
         if 𝑓0 = ∫𝐾 𝑓(𝑥)𝑑𝑥, and ∫0 𝑓𝑖1 ,…,𝑖𝑠 𝑑𝑥𝑖𝑝 = 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 ,…,𝑖𝑠 .
         The global sensitivity indices are the dispersion ratios
                                               𝑆𝑖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 ,…,𝑖𝑠 ,      𝑆𝑦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 ,
                     1
             𝐷𝑦tot = ∫[𝑓(𝑥) − 𝑓(𝑦 ′ , 𝑧)]2𝑑𝑥𝑑𝑦 ′ = 𝐷 + 𝑓02 − ∫ 𝑓(𝑥 ′ )𝑓(𝑦, 𝑧 ′ )𝑑𝑥 ′𝑑𝑦.
                     2
       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
1 ≤ 𝑖 ≤ 𝑛.


5. LPτ Sequences
       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 [4]. These strictly
ordered sequences are based on multiple binary division of the domain by all the dimensions. They
allow one to construct a mesh on a hypercube so that the mesh nodes fill it as uniformly as possible.



                                                                                                                597
Proceedings of the VIII International Conference "Distributed Computing and Grid-technologies in Science and
             Education" (GRID 2018), Dubna, Moscow region, Russia, September 10 - 14, 2018



        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                                 1




                       0                          1         0                       1
                              Figure 1. Simple mesh (left) and LPτ mesh (right)


6. Implementation and Scalability
          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 [5]: it calculates values of energies and
forces for several dozens of molecules.
          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.




                                                                                                        598
Proceedings of the VIII International Conference "Distributed Computing and Grid-technologies in Science and
             Education" (GRID 2018), Dubna, Moscow region, Russia, September 10 - 14, 2018



         Figure 3 shows the measured dependence                             84

for both programs (CC and SI), since the values                             77

                                                                            70
of acceleration match with high accuracy. Linear
                                                                            63
approximation is done with the least squares                                56
technique. The efficiency of parallelism (the ratio




                                                             Acceleration
                                                                            49

of acceleration increase and cores number                                   42

increase) is (71.2±1.3) %.                                                  35

                                                                            28
         The time of operation of the serial part of
                                                                            21
the programs, i. e. the time of SI or CC                                    14
computation is less than 1 second for 112 points                             7

and is about 30 seconds for 130000 points, which                             0
                                                                                 0   7   14   21   28   35   42   49   56   63    70   77   84   91   98 105 112 119
is significantly smaller than the time of the                                                                               P
parallel part of the programs.                                              Figure 3. Accel. dependence on the cores number

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
                          Relative Sensitivity Indices                      Relative Correlation Coefficients
                       Figure 4. SI and CC for ReaxFF force field parameters for ZnO
        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


                                                                                                                                                                599
Proceedings of the VIII International Conference "Distributed Computing and Grid-technologies in Science and
             Education" (GRID 2018), Dubna, Moscow region, Russia, September 10 - 14, 2018



1.5 passes. Thus, sorting the parameters using the Sobol’s SI gives for this problem a better result than
using the correlation analysis.


8. Conclusion
         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.


References
[1] K. Nomura, R. K. Kalia, A. Nakano, P. Vashishta, J. L. Landa. A scalable parallel algorithm for
large-scale reactive force-field molecular dynamics simulation // Comp. Phys. Comm. — 2008. —
Vol. 178(2). — P. 73 – 87.
[2] Stepanova M.M., Shefov K.S., Slavyanov S.Yu. Multifactorial global search algorithm in the
problem of optimizing a reactive force field // Theoretical and Mathematical Physics. — 2016. — Vol.
187, Issue 1. — P. 603 – 617.
[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)
[4] Sobol’ I. M., Statnikov R. B. Vybor optimal’nykh parametrov v zadachakh so mnogimi kriteriyami
[Choise of optimal parameters in problems with multiple criteria]. — M.: «Drofa». 2006. — 175 p. (in
Russian)
[5] LAMMPS package: https://lammps.sandia.gov/ (accessed 30.09.2018)




                                                                                                        600