Numerical Research of Materials Crystal Lattice Parameters Based on Rare-Earth Metals

Geometrical parameters (coordinates and angles) of CeO2 crystal lattice by molecular dynamics method are calculated. Calculated parameters of crystal lattice are applied for definition the energy band structure via Hartree-Fock method in an approximation to CO LCAO (crystal orbitals as linear combination of atomic orbitals) and using the model of cyclic cluster. Calculated minimum energy band p-d is within the value range of experimental data. Valence band maximum is 4.2 while minimum energy band p-d width is 2.8 eV. Quantum-chemical calculations are accelerated by Schwarz inequality and direct inversion method in iterative subspace. The obtained mathematical model is implemented into software package for calculating material properties.


Introduction
In view of the fact that the computer capacity is being improved, the calculation of crystal physical properties is becoming a highly topical problem at present.This, in itself, generates interest in the possible reduction of experiments, as such calculations could present the preliminary characteristics of investigated sample and finalize the relevant experimental proof for these characteristics.
In this case such updated methods as Hartree-Fock method and Ab Initio method, molecular dynamics method (MDM) and hybrid methods are applied [1,2].These methods having different mathematical framework are focused on various research areas.However, it should be noted that Ab Initio method exhibits high calculation accuracy, while MDM-rapid rate.
It is known that about one-third of computer time spent supercomputers to solve problems in the area of atomistic modeling [3].In this way, currently known mathematical models that describe atomic systems from first principles are resource-intensive.Therein one has the problem of finding new materials -high need of mathematical models in computational resources severely limits the possibility of modeling the atomic structures, in particular -possible number of simulated atoms.
To solve this problem has been developed the software package for materials modeling that supports the methods of molecular dynamics and quantumchemistry models.
Presently for crystals there are a restricted number of software packages and mathematical models which can be used for description their structures from the first principles.With respect to quantum chemistry, the calculation of molecule reaction properties and electronic crystal properties is more challenging than the calculation of crystal physical properties.In this case, there is a set of mathematical models suitable for posed tasks.
By applying the band theory for solids the following energy bands were defined: maximum and minimum of conduction bands, as well as valence and energy bands.These properties were identified in band diagram (Fig. 1).
The aim of this paper is to demonstrate the operability of software package for materials modeling both in terms of molecular dynamics and Ab Initio simulations by the example of CeO 2 crystal.In the present paper geometric coordinates of primitive cells and perfect CeO 2 crystals are examined by molecular dynamics method.The coordinates are calculated on the basis of the energy band structure via Hartree-Fock method in an approximation to CO LCAO (crystal orbitals as linear combination of atomic orbitals) by applying cyclic -cluster model.

Mathematical model
One of the most commonly-accepted mathematical models in first-principles calculations of crystal properties is Hartree-Fock method (or Hartree-Fock method in an approximation to CO LCAO -crystal orbitals as linear combination of atomic orbitals).
Regardless of several attempts to accelerate the calculations by Hartree-Fock method [4][5][6], there are still reasons according to which the calculations of real crystals are difficult.It is a well-known fact that even 1 mm³ volume crystal contains an immense amount of atoms, which excludes the possible calculations of realworld samples.This is due to a lack of computing resources for short − time interval calculations [3].
There are several approaches in solving such tasks.For example, a primitive cell is closed by pseudo-atoms on its boundaries and the derived system is calculated.This would be an optimal solution only if the basic criterion involves the shortest calculation time.However, the calculation accuracy, in this case, remains low.
Another approach is based on cyclic-primitive cell and / or cyclic -cluster method.The calculation difficulty of crystals involves the formal infinite number of electrons in the system, while density matrix construction is possible only in case of the finite number of electrons This problem could be solved by imposition of cyclic boundary conditions [7].It is this approach that is used in the following paper.
Mathematical calculations were taken from famous LCAO-CO method (crystal orbital as a linear combination of atomic orbitals) that is an analogue of the Hartree-Fock-Roothaan method for crystals.On the basis of method results it was possible to connect atom properties which form a crystal with the crystal itself.
Calculation via Hartree-Fock method involves the following sequence: The first step includes calculating electron kinetic energy: further, external potential value: Then, the single-electron overlap integrals and electron-electron repulsion integrals are calculated: Hamiltonian operators are given by the following matrix expression: Based on the variation principle C V is determined from the equation: Unconstrained Hartree-Fock method is applied including the self-coupling procedure: 1. determining occupied orbital 2. designing density matrix: 3. designing Fock matrix: where, , 4. Integrals based on total energy calculation: i j i j i j i j i j i j i j 5. Verification of E -E old < const, where constconstant value, i.e. convergence criterion of selfcoupling for each integral of calculated total energy.

LCAO-CO approximation comparable to LCAO MO (old version) that crystal orbitals
as Bloch atomic function totals are introduced into the systems with translation symmetry: ( , ) exp( ) ( ) In LCAO-CO one-electron functions for the crystal are constructed as linear combinations of atomic orbitals sums in blocks: Coefficients ( ) np C k are found with the help of variational principle from the Hartree-Fock-Roothaan equations: and instead of one such system for a molecule we obtain N systems for crystals.Matrix ( ) pq F k is found as a sum of matrices

REE-2016
Density matrix for crystal is calculated with the formula: Lattice sums of matrix elements calculated with functions from the zero (central cell): ) While solving the Hartree-Fock-Roothaan equations for crystals it is necessary to know ( ) np C k of the whole reduced Brillouin zone for fixed value because the specified form of matrix elements as follows: pq H k − Hamiltonian elements for the crystal which contains energy of electron interaction with cores V(r), and also it contains kinetic energy.As a result, for self-consistent calculation of the electron crystal structure it is necessary to summarize a large number of occupied electronic states on every stage of integration process.
In the generalized eigenvalue problem the matrices F and S are calculated in the whole reduced Brillouin zone (for fixed value), while Brillouin zone integration in calculating Fock matrix is substituted by summation over spatial points.
Highly precise short-time calculations are one of the major problems in modern quantum chemistry.Applying Hartree-Fock method, including Roothaan equation, (for example, for both MO LCAO and CO LCAO approximations) calculation time is relatively more than in the case if applying density functional theory (proportional N 4 , where N -number of basic functions used in calculations).In this case, Roothaan equation is more rarely applied than density functional theory.This problem for Hartree-Fock method is partially solved by direct inversion in iterative subspace and Schwarz inequality application.Thus, total calculation decreases from N 4 to N 2 but calculation precision does not decrease [8].
Schwarz inequality: This, in its turn, makes it possible to calculate upper value limits (threshold) for all four two-electron central integrals.This threshold could be ignored in the calculation itself.
In this paper the authors applied DIIS (direct inversion in the iterative subspace) to accelerate selfcoupling convergence in Hartree-Fock domain.Its specific feature is plotting on current iteration linear combination of approximated error vectors from the previous iterations (W σ '): Defined coefficients on current iteration are applied for extrapolation of variable functions in the following iteration [9].
Likewise, mathematical manipulations from [7,10] are used in the calculations, which would provide the calculation of a perfect crystal by applying cyclic-cluster method.

Numerical experiment
The first modeling step includes the calculation of CeO 2 crystal geometry by molecular dynamics method.
Molecular dynamics method is applied in automation of input data implementation and in calculations based on Hartree-Fock method.This method defines more precisely the atomic coordinates in the unit cell, which is then used in Hartree-Fock method for further calculations.
Based on molecular dynamics method substance atoms are presented as material points, the interaction of which is described by potential interaction functions.Function types and coefficients are selected for each substance according to experimental thermodynamic data [11].
In molecular dynamics method the time evolution of atom systems is traced by integrating motion equation [12]: Solution algorithm basis includes the following: x coordinates and the velocity of all particles are given at the initial time instant; x in each successor step new coordinates and particle velocity are calculated for all acting energy.The following steps are included within the framework of simulating investigated crystalline structure: x define initial coordinates and atom velocity of modeled system; x select ensemble; x determine the heating procedure and select relevant thermostat; x select simulation time step; x select used integrating unit; x select interaction potential.
Interatomic energy for CeO 2 could be the combination of short-acting and short-range energy [13][14][15].In our case, Buckingham potential is used to simulate short-range energy in a crystal: Coefficients for Buckingham potential is taken from [16].
Long-range interaction energy could be written by Coulomb law: By applying the above-mentioned potentials a more exact simulation of CeO 2 crystal atom systems could be conducted.
The next MD simulation step is integrating motion equation to determine instantaneous rates and atom position in the system.The most common integrating methods are Verlet integration, accelerated Verlet integration, leapfrog integration and linear multistep method.In our case, Verlet algorithm is used in the integrating motion equation (step equals 1 fs) [17]: CeO 2 crystal cluster model (space group Fm3m) was plotted via software package.The number of N atoms in grid block was 1280.Periodic boundary conditions were also applied.Initial atom rate was described by Maxwell distribution at a given temperature.
In simulation the NVE ensemble is used [18].The system is connected to thermostat Nose-Hoover [19], which, in its turn, keeps constant temperature close to 295K in the system: where modeled temperature 295K; pressure 0 Pa; relaxation of a system -during 1 ns.
After calculating the geometrical parameters of crystal lattices, Hartree-Fock method is used to calculate CeO 2 crystal (transformation group Fm3m) in an approximation to CO LCAO.
Quantum-chemical calculation results are depicted in figure 1 and table 2. Table 2 also includes results of previous studies and experimental data.Calculation method applied in this paper is denoted as CO LCAO CM (crystal orbitals as linear combination of atomic orbitals via cyclic-cluster method).The calculations are performed over the coordinates obtained by molecular dynamics method, which, in its turn, introduce an error into the calculation itself.
Minimum energy band p-d is 0.2eV less than experimental range limit.Maximum band-gap deviates from the band by 0.3 eV.This could be considered a good result as the coordinates in the calculation are not true.It should be noted that the experiment results are not single numerical values, but a band.This is conditioned by the complexity of the experiment procedure, as the investigated material is fine powder.

Conclusion
Operability of software package was confirmed by carrying out numerical experiment.The results of calculation the band structure of CeO 2 ideal crystal were compared with experimentally obtained data for a real crystal, which lattice has defects.Despite this, calculation results are in good agreement with the experimental data.Currently, there are models providing the possible calculation of real crystals [28], which would be further implemented into developed software package for simulation materials properties.
At the present time there is quite a large number of software packages for the simulation of atomic structures.Taking into account the problem, which has been indicated in introduction, the software package will be developed further in the area of scalability calculation resources and optimization the algorithms of mathematical models.

Table 1 .
Calculated coordinate results by molecular dynamics method.Energy band diagram for CeO 2 is illustrated in Fig. 1: X-axis indicates wave vector values in Brillouin zone for several directions, while Y-axis -electron energy state magnitude.Dotted line is Fermi energy.Upper group Ce is 5d orbital, middle Ce 4f orbital and lower O 2p orbital.

Table 2 .
Calculation results of CeO 2 energy band structure.