Hopping conductivity in a system with ZnS crystal lattice by non-constant force field molecular dynamics

In the paper non-constant force field molecular dynamics was used to study conductivity behavior on ZnS crystal lattice. The considered conductivity provided by electron hopping between localization centers placed randomly according to ZnS geometry. It was shown that the conductivity behavior depends on the maximal hopping distance. For the small distances the conductivity passes through the maximum around equimolar concentrations of electron donors and acceptors. Increasing in the maximal hopping distance leads to increasing in conductivity values and change shape of its concentration dependence.


Introduction
Concentration dependences of conductivity in different systems has already for a long time obtained the attention of researchers [1][2][3].Such questions are often referred as percolation problems, including conductivity provided by polaron hopping between localization cen ters [2,3].At the moment a big amount of equations were suggested for describing concentration dependences of conducti vity in various systems.However, analytical solutions can not account all multiplicity of factors acting on conductivity: thermal motion of ions, their mobilities, sizes and etc.An attempt to consider all microscopic picture of the conductivity process can be made with numerical methods which began to develop with progress in computer engineering.One of the most widespread me-thods for studying of many-boded systems on a microscopic level is molecular dynamics simulation (MD) [4,5].Unfortunately, classical MD can not deal with polaron hopping conductivity, because this phenomenon has quantum-mechanical nature.Early we have suggested scheme which allows including hopping conducti vity into MD [6].This approximation is possibility of particles to change their oxidation degree (and, consequently, properties and interaction laws) runtime.Thus force field becomes variable and method can be called as "Non-constant force field molecular dynamics" (NCFFMD).The method is implemented in "azTotMD" software which is available on website http://ncffmd.ru/ [7].In our previous work we have demonstrated possibilities of the method and the software for simulation of redox processes in liquid media [8].The aim of the present paper is to study percolation be-havior of a system with polaron hopping between localization centers placed according to ZnS crystal lattice.

Simulation details
MD simulations were performed with "azTotMD" software [7] in canonical (NVT) ensemble.Newton's equations of motion were integrated by Velocity Verlet algorithm [9] with timestep of 1 fs during 100'000 steps (0.1 ns).Equilibration time was 1 ps (1000 timesteps).Electrostatic interactions were accounted using the Ewald summation.Nosé-Hoover thermostat [10] with relaxation time of 0.2 ps was used for maintaining the temperature around 298 K.The considered system consists of electron donors (A + ), electron acceptors (A 2+ ) and counterions (X -).The number of X -ions was chosen equal to 500 for all studied systems.Amounts of A + and A 2+ cations were given in such way to keep electroneutrality and obtain desired ratio of A + / A 2+ concentrations.Short range interactions were given by pair potential suggested for CuCl-CuCl 2 binary system [11] since solid CuCl has a ZnS structure.Initial configurations were generated with the abovementioned site [7].The starting ion coordinates correspond to ZnS crystalline structure with some vacant sites.The box was cubic with the edge length of 25.7 Å for all studied systems.This box length roughly corresponds to a cell parameter of CuCl.The sizes of boxes was the same for simplification and, in addition, ionic radii of Cu + and Cu 2+ (which are prototypes of ions A + and A 2+ ) are close to each other and equal to 0.77 and 0.73 Å, correspondingly [12].
For activation of electron transfer routine during the simulation value of ejump directive was set as 1 in control.txtfile (one of input files for the program).The program performs electron transfer if the system decreases energy by this transfer.The difference in system energy before and after electron transfer is determined according to formula [6]: where ∆U ij is energy difference after electron transfer from i-th particle to j-th particle, V ik is the Van der Waals energy of interaction between i-th and k-th par ticles, provided by corresponding pair potentials, q is the electric charge of the particle, r is the distance between par ticles, C is the constants in Coulomb's law, ε is the relative permittivity of the media, ΔE x is the voltage drop on the X axis, a is the box length and upper indices I and II mean states of particles before and after electron transfer, correspondingly.Electronic current (I) was determined through a time derivative of difference in the number of electrons transferred in positive and negative directions: where e is the electron's charge, n + and n - are the numbers of electron hops through Oyz edge in positive and negative directions along x-axis.

Results and discussion
At first let us discuss a structure of the system.The numerical experiments showed that with the presented here force field the crystal lattice is stable only in the case of the even numbers of both A + and A 2+ cations.In this case one can observe crystal lattice with ZnS structure independent of the time of the system evolution, Fig. 1.One can see some vacant sites in cation sublattice, because every A 2+ ion demands one cation vacancy to save electroneutrality.Despite of this voids the system keeps own crystal lattice right up to composition of 0.2AX•0.8AX 2 .
Examples of radial distribution functions (RDF) are given in Fig. 2. The RDFs are characteristic for crystalline systems; there are a number of well-resolved maxima.Distance for the first maximum determines more probable distance between ions in the first coordination sphere for a chosen pair.The position of the first maximum for A + -A 2+ pair does almost not depend on composition and electric field and equals to ~3.7 Å.The distance at which RDF for A + -A 2+ pair starts to exceed 0.1 is also almost constant and equal to ~2.7 Å.This means that if length of electron hopping is lesser than 2.7 Å, the system will not have electronic conductivity.For this reason, we set maximal hopping distance to be much higher than 2.7 Å by specifying of rElec directive in the control.txtfile.
Usually, the number of electron hops through some plane is a linear function of time, Fig. 3. Without external electric field these numbers are the same for positive and negative directions, but in the case of the field the slopes of these lines differs from each other (Fig. 3).The corresponding value of current can be obtained from expression (2).Note that for observation of noticeable current we need to apply an external electric field with a colossal voltage, because current density of 1 A•cm -2 is approximately equal to 0.625•10 -9 singlecharged particles per picosecond, per square angstrom (in units more conve nient for MD).To see the noticeable number of particles for the simulation time we need to obtain high current which requires high Fig. 1.The structure of the simulated system 0.5AX•0.5AX 2 obtained after simulation under 1V electric field during 100 000 timesteps.Big green spheres denote X -anions, small red and orange ones -A 2+ and A + cations, correspondingly Fig. 2. The radial distribution functions, g(r), of the simulated system 0.5AX•0.5AX 2 for A + -X -and A 2+ -A + ion pairs 0,0 voltages.However, extremely high voltages can lead to the destruction of a system.In our simulations we used voltage of 1 V, higher voltages led to breaking of crystal lattice and liquid-like structure of the system.Electronic current as a function of composition for different values of the maximal hopping distance is presented in Fig. 4. As expected, the current grows with this parameter.At small values of the distance the current passes through the maximum around x = 0.5.From statistical point of view the current will be maximal if pro bability of finding the electron acceptor (A 2+ ) near the electron donor (A + ) is maximal.This probability is proportional to the product of their concentrations which is maximal at x = 0.5.The difference of observed position of maximum from 0.5 can be caused by some reasons: asymmetry in pair potentials, different mobility of A + and A 2+ ions and etc.At high values of the maximal hopping distance the current grows with concentration of A + species.This implies that if electron can hop on enough long distances the conductivity will be limited by the concentration of electron donors.In other words, if there is an electron donor, an electron acceptor always can be found.So the maximal length of the electron hop determines the shape of the concentration dependence of the conductivity.

Conclusions
Non-constant force field molecular dynamics is able to simulate electron hops between electron donors and acceptors affected by thermal movement.Thus, it is possible to study polaron conductivity of a given system by this method.In this work concentration dependence of conductivity was considered in the case of ZnS geometry for electron localization cen ters.It was shown that the position of the conductivity maximum depends on the maximal hopping distance.i, e • ps -1

Fig. 3 .Fig. 4 .
Fig.3.The number of electron hops through the Oyz plane in positive and negative directions (and their difference) as a function of time.The simulated system is 0.5AX•0.5AX 2 , the maximal hopping distance is 4.5 Å, the applied voltage is 1 V