Iteration Scheme for Solving the System of Coupled Integro-Differential Equations for Excited and Ionized States of Molecular Systems

Investigation of the interaction of electromagnetic radiation with molecular systems provides most of the information on their structure and properties. Interpretation of experimental data is directly determined by the knowledge of the structure of energy levels and its change in the transition of these systems to an excited state. A key task of the methods for calculating the molecular orbitals of excited states is to accurately describe the emerging vacancies of the molecular core, leading to radial relaxation of the electron density. We propose an iterative scheme for solving a system of coupled integro-differential equations for obtaining molecular orbitals of electron configurations with excited/ionized deep and subvalent shells in a single-center representation. The numerical procedure of the iterative scheme is reduced to solving a boundary value problem based on a combination of the three-point difference scheme of Numerov and Thomas algorithm. To increase the rate of convergence of the computational procedure, an accurate account is taken of the behavior of the electron density near the nuclei of the molecular system. The realization of the algorithm of the computational scheme is considered on the example of a diatomic hydrogen fluoride molecule. The energy characteristics of the ground and ionized states of the molecule are estimated, and also the spatial distribution of the electron density is presented for the example of the σ-symmetry shell.


Introduction
Theoretical study of the absorption processes [1], photoionization [2], elastic [3], and inelastic [4] scattering of an X-ray photon by a molecular system is a complex problem that requires an exact description of the excited states with emerging vacancies in deep and subvalent electronic shells.Existing methods for calculating wave functions can be classified according to the form of their representation as multicenter (MC) [5] and single-center (SC) [6].
In the MC representation, molecular orbitals (MO) in different regions of space are described with respect to the centers associated with the nuclei of the molecular system.In this approach, the main limitation is related to the choice of the basic set of wave functions affecting the accuracy and amount of calculation [7,8].In the calculations, the construction of basic functions in the form of a superposition of Slater [9] or Gaussian [10] orbitals are widely used.The most well-known realization of the MC approach is the method of linear combination of atomic orbitals (MO LCAO) [11,12].Its equations were first written by Roothaan [5].In this case, the problem of minimizing the energy functional is reduced to finding the coefficients for the expansion of MO in basic atomic orbitals.A fairly complete overview on the use of the MO LCAO presented in [13], shows that the main drawback of the method is the principal difficulties in obtaining excited states of molecular systems caused by the problem of choosing the optimal basis set.
Another approach is to extend the methods for calculating the wave functions of atomic electrons to molecular systems that have been developed in the SC method of MO calculation [6,14,15].The SC approach is based on the MO representation in the form of an expansion in the basis of functions centered relative on one point of space.The application of the SC method has a number of advantages, especially when studying electron-photon interaction processes in many-electron systems [15][16][17].From the computational point of view, the simplicity of using SC in comparison with MC is due to the developed theory and methods of atomic calculations.The transition to molecular systems is determined by introducing summation over the orbital angular momenta, since one-electron atomic wave function is formally the SC component of the MO decomposition.This approach avoids the problems associated with the redundancy of the basis, or at least controls them, in contrast to the MC approach.In addition, the absence of multicenter integrals makes it possible to substantially simplify the procedure for performing calculations.
Historically, the first way to implement the SC method was to represent the partial harmonics of the SC expansion in the form of a linear combination of Slater-type orbitals and directly minimize the expressions for the total energy [6,11].In a different approach, described in detail in [15,16], another variant of the SC method is proposed, which reduces to solving the secular equation.In this case, both the discrete function and the continuous spectrum function are used as basis functions, which are a solution of the Hartree-Fock equations.This realization of the SC approach is not free from the problems of choosing optimal basis sets, as well as from the slow convergence of the SC expansion series of the MO.Moreover, it has a strong influence on the energy of the state of the partial SC harmonics of the MO decomposition with large orbital moment values under low contribution of these harmonics to the function norm.The last two drawbacks can also be attributed to the SC method in general.
The application of the variational principle for obtaining the wave functions of electrons in the case of an atom leads to the well-known system of integro-differential Hartree-Fock equations.In molecular SC systems, the representation allows to obtain a system of integro-differential equations connecting different partial harmonics of MO.This approach is most accurate due to the elimination of the problem of selecting basic sets for describing the MO.In works [17][18][19][20][21][22], this method was realized when calculating the wave function of an excited photoelectron of a discrete and continuous spectrum on the basis of a system of coupled integro-differential equations.Herewith, despite the application of the widely known method of solving the differential equations of Numerov [23], the developed approach has a number of drawbacks related to the use of the approximation of the frozen core, and also the restriction of the SC expansion of the MO to the harmonics.It makes an essential contribution to the expansion of the wave function, due to its slow convergence.In particular, the accuracy of describing vacancies in the excited and ionized states restricts the application of the approximations of "Z + 1" and "frozen core".In calculating the weakly excited states of molecular systems, the delta-SCF [24] method based on the DFT theory has been developed, requiring the search for a suitable exchange-correlation potential for an accurate description of the electron-electron interaction.At the same time, the application of the method is limited when solving a wide class of problems of elastic and inelastic X-ray photon scattering, as well as photoexcitation processes, which is due to the need to obtain highly-excited electronic states.At the same time, the slow convergence of the MO in the SC representation has a significant effect on the accuracy of the calculated values.According to the estimate made in [18], in order to describe the one-electron energies of MO molecules with ligands of the second period of the Mendeleyev periodic table with an accuracy of ~3%, one must take into account the partial harmonics with orbital angular momentum l ≥ 50 in the SC expansion of MO.The increase in the number of terms of the SC expansion leads to an increase in the system of differential equations, which causes more computational capacity giving an insignificant contribution of the obtained higher harmonics of the expansion to the observed quantities, but important for the energy characteristics of the systems [22].In [17,25,26] the problem of slow convergence of the SC of MO expansion was solved by "effective" accounting of the influence of partial harmonics with large values of the orbital angular momentum, by increasing the participation coefficient of the last from the accounted partial harmonics in the limited SC expansion when calculating the photoelectron MO.At the same time, such an "effective" technique in calculating individual spectroscopic quantities does not allow to correctly describe the electron-electron Coulomb interaction, including the photoabsorption cross-section.
The paper proposes an iterative scheme for calculating the excited and ionized states of molecular systems obtained outside the framework of the widely used approximations of "Z + 1" and "frozen core".The numerical procedure is based on solving a system of coupled integro-differential equations for all MOs included in the electronic configuration of the state.The higher harmonics of MO are represented as a composition of the wave functions of the electron shells of atoms expanded in a series of spherical harmonics.This ensures the required description accuracy of electron density near the nuclei of the molecular system.The proposed computational procedure is based on a combination of three-point difference Numerov scheme and Thomas algorithm.The theoretical basis and details of the iteration scheme are described.The efficiency of higher harmonics consideration in the MO decomposition is illustrated by the example of calculating the energy characteristics of the main electronic configuration of a hydrogen fluoride molecule with an increase in the ligand charge.The space expansion of electron density is also presented.The energy characteristics of the ionized states of the σ-symmetry molecule are estimated.

Basic Equation of Iteration Scheme
Integro-differential equations underlie the iterative scheme of obtaining radial parts of MO.They are the generalization of Hatree-Fock equation first recorded for an atom [27].The equations were obtained on the basis of minimizing the energy functional in the Born-Oppenheimer approximation [28] for the wave function on N-electron molecular system Ψ = (N!) −1/2 det φ nγ , represented as antisymmetrized product of one-electron MO φ nγ .After the variation procedure we obtain a set of coupled integro-differential Hartree-Fock equations systems each of which corresponds to a certain MO nγ of molecule electron configuration.In this case, the single-electron MO φ nγ , in the nonrelativistic SC representation, has the form [6]: where P nγ l (r) is the radial and Y µ l (ϑ, ϕ) is the angular part of l, the symmetry partial harmonic of MO with the fixed (for a linear molecule) value of the µ, the projections of angular momentum to the quantization axis OZ (molecule axis); n is the main quantum number; γ is the symmetry of MO; and r, ϑ, ϕ are the spherical coordinates.
Considering Equation (1) the system of coupled integro-differential equations for MO nγ has the form: where Ω nγ l,l (r) is the local potential, describing the electron interaction of MO nγ with the nuclei of the molecular system and the direct Coulomb interaction with the electrons of other MO; X nγ l (r) is the nonlocal potential, describing the exchange Coulomb interaction of MO electrons nγ with electrons of other MO.
Local potential Ω nγ l,l (r) has the following form: where Z 0 is the charge of molecular nucleus, located at the chosen origin of coordinates, W nγ l,l (r) is the potential of crystal field, describing the electron interaction of MO nγ with the nuclei of molecular system; V nγ l,l (r) is the potential of the direct Coulomb interaction of MO nγ electrons with the electrons of other MO; and ε nγ is the single-electron energy nγ of MO.
The expression for the potential of crystal field in Equation (3) has the form: where λ nγ µ is the number of nγ electrons of MO with µ, the projections of angular momentum to the quantization axis OZ (molecule axis); N L is the number of ligands in the molecule; Z α is the charge of the α-th ligand; 4 of 14 nteraction of МО electrons y of МО. as the form: rojections of angular momentum to the the molecule; is the charge of the α-th =max( , ) и ), , , are the f the spherical harmonic, the calculation rons with the electrons from oulomb interaction between iguration of a molecule; = tion of electrons of МО with (8) e Coulomb interaction between on configuration of a molecule; e expressions of the direct and artial harmonics of SC expansion over the μ-projection of angular he molecular system nuclei and O for calculating the system oefficients of the wave function are the coordinates of the α-th ligand; and C k lµ l µ is the matrix element of the spherical harmonic, the calculation methods of which are presented in [28].
The potential, describing the direct Coulomb interaction of electrons nγ with the electrons from other MO (Equation ( 3)) has the form: where α nγ n γ is the coefficient for the potential describing the direct Coulomb interaction between nγ and n γ MO depending on the quantum number of electron configuration of a molecule; Nonlocal potential describing the interchange Coulomb interaction of electrons nγ of MO with the electrons of other MO in Equation (2) has the form: where β nγ n γ is the coefficient for the potential describing interchange Coulomb interaction between nγ and n γ of MO depending on the quantum number of electron configuration of a molecule; The radial part of the electronic potential is determined in the expressions of the direct and interchange potentials.It is calculated on the basis of radial parts of partial harmonics of SC expansion of MO: In order to simplify Equations ( 6)-( 9), the summing is omitted over the µ-projection of angular momentum of the electrons participating in Coulomb interaction.

The Inclusion of Higher Spherical Harmonics in SC Expansion of MO
To accurately describe the electron density in the vicinity of the molecular system nuclei and solve the problem of slow convergence of the functional series (1), MO nγ for calculating the system of coupled differential equations has the following form: where A nγµ is the normalization factor; Λ nγL µα is the contribution coefficients of the wave function ψ nγL µl 0 (r) with the main quantum number n and symmetry γ of the L-th ligand of the molecular system in the form: The functions ψ nγL µl 0 (r) are the wave functions of the electron shells of atoms entering the molecular system as ligands (atoms are not at the origin) expanded in a series of spherical harmonics.
In this case, the partial harmonics of MO in Equation ( 11) with the values of the orbital angular momentum l > l 0 are obtained by solving the system of Equation ( 2).The terms of the expansion Equation (11) with orbital angular momenta l > l 0 are obtained on the basis of the linear combination ψ nγL µl 0 (r) with the contribution coefficients Λ nγL µα .Such a representation of MO is due to the fact that the shape of the partial spherical harmonics with large values of the orbital angular momentum is determined to within a coefficient by the ligand potential of the molecular system.
The choice of the value of l 0 in Equation ( 11) is a result of a numerical experiment and depends on the nature of the molecular system and the number of partial harmonics involved in the formation of the chemical bond of the compound.
To determine the contribution coefficients Λ nγL µα in Equation (11) the partial harmonics are linked in using the equation: That is, radial part of the l 0 -th partial harmonic of the nγ molecular orbital is represented as a superposition of radial partial harmonics of ligand functions with weighting coefficients Λ nγL µα .Equation ( 13) allows compiling a system of algebraic equations for determining the contribution coefficients Λ nγL µα of ligand functions with γ symmetry in the form of: Thus, the determined MO is normalized to unity.

The Expansion of the Wave Functions of the Electron Shells of an Atom in a Series of Spherical Harmonics
The procedure of obtaining the ligand functions ψ nγL µl 0 (r) used in Section 2.2 is based on the recalculation of the wave functions of the electron shells of atoms P A nl (r), entering the molecular system in a new unified system of coordinates.
In the transition to a uniform coordinate system of the wave function of an electron shell related to a certain atom is represented in the form of an expansion in a series of spherical harmonics [29]: where Φ A nlm (r, ϑ, ϕ) is the wave function of nl electron shell in the atomic system of coordinates; r, ϑ, ϕ are the spherical coordinates of the atomic system; Γ nlµ LM (R) is the coefficient of expansion in spherical harmonics in the molecular spherical coordinates system; L is the orbital angular momentum of expansion of the wave function in the molecular coordinates system; M is the projection of the orbital angular momentum in the molecular coordinates system; and R, Ξ, φ are the spherical coordinates in the molecular system.
The coefficient of expansion of the series in Equation ( 15) can be represented in the form [29]: where ( , , ) is the wave function of electron , , are the spherical coordinates of the atomic system spherical harmonics in the molecular spherical coord momentum of expansion of the wave function in the projection of the orbital angular momentum in the molecu spherical coordinates in the molecular system.
The coefficient of expansion of the series in Equation Υ are the origin coordinates of the molecular system in coefficients of the theta-function representation in the trig to a certain atom is represented in the form of an expansion in a series of where ( , , ) is the wave function of electron shell in the ato , , are the spherical coordinates of the atomic system; Γ ( ) is the spherical harmonics in the molecular spherical coordinates system; momentum of expansion of the wave function in the molecular coo projection of the orbital angular momentum in the molecular coordinates spherical coordinates in the molecular system.
The coefficient of expansion of the series in Equation ( 15 where R M , In the transition to a uniform coordinate system of the wave function of an electron shell related to a certain atom is represented in the form of an expansion in a series of spherical harmonics [29]: where ( , , ) is the wave function of electron shell in the atomic system of coordinates; , , are the spherical coordinates of the atomic system; Γ ( ) is the coefficient of expansion in spherical harmonics in the molecular spherical coordinates system; is the orbital angular momentum of expansion of the wave function in the molecular coordinates system; is the projection of the orbital angular momentum in the molecular coordinates system; and , Ξ, are the spherical coordinates in the molecular system.
The coefficient of expansion of the series in Equation ( 15) can be represented in the form [29]: Θ ( ) = sin( ) ∑ ( ) Υ are the origin coordinates of the molecular system in the atomic system, are the coefficients of the theta-function representation in the trigonometric form.
, F are the origin coordinates of the molecular system in the atomic system, T lm s are the coefficients of the theta-function representation in the trigonometric form.
Integration in Equation ( 16) is performed on the basis of Equation ( 17) and also relations in Equations ( 18) and ( 19), which connect the coordinates in different systems.It enables the expansion of the wave function of the electron shells of atoms with respect to the chosen molecular center.

Scheme for Numerical Solution for the System of Equations for MO
The numerical procedure of MO finding is reduced to solving a system of coupled differential Equation (2) for the radial parts of partial harmonics P nγ l (r) of Equation ( 1).The Equation (2) in general can be represented in a matrix form: can be obtained «on the left» K L n , M L n and «on the right» K R n , M R n from the crosslinking point Z of the solutions.
Boundary conditions are used in order to obtain take-off ratios of recurrence relationships (Equation ( 24)).The conditions are imposed on the radial parts of partial harmonics of MO and are determined by their exponential attenuation if r → 0 ( r → ∞ ).Application of such conditions to Equation (24)

results in matrix equations
To obtain the starting coefficients of the recurrence relations (Equation ( 24)), the boundary conditions are applied on the radial parts of the MO partial harmonics, which are due to their exponential attenuation as r → 0 (r → ∞).Applying these conditions to (24) leads to matrix equations in the form: The iterative solution of Equations ( 25) and ( 26) allow restoring the coefficient values K L n , M L n , and K R n , M R n at all points of the coordinate grid.The solution of Equation ( 21) must satisfy the smoothness requirement corresponding to the equality of the first-order derivatives of a function obtained by sweeping the values "on the left" and "on the right" at the crosslinking point Z: This condition leads to a system of linear inhomogeneous equations with respect to the value of the partial harmonic at the crosslinking point P Z : The algorithm for implementing the numerical procedure of solving the system of Equation ( 21) consists in determining the eigenvalues ε nγ and corresponding MO of the γ-type symmetry.
The solution of the homogeneous system of equations obtained from Equation ( 21) by zeroing out the nonlocal summand X = 0 is performed at the stage of initialization of the iterative procedure using the initial approximation (see Section 3.2) for MO.This is, in fact, the solution of the Hartree equation, which does not take into account the exchange terms.When the numerical procedure is implemented, the homogeneous Equation ( 28) is solved, which right-hand side is equal to zero.Moreover, the existence of non-zero solutions is determined by the criterion: The selection of solutions on the energy scale until reaching the corresponding MO nγ.For example, the σ symmetry states for diatomic molecules can be selected by the principal quantum number.In this case, the search for highly excited states of the discrete spectrum involves the sequential determination of all lower one-electron energies.The implemented procedure for finding solutions accounts the obtained spectrum of state energies and uses them as "reference points" on the one-electron energy scale ε nγ of MO within the range of [−E, 0] at subsequent iterations of the computational procedure.This approach ensures time saving for calculations and increases the efficiency of the algorithm for finding solutions of Equation (21).
The final set of radial parts of the partial harmonics P nγ l (r) c l ≤ l 0 is found for MO when solving Equation (21).Following the procedure given in Section 2.2, the partial harmonics with orbital angular momenta l > l 0 are determined by the superposition of ligand functions.The participation coefficients Λ nγL µα of the ligand wave functions Ψ nγL µl 0 (r) are determined by the solution of the system of algebraic Equation ( 14).Here, the functions Ψ nγL µl 0 (r) are preformed according to the numerical procedure given in Section 2.3 on the basis of the wave functions of electron shells of atoms, obtained by solving the Hartree-Fock equations.
Equation ( 21) has a unique solution for any value of one-electron energy ε nγ when accounting nonlocal terms X.Therefore, the selection of solutions is based on the overlap integral, which is obtained according to the expression: The criterion for selecting solutions for MO at the intermediate iteration step of the computational procedure is the correspondence of the overlap integral of the wave function (Equation ( 30)) to the value obtained at the first iteration of the matching process.

Starting Wave Function MO and Stopping Criteria
Initial approximation for the radial parts of MO partial harmonics of molecule electronic configuration is formed on the basis of the procedure described in Section 2.2.The wave functions of the electron shells of the atom located in the origin of the molecular coordinate system are used for the radial parts of MO partial harmonics with l ≤ l 0 .Ligand functions ψ nγL µl 0 (r) are calculated for the formation of MO partial harmonics with l > l 0 using the procedure described in Section 2.3.
The calculation is based of preliminary obtained wave functions for the electron shells of neutral atoms in the Hartree-Fock approximation.
The iterative scheme for calculating the MO includes two basic computational cycles.The cycle of recalculation of MO wave functions is being performed till the defined accuracy of the found solution of Φ i nγ on the i-th step is reached in comparison with the i-1-th of the process.Here, an accuracy estimation of the matching process starting is carried out on each step starting from the second one according to the formula: The accuracy ∆ε τ of determining ε nγ is used for recalculating the wave function of MO in the new potential when solving Equation (28).

Summary of the Iterative Scheme
A general iterative scheme for implementing a numerical procedure of calculating the excited and ionized states of molecular systems is represented in Scheme 1.
The accuracy Δ of determining is used for recalculating the wave function of MO in the new potential when solving Equation (28).

Summary of the Iterative Scheme
A general iterative scheme for implementing a numerical procedure of calculating the excited and ionized states of molecular systems is represented in Scheme 1.The construction of the MO wave functions is performed on the preparatory stage according to the Section 2.2.In particular, the ligands functions Ψ nγL µl 0 (r) are formed on the basis of the wave functions of the electron shells of atoms, expanded into series of spherical harmonics (see Section 2.3).The local potential Ω nγ l,l (r) is calculated from Equation (3) for the orbital angular momenta l ≤ l 0 , which includes the potential of crystal field W nγ l,l (r), and direct Coulomb interaction V nγ l,l (r) (see Section 2.1).At the first iteration, a homogeneous system of Equation ( 21) is solved for every MO with X = 0. One-electron energies ε nγ are searched for in accordance with condition (29) within the interval [−E, 0].There is a direct correspondence between the solution number and the principal quantum number of MO with the chosen symmetry type.The found one-electron energy ε nγ provides the radial parts P nγ l (r) of MO with the use of Equation (28).Wave functions of MO with higher harmonics l > l 0 are formed via Equation (11) on the basis of the partial harmonics MO nγ with l ≤ l 0 (see Section 2.2).
Participation coefficient Λ nγL µα of ligand function Ψ nγL µl 0 (r) is determined by the solution of the system of Equation ( 14).The MO is normalized to unity for determining A nγµ in Equation (11).The resulting normalization integrals { Φ nγ |Φ nγ } are used as selection criteria of the solution on the subsequent steps of the iteration scheme.The found single-electron energies {ε nγ } allow to determine intervals on the energy scale for searching the energy values of the new iteration, thereby increasing the efficiency of the computational scheme.
Hereafter, an iterative procedure is realized aiming to minimize the energy of molecular system.The criterion of achieving a positive result is a decrease in the deviation of MO wave function calculated in accordance with Equation (31) below a specified value.The cycle of the numerical procedure on the i-th iteration involves recalculation of the electronic potentials V nγ l,l (r) and X nγ l (r) using the MO wave functions Φ i−1 nγ obtained on the i-1-st iteration step.The solution of Equation ( 21) is searched.On the first step, the value nγ is set in accordance with the criterion of Equation (29).The next step is searching the value of nγ in the energy range nγ ± ∆ for which the overlap integral Φ nγ |Φ nγ i of the given solution Φ i nγ (r) coincides with the normalization integral determined at the first iteration of the process.

Results and Discussions
The program code realizing the above iterative scheme, as well as the results of calculating the ground and excited states are available at the electronic address [33].
The hydrogen fluoride (HF) molecule was chosen for testing the numerical procedure due to the following circumstances: A relatively small number of molecule shells allows the exploration of the possibility of the developed iterative scheme and the procedure for taking into account the higher harmonics in the SC expansion of MO when they are implemented on the ECM with limited CPU time; 2.
The spectrum of free states and the structure of the wave functions have been thoroughly investigated [34].
In the calculations, the reference system associated with the F atom is chosen.To systemize the one-electron states, the quantum numbers characteristic for molecules of the point symmetry group C ∞ν are used.The MO calculations for the HF molecule are performed for the equilibrium distance between the nuclei of hydrogen atoms and fluorine R 0 = 1.733 a.u.[35].
Comparative calculations of the energy characteristics (one-electron MO energies and total energy) of the ground state 1σ 2 2σ 2 3σ 2 1π 4 of the HF molecule are given in Table 1.In section (a) of Table 1, the calculation was performed via the MC method of MO LCAO using the General Atomic and Molecular Electronic Structure System [12].The calculation results in section (b) of Table 1 were obtained in this work using the MO expanded into a series (1) and limited l ≤ 8.They show sufficient accuracy which correlate with the calculation methods.The use of higher harmonics of the ligand wave functions for an accurate description of MO in the molecular nuclei region (see Section 2.2) gives an additional contribution to the total energy and regularly changes the energy level of the 1σ-state of the HF molecule.The changes in the energy quantities observed in section (c) of Table 1 are primarily due to the contribution of the 1s state of the H atom, which is a ligand in this case.Construction of MO is determined by the selection of l 0 which includes the ligand functions in the MO structure.Figures 1 and 2 illustrate the procedure of determining l 0 with respect to the calculation (section (c), Table 1).Figure 1 displays SC expansion of MO 2σ of the ground state of HF molecule (calculation is performed for the SC series (1) with l ≤ 60). Figure 2 displays the SC expansion of the wave function of H atom 1s state (the solution of the Hratree-Fock equations), decomposed in the molecular coordinate system (see Section 2.3).
harmonics of the ligand wave functions for an accurate description of MO in the molecular nuclei region (see Section 2.2) gives an additional contribution to the total energy and regularly changes the energy level of the 1σ-state of the HF molecule.The changes in the energy quantities observed in section (c) of Table 1 are primarily due to the contribution of the 1s state of the H atom, which is a ligand in this case.Construction of MO is determined by the selection of which includes the ligand functions in the MO structure.Figures 1 and 2 illustrate the procedure of determining with respect to the calculation (section (c), Table 1).Figure 1 displays SC expansion of MO 2σ of the ground state of HF molecule (calculation is performed for the SC series (1) with ≤60). Figure 2 displays the SC expansion of the wave function of H atom 1s state (the solution of the Hratree-Fock equations), decomposed in the molecular coordinate system (see Section 2.3).The example of MO 2σ shows (see Figure 1b) that the form of radial parts of the partial harmonics of SC expansion of MO starting with =8 ( is the partial harmonic denoted in the figure) is in good agreement with the form of radial parts of the partial harmonics of SC expansion of the wave function of H atom 1s state (see Figure 2b).It can be concluded that the influence of the ligand field with the charge = 1 on radial parts of the partial harmonics of SC expansion of 2σ-МО with orbital moments ≥ 8 is decisive.Thus, in accordance with the procedure proposed above (see Section 2.2) the crosslinking of MO (1) with ligand functions ( ) is done starting from = 8.
A significant influence of the higher harmonics of the SC expansion on the energy characteristics of the molecular system is observed with the growth of the ligand charge.To estimate this effect, the calculation of the total energy of the electron configuration 1 2 3 1 is performed The example of MO 2σ shows (see Figure 1b) that the form of radial parts of the partial harmonics of SC expansion of MO starting with l = 8 (l is the partial harmonic denoted in the figure) is in good agreement with the form of radial parts of the partial harmonics of SC expansion of the wave function of H atom 1s state (see Figure 2b).It can be concluded that the influence of the ligand field with the charge Z L = 1 on radial parts of the partial harmonics of SC expansion of 2σ-MO with orbital moments l ≥ 8 is decisive.Thus, in accordance with the procedure proposed above (see Section 2.2) the crosslinking of MO (1) with ligand functions ψ nγL µl 0 (r) is done starting from l 0 = 8.A significant influence of the higher harmonics of the SC expansion on the energy characteristics of the molecular system is observed with the growth of the ligand Z L charge.To estimate this effect, the calculation of the total energy of the electron configuration 1σ 2 2σ 2 3σ 2 1π 4 is performed depending on the charge Z L which varied within the range of Z L = 2 ÷ 5 (instead of case Z L = 1).The results obtained for the total energies are given in Table 2.
The values of the total energies on the basis of MO with SC expansion limited to l ≤ 4 are given in column (a) of Table 2.The calculation results accounting the ligand functions with expansion up to l = 60 are given in the column (b) of Table 2.It is obvious that higher harmonics of the SC expansion of MO begin to significantly affect the energy characteristics of the molecular system.The calculation results of the MO electron state 1σ 2 2σ 2 3σ 2 1π 4 show that this iterative scheme solves the main problem of the SC methods.That is the description of MO behavior in the area of ligand nuclei.
It is also evident that with increasing the ligand charge, the contribution of the higher harmonics of MO, taken into account by this SC approach, increases substantially.Thus, the represented iterative scheme of the SC calculation of MO becomes effective when studying molecules with nonhydrogen ligands.The results of the research show that taking into account the higher harmonics of the expansion of the wave functions of the ligands electron shells is a necessary condition in performing calculations in the SC representation for obtaining accurate energy characteristics of molecular systems.
The possibilities of obtaining the excited states represented by the iterative scheme are considered on the example of ionization of σ symmetry shells of HF molecule ground state.
Table 3 displays the results of accounting the ligand functions with MO expansion up to l = 60.Energy characteristics of the electron configurations: 1σ ≡ 1σ 1 2σ 2 3σ 2 1π 4 ; 2σ ≡ 1σ 2 2σ 1 3σ 2 1π 4 ; 3σ ≡ 1σ 2 2σ 2 3σ 1 1π 4 are given as the results.Obtaining the ionized states in the «Z + 1» approximation of MC approach requires a special charge increase of the chosen core of the molecular system, which corresponds to the creation of a vacancy.Such an artificial technique is used in connection with the absence of the possibility of directional ionization of the electron shell.Unified approach to the calculation of any states is in the proposed iterative scheme.It ensures the possibility of ionization of any electronic shell.

Conclusions
The iterative scheme proposed in the paper is intended for calculating the wave functions of excited and ionized states of molecular systems.The procedures is based on a numerical solution of a system of coupled integro-differential equations for the radial parts of MO partial harmonics that determine the formation of the chemical bond of the many-electron compound and take into account the higher harmonics of the SC series.The harmonics describe the electron density strongly connected with the molecular system cores.The realized calculation scheme lacks the drawbacks of the existing SC calculation methods associated with the slow convergence of the expansion.The problem of calculating the excited states was overcome.It was connected with the use of "artificial methods" of creating vacancies in electron shells (approximation «Z + 1»).The proposed computation scheme enables performing calculations for the excited states with any number of unfilled electron shells.The calculations carried out for the HF molecule showed good agreement with the results of existing MC approbated software systems.
Ξ) cos(Υ) + sin(Ξ)sin(Υ)cos( − Υ are the origin coordinates of the molecular system in the atomic sys coefficients of the theta-function representation in the trigonometric form ) + sin(Ξ) sin( to a certain atom is represented in the form of an expans the spherical coordinates of the atomic system spherical harmonics in the molecular spherical coor momentum of expansion of the wave function in the projection of the orbital angular momentum in the molec spherical coordinates in the molecular system.The coefficient of expansion of the series in Equatio Γ Ξ) cos(Υ) + sin(Ξ Υ are the origin coordinates of the molecular system i coefficients of the theta-function representation in the tri

Scheme 1 .Scheme 1 .
Scheme 1. Iteration scheme of wave functions of MO of molecular systems electron configurations.Scheme 1. Iteration scheme of wave functions of MO of molecular systems electron configurations.

Figure 1 .
Figure 1.Radial parts of partial harmonics of SC expansion of 2σ-MO of the ground HF molecule state.Calculation of the ground state is performed for the equilibrium distance ( = 1.733 a.u.[35]) between the atom nuclei of H and F: (a) SC expansion of MO into a series (1) с ≤60; and (b) SC expansion of MO into a series (1) с ≥ 8.

Figure 1 .
Figure 1.Radial parts of partial harmonics of SC expansion of 2σ-MO of the ground HF molecule state.Calculation of the ground state is performed for the equilibrium distance (R 0 = 1.733 a.u.[35]) between the atom nuclei of H and F: (a) SC expansion of MO into a series (1) c l ≤ 60; and (b) SC expansion of MO into a series (1) c l ≥ 8. Algorithms 2018, 11, 1 11 of 14

Figure 2 .
Figure 2. Radial parts of partial harmonics of SC series of H atom 1s wave function found via a solution of Hatree-Fock equations in the molecular axis of reference: (a) SC expansion into series (1) with ≤60; and (b) SC expansion into series (1) with ≥ 8.

Figure 2 .
Figure 2. Radial parts of partial harmonics of SC series of H atom 1s wave function found via a solution of Hatree-Fock equations in the molecular axis of reference: (a) SC expansion into series (1) with l ≤ 60; and (b) SC expansion into series (1) with l ≥ 8.

Table 1 .
One-electron MO energy and the total energy of the HF molecule (in atomic units).

Table 3 .
One-electron energy MO and the total energy of the ionized states of the HF molecule and (in atomic units).