`ensemble` (GLSD)

GLSD evolves both direction and magnitude of spin:R:3 magnetic moments. It implements the generalized Langevin equation of Ma and Dudarev, Phys. Rev. B 86, 054416 (2012), Eqs. (18)–(19) and (33), with a stochastic Heun discretization. It requires a Spin potential and uses its existing magnetic force, atomic force and virial interfaces. For fixed-magnitude dynamics use SIB.

Quick start

The following examples use dimensionless damping 0.01 as a starting setting for pilot runs, not a default or a calibrated material parameter. No equilibrium-moment scan or element table is required to set damping. Read the guidance below before interpreting physical relaxation times. After loading a Spin potential and a model.xyz with magnetic moments:

# Deterministic dissipative spin dynamics; fixed atomic positions
ensemble nve_glsd damping 0.01 stemp -1
time_step 0.1
dump_xyz 100 spin.xyz precision double spin mforce
dump_restart 1000
run 1000

Other ensemble choices:

# Fixed lattice; spins coupled to a 300 K bath
ensemble nvt_glsd 300 300 100 lattice off damping 0.01 seed 731

# Coupled spin and lattice dynamics (lattice NHC)
ensemble nvt_glsd 300 300 100 damping 0.01 seed 731

# Coupled spin, lattice and cell dynamics (MTTK)
ensemble npt_glsd temp 300 300 iso 0 0 tperiod 100 pperiod 1000 damping 0.01 seed 731

The example damping is not a fitted Fe parameter, and the timestep is not a stability guarantee. Check timestep convergence and the potential’s coverage of the sampled moment magnitudes. A potential trained at one fixed magnitude need not describe longitudinal energetics.

Syntax and parameters

ensemble nve_glsd damping <value> [gamma <value>] [stemp <-1|value>] [seed <value>]
ensemble nvt_glsd <T_1> <T_2> <T_coup> [lattice on|off] \
         damping <value> [gamma <value>] [stemp <-1|0|value>] [seed <value>]
ensemble npt_glsd temp <T_1> <T_2> <pressure_terms> \
         [tperiod <tau_temp>] [pperiod <tau_press>] \
         damping <value> [gamma <value>] [stemp <-1|0|value>] [seed <value>]

damping: Required, finite, non-negative and dimensionless. It is \(\eta=\hbar\gamma'_s\), where \(\gamma'_s\) is the paper’s Cartesian GLSD dissipation coefficient. This is a new input normalization of that coefficient, not another physical term and not the paper’s \(\gamma_s\) or Gilbert alpha. There is no material-independent physical default. Zero disables damping and noise. The code converts it once per ensemble specification to the constant internal mobility \(\lambda=\eta g_L^2/\hbar=\eta g^2\hbar\), in muB^2/(eV ps). It is never inferred from element identity, initial or instantaneous magnetic moments, volume or force curvature.
gamma: Positive, finite numerical precession coefficient \(g=g_L/\hbar\); default \(2/\hbar \simeq 3038.54\) in the eV/ps convention. GLSD interprets an explicit value as \(g_L=g\hbar\), so changing gamma at fixed damping also changes \(\lambda\) quadratically and the noise amplitude accordingly. Keyword order has no effect. Keep gamma at its default unless a different precession coefficient is intended. This is not the paper’s damping parameter.
stemp: -1 disables noise while retaining deterministic damping; 0 follows the linearly ramped ensemble target; a positive value fixes the spin bath temperature in K. Default: -1 for nve_glsd, 0 otherwise. nve_glsd rejects stemp 0 because it has no ensemble temperature.
seed: Positive integer; default 12345. Keys the counter-based Gaussian stream.
lattice: Only for nvt_glsd, defaults to on. off fixes positions and velocities. nve_glsd always fixes them; npt_glsd always evolves atoms and the cell. The name nve_glsd denotes the fixed-lattice driver: with positive mobility the spin subsystem dissipates energy, so its energy is not conserved even when noise is disabled.

The NVT temperature and lattice coupling parameters have the same meaning as NHC. T_coup is required even with lattice off, when it has no effect on spins. The NPT prefix has the syntax of MTTK; all GLSD options follow the MTTK options. GLSD uses dimensionless damping; SIB uses Gilbert alpha. These keywords select different dissipative equations and are not interchangeable.

Choosing damping

For equilibrium pilot runs, damping 0.01 is a convenient numerical starting choice, not a demonstrated universal optimum or a material fit. Compare, for example, 0.005, 0.01 and 0.02 after equilibration; check correlation times and shrink the timestep before comparing statistical averages. This suggested bracket is not a calibrated range or a confidence interval. More damping is not guaranteed to improve mixing or stability. For quantitative time-dependent predictions, use independently calibrated kinetics instead of this pilot-run choice.

The transformation removes units from the input; it does not remove physical material dependence. A common damping across elements defines a common model bath coupling, not identical relaxation times for all elements.

A Spin potential specifies energy and forces; these do not by themselves specify coupling to the dissipative bath. The fluctuation–dissipation relation fixes the noise amplitude given damping, gamma and temperature; it does not determine damping. Merely obtaining a stable run or the correct equilibrium temperature cannot identify a physical relaxation rate.

For equilibrium sampling

At a common, constant bath temperature, positive constant mobilities have the same continuum Cartesian Gibbs invariant distribution, provided it is normalizable and the dynamics equilibrates. Mobility changes equilibration and correlation times. Choose it by pilot runs: check stationarity and mixing of energy, magnetization and moment-length statistics; repeat at smaller timestep and several positive mobilities. Compare equilibrium observables within sampling uncertainty. Do not infer physical relaxation times from an uncalibrated sampling coefficient. With unequal spin/lattice bath temperatures or driven fields, even steady observables can depend on mobility; the equilibrium argument does not apply.

For physical time-dependent dynamics

Calibrate against a relaxation measurement or an independently justified dynamical calculation for the same material, temperature, state and bath. A coefficient from another Hamiltonian or a particular Fe experiment is only a candidate, not a transferable material constant. Report its source and test sensitivity. Do not increase it just to accelerate the simulation.

For an initial local estimate

In a fixed-environment radial mode near a stable nonzero moment, \(E(r)\simeq E(r_0)+\kappa(r-r_0)^2/2\), with positive curvature \(\kappa\) in eV/muB^2, the linear deterministic radial equation gives

\[\frac{d\, \delta r}{dt}=-\lambda\kappa\, \delta r,\qquad \tau_\parallel=\frac{1}{\lambda\kappa},\qquad \lambda=\frac{1}{\kappa\tau_\parallel}.\]

This is a local estimate derived from the implemented equation, not a general fit to a coupled finite-temperature relaxation spectrum. Curvature alone is insufficient: a target relaxation time is also needed. For example, hypothetical \(\kappa=0.4\) eV/muB^2 and \(\tau_\parallel=0.5\) ps give mobility 5. These numbers are an illustration, not Fe data. A sampling-only chosen time is a numerical choice, not a measured material timescale.

Knowing a Gilbert damping alpha alone does not uniquely calibrate this full variable-magnitude model. One scalar mobility constrains both transverse and longitudinal relaxation; if the measured rates cannot be matched together, this isotropic GLSD model is insufficient. Relabeling alpha or adding a default would not resolve that physical limitation.

For each candidate, check timestep convergence separately from parameter sensitivity. In the simple radial harmonic model the explicit-Heun stability condition is \(\lambda\kappa\Delta t<2\) (time in ps), but accuracy needs smaller steps and precession or other coupled modes may impose tighter limits. Without kinetic data, users can study equilibrated static statistics under the stated assumptions, or report a mobility sensitivity study; the potential alone cannot justify quantitative relaxation times.

Paper-to-code conventions

The reference is Ma and Dudarev (2012), specifically Eqs. (18)–(19) for drift, Eqs. (20)–(22) for transverse projection, and Eq. (33) for noise. No LAMMPS implementation is used as the specification.

Write the paper’s dimensionless angular-momentum vector as \(\mathbf s\) and its Landé factor as \(g_L\). The paper defines the magnetic moment \(\mathbf M=-g_L\mu_B\mathbf s\). GPUMD stores the numerical magnetic moment \(\mathbf S=\mathbf M/\mu_B=-g_L\mathbf s\). Hence \(\mathbf H_s=-g_L\mathbf H\), where \(\mathbf H=-\partial E/\partial\mathbf S\). The paper’s equation is

\[d\mathbf s = [\hbar^{-1}\mathbf s\times\mathbf H_s +\gamma'_s\mathbf H_s]dt +\sqrt{2\gamma'_s k_BT}\,d\mathbf W_s.\]

Changing variables gives the equation below with numerical coefficients \(g=g_L/\hbar\) (input gamma) and \(\lambda=g_L^2\gamma'_s\) (internal mobility), when time is in ps. The user input is \(\eta=\hbar\gamma'_s\) (damping). For dimensional moments the latter is \(\lambda_M=(g_L\mu_B)^2\gamma'_s\). The Wiener process changes sign, which does not change its distribution. Thus the negative precession cross product in code is the paper’s positive cross product expressed in magnetic-moment coordinates, not a sign change in the physical model. The default uses \(g_L=2\). The numerical Fe example needs care. Section III.C quotes \(\gamma'_s=5.88\times10^{13}\) eV^-1 s^-1, which literally maps to mobility 235.2 at \(g_L=2\). However, Eq. (28) states \(\gamma'_s=s_0^2\gamma_s/\hbar\). Combining the paper’s \(M_0=2.2\mu_B\), \(s_0=M_0/(g_L\mu_B)=1.1\), and \(\gamma_s=0.008\) gives \(\gamma'_s=1.47065\times10^{13}\) eV^-1 s^-1 and mobility 58.826. These printed definitions and numbers leave a factor-of-four normalization discrepancy. Do not treat the quoted 235.2 as an unambiguous experimental calibration. The validation case using it checks coordinate and time-unit conversion for that specified coefficient only. In the dimensionless convention the two routes give approximately 0.03870 and 0.00968, respectively; neither equals 0.008.

The damping 0.008 was obtained by fitting laser-driven Fe-film magnetization in the coupled spin–lattice–electron model of Ma, Dudarev and Woo (2012), Sec. III.B. That fit also involved lattice–electron coupling and absorbed pulse energy. Transferring it to a new potential is an assumption, not a new experimental fit. If this prior is adopted explicitly, matching transverse damping at a chosen reference moment \(m_0=|S_0|\) gives \(\lambda=m_0^2\gamma_s/\hbar\) and \(\eta=(m_0/g_L)^2\gamma_s\). Keep the resulting mobility constant; recomputing it from the instantaneous moment would change the model and require a different noise treatment.

Heun is the numerical discretization chosen here, not a claim to reproduce the paper’s particular numerical implementation. The paper’s illustrative Heisenberg–Landau Hamiltonian, fitted Fe coefficients and temperature estimator involving energy Hessians are not added to GPUMD: the supplied Spin potential provides the Hamiltonian. Standard thermo temperature remains the lattice kinetic temperature, not the paper’s Eq. (34) spin temperature.

Converting an earlier input

The mobility corresponding to a dimensionless damping is \(\lambda=\eta g^2\hbar\). GPUMD uses \(\hbar=0.06465412\times10.18051/1000\) eV ps. Startup output reports both damping and converted mobility.

Equation and integration

Let \(\mathbf S\) denote the numerical magnetic moment in units of \(\mu_B\), \(\mathbf H=-\partial E/\partial\mathbf S\) the public magnetic force, \(g\) the existing gamma numerical coefficient and \(\lambda\) the constant internal mobility derived from damping. Time is converted from fs to ps:

\[d\mathbf S_i = \mathbf a_i\,dt + \sqrt{2\lambda k_BT}\,d\mathbf W_i, \qquad \mathbf a_i=-g\,\mathbf S_i\times\mathbf H_i+\lambda\mathbf H_i.\]

The precession sign agrees with SIB. There is no Gilbert \(1/(1+\alpha^2)\) factor in this mode. With dimensional magnetic moments, the precession coefficient has units \(\mu_B/(\mathrm{eV\,ps})\); GPUMD’s historical gamma input absorbs the numerical Bohr-magneton unit and is printed in \((\mathrm{eV\,ps})^{-1}\). add_mfield still takes tesla and performs the Zeeman conversion before integration.

For \(r=|\mathbf S|>0\), the deterministic longitudinal and transverse parts are

\[\dot r=\lambda\,\mathbf e\cdot\mathbf H,\qquad \dot{\mathbf S}_\perp=-g\,\mathbf S\times\mathbf H -\frac{\lambda}{r^2}\mathbf S\times(\mathbf S\times\mathbf H), \quad \mathbf e=\mathbf S/r.\]

Both transverse cross-product relations are retained. The transverse direction damping coefficient is \(\lambda/r\), rather than the fixed-magnitude Gilbert value \(g\alpha/(1+\alpha^2)\). One isotropic Cartesian mobility controls longitudinal and transverse relaxation. Independent rates and temperature-dependent LLB coefficients are not implemented. The code evaluates \(\lambda\mathbf H\) directly, without divisions by \(r\). At zero moment a nonzero magnetic force or thermal noise can create a moment; zero force at zero temperature leaves it zero.

For fixed positions and a time-independent conservative energy, the continuum zero-temperature equation satisfies

\[\frac{dE}{dt}=-\lambda\sum_i |\mathbf H_i|^2\leq0.\]

Precession does no work. This does not imply strict energy decrease at every finite numerical step, or for driven fields, a thermostat or moving atoms.

Each step uses one Gaussian vector \(\boldsymbol\xi_i\) with covariance \(2\lambda k_BT\Delta t\,I\):

\[\mathbf S^*=\mathbf S_n+\Delta t\,\mathbf a_n+\boldsymbol\xi,\qquad \mathbf S_{n+1}=\mathbf S_n+ \frac{\Delta t}{2}(\mathbf a_n+\mathbf a^*)+\boldsymbol\xi.\]

The first force is cached at the starting endpoint. After the lattice first half-step and spin predictor, the potential and external field are recomputed at the predictor. The lattice second half-step and spin corrector follow, then a full endpoint evaluation refreshes energy, force, virial and magnetic force. Stateful mechanical modifiers are sampled once, as in SIB. There are two potential evaluations per step after initialization. Heun is second order for smooth deterministic drift, and strong order one for additive noise. It is explicit and conditionally stable. A radial harmonic well requires \(\lambda k\Delta t<2\); accuracy generally needs a much smaller step. Undamped precession has finite-step norm drift with Heun, unlike Cayley SIB. Use timestep convergence and SIB when exact fixed-length preservation is required.

Noise, equilibrium measure and limits

Independent Cartesian Wiener increments obey \(\langle dW_{ia}dW_{jb}\rangle=\delta_{ij}\delta_{ab}dt\). The diffusion coefficient is \(D=\lambda k_BT\), with \(k_B=8.617343\,10^{-5}\) eV/K. Noise is additive: Ito and Stratonovich interpretations coincide in Cartesian coordinates. The same increment is reused in predictor and corrector, never redrawn or counted twice.

For a smooth, confining, time-independent energy at constant temperature, the continuum invariant measure is

\[P(\{\mathbf S\})\prod_i d^3S_i \propto \exp[-E/(k_BT)]\prod_i d^3S_i.\]

The dissipative probability current vanishes for this density; precession is divergence-free and tangent to constant energy. In radial variables the volume element contains \(r^2dr\,d\Omega\), and the Ito radial drift contains \(2\lambda k_BT/r\). No extra radial entropy force is added in Cartesian coordinates. A flat \(dr\,d\Omega\) measure would be a different physical model. This is a classical bath, not an automatic quantum or material-specific longitudinal entropy correction. Finite-step sampling has bias and is not an exact Boltzmann sampler. Thermalization requires positive mobility and temperature; stemp -1 has deterministic damping only.

A bounded neural-network energy or pure Zeeman energy need not define a normalizable distribution over unbounded moments. Short trajectories with such a model establish only local stability, not global thermal equilibrium or extrapolation validity. Different spin measures or material-specific kinetics require separate specification and validation.

Output, restart and multiple GPUs

Use dump_xyz ... precision double spin mforce. GLSD has no independent spin_vel state; noisy paths have no instantaneous differentiable velocity. Entering GLSD clears inherited spin velocity and marks it uninitialized; requesting spin_velocity output then fails explicitly. Switching later to TSPIN lets that integrator initialize its own velocity. No spin kinetic energy is included in the lattice temperature.

dump_restart saves the cell and per-atom state, including magnetic moments. Copy restart.xyz to model.xyz to start a new simulation from that state. Thermostats, barostats and the noise counter are initialized for each new run; seed selects its random sequence. Reusing a seed restarts the same sequence. Neither the original ramp progress nor the random stream is restored, so a stochastic restart need not follow an uninterrupted trajectory. See Spin restart for the input procedure.

Integration and Gaussian generation run on GPU 0, keyed by seed, original atom index, step and component. Multi-GPU NEP distributes predictor/endpoint spins, computes owner forces with complete halos and gathers owner outputs into original atom order. It neither generates per-device noise nor sums ghost magnetic forces. Identical input order and seed therefore give identical random increments on one and four GPUs. Neighbor ordering and floating-point reductions can still differ; compare deterministic short-time errors and stochastic distributions, not long noisy trajectories bitwise. Only elements listed in spin_type evolve and receive Cartesian thermal noise, including active atoms whose current moment is zero. Other elements retain zero effective moments.

ensemble (GLSD)