MD simulation of the radiation influence on the thermophysical properties and the structure of water L. A. Bulavin,1 K. V. Cherevko,1 D. A. Gavryushenko,1 V. M. Sysoev,1 and T. S. Vlasenko2,∗ 1Physics Faculty, Taras Shevchenko National University of Kyiv, 4 Glushkova Av., Kyiv, 03022, Ukraine 2Institute for Safety Problems of Nuclear Power Plants, 7 National Academy of Sciences of Ukraine, 12 Lysogirska St., Kyiv, 03028, Ukraine 1 (Dated: January 25, 2017) 0 2 The results of themolecular dynamicssimulation of theradiation influenceon the structureand n thermophysicalproperties of water are presented. Thechanges in theradial distribution functions, a momentum distribution function and the selfdiffusion coefficients are quantified. It is shown that J the irradiation causes the structural and thermodynamics properties of water. The ”effective tem- perature” of the stationary nonequilibrium water system under the irradiation allowing to define 4 2 the correspondent equilibrium system with the same structural and thermodynamic properties is calculated. It is confirmed that the structural changes in the liquid systems under irradiation are ] causedbythechangesinthecoefficientsoftheMaxwelldistributionfunctionduetothemomentum h exchangebetweentheactiveparticles andtheparticles formingtheliquid. Toexplain thephenom- c ena observed in the molecular dynamics simulation the results are quantitatively compared to the e predictions of the theoretical model of the phenomenon that is based on the Bogolyubov chain of m equations and with theexperimental data. - t a t Recently, quite a number of works devoted to the ra- sponsibleforthechangesinthestructuralandthermody- s . diation interaction with the matter appeared in the lit- namicpropertiesofwaterunderirradiationbyαparticles t a erature [1–6]. That field is important both from funda- (He). m mental point of view and for the different applications. Earlier,basedonthefundamentalBogolyubovchainof - The studies may reveal the underlaying physical mecha- equations[31]therewassuggestedthemodelrelatingthe d nisms of the interactions and their results may be used structuralandthermophysicalpropertiesofthe nonequi- n for developing new technologies in medicine, nuclear en- libriumliquidsystemsunderirradiationinthestationary o c ergy, etc. Nowadays most of the studies in the field deal state [32]. Within that model it was suggested that the [ with the solid state [7–9], ionic liquids [10, 11] or biolog- irradiation changed the structure of the liquid system. 1 ical structures [12, 13]. For the above objects computer That obtained new structure of the nonequilibrium liq- v simulation studies has shown to be on of the most ef- uidsysteminthestationarystatewascharacterizedbya 5 fective instruments [14–16] as the experiments with the new parameter that is the “effective temperature” T . eff 6 irradiation influence on matter are difficult to perform That allowed writing down the classical BGY equation 9 and to interpret. At the same time there are quite few [33]forthenonequilibriumsysteminthestationarystate: 6 works devoted to the irradiation influence on liquid sys- 0 1. tweimdesl.yIutssehdotuoldstbuedynothteedeqthuailtibcroiummpuatnedr sniomnuelqautiiloinbsriuarme kTeff∂F2∂(rr1,r2) + ∂Φ(|r∂1r−r2|)F2(r1,r2) 0 1 1 17 lcioqmuipdustyesrtseimmsu[l1a7ti–o2n2s].aTrehneoretfuorsee,dittolosotkusdysutrhperirsaindgiatthioant +ρ ∂Φ(|r∂1r−r3|)F3(r1,r2,r3)dr3 =0, (1) 1 : influenceonliquids. Furthermore,atpresent,mostofthe Z v i works in the field deal with the physical-chemical stage with Teff standing instead of T. Here F2(r1,r2) and X and describe the radiolysis [23–25] but only few can be F3(r1,r2,r3) are the 2-nd and 3-d order distribution r found that search for the physical mechanisms involved functions respectively depending on space coordinates a in the process. One should also mention the existing r1,r2,r3, Φ(ri rj )⁀is the potential of interaction be- | − | experimental studies of the radiation influence on liquid tweeni-thandj-thparticles,ρ= N ⁀isthenumericalden- V matter [26–30]. Those experiments show the strong de- sity. The importantthingto mentionis the factthatthe pendence of the thermodynamic properties of liquids on introducedparameterTeff was,actually,theanalogueof the radiation. Therefore, accounting for the importance the thermodynamic temperature for the nonequilibrium ofunderstandingthebehaviorofliquidsunderirradiation system in focus. It was suggested to be equal to the for the different applications it seems to be attractive to thermodynamic temperature of the corresponding equi- study the physical nature of the interaction. librium system with the same structural properties. It wasshownthatthe“effectivetemperature”couldbecal- The aim of this work is to use the molecular dynam- culatedfromtheperturbedmomentumdistributionfunc- ics (MD) methods to study the physical mechanisms re- tions of the systems caused by the momentum exchange 2 betweentheactiveparticlesandtheparticlesformingthe TABLE II. Masses and intermolecular potential parameters liquid. of water and He[39–41] kT ∂f2(p1,p2)dp dp = p1f (p ,p )dp dp . Atom m,a.u.m. q,e ε, kJ·mol−1 σ, ˚A eff ∂p1 1 2 − m 2 1 2 1 2 H 1 +0.41 0 0 Z Z (2) O 16 -0.82 0.65 3.1656 Therefore, the developed theoretical model of the pro- He 4.0026 +2 0.084 2.6 cess[32]allowedcalculatingthermophysicalpropertiesof the stationary nonequilibrium liquid systems under irra- diation from the thermophysical properties of the corre- order to reach the stationary state of the system we add spondentequilibriumliquid systemsconsideringthe new energyinadiscretewaywiththestep0.05keV.Itisdone parameter T . byacceleratingHeevery2pstohavethetotalirradiation eff In this work we present the results of the MD simu- energy ranging from 0.05 keV to 0.25 keV. lation of the radiation influence on water and compare To justify the choice of the NVE ensemble one can them with the predictions of the theoretical model sug- mention the fact that in all the other ensembles inclu- gested earlier. sion of the thermostats makes the momentum distri- The MD simulations using the rigid simple point bution function move towards the Maxwellian one and, charge (SPC) model of water [34] are performed for a hence, disturbs the real physical picture in the momen- cubic simulation box containing 16384 water molecules tumspace. Therefore,theyarnotacceptableforourcase and 1 He particle at the temperature 300 K, and at am- asoneexpectthe changesinthemomentumdistribution bient pressure (0.1 MPa). The geometric parameters of functiontobeoneofthemainphysicalmechanismsofthe the SPC model used are given in Tab I. The DL POLY radiation influence on he structural and thermodynamic propertiesofwater. Alltheparametersofthe simulation process are given in Tab III. TABLEI.GeometricparametersoftherigidSPCmodelused [35] θo 109.47 TABLE III. Simulation parameters 6 HOH rOH, ˚A 1 Numberof water molecules 16384 rHH,˚A 1.63 Box size,˚A 79 Density ρ, kg·m−3 [42] 992,27 Temperature T, K 300 package (4.06 version) is used [36], with the Ewald par- Equilibration time (NVT ensemble), ps 2000 ticle mesh method for the evaluation of the Coulombin- Simulation time (NVE ensemble), ps 2-10 teractions. The intermolecular potential energy between R , ˚A 14 cutoff atomic sites is calculated in a standard way by a sum Boundary conditions periodic of the Lennard-Jones (12-6) potential and the Coulomb Simulation software DL-POLY4.06 electrostatic interaction TheSPC rigidmodelwasusedinoursimulationsasit U(rij)=ULJ(rij)+UCoul(rij) givesstructuralanddynamic parametersofthe bulk wa- 12 6 terthatareinaccordwiththeexperimentaldata[39,43] σ σ q q ij ij i j = 4εij + (3) and at the same time it gives the results for the radial i≺j (cid:18)rij(cid:19) −(cid:18)rij(cid:19) ! i≺j 4πε0rij distribution functions (RDF) and momentum distribu- X X tion functions that are not blurred by inclusion of the with the off-diagonal interaction parameters calculated internal degrees of freedom. It seems to be attractive as using the Lorentz-Berthelot mixing rules [37, 38] we suggest that one of the main mechanisms responsi- σ +σ ble for the changes of the thermodynamic properties of i j σ = , (4) ij 2 the liquid systems under irradiation is the change in the εij =√εiεj. (5) momentum distribution function [32]. Therefore, using SPC model seems to be a reasonable choice in the first The parametersusedaregiveninTab. II.Allthe sim- approximation. ulations are performed in two steps. First, the system is Inordertohavethereliabledataforthenoneqilibrium equilibratedintheNVT ensembleusingBerendsenther- system in the stationary state we have averaged the re- mostat and the box is fixed. In all the future simulation sults of the three runs with the same initial equilibrated the box size remains unchanged. Next, the radiation is system. The g , g , and g RDFs are shown at OO OH HH switched on. At this stage the simulation is performed Figs. 1, 2, and 3 respectively. in the NVE ensemble. Irradiation is included by accel- It can be seen from the figures that all the changes in erating the He atom originally present in the system. In the RDFs are quite small. It can be explained by the 3 0keV 1,4 0,05keV 3,0 1,2 0,1keV g(r)2,5 g(r) 1,2 00,,125kekeVV 2,5 1,0 0,25keV 1,0 1,3 g(r)OO12,,50 2,02,6 2r ,(8A) 3,0 3,2 r3 (,A6) 4,0 g(r)HH 0,8 0,88 g(r)11,,12 0,6 0keV g(r)0,84 1,0 1,0 0,05keV 0,1keV 0,4 0,80 0,9 0,15keV 2,2 2,4 2,6 2,8 0,5 0,2keV 0,2 0,76 r (A) 0,25keV 2,8 3,0 3,2 0,0 0,0 r (A) 0 2 4 6 8 10 12 14 16 0 2 4 6 8 10 12 14 16 r (A) r (A) FIG.1. Oxygen-oxygenRDFsfortheradiationenergiesrang- FIG. 3. Hydrogen-hydrogen RDFs for the radiation energies ing from 0 KeV to0.25 KeV ranging from 0 KeV to 0.25 KeV 1,5 seen (Fig.3). Still, the tendency of the first peak to de- 0keV crease with the increasing radiation energy is also ob- 1,6 g(r) 00,,015kekeVV served at the graph. 0,15keV 1,4 0,2keV To find the physical mechanism of the radiation influ- 0,25keV ence on water causing the observed structural changes 1,2 1,01,6 1,7 1,8 1,9 2,0 2,1 the obtainedsimulation results are comparedto the the- r (A) g(r)OH 1,0 0,6 oorredteircatol mcaolcduellaotef tthhee“perffoeccetsisvedteevmelpoepreadtueraersl”ietrh[a3t2c].haIrn- 0,8 acterize the structural and thermophysical properties of 0,6 g(r)0,4 the stationary nonequilibrium liquid system under irra- diation the momentum distribution functions have been 0,4 extracted from the simulation data for the different ir- 0,2 0,2 2,2 2,4 2,6 2,8 radiation energies. The obtained curves approximated 0,0 r (A) withthe Maxwelltype functions f(p)=Aexp( φp2)are 0 2 4 6 8 10 12 14 16 shown at Fig. 4. One can see that the positi−on of the r (A) maximumshiftstowardthehighervelocitiesandthepeak FIG. 2. Oxygen-hydrogen RDFs for the radiation energies broadens with the increasing radiation energy. Such a ranging from 0 KeV to 0.25 KeV behavior qualitatively confirms our hypothesis that one ofthe important mechanismsof the irradiationinfluence on the properties of water is the change of the velocity averaging over the three different runs. Still, it can be distributionfunctionduetothemomentumexchangebe- easily seen that there is the correlation in between the tween the active particles and the particles forming the radiation energy and the changes in the RDFs. At the liquid. inset of Fig. 1 one can see the decrease of the height of In order to find the thermophysical properties of the the first peak with the increasing radiationenergy. Such system we calculate the “effective temperature” T eff a behavior suggests the tendency to blurring of the first from the equation (2) with the correspondent momen- and second coordination spheres. At the same time its tum distribution functions f( p ) shown at Fig. 4. The position remains the same at about 2.78 ˚A that corre- resulting values are givenin T|ab.|IV As our systemis in sponds to the geometric criteria for the hydrogen bond that is R 3.3˚A [21]. OO ≤ TABLE IV. Effective temperatures for the different irradia- For g (r) the decrease of the height of the first peak OH tion energies with the increasing radiation energy(the inset at Fig.2) is seen again. Its position at 1.8 ˚A corresponds to the Irradiation energy,keV 0 0.05 0.1 0.15 0.2 0.25 geometric criteria for the hydrogen bond that is ROH Teff, K 300.2 303.4 308.3 313.5 315.2 322.8 2.6˚A [21]. Such a behavior suggests the tendency of th≤e Tirrad, K 300 312 324 335 347 359 radiation to destroy the net of hydrogen bonds. Unfortunately, the changes in the g (r) are hardly the stationary nonequilibrium state it is not possible to HH 4 SPC 0,00077 SPC/E TIP4P 10 TIP4P/2005 Experiment This work Teff 0,00073 s) This work Tirrad 2)/ ^ m 9)( 0^(- 5 0,00069 1,5·10-20 1,7·10-20 1,9·10-20 D (1 0,0007 0,0006 0,00020 0 260 280 300 320 340 360 380 0,0005 T (K) y sit FIG. 5. Selfdiffusion coefficient dependence on temperature: n de 0,0004 0,00005 Datafrom literature for thedifferent models and experiment y [44]; Dependence obtained in this work from T (circles) babilit 0,0003 3,0·10-20 4,0·10-20 and from Tirrad (diamonds). eff o Pr 0,0002 tedagainsttheT areinaccordwiththedataavailable eff in the literature for the SPC rigid model. At the same 0,0001 time calculating the “real” temperature for the system in the nonequilibrium state is not possible and plotting 0 D against Tirrad (Fig. 5) doesn’t show correspondence 0 1,0·10-20 3,0·10-20 5,0·10-20 to the existing data. Therefore, one can come to the p (kg·m·s-1) conclusionthatitisthe“effective”temperaturethatcan explain the observed behavior of the selfdiffusion coeffi- FIG. 4. Momentum distribution functions for the radiation cient. energies ranging from 0 KeV to 0.25 KeV In this work the results of the MD simulations of the α-particles radiation influence on water are presented. It is shown that irradiation causes the changes in the define the “real” temperature of the system. Therefore, structural and thermodynamic properties of water. The to havesomedata forcomparisonin the lastrawofTab. main observed changes in the structure are the blurring IV there are the temperatures that should have the sys- of the first and the second coordination spheres as well tem if all the radiation energy goes for heating only. It as destruction of the net of hydrogen bonds. canbe seenthatthe numbersaredifferent. Atthispoint it is important to mention that the RDFs of the equi- The comparison of the obtained results with the pre- librium system calculated in the NVT ensemble at real dictions of the earlier developed theoretical model of temperatures equal to T given in Tab. IV coincide the process based on the Bogolyubov chain of equations eff with the RDFs shown at Figs. 1-3 for the corresponding showsthatthe observedchangesin structurearecharac- irradiation energies. terizedbyanewparameterthatisthe“effectivetemper- To check the hypothesis that it is the effective tem- ature”. It means that the structure of the liquid system perature but not the “real” or T that describes the under irradiation in the stationary nonequilibrium state irrad thermophysical properties of the system under the irra- isthe sameasthestructureofthecorrespondingequilib- diation it is intereesting to compare simulation results rium system with its thermodynamic temperature equal withtheexperimentaldata. Theeffectsmightbeseenin to the “effective temperature”. Therefore, the “effective allthe thermophysicalproperties thataredefined by the temperature” might be treated as the “structural tem- structure of the system like surface tension coefficient, perature”ofthenonequilibriumsysteminthestationary selfdiffusioncoefficient, etc.. Wehaveanalyzedtheselfd- state. iffusion coefficient dependence on the effective tempera- The changes of the selfdiffusion coefficient under irra- ture. The corresponding graph is shown at Fig. (5 ). At diation are shown to be caused by the changes in the thesamefigurethereareshowngraphsoftheselfdiffusion structure of the liquid system that are described by the coefficient dependence on real temperature that can be growth of the “effective temperature”. The dependence foundintheliterature[44]. OnecanseefromFig. 5that of the selfdiffusion coefficient on the “effective temper- the obtained results for the selfdiffusion coefficient plot- ature” shows the good correspondence to the existing 5 data for the temperature dependence of the selfdiffusion [12] E. Alizadeh, A. G. Sanz, G. Garcia, and L. Sanche, coefficient. Fromthe obtained results it is seen that self- Phys. Chem. Lett. 4 (5), 820 (2013). diffusion coefficient grows with the increasing radiation [13] M. Spotheim-Maurizot and M. Davidkova, J. Phys.: Conf. Ser. 261, 012010 (2011). energy. At the same time the changes in the “real tem- [14] H. Christie, M. Robinson, D. Roach, D. Ross, I. Suarez- perature ” of the nonequilibrium system can not explain Martinez, and N. Marks, Carbon 81, 105114 (2015). the observed growth of the selfdiffusion coefficient in a [15] G. R. Lumpkin, K. L. Smith, M. G. Blackford, B. S. correctway. Therefore,itisshownthatthe thermophys- Thomas,K.R.Whittle,N.A.Marks, andN.J.Zaluzec, ical properties of the liquid system under irradiation in Phys. Rev.B 77, 214201 (2008). the stationary nonequilibrium state are the same as the [16] K. Trachenko, M. T. Dove, E. Artacho, I. T. Todorov, thermophysical properties of the corresponding equilib- and W. Smith,Phys. Rev.B 73, 174207 (2006). [17] L. Xu, S. V. Buldyrev, H. E. Stanley, and G. Franzese, rium system with its thermodynamic temperature equal Phys. Rev.Lett. 109, 095702 (2012). to the “effective temperature”. [18] E. A. Elfimova, A. O. Ivanov, and P. J. Camp, It is shown that the changes in the coefficients of the Phys. Rev.E 88, 042310 (2013). Maxwelldistribution function due to the momentum ex- [19] S. Toxvaerd,Phys. Rev.E 58, 704 (1998). change between the active particles and the particles [20] D.vanderSpoel,P.J.vanMaaren, andH.J.C.Berend- formingtheliquidisoneoftheimportantphysicalmech- sen, J. Chem. Phys. 108, 10220 (1998). anisms of the radiation influence on liquid systems. The [21] G.W.Robinson,S.B.Zhu,S.Singh, andM.W.Evans, Water in Biology, Chemistry and Physics: Experimental knowledgeofthe distortedmomentumdistributionfunc- Overviews and Computational Methodologies, WorldSci- tionoftheliquidsystemunderirradiationinthestation- entific Series in Contemporary Chemical Physics, Vol. 9 arynonequilibriumstateallowscalculatingthe“effective (World Scientific, Singapore, 1996). temperature” of such a system and, hence, its structural [22] M.P.AllenandD.J.Tildesley, Computer Simulationof and thermodynamic properties. Liquids (Oxford UniversityPress Inc., 1987). TheobtainedintheworkresultsoftheMDsimulation [23] I. Draganic, Rad. Phys.Chem. 72, 181 (2005). oftheradiationinfluenceonwaterquantitativelyconfirm [24] W. G. Burns, Nature 339, 515 (1989). [25] H. Sims, Rad. Phys.Chem. 75, 1047 (2006). the predictions of the previously introduced theoretical [26] Y. Kolesnichenko, Nuclear Fusion 15, 35 (1975). modeloftheprocessandareinaccordwiththe available [27] W. Martino, F. de la Mora, and Y. Yoshida, Green in the literature experimental data. Chemistry 8, 390 (2006). [28] M. Zenkievicz, J. Achiev. Mat. Manufact. Eng. 2, 43 (2007). [29] M. Byung and H. Jung, Appl. Phys. Lett. 93, 244105 (2008). ∗ [email protected] [30] B. M. Weon, J. H. Je, Y. Hwu, and G. Margaritondo, [1] E. Zarkadoula, S. L. Daraszewicz, D. M. Phys. Rev.Lett. 100, 217403 (2008). Duffy, M. A. Seaton, I. T. Todorov, K. Nord- [31] N. Bogolyubov, “Studies in statistical mechanics,” lund, M. T. Dove, and K. Trachenko, (North-Holland,1962)Chap.Problemsofdynamicalthe- J. Phys.: Condens. Matter 25, 125402 (2013). ory in statistical physics, p. 5. [2] L. Yuan, J. Peng, L. Xu, M. Zhai, J. Li, and G. Wei, [32] L. A. Bulavin, K. V. Cherevko, D. A. Gavryushenko, J. Phys. Chem. B 113, 89488952 (2009). V. M. Sysoev, and T. S. Vlasenko, Phys. Rev. E 93, [3] H. Luna and E. C. Montenegro, 032133 (2016). Phys.Rev.Lett. 94, 043201 (2005). [33] H. N. V. Temperley, J. S. Rowlinson, and G. S. Rush- [4] W.J.Weber,R.C. Ewing, C.R.A.Catlow, T. D.dela brooke, eds., Physics of Simple Liquids (North-Holland Rubia, L. W. Hobbs, C. Kinoshita, H. Matzke, A. T. Publishing Company,Amsterdam, 1968). Motta, M. Nastasi, E. K. H. Salje, E. R. Vance, and [34] H.Berendsen,J.Postma,W.vanGunsteren, andJ.Her- S.J. Zinkle, JMR 13, 1434 (1998). mans,“Interactionmodelsforwaterinrelationtoprotein [5] E. M. Fielden and P. O’Neill, eds., hydration,” (Reidel, Dordrecht,1981) pp.331–342. The Early Effects of Radiation on DNA, NATO ASI [35] Y. Wu, H. L. Tepper, and G. A. Voth, Series H:Cell Biology, Vol.54 (Springer-Verlag, 1991). J. Chem. Phys. 124, 024503 (2006). [6] G. O. Phillips, G. J. Moody, and G. L. Mattok, [36] I. Todorov, W. Smith, K. Trachenko, and M. Dove, J. Chem. Soc. ,3522 (1958). J. Mater. Chem. 16, 1911 (2006). [7] K. Trachenko, J. M. Pruneda, E. Artacho, and M. T. [37] H. A.Lorentz, Ann.Phys. 248, 127136 (1881). Dove,Phys. Rev.B 71, 184104 (2005). [38] D. Berthelot, C. R. Acad.Sci. 126, 1703 1855 (1898). [8] K. Trachenko, M. T. Dove, and E. K. H. Salje, [39] J. Zielkiewicz, J. Chem. Phys. 123, 104501 (2005). Phys.Rev.B 65, 180102(R) (2002). [40] A.P.LyubartsevandA.Laaksonen,J.Phys.Chem.100, [9] K. Nordlund, M. Ghaly, R. S. Averback, 16410 (1996). M. Caturla, T. Diaz de la Rubia, and J. Tarus, [41] H. Y. Tang, , and I. J. Ford, Phys.Rev.B 57, 7556 (1998). J. Chem. Phys. 125, 144316 (2006). [10] L. Yuan, J. Peng, M. Zhai, J. Li, and G. Wei, [42] E. Lemmon, M. Huber, and M. McLinden, REFPROP: Rad.Phys. Chem. 78, 737739 (2009). Reference fluid thermodynamic and transport properties. [11] M. Qi, G. Wua, Q. Li, and Y. Luo, Version 8.0, NIST standard reference database (NIST, Rad.Phys. Chem. 77, 877883 (2008). Gaithersburg, 2007). 6 [43] B. Guillot, J. Mol. Liq.101, 219 (2002). [44] G. Guevara-Carrion, J. Vrabec, and H. Hasse, J. Chem. Phys. 134, 074508 (2011).