CC BY
17
UDC 519.632, 544.139, 538.911
Fitting of 3D Reaction Surface of Three-Atomic System in the Natural Collision Coordinates
A. Gevorkyan*, T. Grigoryan*, M. Nalbandyan*, V. Hovhannisyan^, G. Nyman*
* Laboratory of Quantum Systems and Processes Institute for Informatics and Automation Problems P. Sevak 1, 0014 Yerevan, Armenia
^ Yerevan State University Department of Applied Mathematics 0049, Alex Manoogian 1, Yerevan, Armenia * Department of Chemistry Göteborg University SE-412 96, Goteborg, Sweden
The objective of this investigation is to prove the possibility of representation of the 3D interaction between particles in various reacting three-body systems by analytical function with a set of adjusting parameters in natural collision coordinates. Realization of this idea implies the procedure of fitting 3D numerical data by 3D analytical function or more precisely, calculation of adjusting parameters in mentioned analytical function. In the work LevenbergMarquardt algorithm is used on the basis of which the numerical method is developed. The possibility of implementation of 3D fitting with the big accuracy, on the example of reaction H + H2 is shown.
Key words and phrases: nonlinear optimization, parallel computation, reaction surface of atom-diatom collision.
1. Formulation of the Problem
Recently as was shown by authors [1], the three-body quantum reactive scattering problem in the curvilinear Natural Collision Coordinates (NCC) system may be reduced to the inelastic single-arrangement problem. Mathematically the problem consists of solution of a first-order ordinary differential equation's system. For numerical investigations of this system it is necessary to define a full interaction potential between particles in the NCC system (u, v, $). Here u is a coordinate along the curve of coordinate of reaction 0^/, which connects together (in) and (out) scattering asymptotic subspaces, v is a normal to the curve 0^/ coordinate, along which the full wavefunc-tion is localized, $ is a scattering angle. Recall that usually the reaction potential is constructed by means of ab — initio quantum calculations, after which this numerical data are used for fitting and reconstructing the analytical form of the interaction potential in terms of scaled Jacobi coordinates (qo,qi,$). Now the 3D analytical forms are well known for many reaction potentials V(qo,qi,&). For definition of reaction potential in the NCC system, the coordinate transformations (qo,qi,$) ^ (u,v,$) in the expression of potential V(qo,qi,$) are carried out. We can organize a one-to-one mapping between coordinate systems (q0, qi,$) ■ (u, v, $) in some subspace of intrinsic 3D configuration down form the curve 0. Following the work [2,3], we can define the curve 0, which connects (in) and (out) asymptotic channels in plane (qo,qi):
q0 = —^ + H + q+q, q- <qt < (1)
Wl — qeq)
where a and b are some constants. In eq. (1) q- and q+q are the mass-scaled equilibrium bond length of molecules in (in) and (out) channels correspondingly. Note that
Received 28th November, 2009.
This work was partially supported by ISTC-1451 project.
the variable ql is defined on a part of axis q1 £ (qeq, and can have only positive values.
Now we can write the inverse transformations from (q0,q1) to (u,v):
qo(u,v) = q^(u) — v sin <p(u), q1 (u, v) = q{(u) + v cos <p(u), (2)
where the angle ip(v) is determined from the requirement that the coordinate system (u, v) should be orthogonal:
dqc0
dql
= cot <p(u), lim cot <p(u) = {mAm,C'\ , (3)
I mBM J
where mA,ms and mc are masses of scattering particles, M is its sum.
The coordinate u describes the translational motion of three-body system between reactant and product channels and is changed along the curve 3if measured from an initial point u0. It in particularly can be determined by equation:
u = uo — c a + b(q'i — q~q). (4)
(Ql — Qeq)
Under the numerical modelling of the system of differential equations near the subspace borders a computation error appears. This problem we can solve by fitting numerical data in the NCC system. In the limits of the NCC system, the full interaction potential may be represented in the following form:
V(qo,qw&) = u(u,v,ti) = ^ Uj(u,v)Pj(costi), (5)
3
where Pj(x) is Legendre polynomial and m < Taking into account the orthogonal property of Legendre polynomials
+i
2
Pj (x)Pj' (x)dx = àjf
we can find that the full interaction potential between particles may be represented in the following view:
m
V (go,gi,â) = U (u,v,iï) = ^ Uj (u,v)Pj (cos ■&). (6)
3
Taking into account the orthogonal property of Legendre polynomials
+i
IP3 (x)Pj' (x)dx = ^^ Ô, -i
we can find:
n
u3 = (j + 1 /2) J u(u, v, ïï)Pj (cos ■&) sin tfdtf, (7)
It is obvious that if the V(q0, q1,r&) are known in the kind of an analytical function or are specified in the form of an numerical array, we can generate 2D numerical arrays (databases) and look for its analytical approximation. Based on our experience, it is
■LJ
convenient to present these functions as follows:
2
Uj (u,v) = ( ^ Ajk)(u)vk) e-2aj (u)v - Bj (u)e-13!(u)v, (8)
' k=0
where the functions Ajk\u) and Bj (u) provide a smooth transition from the bound state (AC) in the (in) channel to the bound state (AB) in the (out) channel, A, B and C describe the reacting particles. Analyzing the geometrical and topological features of the different energetic surfaces of reactions shows that we can use the following analytical form for these functions:
_ F(1) - F(0) F(2)^2 F-(u) = F(0) + -j_+__j_-__(9)
Cj (U) rj + 1+ e-2-nu + (e"fjU + e-7i«)2 ' (9)
where f(0) , f(1) , f(2) , and Tj are some adjusting parameters and
Fj(U) = (Ak(u), Bj(u), aj(u), fa(u)).
Thus now the main problem is the elaboration of numerical method for computation of adjusting parameters which would give it a possibility to carry out approximation of a numerical array with the given accuracy.
2. Chi-Square Fitting Method
There are various methods for fixing adjusting parameters in expressions (7)-(8), the relaxation method, the Newton method, and the modified Newton method [4,5] et etc. All the mentioned methods are based on the procedure of inverse Jacobian matrix computation with respect to the adjusting parameters. However the direct calculation of the Hessian by means of minimization methods is impossible. For a solution of considered problem in this paper we use nonlinear minimization method.
The minimization methods permit the iterative evaluation of a function and of its gradients. The second essential difference from all previous methods is that all the above methods are based on the linear search or linear minimization. The methods in issue are based on nonlinear procedures with the use of least-squares formalism. The calculation procedure of the gradient and Hessian in the Levenberg-Marquardt method is described in [6].
The j component of decomposition in the expression of full 3D reaction potential (3) under the consideration can be represented as: Uj = Uj(u, v, C), where C is the vector of adjusting parameters having the dimensions M. The mean-square deviation of the function %2 can be determined the following way:
N
x2(u, v) = £
=1
Uj(ui, Vi) - Uj(ui, Vi, C)
Oi
(in literature this method is called the "chi-square") ai-described standard deviation. The gradient of function is:
M = _2 v M-M^i, vi, C)] 9Uj(Ui,Vi, C)
dCk a.
=1
2
dCk
where Ujj = Uj(ui,vi).
The second derivative is computed as follows:
2
d2x*
^ 1 rUj(Ui,Vi, C) Uj(Ui,Vi, C)
1 Vj (u■i, Vi, C) Uj (U,i, Vi
¿¡[ dCk dCi
dCk dCi
Here we introduce some designations:
i ^2 « _i ^2
(10)
2 dCk 2 dCk dCi'
From (10) we obtain the relation:
M
at = ^2 ¡3k i SCi, i=i
where Ci satisfy a system of linear equations:
5Ci = const ■ ai, The second partial derivatives are computed by formula:
So, here all necessary formulas for a closed computation of adjusting parameters in 3D model potential (4)-(9) are presented.
The adjusting parameters in our previous article [7] about a reaction surface of the collinear collision, have been computed using the Levenberg-Marquardt's nonlinear optimization method [6].
3. Fitting of 3D Reaction Surface H + H2 System
The approximation of numerical data with an analytical surface is done by applying several numerical methods, mainly the Levenberg-Marquardt algorithm and Fourier transformation. The first step in this process is the analysis of the numerical data and finding an analytical function that visually matches the numerical data. In the case of collinear surface fig. 1 the following scenario has been followed.
Figure 1. Reaction Surface of Collinear Collision H + H2 by using Quantum-Chemical ab
Initio Calculations
First, a pre-calculated numerical array of reaction surface points is considered as a set of data points over (u, v) plain. The surface is sliced over the u axis, and the graphs of slices are analyzed. In the process of analyzing of the slices it was pointed out, that each slice can be approximated using the simple version of generalized Mors potential (see (8)):
Ugm (u,v) = A(u)
(—2a.(u)(v — vo(u))_ 2e-a(u)(v-vo(u))
(11)
For the fitting to be more accurate the following modifications on the numerical surface are applied. First, the surface is analyzed to find the optimal value of Az shift along the z axis (potential energy axis), for which z ^ 0 when v ^ _<x for any value of u. After that, on each slice, only negative values are considered for fitting (as it was shown by a series of numerical experiments, this gives the most accurate results). It's also worth mentioning, that the points of our interest lye in the "flute", and are localized near the "peak" in fig. 2 (that is, approximately, in the range u G [-7, +7]), that's why the numerical experiment aims to produce the best results in that area. The result of fitting the slices produces a 2D numerical array, that contains values for parameters A(u), and for each value of u.
A( u)
Figure 2. The Behavior of Modified Echart Function for a Reaction Surface of Collinear
Collision
Further, these sets of points are fitted (again using Levenberg-Marquardt algorithm) to the following modified Echart function:
J 2
A(u) = a +- _ a , + de-a(u-u0)2 + C1
K ' 1 I „''VI II-tl.n I
1 + e^(u-no) (el(u-uo) _ e—y(u-uo)j
2
where a, b, c, d, u0 and 7 are some adjusting parameters. Finding an approximated analytical representation of Mors potential parameters gives us an analytical approximation of reaction surface. This analytical representation of surface (fig. 3) gives a value of relative error of not more than 5 percent, at maximum in the region which important for a elementary atom-molecular processes (fig. 4). The analysis of input data shows that some inaccuracy in it can lead to this error. Hence, further reducing of relative error can be done by further refinement of input data in the first step. The non-collinear problem requires finding an analytical representation of (7) surface (using numerical values of that we already have), which is later used in (5) formula. As the numerical experience has proved, the values of m that are bigger than 8 give a variance in potential values less than 1 percent, the potential surface for non-collinear reaction can be calculated using above mentioned decomposition by Legendre polynomials with high degree of accuracy taking into account the first 8 Legendre polynomials.
For analytical approximation of numerically calculated Uj surface, again, the surface is sliced along its u axis, and the slices are first analyzed visually. During that
Figure 3. The Reaction Surface of Collinear Collision H + H2 Constructed by Analytical Formula (11) After Fixing of Adjusting Parameters
Figure 4. The Relative Error Between Figures 1 and 3
process, it was inspected, that some of the slices are "well behaved" while others follow an irregular pattern. To be able to deal with these irregularities, it was decided to approximate the slices using Fourier transformations. The Fourier polynomial of degree 40 showed a good approximation in the surface area of our interest (the peak and its surroundings).
References
1. Novel Algorithms for Quantum Simulation of 3D Atom-Diatom Reactive Scattering / A. S. Gevorkyan, G. G. Balint-Kurti, A. V. Bogdanov, G. Nyman // ICCSA. — 2007. — Vol. 3. — Pp. 1114-1123.
2. Marcus R. A. On the Analytical Mechanics of Chemical Reactions. Quantum Mechanics of Linear Collisions // J. Chem. Phys. — 1966. — Vol. 45. — P. 4493.
3. Marcus R. A. Analytical Mechanics of Chemical Reactions. III. Natural Collision Coordinates // J. Chem. Phys. — 1968. — Vol. 49. — P. 2610.
4. Самарский А. А. Численные методы. — Наука, 1997.
5. Бахвалов Н. С. Введение в численные методы. — М.: Наука, 1987.
6. William H. Numerical Recipes in C. — Cambridge University Press, 1992.
7. Gevorkyan A. S. et al. Modeling of the Potential Energy Surface of Regrouping Reaction in Collinear Three-Atom Collision System // Computational Science-ICCS 2003 / Ed. by A. V. Bogdanov, P. Sloot et al. — St-Petersburg, Russia: SpringerVerlage, 2003. — Pp. 567-577.
УДК 519.632, 544.139, 538.911 Подгонка трёхмерной поверхности реакции системы трёх атомов в естественных координатах столкновений
А. Геворкян*, Т. Григорян*, М. Налбандян*, В. Оганнисян^",
Г. Ньюман*
* Лаборатория квантовых систем и процессов Институт информатики и проблем автоматизации П. Севак 1, 0014 Ереван, Армения ^ Ереванский университет Факультет прикладной математики, 0049, Алекса Манукяна 1, Ереван, Армения
* Факультет химии Готеборгский университет БЕ-412 96, Готеборг, Швеция
Цель данного исследования — доказать возможность представления трёхмерного взаимодействия между частицами в различных реагирующих системах трех тел с помощью аналитической функции с набором параметров подгонки в естественных координатах столкновения. Реализация этой идеи подразумевает процедуру подгонки трёхмерных числовых данных с помощью 3-0 аналитической функции или более точно, вычисление параметров подгонки в данной аналитической функции. В работе используется алгоритм Левенберга—Марквардта, на основе которого был развит численный метод. Возможность реализации трёхмерной подгонки с большой точностью продемонстрирована на пример реакции Н + Н2.
Ключевые слова: нелинейная оптимизация, параллельные вычисления, поверхность реакции атом-молекулярного столкновения.
