polymerization process Research
IN OXIDE MELTS BY MOLECULAR DYNAMICS
AND STATISTICALGEOMETRICAL METHODS
Voronova L.I., Voronov V.I., Gluboky J.V.
The method of forecasting of structure and physicalchemical properties of oxide melts on a base of the complex simulation by a quantumchemical method MNDO, moleculardynamic and statisticalgeometrical methods realized by the authors by the way of a program complex is described. The methods of simulation of structure with different levels of an approximation are considered: the shortrange order, nanostructure, microinhomogeneity. The method of сovalent bonds network covering developed by the authors for exploration of oxide melts polymerization processes is explicitly described. Some results of structural simulation of the systems SiO_{2}CaO and SiO_{2}Na_{2}O are shown.
Introduction
The oxide melts are disordered stronginteracting polymerized systems. One of the basic problems of physical chemistry of oxide melts is the exploration of correlation of structure and physicalchemical properties of these objects.
The absence of the sequential analytical theory of slag melts and heavy complexities of their experimental study stimulated active development of a new scientific direction "computer science of materials". The results of computer experiments are used for design of materials with required properties.
Among computer methods used for study of stronginteracting systems of many particles, methods of computer simulation have widely spread, such as the quantumchemical semiempirical methods, in particular, the method MNDO (modified neglect by diatomic differential overlapping); Monte Carlo method (MC); a method of molecular dynamics (MD). These methods, practically from "of the first principles", allow to obtain various physicalchemical properties. The adequacy of obtained results is defined by accuracy of used mathematical models.
One of actual problems, successfully explored within the framework of "computer science of materials", is the study of the polymerization processes of multicomponent oxide melts. It is possible to use both the MC and the MDmethods for these purposes, however the MDmethod is more preferable, because with its help it is possible not to only define parameters of structure formation, but to explore multilaterally a fundamental problem of “structureproperty” correlations.
It is possible by results of moleculardynamic simulation, with engaging statisticalgeometrical (SG) methods to investigate system structure with different levels of detailing (shortrange order, nanostructure, microinhomogeneity), to define features of shortrange and extended structure, to pick out characteristic building blocks and regularities of their relative location.
Complex MNDOMDSg simulation
For forecasting oxide melts structure and their physicalchemical properties we use a complex method that includes a semiempirical quantumchemical MNDO (modified neglect by diatomic differential overlapping) method, molecular dynamics (MD) method and statisticalgeometrical (SG) methods.
The complex method simulation phases 
MNDO simulation in a cluster approximation, obtaining the body of basic data;

Superimposed potentials parametrization on the basic data processed by multivariate optimization methods;

Molecular dynamics experiment, arrays of physical chemical properties as result;

Structure investigation by the method of сovalent bonds network covering

Structure investigation by statisticalgeometrical methods;

Investigation of “compositionstructureproperty” correlation;
The simulation outcomes 
Parameters of potential functions in an ionic and an ioniccovalent approximation;

Energy parameters;

Thermodynamic parameters: molar heat capacities, adiabatic and isothermal compressibilities, temperatureexpansion factor;

Transport properties: autocorrelation functions, rootmeansquare displacements, diffusion factors, volume and shear viscosity, electric conductivity, thermal conductivity;

Spectrum characteristics: a density function of oscillatory states, infrared absorption spectrums;

The averaged structural parameters of the shortrange order: partial functions of radial and angle distribution, bond lengths, average coordination numbers and their allocations;

Estimation of the system polymerization degree: the distribution functions of complex anions on some parameters, lifetime of polyanions, portion of Oxygenium of a different type, portion of plane rings;

Statistical geometrical parameters of extended and longrange structure: allocation of number of sides, volumes and squares of polyhedrons, squares of sides, spherical tetrahedral and octahedral factors;

Correlation dependencies of physicalchemical properties and basis nanostructure of melts as tables and diagrams;

Visualization of the system evolution
Ionic covalent model (ICM)
For forecasting properties of the multicomponent oxide melts containing networkforming Si, B, Al ions and improving results adequacy we use perspective ioncovalent model (ICM) applied for systems, containing stable longliving clusters with a high share of covalent connections [1]. The pair sphericalsymmetrical longrange Coulomb interaction, the two and threeparticle covalent interactions connected to quantum mechanical dipole and quadrupole moment are taken into account in ICM.
The particles of modeling system – cations networkformers(CNF), cationsmodifiers (CM) and anions of oxygen(O) have the following attributes in ICM: m_{i} mass, q_{i}–effective charge, _{i} – effective radius,  radiusvector,  velocity; besides there are set: d_{0} – CNFO bond length, _{0} equilibrium OCNFO bond angle and force constants: _{kit} – twoparticle covalent interaction, _{kit} – threeparticle covalent interaction.
The particle’s potential is defined by particle’s attributes and its implementation of elementary structural unit “implication condition“. Since all particles have a charge as attribute, The Coulomb interaction is described by pair ionic spheresymmetric potential
=+, (1)
that have terms in the Pauling's form, most often used at MDsimulation of oxide melts [2].
, (2)
where r_{ij}_{ }is a particle distance, e is the electron charge, n is repulsive curve steepness parameter.
The implication condition: In case the particle belongs to a "regular" elementary complex that includes atoms of oxygen in number of atom networkcreator valency, potential (1) is supplemented by two and threepartial covalent contributions, associated with the quantummechanical dipole and quadrupole moments:
twoparticle potential: , (3)
threepartial potential: , (4)
where m is the index of silicon atom in center of an elementary structural unit; k and k' are elementary structural unit oxygen atom indexes; q_{kmk} is OCNFO equilibrium angle; 1.5d_{0}, 1.5_{0} are maximum covalent two and threeparticle forces action radius and angle respectively.
The introduction of the last terms in the equations (3) and (4) allows to explicitly describe an advantage in energy under considering the covalent interactions, saving continuity of potential function.
Covalent contributions are described in the Keating's approximation [3]:
, (5)
. (6)
Potential energy _{m} of covalent interaction of belonging to a melementary complex particles:
. (7)
Complete potential energy of a simulated system, considering the covalent interactions inside an elementary complex:
. (8)
Basic MNDO accounts and parametrization of potential functions
The parameters of superimposed potentials (26) are determined with the help of semiempirical quantumchemical method MNDO, based on the electronic structure computation, which gives the most stabilized and close to the experiment results.
The basis of all semiempirical quantumchemical methods is the cluster approximation. Representative cluster or complex (fig. 1) retaining main characteristics of the system is cut out and its Schrodinger equation is solved. The series of used approximations and limitations are described in the special literature [4,5].
Figure 1. The representative clusters for MNDOsimulation of a system SiO_{2}CaO.
Some characteristics related to multiparticle interaction are defined by the MNDOmethod. Namely they are: total energies of complexes, heats of formation, energies of atomization, effective atom charges, orbital population, equilibrium bond lengths valence angles, force constants. However, it is impossible to use directly these basic MNDOdata in MDsimulation for some reasons. We have developed special technique that allows to interpret the results of basic MNDOcomputations in the ioniccovalent model terms [1].The parametrization of a superimposed potential functions in ICM is made with the help of the program of multivariate optimization. Here the force constants of two and threeparticle covalent interaction are varied parameters, that are searched by minimization of the criterion function. This function is built as the sum of squared forces deflections calculated through potential gradient (26) and corresponding forces in a linear approximation obtained in MNDOcalculations.
MD simulation
The “particles model” is used for molecular dynamics simulation, i.e. there is the mutual unique dependence between physical particles and computer model particles. The classical differential equations of motion based on Newtonian laws are solved for each particle:
, (9)
which we approximate by finitedifference equations on Beeman algorithm [7]:
,
(10)
Where m_{i}, _{ }  mass and acceleration of particle i,  complete force affected at particle i,  simulation time step (is usual~10^{15} c). The periodic boundary conditions are used [8], the real volume of cube depends on the melt density.
Structure research techniques
One of obtained by the MD computer experiment results is collection of coordinates and velocities of the particles of simulated system for each configuration, that allows to investigate both structure and dynamics of a system in detail.
The averaged structural parameters of the shortrange order
The averaged structural parameters of the shortrange order determined traditionally on the base of these data are described below.
The radial distribution function (RDF) of particles of type in an spherical layer of radius r around particles of type :
, (11)
where is an average number of particles of type, located at distance from r till from particles of a type; V  volume of the system.
Coordination number is an average number of particles located inside a sphere with the radius equal to distance up to the first minimum of radial distribution function.
, (12)
where  distance up to the first minimum of RDF: an average distance between particles of type (bond length or radius of the first coordination sphere).
The angles between particles are defined for particles of type that located inside of the first coordination sphere of particle of type . The cosine law is applied to calculate angles on known coordinates of particles.
The angle distribution function (ADF) are calculated by the formula:
, (13)
where is number of angles between particles in limits from to for particles of type, located in the first coordination sphere of a particle of type, at step k, is number of all angles at step k.
The structural parameters of the shortrange order for the melt with the mole fraction 0.3 N_{SiO2 }of binary system SiO_{2}CaO are indicated in fig.2 and fig.3. The square of the first peak of the curve 1 in fig. 2a defines the coordination number of silicium on oxygen, r_{0} – an average distance between atoms of two types (bond length).
However the listed averaged shortrange order parameters do not contain information about structure features, because they reflect the averaged pattern on the whole volume of a sample. They say nothing about extended and longrange structure, which is defining in forming of physicalchemical characteristics of slag melts, that’s why the other approaches for exploration of structure are actively developed.
The computation of distributions of coordination numbers for atoms – networkcreator is informative enough [10]. It gives an opportunity to obtain parameters defining features of the shortrange order having a good agreement with the experimental data.
The basis of this technique, as well as in a classical case, is the radial distribution functions of atoms. The tables, combining information about relative numbers of ions of different types of this or that coordination; share of contents of elementary structural units as function of molar ratios of slag components; contents of nonbridging oxygen as function of glass composition are used as estimated parameters.
a ) b)
Figure 2. a) The radial distribution function (RDF) for 1)SiO, 2)OO, 3)CaO,
b) The angle distribution function (ADF) 1)OSiO, 2)SiOSi
a ) b)
Figure 3. Temperature dependences of a) average distances (bond length), b) coordination numbers, for 1)SiSi, 2)SiO, 3)OO, 4)CaO.
The method of Voronoy’s polyhedrons and Delone’s simplexes
For the structure detailing it is expedient to use statisticalgeometrical methods, in which with help of a special fragmentation of a system on polyhedrons it is possible to get both regularities and features of shortrange and extended structure of the explored system. There are number of fragmentation algorithms, some of them are Tanemura's algorithm, the method of circumscribed spheres, facets bypass method.
The most potent and deeply developed is the method of Voronoy's polyhedrons and Delone's simplexes [11]. The fragmentation is the set of the atom centers of a system {A}. For each center {A} it is always possible to indicate the range, which points are closer to the given center, than to any other center of a system. These are Voronoy's polyhedrons (fig. 4 a) the fragmentation of a system on tetrahedrons (Delone's simplexes) is simultaneously created, the apexes of tetrahedrons are the atoms, and the sphere circumscribed around of a tetrahedron, does not include other atoms of a system (fig. 4b). A fragmentation of a system on Voronoy's polyhedrons and the Delone's simplexes are unambiguous and fill in the whole space without hollows.
a) b)
Figure 4. Example of construction a) of Vononoy’s area around the arbitrary centre of a system, b) Delone’s simplex (123) the circle, circumscribed around of atoms, does not include other atoms
There is an unambiguous correlation between a Voronoy’s fragmentation and Delone’s fragmentation since each apex of a Voronoy’s polyhedron is the center of circumscribed around of a particular Delone’s simplex sphere. The explorations of structure through Voronoy’s and Delone’s fragmentations are based on analysis of metrical and topological performances of obtained formations. For Voronoy’s polyhedrons the following characteristics are used:  an average number of facets on a Voronoy’s polyhedron; the sphericity coefficient , where S the area of the surface of a polyhedron, V its volume; a topological index n_{3}n_{4}n_{5}n_{6}n_{7} …, where n_{i} number of icornered facets at this polyhedron; the matrix of neighbourhoods {n_{ij}}, where each element n_{ij} is number corresponding amount of jcornered facets around of all icornered facets on this polyhedron; an index of neighbourhoods  diagonal elements of neighbourhoods matrix.
As against Voronoy’s polyhedrons, the Delone’s simplexes have the same topological type, therefore they can be distinguished among themselves only metric – by size and form. In these purposes such characteristics are used , as tetrahedricity
,
where l_{j} is edge lengths of the given simplex, and l_{0} – an average length of its edges;
and octahedricity
,
where l_{m} is the length of the longest edge of a simplex. The less values T and O, the closer simplex to a regular tetrahedron and octahedron accordingly. The method allows to get information about current structure of a system, various structural characteristics, as system polymerization degree, free ions presence, passability of a system for particles of a particular type, diffusive characteristics.
The covalent bonds network covering method
We have developed a “covalent bonds network covering” method for exploration of structure of polymerizing systems [9].
The method is that for all simulated structures the complex (i.e. complex anions) distribution functions (CDF) are obtained on various characteristic parameters рarm of complexes T, existing in a system; рarm = Tyрe, N, , , , , where Tyрe  the type of a complex, is defined by number of varied particles which are included in a complex T (Tyрe), N  total number of particles in a complex T (N),  sum of cations networkcreators in a complex T(),  sum of free oxygen in a complex T(),  sum of nonbridging oxygen in a complex T (),  sum of bridging oxygen in a complex Т (). It is possible to define degree of polymerization of a system on CDF of the following kinds:
1. , where is the portion of T(рarm) complexes; is the number of configurations, on which the complex exists, the sum in the nominator is taken over all the complexes of this type; is the number of configurations on which any complex T_{j} exists, the sum in the denominator is taken over all j existing complexes.
2. The lifetimes for various complexes are defined: , where is average lifetime of a complex T(рarm), is lifetime of a complex T(рarm)_{i}, the sum is taken over complexes with the same value of the parameter.
3. The portions of differenttypes oxygen atoms are calculated:
, ,
where is the portion of bridging oxygen, is the portion of nonbridging oxygen and is the portion of free oxygen in a system, , are quantities of bridging and nonbridging oxygen atoms respectively, N_{O} is the amount of oxygen atoms in a system, K is a total amount of system configurations in a phase.
4. The portions of closed structural units D(grouр)_{n} of a different degree of complexity are
calculated: , where N gr(n, k) is a number of closed structural
units of a ndegree of complexity on k configuration.
Similarly it is possible to define portions of plain rings of different dimensionality:
, where N ring(n,k)  number of rings
of n^{th} dimensionality on a kconfiguration, and complexes charge:
, where n is the number of links of those particle types, which are indicated in inferior indexes of the sums.
a) b)
Figure 5. а) portions D_{T(cat)} of complexes of a different type, b) configurational lifetime of complexes of the system SiO_{2}Na_{2}O
Some results of investigation of the process of the polymerization in the system SiO2Na2O is shown in fig. 6. Seven melts were simulated in the whole range of compositions. At increase of a molecular fraction SiO2 the process of large silicooxygen groupings forming, leading to networking, is observed. Simultaneously, there is a process of their constant regenerating, as the lifetime even of large classifications (more than 50 atoms of silicium) in an middle region of compositions does not exceed 400 configurations. It is connected with the fact that both the atoms of oxygen, and smallsized silicooxygen anions having 12 atoms of silicium, permanently associate and abandon large complexes, owing to the migration of atoms of sodium.
a) b)
F igure 6. Portions a) of plain rings D(ring)_{n} with different number of cations networkformers n, b)of closed structural units D(group)_{n} with number of cations n;
n= 1) 4, 2) 5, 3) 6, 4) 7, 5) 8.
In fig. 6 the portions of plain rings and portions of closed structural units of different complexity in dependence on composition of a melt of a system SiO_{2}CaO are shown. The plain rings exist only in range of compositions (0.70.9) N_{SiO2} and include from 4 up to 8 cations of silicium. For structure 0.7 N_{SiO2} the pentatomic rings, and in structure 0,9 N_{SiO2}  fouratomic ones dominate. The portions of rings with the major contents of cations networkcreator are minor. The closed structural units differ from rings in that the atoms, included in them, do not belong to one plane. As follows from fig. 6, such structural formations are present already at structure 0.5 N_{SiO2} and have 4 or 5 silicium atoms. For a range of structures (0.70.9) N_{SiO2} their portion decreases and larger closed units containing 6,7,8 atoms of silicium are formed.
With usage of the stated approaches we obtained series of results for oxide binary mixtures containing atoms of Si, Al, Na, K, Ca, Mg, which were the basis for creating the structure  property correlation dependencies. The analysis of the obtained results shows efficiency of usage of the different approaches for structure analysis, so the obtaining of multivariate characteristics of structural parameters allows to penetrate deeply into the nature of polymerized oxide melts and to obtain more reliable prognoses of properties.
Acknowledgments
The work is executed by support RFBR, grant № 970332531 REFERENCES

Voronova L.I., Gluboky Y.V., Voronov V.I., Grokhovetsky R.V: 'Raschet samosoglasovannogo nabora potentzialnyh parametrov dlya MNDOMD modelirovaniya binarnyh oksidnyh rasplavov'. Rasplavy 1999 2 6674.

Mitra S.K: 'Molecular dynamics simulation of silicon dioxide glass'. Phyl.Mag. B 1982 45(5) 529548.

Gaskell P.H., Tarrant I.D: 'Refinement of random network model for vitreous silicon dioxide'. Phil.Mag. 1980 42(2) 256286

Dewar M.J.S. and Thiel W. J: 'Ground States of molecules. 38. The MNDO Method. Approximations and Parameters'. J.Am. Chem. Soc. 1977 99(15) 48994907.

Semiempirical methods of electronic structure calculation. Edited by G.A.Segal, Moskow, Mir, 1980.

Bliznuk A.A., Voytuk A.A: 'Complex software MNDO85 for calculation of electronic structure, physicalchemical properties and reactivity of molecular system by semiempirical methods MNDO, MNDOC and AM1'. Journ.Struct.Chem. 1986 27(4) 190191

Alder B. J, Wainwright T. E: 'Studies in molecular dynamics. I. General Method'. J. Chem. Phys. 1959 31 459466.

Beeman D: 'Some multistep methods for use in molecular dynamics calculations'. J.Comput.Phys. 1976 20 130139.

Voronova L.I., Voronov V.I., Gluboky Y.V., Grokhovetsky R.V.: 'Ocenka stepeni polimerizatzyi sistemy SiO_{2}CaO po rezultatam molekulyarnodinamicheskogo modelirovaniya 1. Matematicheskaya model i ee programnaya realizatziya'. Sbornik nauchnyh trudov “Matematicheskoe i programnoe obespechenie nauchnyh issledovanii i obucheniya”, Kurgan State University, 1998 4954.

T.F.Souls, A.K.Varshneya: 'Molecular Dynamic Calculations of A Sodium Borosilicate Glass Structure'. Journal of the Amer. Ceram. Soc. 1987 64(3) 145150.

Medvedev N.N: `Metod VoronogoDelone v issledovanii struktury nekristallicheskih upakovok`, Novosibirsk, NSU, 1994.
