text:
ArticlePDF Available

# Aging and non-equilibrium critical phenomena in Monte Carlo simulations of 3D pure and diluted Ising models

Authors:

## Abstract

We investigate the non-equilibrium critical evolution of statistical systems and describe of its some features, such as aging and violation of the fluctuation-dissipation theorem. We consider some theoretical results of computations for universal quantities that have been obtained in recent years, such as the exponents determining the scaling behavior of dynamic response and correlation functions and the fluctuation-dissipation ratio, associated with the non-equilibrium critical dynamics. We derive our original Monte Carlo simulation results for 3D pure and diluted Ising models with Glauber and Metropolis dynamics. For these models, we analyze the influence of critical fluctuations, different non-equilibrium initial states, and the presence of nonmagnetic impurities on the two-time dependence of correlation and response functions on characteristic time variables, such as waiting time $t_{\rm w}$ and time of observation $t - t_{\rm w}$ with $t > t_{\rm w}$ . We discuss the obtained values of the non-equilibrium exponents for the autocorrelation and response functions and values of the universal long-time limit of the fluctuation-dissipation ratio $X^{\infty }$ . Our simulation results demonstrate that the insertion of disorder leads to new values of $X^{\infty }$ with $X_{\rm diluted}^{\infty } > X_{\rm pure}^{\infty }$ .
Prog. Theor. Exp. Phys. 2015, 053A01 (20 pages)
DOI: 10.1093/ptep/ptv055
Aging and non-equilibrium critical phenomena
in Monte Carlo simulations of 3D pure and
diluted Ising models
Pavel V. Prudnikov, Vladimir V. Prudnikov, Evgeny A. Pospelov, Petr N. Malyarenko, and
Andrey N. Vakilov
Department of Theoretical Physics, Omsk State University, Omsk 644077, Russia
E-mail: prudnikp@univer.omsk.su
Received September 18, 2014; Revised January 26, 2015; Accepted March 19, 2015; Published May 1, 2015
...............................................................................
We investigate the non-equilibrium critical evolution of statistical systems and describe of its
some features, such as aging and violation of the ﬂuctuation-dissipation theorem. We consider
some theoretical results of computations for universal quantities that have been obtained in recent
years, such as the exponents determining the scaling behavior of dynamic response and correla-
tion functions and the ﬂuctuation-dissipation ratio, associated with the non-equilibrium critical
dynamics. We derive our original Monte Carlo simulation results for 3D pure and diluted Ising
models with Glauber and Metropolis dynamics. For these models, we analyze the inﬂuence of
critical ﬂuctuations, different non-equilibrium initial states, and the presence of nonmagnetic
impurities on the two-time dependence of correlation and response functions on characteristic
time variables, such as waiting time twand time of observation ttwwith t>tw. We discuss
the obtained values of the non-equilibrium exponents for the autocorrelation and response func-
tions and values of the universal long-time limit of the ﬂuctuation-dissipation ratio X.Our
simulation results demonstrate that the insertion of disorder leads to new values of Xwith
X
diluted >X
pure.
...............................................................................
Subject Index A22, A40, A56
1. Introduction
The collective behavior of statistical systems close to critical points is characterized by an extremely
slow dynamics that, in the thermodynamic limit, does not give them time to relax to an equilibrium
state after a change in some thermodynamic parameters and conditions. The non-equilibrium evolu-
tion displays in this case some of the features typically observed in glassy materials, such as aging
and violation of the ﬂuctuation-dissipation theorem (FDT) [15]. It can be detected through deter-
mination of dynamic susceptibilities and correlation functions of the order parameter, the scaling
behavior of which is characterized by universal exponents, scaling functions, and amplitude ratios.
This universality allows one to calculate these quantities in different statistical models [6,7]and
Monte Carlo methods are a natural approach for this analysis.
We consider here some of the theoretical results of computations that have been obtained in recent
years for universal quantities, such as the ﬂuctuation-dissipation ratio, associated with the non-
equilibrium critical dynamics, with particular focus on our original Monte Carlo simulation results
for the 3D pure and diluted Ising models with Glauber and Metropolis dynamics.
© The Author(s) 2015. Published by Oxford University Press on behalf of the Physical Society of Japan.
which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.
PTEP 2015, 053A01 P. V. Prudnikov et al.
The investigation of the inﬂuence of structural disorder on the critical behavior of different statisti-
cal systems remains one of the main problems in condensed-matter physics. The disorder breaks the
translational symmetry of the crystal and thus greatly complicates the theoretical description of the
material. It is known that the critical behavior of a system is characterized by an anomalously large
response even for any weak perturbation; therefore, the presence of even low defect concentration
might lead to drastic changes of this behavior and its characteristics. The description of such systems
requires the development of special analytical and numerical methods [8,9].
Research into the effects produced by weak quenched disorder on critical phenomena has been
carried out for many years [1017]. This work determined the criterion [10] under which the disorder
affects the critical behavior only if the speciﬁc heat exponent of the pure system αis positive. In
this case, a new universal critical behavior, with new critical exponents, is observed. Otherwise, for
α<0, the presence of disorder is irrelevant for the critical behavior. Only systems with an effective
Hamiltonian that is isomorphic near the critical point of the Ising model satisfy this Harris criterion,
and a drastic change in the universality class of the critical behavior is realized with its introduction
to a system of whatever concentration of short-range correlated quenched defects [9,18,19].
Studies of the critical behavior of diluted Ising-like magnets have been carried out by
renormalization-group (RG) methods, by numerical Monte Carlo methods, and experimentally, and
the results of these investigations have been presented in a large number of publications (for a review,
see Ref. [20]). In Refs. [21,22], the authors declared interesting ideas about replica symmetry break-
ing in the systems with quenched disorder. A ﬁeld-theoretic description of the critical behavior of
weakly disordered systems with introduced replica symmetry-breaking potentials, which was car-
ried out in the two-loop approximation [23,24], has shown the stability of the critical behavior of
disordered 3D systems with respect to the replica symmetry-breaking effects. The results of the
investigations conﬁrm the existence of a new universal class of critical behavior for diluted Ising-
like systems. However, it remains unclear whether the values of critical exponents are independent
of the rate of dilution of the system, how the crossover effects change these values, and whether two
or more regimes of the critical behavior exist for weakly and strongly disordered systems. These
questions are of interest to researchers and have been discussed before now [25,26].
This study is devoted to numerical investigation of the non-equilibrium critical dynamics with a
non-conserved order parameter (model A) [27] in 3D pure and site-diluted Ising systems with spin
concentrations p=0.8and 0.6. In the following section, we derive the main concepts and charac-
teristics of non-equilibrium relaxation in statistical systems with slow dynamics and peculiarities in
their non-equilibrium behavior at critical points. In Sect. 2, we also review some theoretical results
of computations in the range of non-equilibrium phenomena that have been obtained in recent years
for some statistical models. In Sects. 3 and 4, we introduce the 3D pure Ising model and 3D Ising
model with quenched point-like defects, and derive the results of our original Monte Carlo inves-
tigations into two-time autocorrelation and response functions and the ﬂuctuation-dissipation ratio
for systems starting from high-temperature disordered and low-temperature ordered initial states.
The critical exponents and the asymptotic values of the universal ﬂuctuation-dissipation ratio X
obtained under these two conditions are compared. The ﬁnal section contains an analysis of the main
results and our conclusions.
2. Non-equilibrium critical dynamics and its main peculiarities and characteristics
Statistical systems with slow dynamics have recently attracted considerable theoretical and exper-
imental interest, in view of the rich scenario of phenomena they display: extreme slowing down
2/20
PTEP 2015, 053A01 P. V. Prudnikov et al.
of relaxation processes, memory effects, hysteresis, etc. After a perturbation, a system with slow
dynamics does not achieve equilibrium, even after a long time, and its dynamics is not invariant
either under time translations or under time reversal, as it should be in thermal equilibrium. Dur-
ing this never-ending relaxation, aging occurs: two-time quantities such as response and correlation
functions depend on characteristic time variables such as waiting time twand time of observation
ttwwith t>twnot via ttwonly. It is important that decays for these two-time quantities as
functions of time of observation ttware slower for larger waiting times tw.
At variance with one-time quantities like the order parameter converging to asymptotic values in
the long-time limit, two-time quantities clearly demonstrate aging. Aging is known to occur in such
disordered and complex systems as glassy materials [13] and only in the last ten years has attention
been focused on simpler systems such as critical ones, whose universal features can be rather easily
investigated by using different methods, and which might provide insight into more general cases
[2830].
Consider a system with a critical point at temperature Tc, order parameter S(x,t), and prepare
it in some initial conﬁguration that might correspond to an equilibrium state at a given tempera-
ture T0. At time t=0, bring the system into contact with a thermal bath with a temperature Tsnot
equal to T0. The relaxation process is expected to be characterized by some equilibration time teq(Ts)
such that, for tteq(Ts), equilibrium is attained and the dynamics is stationary and invariant under
time reversal, whereas, for 0<tteq(Ts), the evolution depends on the speciﬁc initial condition.
Upon approaching the critical point Ts=Tc, the equilibration time diverges as teq τzν,where
τ=(TTc)/Tcis the reduced critical temperature, zis the dynamic critical exponent, νis the
exponent for correlation length, and therefore equilibrium is never achieved. To monitor the time
evolution, we consider the average order parameter M(t)=1
VddxS(x,t), the time-dependent
correlation function of the order parameter
C(t,tw)=1
VddxS(x,t)S(x,tw)S(x,t)S(x,tw),(1)
where < ... > stands for the mean over the stochastic dynamics, and the linear response (suscepti-
bility) Rx(t,tw)to a small external ﬁeld, applied at time tw, which is deﬁned by the relation
R(t,tw)=1
VddxδS(x,t)
δh(x,tw)h=0
.(2)
Note that causality implies that R(t,tw>t)=0. According to the general picture of the relax-
ation process, one expects that, for t>twteq(Ts),C(t,tw)=Ceq(ttw)and R(t,tw)=Req
(ttw),whereCeq and Req are the corresponding equilibrium quantities, related by the ﬂuctuation-
dissipation theorem (FDT):
Req(t)=−1
Ts
dCeq(t)
dt .(3)
The FDT suggests the deﬁnition of the so-called ﬂuctuation-dissipation ratio (FDR):
X(t,tw)=TsR(t,tw)
twC(t,tw)(4)
with t>tw.Fort>twteq (Ts), the FDT yields X(t,tw)=1. This is not generically true in the
aging regime. The asymptotic value of the FDR,
X=lim
tw→∞ lim
t→∞ X(t,tw), (5)
is a very useful quantity in the description of systems with slow dynamics, since X=1when-
ever the aging evolution is interrupted and the system crosses over to equilibrium dynamics, i.e.,
3/20
PTEP 2015, 053A01 P. V. Prudnikov et al.
Tabl e 1. Values of the ﬂuctuation-dissipation ratio Xfor some systems for a high-tem-
perature initial state with m01.
Model T<TcT=TcT>Tc
Free Gaussian ﬁeld [33] Exact – 1/21
dD spherical model [34] Exact 012/d1
1D Ising model [35] Exact – 1/21
2D Ising model [34]MC 0.26(1)
[36]MC 0.340(5)
[37,38]MC 0.330(5)
2D 3-state Potts model [37,38]MC 0.406(1)
2D 4-state Potts model [37,38]MC 0.459(8)
3D Ising model [34]MC 00.40a
3D XY model [39]MC 00.43(4)
aThe value X0.40 was announced in Ref. [34] as result of preliminary simulations
on the 3D Ising model without demonstration of any obtained data either in Ref. [34]or
later publications.
teqTs<. Conversely, X= 1is a signal of an asymptotic non-equilibrium dynamics. More-
over, Xcan be used to deﬁne an effective non-equilibrium temperature Teff =T/X, which might
have some features of the temperature of an equilibrium system, e.g., controlling the direction of heat
ﬂows and acting as a criterion for thermalization [31].
What is known in general about X? As a consequence of the ﬂuctuation-dissipation theorem,
X(T>Tc)=1. On the other hand, on the basis of general scaling arguments for the phase-
ordering regime [32], it has been shown that X(T<Tc)=0. These results are expected to be
actually independent of the speciﬁc system and of its microscopic details. For T=Tc,thereareno
general arguments constraining the value of Xand therefore it has to be determined for each spe-
ciﬁc model. In Table 1we report some of the values that have been found either by exact solutions
or by means of Monte Carlo (MC) simulations (a more complete table can be found in Ref. [28]).
Clearly, X(T=Tc)depends on the model and on the space dimensionality d. Nevertheless, it
has been argued on the basis of scaling arguments [32,34]thatX(T=Tc)should be a universal
quantity associated with the critical dynamics.
At the present time, it is well known that the two-time dependence for autocorrelation and response
functions for systems starting from a high-temperature initial state with m0=0(or m01) satisﬁes
the following scaling forms:
C(t,tw)=AC(ttw)a+1d/z(t/tw)θ1fC(tw/t),
where fC(tw/t)and fR(tw/t)are ﬁnite for tw0,a=(2ηz)/z,θ=θ(2zη)/z,
and θis the initial slip exponent [40]. ARand ACare non-universal amplitudes that are ﬁxed by the
condition fR,C(0)=1. With this normalization, fR,Care universal. From these scaling forms, the
universality of Xfollows as an amplitude ratio X=AR/[(1θ)AC][28,32,34].
A remarkable property of the non-equilibrium relaxation process in this case is the increase
of magnetization m(t)tθfrom a non-zero initial magnetization m01at short times t<tcr
m1/θ+β/zν
0(see, e.g., Fig. 1). The initial rise of magnetization is changed to the well known decay
m(t)tβ/zνfor ttcr. The critical exponents θand θdepend on the dynamic universality class
[27] and have been calculated by the RG method for a number of dynamic models, such as a model
4/20
PTEP 2015, 053A01 P. V. Prudnikov et al.
m(t)
t (MCS/s)
t
tvz
Fig. 1. Time evolution of the magnetization m(t)from the initial state with magnetization m0=0.03 1at
critical temperature Tc=3.499 48 as a result of Monte Carlo simulation of samples with spin concentration
p=0.8and linear size L=128.
with a non-conserved order parameter [26,4042] (model A), a model with an order parameter cou-
pled to a conserved density [43] (model C), and models with reversible mode coupling [44] (models
E, F, G, and J).
Three stages of the non-equilibrium relaxation process can be distinguished. The ﬁrst quasi-
equilibrium stage is observed for small time separation ttwtwwith tw1, where aging does
not exist and the dynamic evolution of the correlation and response functions exhibits a station-
ary part and does not depend on waiting time with C=C(ttw)and R=R(ttw). In the limit
tw/(ttw)→∞, the scaling functions fC(tw/t)and fR(tw/t)in (7) are characterized by power
law dependences fC(tw/t)(tw/t)θ1and fR(tw/t)(tw/t)θthat lead to their quasi-equilibrium
behavior. The second aging regime is realized for times ttwtw1, where the correlation and
response functions are derived by relations
C(t,tw)t2β/(νz)
wˆ
FC(t/tw),
R(t,tw)t2β/(νz)1
wˆ
FR(t/tw),
(7)
and therefore at different waiting times do not superpose and are characterized by dif-
ferent slopes for each tw(in (7) the relation 2β/(ν z)=d/za1was used and the
scaling functions ˆ
FC(t/tw)=[(ttw)/tw]a+1d/z(t/tw)θ1fC(tw/t)and ˆ
FR(t/tw)=[(t
tw)/tw]ad/z(t/tw)θfR(tw/t)were introduced). At long time separations with ttwtw1,
the scaling functions in (7) decay as power laws:
ˆ
FC(t/tw)(t/tw)ca,ˆ
FR(t/tw)(t/tw)cr,(8)
where the exponent ca=d/zaθ=d/zθis the same; it describes the time dependence of
the autocorrelation function in the short-time regime of non-equilibrium behavior [26,45,46]. At this
short-time dynamics stage, aging effects are not developed either. Scaling analysis of the response
function R(t,tw)behavior in the short-time dynamics regime predicts that cr=ca.
Renormalization-group investigations of non-equilibrium critical behavior in d-dimensional pure
systems with n-component order parameter and a weakly dilute random Ising model (RIM) for purely
dissipative dynamics (model A) were carried out in Refs. [47]and[48]. The asymptotic values of
the FDR Xwere calculated with the use of the ε-expansion method (ε=4d) in the two-loop
approximation for pure systems [47],
(X)1
2=1+n+2
4(n+8)ε+ε2n+2
(n+8)2n+2
8+3(3n+14)
4(n+8)+c+Oε3,(9)
5/20
PTEP 2015, 053A01 P. V. Prudnikov et al.
where c=−0.0415 ··· (its analytic expression is given in Ref. [47]), and in the one-loop approxi-
mation for the diluted Ising model [48],
X=1
21
46ε
53 +O(ε). (10)
The following results were obtained: X
3DIs =0.429(6)for the 3D Ising (ε=1,n=1) model,
X
3DXY =0.416(8)for the XY (ε=1,n=2) model, and X
2DIs =0.30(5)for the 2D Ising (ε=2,
n=1) model, which are in good agreement with the Monte Carlo results for these models, given
in Table 1. For the 3D diluted Ising model (which is the only physically relevant case), the value
X
3DRIM 0.416 was obtained. To this order, it is not clear whether disorder really changes X
in a sensible way or not. In any case, this could not be safely stated from low-order computations
since the ε-expansion is known to be not well behaved at d=3[20,24,49]. The results of Monte
Carlo investigations into the non-equilibrium critical dynamics for the 3D diluted Ising model will
be presented below.
If the initial state of a system is characterized by magnetization m0= 0, the renormalization-group
analysis of non-equilibrium dynamics for model A predicts that the correlation C(t,tw)and response
R(t,tw)functions display the following scaling behaviors after a quench to T=Tc[40,50]:
C(t,tw)=aC(ttw)a+1d/z(t/tw)θ1FC(tw/t,t/tm),
where modiﬁcation of these scaling relations in comparison with (6) is connected with a new
timescale tm, set by the initial value of the magnetization m0, and which displays a universal
dependence on it:
tm=Bmm1
0,(12)
where the universal scaling exponent κ>0is given, in terms of static and dynamic equilibrium and
non-equilibrium exponents, by κ=θ+a+β/(νz).
The two-time quantities C(t,tw)and R(t,tw)are homogeneous functions of the three timescales
ttw,tw,andtm. In particular, when tw<ttm, which is always the case with m0=0, the scaling
forms of Cand Rbecome as they are in Eq. (6) with FC,R(x,0)=fC,R(x). In the opposite case, with
large times compared to tm, i.e., tmtw<t, the scaling forms of C(t,tw)and R(t,tw)become [50]
C(t,tw)aC(ttw)a+1d/z(t/tw)¯
θ1¯
FC(tw/t),
θ¯
FR(tw/t), (13)
where the new exponent ¯
θ=−βδ/(νz)=−
1+a+β/(ν z)and ¯
FC,Rare universal scaling func-
tions related to the large-ybehavior of FC,R(x,y). In the aging regime that is realized for times
ttwtwtm, the correlation and response functions are derived by the relations
C(t,tw)t2β/(νz)
w˜
FC(t/tw),
R(t,tw)t2β/(νz)1
w˜
FR(t/tw)
(14)
with the scaling functions ˜
FC,R(t/tw), which decay at a long time separation limit with ttw
twtmas a power law:
˜
FC,R(t/tw)(t/tw)φ,(15)
where the exponent φ=d/za+βδ/(ν z).
6/20
PTEP 2015, 053A01 P. V. Prudnikov et al.
In Ref. [50], the non-equilibrium behavior of the d-dimensional Ising model with purely dissipative
dynamics was investigated during its critical relaxation from a magnetized initial state. The uni-
versal scaling forms of the two-time response and correlation functions were derived within the
ﬁeld-theoretical approach and the associated scaling functions were computed in ﬁrst order of the
ε-expansion. It was shown that aging behavior is clearly displayed and the asymptotic universal
ﬂuctuation-dissipation ratio is characterized by relation
X=4
573
600 π2
100ε+Oε2,(16)
which gives X
3DIs 0.78 for the 3D Ising (ε=1,n=1) model and X
2DIs 0.75 for the 2D Ising
(ε=2,n=1) model. The obtained results were conﬁrmed by Monte Carlo simulations of the 2D
Ising model with Glauber dynamics, from which it was found that X
MC =0.73(1).
3. Monte Carlo simulations of 3D pure Ising models with relaxation from
high-temperature and low-temperature initial states
One of the simplest non-trivial models in which aging occurs is the Ising model in d-dimensions
evolving with a purely dissipative dynamics after a quench to the critical point. Its Hamiltonian on a
hypercubic lattice is given by
H=−J
<i,j>
SiSj,(17)
where J>0is the short-range exchange interaction between spins Siﬁxed at the lattice sites, and
assuming values of Si1. We performed our Monte Carlo simulations using the heat-bath updat-
ing rule [51], simulating a large Ising spin system on a cubic lattice with linear size L=128 with
periodic boundary conditions at Tc=4.5114(1). We compute the magnetization,
M(t)=1
L3
L3
i=1
Si(t),(18)
and the two-time autocorrelation function,
C(t,tw)=1
L3
L3
i=1
Si(t)Si(tw)M(t)M(tw), (19)
where the angle brackets stand for an average over the initial conﬁgurations and MC realizations.
The averaging of C(t,tw)is carried out on 3000 MC runs for every tw.
In the case when we simulated the dynamics of systems starting from the high-temperature initial
state with magnetization value m0=0.02, the response function and ﬂuctuation-dissipation ratio
were calculated using the relations [37,38]:
R(t,tw)=1
TL
3
L3
i=1Si(t)Si(tw+1)SW
i(tw+1),(20)
where SW
i=tanh(Jm=iSm/T),and
X(t,tw)=TR(t,tw)
twC(t,tw)=N
i=1Si(t)Si(tw+1)SW
i(tw+1)
N
i=1Si(t)[Si(tw+1)Si(tw)].(21)
The averaging of R(t,tw)and X(t,tw)is carried out on 90 000 MC runs for every tw.
7/20
PTEP 2015, 053A01 P. V. Prudnikov et al.
t
t
t
t
t
t
t
t
tCt
tCt
ttt
(a) (b)
Fig. 2. Time dependence of correlation function C(t,tw)as a function ttwfor different initial non-equilib-
rium states (a)m01;(b)m0=1.
We must note that the non-equilibrium behavior of R(t,tw)is characterized by very large ﬂuctua-
tions at the critical point and these ﬂuctuations drastically increase in the limit m00. Therefore,
for systems starting from the high-temperature initial state, we used m01, but not m0=0,and
realized the calculations of R(t,tw)with greater statistics than calculations of the autocorrelation
function with the same m0.
In the case when we considered the dynamics of systems starting from the low-temperature initial
state with magnetization value m0=1, we calculated the integrated response function [50,52]:
χ(t,tw)=tw
0
dtR(t,t)=1
TcN
N
i=1Si(t)Si(tw)(22)
with response function determined by (2) and the function Si(tw), which is computed during the
simulation from time t=0to twand is deﬁned by the relation
Si(tw)=
tw
s=0Si(s)SW
i(s).(23)
In the large time limit,
Tcχ(C)=C
0
X(q)dq,
and the ﬂuctuation-dissipation ratio can be deﬁned as
X(t,tw)=lim
C0
Tcχ(t,tw)
C(t,tw).(24)
The MC-obtained time dependences of the autocorrelation function C(t,tw)and the response func-
tion R(t,tw)from observation time ttwfor different initial non-equilibrium states and twvalues
are presented in Figs. 2and 3. The curves of C(t,tw)and R(t,tw)demonstrate the aging effects, i.e.,
the slowing down of time correlations and decreasing of response with increasing system age tw.
In the aging regime, the time dependence of correlation and response functions is character-
ized by scaling relations (7)and(8). These scaling relations are rather well displayed by our data,
as shown in Fig. 4, where we used the exponent values β=0.325(1),ν=0.630(1)[53], and
z=2.024(6)[54]. For time intervals with (ttw)/tw1, we determined the values of exponents
ca=1.333(40)and cr=1.357(18), which demonstrate a good agreement with each other and with
value ca=1.362(19), obtained by the short-time dynamics method in Ref. [45].
8/20
PTEP 2015, 053A01 P. V. Prudnikov et al.
tRt
tt
t
t
t
Fig. 3. Time dependence of response function R(t,tw)as a function ttwfor a high-temperature initial state
with m01and for different values of tw.
t
ttCt ttRt
ttt
zz
(a) (b)
Fig. 4. Scaling collapse of correlation C(t,tw)(a)and response R(t,tw)(b)functions for a high-temperature
initial state with m01.
For a completely ordered initial state with m0=1, the short-time dynamics regime for autocorrela-
tion function C(t,tw)(Fig. 5) is characterized by the exponent φ(15) and we obtained the value φ=
2.742(32), which is in good agreement with the theoretical value φ=1+d/z+β/z)=2.737(8).
The obtained data for the ﬂuctuation-dissipation ratio for the case with m1are presented in
Fig. 6. The asymptotic value of the FDR, X=0.380(13), was obtained by linear extrapolation of
data for X(tw/(ttw)) to the limit tw/(ttw)0.
For a completely ordered initial state with m0=1, the asymptotic value of the FDR,
X=0.77(6), was obtained as a realization of the limit dependence Tcχ(C)(Fig. 7), in accor-
dance with relation (24). This value of Xis in excellent agreement with the theoretical ﬁeld value
X0.78 [50].
It should be noted that the non-equilibrium evolution of a system from a completely ordered initial
state is accompanied by a change in the magnetic domain structure from a single-domain initial
structure to a multidomain ﬂuctuation structure at the critical point. In diluted systems, this process is
additionally characterized by pinning of domain walls on defects. These phenomena have an essential
inﬂuence on the two-time dependence of the correlation and response functions and lead to a new
9/20
PTEP 2015, 053A01 P. V. Prudnikov et al.
tt
tt
tt
Ct
z
Fig. 5. Scaling collapse of correlation function C(t,tw)for a low-temperature initial state with m0=1.
tX
ttt
t
Fig. 6. Fluctuation-dissipation ratio X(t,tw)as a function of tw/(ttw)for ttwtw.
t
Tt
t
X
Ct
Fig. 7. Dependence of Tcχ(t,tw)on C(t,tw)for different values of twwith a demonstration of its deviation
from a line with X=1.
asymptotic value of the FDR. At present, our investigations into the inﬂuence of disorder on the non-
equilibrium dynamics of the 3D Ising model with evolution from a completely ordered initial state
are in progress.
10/20
PTEP 2015, 053A01 P. V. Prudnikov et al.
4. Monte Carlo simulations of 3D diluted Ising models with relaxation from
a high-temperature initial state
The Hamiltonian of the ferromagnetic Ising model diluted by nonmagnetic impurity atoms is given
by
H=−J
<i,j>
pipjSiSj,(25)
where random occupation numbers pitake the values 0 or 1, and piequals 1 if the site of the lattice
contains spin and 0 otherwise. We considered the cubic lattice with periodic boundary conditions.
Let us denote Ns=pL3as the number of spins in a lattice with linear size Lwhere pis the spin
concentration.
We also performed our Monte Carlo simulations with the use of the heat-bath updating rule for Ising
spin systems with spin concentrations p=0.8and 0.6 on a cubic lattice with L=128. We computed
the magnetization
M(t)= 1
Ns
Ns
i=1
piSi(t)av
,(26)
and the two-time autocorrelation function
C(t,tw)= 1
Ns
Ns
i=1
piSi(t)Si(tw)av
,(27)
where the square brackets [···]av stand for an additional average over disorder conﬁgurations. The
averaging of C(t,tw)was carried out on over 1000 samples with different disorder conﬁgurations
with 15 MC runs for each sample.
In this part of the investigation, we have simulated the non-equilibrium critical dynamics of the
systems starting from a high-temperature initial state only with magnetization value m0=0.01 for
p=0.8and m0=0.005 for p=0.6. The response function and ﬂuctuation-dissipation ratio were
calculated using the relations
R(t,tw)=1
Tc1
Ns
Ns
i=1Si(t)Si(tw+1)SW
i(tw+1)av
(28)
with SW
i=tanh(Jm=iSm/T),and
X(t,tw)=Ns
i=1Si(t)Si(tw+1)SW
i(tw+1)av
N
i=1Si(t)Si(tw+1)Si(tw)av
.(29)
The averaging of R(t,tw)and X(t,tw)was carried out on over 5000 samples with different disorder
conﬁgurations with 15 runs for each sample.
We give in Fig. 8the obtained curves for the time evolution of the magnetization M(t)from the
initial states with magnetization m0=0.02 at Tc=4.5114(1)for the pure Ising model, m0=0.01
at Tc=3.4995(2)for samples with spin concentration p=0.8,andm0=0.005 at Tc=2.4241(1)
for samples with p=0.6[25,26], which are characterized by a power-law dependence in the form
M(t)tθin the initial stage of evolution. The values of the exponent θcalculated with the use
of corrections to scaling are as follows: θ=0.106(4)for p=1,θ=0.127(16)for p=0.8,and
θ=0.186(39)for p=0.6[26,55,56]. The timescale of a critical initial increase of the magneti-
zation is tcr m1/(θ+β/zν)
0. The initial rise of magnetization is changed to the well known decay
11/20
PTEP 2015, 053A01 P. V. Prudnikov et al.
tM
t
p
p
p
Fig. 8. Time dependence of magnetization M(t)for different spin concentrations p.
ttC
ttR
tt tt
t
t
p
p
(a) (b)
Fig. 9. Time dependence of correlation C(t,tw)(a)and response R(t,tw)(b)functions from ttwfor
different spin concentrations p.
m(t)tβ/zνfor ttcr. However, in the limit of m00, the timescale tcr goes to inﬁnity. The
times with t<tcr can be considered as the time intervals of non-equilibrium behavior of the systems
considered.
The curves in Fig. 8demonstrate that time intervals of non-equilibrium behavior increase with
increasing dilution of spin systems. Therefore, if, for investigation of different non-equilibrium
regimes in the pure system with tcr 103MCS/s, we considered the values of waiting time tw<102
MCS/s, then, for diluted systems with tcr 104MCS/s, we can use the values of tw<103MCS/s.
The curves for the time dependence of the autocorrelation C(t,tw)and response R(t,tw)func-
tions are plotted in Fig. 9for different values of waiting time twand spin concentrations p,which
demonstrate that the aging effects are increased with increasing defect concentrations.
To check the scaling predictions for C(t,tw)and R(t,tw)given by relations (7), we plot the depen-
dences of t2β/(νz)
wC(t,tw)and t1+2β/(νz)
wR(t,tw)versus t/twin Fig. 10 with the use of exponent
values: z=2.191(21),2β/ν =1.016(32)for p=0.8[26]andz=2.663(30),2β/ν =0.924(80)
for p=0.6[55,56]. These functions demonstrate the collapse of curves for different twwith ﬁxed
spin concentration pinto a single curve with universal scaling dependence. Note that systems with
different pare characterized by different scaling functions FC,R(t,tw,p).
For an evolution stage with long time separations ttwtw1, the scaling functions
FC,R(t,tw,p)have power-law dependences (8), which are characterized by exponents caand cr.
12/20
PTEP 2015, 053A01 P. V. Prudnikov et al.
ttCt
tt tt
zttRt
z
pp
p
p
p
p
(a) (b)
Fig. 10. Scaling collapse of correlation C(t,tw)(a)and response R(t,tw)(b)functions for different spin
concentrations p.
Xtt
ttt
pt
t
t
t
t
t
t
pp
p
p
p
p
Fig. 11. Dependence of the FDR X(t,tw)on tw/(ttw)for ttwtwand for different p, where
p=1.0,tw=15 (); p=1.0,tw=30 (); p=0.8,tw=30 (); p=0.8,tw=50 (); p=0.6,tw=30 ();
p=0.6,tw=50 (); p=0.6,tw=150 ().
We determined the values of these exponents: ca=1.237(22),cr=1.251(22)for p=0.8and
ca=0.982(30),cr=0.950(8)for p=0.6. They demonstrate very good agreement with each
other for ﬁxed spin concentration p, but differences between values caand crfor different pand
ca=1.333(40)and cr=1.357(18)for the pure Ising model exceed the statistical errors of their
determination. These values of caare in good agreement with ca=1.242(10)for a weakly diluted
system with p=0.8and ca=0.941(21)for a strongly diluted system with p=0.6, obtained in
Refs. [26,55,56] by the short-time dynamics method, but are in poor agreement with ca=1.05(3),
determined in Ref. [57] as a characteristic, independent of p, of the non-equilibrium critical behavior
of the autocorrelation function. The reasons for this discrepancy have been discussed in detail in our
earlier paper [26].
In later investigations, we computed the ﬂuctuation-dissipation ratio in compliance with relation
(29). The data obtained are plotted in Fig. 11 for different spin concentrations p. For the analysis, we
used a time interval from ttwtwto ttwtw. In contrast to the pure Ising model, the data
for X(t,tw)for the diluted Ising model are characterized by an explicit dependence from tw.Atthe
beginning, we calculated the asymptotic values of X(tw,p)for different twfrom the plot in Fig. 11 in
the limit tw/(ttw)0, and then, using the obtained values of X(tw,p), we made an extrapolation
of 1/tw0to gain the asymptotic ﬂuctuation-dissipation ratios X(p). The results are presented
13/20
PTEP 2015, 053A01 P. V. Prudnikov et al.
Xt
p
p
p
t
Fig. 12. Dependence of the FDR X(tw,p)on 1/twfor different p. The asymptotic value of X(p)can be
obtained in the limit 1/tw0.
Tabl e 2 . Values of the FDR Xfor spin concentra-
tions p=1.0,0.8, and 0.6.
X
twp=1.0p=0.8p=0.6
10 0.361(2)
15 0.371(4)
20 0.365(10)0.373(9)
25 0.369(9)
30 0.374(14)0.384(5)0.382(1)
50 0.379(10)0.397(4)0.407(3)
100 0.406(6)0.427(6)
150 0.412(9)0.437(9)
→∞ 0.380(13)0.413(11)0.446(8)
in Fig. 12 and in Table 2. The ﬁnal values of the asymptotic ﬂuctuation-dissipation ratio for the
systems considered are X(p=1)=0.380(13),X(p=0.8)=0.413(11),andX(p=0.6)=
0.446(8).
The obtained value X(p=0.8)=0.413(11)for the 3D weakly diluted Ising model is in very
good agreement with the result X
3DRIM 0.416 from low-order renormalization-group computa-
tions by the ε-expansion method [48], but the value X(p=0.6)=0.446(8)for the 3D strongly
diluted Ising model does not agree with the above-denoted values in the limits of errors. Also, we
must mention that the obtained value X(p=1)=0.380(13)for the 3D pure Ising model does
not agree with X
3DIs =0.429(6), which has been calculated with the use of the ε-expansion method
in the two-loop approximation for pure systems with a non-conserved order parameter dynamics
[47], but it is characterized by close agreement with X(p=1)0.40, noted in Ref. [34]asthe
result of preliminary simulations on the 3D pure Ising model, without, however, demonstration of
any obtained data in either Ref. [34] or later publications.
Simulations with probing magnetic ﬁeld
Another way to compute the asymptotic ﬂuctuation-dissipation ratio is the application of random
probing magnetic ﬁelds hiwith small amplitude hin the simulation process after tw[58]. For this
14/20
PTEP 2015, 053A01 P. V. Prudnikov et al.
tt
ttt
ttt
ttt
p
p
p
t
tt
t
Fig. 13. The integrated response function time dependences χ(t,tw)for different values of twand different
spin concentrations p.
case, the Hamiltonian for the diluted Ising model, taking into account a local magnetic ﬁeld hi,is
given by
H=−J
<i,j>
pipjSiSj
i
hipiSi,(30)
and the two-time integrated response function χ(t,tw)is introduced through the relation
χ(t,tw)=t
tw
dtRt,t=
1
h2pL3
pL3
i=1
pihi(tw)Si(t)
av
,(31)
where the line stands for an average over the random ﬁeld realizations with hi=0and hihj=h2δij.
For simplicity, the random hiwere taken from a bimodal distribution (hih)with small amplitude
h=0.01 in order to avoid nonlinear effects.
In the large time limit,
Tcχ(C)=1
C
X(q)dq,
and the ﬂuctuation-dissipation ratio can be deﬁned as
X(t,tw)=−lim
C0
Tcχ(t,tw)
C(t,tw)(32)
through dependence on the integrated response function χ(t,tw)computed from the autocorrelation
function C(t,tw).
In this part of the investigation, we determined the time dependence of the integrated response
function using the Metropolis algorithm. The simulations were also performed on lattices with linear
size L=128 with spin concentrations p=1.0,0.8,and0.6. The spin systems considered were
quenched in the critical point from high-temperature initial states with m0=0.02 for pure samples,
m0=0.01 for samples with spin concentration p=0.8,andm0=0.005 for samples with p=
0.6. The averages were taken with the use of 5000 samples characterized by different independent
conﬁgurations of defects and 10 different realizations of random ﬁelds.
The computed time dependences for the integrated response function are plotted in Fig. 13 for
different values of waiting time twand spin concentrations p.
To obtain the FDR on the basis of relation (32), we analyzed the dependences Tcχ(t,tw)from
C(t,tw), found the slopes of curves for different tw, and then made an extrapolation, tw→∞.
15/20
PTEP 2015, 053A01 P. V. Prudnikov et al.
t
t
tCt
tCt
Ttt
Ttt
p
p
p
p
p
Fig. 14. Tcχ(t,tw)versus C(t,tw)for different p. The error bars are smaller than the symbols.
X
t
p
p
p
p
p
Fig. 15. Dependence of the FDR X(tw,p)on 1/twfor p=0.8and 0.6. The asymptotic value of X(p)can
be obtained in the limit 1/tw0. In the inset, the dependence of X(tw)on 1/twis given separately for the
pure Ising model.
Tabl e 3. The values of the FDR Xfor spin concen-
trations p=1.0,0.8, and 0.6computed with applica-
tion of a probing magnetic ﬁeld.
X
twp=1.0p=0.8p=0.6
10 0.586(24)
25 0.460(22)
50 0.437(26)
250 0.708(16)0.726(13)
500 0.553(17)0.583(14)
1000 0.494(17)0.519(29)
→∞ 0.391(12)0.419(11)0.443(10)
The results for these pure and diluted systems are given in Figs. 14 and 15 and also presented in
Tabl e 3. It should be noted that the values of X(tw), in accordance with relation (32), are calcu-
lated in the limit C(t,tw)0, which corresponds to the stage in the interval ttw1. The inset
16/20
PTEP 2015, 053A01 P. V. Prudnikov et al.
in Fig. 14 highlights the sections of the dependences of Tχ(t,tw)on C(t,tw), which satisfy these
criteria and on which the values of X(tw)were determined.
The comparison shows that the ﬁnal values of the asymptotic ﬂuctuation-dissipation ratio,
X(p=1.0)=0.391(12),X(p=0.8)=0.419(11),andX(p=0.6)=0.443(10), obtained
with application of a probing magnetic ﬁeld, are in very good agreement with the values
X(p=1)=0.380(13),X(p=0.8)=0.413(11),andX(p=0.6)=0.446(8), computed
with realization of heat-bath dynamics simulations.
5. Conclusions
In the introduction to this paper, we reviewed some theoretical results of computations in the range of
non-equilibrium phenomena that have been obtained in recent years for universal quantities, such as
the exponents determining the scaling behavior of dynamic response and correlation functions and
the ﬂuctuation-dissipation ratio, associated with the non-equilibrium critical dynamics. It was noted
that the value of the asymptotic ﬂuctuation-dissipation ratio Xfor quenches from the disordered
initial state with m01to the critical temperature Ts=Tcdepends on the universality class of
critical behavior to which the speciﬁc statistical system belongs. In the original part, we considered
the results of our MC simulations for 3D pure and diluted Ising models with Glauber and Metropolis
dynamics for quenches to their critical temperatures as functions of spin concentration. The inﬂuence
of critical ﬂuctuations, different non-equilibrium initial states, and site-quenched disorder on the two-
time dependence of the correlation and response functions in a non-equilibrium critical regime was
investigated for these models.
Analysis of the time dependences for the autocorrelation and response functions in the aging
regime of evolution showed that the aging effects increase with increasing defect concentration.
The obtained values of the autocorrelation and response function critical exponents ca,crand the
ﬂuctuation-dissipation ratio Xdemonstrate that statistical systems described by 3D pure and
diluted Ising models with a non-conserved order parameter dynamics belong to different classes
of non-equilibrium critical behavior.
Analysis of the simulation results shows that the insertion of disorder leads to new universal FDRs
with X
strong diluted >X
weak diluted >X
pure.
As mentioned above, the values computed for the 3D weakly diluted Ising model with p=0.8,
X=0.413(11)and 0.419(11), are in very good agreement with the result X
3DRIM 0.416
from low-order renormalization-group computations by the ε-expansion method [48], but, for
the 3D pure Ising model, the obtained values, X=0.380(13)and 0.391(12), do not agree with
X
3DIs =0.429(6), which has been calculated with the use of the ε-expansion method in the two-loop
approximation for pure systems with a non-conserved order parameter dynamics [47]. We think that
this difference of results is because the series of ε-expansion obtained in two-loop order (9) is not
an oscillating series. Therefore, it is impossible for summation of this series to apply the more exact
simplest Pade approximant method accounting for the direct [2,0] and inverse [0,2] approximants (the
[1,1] approximant has a pole in the important range). To obtain more reliable renormalization-group
results, one needs to calculate Xin the next three-loop order by the ε-expansion method or to
apply this more accurate method of ﬁeld-theoretic renormalization-group description directly to 3D
systems with the use of different resummation methods [24,54,59].
When discussing the inﬂuence of disorder on the universality of the critical behavior in the
site-diluted 3D Ising model, we must note that the results of Monte Carlo investigations into this
17/20
PTEP 2015, 053A01 P. V. Prudnikov et al.
model are quite contradictory. The results of some researchers are aimed at proving the concept
that the values of the critical exponents do not depend on the concentration of defects down to
the percolation threshold with ν=0.684(5),β=0.355(3),γ=1.342(10)[60], z=2.62(7)[61],
z=2.35(2)[62], and θ=0.10(2)[63] obtained by a speciﬁc procedure for ﬁtting intermediate
values of the effective exponents and amplitudes in the scaling dependence of the calculated thermo-
dynamic characteristics for different spin concentrations, using a ﬁtted exponent of the corrections
to the scaling ω=0.370(63)[60], ω1=0.50(13)[61], and ω2=0.82(8)[62]. The results of other
researchers indicate that there are two universal classes of critical behavior for weakly diluted sys-
tems with ν=0.68(2),β=0.34(2)[64], z=2.38(1)[65], ν=0.682(3),β=0.344(3)[66,67],
ν=0.683(4),β=0.310(3),γ=1.299(3)[68], ν=0.696(3),γ=1.345(4),ω=0.23(13)[25],
z=2.20(7)[14,15], z=2.191(21),ω=0.256(55),θ=0.127(16)[26], and z=2.28(7)[69],
and for strongly diluted systems with ν=0.72(2),β=0.33(2),γ=1.51(3)[64], z=2.53(3)
[65], ν=0.717(7),β=0.313(12)[66,67], ν=0.725(6),β=0.349(4),γ=1.446(4)[68], ν=
0.725(4),γ=1.415(11),ω=0.28(15)[25], z=2.58(9)[14,15], z=2.663(30),ω=0.286(10),
θ=0.167(18)[55,56], and z=2.67(8)[69].
In this discussion, we consider the results obtained in Ref. [70] very important, where the ﬁxed-
point structure of the 3D site- and bond-diluted and ±JIsing models using the numerical domain-
wall renormalization-group method has been studied. It was shown that the observed random ﬁxed
points are characterized by close values. This fact allows us to suggest that there exists a universal
ﬁxed point characterizing the 3D disordered ferromagnetic Ising model, irrespective of the type of
disorder. Unfortunately, this global conclusion was made on the basis of simulation results of systems
with very small linear sizes, L=8and 12. For the 3D site-diluted Ising model, the inﬂuence of the
correction to the ﬁnite-size scaling on the critical temperature, the values of the critical exponents, and
other critical characteristics is very important, especially for strongly diluted systems. In addition,
we must note that, in the renormalization-group ﬁeld-theory description of the random Ising model
(RIM), a model with two vertices characterizing the interaction of the order parameter ﬂuctuations is
applied [16,7173]. In this case, the critical behavior of RIM is determined by a stable ﬁxed point for
these two vertices with negative values for both exponents for the eigenvalues of the stability matrix.
In Ref. [70], the author uses a method with a single effective renormalization parameter r, responsible
for the inﬂuence of disorder and stability of a ﬁxed point (another parameter tis connected with the
effective temperature of the system). Such a transition from two-parametric phase space to one-
parametric space can be characterized by inaccuracy in the obtained results. However, the method
applied in Ref. [70] is very interesting, and it is desirable to see its development to 3D RIM with the
largest system sizes and with possible two-parametric modiﬁcation.
Acknowledgements
This work was supported by the Russian Scientiﬁc Fund through project No. 14-12-00562. The simula-
tions were supported by the Supercomputing Center of Lomonosov Moscow State University and Joint
Supercomputer Center of the Russian Academy of Sciences.
References
[1] L. Berthier and J. Kurchan, Nature Phys. 9, 310 (2013).
[2] E. Vincent, J. Hammann, M. Ocio, J.-P. Bouchaud, and L. F. Cugliandolo, Complex behavior of glassy
systems. In Lecture Notes in Physics (Springer, Berlin, 1997), Vol. 492, p. 184.
18/20
PTEP 2015, 053A01 P. V. Prudnikov et al.
[3] J.-P. Bouchaud, L. F. Cugliandolo, J. Kurchan, and M. Mezard, Spin glasses and random ﬁelds.
In Directions in Condensed Matter Physics, ed. A. P. Young (World Scientiﬁc, Singapore, 1998),
Vol. 12., p. 443.
[4] L. F. Cugliandolo, Slow relaxation and nonequilibrium dynamics in condensed matter. In Les Houches:
Ecole d’Et´
e de Physique Theorique, eds. J.-L. Barrat (Springer, Berlin, 2003), Vol. 77, p. 371.
[5] M. Henkel and M. Pleimling, Non-Equilibrium Phase Transitions, Volume 2: Ageing and Dynamical
Scaling far from Equilibrium (Theoretical and Mathematical Physics) (Springer, Heidelberg, 2010),
p. 544.
[6] N. Afzal and M. Pleimling, Phys. Rev. E 87, 012114 (2013).
[7] G. Parisi, F. Ricci-Tersenghi, and J. J. Ruiz-Lorenzo, Eur. Phys. J. B 11, 317 (1999).
[8] Y. Ozeki and N. Ito, J. Phys. A: Math. Theor. 40, R149 (2007).
[9] Y. Ozeki, S. Yotsuyanagi, and N. Ito, J. Phys. Soc. Jpn. 81, 074602 (2012).
[10] A. B. Harris, J. Phys. C: Solid State Phys. 7, 1671 (1974).
[11] D. E. Khmelnitskii, Sov. Phys. JETP 41, 981 (1976).
[12] G. Grinstein and A. Luther, Phys. Rev. B 13, 1329 (1976).
[13] V. V. Prudnikov and A. N. Vakilov, Sov. Phys. JETP 74, 990 (1992).
[14] V. V. Prudnikov and A. N. Vakilov, JETP Lett. 55, 741 (1992).
[15] V. V. Prudnikov and A. N. Vakilov, Sov. Phys. JETP 76 (3), 469 (1993).
[16] V. V. Prudnikov, S. V. Belim, A. V. Ivanov, E. V. Osintsev, and A. A. Fedorenko, JETP 87, 527 (1998).
[17] J. Q. Yin, B. Zheng, V. V. Prudnikov, and S. Trimper, Eur. Phys. J. B 49, 195 (2006).
[18] V. V. Prudnikov, P. V. Prudnikov, and A. A. Fedorenko, J. Phys. A: Math. Gen. 32, 8587 (1999).
[19] V. V. Prudnikov, P. V. Prudnikov, and A. A. Fedorenko, Phys. Rev. B 62, 8777 (1999).
[20] A. Pelissetto and E. Vicari, Phys. Rep. 368, 549 (2002).
[21] V. S. Dotsenko, A. B. Harris, D. Sherrington, and R. B. Stinchcombe, J. Phys. A: Math. Gen.
28, 3093 (1995).
[22] V. S. Dotsenko and D. E. Feldman, J. Phys. A: Math. Gen. 28, 5183 (1995).
[23] V. V. Prudnikov, P. V. Prudnikov, and A. A. Fedorenko, Phys. Rev. B 63, 184201 (2001).
[24] V. V. Prudnikov and P. V. Prudnikov, JETP 95, 550 (2002).
[25] V. V. Prudnikov, P. V. Prudnikov, A. N. Vakilov, and A. S. Krinitsyn, JETP 105, 371 (2007).
[26] V. V. Prudnikov, P. V. Prudnikov, A. S. Krinitsyn, A. N. Vakilov, E. A. Pospelov, and M. V. Rychkov,
Phys. Rev. E 81, 011130 (2010).
[27] P. C. Hohenberg and B. I. Halperin, Rev. Mod. Phys. 49, 435 (2010).
[28] P. Calabrese and A. Gambassi, J. Phys. A: Math. Gen. 38, R133 (2005).
[29] L. Berthier, P. C. W. Holdsworth, and M. Sellitto, J. Phys. A: Math. Gen. 34, 1805 (2001).
[30] A. Gambassi, Eur. Phys. J. B 64, 379 (2008).
[31] L. F. Cugliandolo, J. Kurchan, and L. Peliti, Phys. Rev. E 55, 3898 (1997).
[32] C. Godr´
eche and J. M. Luck, J. Phys.: Condens. Matter 14, 1589 (2002).
[33] L. F. Cugliandolo, J. Kurchan, and G. Parisi, J. Phys. (France) 4,1641 (1994).
[34] C. Godr´
eche and J. M. Luck, J. Phys. A: Math. Gen. 33, 9141 (2000).
[35] E. Lippiello and M. Zannetti, Phys. Rev. E 61, 3369 (2000).
[36] P. Mayer, L. Berthier, J. P. Garrahan, and P. Sollich, Phys. Rev. E 68, 016116 (2003).
[37] C. Chatelain, J. Phys. A: Math. Gen. 36, 10739 (2003).
[38] C. Chatelain, J. Stat. Mech. Theor. Exp. 2004, P06006 (2004).
[39] S. Abriet and D. Karevski, Eur. Phys. J. B 41, 79 (2004).
[40] H. K. Janssen, B. Schaub, and B. Schmittmann, Z. Phys. B 73, 539 (1989).
[41] V. V. Prudnikov, P. V. Prudnikov, I. A. Kalashnikov, and S. S. Tsirkin, JETP 106, 1095 (2008).
[42] K. Oerding and H. K. Janssen, J. Phys. A: Math. Gen. 28, 4271 (1995).
[43] K. Oerding and H. K. Janssen, J. Phys. A: Math. Gen. 26, 3369 (1993).
[44] K. Oerding and H. K. Janssen, J. Phys. A: Math. Gen. 26, 5295 (1993).
[45] A. Jaster, J. Mainville, L. Schulke, and B. Zheng, J. Phys. A: Math. Gen. 32, 1395 (1999).
[46] V. V. Prudnikov, P. V. Prudnikov, B. Zheng, S. V. Dorofeev, and V. Yu. Kolesnikov, Prog. Theor. Phys.
117, 973 (2007).
[47] P. Calabrese and A. Gambassi, Phys. Rev. E 66, 066101 (2002).
[48] P. Calabrese and A. Gambassi, Phys. Rev. B 66, 212407 (2002).
[49] R. Folk, Yu. Holovatch, and T. Yavorskii, Phys. Usp. 46, 169 (2003).
[50] P. Calabrese, A. Gambassi, and F. Krzakala, J. Stat. Mech. 6, 2 (2006).
19/20
PTEP 2015, 053A01 P. V. Prudnikov et al.
[51] W. Janke, Ageing and the glass transition. In Lecture Notes in Physics, eds. M. Henkel, M. Pleimling,
and R. Sanctuary (Springer, Berlin, 2007), Chap. 5, p. 207.
[52] F. Ricci-Tersenghi, Phys. Rev. E 68, 065104 (2003).
[53] R. Guida, and J. Zinn-Justin, J. Phys. A, Math. Gen. 31, 8103 (1998).
[54] V. V. Prudnikov, P. V. Prudnikov, and A. S. Krinitsyn, Theor. Math. Phys. 147, 559 (2006).
[55] V. V. Prudnikov, P. V. Prudnikov, A. N. Vakilov, A. S. Krinitsyn, and M. V. Rychkov, Proceedings of
workshop on computational technologies in natural sciences. In Computational Physics, ed. R. R.
Nazirov (KDU, Moscow, 2009), Vol. 1, p. 240 [in Russian].
[56] V. V. Prudnikov, P. V. Prudnikov, and A. N. Vakilov, Description of Nonequilibrium Critical Behavior in
Disordered Systems by the Theoretical Methods (Fizmatlit, Moscow, 2014) [in Russian].
[57] G. Schehr and R. Paul, Phys. Rev. E 72, 016105 (2005).
[58] A. Barrat, Phys. Rev. E 57, 3629 (1998).
[59] P. V. Prudnikov and V. V. Prudnikov, Phys. Rev. B 80, 024115 (2009).
[60] H. G. Ballesteros, L. A. Fern´
andez, V. Mart´
in-Mayor, M. Sudupe, G. Parisi, and J. J. Ruiz-Lorenzo,
Phys. Rev. B. 58, 2740 (1998).
[61] G. Parisi, F. Ricci-Tersenghi, and J. J. Ruiz-Lorenzo, Phys. Rev. E 60, 5198 (1999).
[62] M. Hasenbusch, A. Pelissetto, and E. Vicari, J. Stat. Mech.: Theory Exp. 2007, P11009 (2007).
[63] G. Schehr and R. Paul, J. Phys.: Conf. Ser. 40, 27 (2006).
[64] H.-O. Heuer, Phys. Rev. B 42, 6476 (1990).
[65] H.-O. Heuer, J. Phys. A: Math. Gen. 26, L341 (1993).
[66] S. Wiseman and E. Domany, Phys. Rev. Lett. 81, 22 (1998).
[67] S. Wiseman and E. Domany, Phys. Rev. E 58, 2938 (1998).
[68] A. K. Murtazaev, I. K. Kamilov, and A. B. Babaev, JETP 99, 1201 (2004).
[69] V. V. Prudnikov, A. N. Vakilov, and D. V. Talashok, JETP Lett. 100, 675 (2014).
[70] K. Hukushima, J. Phys. Soc. Jpn. 69, 631 (2000).
[71] G. Jug, Phys. Rev. B 27, 609 (1983).
[72] I. O. Mayer, J. Phys. A: Math. Gen. 22, 2815 (1989).
[73] A. Pelissetto and E. Vicari, Phys. Rev. B 62, 6393 (2000).
20/20
Article
Full-text available
Monte Carlo study of non-equilibrium critical behavior of three-dimensional Heisenberg model in isotropic case and with anisotropy of easy axis type is carried out. Relaxational Glauber-like dynamics of these models with evolution from high-temperature initial state is investigated. Realization of aging is demonstrated for two-time dependence of the autocorrelation function and dynamical susceptibility. Asymptotic fluctuation-dissipation ratios are determined for isotropic X ∞ = 0.383(6) and anisotropic X ∞ = 0.392(7) Heisenberg models. Significant influence of easy-axis anisotropy on non-equilibrium critical behavior of the 3D Heisenberg model leading to characteristics typical for the 3D Ising model is revealed.
Article
Full-text available
A Monte Carlo simulation of the non-equilibrium behavior of multilayer magnetic nanostructure Co/Cu(100)/Co consisting of alternating magnetic and nonmagnetic nanolayers is carried out. Analysis of calculated two-time autocorrelation function for structure relaxing from both high-temperature and low-temperature initial states reveals aging characterized by a slowing down of correlation characteristics with increase of the waiting time. It is shown that, in contrast to bulk magnetic systems, the aging effects in nanostructure arise not only at the ferromagnetic ordering temperature T c but also within a wide temperature range at T ≤ T c . For evolution from high-temperature initial state, the study of dependence of aging characteristics on thickness N of cobalt films reveals a weakening of the aging with increasing N at the critical temperatures T c (N) and an opposite tendency at temperatures T < T c (N) with strengthening of aging with increasing N for considered N ≤ 9 ML. This phenomenon is connected with increasing correlation and relaxation times in nanostructures when temperature is decreased. For case of the low-temperature initial state, it is shown that correlation times are two-three orders of magnitude smaller than those in the evolution from a high-temperature initial state at the same t w values. In this case, time behavior of the autocorrelation function doesn't depend considerably on temperature for T s ≤ T c and thickness N of cobalt films. Simulation of transport properties in Co/Cu(100)/Co structure permitted to calculate temperature dependence of its equilibrium magnetoresistance values. For the first time, it was revealed influence of non-equilibrium behavior on the magnetoresistance with demonstration of nontrivial aging effects. It has been shown that the magnetoresistance reaches plateau in asymptotic long-time regime with values , which depend on type of initial state, thickness of cobalt films, and temperature.
Article
The nonequilibrium behavior of the magnetic superstructures consisting of alternating magnetic and nonmagnetic nanolayers is numerically simulated by Monte Carlo methods. An analysis of the calculated two-time dependence of an autocorrelation function during the evolution from a high-temperature initial state has revealed aging effects, which are characterized by slowing down of the correlation effects in the system when the waiting time increases. In contrast to bulk magnetic systems, the aging effects are shown to appear in magnetic superstructures both near the critical ferromagnetic ordering temperature Tc in films and in a low-temperature phase at T ≤ Tc. The aging effects in the correlation processes in a magnetic multilayer structure weaken when ferromagnetic film thickness N increases at T = Tc(N), and these effects increase with film thickness N at temperatures T = Tc(N)/2. When simulating the transport properties of the Co/Cu(001)/Co structure, we calculated the temperature dependence of equilibrium magnetoresistance and were the first to reveal the influence of nonequilibrium behavior of the structure on the magnetoresistance and the manifestation of the aging effects in it.
Article
Full-text available
A Monte Carlo simulation of the non-equilibrium behavior of multilayer magnetic structures consisting of alternating magnetic and nonmagnetic nanolayers is realized. The calculated two-time autocorrelation function for the structure during its evolution starting from high-temperature initial state is analyzed. The analysis reveals aging effects characterized by a slowing down of correlation characteristics in the system with increase of the waiting time. The dependence of aging characteristics on thickness of ferromagnetic films is investigated. It is shown that, in contrast to bulk magnetic systems, the aging effects in magnetic superstructures arise not only at the ferromagnetic ordering temperature T c in the films but also within a wide temperature range at T ≤ T c .
Article
The effect of various initial magnetizations m0 and structural defects the nonequilibrium critical behavior of the three-dimensional Ising model is numerically studied. Based on an analysis of the time dependence of the magnetization and the two-time dependence of the autocorrelation function and dynamic susceptibility, the significant effect of initial states on relaxation magnetizations and aging effects characterized by anomalous relaxation inhibition and correlation in the system with increasing waiting time was revealed. The fluctuation–dissipation theorem violation was studied, and the values of the limit fluctuation–dissipation ratio (FDR) are calculated. It is shown that two universality subclasses can be distinguished in the nonequilibrium critical behavior of the three-dimensional Ising model with random initial magnetization m0 These subclasses correspond to the system evolution from the high-temperature (m0 = 0) and low-temperature (m0 = 1) initial states with limit FDRs characteristic of these states.
Article
Full-text available
Comprehensive Monte Carlo simulations of the short-time dynamic behaviour are reported for the three-dimensional Ising model at criticality. Besides the exponent 0305-4470/32/8/008/img7 of the critical initial increase and the dynamic exponent z, the static critical exponents 0305-4470/32/8/008/img8 and 0305-4470/32/8/008/img9 as well as the critical temperature are determined from the power-law scaling behaviour of observables at the beginning of the time evolution. States of very high temperature as well as of zero temperature are used as initial states for the simulations.
Article
Full-text available
In the past few years, systems with slow dynamics have attracted considerable theoretical and experimental interest. Ageing phenomena are observed during this everlasting non-equilibrium evolution. A simple instance of such a behaviour is provided by the dynamics that takes place when a system is quenched from its high-temperature phase to the critical point. The aim of this review is to summarize the various numerical and analytical results that have been recently obtained for this case. Particular emphasis is put on the field-theoretical methods that can be used to provide analytical predictions for the relevant dynamical quantities. Fluctuation-dissipation relations are discussed and in particular the concept of fluctuation-dissipation ratio (FDR) is reviewed, emphasizing its connection with the definition of a possible effective temperature. The renormalization-group approach to critical dynamics is summarized and the scaling forms of the time-dependent non-equilibrium correlation and response functions of a generic observable are discussed. From them, the universality of the associated FDR follows as an amplitude ratio. It is then possible to provide predictions for ageing quantities in a variety of different models. In particular, the results for models A, B and C dynamics of the O(N) Ginzburg-Landau Hamiltonian, and model A dynamics of the weakly dilute Ising magnet and of the phiv3 theory are reviewed and compared with the available numerical results and exact solutions. The effect of a planar surface on the ageing behaviour of model A dynamics is also addressed within the mean-field approximation.
Chapter
Results of Monte Carlo simulations in three and four spatial dimensions of a simple model that seems to have the necessary ingredients for disordered type-II superconductor behavior in an external magnetic field are reported. The data suggest that in d = 3 dimensions there is a finite temperature phase transition at T ≈ 0.45 into a truly superconducting vortex glass phase with infinite d.c. conductivity The (effective) correlation length exponent v and the dynamic critical exponent z at this transition are in good agreement with experiments. In d = 4 dimensions the gauge glass transition is located at T ≈ 0.95. It is concluded that the lack of time reversal symmetry in the model places it into a universality class different from that of the XY spin glass.
Article
The computer simulations of the critical dynamics of the structurally disordered three-dimensional Ising model are performed by the damage spreading method. For the systems with spin densities p = 0.6 and 0.8, we calculate the critical temperature and dynamic critical exponent z characterizing the behavior of relaxation properties near the critical point. The analysis of the results demonstrates the nonuniversality of the critical dynamics in the disordered Ising model. To interpret such dynamics, we introduce the concept of two universal types of critical behavior corresponding to weak and strong disorder, respectively.
Article
The ferromagnetic critical exponents are estimated along the phase boundary for the gauge glass model in three dimensions. The nonequilibrium relaxation method is applied to estimate the transition temperature and critical exponents. Together with the previously obtained results for the ± J Ising model in two and three dimensions, universality classes are examined for these disordered models. The result reveals a direct numerical confirmation of the Harris criterion for typical models with positive, zero and negative values of α0, the exponent for the specific heat in the regular case. In addition, we show a nonuniversal behavior of the dynamical exponent z with respect to the randomness for these models irrespective of α0.
Article
Ageing phenomena and dynamical scaling behaviour are studied both in ferromagnets and in critical spin glasses. For the ferromagnetic systems we find that dynamical scaling functions measured after a quench into the ordered low temperature phase are in complete agreement with predictions coming from the recently proposed theory of local scale invariance. For a quench to the critical point corrections to these predictions are s hown to exist. Critical spin glasses are found to display the same ageing phenomenology as critical ferromagnets. Numerical evi- dence indicates that the concept of universality is weaker i n critical spin glasses than in critical ferromagnets.


text:
Process affinity (or CPU pinning) means to bind each MPI process to a CPU or a range of CPUs on the node. It is important to spread MPI processes evenly onto different NUMA nodes.

Thread affinity means to map threads onto a particular subset of CPUs (called "places") that belong to the parent process (such as an MPI process) and to bind them to these places so the OS cannot migrate them to different places. It helps to take advantage of the local process state and to achieve better memory locality.

Memory locality is the degree to which data resides in memory that is close to the processors/threads working with the data.

Modern processors have multiple sockets and NUMA (Non-Uniform Memory Access) domains. Threads accessing memory in a remote NUMA domain is slower than accessing memory in a local NUMA domain.

Improper process and thread affinity could slow down code performance significantly. A combination of OpenMP environment variables and runtime flags are needed for different compilers and for the batch scheduler used on the system.

## Node architecture¶

### Cori Haswell¶

Figure 1 below illustrates a Haswell compute node. Each node contains 2 processors; There is 1 socket per processor, thus 2 sockets per node. Each processor has 16 cores, and each core has 2 hyperthreads. Socket 1 has physical cores 0 to 15, and socket 2 has physical cores 16 to 31. Core 0 has 2 hyperthreads, with the logical CPUs numbered as 0 and 32; Core 1 has logical CPUs of 1 and 33, and so on. When OMP_PLACES is set to "cores", each OpenMP thread binds to one core; And when OMP_PLACES is set to "threads", each OpenMP thread binds to one hyperthread.

Fig 1: CPUs, cores, and sockets on a Cori Haswell node.

Below is the numactl -H result from a Haswell compute node:

$numactl -H available: 2 nodes (0-1) node 0 cpus: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 node 0 size: 64430 MB node 0 free: 63002 MB node 1 cpus: 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 node 1 size: 64635 MB node 1 free: 63395 MB node distances:node 0 1 0: 10 21 1: 21 10  It reports that each Haswell compute node is composed of 2 NUMA "nodes" (i.e., NUMA domains 0 and 1). NUMA node 0 has logical CPUs 0-15 and 32-47. NUMA node 1 has logical CPUs 16-31 and 48-63. The tool also reports the total and available memory sizes on each NUMA domain and the relative distances from CPUs on each NUMA domain to access memory in both NUMA domains. Memory distance from NUMA domain 0 to 0 or 1 to 1 is closer than between NUMA domains 0 and 1. Note The values 10 and 21 are relative numbers; they do not necessarily mean that the ratio of the distances is exactly 21 to 10. ### Cori KNL¶ A Cori KNL node has 68 cores, 4 hyperthreads (cpus) per core and an additional level of hierarchy: each pair of cores share an L2 cache on a tile. The node is a single socket. Depending on the KNL mode a node is configured, its tiles can be organized as 1, 2 or 4 NUMA nodes. Moreover, the MCDRAM can be used as an invisible-to-the-OS cache, or as one or more NUMA nodes. Figure 2 below illustrates a KNL node in the "quad,cache" mode, in which all cpus share a single NUMA node with the DDR memory, and the MCDRAM is used as cache so is not considered as a separate NUMA node. (Note that tile numbering in the figure is illustrative, not accurate). Fig 2: CPUs, cores, and tiles on a Cori KNL node in "quad,cache" mode Below is the numactl -H result from a KNL "quad,cache" compute node: $ numactl -H
available: 1 nodes (0)
node 0 cpus: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76
77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101
102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121
122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141
142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161
162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181
182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201
202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221
222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241
242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261
262 263 264 265 266 267 268 269 270 271
node 0 size: 96527 MB
node 0 free: 93656 MB
node distances:
node   0
0:  10

Each KNL "quad,cache" compute node has only a single NUMA "node" (i.e., NUMA domains 0), which has logical CPUs 0 to 271, i.e., all the logical CPUs on the 68 physical cores. The KNL compute node has 96 GB of memory, and the 16GB MCDRAM is "invisible" as cache memory.

## Slurm Options¶

srun -n <total tasks> -c <logical CPUs per task> --cpu-bind <binding option> <executable>

The srun option -c sets the "number of logical cores (CPUs) per task" for the executable. For MPI and hybrid MPI/OpenMP applications one task corresponds to one MPI process.

Tip

The -c flag is optional for jobs that use one task per physical core.

### Cori¶

Tip

A general recommendation is to have at least one MPI rank per NUMA domain, and use OpenMP thread parallelism inside a NUMA domain so the data each thread works on are within the same NUMA domain.

#### Cori Haswell Slurm affinity¶

On Haswell, there are a total of 32 physical cores (each with 2 hyperthreads, so 64 logical cpus total), the value of -c should be set to 64 / tasks_per_node.

Example

To use 16 tasks per node, the -c value should be set to $64/16 = 4$.

In the general case when tasks per node is not a divisor of 64:

Example

To run with 12 tasks per node, the -c value should be set to $2*\left \lfloor{32/12}\right \rfloor = 4$.

Warn

srun option --cpu-bind=cores is needed when tasks per node is not a divisor of 64 on Haswell. In the case of more than 32 tasks per node (not recommended) use --cpu-bind=threads.

#### Cori KNL Slurm affinity¶

On Haswell, there are a total of 68 physical cores (each with 4 hyperthreads, so 272 logical cpus total), the value of -c should be set to 272 / tasks_per_node.

Example

To use 34 tasks per node, the -c value should be set to $272/34=8$.

In the general case when tasks per node is not a divisor of 272:

Example

To run with 8 tasks per node, the -c value should be set to $4*\left \lfloor{68/8}\right \rfloor = 32$.

Warn

srun option --cpu-bind=cores is needed when tasks per node is not a divisor of 272 on KNL. In the case of more than 68 tasks per node (not recommended) use --cpu-bind=threads.

### Perlmutter¶

#### AMD - EPYC 7763¶

The AMD EPYC processor has 64 physical cores each with 2 threads, for a total of 128 logical cpus.

The -c value should be:

#### GPUs¶

When allocating CPUs and GPUs to a job in Slurm, the default behavior is that all GPUs on a particular node allocated to the job can be accessed by all tasks on that same node:

// gpus_for_tasks.cpp
#include <iostream>
#include <string>
#include <cuda_runtime.h>
#include <mpi.h>

int main(int argc, char **argv) {
int deviceCount = 0;
int rank, nprocs;

MPI_Init (&argc, &argv);
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
MPI_Comm_size(MPI_COMM_WORLD, &nprocs);

cudaGetDeviceCount(&deviceCount);

printf("Rank %d out of %d processes: I see %d GPU(s).\n", rank, nprocs, deviceCount);

int dev, len = 15;
char gpu_id[15];

for (dev = 0; dev < deviceCount; ++dev) {
cudaSetDevice(dev);
cudaGetDeviceProperties(&deviceProp, dev);
printf("%d for rank %d: %s\n", dev, rank, gpu_id);
}

MPI_Finalize ();

return 0;
}

CC -cuda -o gpus_for_tasks gpus_for_tasks.cpp

$srun -C gpu -N 1 -n 2 -c 64 --cpu-bind=cores -G 4 ./gpus_for_tasks Rank 0 out of 2 processes: I see 4 GPUs. Their PCI Bus IDs are: 0 for rank 0: 0000:02:00.0 1 for rank 0: 0000:41:00.0 2 for rank 0: 0000:81:00.0 3 for rank 0: 0000:C1:00.0 Rank 1 out of 2 processes: I see 4 GPUs. Their PCI Bus IDs are: 0 for rank 1: 0000:02:00.0 1 for rank 1: 0000:41:00.0 2 for rank 1: 0000:81:00.0 3 for rank 1: 0000:C1:00.0  Likewise, if 1 GPU is needed for each task in a 2-task job using a node, the resulting 2 GPUs are visible to both tasks: $ srun -C gpu -N 1 -n 2 -c 64 --cpu-bind=cores --gpus-per-task=1 ./gpus_for_tasks
Rank 1 out of 2 processes: I see 2 GPUs. Their PCI Bus IDs are:
0 for rank 1: 0000:81:00.0
1 for rank 1: 0000:C1:00.0
Rank 0 out of 2 processes: I see 2 GPUs. Their PCI Bus IDs are:
0 for rank 0: 0000:81:00.0
1 for rank 0: 0000:C1:00.0

For some applications, it is desirable that only certain GPUs can be accessed by certain tasks. For example, a common programming model for MPI + GPU applications is such that each GPU on a node is accessed by only a single task on that node.

Such behavior can be controlled in different ways. One way is to manipulate the environment variable CUDA_VISIBLE_DEVICES, as documented in the Environment variables section of the CUDA C programming guide. This approach works on any system with NVIDIA GPUs. The variable must be configured per process, and may have different values on different processes, depending on the user's desired GPU affinity settings.

If the number of tasks is the same as the number of GPUs (4 on the GPU nodes) on each node and each task is to use one GPU only, a simple way of initializing CUDA_VISIBLE_DEVICES is to set it to $SLURM_LOCALID. The SLURM_LOCALID variable is the local ID for the task within a node. Since the local ID is defined after launching an srun command, you will need to wrap the environment variable setting as well as the executable run in one shell as in the following example: #!/bin/bash # select_cpu_device wrapper script export CUDA_VISIBLE_DEVICES=$SLURM_LOCALID
exec $*  $ srun -C gpu -N 1 -n 2 -c 64 --cpu-bind=cores --gpus-per-task=1 ./select_local_device ./gpus_for_tasks
Rank 1 out of 2 processes: I see 1 GPUs. Their PCI Bus IDs are:
0 for rank 1: 0000:C1:00.0
Rank 0 out of 2 processes: I see 1 GPUs. Their PCI Bus IDs are:
0 for rank 0: 0000:81:00.0

If CUDA MPS (Multi-Process Service) is being used and the number of tasks is larger than the number of GPUs on each node, you can modify the wrapper script so that GPUs can be evenly distributed in a round-robin fashion, unless it is already done by your application. Below is an example of the script for using 4 GPUs per node:

#!/bin/bash
# select_cpu_device wrapper script
export CUDA_VISIBLE_DEVICES=$(( SLURM_LOCALID % 4 )) exec$*

Another way to achieve a similar result is to use Slurm's GPU affinity flags. In particular, the --gpu-bind flag may be supplied to either salloc, sbatch, or srun in order to control which tasks can access which GPUs. A description of the --gpu-bind flag is documented here and via man srun. The --gpu-bind=map_gpu:<gpu_id_for_task_0>,<gpu_id_for_task_1>,... option sets the order of GPUs to be assigned to MPI tasks. For example, adding --gpu-bind=map_gpu:0,1 results in assigning GPUs 0 and 1 to MPI tasks in a round-robin fashion:

$srun -C gpu -N 1 -n 2 -c 64 --cpu-bind=cores --gpus-per-task=1 --gpu-bind=map_gpu:0,1 ./gpus_for_tasks Rank 0 out of 2 processes: I see 1 GPUs. Their PCI Bus IDs are: 0 for rank 0: 0000:81:00.0 Rank 1 out of 2 processes: I see 1 GPUs. Their PCI Bus IDs are: 0 for rank 1: 0000:C1:00.0  such that each task on the node may access a single, unique GPU. If a task needs to access only one GPU, there is a convenient option for --gpu-bind flag. With --gpu-bind=single:<tasks_per_gpu>, the first <tasks_per_gpu> tasks are bound to the first GPU available, the second <tasks_per_gpu> tasks are bound to the second GPU available, etc. The following is a hypothetical example for 6 MPI tasks on a node, sharing 4 GPUs: $ srun -C gpu -N 1 -n 6 -c 20 --cpu-bind=cores --gpus-per-node=4 --gpu-bind=single:2,2,1,1 ./gpus_for_tasks
Rank 4 out of 6 processes: I see 1 GPUs. Their PCI Bus IDs are:
0 for rank 4: 0000:81:00.0
Rank 1 out of 6 processes: I see 1 GPUs. Their PCI Bus IDs are:
0 for rank 1: 0000:02:00.0
Rank 3 out of 6 processes: I see 1 GPUs. Their PCI Bus IDs are:
0 for rank 3: 0000:41:00.0
Rank 2 out of 6 processes: I see 1 GPUs. Their PCI Bus IDs are:
0 for rank 2: 0000:41:00.0
Rank 5 out of 6 processes: I see 1 GPUs. Their PCI Bus IDs are:
0 for rank 5: 0000:81:00.0
Rank 0 out of 6 processes: I see 1 GPUs. Their PCI Bus IDs are:
0 for rank 0: 0000:02:00.0

To run a job across 8 GPU nodes using 32 tasks total, with each task bound to one GPU, one could run the following:

$srun -C gpu -N 8 -n 32 -c 32 --cpu-bind=cores --gpus-per-task=1 --gpu-bind=map_gpu:0,1,2,3 ./gpus_for_tasks Rank 31 out of 32 processes: I see 1 GPUs. Their PCI Bus IDs are: 0 for rank 31: 0000:C1:00.0 Rank 7 out of 32 processes: I see 1 GPUs. Their PCI Bus IDs are: 0 for rank 7: 0000:C1:00.0 Rank 3 out of 32 processes: I see 1 GPUs. Their PCI Bus IDs are: 0 for rank 3: 0000:C1:00.0 ... Rank 20 out of 32 processes: I see 1 GPUs. Their PCI Bus IDs are: 0 for rank 20: 0000:02:00.0 Rank 24 out of 32 processes: I see 1 GPUs. Their PCI Bus IDs are: 0 for rank 24: 0000:02:00.0 Rank 8 out of 32 processes: I see 1 GPUs. Their PCI Bus IDs are: 0 for rank 8: 0000:02:00.0  ## OpenMP Environment Variables¶ ### OMP_PLACES¶ OMP_PLACES defines a list of places that threads can be pinned on. It is used for complex layouts of threads. The possible values are: • threads: Each place corresponds to a single hardware thread on the target machine. • cores: Each place corresponds to a single core (having one or more hardware threads) on the target machine. • sockets: Each place corresponds to a single socket (consisting of one or more cores) on the target machine. • A list with explicit place values: such as: "{0,1,2,3},{4,5,6,7},{8,9,10,11},{12,13,14,15}” or “{0:4},{4:4},{8:4},{12:4}” can also be used. It has the following form "{lower-bound:length:stride}" Thus, specifying "{0:3:2}" is the same as specifying "{0,2,4}". Multiple locations can be included in a place. On Cori, it is recommended to set OMP_PLACES=threads or OMP_PLACES=cores. ### OMP_PROC_BIND¶ OMP_PROC_BIND sets the binding of threads to processors. The options are: • true: Thread affinity is enabled with an implementation-defined default place list. • false: Thread affinity is disabled. • spread: Bind threads as evenly distributed (spread) as possible. • close: Bind threads close to the master thread while still distributing threads for load balancing. • master: Bind threads to the same place as the master thread. On Cori, it is recommended to set OMP_PROC_BIND=true or OMP_PROC_BIND=spread. ### OMP_NUM_THREADS¶ OMP_NUM_THREADS sets the number of threads to be used for the OpenMP parallel regions. ## Job Script Generator¶ An interactive Job Script Generator is available at MyNERSC to provide some guidance on getting optimal process and thread binding on Cori Haswell and Cori KNL. It is a GUI tool that allows you to input the Slurm parameters of your intended application run, and get a job script template. ## Methods to Check Process and Thread Affinity¶ ### Use NERSC Prebuilt Binaries¶ Pre-built binaries from a small test code xthi.c with pure MPI or hybrid MPI/OpenMP can be used to check affinity. Binaries are in users default path, and named as such: check-mpi.<compiler>.<machine> (MPI), or check-hybrid.<compiler>.<machine> (MPI/OpenMP), for example: check-mpi.intel.cori, check-hybrid.intel.cori, check-mpi.gnu.cori, check-hybrid.gnu.cori, etc. The Cori binaries can be used to check both Haswell or KNL. It is recommended that you replace your application executable with one of the small test binaries, and run with the exact same number of nodes, MPI tasks, and OpenMP threads as those your application will use and check if the desired binding is obtained. Below is sample output (with interleaved notes) on Cori Haswell with an interactive batch session with various tests from different settings. The test results report that for each MPI rank, which node it runs on; and for each OpenMP thread that belongs to each MPI rank, which logical CPUs it binds to. #### MPI example¶ 8 MPI ranks on one node. Using the correct -n, -c, and --cpu-bind=cores options, the MPI tasks are spread out, and bind to both sockets on the Haswell node. MPI ranks 0,2,4,6 are on the first socket, and MPI ranks 1,3,5,7 are on the second socket. Each MPI rank binds to 4 physical CPUs (which has 8 logical CPUs total). For example, MPI rank 0 binds to physical cores 0-3, which includes logical CPUs 0-3,32-25. cori$ salloc -N 1 -C haswell --qos=interactive -t 20:00
salloc: Pending job allocation
21541802 salloc: job 21541802 queued and waiting for resources salloc:
job 21541802 has been allocated resources salloc: Granted job
allocation 21541802

elvis@nid00224:~> srun -n 8 -c 8 --cpu-bind=cores check-mpi.intel.cori
Hello from rank 0, on nid00224. (core affinity = 0-3,32-35)
Hello from rank 1, on nid00224. (core affinity = 16-19,48-51)
Hello from rank 2, on nid00224. (core affinity = 4-7,36-39)
Hello from rank 3, on nid00224. (core affinity = 20-23,52-55)
Hello from rank 4, on nid00224. (core affinity = 8-11,40-43)
Hello from rank 5, on nid00224. (core affinity = 24-27,56-59)
Hello from rank 6, on nid00224. (core affinity = 12-15,44-47)
Hello from rank 7, on nid00224. (core affinity = 28-31,60-63)

#### MPI/OpenMP example¶

8 MPI tasks, 4 OpenMP threads per MPI task. You could see the MPI tasks are again spread out on both sockets. Setting OMP_PLACES=threads will result in each thread bind to a specific hyperthread (a unique logical CPU as reported by Slurm). For GNU compiler, we set OMP_PROC_BIND=true, it spreads out the threads onto CPUs.

elvis@nid00224:~> export OMP_NUM_THREADS=4
elvis@nid00224:~> export OMP_PROC_BIND=true
elvis@nid00224:~> srun -n 8 -c 8 --cpu-bind=cores check-hybrid.gnu.cori |sort -k 4,6
Hello from rank 0, thread 0, on nid00224. (core affinity = 0)
Hello from rank 0, thread 1, on nid00224. (core affinity = 1)
Hello from rank 0, thread 2, on nid00224. (core affinity = 2)
Hello from rank 0, thread 3, on nid00224. (core affinity = 3)
Hello from rank 1, thread 0, on nid00224. (core affinity = 16)
Hello from rank 1, thread 1, on nid00224. (core affinity = 17)
Hello from rank 1, thread 2, on nid00224. (core affinity = 18)
Hello from rank 1, thread 3, on nid00224. (core affinity = 19)
Hello from rank 2, thread 0, on nid00224. (core affinity = 4)
Hello from rank 2, thread 1, on nid00224. (core affinity = 5)
Hello from rank 2, thread 2, on nid00224. (core affinity = 6)
Hello from rank 2, thread 3, on nid00224. (core affinity = 7)
Hello from rank 3, thread 0, on nid00224. (core affinity = 20)
Hello from rank 3, thread 1, on nid00224. (core affinity = 21)
Hello from rank 3, thread 2, on nid00224. (core affinity = 22)
Hello from rank 3, thread 3, on nid00224. (core affinity = 23)
Hello from rank 4, thread 0, on nid00224. (core affinity = 8)
Hello from rank 4, thread 1, on nid00224. (core affinity = 9)
Hello from rank 4, thread 2, on nid00224. (core affinity = 10)
Hello from rank 4, thread 3, on nid00224. (core affinity = 11)
Hello from rank 5, thread 0, on nid00224. (core affinity = 24)
Hello from rank 5, thread 1, on nid00224. (core affinity = 25)
Hello from rank 5, thread 2, on nid00224. (core affinity = 26)
Hello from rank 5, thread 3, on nid00224. (core affinity = 27)
Hello from rank 6, thread 0, on nid00224. (core affinity = 12)
Hello from rank 6, thread 1, on nid00224. (core affinity = 13)
Hello from rank 6, thread 2, on nid00224. (core affinity = 14)
Hello from rank 6, thread 3, on nid00224. (core affinity = 15)
Hello from rank 7, thread 0, on nid00224. (core affinity = 28)
Hello from rank 7, thread 1, on nid00224. (core affinity = 29)
Hello from rank 7, thread 2, on nid00224. (core affinity = 30)
Hello from rank 7, thread 3, on nid00224. (core affinity = 31)

#### MPI/OpenMP example 2¶

Below is another hybrid MPI/OpenMP example. 4 MPI tasks, 4 OpenMP threads per MPI task. You could see the MPI tasks are again spread out on both sockets. Setting OMP_PLACES=cores will result in each thread bind to a specific core (on Haswell, it has 2 logical CPUs on a physical core as reported by Slurm). For Intel compiler, we set OMP_PROC_BIND=spread.

elvis@nid00224:~> export OMP_NUM_THREADS=4
elvis@nid00224:~> export OMP_PLACES=cores
elvis@nid00224:~> srun -n 4 -c 16 --cpu-bind=cores check-hybrid.intel.cori |sort -k 4,6
Hello from rank 0, thread 0, on nid00224. (core affinity = 0,32)
Hello from rank 0, thread 1, on nid00224. (core affinity = 2,34)
Hello from rank 0, thread 2, on nid00224. (core affinity = 4,36)
Hello from rank 0, thread 3, on nid00224. (core affinity = 6,38)
Hello from rank 1, thread 0, on nid00224. (core affinity = 16,48)
Hello from rank 1, thread 1, on nid00224. (core affinity = 18,50)
Hello from rank 1, thread 2, on nid00224. (core affinity = 20,52)
Hello from rank 1, thread 3, on nid00224. (core affinity = 22,54)
Hello from rank 2, thread 0, on nid00224. (core affinity = 8,40)
Hello from rank 2, thread 1, on nid00224. (core affinity = 10,42)
Hello from rank 2, thread 2, on nid00224. (core affinity = 12,44)
Hello from rank 2, thread 3, on nid00224. (core affinity = 14,46)
Hello from rank 3, thread 0, on nid00224. (core affinity = 24,56)
Hello from rank 3, thread 1, on nid00224. (core affinity = 26,58)
Hello from rank 3, thread 2, on nid00224. (core affinity = 28,60)
Hello from rank 3, thread 3, on nid00224. (core affinity = 30,62)

### Use OpenMP 5.0 Environment Variables¶

A runtime thread affinity display feature is now included in the OpenMP 5.0 standard. This feature allows more portable and standard thread affinity display info across all compilers. It is available now for Intel compiler (18.0.5 and newer), and will be available soon for gcc (9.0 and newer) and CCE (9.0 and newer) on Cori:

The runtime display OpenMP thread affinity feature is enabled with two environment variables:

• OMP_DISPLAY_AFFINITY: set to TRUE or FALSE. Setting this to true will cause the system to display affinity information for all OpenMP threads when entering the first parallel region and when any thread affinity information changes in subsequent parallel regions.

• OMP_AFFINITY_FORMAT: set to a string that defines the output affinity values that will be output and the format used when displaying them.

The format for each field is: %[[[0].]size]type, where size defines the number of characters used for an output field and type indicates the information to output. The period indicates that the values are to be right-justified (the default is left-justified) and the 0 indicates that you want leading zeros to be included.

Some sample OMP_AFFINITY_FORMAT strings are:

OMP_AFFINITY_FORMAT="host=%H, pid=%P, thread_num=%n, thread affinity=%A"

OMP_AFFINITY_FORMAT="Thread Level=%0.3L, Parent TLevel=%5a, thread_num=%5n, thread_affinity=%15A, host=%10H”

Not specifying a size for each field allows it to expand the display result as needed. Notice this is an OpenMP only feature, so it does not have specific info on MPI ranks, etc. The host and pid fields are useful in identifying all the OpenMP threads belong to the same MPI task.

There are also omp_display_affinity() and omp_capture_affinity() APIs available that you could call from specific threads. Details for the available display affinity fields (with short and long names) and runtime APIs are in OpenMP 5.0 Specification and in the OpenMP 5.0 Examples document.

Below is sample output on Cori KNL (using intel/19.0.3.199, with a sample code hybrid-hello.f90 with the code output commented out). You could see that with OMP_PROC_BIND=spread, the OpenMP threads are evenly spread out, and with OMP_PLACES=threads, each thread binds to a specific hyperthread on a core; and with OMP_PLACES=cores, each thread binds to a specific core, and is allowed to freely migrate within the core on any of the hyperthreads (as reported as logical CPUs by Slurm. For example, logical CPUs 0,68,136,204 are all on physical core 0).

cori$module swap intel intel/19.0.3.199 cori$ ftn -qopenmp -o hybrid-hello hybrid-hello.f90

Then either in the batch script of in an interactive session set the appropriate variables before the srun command.

export OMP_NUM_THREADS=2
export OMP_DISPLAY_AFFINITY=true

$srun -n 8 -c 32 --cpu-bind=cores ./hybrid-hello |sort -k1,3 host=nid02496, pid=150147, thread_num=0, thread affinity=0 host=nid02496, pid=150147, thread_num=1, thread affinity=4 host=nid02496, pid=150148, thread_num=0, thread affinity=8 host=nid02496, pid=150148, thread_num=1, thread affinity=12 host=nid02496, pid=150149, thread_num=0, thread affinity=16 host=nid02496, pid=150149, thread_num=1, thread affinity=20 host=nid02496, pid=150150, thread_num=0, thread affinity=24 host=nid02496, pid=150150, thread_num=1, thread affinity=28 host=nid02496, pid=150151, thread_num=0, thread affinity=32 host=nid02496, pid=150151, thread_num=1, thread affinity=36 host=nid02496, pid=150152, thread_num=0, thread affinity=40 host=nid02496, pid=150152, thread_num=1, thread affinity=44 host=nid02496, pid=150153, thread_num=0, thread affinity=48 host=nid02496, pid=150153, thread_num=1, thread affinity=52 host=nid02496, pid=150154, thread_num=0, thread affinity=56 host=nid02496, pid=150154, thread_num=1, thread affinity=60  $ export OMP_PLACES=cores
\$ srun -n 8 -c 32 --cpu-bind=cores ./hybrid-hello |sort -k1,3

### Slurm cpu-bind flag¶

The srun flag --cpu-bind=verbose can be used to report process and thread binding. This option is recommended for advanced users only to interpret the "mask" output for CPUs on the node, such as:

cpu-bind=MASK - nid02496, task  0  0 [151947]: mask 0xff000000000000000ff000000000000000ff000000000000000ff set

### Compiler Specific Environment Variables¶

Or you can set the following runtime environment to obtain affinity information as part of the job stdout:

#### Intel¶

export KMP_AFFINITY=verbose

#### Cray CCE¶

export CRAY_OMP_CHECK_AFFINITY=TRUE

The results from each compiler have different formats and are not portable.

format: plain