Hopping conductivity in a system with ZnS crystal lattice by non-constant force field molecular dynamics 153 D O I: 1 0. 15 82 6/ ch im te ch .2 01 8. 5. 3. 05 A. A. Raskovalov Institute of High-Temperature Electrochemistry, Ural Branch of the Russian Academy of Sciences, 20 Akademicheskaya St., Ekaterinburg, 620137, Russian Federation E-mail: other@e1.ru 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. Keywords: molecular dynamics; non-constant force field; polaron hopping; ZnS lattice. Received: 08.10.2018. Accepted: 26.10.2018. Published: 31.10.2018. © Raskovalov A. A., 2018 Raskovalov A. A. Chimica Techno Acta. 2018. Vol. 5, No. 3. P. 153–157. ISSN 2409–5613 Introduction Concentration dependences of  conductivity in  different systems has already for a long time obtained the atten- tion of researchers [1–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 concentra- tion dependences of conducti vity in vari- ous 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 en- gineering. One of the most widespread me- thods for studying of many-boded systems on a microscopic level is molecular dyna- mics simulation (MD) [4, 5]. Unfortunately, classical MD can not deal with polaron hopping conductivity, because this phe- nomenon has quantum-mechanical nature. Early we have suggested scheme which al- lows including hopping conducti vity into MD [6]. This approximation is possibili- ty 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 dy- namics” (NCFFMD). The method is im- plemented in “azTotMD” software which is available on website http://ncffmd.ru / [7]. In our previous work we have de- monstrated possibilities of the method and 154 the software for simulation of redox pro- cesses 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 accord- ing to ZnS crystal lattice. Simulation details MD simulations were performed with “azTotMD” software [7] in  canoni- cal (NVT) ensemble. Newton’s equations of motion were integrated by Velocity Ver- let algorithm [9] with timestep of 1 fs du- ring 100’000 steps (0.1 ns). Equilibration time was 1 ps (1000 timesteps). Electro- static interactions were accounted using the Ewald summation. Nosé-Hoover ther- mostat [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 (A2+) and counterions (X-). The number of X- ions was chosen equal to 500 for all studied systems. Amounts of A+ and A2+ cations were given in such way to keep electroneutrality and obtain desired ratio of A+ / A2+ concentrations. Short range in- teractions were given by pair potential sug- gested for CuCl-CuCl2 binary system [11] since solid CuCl has a ZnS structure. Ini- tial 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 Cu2+ (which are prototypes of ions A+ and A2+) are close to each other and equal to 0.77 and 0.73 Å, correspon- dingly [12]. For activation of electron transfer rou- tine during the simulation value of ejump directive was set as 1 in control.txt file (one of input files for the program). The pro- gram performs electron transfer if the sys- tem decreases energy by this transfer. The difference in system energy before and after electron transfer is determined according to formula [6]: �U V V V V C q q q r q ij ik II ik I jk II jk I k i j k i II i I ik j I � � � ��� �� � � � � � � , 1 � II j I jkk i j x i j q r E a x x �� � � � � � � � � � �� � � , , � (1) where ∆Uij is energy difference after elec- tron transfer from i-th particle to j-th par- ticle, Vik is the Van der Waals energy of in- teraction 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 rela- tive permittivity of the media, ΔEx 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 cur- rent (I) was determined through a  time derivative of  difference in  the number of  electrons transferred in  positive and negative directions: I e d n n dt � �� �� � , (2) where e is the electron’s charge, n+ and n- are the numbers of electron hops through Oyz edge in positive and negative direc- tions along x-axis. 155 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 A2+ cations. In this case one can observe crys- tal 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 A2+ ion demands one cation vacancy to save electroneutra- lity. Despite of this voids the system keeps own crystal lattice right up to composition of 0.2AX·0.8AX2. Examples of radial distribution func- tions (RDF) are given in Fig. 2. The RDFs are characteristic for crystalline systems; there are a number of well-resolved ma- xima. Distance for  the first maximum determines more probable distance be- tween ions in the first coordination sphere for a chosen pair. The position of the first maximum for  A+–A2+ pair does almost not depend on  composition and elec- tric field and equals to  ~3.7 Å.  The dis- tance at which RDF for A+–A2+ 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 con- ductivity. For this reason, we set maximal hopping distance to be much higher than 2.7 Å by specifying of rElec directive in the control.txt file. 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 posi- tive and negative directions, but in the case of the field the slopes of these lines differs from each other (Fig. 3). The correspond- ing 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 vol- tage, because current density of 1 A·cm–2 is approximately equal to 0.625·10–9 single- charged 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.5AX2 obtained after simulation under 1V electric field during 100 000 timesteps. Big green spheres denote X- anions, small red and orange ones — A2+ and A+ cations, correspondingly Fig. 2. The radial distribution functions, g(r), of the simulated system 0.5AX·0.5AX2 for A+–X– and A2+–A+ ion pairs 1 2 3 4 5 6 7 8 9 0,0 0,5 1,0 1,5 2,0 2,5 3,0 3,5 4,0 4,5 5,0 5,5 g( r) , Å r, Å pair A+ – X– pair A2+ – A+ 156 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 lat- tice and liquid-like structure of the system. Electronic current as a function of com- position for different values of the maxi- mal 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 (A2+) 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 posi- tion of maximum from 0.5 can be caused by some reasons: asymmetry in pair po- tentials, different mobility of A+ and A2+ 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 dy- namics 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  con- ductivity was considered in the case of ZnS geometry for electron localization cen ters. It was shown that the position of the con- ductivity maximum depends on the maxi- mal hopping distance. 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.5AX2, the maximal hopping distance is 4.5 Å, the applied voltage is 1 V Fig. 4. Electronic current (i) as a function of composition (x) in the xAX–(1–x)AX2 system under external electric field of 1 V and different maximal hopping distance (re) 0 10 20 30 40 50 60 70 80 90 100 0 500 1000 1500 2000 2500 di�erence in negative direction n τ, ps in positive direction 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0 1 2 3 4 5 6 7 x(AX) re = 3.5Å re = 3.8Å re = 4.0Å re = 4.5Å i, e · p s– 1 157 Acknowledgements The reported study was funded by Russian Foundation for Basic Research (RFBR), according to the research project No. 16-33-60095 mol_а_dk. References 1. Kirkpatrick S. Percolation and Conduction. Rev Mod Phys. 1973; 45(4):574–88. DOI: 10.1103 / RevModPhys.45.574. 2. Kurkijarvi J. Conductivity in random systems. II. Finite-size-system percolation. Phys Rev B. 1974; 9(15):770–4. DOI: 10.1103 / PhysRevB.9.770. 3. Shklovskii BI, Éfros AL. Percolation theory and conductivity of  strong- ly inhomogeneous media. Sov Phys  Usp. 1975; 18(11):845–62. DOI: 10.1070 / PU1975v018n11ABEH005233. 4. Allen MP, Tildesley DJ. Computer Simulation of Liquids. Oxford: Clarendon Press; 1987. 385 p. 5. Frenkel D, Smit B. Understanding Molecular Simulation From Algorithms to Ap- plications. San Diego, California: Academic Press; 1996. 638 p. 6. Raskovalov AA. A new extension of classical molecular dynamics: an electron transfer algorithm. J Comp Chem. 2017; 38(5):926–32. DOI: 10.1002 / jcc.24755. 7. Non-constant force field molecular dynamics. Available from: http://ncffmd.ru / set. php?lang=en (accessed: 22.02.2017). 8. Raskovalov AA, Latypov AA. Simulation of Red-Ox Reactions in Liquid Media with Non-Constant Force Field Molecular Dynamics. Rasplavy [Melts]. 2018; 6. In press. DOI: 10.1134 / S0235010618060063. 9. Verlet L. Computer “Experiments” on Classical Fluids. I. Thermodynamical Proper- ties of Lennard-Jones Molecules. Phys Rev. 1967; 159(1):98–103. DOI: 10.1103 / Phys- Rev.159.98. 10. Martyna GJ, Tuckerman ME, Tobias DJ, Klein ML. Explicit Reversible Inte- grators for  Extended Systems Dynamics. Mol Phys. 1996; 87(5):1117–57. DOI: 10.1080 / 00268979600100761. 11. Raskovalov AA, Shevelin PYu. Physico-chemical Properties of the Molten CuCl — CuCl2 System: Experiment, Thermodynamics and Molecular Dynamics Simulations. J Sol Chem. 2018. In press. DOI: 10.1007 / s10953-018-0817-x. 12. Shannon RD. Revised effective ionic radii and systematic studies of  interato- mic distances in halides and chalcogenides. Acta Cryst. 1976;A32:751–67. DOI: 10.1107 / S0567739476001551.