Magnetostatic Dipolar Energy of Large Periodic Ni fcc Nanowires, Slabs and Spheres

Magnetostatic Dipolar Energy of Large Periodic Ni fcc Nanowires, Slabs and Spheres I. Cabriaa,∗ aDepartamento de F´ısica Te´orica, At´omica y ´ Optica, Universidad de Valladolid, 47011 Valladolid, Spain Abstract The computational effort to calculate the magnetostatic dipolar energy, MDE, of a periodic cell of Nmagnetic moments is an O(N2) task. Compared with the calculation of the Exchange and Zeeman energy terms, this is the most computationally expensive part of the atomistic simulations of the magnetic properties of large periodic magnetic systems. Two strategies to reduce the computational effort have been studied: An analysis of the traditional Ewald method to calculate the MDE of periodic systems and parallel calculations. The detailed analysis reveals that, for certain types of periodic systems, there are many matrix elements of the Ewald method identical to another elements, due to some symmetry properties of the periodic systems. Computation timing experiments of the MDE of large periodic Ni fcc nanowires, slabs and spheres, up to 32000 magnetic moments in the periodic cell, have been carried out and they show that the number of matrix elements that should be calculated is approximately equal to N, instead of N2/2, if these symmetries are used, and that the computation time decreases in an important amount. The time complexity of the analysis of the symmetries is O(N3), increasing the time complexity of the traditional Ewald method. MDE is a very small energy and therefore, the usual required precision of the calculation of the MDE is so high, about 10−6 eV/cell, that the calculations of large periodic magnetic systems are very expensive and the use of the symmetries reduces, in practical terms, the computation time of the MDE in a significant amount, in spite of the increase of the time complexity. The second strategy consists on parallel calculations of the MDE without using the symmetries of the periodic systems. The parallel calculations have been compared with serial calculations that use the symmetries. ∗Corresponding author. Tel.: +34 983 423141; Fax: +34 983 423013 Email address: [email protected] (I. Cabria) 1. Introduction and Motivation Magnetic anisotropy is one of the most important properties of magnetic materials, from a scientific and also from a technological point of view. Some of the applications of the magnetic anisotropy are permanent magnets, magnetic memories, electric motors and magnetic field sensors. The magnetic anisotropy energy, MAE, is the energy change due to a change of the magnetization direction. There are two contributions to the magnetic anisotropy energy: The electronic band structure (or simply, electronic) anisotropy energy and the shape anisotropy. The electronic contribution is due to the simultaneous occurrence of the electron relativistic interaction (spin-orbit coupling) and spin-polarization in the electronic structure of the magnetic systems. The magnetostatic anisotropy energy results from the classical magnetic dipolar interactions and therefore, is also called magnetostatic dipolar anisotropy energy, MDAE. Due to the long range character of the magnetic dipolar interactions, the magnetostatic dipolar anisotropy energy depends, in general, on the shape of the magnetic system and hence, the shape anisotropy is usually adscribed to this anisotropy. The magnetostatic dipolar anisotropy energy is zero for cubic systems and negligibly small for weak anisotropic systems such as cobalt. However, for systems with a large anisotropy, such as layered materials and nanowires of ferromagnetic atoms, the magnetostatic dipolar anisotropy energy can not be neglected and is comparable with the electronic band anisotropy energy or even larger [1–6]. In the case of nanowires, the elongated shape enhances the magnetostatic dipolar anisotropy of these materials. A flip of the orientation of the magnetic moments from out-of-plane to in-plane as the number of layers or the thickness of magnetic layered materials increases, is observed in the theoretical calculations [1–10] and in the experiments [11–17]. The electronic or band contribution causes an out-of-plane or perpendicular orientation of the magnetic moments of layered systems, while the magnetic dipolar anisotropy causes and in-plane or parallel orientation. The explanation of the flip of the orientation is that the dipolar interactions increase as the thickness of the layered materials increases and are larger than the spin-orbit interactions, which are surface terms, if the thickness is large enough. For thick layered materials, the preferred orientation will be in-plane and the MDAE will depend linearly on the number of layers. Hence, in the simulations of large thick layered magnetic systems, the approach of considering the MAE composed only by the magnetostatic dipolar anisotropy energy is usually adopted. The calculation of the MAE as the electronic part plus the dipolar-dipolar part is a hybrid, relativistic-classical, approach. A consistent treatment of the MAE should consist on a fully relativistic calculation of the system. Jansen [18, 19] proved that the shape anisotropy is caused by the Breit interaction [20, 21], a relativistic correction of the Coulomb interaction between electrons. Bornemann et al. [22] did fully relativistic band structure calculations of magnetic layered systems, accounting simultaneously for spin-orbit coupling and the Breit interaction. 2 They compared numerical results of the Breit interaction and the magnetostatic dipolar contributions to the MAE and they found that they were very close. The relativistic calculations are computationally much more expensive than the classical dipole-dipole calculations. Therefore, it makes sense, from a practical point of view, to calculate the shape anisotropy as the classical magnetostatic dipolar anisotropy, instead of carrying out relativistic calculations. The most expensive part of the atomistic simulations of the magnetic properties of periodic magnetic systems of certain thickness, such as nanowires and films of ferromagnetic atoms, is the calculation of the magnetostatic dipolar energy, MDE [9, 23–28]. To simulate these materials with realistic models, it is necessary to consider a large number of atoms and the details of the geometric structure. According to experiments, magnetic nanowires have diameters of the order of 10-100 nm [29–33]. The smallest cells to simulate nanowires of 10 and 35 nm contain about 2500 and 31000 atoms, respectively. In the case of arrays of magnetic nanowires, it is important to consider the structure in the edges or surface of the nanowires and the distances between the walls of the nanowires in the array. However, calculations of the MDE of systems with a large number of atoms are very expensive. The present paper is devoted to reduce the computation time of the calculation of the MDE and MDAE with high precision by using the symmetries of the periodic magnetic systems, and in doing so, to reduce the computation time of the simulations of the magnetic properties of large periodic layered magnetic materials. An analysis of the Ewald method in its traditional form [34–36], to calculate the MDE of periodic magnetic systems, whose time complexity is O(N2), has been carried out in the present research, finding that many matrix elements of the Ewald summation method are identical to others, depending on the type of Bravais lattice cell and if the basis atoms of the cell of the periodic magnetic system satisfy certain conditions or symmetries. When these symmetries are applied, the number of matrix elements that should be calculated is approximately or even equal to the number Nof magnetic moments of the periodic cell. Periodic layered magnetic systems such as nanowires, slabs and multilayers satisfy the symmetries. The usual required precision of the MDEs is high, about 10−6eV/cell, and the computation of the matrices to obtain MDEs with that precision is very expensive and hence, the application of these symmetries reduces drastically the computing time of the calculation of the MDEs of large magnetic periodic systems. The novelty of the analysis and application of the symmetries of periodic magnetic systems is that this analysis was not performed in former forms of the Ewald summation method: The method developed by Perram et al. [37], which is an O(N3/2) method and is based on the linkedcell spatial decomposition technique [38, 39], the Particle Mesh Ewald, PME, method [40, 41], based on using fast Fourier transform, FFT, techniques to evaluate the recip- 3 rocal space part of the Ewald method and has a complexity O(NlogN), and the fast multipole method, FMM [42], which is an O(N) method. The method devised by Perram et al. [37] reduces the time complexity of the traditional Ewald summation without approximations. This paper is organized as follows. Section 2 is devoted to the theory of the magnetostatic dipolar interaction energy of a lattice of magnetic moments or dipoles. Section 3 explains the analysis of the symmetries of periodic magnetic systems to reduce the computation time of the MDE. Section 4 is a brief description of the three Ni fcc periodic magnetic systems studied: Nanowires, slabs and spheres. The next section is the discussion of the computation timing results of the calculations of the MDE of Ni fcc nanowires, slabs and spheres up to 32000 magnetic moments in the periodic cell, using and not using the symmetries. The last section is a comparison of the two strategies to reduce the computation time: Parallel calculations not using the symmetries, serial calculations using the symmetries, and the combination of parallelization and use of the symmetries. 2. Theory of the Magnetostatic dipolar energy of a lattice of magnetic moments 2.1. Magnetic dipolar interaction energy between two magnetic dipoles In classical electromagnetism, the magnetic vector potential ~ Aat point #” r, due to a magnetic moment ~ mlocated at the origin of coordinates, is given by ~ A(#” r)=µ0 4π ~ m×#” r r3.(1) The magnetic field ~ Bat point #” rproduced by the magnetic moment ~ m, located at the origin, is calculated from the above magnetic vector potential, Eq. 1 and is given by ~ B(#” r)= #” ∇ × ~ A(#” r)=−µ0 4π~ m∇21 r− #” ∇(~ m· #” ∇)1 r =µ0 4π~ m8π 3δ(#” r)+3#” r(~ m·#” r) r5−~ m r3.(2) The contact term is proportional to the Dirac delta function in three dimensions, δ(#” r). This term cancels out if #” r,0. Therefore, this term is not usually considered in the magnetic field due to a magnetic moment. If the magnetic moment ~ mis located at the point #” r, then the magnetic field produced at the point #” r′due to the magnetic moment ~ mlocated at #” ris given by ~ B(#” r′)=µ0 4π3(#” r′−#” r)(~ m(#” r)·(#” r′−#” r)) |#” r′−#” r|5 −~ m(#” r) |#” r′−#” r|3.(3) The magnetic dipolar interaction energy between the magnetic moment ~ mlocated at #” rand the magnetic moment ~ m′located at #” r′is given by: Em,m′=−~ m′(#” r′)·~ B(#” r′)=µ0 4π~ m′(#” r′)·~ m(#” r) |#” r′−#” r|3 −3(~ m′(#” r′)·(#” r′−#” r))(~ m(#” r)·(#” r′−#” r)) |#” r′−#” r|5,(4) where the expression of the magnetic field at the point #” r′, Eq. 3, has been used. 4 2.2. Magnetostatic dipolar energy of a lattice of magnetic moments The magnetostatic dipolar energy of a lattice of magnetic moments consists on the summation of the magnetic dipolar interaction energies, Eq. 4, between the magnetic moments of the lattice. This summation is given by: Ed=1 2 µ0 4πX iX jX n~ mi·~ mj |~ Rn+~ ri−~ rj|3 −3(~ mi·(~ Rn+~ ri−~ rj))(~ mj·(~ Rn+~ ri−~ rj)) |~ Rn+~ ri−~ rj|5,(5) where iand jdenote the atoms in the cell, ~ riis the position of atom iin the cell, ~ miis the magnetic moment of atom iand the vector ~ Rn+~ ri−~ rjconnects the magnetic moments ~ miand ~ mj, located at ~ Rn+~ riand~ rj, respectively. ~ Rnis a lattice site: ~ Rn=na~ a+nb~ b+nc~ cand nstands for n=(na,nb,nc). The sum runs over all the lattice sites ~ Rnexcept over that for which the denominator in Eq. 5 is zero. If all the magnetic moments ~ miand ~ mjof the cell are parallel to the direction b n, i.e., it is a ferromagnetic system, then ~ mi=mib n, with i=1−N, and the magnetostatic dipolar energy is given by: Ed(b n)=1 2 µ0 4πX iX j mimjMi j(b n),(6) where the quantities Mi j(b n)=Mi j are called the ferromagnetic dipolar Madelung constants and are given by Mi j(b n)=X n1 |~ Rn+~ ri−~ rj|3−3(b n·(~ Rn+~ ri−~ rj))2 |~ Rn+~ ri−~ rj|5.(7) These constants can be further developed, taking into account the angle θ′ ni j between the magnetic moments and the vector ~ Rn+~ ri−~ rj: Mi j =X n 1−3(cosθ′ ni j)2 |~ Rn+~ ri−~ rj|3,(8) where the cosine of the angle θ′ ni j is given by: cosθ′ ni j =b n·(~ Rn+~ ri−~ rj) |~ Rn+~ ri−~ rj| .(9) The magnetic moment, the vector ~ Rn+~ ri−~ rjand the angle θ′ ni j are depicted in Fig. 1). x y z Rn +i-j θnij ϕnij mθ| nij Figure 1: Vector ~ Rn+~ ri−~ rj, the spherical angles θni j and φni j of this vector with respect to the Cartesian reference system, the magnetic moment ~ mand the spherical angle θ′ ni j between the vectors ~ Rn+~ ri−~ rjand ~ m. The magnetic moment ~ m=mb n. Using the complex spherical harmonic for l=2 and m=0 [43–45], given by: Ycomplex 2,0=r5 16π(3cos2θ−1) ,(10) 5 the Madelung constants are written as: Mi j =−r16π 5X n Ycomplex 2,0(θ′ ni j, φ′ ni j) |~ Rn+~ ri−~ rj|3.(11) The complex spherical harmonic Ycomplex 2,0(θ′ ni j, φ′ ni j) can be written as: Ycomplex 2,0(θ′ ni j, φ′ ni j)= 2 X m=−2 D2,m,0(α, β, γ)Ycomplex 2,m(θni j, φni j),(12) where α,βand γare the Euler angles that define the direction of the magnetic moments with respect to a Cartesian reference system, D2,m,0are the Wigner rotation matrix elements [46–48], and θni j and φni j are the spherical angles of the vector ~ Rn+~ ri−~ rjwith respect to the Cartesian reference system (See Fig. 1). Inserting Eq. 12 into Eq. 11, the Madelung constants turn into: Mi j =−r16π 5 2 X m=−2 D2,m,0(α, β, γ) X n Ycomplex 2,m(θni j, φni j) |~ Rn+~ ri−~ rj|3.(13) If the magnetic moments are in units of the Bohr magneton µB, then: Ed(b n)=µ2 B 2 µ0 4πX iX j mimjMi j .(14) The quantity µ2 Bµ0/8πis equal to 1/c2in atomic Rydberg units. Therefore, the MDE in atomic Rydberg units is given by: Ed(b n)=1 c2X iX j mimjMi j .(15) The Madelung constants Mi j can be written as a combination of real spherical harmonics, using the relationship between the real and complex spherical harmonics (See Eq. 35 in the Appendix A) [43, 45]: Mi j =khD2,0,0X n Yreal 2,0(θni j, φni j) |~ Rn+~ ri−~ rj|3+ 2 X m=1 D2,m,0X n Yreal 2,m(θni j, φni j)+iYreal 2,−m(θni j, φni j) √2(−1)m|~ Rn+~ ri−~ rj|3+ 2 X m=1 D2,−m,0X n Yreal 2,m(θni j, φni j)−iYreal 2,−m(θni j, φni j) √2|~ Rn+~ ri−~ rj|3i,(16) where k=−r16π 5. Let’s define the matrix elements Sm(i,j): Sm(i,j)=X n Yreal 2,m(θni j, φni j) |~ Rn+~ ri−~ rj|3,(17) with m=−2,−1,0,1,2. Using Eq. 16 and the quantities Sm(i,j) defined in Eq. 17 and with some algebra calculations, the Madelung constants can be written as: Mi j =khS0(i,j)D2,0,0+S1(i,j) √2(−D2,1,0+D2,−1,0)+ iS −1(i,j) √2(−D2,1,0−D2,−1,0)+S2(i,j) √2(D2,2,0+D2,−2,0)+ iS −2(i,j) √2(D2,2,0−D2,−2,0)i.(18) The Wigner rotation matrix elements have some properties [46–48] that can be used to simplify the Madelung constants Mi j (See Appendix B): D2,−1,0=−D∗ 2,1,0and D2,−2,0=D∗ 2,2,0. Taking into account these properties and 6 Eq. 18, and with some additional algebra, the Madelung constants Mi j can be finally written as: Mi j =khS0(i,j)D2,0,0−S1(i,j)√2 Real(D2,1,0) +S−1(i,j)√2 Imag(D2,1,0)+S2(i,j)√2 Real(D2,2,0) −S−2(i,j)√2 Imag(D2,2,0)i.(19) The MDE is calculated using Eqs. 15, 17 and 19. The matrix elements Sm(i,j) in Eq. 17 are computed by means of the Ewald summation method [34, 35]. 2.3. Magnetostatic dipolar anisotropy energy The magnetostatic dipolar anisotropy energy, MDAE, is the difference between the magnetostatic dipolar energies for two different magnetization directions. For instance, in the case of magnetizations ~ Mparallel and perpendicular to the c-axis b cof a layered system (this axis is perpendicular to the plane of the layers), the magnetostatic dipolar anisotropy energy is given by: MDAE(k,⊥)=Ed(b nkb c)−Ed(b n⊥b c),(20) where b n=~ M/Mis a unitary vector along the magnetization, b cis a unitary vector along the c-axis and the magnetostatic dipolar energies Ed’s are given by Eq. 15, with the corresponding orientations of the magnetizations. In the study of the MDAE of layered magnetic systems, the directions of interest are the axis perpendicular and parallel to the plane of the layers. The parallel axis is not well defined, because there are many axes lying in the plane of the layers. Usually the perpendicular axis is denoted as the zaxis and the parallel axis could be any axis lying in the xy plane. This is the convention that has been followed in this paper, unless otherwise noted. 3. Analysis of the Symmetries of the S matrices The magnetostatic dipolar energy, MDE, is a longrange interaction and hence, in a periodic system of N magnetic moments, the interaction of each magnetic moment iwith every other magnetic moment jmust be calculated. The MDE of periodic magnetic systems is calculated by means of the Ewald’s lattice summation method [34–36]. This method is used to calculate the five matrix elements Sm(i,j) (m=-2,-1,0,1,2), Eq. 17, related to the magnetostatic dipolar interaction between the magnetic moments iand jin all the cells (the real cell and the replicated cells). These five matrix elements are then, used to calculate the matrix element Mi j through Eq. 19. The Madelung constants or matrix elements Mkk, with k=1−N, are all identical and hence, only one of these matrix elements should be calculated. On other hand, Mi j =Mji. Therefore, there are only N(N−1)/2+1 different Madelung constants Mi j in the summation of Eq. 6. This means that the time complexity of the calculation of the MDE, Eq. 6, of periodic magnetic systems using the traditional Ewald method is O(N2), because there are N(N−1)/2+1 different matrix elements Mi j in that equation, or five times N(N−1)/2+1 different matrix elements Sm(i,j), if Eq. 19 is considered. A detailed analysis of 7 the time complexity of the traditional Ewald method was published by Petersen [41] and Wang and Holm [36]. Each Madelung constant Mi j (or equivalently each of the five Sm(i,j) matrix elements) is a summation over the infinite number of lattice sites ~ Rnof the magnetic periodic system (See Eq. 7). The summation to calculate Mi j in Eq. 7 is obtained by applying cutoffdistances in real and reciprocal spaces and it converges rapidly. The MDEs and MDAEs are very small energies. The Madelung constants Mi j must be calculated with enough precision to ensure MDEs, and especially MDAEs, with a precision of at least 10−6eV/cell. The MDAE is the difference between two MDEs and both must be enough accurate, to obtain the MDAE as an accurate difference, without effects due to the compensation of errors. A strategy to reduce the computation time of the calculation of the MDE with high precision, without changing the cutoffdistances, consists on using the symmetries of the periodic magnetic system. To use those symmetries, one should consider and analyze the S matrices in more detail. The S matrices of the Ewald method applied to the calculation of the magnetostatic dipolar energy are given by Eq. 17, where Yreal 2,mis a real spherical harmonic of l=2 and m=−2,−1,0,1,2, ~ riand ~ rjare the positions of the iand jatoms in the cell, respectively, and ~ Rnis a Bravais lattice vector or lattice site, i.e., ~ Rn=na~ a+nb~ b+nc~ c, with ~ a,~ band~ cequal to the lattice vectors of the cell, and na,nb and ncare integer numbers. The position vector of atom i is given by ~ ri=(xi,yi,zi). The real spherical harmonics in the definition of Sm(i,j), Eq. 17, depend on the spherical angles θni j and φni j and are obtained from the Eqs. 37 in Appendix A, by making the following replacements in those equations: x replaced by Xn+xi−xj,yreplaced by Yn+yi−zjand z replaced by Zn+zi−zjand r=|~ Rn+~ ri−~ rj|. For instance, the real spherical harmonic Yreal 2,1is given by: Yreal 2,1(θni j, φni j)=r15 4π (Xn+xi−xj)(Zn+zi−zj) |~ Rn+~ ri−~ rj|2. (21) The S matrices are symmetric, i.e., Sm(i,j)=Sm(j,i). This is taken into account in all the calculations and this does not depend on the type of Bravais lattice cell, nor in the values of the vectors ~ ri−~ rjof the basis atoms of the cell. If the vectors ~ ri−~ rjand ~ rk−~ rlof the basis atoms of the cell and the Bravais lattice cell satisfy certain conditions, then the matrix elements Sm(i,j) are equal to ± Sm(k,l), with m=−2,−1,0,1,2. These symmetries or conditions allow us to reduce the number of matrix elements that should be calculated. The general symmetry or condition that must be satisfied is as follows: If any vector ~ Rn+~ ri−~ rjis equal to the vector ~ T+~ rk−~ rl, such as Yreal 2,m(θni j, φni j) |~ Rn+~ ri−~ rj|3=± Yreal 2,m(θtkl, φtkl) |~ T+~ rk−~ rl|3,(22) and the vector ~ Tis a Bravais lattice vector, i.e., ~ T=~ Rp= pa~ a+pb~ b+pc~ c, with pa,pband pcbeing integer numbers, then Sm(i,j)±Sm(k,l). 8 If ~ T=~ Rp, then Eq. 22 implies a reordering of the sums in the summation that defines Sm(i,j), Eq. 17, but the value of the summation does not change, except for a sign in some cases, depending on the value of m. If ~ Tis not a Bravais lattice vector, then Eq. 22 is not satisfied and the absolute value of the summation in Eq. 17 changes. Eq. 22 will be satisfied depending on the values of~ ri−~ rj and ~ rk−~ rl, and on the type of Bravais lattice cell. There are at least eight symmetries or conditions of ~ ri−~ rjand ~ rk−~ rlthat could lead to the fulfillment of Eq. 22. The first and second conditions satisfy Eq. 22 for any of the 14 Bravais lattice cells: 1) If~ ri−~ rjis equal to~ rk−~ rlthen Sm(k,l)=Sm(i,j) for any value of m: If~ ri−~ rj=~ rk−~ rl, then ~ T+~ rk−~ rl=~ Rn+~ rk−~ rl=~ Rn+~ ri−~ rj, θtkl =θnkl =θni j and φtkl =φnkl =φni j, which implies that Sm(i,j)=Sm(k,l). An obvious and particular case of this symmetry is ~ r1−~ r1=~ r2−~ r2=...=~ rN−~ rN. This means that all the elements of the diagonal of the corresponding Smmatrices are identical: Sm(i,i)=Sm(1,1) for any value of i, with m=−2,−1,0,1,2. Hence, to calculate the elements of the diagonal of Sm, only one element, Sm(1,1), has to be calculated. Notice that S0(1,1) is different from S2(1,1) and so on for m=−2,−1,0,1,2. Only five matrix elements are necessary to calculate the corresponding diagonals of the five Smmatrices. From physical arguments, without math calculations, it can be also derived that Sm(i,i)=Sm(1,1) for any value of i: The interaction of atom iof the cell with all the atoms iof the replicated cells, is the same that the interaction of atom jwith all the atoms jof the replicated cells. The fact that Sm(i,i)=Sm(1,1) for any value of i, is applied in all the calculations, not only on the calculations that use the symmetries of the periodic magnetic system. 2) If ~ ri−~ rjis equal to -(~ rk−~ rl) then Sm(k,l)=Sm(i,j) for any value of m: This symmetry comes from the fact that the S matrices are symmetric: If ~ ri−~ rj=−~ rk−~ rl→~ Rn+~ ri−~ rj=~ Rn+ ~ rl−~ rk, which means that Sm(i,j)=Sm(l,k). The matrices Smare symmetric matrices, therefore Sm(l,k)=Sm(k,l), and Sm(i,j)=Sm(k,l). The following six symmetries or conditions, 3-8, do not fulfill Eq. 22 for all the Bravais lattice cells. 3) If xi−xj=-(xk−xl), yi−yj=yk−yland zi−zj=zk−zl, then, taking into account the dependence on xi−xjand xk−xlof Y2,m: S−2(k,l)=−S−2(i,j) S−1(k,l)=S−1(i,j) S0(k,l)=S0(i,j) (23) S1(k,l)=−S1(i,j) S2(k,l)=S2(i,j). The Eqs. 23 can be proved as follows. If xi−xj=-(xk− 9 metries of these periodic systems. For instance, the calculation of a Ni fcc nanowire of 5025 atoms takes about 27000 seconds not using the symmetries, and about 130 seconds using the symmetries, in the mentioned cluster and with the same cutoffdistances. The reduction factor is about 200 for nanowires, 370 for slabs and about 390 for spheres, for large values of the number Nof atoms (or magnetic moments). Another way to realize the reduction of the computation time is to fix the amount of the computation time of the calculations and to find out the number of atoms of the nanowires calculated in that same fixed amount of time. For instance, a calculation of a nanowire of 19000 atoms using the symmetries and another calculation of a nanowire of 2900 atoms not using the symmetry, will take approximately the same amount of time, about 6000 seconds. The computation time of the calculations not using the symmetries grows faster than the computation time of the calculations using the symmetries. This can be noticed in Fig. 7. Another interesting fact is that the use of the symmetries has a much larger impact on the calculations of large systems than on the calculations of small systems: For the smallest nanowire studied without using the symmetries, the reduction factor of the computation time is about six and for the largest nanowire studied without using the symmetries, which has 5025 atoms, the reduction factor is about 200. A similar behaviour has been found in slabs and spheres. 0 5000 10000 15000 20000 25000 30000 Number of atoms 0 10000 20000 30000 40000 Computation time (seconds) No symmetries used Symmetries used 0 5000 10000 15000 20000 25000 30000 Number of atoms 0 10000 20000 30000 40000 Computation time (seconds) No symmetries used Symmetries used 0 5000 10000 15000 20000 25000 30000 Number of atoms 0 10000 20000 30000 40000 Computation time (seconds) No symmetries used Symmetries used Figure 7: (Color online) Computation time vs number Nof atoms of the calculations of Ni fcc nanowires, slabs and spheres (top, central and bottom panel, respectively), not using and using the symmetries of the S matrix. 16 It can be also noticed in Fig. 7 that the computation time of the calculations using the symmetries has approximately the same dependence on the number Nof atoms in the three types of geometries studied, which implies that the number of atoms is much more relevant than the type of geometry. 5.3. Analysis of the Computation Time of the Calculations done Using the Symmetries of Ni fcc Systems To analyze the dependence on Nof the computation time of the calculations done using the symmetries, the two main contributions to the computation time have been considered: The time to find and analyze the symmetries of the periodic magnetic system and to determine which matrix elements should be calculated, ta, and the time to calculate the matrix elements that should be calculated, tm. These two times are plotted in Fig. 8. The computation time tais larger than tm, and taincreases faster than tmas the number Nof atoms increases. The computation time to find and analyze the symmetries and to calculate the matrix elements of Ni fcc nanowires, slabs and spheres are plotted in Fig. 8, respectively. The time to find and analyze the symmetries, ta (See Fig. 8), has the same dependence on Nas the total computation time of the calculations done not using the symmetries (See Fig. 7). If the symmetries of the periodic magnetic system are not used, then tmis proportional to N2for large values of N. If the symmetries are used, then tmis proportional to 0 5000 10000 15000 20000 25000 30000 Number of atoms 0 10000 20000 30000 40000 Computation time (seconds) Symmetry analysis Calculation matrix elements 0 5000 10000 15000 20000 25000 30000 Number of atoms 0 10000 20000 30000 40000 Computation time (seconds) Symmetry analysis Calculation matrix elements 0 5000 10000 15000 20000 25000 30000 Number of atoms 0 10000 20000 30000 40000 Computation time (seconds) Symmetry analysis Calculation matrix elements Figure 8: (Color online) Computation time to analyze the symmetries, ta, and to calculate the matrix elements, tm, vs number of atoms of the calculations of Ni fcc nanowires, slabs and spheres (top, central and bottom panel, respectively), when the symmetries of the periodic magnetic system are used. 17 N. To understand why tmis proportional to Nif the symmetry is used and to N2if the symmetries are not used, tmhas to be further analyzed. The time to calculate the S matrix elements is proportional to the number Mof matrix elements: tm=aM. If the symmetries are not used, then the number Mof matrix elements of S is not reduced and Mis equal to N(N−1)/2+1, which means that M is proportional to N2for large values of the number Nof atoms of the cell. If the symmetries are used, then Mis approximately equal to N. This fact can be noticed in Fig. 9. The number Mof matrix elements that should be calculated vs the number of atoms of the nanowires, slabs and spheres, when the symmetries are used, is plotted in Fig. 9. It can be noticed in that Figure that the number of matrix elements that should be calculated is practically equal to the number Nof atoms. For instance, the rightmost point in Fig. 9 corresponds to a Ni fcc nanowire with 28917 atoms and 29066 matrix elements. In the case of nanowires, M is very close to N, but not exactly equal to N.Mis exactly equal to the number Nof atoms of the slabs, for any value of N.Mis slightly higher than the number Nof atoms of the spheres, for any value of N. The fact that M=Nfor slabs is probably due to the higher symmetry of the slabs, compared to nanowires and spheres. This dependence of the number Mof matrix elements that should be calculated on the number Nof magnetic moments explains the dependence on Nof the computation time to calculate the matrix elements, tm. That com- 0 5000 10000 15000 20000 25000 30000 Number of atoms 0 10000 20000 30000 No. matrix elements that should be calculated Symmetries used M=N 0 5000 10000 15000 20000 25000 30000 Number of atoms 0 10000 20000 30000 No. matrix elements that should be calculated Symmetries used M=N 0 5000 10000 15000 20000 25000 30000 Number of atoms 0 10000 20000 30000 No. matrix elements that should be calculated Symmetries used M=N Figure 9: (Color online) Number Mof matrix elements that should be calculated vs number Nof atoms of the Ni fcc nanowires, slabs and spheres (top, central and bottom panels, respectively), when the symmetries of the periodic magnetic systems are used. 18 putation time is proportional to M:tm=bM. If the symmetries of the periodic magnetic system are used, then tm=bM ≈bN, and if the symmetries are not used or the periodic magnetic system has not symmetries, then tm=bM =[N(N−1)/2+1] ≈cN2. 6. Parallelization of the calculations The second strategy to decrease, in practical terms, the computation time of the MDE of large magnetic periodic systems is to implement and run parallel calculations. The calculation of the elements of the S matrix has been parallelized in a simple way in the present research: If the parallel calculation is carried out by pprocessors, then every processor calculates M/pmatrix elements, where Mis the number of matrix elements that should be calculated. Parallel calculations of the MDE of the following Ni fcc periodic systems have been carried out, using up to 96 processors and not using the symmetries: A Ni fcc nanowire of radius 19a(4513 magnetic moments), a Ni fcc slab of 2600 atomic layers (5200 magnetic moments) and a Ni fcc sphere of radius 5.2a(5000 magnetic moments), with a=3.52 Å. The computation times of these parallel calculations as a function of the number pof processors are plotted in Fig. 10. A serial calculation, i.e., using only one processor, and not using the symmetries takes about 20000 seconds in the case of the nanowire and slab and 21600 seconds in the case of the sphere. Using the 96 processors, the larger number available in our computer resources, and 012 24 36 48 60 72 84 96 Number of processors 100 1000 10000 Computation time (seconds) No symmetries used Amdahl law s=0.005 012 24 36 48 60 72 84 96 Number of processors 100 1000 10000 Computation time (seconds) No symmetries used Amdahl law s=0.007 012 24 36 48 60 72 84 96 Number of processors 100 1000 10000 Computation time (seconds) No symmetries used Amdahl law s=0.008 Figure 10: Computation time, in logarithmic scale, vs number of processors of the calculations of the following periodic magnetic systems, not using the symmetries: A Ni fcc nanowire of radius 19a, a Ni fcc slab of 2600 atomic layers and a Ni fcc sphere of radius 5.2a, a=3.52 Å (top, central and bottom panels, respectively). 19 not using the symmetries, the computation times of the nanowire, the slab and the sphere are about 820, 900 and 710 seconds, respectively (See Fig. 10 and Table 1). Hence, the reduction factors due to the parallelization are between 22 and 30. 6.1. Comparison of Parallelism and Use of the Symmetries Parallelization of the calculations reduces, obviously, the computation time, but much less than the use of the symmetries. As it has just been indicated above, parallel calculations using 96 processors and not using the symmetries of a nanowire, a slab and a sphere take about 820, 900 and 710 seconds, respectively. Those computation times are longer than the corresponding computation times of serial calculations using the symmetries: 100, 140 and 55 seconds, respectively (See Table 1). The reduction factors due to the use of the symmetries are between 143 and 390, about 6-13 times larger than the reduction factors due to the parallelization, which are between 22 and 30 (See Table 1). Therefore, it is much more efficient (less computation time and also less computer resources) to run serial calculations using the symmetries than to run parallel calculations without using the symmetries. The combination of parallelism and the analysis of the symmetries is also possible. This type of parallel calculations are based on the parallelization of the algorithm to calculate the matrix elements and the algorithm to analyze Table 1: Computation times and reduction factors of the calculations of a nanowire of radius 19a(up), a slab of 2600 atomic layers (center) and a sphere of radius 5.2a(down), a=3.52 Å, as a function of the number of processors and the use of the symmetries. Number of Use of the Time Reduction processors symmetries (seconds) factor 1 no 20000 – 96 no 820 24 1 yes 100 200 4-6 yes 40 500 1 no 20000 – 96 no 900 22 1 yes 140 143 4-6 yes 60 333 1 no 21600 – 96 no 710 30 1 yes 55 390 4-6 yes 30 720 the symmetries. The parallel version of the calculation of the matrix elements distributes evenly these calculations among the processors. The parallel version of the analysis of the symmetries also distributes evenly the calculations, but it is more complex. As it was explained before, the serial algorithm to analyze the symmetries of the magnetic system consists on a conditioned comparison of the pairs of vectors ~ ri−~ rj and ~ rk−~ rlof the basis atoms of the cell. The pairs that 20 satisfy some of the symmetries or conditions are not compared anymore. The parallel version of that algorithm distributes the comparisons of the vectors as follows. Each processor compares a=nv/pvectors, where nv is the total number of vectors~ ri−~ rjthat will be compared and pis the number of processors. Processor kcompares the vectors from 1+ak to a(k+1)−1. The index kruns from 0 to p−1. The last processor, k=p−1, runs from 1+(p−1)a to nv. Finally, the master node gathers the results. The parallel algorithm to analyze the symmetries is an O(N3/p) algorithm. It reduces in an important amount the computation time of the analysis, but with a price: The result of the parallel version of the algorithm to analyze the symmetries is that the total number of matrices Mpthat should calculated using pprocessors is approximately equal to pN, where Nis the number of atoms, if the magnetic system satisfies the conditions and symmetries. In a serial calculation, the number of matrices that should be calculated is approximately N. The number of matrices that should be calculated, Mp, is therefore, larger than in a serial calculation, although of the same order of magnitude. This increase of the number of matrices that should be calculated has not an important impact on the computation time to calculate the matrix elements, because the calculation of the matrix elements is an O(Mp/p)=O(N) task in a parallel calculation, the same as in a serial calculation. Hence, the result of the parallelization of both algorithms is an important reduction of the total computation time, compared with the other types of calculations, as can be noticed in Table 1. Parallel calculations of the nanowire, slab and sphere take about 40, 60 and 30 seconds, respectively, using between four and six processors and the symmetries, and about 100, 140 and 55 seconds, respectively, using one processor (a serial calculation) and the symmetries (See Table 1). The parallel calculations with 4-6 processors are the optimal ones: Calculations with a larger or a smaller number of processors take longer. These calculations with 4-6 processors are about two times faster than the serial calculations using the symmetries. 6.2. Comparison with the Amdahl law The dependence of the computation time of the calculation of the MDE, not using the analysis of the symmetries, on the number pof processors has been compared with the Amdahl law [49]. This law states that the time of a calculation using pprocessors is given by: Tp=Ts(1 +(p−1)s) p,(33) where sis between 0 and 1 and is the proportion of the code that remains serial, because is not parallelized or can not be parallelized, and Tsis the time of a calculation with one processor (serial run). The results of the parallel calculations of the nanowire, the slab and the sphere have been fitted to the Amdahl law, Eq. 33, obtaining a value of sequal to 0.005, 0.007 and 0.008 for the nanowire, slab and sphere, respectively. The fitting functions are plotted 21 as solid lines in Fig. 10. These values of smeans that a 0.5-0.8 % of the code is serial and a 99.5-99.2 % is parallelized. According to the Amdahl law, Tpshould be about 300 seconds for the nanowire, slab and sphere, using 96 processors. However, the computation time using 96 processors is 820, 900 and 710 seconds, respectively (See Table 1). This is an expected behaviour: The predictions of the Amdahl law are not at all correct for large values of p. This can be better noticed in the plots of the speedup, Ts/Tp, of the parallel calculations of the nanowire, slab and sphere in Fig. 11. The real speedup matches very well the speedup predicted by the Amdahl law for p<=20, but it deviates largely from the predictions for p>20. The speedup is approximately constant above p>20. This is the expected behaviour of the speedup when the size of the problem is relatively small. In the present case, the size of the problem is the number Nof magnetic moments, which is about 5000 for the studied nanowire, slab and sphere. Larger values of Nwill improve the real speedup. Finally, it should be considered that the basis atoms or magnetic moments of a periodic cell could be such that their position coordinates do not satisfy the conditions or symmetries explained in section 3 of this paper. In that case, the whole periodic system (lattice cell +basis atoms) would be a low symmetry system. If the periodic magnetic system has a low symmetry, then the use of the symmetries does not reduce the number of matrix ele- 012 24 36 48 60 72 84 96 Number of processors 0 10 20 30 40 50 60 70 Speedup = Ts/Tp No symmetries used Amdahl law s=0.005 012 24 36 48 60 72 84 96 Number of processors 0 10 20 30 40 50 60 70 Speedup = Ts/Tp No symmetries used Amdahl law s=0.007 012 24 36 48 60 72 84 96 Number of processors 0 10 20 30 40 50 60 70 Speedup = Ts/Tp No symmetries used Amdahl law s=0.008 Figure 11: Speedup vs number of processors of the calculations of the following periodic magnetic systems, not using the symmetries: A Ni fcc nanowire of radius 19a, a Ni fcc slab of 2600 atomic layers and a Ni fcc sphere of radius 5.2a, a=3.52 Å (top, central and bottom panels, respectively). 22 ments or the reduction is very small and hence, the reduction of the computation time is very small. Therefore, for periodic magnetic systems with a low symmetry, the parallel calculations without using the symmetries are faster than the serial calculations using or not using the symmetries. 7. Conclusions Two strategies to reduce the computational effort of the calculation of the MDE of large magnetic periodic systems have been studied. The first strategy consists on an analysis of the symmetry properties of periodic magnetic systems of Nmagnetic moments, in order to reduce the number of matrix elements that should be calculated in the traditional Ewald method used to calculate the MDE. The number of matrix elements of this method is N2/2 and hence, its time complexity is O(N2). It has been shown that if the periodic magnetic system has certain symmetries, then there are many matrix elements are identical to other elements, except the sign of some matrix elements. This reduces the number Mof matrix elements that should be calculated to approximately N, instead of N2/2, according to computation timing experiments carried out in large periodic magnetic systems, such as large Ni fcc nanowires up to 32000 magnetic moments. This decreases considerably the computation time of the MDE. This reduction is in contrast with the fact that the analysis of the symmetries is an O(N3) task, which increases the time complexity of the traditional Ewald method. The origin of this contrast is that the MDE and MDAE are very small energies and therefore, the usual required precision to calculate these energies is so high, 10−6eV/cell, that the calculation of the matrix elements is very expensive and, in practice, the computations carried out using the analysis of the symmetries are much faster, in spite of the larger time complexity of the analysis of the symmetries. The second strategy to reduce the computation time of the calculations of the MDE is the parallelization of the calculations, without using the symmetries of the system. For periodic magnetic systems with high symmetry, the parallelization of the calculations of the MDE reduces the computation time, but much less than the use of the symmetries in a serial calculation and using more computational resources. However, for periodic magnetic systems with low symmetry, the use of the symmetries reduces very little the computation time of a serial calculation and running parallel calculations without using the symmetries is faster. Finally, the use of both, the parallelization and the symmetries of the periodic magnetic system, is the fastest procedure. There are several future lines of improvement of the present research. The most important one consists on finding and studying more symmetries or conditions of the periodic magnetic system that reduce the number of matrix elements of the Ewald summation method that should be calculated. The reduction of the time complexity of the analysis of the symmetries of periodic magnetic systems and the application of the proposed analysis of 23 the symmetries to the non-traditional forms of the Ewald summation method are also important research lines. Finally, the derivation of the mathematics involved in the calculation of the MDE and MDAE of periodic cells of non-collinear magnetic dipoles is underway. Acknowledgments This work was supported by MINECO of Spain (Grant MAT2014-54378-R), Junta de Castilla y Le´ on (Grants VA050U14 and VA124G18) and the University of Valladolid. The facilities provided by Centro de Proceso de Datos - Parque Cient´ ıfico of the University of Valladolid are acknowledged. Appendix A. Complex and Real Spherical Harmonics The complex spherical harmonics are defined by [43– 45]: Ycomplex l,|m|= Θl,|m|(θ)ei|m|φ Ycomplex l,−|m|=(−1)|m|Θl,|m|(θ)e−i|m|φ=(−1)|m|Ycomplex∗ l,|m|(34) The complex spherical harmonics can be also written as a combination of real spherical harmonics [43–45]: Ycomplex l,0=Yreal l,0 Ycomplex l,|m|=(−1)|m| √2Yreal l,|m|+iYreal l,−|m|(35) Ycomplex l,−|m|=1 √2Yreal l,|m|−iYreal l,−|m| The real spherical harmonics as a combination of complex spherical harmonics are obtained from Eqs. 35: Yreal l,0=Ycomplex l,0 Yreal l,|m|=1 √2Ycomplex l,−|m|+(−1)|m|Ycomplex l,|m|(36) Yreal l,−|m|=i √2Ycomplex l,−|m|−(−1)|m|Ycomplex l,|m| The real spherical harmonics of l=2 are given by [43– 45]: Yreal 2,0=r5 16π(3cos2θ−1) =r5 16π 3z2−x2−y2−z2 r2 Yreal 2,1=r15 4πsinθcosθcosϕ=r15 4π xz r2 Yreal 2,−1=r15 4πsinθcosθsinϕ=r15 4π yz r2(37) Yreal 2,2=r15 16πsin2θcos2ϕ=r15 16π x2−y2 r2 Yreal 2,−2=r15 16πsin2θsin2ϕ=r15 4π xy r2. Appendix B. Some Properties of the Rotation Matrix Elements The Wigner rotation matrix elements are given by [46– 48]: Dl,m′,m(α, β, γ)∗=(−1)m′−mDl,−m′,−m(α, β, γ) Dl,m′,m(α, β, γ)=e−im′αdl,m′,m(β)e−imγ(38) dl,m′,m(β)=(−1)m′−mdl,−m,−m′(β) If the above definition is applied to the specific cases l= 2, m′=±1, ±2 and m=0, the following matrix elements 24 are obtained: D2,1,0=e−iαd2,1,0(β) D2,−1,0=−eiαd2,1,0(β)=−D∗ 2,1,0(39) D2,−2,0=D∗ 2,2,0. References [1] G. Y. Guo, H. Ebert, W. M. Temmerman, First principles determination of the magnetisation direction of Fe monolayers in noble metals, J. Phys.: Condens. Matter 3 (1991) 8205–12. [2] L. Szunyogh, B. ´ Ujfalussy, P. Weinberger, Magnetic anisotropy of iron multilayers on Au(001): Firstprinciples calculations in terms of the fully relativistic spin-polarized screened KKR method, Phys. Rev. B 51 (1995) 9552–9. [3] M. T. Johnson, P. J. H. Bloemen, F. J. A. den Broeder, J. J. de Vries, Magnetic anisotropy in metallic multilayers, Rep. Prog. Phys. 59 (1996) 1409–58. [4] I. Cabria, H. Ebert, A. Y. Perlov, Microscopic origin of the magnetocrystalline anisotropy energy of ferromagnetic-semiconductor multilayers, Europhys. Lett. 51 (2000) 209–15. [5] I. Cabria, A. Y. Perlov, H. Ebert, Magnetization profile and magnetocrystalline anisotropy of ferromagnet-semiconductor heterostructure systems, Phys. Rev. B 63 (2001) 104424. [6] R. G´ omez-Abal, A. M. Llois, Magnetic anisotropy of extended defects and vicinal surfaces of 3d transition metals, Phys. Rev. B 65 (2002) 155426. [7] R. Lorenz, J. Hafner, Magnetic structure and anisotropy of thin Fe films on Cu(001) substrates, Phys. Rev. B 54 (1996) 15937–49. [8] L. Szunyogh, B. ´ Ujfalussy, P. Weinberger, Magnetic structure and anisotropy in Fe/Cu(001) over- and interlayers with antiferromagnetic interlayer coupling, Phys. Rev. B 55 (1997) 14392–6. [9] M. Rapini, R. A. Dias, B. V. Costa, Phase transition in ultrathin magnetic films with long-range interactions: Monte Carlo simulation of the anisotropic Heisenberg model, Phys. Rev. B 75 (2007) 014425. [10] H. Leon, E. Estevez-Rams, Magnetic dipolar interaction in ultrathin films: Structural and size effects, J. Magn. Magn. Materials 321 (2009) 2150–9. [11] C. Liu, S. D. Bader, Perpendicular surface magnetic anisotropy in ultrathin epitaxial Fe films, J. Vac. Sci. Technol. A 8 (1990) 2727–31. [12] R. Allenspach, M. Stampanoni, A. Bischof, Magnetic domains in thin epitaxial Co/Au(111) films, Phys. Rev. Lett. 65 (1990) 3344–7. [13] P. Krams, F. Lauks, R. L. Stamps, B. Hillebrands, G. G¨ untherodt, Magnetic anisotropies of ultrathin 25