PaperlayerView exact revision on arXiv

Source-linked AI summary

Massively parallel symplectic algorithm for coupled magnetic spin dynamics and molecular dynamics

J. Tranchida, S. J. Plimpton, P. Thibaudeau, A. P. Thompson

arXiv:1801.10233v1cond-mat.stat-mechcond-mat.mes-hallcond-mat.mtrl-sci

TL;DR

Accurate simulations of magneto-elasticity matter for testing theories and interpreting experiments, while large-scale spin–lattice dynamics requires coupled spin and atomic motion. The paper implements SD–MD in LAMMPS with a symplectic integrator and parallel sectoring algorithm. The parallel implementation was statistically equivalent to the serial algorithm and performed well above 500 atoms per processor, at about 4–5 times the cost of analogous MD.

  • Problem

    Accurately simulating magneto-elasticity is important for testing competing theories and predicting and interpreting experiments.

    Source: §1 · p. 1

  • Method

    The paper implements coupled spin–lattice dynamics in LAMMPS using a symplectic Suzuki–Trotter integration algorithm and synchronous sub-lattice sectoring for parallelization.

    Source: §1 · p. 2

  • Results

    The parallel algorithm was statistically equivalent to the serial algorithm and performed well for cases with more than 500 atoms per processor, at 4 to 5 times the cost of an analogous NVE MD LAMMPS run.

    Source: §6 · p. 11

  • Takeaways & Limitations

    The sectoring method works for ordered and disordered particle configurations when each processor domain is at least twice the cutoff distance for short-range spin–spin interactions.

    Source: §6 · p. 11

    The study assumes uncorrelated white noise; when dynamics approach the simulated relaxation timescales, the Markov hypothesis may break down and colored noise may be more appropriate.

    Source: §4.3 · p. 7

Abstract

from arXiv · show

A parallel implementation of coupled spin-lattice dynamics in the LAMMPS molecular dynamics package is presented. The equations of motion for both spin only and coupled spin-lattice dynamics are first reviewed, including a detailed account of how magneto-mechanical potentials can be used to perform a proper coupling between spin and lattice degrees of freedom. A symplectic numerical integration algorithm is then presented which combines the Suzuki-Trotter decomposition for non-commuting variables and conserves the geometric properties of the equations of motion. The numerical accuracy of the serial implementation was assessed by verifying that it conserves the total energy and the norm of the total magnetization up to second order in the timestep size. Finally, a very general parallel algorithm is proposed that allows large spin-lattice systems to be efficiently simulated on large numbers of processors without degrading its mathematical accuracy. Its correctness as well as scaling efficiency were tested for realistic coupled spin-lattice systems, confirming that the new parallel algorithm is both accurate and efficient.

1. Introduction

Magneto-elasticity links magnetic moments and mechanical strain, but its interacting scales make accurate, large-system simulations challenging. The paper presents a general LAMMPS implementation of coupled spin–lattice dynamics with a symplectic integrator and parallel algorithm.

1. Introduction

  • Magneto-elasticity couples stress-driven magnetic rotation with field-driven rotation that generates strain, and simulations can test theories and interpret experiments.These effects also underpin actuator and sensor capabilities in smart devices [1].Source: §1 · p. 1
  • Local quantum methods offer high accuracy for hundreds of atoms over picoseconds, while semiclassical molecular or spin dynamics can simulate millions over nanoseconds.The approaches therefore fill complementary roles in studying interactions across length scales.Source: §1 · p. 2
  • SD–MD uses effective parametrized potentials to separate spin and lattice variables from electrons and derive coupled equations of motion.The approach is motivated by spin and lattice timescales being slow relative to electronic timescales.Source: §1 · p. 2
  • The paper implements SD–MD in LAMMPS, combining a Suzuki–Trotter-based symplectic algorithm with synchronous sub-lattice sectoring for parallel simulations.It presents algorithm accuracy and analyzes parallel accuracy and scalability.Source: §1 · p. 2

2. Spin–lattice dynamics in the microcanonical ensemble

The microcanonical treatment describes spin and lattice evolution from a shared Hamiltonian, with exchange and anisotropic interactions shaping their coupling. The resulting equations propagate atomic positions, momenta, and spins as a coupled system.

2. Spin–lattice dynamics in the microcanonical ensemble

  • In the microcanonical ensemble, the system is treated as closed with fixed total energy, composition, volume, and shape.The section introduces equations of motion for both magnetic and mechanical degrees of freedom in this setting.Source: §2 · p. 2

2.1. Equations for atomistic spin dynamics

  • Atomistic spin dynamics represents each atomic site with a classical spin vector whose adiabatic motion follows Thomas precession, with finite-temperature damping modeled by Landau–Lifshitz–Gilbert dynamics.The spin vector is ℏg_i times the dimensionless unit vector s_i; damping is connected to the fluctuation-dissipation theorem.Source: §2.1 · p. 3
  • Distance-dependent exchange naturally couples spins and lattice, but radial exchange alone cannot model anisotropic magnetostriction from magneto-crystalline anisotropy.Accurate, realistic spin–lattice simulations therefore need spin–orbit interactions; evaluating their very small energies from first principles requires high accuracy and is difficult.Source: §2.1 · p. 3

2.2. Equations for spin–lattice dynamics

  • The coupled spin–lattice Hamiltonian combines spin energy, particle kinetic energy, and a mechanical binding potential, yielding equations for atomic positions, momenta, and spins.Momentum propagation includes both mechanical forces and magnetic forces dependent on interatomic distance and spin configuration; the system has 9 × N coupled equations.Source: §2.2 · pp. 3–4

3. Integration algorithm

The integration algorithm splits the coupled spin–lattice evolution into non-commuting operators using a Suzuki–Trotter decomposition. Its single-spin propagator preserves spin norm, and tests report good conservation of energy and total magnetization up to a timestep of 10−2 ps.

Operator splitting

  • A Suzuki–Trotter decomposition advances the coupled variables by splitting the Liouville operator into separate position, momentum, and spin operators.The operators need not commute, and the displayed decomposition is accurate to O(...).Source: §3 · p. 4

Spin propagation

  • Spin rotations are non-commuting, so each spin’s precession vector must be computed immediately before updating it with the latest neighboring-spin orientations.This ordering makes straightforward concurrent updates unsuitable; the proposed synchronous sublattice or sectoring algorithm supports nonuniform or disordered systems while preserving the Suzuki–Trotter properties.Source: §3 · p. 5
  • The single-spin propagator uses low-order trigonometric approximants that preserve spin norm without norm rescaling.The approximants satisfy cos²(∥ω_i∥Δt) + sin²(∥ω_i∥Δt) = 1; avoiding rescaling also avoids breaking symplecticity and time-reversal properties.Source: §3 · p. 5

Numerical accuracy

  • The algorithm remains reasonably accurate at timesteps up to 10−2 ps, while scaling begins to break down above that threshold.The simulations measured energy and total-magnetization deviations; the stated breakdown indicates the numerical-stability limit is being approached.Source: §3 · pp. 5–6
  • In a closed 500-atom cobalt simulation, energy and total magnetization were very well conserved, and the algorithm remained reasonably accurate up to 10−2 ps.The figure plots energy per atom and |M| for several timestep sizes, with inset deviations against timestep size.Source: §3 · p. 5, Figure 2 · p. 6

4. Spin–lattice dynamics in a canonical ensemble (NVT)

The canonical-ensemble approach couples spin and lattice degrees of freedom to thermal reservoirs using Langevin thermostats and measures their temperatures. Simulations show temperature relaxation through spin–lattice coupling, while the stochastic treatment assumes uncorrelated white noise.

4. Spin–lattice dynamics in a canonical ensemble (NVT)

  • The NVT approach connects the spin–lattice system to an external energy reservoir, using Langevin random forces and damping in accord with the fluctuation–dissipation theorem.The authors note that damping contracts phase-space volume and makes the dynamics non-Hamiltonian.Source: §4 · p. 6

4.1. Connecting spins to a random bath

  • Spin-bath fluctuations enter the stochastic Landau–Lifshitz–Gilbert equation as a multiplicative random torque, with noise amplitude related to thermostat temperature through a fluctuation–dissipation relation.The midpoint Stratonovich prescription is used because it preserves individual-spin norm and time micro-reversibility.Source: §4.1 · p. 6

4.2. Measuring lattice and spin temperatures

  • For spin temperature, the authors choose Nurdin et al.’s definition because it has the same stated domain of validity as the kinetic lattice-temperature expression and has been used for spin thermalization.The lattice kinetic-temperature expression is assumed sufficient for this work despite relying on approximations including equivalence of ensembles and ergodicity.Source: §4.2 · p. 7

4.3. Thermalizing the spin–lattice system

  • The stochastic treatment uses uncorrelated white noise under the Markov hypothesis, which may break down when dynamical timescales approach those of the simulated relaxation processes.The authors identify colored noise, such as an Ornstein–Uhlenbeck process, as a better approximation in that regime.Source: §4.3 · p. 7
  • Thermostatted degrees of freedom stayed at the target temperature, while the others relaxed to it in a few hundred picoseconds.This occurred whether the thermostat was applied to spins or to lattice motion, consistent with the spin–lattice coupling time shown in Fig. 1.Source: §4.3 · p. 7, Figure 3 · p. 8
  • In a 1,372-atom NVE simulation, spin and lattice temperatures equilibrated while total energy remained constant as energy repartitioned among kinetic, atomic-potential, and spin-potential components.The spins began in a high-temperature paramagnetic state, the lattice at 200 K, and the atomic potential energy at a minimum.Source: §4.3 · pp. 7–8

5. Parallel implementation of the spin–lattice algorithm

The sectoring algorithm parallelizes spin updates while preserving the requirement that each spin uses current neighbor values. Tests report serial-like energy conservation, statistical equivalence of spin trajectories, and efficient scaling with a roughly 4–5× cost relative to standard MD.

5.1. Synchronous sub-lattice algoritm

  • Sectoring requires only finite-cutoff spin interactions, extending beyond regular lattices to vibrating, disordered, and diffusive systems.Sectors are chosen so their spatial extent exceeds the interaction cutoff; processors can then update the same sector concurrently without interprocessor spin communication.Source: §5.1 · pp. 8–9
  • Before each sector update, processors exchange current ghost-spin values needed for neighboring interactions, then update sectors in a fixed order.The implementation updates all adjacent ghost spins for convenience, including some that have not changed; atoms are assigned to sectors once at the start of each timestep.Source: §5.1 · p. 9

5.2. Accuracy of the parallel algorithm

  • Parallel and serial trajectories have indistinguishable energy fluctuations; their magnetization differs slightly by timestep but shows the same step-to-step variation.The authors attribute numerical trajectory differences to the different spin-update order and describe the trajectories as statistically equivalent.Source: §5.2 · p. 10, Figure 7 · p. 10

5.3. Scaling results

  • In strong scaling, communication costs grow as atoms per node shrink, and the smaller problem is limited to 16 nodes by the sector-width cutoff requirement.The communication effect is less pronounced for the larger problem.Source: §5.3 · pp. 10–11, Figure 8 · p. 10
  • Weak-scaling simulation rate decreases very little with node count, while SD–MD costs about 4–5× standard MD.The reported cost comparison is with standard NVE MD using the same EAM mechanical potential.Source: §5.3 · pp. 10–11, Figure 8 · p. 10, Figure 9 · p. 11

6. Conclusion

The LAMMPS SD–MD implementation combines symplectic spin propagation with coupled or uncoupled spin–lattice simulations. Its serial accuracy and parallel equivalence and scaling were tested, with the coupled algorithm running 4 to 5 times slower than lattice-only NVE MD.

6. Conclusion

  • The LAMMPS package supports coupled or uncoupled simulations and uses statistically equivalent symplectic Suzuki–Trotter spin propagation in serial and parallel, with different spin-update orders.The sectoring method works for ordered and disordered configurations when each processor domain is at least twice the short-range spin-interaction cutoff width.Source: §6 · p. 11
  • Serial coupled dynamics conserved total internal energy and magnetization norm up to second order in timestep size in NVE tests, while NVT relaxation tests thermostatted either spins or lattice.The reported serial accuracy analysis covered both microcanonical and canonical ensembles.Source: §6 · p. 11
  • 4 to 5 times slower than lattice-only NVE LAMMPS MD, the coupled SD–MD algorithm showed this reported runtime comparison.The comparison is with an analogous run treating only atomic-lattice motion.Source: §6 · p. 11
  • Parallel SD–MD matched the serial algorithm’s conservation of total energy and magnetization norm as processor counts varied, with good strong- and weak-scaling performance above 500 atoms per processor.The parallel integration uses sectoring, also called the synchronous sublattice method.Source: §6 · p. 11
  • Implementing the methods in open-source LAMMPS makes them available for coupled spin–lattice simulations and for testing new models and algorithms.Applications to magnetoelastic effects in magnetic alloys were deferred to subsequent publications.Source: §6 · p. 11

Appendix A. Definition of the spin Hamiltonian

The appendix defines the spin Hamiltonian as the magnetic system’s total energy and presents it as a sum of interaction terms. In the equations of motion, spin-derived mechanical forces are computed from its partial derivative.

Appendix A. Definition of the spin Hamiltonian

  • The spin Hamiltonian represents the magnetic system’s total energy, and its partial derivative supplies the mechanical force from spins in the equations of motion.The passage introduces the Hamiltonian’s general form as a sum of terms.Source: Appendix A · p. 13
  • The Hamiltonian’s listed terms include Zeeman and exchange interactions, along with two magnetic-anisotropy terms.Zeeman and exchange interactions are defined in Section 2.Source: Appendix A · p. 13
  • The remaining listed terms are Dzyaloshinskii–Moriya, magneto-electric, and dipolar interactions.The appendix subsections provide definitions of these interactions.Source: Appendix A · p. 13

Appendix A.1. Parametrization of the exchange interaction

The appendix models exchange with an effective parametrization that uses a short-range cutoff and a continuous isotropic function. It illustrates the approach with exchange data for Fe, Ni, Co, and BFO.

Appendix A.1. Parametrization of the exchange interaction

  • The appendix situates exchange in a quantum-mechanical account of magnetism, attributing effective coupling between magnetic moments to electronic properties including the Pauli principle, Coulomb repulsion, and electron hopping.The passage notes that this mechanism has no classical analogue.Source: Appendix A.1 · p. 13
  • The effective exchange parametrization assumes intensity decreases rapidly with few sign oscillations, allowing a rigid cutoff radius Rc for neighboring spins.The plotted Rc accounts for the three first neighbor shells.Source: Appendix A.1 · p. 13, Figure A.10 · p. 14
  • For 3d atoms, exchange intensity is approximated by a continuous isotropic function based on the Bethe–Slater curve and parameterized by α, δ, and γ.The parameters are α in eV, δ in Å, and dimensionless γ.Source: Appendix A.1 · p. 13
  • The examples cover bcc-Fe, fcc-Ni, fcc-Co, and antiferromagnetic BFO, using ab initio exchange data for the metals and inelastic neutron-scattering measurements for BFO.The appendix identifies fitted coefficients for these four materials in Table A.1.Source: Appendix A.1 · pp. 13–14, Table A.1 · p. 14

Appendix A.2. Spin orbit coupling

The appendix describes spin-orbit-related interactions used to represent anisotropy and noncollinear spin states. It notes that the simple single-site anisotropy model does not capture the more complex forms found in most magnetic crystals.

Appendix A.2. Spin orbit coupling

  • Spin-orbit coupling is important for magnetocrystalline anisotropy and Dzyaloshinskii–Moriya interactions, despite energies orders of magnitude below exchange.The passage links technological applications to effective interactions arising from spin-orbit coupling.Source: Appendix A.2 · p. 14
  • A simple one-site magnetocrystalline-anisotropy model shapes the energy’s angular dependence with spherical harmonics, with a uniaxial form based on the squared spin projection onto an axis.The anisotropy constant’s sign determines easy-axis versus easy-plane behavior.Source: Appendix A.2 · p. 14
  • The simple anisotropy model does not represent the more complex forms common in magnetic crystals, such as cubic anisotropy, and lacks clear dependence on lattice parameters.A more sophisticated two-site, non-local Néel anisotropy model is limited to the pseudo-dipolar term.Source: Appendix A.2 · p. 14
  • A fast-decreasing two-site term can be omitted by redefining the exchange function J(r_ij), and the implementation uses the same three-parameter function as for exchange.This is the g1(r) term described in the appendix.Source: Appendix A.2 · p. 14
  • Exchange combined with spin-orbit coupling can produce noncollinear spin states; exchange coupled to Dzyaloshinskii–Moriya interaction is a common simulation approach, and that interaction is key to stabilizing skyrmions.The Dzyaloshinskii–Moriya vector specifies the interaction’s intensity and direction.Source: Appendix A.2 · p. 14

Appendix A.3. Magneto–electric interaction

For materials such as BFO, the paper represents atomic-scale interplay between magnetism and permanent dielectric polarization through an antisymmetric spin–spin effective interaction. The polarization sets its direction and intensity and may also be induced by a sufficiently strong external electric field.

Appendix A.3. Magneto–electric interaction

  • In BFO, magnetism and permanent dielectric polarization coexist and interact at the atomic scale.Source: Appendix A.3 · p. 14
  • The interaction is modeled by an antisymmetric spin–spin effective interaction, a particular case of the Dzyaloshinskii–Moriya vector.Source: Appendix A.3 · p. 15
  • The screened dielectric atomic polarization E specifies the interaction’s direction and intensity and can be induced by a sufficiently intense external electric field.Source: Appendix A.3 · p. 15

Appendix A.4. Dipolar interaction

The dipolar interaction is a long-range, non-local contribution associated with domain-wall stabilization and magnetic-domain nucleation, but its computational cost limits large magnetic simulations. The paper discusses two ways to compute it, neither yet implemented in the LAMMPS SD-MD package, and notes conditions under which it can be omitted.

Appendix A.4. Dipolar interaction

  • The long-range dipolar interaction stabilizes domain walls and contributes to magnetic-domain nucleation, so it is also called a demagnetizing field.Source: Appendix A.4 · p. 15
  • Because it has no finite cutoff, the dipolar interaction scales computationally as N^2 and is a main cost limit for large magnetic simulations.Source: Appendix A.4 · p. 15
  • Two approaches considered for computing the dipolar effective field and force are Suzuki–Trotter-based r-RESPA and periodic-boundary Ewald sums.Source: Appendix A.4 · p. 15
  • Neither approach is implemented in the LAMMPS SD-MD package yet; dipolar effects can be omitted for sufficiently small systems below the paramagnetic limit that avoid domain-wall nucleation.Source: Appendix A.4 · p. 15

Appendix B. Parametrization of a mechanical potential

The appendix describes how magnetic interactions are handled alongside mechanical potentials, including the standard use of EAM for metals and the subtraction of magnetic ground-state energy from fitted mechanical potentials. It contrasts this approach with a more detailed magnetic potential and identifies improved empirical magneto-mechanical potentials as future work.

Appendix B. Parametrization of a mechanical potential

  • The standard mechanical interaction for metals in MD uses an empirical many-body potential such as the embedded atom method (EAM).Source: Appendix B · p. 15
  • Dudarev et al.’s semi-empirical potential for magnetic α-iron derives its embedding form from Stoner and Ginzburg–Landau models, but does not naturally capture spin-orbit anisotropy.Source: Appendix B · p. 15
  • Unlike that detailed potential, this paper fixes each atom’s magnetic-moment magnitude and evolves only its direction, while noting that implementing the detailed potential is conceptually possible.Source: Appendix B · p. 15
  • Because EAM potentials are fitted to experimental or ab initio data, magnetic interactions are already implicitly included and are subtracted from the mechanical potential.Source: Appendix B · p. 15
  • For a ferromagnet, the appendix identifies exchange energy as the dominant magnetic contribution to the ground-state energy.Source: Appendix B · p. 15
  • Combining experimental and ab initio results to derive better empirical magneto-mechanical potentials is left for future work.Source: Appendix B · p. 15
Loading 1801.10233v1…