Combined Molecular and Spin Dynamics Simulation of BCC Iron with Vacancy Defects
Mark Mudrick, Markus Eisenbach, Dilina Perera, David P. Landau
aa r X i v : . [ phy s i c s . c o m p - ph ] M a r Combined Molecular and Spin Dynamics Simulation of BCC Iron with VacancyDefects
Mark Mudrick, ∗ Markus Eisenbach, Dilina Perera, and David P. Landau Center for Simulational Physics, The University of Georgia, Athens, Georgia 30602, USA Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831, USA Department of Physics and Astronomy, Texas A&M University, College Station, Texas 77843, USA (Dated: September 21, 2018)Utilizing an atomistic computational model which handles both translational and spin degreesof freedom, combined molecular and spin dynamics simulations have been performed to investigatethe effect of vacancy defects on spin wave excitations in ferromagnetic iron. Fourier transformsof space and time-displaced correlation functions yield the dynamic structure factor, providingcharacteristic frequencies and lifetimes of the spin wave modes. Comparison of the system with a5% vacancy concentration with pure lattice data shows a decrease in frequency as well as a decreasein lifetime for all transverse spin wave excitations observed. Additionally, a rugged spin wave lineshape for low- q spin waves indicates the presence of multiple localized excitations near defect sitesresulting in reduced excitation lifetimes due to increased magnon-magnon scattering. We observefurther evidence of increased magnon-magnon scattering as the peaks in the longitudinal spin wavespectrum become less distinct. I. INTRODUCTION
Magnetic properties of ferromagnetic metals at varyinglevels of impurity have previously been investigated usingneutron scattering techniques . Notably, the introduc-tion of substitutional defects into a magnetic system hasbeen shown to distort the spectrum of the system’s char-acteristic excitations . This effect has been investigatedboth theoretically as well as experimentally for lat-tice vibrations and spin waves. These experiments haveprovided spin wave dispersion curves and stiffness con-stants. In particular, the spin wave stiffness parameter D has shown sensitivity to variations in the concentrationof magnetic defects . Additionally, both spin wave ener-gies and lifetimes decrease as the impurity concentrationgrows.Computational techniques have long been utilized toinvestigate collective excitations in systems spanning abroad range of practical interest. Molecular dynamics(MD) methods have been used to model vibrationalbehavior in metals, alloys, biological systems, etc. whilespin dynamics (SD) simulations have investigated mag-netic properties of classical lattice-based spin models .SD simulations have replicated experimental findings insimple spin systems such as RbMnF and have success-fully predicted the existence of two-spin-wave modes .With the use of coordinate dependent exchange inter-actions, this method has also had success investigatingsystems with more complex magnetic properties such aspropagating spin wave modes and external magneticfield effects in bcc iron .Each of these methods numerically solves the classi-cal equations of motion which describe the dynamicalevolution of the model under consideration. Despite theindividual capabilities of MD and SD simulation meth-ods, the coupling of lattice and spin degrees of freedomis inherently neglected by both. In magnetic materials, atomic motion affects magnetic moments which dependnon-trivially on the local atomic environment. Likewise,magnetic interactions have been shown to contribute tostructural properties of these materials, including theBCC structure of iron . Therefore, these degreesof freedom, atomic and magnetic, should be consideredjointly in any model aiming to investigate the excitationsof such a system.The combined molecular dynamics - spin dynamics(MD-SD) method treats the spin subsystem using anextension of the Heisenberg model where the exchangeinteraction is a coordinate-dependent, pairwise functionof inter-atomic distance . Inter-atomic interactions arehandled with a non-magnetic many-body potential. Thetime integration of the coupled equations of motion is cal-culated using an algorithm based on the Suzuki-Trotterdecomposition of the exponential time evolution oper-ator. This algorithm is time reversible, efficient, andis known to conserve phase-space volume . To han-dle systems of practical interest, such a model must uti-lize empirical many-body inter-atomic potentials and ex-change interactions parameterized by experimental dataand first-principles calculations . Using this compress-ible magnetic model parameterized for BCC iron, the in-terplay between spin waves (magnons) and lattice vibra-tions (phonons) has previously been investigated, show-ing spin wave dampening as well as magnon-phonon cou-pling in pure systems .This paper will extend the current MD-SD frameworkto investigate BCC iron with vacancies. We show theeffect of these impurities on magnon dispersion curves aswell as on individual spin wave excitation peaks. Thespin waves are dampened and broadened in the presenceof vacancy defects, in agreement with experimental neu-tron scattering data. II. MODELA. MD-SD Algorithm
The MD-SD method extends the traditional MD ap-proach through the addition of a third phase variable- the spin angular momentum S i - to the position andvelocity degrees of freedom. This spin variable is incor-porated into the MD-SD Hamiltonian such that H = N X i =1 mv i U ( { r i } ) − X i < j J ij ( { r k } ) S i · S j (1)for a system of N magnetic atoms of mass m withpositions { r i } , velocities { v i } , and atomic spins { S i } .The first term in the MD-SD Hamiltonian representsthe kinetic energy of the atoms, and U ( { r i } ) is a non-magnetic interatomic potential. The third term describesa Heisenberg-like exchange interaction which includesa coordinate-dependent exchange parameter J ij ( { r k } ).The introduction of coordinate dependence into J ij ( { r k } )allows for the exchange interaction to depend on the lo-cations of all atoms in the vicinity of atoms i and j . Thisframework is general and may be utilized for any mag-netic material given proper parameterization of the inter-atomic potential and exchange parameter. To investigateBCC iron, we employ the embedded atom interatomicpotential U ( { r i } ) developed and parameterized for ironby Dudarev and Derlet . The exchange interaction usedin this work is a pairwise function J ( r ij ) parameterizedby first principles calculations for iron .This MD-SD Hamiltonian has true dynamics accordingto the classical equations of motion d r i dt = v i (2a) d v i dt = f i m (2b) d S i dt = 1¯ h H effi × S i (2c)where f i and H effi represent the interatomic force andeffective field acting on the i th atom, respectively.In order to investigate the collective excitations of thesimulated system, space-displaced, time-displaced corre-lation functions are computed. By performing Fouriertransformations of these correlation functions, we ob-tain information about the frequency and lifetimes ofthese excitations. During a simulation run, the spatialFourier transform of the space-displaced, time-displacedspin-spin correlation function, or the intermediate scat-tering function, is calculated on the fly F kss ( q , t ) = 1 N h ρ ks ( q , t ) ρ ks ( − q , i , (3) where k represents the real-space directions { x, y, z } and ρ s ( q , t ) represents the microscopic spin density, definedas ρ s ( q , t ) = X i S i ( t ) e − i q · r i ( t ) . (4)During a microcanonical simulation of a ferromagneticmaterial, the total magnetization vector is a constant ofmotion. This property allows for the differentiation ofspin excitations which propagate parallel and perpendic-ular to the magnetic symmetry axis through the choiceof a Cartesian coordinate system such that the z -axis isparallel to the magnetization vector. The components of F ss ( q , t ) are then regrouped into a longitudinal compo-nent F Lss ( q , t ) = F zss ( q , t ) (5)and a transverse component F Tss ( q , t ) = 12 [ F xss ( q , t ) + F yss ( q , t )] . (6)A temporal Fourier transform of the intermediate scatter-ing function yields the spin-spin dynamic structure factor S L,Tss ( q , ω ) = 12 π Z + ∞−∞ F L,Tss ( q , t ) e − iωt dt (7)for momentum transfer q and frequency ω . The dynamicstructure factor obtained from MD-SD simulations maybe compared to the measurements made by inelastic neu-tron scattering experiments . B. Simulation Details
Time integration of the equations of motion shown inEq. (2) is performed using a second-order Suzuki-Trotterdecomposition of the time evolution operator . Weequilibrate the position and spin subsystems of multipleinitial configurations to the desired temperature T usingthe Metropolis Monte Carlo method. Initial velocities arethen drawn from the Maxwell-Boltzmann distribution atthis temperature. A short microcanonical MD-SD run isperformed in order to equilibrate the three sub-systems.Once properly equilibrated, these configurations are inde-pendently time-evolved using the microcanonical MD-SDintegration scheme. For this time integration, a time-stepof δt = 1 fs was found to adequately conserve energy andmagnetization. Simulation data shown in this paper arefrom MD-SD runs of 1 ns (10 time-steps) in duration.During the time integration, the intermediate scatter-ing function is recorded on the fly for each independentconfiguration. Once all runs have completed we averagethese data, yielding an estimate of the canonical ensem-ble average at temperature T .In order to include vacancy type defects, we leave ran-domly selected sites devoid of an atom in each initial ω (meV) S ss T ( q , ω ) Pure Lattice5% Vacant Sites q = (0.137 Å -1 , 0, 0) (a)
75 85 95 10500.010.020.03 S ss T ( q , ω ) Pure Lattice5% Vacant Sites (b) q = (0.581 Å -1 , 0.581 Å -1 , 0)
300 350 400 450 500 ω (meV) S ss T ( q , ω ) Pure Lattice5% Vacant Sites (c) q = (3.32Å -1 , 3.32 Å -1 , 3.32 Å -1 ) FIG. 1. Transverse spin-spin dynamic structure fac-tor obtained from MD-SD simulations for L = 16 at T = 300 K for both the pure lattice and the sys-tem with 5% of lattice sites left vacant. (a) q =(0 . − , , q = (0 . − , . − , q =(3 . − , . − , . − ). Symbols represent simulationdata while solid lines show multi-peak Lorentzian curve-fitting. Dashed lines indicate the constituent peaks whichmake up the solid curve for the systems with 5% vacancy.
40 42 44 ω (meV) S ss T ( q , ω ) Pure Lattice2 Distant Vacancies4 Distant Vacancies
FIG. 2. Transverse spin-spin dynamic structure factor ob-tained from MD-SD simulations for L = 6 at T = 300 K for q = (0 . − , , configuration. For the data presented here, simulationshave been performed for a body-centered cubic system ofedge length L = 16 (8192 sites) with periodic boundaryconditions at T = 300 K. For a 5% vacancy concentra-tion, 7783 sites contain an atom while the remaining 409sites remain vacant. Each independent configuration in-cluded in the overall simulation set has a unique vacancydistribution.Using the canonical ensemble average estimate for thespin-spin intermediate scattering function obtained fromthese simulations, we calculate the dynamic structurefactor. Spin wave excitations are observed through thetransverse spin-spin dynamic structure factor S Tss ( q , ω )as peaks of Lorentzian form, S Tss ( q , ω ) = I Γ ( ω − ω ) + Γ (8)where ω is the characteristic frequency of the excitation, I is the amplitude, and Γ is the half width at half max-imum. Γ is inversely proportional to the lifetime of theassociated excitation. Fitting simulational data to Eq.8 allows us to observe the impact of vacancy defects onspin wave dispersion curves, as well as on individual spinwave line shapes, excitation energies, and lifetimes. III. RESULTSA. Transverse Spin Waves
1. Spin Wave Energies
In order to demonstrate our approach, we performedcalculations with a 5% vacancy concentration, a figureconsistent with the nonmagnetic defect concentration q ( 2 π / a ) ω ( m e V ) Pure Lattice5% Vacancy Γ H P Γ N[q 0 0] [q q q] [q q 0] q ( 2 π / a ) ( ω P u r e - ω % ) / ω P u r e (b) [q 0 0] [q q q] [q q 0] FIG. 3. (a) Spin wave dispersion curve at T = 300 K ob-tained from MD-SD simulations with L = 16. Results areshown for the pure lattice as well as the system with 5% va-cancy concentration. (b) The fractional shift in transversespin wave frequencies due to a 5% concentration of vacancydefects. Results obtained from MD-SD simulations of systemsof size L = 16 at T = 300 K. found in experimental studies. Figure 1 shows the ef-fect of this 5% vacancy concentration on three differenttransverse spin waves in a system of size L = 16. The va-cancy defects cause a dampening of the peaks in S Tss ( q , ω )as well as a shift to lower frequency. These effects on spinwave peaks are observed for all wave vectors, though weonly show a few selected line shapes here. Visible inFig. 1(a), the low- q peak region of the impure systemdisplays a more rugged spectral structure than the purepeaks or any of the higher q peaks from impure systems.While curve fitting these rough peaks to the form of Eq.(8) yields estimates of peak locations and half-widths,applying a fit utilizing the sum of multiple Lorentzianforms more accurately reproduces the asymmetries andadditional structure in these peaks. Included in Fig. 1 are fits to two-peak (Figs. 1(a), 1(b)) and four-peak (Fig.1(c)) Lorentzian functions, describing the distorted sig-nals obtained from the impure system. Multiple peaksare also visible in the small- q region in the ( q, q,
0) and( q, q, q ) symmetry directions, though we only show the( q, ,
0) peak in Fig. 1. The distortion of the these spinwave line shapes is due to the presence of localized spinwave modes near the defect sites. While these multi-peakLorentzian functions provide reasonable fits to our data,the possibility remains of additional structure within thelimits of our resolution.In order to investigate the downward shift of spin wavefrequencies seen in Fig. 1, we simulated a small system( L = 6) from which atoms are removed one by one. Thedecrease in spin wave frequency due to this removal ofatoms is observed in Figure 2. The atoms removed fromthe configurations in Fig. 2 are chosen such that thesites are separated by a distance greater than the cutoffdistance of the interatomic potential. Multiple vacancyconfigurations were generated, though only one configu-ration is shown for each data set in Fig. 2. As the numberof vacancy sites grows, the frequency of the characteristicspin wave spectrum decreases.The full spin wave dispersion curve is shown in Fig-ure 3(a), constructed using the peak locations obtainedthrough Lorentzian curve fitting. To observe the effect ofvacancy defects on spin wave excitations more carefully,we compare the characteristic spin wave frequencies ob-tained from MD-SD simulations for the pure system withthose from systems with a 5% vacancy concentration.We calculate the fractional shift in spin wave frequency,( ω Pure − ω /ω Pure ), shown in Fig. 3(b) directly belowthe dispersion curve. All observed spin wave excitationsdisplay a shift to lower frequency, with this effect beingmost significant near the zone center. The shaded areasin Fig. 3(b) indicate the regions in the small- q S Tss ( q , ω )spectra where we observe multiple peaks. While multi-peak Lorentzian forms fit the high- q signals more accu-rately than the single peak form, identification of theless intense constituent peaks is unreliable. Therefore weinclude only the most intense peak signal from our curve-fitting procedure for high- q excitations in Fig. 3(b).In ferromagnetic materials below the critical tempera-ture T c , spin waves with small- q are isotropic and are ex-pected to approximate a quadratic dispersion relation of the form ¯ hω = D | q | (1 − β | q | ) (9)where D represents the magnetic stiffness constant. Thedispersion curve for low- q is shown in Figure 4 for boththe pure system and that with 5% vacancy concentra-tion. The dispersion parameters D and β in Eq. (9)are quantities accessible through neutron scattering us-ing diffraction methods (DM), triple axis spectrometry(TAS), chopper spectrometry (CS), or other methods.While the dispersion relation shown in Eq. (9) providesan accurate fit to both experimental and computationaldata, the behavior of the parameter β as observed in | q | ( Å -1 ) ω ( m e V ) Pure - [1 0 0]Pure - [1 1 0]Pure - [1 1 1]5% Vacancy - [1 0 0]5% Vacancy - [1 1 0]5% Vacancy - [1 1 1]
FIG. 4. Comparison of low | q | magnon dispersion curves ob-tained from MD-SD simulation (L=16) for the pure systemand that with 5% vacancy defects. Lines in the figure rep-resent a two-parameter curve fit to Eq. 9 with D shown inTable I.TABLE I. Experimental and computational comparison ofspin wave dispersion parameters in impure Fe systems at roomtemperature. D (meV ˚A ) Method Ref.Fe 307 ±
15 CS Loong et al. Fe 281 ±
10 TAS Shirane et al. Fe (4% Si) 270 ± ± ± ± et al. Fe (15% Si) 233 ± experiment is inconsistent. However previous investiga-tion of BCC iron under varying levels of impurity hasshown systematic behavior in the stiffness parameter D .The conditions most notable to this study are the sub-stitutional inclusion of Si, a nonmagnetic atom, into thecrystal. Due to the lack of a magnetic moment and lowermass in comparison to Fe, the Si atom mimics a magneticvacancy. The lines in Fig. 4 represent curve fits to Eq.(9), and the values of the stiffness fit parameter D areshown along with previous experimental data in Table I.While our results overestimate the value of the stiffnessparameter, they capture the trend of a decrease in D caused by nonmagnetic defects observed in experimentaldata.
2. Spin Wave Lifetimes
Spin wave excitation lifetimes are inversely propor-tional to the half width at half maximum of the charac- ( Γ P u r e - Γ % ) / Γ P u r e q ( 2 π / a ) [ q 0 0 ][ q q 0 ][ q q q ] FIG. 5. Fractional shift in half width at half maximum (Γ)of transverse spin waves at T = 300 K obtained from MD-SDsimulations for L = 16 in the [1 0 0], [1 1 0], and [1 1 1]symmetry directions.
14 16 18 ω (meV) S ss L ( q , ω ) Pure Lattice5% Vacancy (2, 0, 0) - (1, 0, 0)
FIG. 6. Longitudinal spin-spin dynamic structure factor ob-tained from MD-SD simulations for q = (0 . − , ,
0) at T = 300 K with L = 16 for both the pure system and that with5% of sites left vacant. Solid lines represent the predictedpeak locations in the pure system, and dashed lines representthe predicted locations of the (2 , , − (1 , ,
0) peaks in thesystem with vacancies. For clarity, only the (2 , , − (1 , , teristic peak in S Tss ( q , ω ) which may be obtained throughcurve fitting to Eq. (8). The effect of vacancy defects onspin wave lifetime is presented in Figure 5, which showsthe fractional change in the half width at half maximumdue to impurities. Fig. 5 shows significant broadeningof spin wave peaks for all q values, indicating a decreasein lifetime for all spin wave excitations. However theeffect is most significant for low q excitations, and thebroadening decreases as q grows larger in all symmetrydirections. This decreased lifetime at low q is due toincreased magnon-magnon scattering caused by the ex-istence of additional spin wave modes, evidenced by therugged structure of low q peaks in S Tss ( q , ω ) as in Fig. 1(a). B. Longitudinal Spin Waves
For classical Heisenberg spin models, Bunker et al. have previously shown that the excitation peaks in thelongitudinal component of S ss ( q , ω ) represent creationand/or annihilation processes resulting from the interac-tion of multiple transverse spin waves . Previous MD-SD simulation of pure BCC iron measured the longitudi-nal component of S ss ( q , ω ), and the frequencies of thesetwo-spin wave peaks were identified using the differenceof transverse spin wave frequencies : ω − ij ( q i ± q j ) = ω ( q i ) − ω ( q j ) (10)where q i and q j are the wave-vectors of the constituentspin waves, and ω is the characteristic frequency of each.Figure 6 shows a portion of the S Lss ( q , ω ) spectrum for q = πLa (1 , ,
0) the L = 16 system. Since the large set ofavailable wave vectors { q i } in an L = 16 system leads tomany two-spin-wave interaction processes, only a smallsection of S Lss ( q , ω ) is presented here. Individual peaksin the pure lattice spectrum are well defined and we haveidentified the two transverse spin waves that combine toproduce each peak in the longitudinal component. An ex-ample of the identification of these peaks is shown by thesolid vertical line in Fig. 6 which represents the differencein frequencies of the πLa (1 , ,
0) and πLa (2 , ,
0) transversespin waves, calculated using Eq. 10. While we have iden-tified the other peaks in the pure lattice data of Fig. 6,only one is shown in the figure for clarity.Fig. 6 also shows the effect of 5% vacancy defectson S Lss ( q , ω ). Qualitatively, the longitudinal spectrumis less well-defined for the impure system, as the indi-vidual excitation peaks seen in the pure system shift tolower frequency and merge into a broad distribution. Asshown previously in Fig. 1(a), small- q spin waves such asthose that contribute to the solid line in Fig. 6 displaymultiple transverse excitation peaks. Therefore to iden-tify the location of longitudinal excitations in the impuresystem, we must consider each combination of the con-stituent transverse spin waves. Both of the πLa (1 , ,
0) and πLa (2 , ,
0) transverse spin wave spectra show dual-peak structure, therefore Eq. 10 has at least four uniquecombinations of spin wave interactions. Each of these re-sultant frequencies are shown as dashed vertical lines inFig. 6, indicating that the loss of individual peak reso-lution in this region of S Lss ( q , ω ) is due to the increase inmultiple spin wave annihilation processes. IV. CONCLUSIONS
The effect of vacancy defects on magnetic excitationsin BCC iron were investigated using combined molec-ular and spin dynamics simulations (MD-SD) at roomtemperature. We calculated the intermediate scatteringfunction on-the-fly during our simulation runs, and usedthis data to obtain the dynamic structure factor S ss ( q , ω )which contains information regarding spin wave energiesand lifetimes. We obtained spin wave energies and halfwidths at half maximum by performing curve-fitting ofthe characteristic spin wave peaks to a Lorentzian func-tion. For all observed spin waves, the introduction ofvacancy defects was shown to decrease the energy of theexcitation with the effect being most significant near thezone center. The S ss ( q , ω ) line shape in low- q spin wavesalso showed rugged structure, indicating the propagationof multiple excitations with different characteristic ener-gies.The presence of multiple spin waves modes increasesmagnon-magnon scattering, which is evident by the de-crease in observed spin wave lifetime for low- q excita-tions. The low q region of the magnon dispersion curvehas been shown experimentally to obey a quadratic dis-persion relation of the form of Eq. (9). From thisquadratic form, the magnetic stiffness constant has beenmeasured experimentally for BCC iron under varyingconditions of magnetic impurity. Our findings using MD-SD simulations are consistent with the observed decreasein the stiffness constant with the introduction of defects,as shown in Fig. 4 and Table I. We investigated longi-tudinal spin wave modes which represent two spin waveinteraction processes. The introduction of defects intothe system resulted in the loss of clearly defined peakswhich are observed in the pure system. We have cal-culated the resultant frequency of interacting transversespin waves, showing that this effect is due to an increasein available spin wave scattering processes. Further quan-titative studies of magnetic materials containing impurityatoms (e.g. Fe − x Si x ) will require reliable interatomicpotentials for those alloys. ACKNOWLEDGMENTS
Part of this work (M.E.) was supported by the U.S. De-partment of Energy, Office of Science, Basic Energy Sci-ences, Material Sciences and Engineering Division. Thiswork was sponsored in part by the Center for DefectPhysics, an Energy Frontier Research Center of the Of-fice of Basic Energy Sciences (BES), U.S. Departmentof Energy (DOE). This study was supported in partby resources and technical expertise from the GeorgiaAdvanced Computing Resource Center, a partnership between the University of Georgia’s Office of the VicePresident for Research and Office of the Vice Presidentfor Information Technology. This research also used re-sources of the Oak Ridge Leadership Computing Facilityat ORNL, which is supported by the Office of Scienceof the U.S. Department of Energy under Contract No.DE-AC05-00OR22725. ∗ [email protected] J. Lynn, Phys. Rev. B , 2624 (1975). H. Mook and R. Nicklow, Phys. Rev. B , 336 (1973). S. Pickart, H. Alperin, V. Minkiewicz, R. Nathans, G. Shi-rane, and O. Steinsvoll, Phys. Rev. , 623 (1967). C.-K. Loong, J. Carpenter, J. Lynn, R. Robinson, andH. Mook, J. Appl. Phys. , 1895 (1984). M. Collins, V. Minkiewicz, R. Nathans, L. Passell, andG. Shirane, Phys. Rev. , 417 (1969). E. Svensson, T. Holden, W. Buyers, R. Cowley, andR. Stevenson, Solid State Comm. , 1693 (1969). S. Takeno, Progress of Th. Phys. , 731 (1963). A. Maradudin, Rep. on Prog. in Phys. , 331 (1965). Y. A. Izyumov and M. Medvedev, Sov. Phys. JETP (1966). B. Antonini, F. Menzinger, A. Paoletti, and A. Tucciarone,Phys. Rev. , 833 (1969). D. C. Rapaport,
The Art of Molecular Dynamics Simu-lation , 2nd ed. (Cambridge University Press, New York,2004). D. P. Landau and M. Krech, J. Phys.: Condens. Matter , R179 (1999). S.-H. Tsai, A. Bunker, and D. P. Landau, Phys. Rev. B , 333 (2000). A. Bunker and D. P. Landau, Phys. Rev. Lett. , 2601(2000). X. Tao, D. Landau, T. Schulthess, and G. Stocks, Phys.Rev. Lett. , 087207 (2005). C. Chui and Y. Zhou, AIP Advances , 037110 (2014). H. Hasegawa and D. G. Pettifor,Phys. Rev. Lett. , 130 (1983). H. C. Herper, E. Hoffmann, and P. Entel,Phys. Rev. B , 3839 (1999). I. P. Omelyan, I. M. Mryglod, and R. Folk, Phys. Rev.Lett. , 898 (2001). M. Krech, A. Bunker, and D. P. Landau, Comp. Phys.Comm. , 1 (1998). P. W. Ma, C. H. Woo, and S. L. Dudarev, Phys. Rev. B , 024434 (2008). D. Perera, D. M. Nicholson, M. Eisenbach, G. M. Stocks,and D. P. Landau, Phys. Rev. B , 014431 (2017). S. L. Dudarev and P. M. Derlet, J. Phys.: Condens. Matter , 7097 (2005). S. W. Lovesey,
Theory of Neutron Scattering from Con-densed Matter (Oxford University Press, Oxford, 1984). S.-H. Tsai, H. Lee, and D. Landau, Am. J. Phys. , 615(2005). G. Shirane, V. Minkiewicz, and R. Nathans, J. Appl. Phys.39