# Nuclear configurational entropy and high energy hadron-hadron scattering
reactions

G. Karapetyan<EMAIL_ADDRESS>Federal University of ABC,
Center of Natural Sciences, Santo André, 09580-210, Brazil  Federal University
of ABC, Center of Mathematics, Santo André, 09580-210, Brazil

###### Abstract

In this work, the high energy hadron-hadron scattering is studied in the
framework of holographic AdS/QCD models, using the Brower-Polchinski-
Strassler-Tan Pomeron exchange kernel and gravitational form factors. We apply
the configurational entropy techniques to estimate the slope of the total
cross section for the total hadron-hadron cross sections at high energies. A
good agreement is derived between our approach and the total cross section for
combinations that include the pion-nucleon, nucleon-nucleon, and pion-pion, as
well as for any high energy data with inclusion of data from the TOTEM
collaboration at the LHC and approximated by the Pomeron exchange. In the case
of pion-nucleon and pion-pion scattering, the agreement for the critical
points to the differential configurational entropy can be reached within 1.1%
even without the involvement of any extra parameters.

##  I Introduction

The QCD theory of strong interactions is a gauge theory, wherein the Color
Glass Condensate (CGC) approach is a strong and convenient tool to investigate
the mechanism of high-energy processes. In this setup, the differential
configurational entropy (DCE) techniques predict the crucial points that help
to understand the massive experimental data obtained by different
collaborative groups daRocha:2021jzn ; daRocha:2021ntm ; Karapetyan:2018yhm ;
daRocha:2021imz ; Fernandes-Silva:2018abr ; Ferreira:2020iry ;
Bernardini:2019stn ; Bazeia:2021stz ; Bazeia:2021jok ; Silva:2021lmy ;
Bazeia:2021pbs ; Bazeia:2020jma ; Braga:2016wzx ; Braga:2021fey ;
daRocha:2021xwq ; Braga:2021zyi ; Braga:2020hhs . At high energy regimes, the
interesting and important information about various aspects of the reaction,
for example, the total hadron-hadron or photon-hadron cross sections, can be
derived by applying the nuclear DCE approach to nuclear systems at extreme
conditions Karapetyan:epjp ; Karapetyan:plb ; Karapetyan:2020epl ;
Karapetyan:2021epjp . During short time of investigation, the comprehensive
studies have been accumulated, in which such approach has been successfully
used to describe the mechanism of the nuclear reactions at high energies
Ferreira:2019inu ; Bernardini:2018uuy ; Braga:2018fyc ; Bernardini:2016hvx ;
Barbosa-Cendejas:2018mng ; Ferreira-Martins:2021cga , glueballs
Bernardini:2016qit , charmonium and bottomonium Braga:2017fsb , the quark-
gluon plasma daSilva:2017jay , baryons Colangelo:2018mrt , and also using the
AdS/QCD model approach Ma:2018wtw . The main details and theory of the CE
concept itself based on Shannon’s information entropy can be found in Refs.
Gleiser:2012tu ; Gleiser:2011di ; Gleiser:2013mga ; Gleiser:2018kbq ;
Sowinski:2015cfa ; Gleiser:2014ipa ; Gleiser:2015rwa . We can also mention
other works, where the critical points of the CE have been used to explain the
phenomena of AdS black holes and their quantum portrait as Bose–Einstein
graviton condensates, also emulating magnetic structures Casadio:2016aum ;
Fernandes-Silva:2019fez ; Braga:2020myi ; Alves:2017ljt ; Alves:2014ksa ;
Bazeia:2018uyg . The measure of the minimum values of CE, or the critical
points, carries all information content about the stability of the complex
nuclear systems. In the framework of the CGC, the nuclear DCE critical points
have been deeply studied by computing the inelastic hadron cross sections
Karapetyan:2018yhm ; Karapetyan:2019epl ; Ma:2018wtw ; Karapetyan:plb ;
Karapetyan:2017edu ; Karapetyan:2016fai . It is possible to compute the
nuclear DCE when the cross section plays the role of a localized function that
is related to the probability of the reaction product for any nuclear spatial
configuration. Hence, one can estimate the critical points of the DCE via the
Fourier transform of the cross section and subsequent calculation of the
related modal fraction Karapetyan:2018yhm ; Karapetyan:plb ;
Karapetyan:2017edu ; Karapetyan:2016fai .

Recently, the inner composition of hadrons, as well as the quark-gluon
pairing, caused a lot of interest because of the numerous experimentally
investigated scattering processes which have been done at high energy regime
as, for instance, the data taken at LHC. Within the framework of QCD, using
the factorization theorem, one can decompose the cross section at high
energies into the hard and the soft components. While the perturbation method
of QCD allows estimating the hard part of the cross section, it is not so
straightforward to calculate the soft part due to its not perturbative origin.
So, it is just possible to try to use some phenomenological approach, as the
parton distribution functions (PDFs) parametrization, which is represented by
the Bjorken scaling variable $x$ and the energy scale $Q^{2}$. At the high
energy regime of hadrons, the nonperturbative approach is one of the most
convenient methods to calculate the process cross section within QCD. In such
a context, the holographic QCD and the AdS/CFT correspondence have been used
to exam the high-energy scattering reactions. One should suggest the Pomeron
exchange (or a multi-gluon exchange) at the small $x$ range to be able to
characterize the whole dynamics of the partons. In Ref. Brower:2006ea  the
deep analysis based on the gauge/string duality, has been done by Brower,
Polchinski, Strassler, and Tan (BPST). This method assumes the contribution of
the Pomeron exchange to the cross sections, which can be caused by a kernel in
AdS space.

Hence, it is possible to use the BPST kernel approach to calculate the total
hadronic cross sections for various scattering processes at a high energy
regime. It is worth remembering that the masses of the participants in the
hadron-hadron scattering process are the only variables in the AdS/CFT
nonperturbative scale. In the AdS space within the BPST approach, the
estimation of the total hadronic cross sections implies the BPST Pomeron
exchange kernel and two density distributions. In the AdS/QCD approach the
density distributions are characterized by the gravitational form factors
Watanabe . For a non-excited state of hadron, one can use the hard-wall model
to describe its behavior. Such a model is valid for AdS geometry in the QCD
scale, where there is a cut-off in the infrared region (IR). In the case of
fixed model parameters based on the given experimental data, the cross section
for hadron scattering reaction can be easily calculated with the assumption of
hadron normalizable mode and its normalized density distribution. Hence, we do
not need to use any other parameters and get a universal description of any
high-energy process within the Pomeron exchange approximation.

We use the DCE paradigm to investigate the pion-nucleon, nucleon-nucleon, and
pion-pion cross section using the Pomeron exchange kernel and gravitational
form factors. Section II is devoted to the Pomeron exchange model in the
framework of holographic QCD for high-energy hadronic interactions. In the
next section, we present the main results of the analysis using DCE approach
and the concept of the critical points based on the experimentally determined
hadron-hadron cross sections. The summary and outlook are given in Sec. IV.

##  II The total hadron-hadron cross sections within the holographic QCD

Let us briefly represent the model description in the framework of the PST
Pomeron exchange kernel, which allows computing the hadron scattering total
cross section ($\chi$). Within this model the scattering amplitude in eikonal
formulation is given by the following expression:

$\displaystyle{\cal A}(s,t)=2is\int d^{2}be^{i\bm{k_{\perp}}\cdot\bm{b}}\int
dzdz^{\prime}P_{13}(z)P_{24}(z^{\prime})\left[1-e^{i\chi(s,\bm{b},z,z^{\prime})}\right].$
(1)

In Eq. (1) the Mandelstam variables are denoted by $s$ and $t$, whereas the
two-dimensional impact parameter is denoted by $\bm{b}$. Finally, $z$ and
$z^{\prime}$ are coordinates in AdS for the incident and target particles,
respectively, and $P_{13}(z)$ and $P_{24}(z^{\prime})$ are the density
distributions of the two hadrons in the AdS space, which are normalized
according to the following conditions:

$\int dzP_{13}(z)=\int dz^{\prime}P_{24}(z^{\prime})=1,$ (2)

Taking into account the optical theorem, in the single-Pomeron exchange model
the total cross section can be expressed as Watanabe

$\sigma(s)=2\int d^{2}b\int
dzdz^{\prime}P_{13}(z)P_{24}(z^{\prime})\Im\chi(s,\bm{b},z,z^{\prime}).$ (3)

The final expression for the total hadron-hadron cross section can be written
as Watanabe

$\displaystyle\sigma(s)=\frac{g_{0}^{2}\rho^{3/2}}{8\sqrt{\pi}}\int
dzdz^{\prime}P_{13}(z)P_{24}(z^{\prime})(zz^{\prime})\Im[\chi_{c}(s,z,z^{\prime})],$
(4)
$\displaystyle\Im[\chi_{c}(s,z,z^{\prime})]\equiv\frac{1}{\sqrt{\tau}}\exp\left[{(1-\rho)\tau-\frac{1}{\rho\tau}\log^{2}\left(\frac{z}{z^{\prime}}\right)}\right],$
(5)

where ${Im}[\chi_{c}(s,z,z^{\prime})]$ i s the imaginary part of the
scattering amplitude and $\tau=\log(\rho zz^{\prime}s/2)$ with $g_{0}^{2}$ and
$\rho$ are the overall factor and the slope of the total cross section,
respectively. It should be mentioned that in the nucleon-nucleon scattering
the main value that describes the process is the nucleon mass ($m_{N}\sim 1$
GeV). This fact supports the idea the leading role of the the low energy
dynamics with a strong coupling in QCD. It means that one can use the modified
BPST kernel with the same functional form as the conformal kernel in the
following form,

$\displaystyle\Im[\chi_{mod}(s,z,z^{\prime})]=\Im[\chi_{c}(s,z,z^{\prime})]$
$\displaystyle\qquad\qquad\qquad\qquad+\mathcal{F}(s,z,z^{\prime})\Im[\chi_{c}(s,z,z_{0}z_{0}^{\prime}/z^{\prime})],$
(6)

where

$\displaystyle\mathcal{F}(s,z,z^{\prime})$
$\displaystyle=1-2\sqrt{\rho\pi\tau}e^{\eta^{2}}\mbox{erfc}(\eta),$ (7)
$\displaystyle\eta$
$\displaystyle=\frac{-\log\left(\frac{zz^{\prime}}{z_{0}z_{0}^{\prime}}\right)+\rho\tau}{\sqrt{\rho\tau}},$
(8)

with the coordinates cut-off parameters, $z_{0}$ and $z^{\prime}_{0}$, in QCD
scale, fixed by hadron masses.

The density distributions, $P_{13}(z)$ and $P_{24}(z^{\prime})$, of the
interacting hadrons in Eq. (1) are described by the gravitational form
factors. Such form factors can be received from the bottom-up AdS/QCD models
using the hadron-Pomeron-hadron three-point functions. As the nucleon can be
expressed as a solution to the five-dimensional Dirac equation Watanabe , the
density distribution of the nucleon is expressed with the left-handed and
right-handed components of the Dirac field ($\psi_{L}$ and $\psi_{R}$) via the
Bessel function in the following form:

$\displaystyle
P_{N}(z)=\frac{1}{2z^{3}}\left[\psi_{L}^{2}(z)+\psi_{R}^{2}(z)\right],$ (9)
$\displaystyle\psi_{L}(z)=\frac{\sqrt{2}z^{2}J_{2}(m_{N}z)}{z_{0}^{N}J_{2}(m_{N}z_{0}^{N})},\
\psi_{R}(z)=\frac{\sqrt{2}z^{2}J_{1}(m_{N}z)}{z_{0}^{N}J_{2}(m_{N}z_{0}^{N})},$
(10)

The fixed cut-off parameter ($z_{0}^{N}=4.081keV$) obeys the conditions of the
normalization $J_{1}(m_{N}z_{0}^{N})=0$ with $m_{p}$ and $m_{n}$ being the
proton and neutron physical masses, respectively. The pion wave function
$\Psi$ is extracted from the bottom-up AdS/QCD model motion equation of mesons
Watanabe . Hence, we can present the density distribution of the pion in the
form of the Bessel function as the following:

$\displaystyle
P_{\pi}(z)=\frac{\left[\Psi^{\prime}(z)\right]^{2}}{4\pi^{2}f_{\pi}^{2}z}+\frac{\sigma^{2}z^{6}\Psi(z)^{2}}{f_{\pi}^{2}z^{3}},$
(11)

where the prime denotes the derivative with respective to $z$ and the pion
wave function reads

$\displaystyle\Psi\left(z\right)=z\Gamma\left(\frac{2}{3}\right)\left({\frac{\alpha}{2}}\right)^{1/3}\Biggl{[}I_{-1/3}\left({\alpha
z^{3}}\right)-I_{1/3}\left({\alpha
z^{3}}\right)\frac{{I_{2/3}\left({\alpha(z_{0}^{\pi})^{3}}\right)}}{{I_{-2/3}\left({\alpha(z_{0}^{\pi})^{3}}\right)}}\Biggr{]}$
(12)

with $z_{0}^{\pi}=3.105{\rm keV}$ and $\alpha=\frac{2\pi}{3}\sigma$, for
$\sigma=(332$ MeV$)^{3}$, and $f_{\pi}$ is the pion decay constant. The cutoff
parameter, $z_{0}^{\pi}$, encodes the mass of $\rho$ meson ($m_{\rho}$) by the
zeroes of the zeroth order Bessel function as $J_{0}(m_{\rho}z_{0}^{\pi})=0$.

##  III Configurational entropy via the BPST Pomeron exchange kernel

To determine the density distribution, one does not need any additional
adjustable parameter. Within the BPST Pomeron exchange kernel approach, the
only two parameters that should be extracted from the experiments are
$g_{0}^{2}$ and $\rho$ (see Eq. (5)). During our analysis, we use the
experimental data of TOTEM collaboration at LHC Watanabe , as well as the
recent hadron-hadron collision data which are given in the Particle Data Group
within the energy range $10^{2}<\sqrt{s}<10^{5}$ GeV. As the check point
values for the two adjustable parameters have been used the data from Ref.
Watanabe , $g_{0}^{2}=6.27\times 10^{2}$ and $\rho=0.824$.

We use the Pomeron exchange kernel to calculate the total hadron-hadron
scattering cross sections at a high energy range. As the using adjustable
parameters of the model do not depend on the intrinsic properties of hadrons,
one can assume the Pomeron exchange approach is a universal tool to describe
the process at high energies. This fact is supporting also by the chosen
normalizable modes of both the nucleon and the pion. In this way, one can
easily compute the cross section for the incident pion, if the parameters for
the nucleon-nucleon interaction are already fixed, and do not include other
adjustable parameters. Besides, the scale of the scattering processes at high
energies is characterized by the constant masses of the participated hadrons.
Therefore, in such a case, we also have the same conditions for all three
processes under study. which leads to the constant value for the total cross
section ratios among them. In this paper, we use the normalized structure
factor for the collective coordinates, which have been obtained by the Fourier
function decomposition into the number of weighted components for the total
hadron-hadron cross sections Bernardini:2016hvx . The cross section for any
collision, as the probability of the nuclear reaction, is governed by the
modal fraction of the CE, and, consequently, contains all the needed
information. On the other hand, the limit number of the critical points that
can be determined using the Shannon entropy and the configurational entropy
concepts, totally characterize an isolated localized nuclear system
Gleiser:2012tu ; Gleiser:2013mga .

If one considers as a localized function the reaction cross section, it is
possible to calculate the latter by via the Fourier transformation
energy–weighted correlation function, given by the cross section,

$\sigma({k,\rho})=\frac{1}{2\pi}\\!\int_{\mathbb{R}}\sigma(z,\rho)\,e^{ikz}dz.$
(13)

Then, the modal fraction of CE can be determined as:

$f_{\sigma({k,\rho})}=\frac{|\sigma({k,\rho})|^{2}}{\int_{\mathbb{R}}|\sigma({k,\rho})|^{2}dk}.$
(14)

The set of the critical points is determined through the suitable expression
for the DCE Gleiser:2012tu :

${\rm CE}(a,b)=-\int_{\mathbb{R}}f_{\sigma({k,\rho})}\log
f_{\sigma({k,\rho})}dk.$ (15)

And finally, one obtains the CE critical points just by applying Eqs. (13 \-
15) employing the total hadronic cross section.

Figure 1: DCE $\times$ $\rho$. Minimum: $\rho=0.814$, DCE = 164.72 nat
(Natural unit of information.)

In Fig. 1 one can see the results of the DCE versus $\rho$ calculations for
the total hadron-hadron scattering cross sections at a high energy regime.
Fig. 1 shows the predominance of the quantum states with the slope of the
total cross section $\rho=0.814$, corresponding to the nuclear DCE = 164.72
nat, with the corresponding global minimum of the nuclear configurational
entropy at such point. One can see that the differential configurational
entropy curve tends to grow on both sides from the observed global minimum,
although with different signs when one looks at the second derivative of the
nuclear DCE with respect to $\rho$, in Fig. 1. The result of our calculation
for hadron-hadron interactions at high energies fits well with the value of
parameter $\rho=0.824$ obtained in Ref. Watanabe , within 1.1%. The critical
point of the nuclear DCE indicates the natural choice of the total hadron
cross section for the parameter, which controls the energy dependence of the
cross sections and confirms the earlier received experimental data with high
enough precision. At $\rho=0.814$, the nuclear system stays at its most stable
configuration. Such minimum has been computed through Shannon’s information
entropy approximation Gleiser:2011di  and contains all information about the
mechanism of the reactions. The fitting procedure of adjustable parameter for
the experimentally determined total hadron-hadron cross sections clearly shows
the conditions for the configurational stability of the excited nuclear system
at high energies. The stability of any localized nuclear configuration is
characterized by the set of limited number degrees of freedom, which determine
the most dominant state(s) after reaching the thermodynamic equilibrium of the
system. And our results once again confirm this fact. Of course, each case of
the reaction requires special and deeper investigation using the comprehensive
results obtained via the recent studies in this area Bazeia:2013usa ;
Correa:2015lla ; Correa:2015vka

##  IV Conclusions

The total cross sections for three hadron-hadron scattering reactions at a
high energy regime have been studied within the Pomeron exchange kernel model
of holographic QCD MarinhoRodrigues:2020yzh ; MarinhoRodrigues:2020ssq . This
model allowed us to calculate the total hadron cross sections using the
expression for the density distributions for the nucleon and for pion, which
are expressed by the gravitational form factors in the framework of bottom-up
AdS/QCD models. The model conditions imply complete independence from the
observed $s$ quantity because of the energy dependence of the cross section
controlled by the BPST kernel. In the considering hadron energy range
($10^{2}<\sqrt{s}<10^{5}$ GeV) the total scattering cross sections are in
satisfactory agreement with the data from Ref. Watanabe . We can conclude
based on the results of this paper that the BPST Pomeron exchange kernel is a
convenient tool in the investigation of high energy hadron interactions with
nonperturbative gluonic dynamics. We applied the nuclear DCE approach and
estimated its global minimum, which characterizes the critical and stable
point of the hadron-hadron reactions. The calculation of the total scattering
cross section within the BPST Pomeron exchange kernel determined the natural
choice for the excited nuclear configuration. We found the satisfactory
agreement of the calculated parameters with those predicted by Watanabe ,
within 1.1%. The stability of the nuclear system can be completely determined
by the critical point of the nuclear DCE. Hence, all the available information
can be reached through the systematic study of the DCE approach.

### Acknowledgments:

GK thanks to The São Paulo Research Foundation – FAPESP (grant No.
2018/19943-6).

## References

  * (1) R. da Rocha, Phys. Lett. B 823, 136729 (2021) [arXiv:2108.13484 [gr-qc]].
  * (2) R. da Rocha, Phys. Rev. D 103, 106027 (2021) [arXiv:2103.03924 [hep-ph]].
  * (3) G. Karapetyan, Phys. Lett. B 786, 418 (2018) [arXiv:1807.04540 [nucl-th]].
  * (4) R. da Rocha, Phys. Lett. B 814, 136112 (2021) [arXiv:2101.03602 [hep-th]].
  * (5) A. Fernandes-Silva, A. J. Ferreira-Martins and R. da Rocha, Eur. Phys. J. C 78, 631 (2018) [arXiv:1803.03336 [hep-th]].
  * (6) L. F. Ferreira and R. da Rocha, Phys. Rev. D 101, 106002 (2020) [arXiv:2004.04551 [hep-th]].
  * (7) A. E. Bernardini and R. da Rocha, Phys. Lett. B 796, 107 (2019) [arXiv:1908.04095 [gr-qc]].
  * (8) D. Bazeia and E. I. B. Rodrigues, Phys. Lett. A 392, 127170 (2021).
  * (9) D. Bazeia, D. A. Ferreira and M. A. Marques, [arXiv:2102.06932 [hep-th]].
  * (10) A. A. Silva, F. M. Andrade and D. Bazeia, [arXiv:2101.05250 [quant-ph]].
  * (11) D. Bazeia, M. A. Liao and M. A. Marques, [arXiv:2101.01691 [hep-th]].
  * (12) D. Bazeia, D. A. Ferreira, F. S. N. Lobo and J. Luís, Eur. Phys. J. Plus 136, 321 (2021) [arXiv:2011.06240 [gr-qc]].
  * (13) N. R. F. Braga, Y. F. Ferreira and L. F. Ferreira, [arXiv:2110.04560 [hep-th]].
  * (14) R. da Rocha, [arXiv:2111.01244 [hep-th]].
  * (15) N. R. F. Braga and O. C. Junqueira, Phys. Lett. B 820, 136485 (2021) [arXiv:2105.12347 [hep-th]].
  * (16) N. R. F. Braga and R. da Mata, Phys. Lett. B 811, 135918 (2020) [arXiv:2008.10457 [hep-th]].
  * (17) N. R. F. Braga and R. da Rocha, Phys. Lett. B 767, 386 (2017) [arXiv:1612.03289 [hep-th]].
  * (18) G. Karapetyan, Eur. Phys. J. Plus 136, 1012 (2021) [arXiv:2105.07546 [hep-ph]].
  * (19) G. Karapetyan, Phys. Lett. B 781, 201 (2018) [arXiv:1802.09105 [hep-ph]].
  * (20) G. Karapetyan, EPL 129, 18002 (2020) [arXiv:1912.10071 [hep-ph]].
  * (21) G. Karapetyan, Eur. Phys. J. Plus 136, 122 (2021) [arXiv:2003.08994 [hep-ph]].
  * (22) L. F. Ferreira and R. da Rocha, Phys. Rev. D 99, 086001 (2019) [arXiv:1902.04534 [hep-th]].
  * (23) A. E. Bernardini and R. da Rocha, Phys. Rev. D 98, 126011 (2018) [arXiv:1809.10055 [hep-th]].
  * (24) N. R. F. Braga, L. F. Ferreira and R. da Rocha, Phys. Lett. B 787, 16 (2018) [arXiv:1808.10499 [hep-ph]].
  * (25) A. E. Bernardini and R. da Rocha, Phys. Lett. B 762, 107 (2016) [arXiv:1605.00294 [hep-th]].
  * (26) N. Barbosa-Cendejas, R. Cartas-Fuentevilla, A. Herrera-Aguilar, R. R. Mora-Luna and R. da Rocha, Phys. Lett. B 782, 607 (2018) [arXiv:1805.04485 [hep-th]].
  * (27) A. J. Ferreira-Martins and R. da Rocha, Nucl. Phys. B 973, 115603 (2021) [arXiv:2104.02833 [hep-th]].
  * (28) A. E. Bernardini, N. R. F. Braga and R. da Rocha, Phys. Lett. B 765, 81 (2017) [arXiv:1609.01258 [hep-th]].
  * (29) N. R. F. Braga and R. da Rocha, Phys. Lett. B 776, 78 (2018) [arXiv:1710.07383 [hep-th]].
  * (30) A. Goncalves da Silva and R. da Rocha, Phys. Lett. B 774, 98 (2017) [arXiv:1706.01482 [hep-th]].
  * (31) P. Colangelo and F. Loparco, Phys. Lett. B 788, 500 (2019) [arXiv:1811.05272 [hep-ph]].
  * (32) C. W. Ma and Y. G. Ma, Prog. Part. Nucl. Phys. 99, 120 (2018) [arXiv:1801.02192 [nucl-th]].
  * (33) M. Gleiser and N. Stamatopoulos, Phys. Rev. D 86, 045004 (2012) [arXiv:1205.3061 [hep-th]].
  * (34) M. Gleiser and N. Stamatopoulos, Phys. Lett. B 713, 304 (2012).
  * (35) M. Gleiser and D. Sowinski, Phys. Lett. B 727, 272 (2013) [arXiv:1307.0530 [hep-th]].
  * (36) M. Gleiser, M. Stephens and D. Sowinski, Phys. Rev. D 97, 096007 (2018) [arXiv:1803.08550 [hep-th]].
  * (37) M. Gleiser and D. Sowinski, Phys. Lett. B 747, 125 (2015) [arXiv:1501.06800 [hep-th]].
  * (38) M. Gleiser and N. Graham, Phys. Rev. D 89, 083502 (2014) [arXiv:1401.6225 [hep-th]].
  * (39) M. Gleiser and N. Jiang, Phys. Rev. D 92, 044046 (2015) [arXiv:1506.05722 hep-th]].
  * (40) R. Casadio and R. da Rocha, Phys. Lett. B 763, 434 (2016) [arXiv:1610.01572 [hep-th]].
  * (41) A. Fernandes-Silva, A. J. Ferreira-Martins and R. da Rocha, Phys. Lett. B 791, 323 (2019) [arXiv:1901.07492 [hep-th]].
  * (42) N. R. F. Braga and R. da Mata, arXiv:2002.09413 [hep-th].
  * (43) A. Alves, A. G. Dias and R. Silva, Braz. J. Phys. 47, 426 (2017) [arXiv:1703.02061 [hep-ph]].
  * (44) A. Alves, A. G. Dias and R. da Silva, Physica A 420, 1 (2015) [arXiv:1408.0827 [hep-ph]].
  * (45) D. Bazeia, D. C. Moreira and E. I. B. Rodrigues, J. Magn. Magn. Mater. 475, 734 (2019) [arXiv:1812.04950 [cond-mat.mes-hall]].
  * (46) G. Karapetyan, EPL 125, 58001 (2019) [arXiv:1901.05349 [hep-ph]].
  * (47) G. Karapetyan, EPL 118, 38001 (2017) [arXiv:1705.1061 [hep-ph]].
  * (48) G. Karapetyan, EPL 117, 18001 (2017) [arXiv:1612.09564 [hep-ph]].
  * (49) R. C. Brower, J. Polchinski, M. J. Strassler, C.- I. Tan, JHEP 0712, 005 (2007).
  * (50) A. Watanabe and M. Huang, Phys. Lett. B 788, 256 (2019) [arXiv:1809.02515 [hep-ph]].
  * (51) D. Bazeia, R. Menezes and R. da Rocha, Adv. High Energy Phys. 2014, 276729 (2014) [arXiv:1312.3864 [hep-th]].
  * (52) R. A. C. Correa and R. da Rocha, Eur. Phys. J. C 75, 522 (2015) [arXiv:1502.02283 [hep-th]].
  * (53) R. A. C. Correa, R. da Rocha and A. de Souza Dutra, Annals Phys. 359, 198 (2015) [arXiv:1501.02000 [hep-th]].
  * (54) D. Marinho Rodrigues and R. da Rocha, Phys. Lett. B 811, 135943 (2020) [arXiv:2009.01890 [hep-th]].
  * (55) D. Marinho Rodrigues and R. da Rocha, [arXiv:2006.00332 [hep-th]].


# Calculating phase diagrams with ATAT

Martin Bäker  Technische Universität Braunschweig, 38106 Braunschweig,
Germany<EMAIL_ADDRESS>
##  1 Introduction

This document is a short and informal tutorial on some aspects of calculating
phase diagrams with the ATAT-tools emc2 and phb and on creating cluster
expansions with maps. It is neither complete, nor in any way an official
document, but mainly a set of collected notes I took during experimentation
with ATAT. I assume that you have read at least [1] and have some basic idea
about ATAT and what it does (but may be stymied by the large number of options
and program variants). For the final sections on phonon and electron
contributions, reading [2] is also helpful. Note that all programs in the
ATAT-toolkit have a -h option which gives a detailed explanation of (most)
parameters you can set and of input and output files.

##  2 Chemical potential

The main variable to govern the composition of the alloy is the chemical
potential. Since this is a semi-grand canonical ensemble (SGCE), the number of
atoms is fixed (if vacancies are present, they are counted as an atom
species); each lattice site contains one atom.

A change of concentration thus means replacing atoms of one sort with one of
the other. (For simplicity, we are looking at two species only.) If we change
one atom, we have (here and in the following, I use $F$ for the free energy
and do not consider $G$) [3, 8.3]

$\frac{\mathrm{d}F}{\mathrm{d}n}=\mu_{A}-\mu_{B}=\mu$ (1)

since we change one atom from A to B. If the chemical potential is defined wrt
the concentrations, replace $n$ by the concentration $x$. Note that in ATAT
(at least in the binary version), the chemical potential used is the
difference between the two species, so there is only one value.

In the semi-grand canonical ensemble, the chemical potential is an external
control variable (like the pressure in the NPT ensemble). So it is defined by
the change in the free energy but it is nevertheless prescribed externally.
The chemical potential can be understood as the thermodynamic force trying to
“push” an atomic species into the system. A chemical potential of $0$ means
that there is no external driving force to replace one species with the other.

In ATAT, there are two conventions for the chemical potential: In the output,
it is given as a physical quantity in the units used (depending, among others,
on the value of $k_{B}$ given with -k); in the input of emc2, the chemical
potential is normalized so that the region where the first phase (from
gs_str.out) is stable corresponds to $\mu\in[0,1]$, the second phase is stable
in the region $\mu\in[1,2]$ and so on. (The disordered phase has negative
$\mu$). Note that the normalization is not done in this way if there are only
two or one ground states (see also section 5).

Therefore, if you call emc2 with
emc2 -gs=1 -mu0=1.5 -mu1=0.5 -dmu=0.04 …
you start at a value of $1.5$ that perfectly stabilizes the second phase
(ground state number 1) and do calculations down to a value of $\mu$ where the
first phase (ground state number 0) is stable. If the specified $\mu$-value
does not correspond to the specified ground state, you may get a warning that
the chemical potential does not stabilize this ground state.

##  3 Phase diagrams

Phase diagrams can be drawn with the chemical potential instead of the
concentration as variable. In a two-phase region, because of the standard
common-tangent construction, the chemical potential is constant. Therefore,
the two-phase region collapses to a line in a $\mu$-$T$ phase diagram. The
advantage of using $\mu$ (and thus the SGCE) instead of $x$ is that only one
phase is stable at any value of $\mu$ (unless we are _exactly_ at a
$\mu$-value that corresponds to a 2-phase region). Therefore, in a Monte-Carlo
simulation, the system should never separate into different phases [4, p. 24].

To calculate phase boundaries, the “thermodynamic function” [1, 4]

$\displaystyle\phi(\beta,\mu)$ $\displaystyle=E-TS-\mu x=-\frac{1}{\beta
N}\ln\sum_{i}\left(\exp(-\beta N(E_{i}-\mu x_{i}))\right)$ (2)
$\displaystyle=F-\mu x$ (3)

is used.

This can be calculated for each phase – the thermodynamically stable phase
will have the lower value of $\phi$; on a phase boundary, the values for both
phases are equal. Therefore, phase boundaries between phase $\alpha$ and
$\gamma$ are defined by

$\Lambda^{\alpha,\gamma}=\\{(\beta,\mu):\phi^{\alpha}(\beta,\mu)=\phi^{\gamma}(\beta,\mu)\\}\,.$
(4)

Since the chemical potential is conjugate to the concentration in a SGCE, the
concentration $x$ in each phase can be calculated from

$x^{\alpha}=-\frac{\partial\phi^{\alpha}(\beta,\mu)}{\partial\mu}\,.$ (5)

Since cluster expansions have their origin in an analogy to the Ising model,
the concentration runs from $-1$ (pure A) to $+1$ (pure B). To get physical
concentrations, use $(x+1)/2$.

To calculate the phase boundary, $\phi$ has to be calculated, but with MC
simulations, only differences (total differentials) can be calculated. The
method can still be employed by starting at a point where $\phi$ is known, for
example from a low-temperature expansion. This is what the program phb does
when the option -ltep is used.

The program emc2 also calculates $\phi$; there, an LTE or HTE is used if no
initial value of $\phi$ (option -phi0) is given.

##  4 Tracking phase boundaries

The basic strategy to follow the phase boundary is from [1, eq. (29)] or [5]:

$\frac{\mathrm{d}\mu}{\mathrm{d}\beta}=\frac{E^{\gamma}-E^{\alpha}}{\beta(x^{\gamma}-x^{\alpha})}-\frac{\mu}{\beta}\,.$
(6)

So we can calculate the change of $\mu$ with temperature from a given point.
So the calculation proceeds by starting at a known $(T,\mu)$-point and goes on
from there, incrementing $\mu$ in finite steps of $\beta$ (or the
temperature).

The parameter -dT in phb thus affects the precision of following the boundary:
Making it smaller means that it is easier to follow a curved boundary.

To check how this works, we can use the example file provided in the mc-folder
of ATAT:

    T          mu        x1         x2         E1           E2
    240  -0.0775028  -0.986175  -0.502882  -0.0495888  -0.0501994
    250  -0.0774447  -0.982053  -0.503663  -0.0493333  -0.0501169¯

The corresponding beta-values are $1./(240\cdot\numprint{8.617e-5})=48.35403$
and $1./(250\cdot\numprint{8.617e-5})=46.41987$. Thus, $\delta\beta$ is
$48.35403-46.41987=1.93416$. With these values, the left-hand side of the
above equation is

$(-0.0775028+0.0774447)/1.93416=\numprint{-3.0039e-5}\,.$

In the current version of ATAT (5/10/18), the last two columns are actually
$E-\mu x$, so the second term of eq. 6 is already included in the “energy”.
The rhs thus becomes

$((-0.0495888+0.0501994)/(-0.986175+0.502882))/48.35402=\numprint{-2.6128E-5}$

which is close to the lhs (although not exactly identical).

A look into the source code of phb.c++ shows in line 416:

          mu+=1.5*dmu-0.5*old_dmu;

so the discretization is not as simple as I assumed here since it uses the
previous value of dmu as well. Nevertheless, the results are close enough to
see that the calculation of $\mu$ does in principle proceed as explained here.

##  5 The simplest possible example

We now try to perform a simple phase diagram calculation from scratch. Use the
following settings to create a binary system where the two species sit on a
simple cubic lattice and separate:111The atoms are called Al and Ni but have
nothing to do with the materials, this is just an example.

    lat.in
    3.5 3.5 3.5 90 90 90
    1. 0 0
    0 1 0
    0 0 1
    0 0 0 Ni, Al

    gs_str.out
    3.500000 0.000000 0.000000
    0.000000 3.500000 0.000000
    0.000000 0.000000 3.500000
    1. 0 0
    0 1. 0
    0 0 1.
    1.000000 1.000000 1.000000 Ni
    end

    3.500000 0.000000 0.000000
    0.000000 3.500000 0.000000
    0.000000 0.000000 3.500000
    1. 0 0
    0 1. 0
    0 0 1.
    1.000000 1.000000 1.000000 Al
    end

These two set up a simple cubic lattice (in lat.in) and the two ground-state
phases (pure “Ni” and “Al” which of course would not form a simple cubic
lattice). Note that the actual value of the lattice constant is totally
irrelevant here since ATAT only looks at cluster configurations.

The clusters are defined in

    clusters.out:
    1
    0.000000
    0

    1
    0.000000
    1
    1.000000 1.000000 1.000000

    6
    3.5
    2
    1.000000 1.000000 1.000000
    1.000 1.0000 0.0000

The first term is a constant term (with no atoms considered), the second is
the single atom cluster (its expansion coefficient gives the energy difference
between a single Ni and a single Al atom), the third is the interaction of
nearest neighbours.

Here are the coefficients to be used with the clusters:

    eci.out
    0.
    0.
    -1

The only one that is really important is the third. It is $-1$, so when
considering two identical atoms (species both $+1$ or $-1$), the energy
contribution is negative, when considering non-identical atoms, the energy is
positive. Therefore, we should expect the system to separate into two phases.

A calculation using phb shows that the chemical potential (as output variable,
see above) is zero at the phase boundary at 0 K, which should be expected
because of symmetry. The energy per atom is $-3$ (each atom has 6 nearest
neighbours, but to avoid double-counting of pair bonds, the energy per atom is
$3\cdot(-1)$).

To calculate the phase diagram, use (note that linebreaks in command lines
should of course not be typed in):

    phb -gs1=0 -gs2=1 -dT=2000 -dx=1e-2 -er=20 -k=8.617e-5
    -ltep=5e-3 -o=ph01.out

The resulting phase diagram, fig. 1, left, shows the miscibility gap. You can
plot it with gnuplot [6] using

    plot [0:1][0:55000] "ph01.out" using (($3+1)/2):1 w lp,\
    "ph01.out" using (($4+1)/2):1 w lp

Figure 1: Phase diagram for the simple system of non-mixing species. Left:
initial run with open gap. Right: Trying to close the gap by changing dT or
er.

The calculation of phb proceeds rapidly up to a temperature of 50000, but then
it stalls. The output from phb shows that a huge number of MC iterations would
be required:

    Phase 1 n_equil= 0 n_avg= 1280000

This is an example of critical slowing down, i. e. the difficulty of finding a
new configuration that is statistically independent from the current one. So
unless you are willing to wait for a very long time, there will be a gap in
your phase diagram at the top of the msicibility gap.

You can try to reduce the temperature-steps to make the gap smaller:

    phb -gs1=0 -gs2=1 -dT=200 -dx=1e-2 -er=20 -k=8.617e-5
     -ltep=5e-3 -o=ph01-200.out

The result is shown in fig. 1, right. Still, at a temperature of 51000, the
number of equilibration steps becomes rather large (72000). You can of course
wait for this to finish, it won’t take too long.

One other way to deal with this that actually helps is to increase the radius
er:

    phb -gs1=0 -gs2=1 -T=46000 -mu=2.9118e-5 -dT=500 -dx=1e-2 -er=50
     -k=8.617e-5 -ltep=5.e-3 -o=ph01-er50.out

Here we start at a reasonable point for $T$ and $\mu$ (taken from the previous
run, of course the phase boundary should be exactly at $\mu=0$ in this
example) but use a larger cell. This run reaches a temperature of 52000 faster
than the other two. (After that, it also stalls) This might seem surprising
because usually critical slowing down becomes worse when the simulation volume
gets larger. I suspect (but do not know) that this is because ATAT calculates
the number of runs needed and there is a trade-off between the needed number
of runs and the simulation volume. So it seems that smaller cells will not
always run faster. The resulting phase diagram is almost closed, but not
quite. The figure also shows that the run with the larger cell is smoother as
should be expected. (And a radius er of 20 is really a bit small for a MC
calculation.)

If this does not work, here is some advice from Axel van de Walle [private
comm.]:

> The way I usually handle this is by running emc2 around the expected top of
> the miscibility gap to explore where the gap may start. emc2 also has the
> option of running for a fixed number of steps (-n option), which bounds the
> stalling time. When you have the top point, you can usually just join it
> with the curves from phb and get a decent-looking phase diagram.

The idea behind this is that emc2 runs at fixed value of $\mu$ so that the
system should always be in a single phase. (emc2 actually contains some checks
on this, see below.) On varying $\mu$, you will then “jump” in the
concentration from one phase to the other. So, following this advice, run emc2
as follows:222emc2 writes several files: ltedat.out, htedat.out, mfdat.out,
and mcsnapshot.out. If you run several instances of emc2 in the same
directory, be aware that these files will be overwritten. The emc2-run itself
is unaffected by this, but if you want to look at configurations
(mcsnapshot.out), it might be better to run each emc2-instance in a separate
directory.

    emc2 -gs=0 -mu0=-0.5 -mu1=0. -dmu=0.1 -T0=2000 -T1=80000 -dT=2000 -keV
       -er=20 -dx=1.e-3 -o=emc-01.out > emc-01-run &
    emc2 -gs=1 -mu0=0.5 -mu1=0. -dmu=0.1 -T0=2000 -T1=80000 -dT=2000 -keV
       -er=20 -dx=1.e-3 -o=emc-10.out > emc-10-run &
    emc2 -gs=0 -mu0=-0.5 -mu1=0. -dmu=0.1 -T0=2000 -T1=80000 -dT=2000 -keV
      -er=40 -dx=1.e-2 -o=emc-01-40.out > emc-01-40-run &

Note that in this case, where we only have two ground states, the
$\mu$-parameter works differently. The help-page of emc2 states:

> If there are only two ground states, the only correction performed is to
> shift mu so that mu=0 stabilizes a two-phase equilibrium between the two
> ground states.

So in our case, the phase transition is at $\mu=0$ (no shift is necessary) and
not at $\mu=1$ (as you might expect, since in the general case, $\mu=1$
stabilizes the boundary between phase 0 and 1). I did one run with a larger
cell size (and larger dx) to check that my cell is not too small. I also piped
the output to a file and sent the process to the background.

We can now plot the $x$-$T$ plane of the states looked at by emc2 and compare
them to the phb-diagram, see fig. 2, left. This plot is created using gnuplot:

    set key bottom
    plot "ph01.out" using (($3+1)/2):1 title "phb run" w lp ls 1,\
     "ph01.out" using (($4+1)/2):1 notitle w lp ls 1 ,\
     "emc-01.out"  using (($4+1)/2):1 title "emc mu0-1" ls 2,\
     "emc-10.out"  using (($4+1)/2):1 title "emc mu1-0" ls 3,\
     "emc-01-40.out"  using (($4+1)/2):1 title "emc mu0-1 40" ls 4

I use columns 3 and 4 of the phb-output (the two concentrations at the phase
boundary) and column 4 of the emc2-output (concentration at the point
currently looked at). Concentration values are converted from $[-1,1]$ to
$[0,1]$ as explained above. The output looks reasonable, although the
boundaries do not agree perfectly (deviation seems to be larger than the
required precision at higher temperatures). This may be due to a slight
overshooting of emc2 due to the hysteresis loop.

Figure 2: Left: Running emc2 on the simple system. The points looked at by
emc2 should all be in the one-phase region, so the boundary of these points
gives a hint of the phase boundary. Right: Adding a run with tstat=0.

You can also plot the potential from emc2, for example

    splot  "emc-01.out" using (($4+1)/2):1:5,\
      "emc-12.out" using (($4+1)/2):1:5

This plots the potential $F-\mu x$ versus the concentration and the
temperature. You expect to see a surface with gaps (in the concentration
region where two phases are in equilibrium) as shown in fig. 3. To see the
phase transition, we can plot $\phi$ vs. $\mu$ and $T$ as shown on the right.
The two surface cross at $\mu=0$ as expected.

Figure 3: Left: potential $\phi$ vs. concentration and temperature; plotted
with gnuplot Right: Plot of $\phi$ vs. $\mu$ and $T$ to show the phase
transition at $\mu=0$.

The states in fig. 2 look reasonable, but they do not follow the phase
boundary very closely. You can estimate that the top of the miscibility gap is
somewhere between 52000 and 60000 which is perhaps not precise enough.

One possibility to study this further is to switch off the checking of phase
transitions in emc2 that causes emc2 to stop at some point. This is done with
the option tstat=0. Since this stalls after some time, I prescribe the number
of equilibrations and calculation steps. This is of course a bit dangerous
since you cannot rely on the statistics to be good enough, so you have to
experiment to see whether things work out (for publishable results, you would
surely use a larger cell and more sweeps):333If you are patient enough, you
can also specify -dx; at least for -gs=0, this finished after about a day of
computing.

    emc2 -gs=0 -mu0=-0.1 -mu1=0.1 -dmu=0.02 -T0=40000 -T1=80000 -dT=2000
      -keV -er=20 -n=5000 -eq=2000 -tstat=0 -o=emc-01-tstat.out
    emc2 -gs=1 -mu0=0.1 -mu1=-0.1 -dmu=0.02 -T0=40000 -T1=80000 -dT=2000
      -keV -er=20 -n=5000 -eq=2000 -tstat=0 -o=emc-10-tstat.out

The result is shown in fig. 2, right. Note that – as should be expected – some
of the data points lie in the miscibility gap, but there are some points
almost exactly at $x=0$ that look like a good candidate for the center of the
miscibility gap. Look into the corresponding output file from emc2 to find
their $\mu$\- and $x$-values:

    emc-01-tstat.out:
    52000.000000¯0.000000¯-1.111069¯0.131531¯
    54000.000000¯0.000000¯-0.899093¯0.004387¯
    56000.000000¯0.000000¯-0.812865¯-0.014910¯
    emc-10-tstat.out:
    52000.000000¯-0.000000¯-1.108399¯0.021828¯
    54000.000000¯-0.000000¯-0.905569¯-0.006905¯
    56000.000000¯-0.000000¯-0.817109¯0.001020

The first concentration value (column 4) is a bit off, all others are close to
a concentration of almost exactly 0. All their $\mu$-values are 0. (This is
the actual, physical $\mu$-value, not the value in the ATAT-input convention.
That it is exactly 0 is due to the symmetry in this example.)

So we can expect a plausible endpoint of the miscibility gap to be at $\mu=0$
somewhere around a temperature of 52000. To check this, we can run phb from
this point _downwards_ to see whether we meet our old phase separation line.
To do so, I increased the size of the cell to 50:

    phb -gs1=0 -gs2=1 -T=52000 -mu=0 -dn -dT=500 -dx=1e-2 -er=50
        -k=8.617e-5 -ltep=5.e-3 -o=ph01-down1.out
    phb -gs1=0 -gs2=1 -T=52500 -mu=0 -dn -dT=500 -dx=1e-2 -er=50
        -k=8.617e-5 -ltep=5.e-3 -o=ph01-down2.out
    phb -gs1=0 -gs2=1 -T=53000 -mu=0 -dn -dT=500 -dx=1e-2 -er=50
        -k=8.617e-5 -ltep=5.e-3 -o=ph01-down3.out

The result of these runs is shown in fig. 4. There is still some noise at the
highest point, but all three runs actually meet with our old phase separation
line. (In a downwards run, phb is only as good as its starting point, so it is
necessary to make sure that you reach a reasonable point.) Thus, things look
quite reasonable.

We can thus use a point on this line to start upwards again (we might also
have done this directly from our first phb-run, but I chose this more
complicated way to illustrate the possibilities), this time with an even
larger cell and a small step size for the temperature:

    phb -gs1=0 -gs2=1 -mu=-2.94e-5 -T=51000  -dT=-100 -dx=1e-2 -er=80
     -k=8.617e-5 -ltep=5.e-3 -o=ph01-upb.out

With this, we get some fluctuations at the top of the miscibility gap, see
fig. 4, right, but they are rather small. Taking all figures together, we can
estimate the temperature of the gap to about $52500\pm 1000$, so the accuracy
is about 2%.

Figure 4: Left: Downwards run using phb from the center point from the
emc2-output with three different starting temperatures Right: Using a starting
point from the downwards runs to run upwards again.

Finally, we can check whether results make sense at all: According to [7], the
miscibility gap should be at

$T_{\text{misc}}=\frac{0.8\Omega_{AB}}{2k_{B}}\text{ with
}\Omega_{AB}=-(z/2)(E_{AA}+E_{BB}-2E_{AB})\,,$ (7)

where $z$ is the coordination number. (According to [8, ex. 5.10], for a
binary system in zeroth order, the temperature is $\Omega/2k_{B}$ regardless
of lattice structure.) For my data ($z=6$, $\Omega_{AB}=12$,
$k=\numprint{8.617e-5}$ in units of eV), the result of this is a temperature
of about 54000, so this looks good enough.

##  6 A more interesting example

A more interesting example is one where a different phase forms. Change the
ECIs, the ground state file and the cluster definition as follows:

    gs_str.out
    3.500000 0.000000 0.000000
    0.000000 3.500000 0.000000
    0.000000 0.000000 3.500000
    1. 0 0
    0 1. 0
    0 0 1.
    1.000000 1.000000 1.000000 Ni
    end

    3.5000 3.5000 0
    3.5000 0 3.5000
    0 3.5000 3.5000
    1. 0. 0.
    0. 1. 0.
    0. 0. 1
    0. 0. 0. Al
    0.5 0.5 0.5 Ni
    end

    3.500000 0.000000 0.000000
    0.000000 3.500000 0.000000
    0.000000 0.000000 3.500000
    1. 0 0
    0 1. 0
    0 0 1.
    1.000000 1.000000 1.000000 Al
    end

The state 0 and 2 are simply the pure elements on a simple cubic lattice. The
second is the primitive unit cell of a checkerboard pattern (NaCl structure).
To visualize it, you can copy this definition into a separate file. Since the
ATAT out-format is a bit special, I use Jesper Kristensen’s tools [9] to
convert it into a POSCAR-format that I can open with the ase tools [10]. To
visualise the clusters, you can use Jesper’s script visualize_clusters.py
[11].

To stabilize the checkerboard pattern, we need to define more clusters than
before:

    eci.out
    0.
    0.
    1
    -0.2

    clusters.out
    1
    0.000000
    0

    1
    0.000000
    1
    1.000000 1.000000 1.000000

    3
    3.5
    2
    1.000000 1.000000 1.000000
    1.000 1.0000 0.0000

    3
    7.00000
    2
    1.00000 1.00000 1.00000
    1.00000 1.00000 -1.00000

Here we have changed the sign of the 2-atom cluster term so that AB as nearest
neighbours is energetically favourable. In order to force a third AB-phase to
form, we need to make an ABA- or BAB- structure favourable, otherwise the
atoms would just dissolve freely and the alloy has full miscibility. The final
cluster in cluster.out contains two atoms on the same “color” of the NaCl-
checkerboard pattern and the negative sign of its ECI makes a sequence AXA or
BXB favourable. So (writing in one dimension) a pattern of type ABABABAAAAA is
better than, for example, ABAABAAAABA because we have more identical second
neighbours if we create a separate AB-phase.

To create these structures, you can use the corrdump-utility on the lat.in-
file; this creates a cluster-file that can be used to define the clusters. To
find all 2-atom clusters that are smaller than a certain size limit (up to
size 7 in this example), simply copy the lat.in-file to a clean directory and
run

    corrdump -2=7.01

This creates a clusters.out-file from which you can pick the clusters you
want.

Of course, usually you do not construct clusters and their ECIs like this (as
we will see below, this method of guessing ECIs is dangerous), but you rather
get them from another calculation (like from ATAT’s maps program, as we will
do later). Nevertheless, this example helps to improve understanding how phb
works and what problems might occur.

To check that the structure definition is correct, we can create
configurations using emc2:

    emc2 -gs=0 -mu0=0.5 -T0=30000  -keV -er=20 -dx=1-e-3
    emc2 -gs=1 -mu0=1.5 -T0=30000  -keV -er=20 -dx=1-e-3
    emc2 -gs=2 -mu0=2.5 -T0=30000  -keV -er=20 -dx=1-e-3

Each run creates a configuration file mcsnapshot.out which we can again plot
using Jesper Kristensen’s tool [9] to convert out-files. Fig. 5 shows the
three structures.

Figure 5: Structures of the three phases.

You can also use emc2 to understand the numbers for the chemical potential.
Simply run emc2 without any MC sweeps whatsoever:

    emc2 -gs=0 -mu0=0.5 -T0=10  -keV -er=20 -n=0 -eq=0 -g2c
    Supercell size: 12 12 12
    10.000000¯-12.000000¯2.400000¯-1.000000 ......¯

and similarly for the other phases. The -g2c-switch tells emc2 to report
quantities for the canonical ensemble; column 3 is thus the energy per atom.
This is $2.4$ (in our units): Each atom has six nearest-neighbour bonds with
energy $+1$ and six next-nearest neighbour bonds with energy $-0.2$, resulting
in an energy of $4.8$ which needs to be divided by 2 to avoid double-counting
of the bonds.

The result of $-12$ for the chemical potential may seem weird. If we are in
the A-rich phase and exchange a single A-atom with B or a B-atom with A, the
total energy change is not $\pm 12$, but only $\pm 9.6$. (6 bonds flip from
$+1$ to $-1$, 6 flip from $-0.2$ to $+0.2$.) ATAT calculates the chemical
potential from the energy difference of the provided ground states. If all
ground states are given correctly, the common-tangent construction then
implies that the chemical potential (change in energy on switching an atom in
the very-low temperature limit) can be calculated this way. (This is also
shown in Fig. 5.3 of the ATAT manual.) As we will see later (spoiler alert),
our list of ground states is incomplete, and so we have a seeming discrepancy
here.

The clusters were designed to create a phase diagram with a 3-phase structure
(A-rich, B-rich, AB “intermetallic”) with a miscibility gap. It should be
absolutely symmetric. Since we have two phase boundaries, we need to do 2 phb
runs, one for the phase boundary between phase 0 and 1 and one for the
boundary between 1 and 2\. Actually, I do a third run just switching 1 and 2
in the gs-parameter to see what happens:

    phb -gs1=0 -gs2=1 -dT=100 -dx=0.001 -er=50 -keV -ltep=1.e-3 -o=phase01.out
    phb -gs1=1 -gs2=2 -dT=100 -dx=0.001 -er=50 -keV -ltep=1.e-3 -o=phase12.out
    phb -gs1=2 -gs2=1 -dT=100 -dx=0.001 -er=50 -keV -ltep=1.e-3 -o=phase21.out

In the following, runs will be abbreviated as “01”, “12” and “21”.

Using these numbers, I would expect the 01-run and the 21-run to be
practically the same. If you compare the numbers in the output files
phase01.out and phase21.out, you see that they are identical for temperatures
up to 2700 where the LTE is valid; the only difference between the two files
is that mu is $-6$ in one case, $+6$ in the other.

At higher temperatures, the runs start to differ. Since these are Monte-Carlo
algorithms using random numbers initialized from the system clock unless
otherwise specified, two runs with the same parameters will never be
identical, but I would still expect the runs to be very similar. This,
however, does not happen. The 01-run proceeds up to a temperature of 63400;
the largest number of MC iterations is 32000. (The exact number may vary
depending on the random numbers in your MC sequence.) The 21-run (which should
be identical) becomes very slow, at a temperature of 13200, we get

    Phase 1 n_equil= 64000 n_avg= 896000
    Phase 2 n_equil= 0 n_avg= 2000

The 12-run shows another problematic behaviour: At a temperature of 13200, we
get the following output:

    13200¯5.96865¯0.450491¯0.746105¯-3.34662¯-3.35492¯
    Phase 1 n_equil= 0 n_avg= 4000
    Phase 2 n_equil= 0 n_avg= 3000
    0¯1
    Looking for phase transition...
    mu           x
    5.992864¯0.768726
    6.016864¯0.782460
    ...
    45.976864¯1.000000
    46.000864¯1.000000

The program stalls and never proceeds to higher temperatures. This may happen
when program “wanders outside of the region of metastability” [1, p. 15] and
tries to find the correct mu-value.

The resulting phase diagram is shown in fig. 6, left.

Figure 6: Left: Calculation of the phase diagram with phb for the 3-phase
alloy. Note that the 12- and 21-runs are stalling whereas the 01-run
calculates everything up to the miscibility gap. Right: Zoom of the
calculation re-done with smaller temperature steps and with larger radius.

Only one of the runs calculates the full 2-phase region up to the high-
temperature state. It is also interesting that there are strange “spikes” at
concentrations of about 0.25 and 0.75 that are possibly due to phb not
correctly tracking the phase boundary. (More on these below.) At the same
temperature, there are smaller discontinuities at concentrations of about 0.12
and 0.88.

To study this further and get a better result in the right-hand part of the
diagram, we can re-run phb with smaller step size for dT. To do so, use values
for $T$ and $\mu$ from the out-file of the 12-run (at temp 12000, the value
for the 21-run differs only in the last digit)

    12000  5.97369   0.380374   0.790879 ....

and start phb from this point (I could of course start from the LTE, but then
the calculation would take much longer, even so, the following runs need
several hours):

    phb -gs1=1 -gs2=2 -T=12000 -mu=5.97369  -dT=10 -dx=0.001 -er=50 -keV
     -ltep=1.e-3 -o=phase12-dt10.out
    phb -gs1=2 -gs2=1 -T=12000 -mu=5.97369 -dT=10 -dx=0.001 -er=50 -keV
     -ltep=1.e-3 -o=phase21-dt10.out

I also run it again with a larger cell (er):

    phb -gs1=0 -gs2=1 -dT=2000 -dx=1e-2 -er=80 -k=8.617e-5 -ltep=5.e-3
    -o=ph01-er80.out

The results are shown in 6, right. The smaller temperature step reproduces the
behaviour from the other side of the concentration region, including the
strange spikes in the curve. They are also a bit strange since the line is not
smooth but seems to oscillate. The run with larger radius seems to close the
phase diagram, but although the two red lines cross, phb does not finish the
calculation.

So things look really strange here. To find out what is going wrong, we can
look at the most critical aspect of the behaviour of the simulation: The spike
at concentrations 0.75/0.25. Perhaps something interesting is happening at
these concentrations?

To find out, we can use emc2. Since ATAT always gets chemical potentials as
input variables, not concentrations, it is difficult to fix a certain
concentration value.444You can specifiy a concentration with -x in emc2, but
this only affects the initialization. We could of course change the source
code (it would probably not be too difficult to use the existing routines to
write a program that does a canonical MC simulation), but we can also simply
try to guess a good value for the chemical potential. Since $\mu=2.5$ (input
variable $\mu$) perfectly stabilizes pure B (concentration 1.0) and $\mu=1.5$
stabilizes the AB-phase at concentration 0.5, a value around 2.0 should be
fine. We run emc2 with a small cell at several $\mu$-values in this region at
a very low temperature. We use a large number of equilibration steps and only
one actual measurement step so that the output file describes the final state.
We start at a disordered state and ignore the phase structure (tstat=0) so
that emc2 can try to simply find the lowest-energy structure. So run several
commands of this type:

    emc2 -gs=-1 -mu0=1.96  -T0=1000  -keV -er=10 -n=1 -eq=500000 -tstat=0

with different values of -mu0 and check the output file. And indeed, at a
value of 1.96, we find a concentration (in ATAT units) of 0.5 for the final
state (physical concentration of 0.75). At -mu0 of 1.97, we find another
structure which has a concentration of 0.75 (0.875 in physical units). Note
that these runs are sensitive to random number initialization, so you may find
slightly different values for -mu0 to stabilize the structures we are looking
at. Depending on the numbers you might also find other stable
structures.555This also explains why we got some seemingly inconsistent result
for the calculation of the chemical potential at the beginning of this
section.

To understand these structures, we can plot the end states using the
mcsnapshot.out files, see fig. 7. Looking at these structures, we see that
each A-atom is surrounded only by B-atoms as nearest neighbours (the optimal
structure for the nearest-neighbour cluster) and that each atom also has
another atom of the same species as next-nearest neighbour along each axis. So
these structures are also stable with our cluster definitions.

Figure 7: Stable structures with concentrations of 0.75 and 0.875.

Thus, the problem in this example lies not with phb, but with my incorrect
definition of stable structures; there are stable structures that were not
included in my gs_str.out-file. The weird behaviour of phb was thus due to
structures interfering that were not included in the definition.

##  7 Calculating a phase diagram from first principles

###  7.1 The Ni-Al phase diagram: First attempt

It seems that guessing cluster coefficients is not so easy – so let’s do a
density functional theory calculation to get a meaningful cluster expansion.
This is done using maps. Since I’m working a lot with nickel alloys, we can
look at the NiAl-system.

The basic idea is explained nicely in the ATAT manual: maps is the control
program that creates new configurations to calculate, pollmach is the program
that manages the calls to the first principles code (VASP in my case) and this
calls runstruct_vasp. runstruct_vasp calls another program called ezvasp that
performs the actual runs.

Before you start, you need a ~/.ezvasp.rc-file in your home directory that
tells ezvasp where the potentials are and how VASP is to be called. My file
looks like this:

    #!/bin/csh
    #enter name of vasp executable here
    set VASPCMD="mpirun -n 12 vasp"
    #enter the directories containing the pseudopotentials here
    set POTPAWLDA="/opt/vasp/potentials/potpaw_lda/"
    set POTPAWGGA="/opt/vasp/potentials/potpaw_gga/"
    set POTPAWPBE="/opt/vasp/potentials/potpaw_pbe.5.4/"
    set POTLDA=$POTPAWLDA
    set POTGGA=$POTPAWGGA

As you can see, I want to call vasp with mpirun and run on 12 of my 16 cores.
I do run in serial mode since I am doing all this on a simple workstation, not
on a massively parallel machine. The VASP potentials are situated in the /opt-
directory. Note that ezvasp -h gives you lots of informations on how to do the
VASP runs, but note also that you do not directly interface ezvasp but only do
it with runstruct_vasp. For example, ezvasp tells you that you can write
magnetic moments in the POSCAR section of the file ezvasp gets, but as far as
I can see you cannot do this if you call it from runstruct_vasp via pollmach.

To do the calculation, you need a lat.in-file as before and a wrapper file.
Here is my vasp.wrap:

    [INCAR]
    PREC = high
    ISMEAR = 1
    SIGMA = 0.07
    NSW=41
    IBRION = 2
    ISIF = 3
    KPPRA = 1000
    USEPOT = PAWPBE
    NPAR=4
    ISPIN=2
    SUBATOM = s/Ni$/Ni_pv/g
    EDIFF = 1e-8
    EDIFFG = 1e-5
    DOSTATIC

I have taken most of the parameters from the ATAT example. Note that I use the
SUBATOM-command that is documented in the ezvasp-help. This allows to replace
an element name with the name of the potential – for Ni, it is good to include
the $p$-electrons and use the Ni_pv-potential. Note also that I use tight
values for EDIFF and EDIFFG. The smaller value for SIGMA is mainly for
consistency with previous calculations I have made with nickel.

Actually, I should also have entered a sensible ENCUT, otherwise, the pure Al
run in directory 1/ has a different cutoff from the other simulations.
Checking afterwards showed that the energy of the Al run was only affected
slightly by this. Still, I strongly recommend to always set ENCUT (which is a
good practice for VASP in general).

Nickel is ferromagnetic, so I also include ISPIN=2. You cannot specify a
MAGMOM in the wrapper-file because this gets directly copied to the later
INCAR-file and as soon as the first 2-particle calculation is done, VASP
complains about the MAGMOM-line. However, for nickel, this should not matter
much because the default initialization of 1 for MAGMOM is fine. (If you do
not include the spin polarization, ATAT finds a stable Ni7Al-phase as
described briefly in [12].)

For completeness, here is the simple lat.in with a reasonable value for the
lattice constant:

    3.52 3.52 3.52 90 90 90
    0   0.5 0.5
    0.5 0   0.5
    0.5 0.5 0
    0.000000 0.000000 0.000000 Ni,Al

To start everything, do (after making sure that your VASP settings are good
enough to give precise results)

    maps -d &
    touch ready
    pollmach runstruct_vasp &

and then maps runs and creates configurations. The maps output is again
explained well in the manual. Note that you can find the ground states of the
system not only in gs_str.out in the same format as before, but also in gs.out
where the structure number (number of the subdirectory) is given. (So if you
want to visualise the ground states, you can use the POSCAR-files in the
directories listed in gs.out.) The file also lists the energy the structure
actually has and the energy you get from the cluster expansion.

If you get errors like (in your system language)

    cp: cannot stat ’OSZICAR’: No such file or directory
    cp: cannot stat ’OUTCAR’: No such file or directory

something is wrong with the VASP call; runstruct_vasp does not directly check
whether VASP has run so the code only notices this when it tries to copy some
files. Check your VASP installation or the interface between maps and VASP.

During the run, maps reports on the performed calculations; don’t worry if you
see lines like

    Finding best cluster expansion...
    1 1 1 1 3.40282e+38

the FLOATMAX is no problem.

After about two days of running (on my machine, using 12 Xeon processors),
ATAT finds four phases (beyond the pure Ni and Al phase): ${\rm Ni_{3}Al}$,
${\rm Ni_{5}Al_{3}}$, ${\rm NiAl}$, ${\rm Al_{4}Ni}$. Comparing this to the
Ni-Al phase diagram, it seems that we get the nickel-side correct, the Al-side
should have an ${\rm Ni_{2}Al_{3}}$\- and an ${\rm Al_{3}Ni}$-phase. A look at
the energies with mapsrep (figure 8) shows that ${\rm Al_{4}Ni}$ is only
barely stable. (The data point lies almost exactly on the common tangent
between Al and NiAl.) We stop the simulation with touch stoppoll. Wait until
the current VASP run is finished and maps has updated the results, then touch
stop.

Figure 8: ATAT calculation result plotted with mapsrep for the NiAl system.

To analyze the results, we can also use the program checkrelax. This
calculates how distorted the relaxed cells of the different runs actually are
(it measures the strain in the cells). Do

    checkrelax > checkrelaxresult

and look at the output file checkrelaxresult. In our case, some of the cells
are very strongly distorted, the file ends with:

    0.2024 54
    0.2069 3
    0.2115 46
    0.2186 8
    0.2456 48
    0.2535 25

If strains are very large, the result is not really useful for the cluster
expansion, because in doing the MC runs later on, we assume that things happen
on our prescribed lattice. Expansion coefficients calculated from a strongly
distorted cell may thus be incorrect.

In principle, we can exclude all those configurations. (See section 7.3 below
on how to do that.) However, in our case there is a problem: Configuration 3
is a ground state as you can see by looking into gs.out which lists all the
ground states. If you look into the corresponding POSCAR in directory 3, you
see that this is the ${\rm NiAl}$-phase, but has an orthorhombic cell instead
of the correct cubic CsCl-structure. The reason for this problem is simple:
The CsCl-structure is a bcc-structure, not an fcc-structure, and thus cannot
be found with the current setting of lat.in. ATAT creates an initial fcc-
lattice that is occupied by alternating Ni and Al atoms and this then gets
distorted to a bcc lattice. (This is reminiscent of the Bain model of the
martensitic transformation in iron [13].)

In any case, if you do this exclusion and calculate a phase diagram (I tried),
you will find that your phase structure looks simply awful and that the
calculation with phb frequently stalls. If you actually need the full Ni-Al
phase diagram, you could do a maps calculation with a bcc lattice as well; you
would then have to do the phase diagram calculations with several different
phases.

However, for this tutorial we choose a simpler way.

###  7.2 Calculating the Ni-rich part of the phase diagram

In maps, you can specify the concentration range you are interested in via the
command-line options -c0=a -c1=b. Let us restrict ourselves to the calculation
of the phase boundary between Ni and ${\rm Ni_{3}Al}$ to avoid trouble with
the bcc phase. Create a new directory, set up lat.in and vasp.wrap as before
and then call maps as

    maps -d -c0=0 -c1=0.25 &
    touch ready
    pollmach runstruct_vasp &

(Instead of doing all the calculations again, it would probably be easier to
copy the old result directories for those structures that are in the right
part of the phase diagram, but I prefer to do a clean start.) After a weekend
of calculations, we can stop the calculation with touch stoppoll and touch
stop.

From this calculation, we get six 2-particle clusters, two 3-particle clusters
and one 4-particle cluster (this is different from [12] where some “special
quasi-random structures”666Special quasi-random structures are supercells that
are designed to correctly reproduce the cluster structure. You can generate
them with the mcsqs-program that is distributed with ATAT. You provide the
maximum size of clusters of 2, 3 or more particles and the cell size you would
like and a Monte Carlo procedure is used to arrange atoms in the way that
optimally represents the cluster coefficients. are used as well to optimize
the cluster calculation.). You get this information most easily from the
clusinfo.out-file:777Be careful with this file; it is not created by maps but
only by mapsrep.

    2 2.489016 6 0.096360
    2 3.520000 3 0.002918
    2 4.311102 12 0.003826
    2 4.978032 6 -0.003782
    2 5.565609 12 -0.001564
    2 6.096819 4 0.004159
    3 2.489016 8 -0.005556
    3 3.520000 12 0.002930
    4 2.489016 2 -0.025176

First number is the number of particles, second the size, third the number of
equivalent configurations, and fourth the ECI. Note that the file eci.out
contains two more entries for the zero- and one-particle “cluster”.

Copy lat.in and all the out-files to a new directory (you can run in the same
directory without trouble, but I would recommend to keep files separated) and
run phb as usual:

    phb -gs1=0 -gs2=1 -dT=10 -dx=0.001 -er=50 -keV -ltep=1.e-3
     -o=phase01.out

This may take a while, but in the end you get a rather nice phase diagram of
the Ni-rich part of the phase diagram, see Fig. 9 (the yellow curve).

Figure 9: Calculated Ni-rich part of the Ni-Al phase diagram. The orange line
marked “all data” describes the first run in this section, the blue line is
discussed in section 7.3. The thick black line is the calculated phase diagram
using the TTNi8 database of Thermocalc [14].

Of course the high-temperature behaviour is wrong (melting is not included),
but overall the diagram does at least look realistic. It differs from a full
calculation done with thermodynamic software (I used Thermocalc [14] in the
figure), but it looks similar to the diagram shown in [12].

###  7.3 Remove distorted configurations

We did not check whether there are distorted configurations that may affect
the results of the computation. I recommend doing this by copying your whole
maps-directory (where all the DFT calculations reside) of the run to a new
name so that you can easily compare old and new runs. Create a
checkrelaxresult\- file with checkrelax > checkrelaxresult . It contains 5
configurations that have a value larger than 0.1:

    0.1062 609
    0.1090 51
    0.1119 13
    0.1666 24
    0.2220 608

To remove these, we create a file called error in the directory of the
corresponding run. For example, touch 609/error removes the result from the
most strongly distorted cell from the calculation. We can exclude all those
directories where checkrelax gives a value larger than 0.1 (the value
recommended in the man page) by simply editing the output from checkrelax,
replacing the first number with the string touch and appending an /error at
the end. I do this with an emacs keyboard macro, but you can also use awk, sed
or whatever your favorite tool for these things is, if the number is
annoyingly large. (However, if it is, it is probable that something is wrong
with your lattice definition as in the case of the NiAl-phase above.) If you
do this, don’t forget to re-run maps (and touch stop after it finishes) after
you touched the error files so that the new ECIs are actually
calculated.888When you re-run maps, the logfile may state that it only looked
at ground states with 0 atoms: The internal database of structures extends at
least up to 0 atoms/unit cell. The cluster calculation is nevertheless done
correctly. If you want to ensure that everything looks o.k., you can use the
-gs=-parameter of maps to specify the number of atoms in the ground state to
be looked at.

We can then do the same phase diagram calculation as before. The result is
also shown in fig. 9. The right-hand phase separation line is much closer to
the Thermocalc prediction; on the left-hand side, the low-temperature region
is a bit worse but the high-temperature behaviour looks quite good. Probably,
we could fine-tune the calculation here by adding some configurations from the
calculation, using SQS.

###  7.4 Phonon and electron influence

In this step, we follow the steps.txt-file that is in ATAT’s tutorial
directory to calculate temperature-dependent ECIs that contain the influence
of phonons and of electronic degrees of freedom. We use the version from the
previous section with removed distorted configurations.

Create a force.wrap-file, for example:

    [INCAR]
    PREC = high
    ISMEAR = 1
    SIGMA = 0.07
    NSW=0
    KPPRA = 1000
    USEPOT = PAWPBE
    NPAR=4
    ISPIN=2
    SUBATOM = s/Ni$/Ni_pv/g
    EDIFF = 1e-8
    EDIFFG = 1e-5

If you forget to create this file, pollmach will warn you: You need a
force.wrap file in one of the directories ., .. , ../.. , etc..

To find out the numbers of the ground states we need to calculate the phonons
for, look into gs.out. There are only two, one for pure Ni, one for ${\rm
Ni_{3}Al}$. For fitsvsl, I evaluate 3 different volumes instead of only 2 and
I set the distance of images to 8 because the manual states “Typically, -er
should be 3 or 4 times the nearest neighbor distance”. Possibly, a larger
value might work better; to get quantitatively reliable results, it would
probably be wise to try different values and check how they influence the
results.

    echo 0 27 > strname.in
    fitsvsl  -ns=3 -er=8
    pollmach -e runstruct_vasp -w force.wrap &
    fitsvsl -f
    gnuplot fitsvsl.gnu

The gnuplot-command plots the force-distance relation, it looks not perfect in
this case (figure 10), but reasonable.

Figure 10: Calculated force-distance relation (linear fit) from fitsvsl.

We continue as in the tutorial

    echo 2000 21 > Trange.in
    foreachfile -e str_relax.out pwd \; svsl -d
    clusterexpand -e svib_ht
    clusterexpand -e fvib

This creates the ECIs for the vibrational DOFs in the desired temperature
range. With

    mkteci fvib.eci

we can create a teci.out-file that contains temp-dependent ECIs.

Now also add electronic dofs:

    foreachfile -e str_relax.out pwd \; felec -d
    clusterexpand -e felec
    mkteci fvib.eci felec.eci

This first uses the foreachfile-utility to do a felec calculation in every
directory that contains str_relax.out. Note that the option -e ignores
directories that contain an error file; if this is not included and error file
are present, the command may stall.

For low temperatures, the temperature-dependent coefficients are similar to
those in eci.out as they should be:

teci.out | eci.out
---|---
-0.392042 | -0.425542
-0.057159 | -0.044139
0.116381 | 0.118159
-0.004468 | -0.004468
0.003976 | 0.003991
-0.00671667 | -0.006715
-0.00670783 | -0.006762
-0.0044705 | -0.004538
-0.00478617 | -0.004812
0.0134192 | 0.013498

The first value differs considerably, I assume that this is due to the zero-
point energy.

Although the phb-help page states that the ECI input file is eci.out, phb does
actually read a teci.out if provided. To calculate the phase diagram, we just
run it:

    phb -gs1=0 -gs2=1 -dT=10 -dx=0.001 -er=50 -keV -ltep=1.e-3
      -o=phase01.out > phase01-run

The result of the calculation is shown in figure 11. Up to a temperature of
1000 K, it agrees well with the previous calculations. However, at larger
temperatures, the phase separation line looks incorrect. It might be possible
that with the additional degrees of freedom, some intermediate phases are
stabilized (for example, the jump at a concentration of about 0.125 may be due
to a ${\rm Ni_{7}Al}$-phase). However, to find out exactly what is happening
here is beyond the scope of this tutorial.

Figure 11: Phase diagram calculated with phb using temperature-dependent
coefficients with phonon and electron contributions.

##  8 Conclusion

ATAT is a powerful tool to calculate phase boundaries and other material
properties. But as the experiments here have already shown, it is not a
program that can just be started and yields publishable results without user
intervention. Critical evaluation of all results, experimentation, and a
mixture of different approaches may be needed to understand even a simple
system like the ones we looked at here.

##  9 Acknowledgements

Thanks to Axel van der Walle for answering a lot of questions on ATAT and
reading the manuscript and to Jesper Kristensen for many helful hints on ATAT.

## References

  * [1] Axel Van De Walle and Mark Asta.  Self-driven lattice-model Monte Carlo simulations of alloy thermodynamic properties and phase diagrams.  Modelling and Simulation in Materials Science and Engineering, 10(5):521–538, 2002.
  * [2] Axel Van De Walle.  Multicomponent multisublattice alloys, nonconfigurational entropy and other additions to the Alloy Theoretic Automated Toolkit.  pages 1–25, 2009.
  * [3] Daan Frenkel and Berend Smit.  Understanding molecular simulation: from algorithms to applications, volume 1.  Elsevier, 2001.
  * [4] Dane Morgan.  Computational Laboratory : Monte Carlo for Phase Stability Calculations Atomistic Modeling Toolbox ( AMTB ).  2005\.
  * [5] A. van de Walle, M. Asta, and G. Ceder.  Talk "Automatic first-principles phase diagrams calculations".  https://www.brown.edu/Departments/Engineering/Labs/avdw/atat/atattalk.pdf, 2002\.
  * [6] T Williams, C Kelley, HB Bröker, J Campbell, R Cunningham, D Denholm, E Elber, R Fearick, C Grammes, and L Hart.  Gnuplot 4.5: An interactive plotting program. 2011.  http://www.gnuplot.info, 2017.
  * [7] F. Soisson and G. Martin.  Monte Carlo simulations of the decomposition of metastable solid solutions: Transient and steady-state nucleation kinetics.  Physical Review B - Condensed Matter and Materials Physics, 62(1):203–214, 2000.
  * [8] David A Porter, Kenneth E Easterling, and Mohamed Sherif.  Phase Transformations in Metals and Alloys, (Revised Reprint).  CRC press, 2009.
  * [9] Jesper Kristensen.  Vasp meets atat.  http://www.jespertoftkristensen.com/JTK/Blog/Entries/2013/11/7_Converting_str.out_to_POSCAR_and_vice_versa.html, 2014\.
  * [10] Ask Hjorth Larsen, Jens Jørgen Mortensen, Jakob Blomqvist, Ivano E Castelli, Rune Christensen, Marcin Dułak, Jesper Friis, Michael N Groves, Bjørk Hammer, Cory Hargus, et al.  The atomic simulation environment–a python library for working with atoms.  Journal of Physics: Condensed Matter, 29(27):273002, 2017.
  * [11] Jesper Kristensen.  How to visualize the clusters found in atat.  http://www.jespertoftkristensen.com/JTK/Blog/Entries/2013/11/8_How_to_visualize_the_clusters_found_in_ATAT.html, 2014\.
  * [12] C. Woodward, Axel Van De Walle, Mark Asta, and D.R. Trinkle.  First-principles study of interfacial boundaries in Ni–Ni3Al.  Acta Materialia, 75:60–70, 2014.
  * [13] Joachim Rösler, Harald Harders, and Martin Bäker.  Mechanical behaviour of engineering materials: metals, ceramics, polymers, and composites.  Springer Science & Business Media, 2007.
  * [14] Jan-Olof Andersson, Thomas Helander, Lars Höglund, Pingfang Shi, and Bo Sundman.  Thermo-calc & dictra, computational tools for materials science.  Calphad, 26(2):273–312, 2002.


provides a cleaner way of sampling subsets according to the Fourier
distribution than the related classical technique we provided in Section 5.
However, the issue arises when the nearest $k$-junta has Fourier mass on lower
levels (below $\log k$ or even a constant, for example). In this case, it is
not clear to us how quantum algorithms provide any advantage over classical
ones. An open question is whether quantum Fourier sampling techniques can be
applied in a more clever way to give faster algorithms in the tolerant testing
paradigm.

Finally, a clear open question is how good of a lower bound one can prove on
the query complexity of the tolerant junta testing problem. Our main result
Theorem 1.2, rules out strictly exponential-in-$k$ query lower bounds for
$k$-junta distance approximation. [PRW20] proved a non-adaptive query
complexity lower bound of $2^{k^{\eta}}$ for
$(k,k,\varepsilon_{1},\varepsilon_{2})$-tolerant junta testing (given a
particular choice of $0<\varepsilon_{1}<\varepsilon_{2}<1/2$), for any
$0<\eta<1/2$. While this is quite close to our upper bound of
$2^{\tilde{O}(\sqrt{k})}$, our algorithm is highly adaptive, while the lower
bound due to [PRW20] applies only to nonadaptive algorithms. Therefore,
another interesting direction would be to explore whether any nontrivial lower
bounds apply to adaptive algorithms for tolerant (junta) testing and distance
approximation.

### Acknowledgements

We thank Anindya De, Shafi Goldwasser, Amit Levi, and Orr Paradise for very
helpful discussions.

## References

  * [Ail+07] Nir Ailon, Bernard Chazelle, Seshadhri Comandur and Ding Liu  “Estimating the distance to a monotone function”  In _Random Structures & Algorithms_ 31.3, 2007, pp. 371–383  DOI: https://doi.org/10.1002/rsa.20167
  * [Amb+16] Andris Ambainis, Aleksandrs Belovs, Oded Regev and Ronald Wolf  “Efficient Quantum Algorithms for (Gapped) Group Testing and Junta Testing”  In _SODA_ SIAM, 2016, pp. 903–922
  * [BGS98] Mihir Bellare, Oded Goldreich and Madhu Sudan  “Free Bits, PCPs, and Nonapproximability-Towards Tight Results”  In _SIAM J. Comput._ 27.3, 1998, pp. 804–915
  * [Bla+19] Eric Blais, Clément L. Canonne, Talya Eden, Amit Levi and Dana Ron  “Tolerant Junta Testing and the Connection to Submodular Optimization and Function Isomorphism”  In _ACM Trans. Comput. Theory_ 11.4, 2019, pp. 24:1–24:33
  * [Bla08] Eric Blais  “Improved Bounds for Testing Juntas”  In _APPROX-RANDOM_ 5171, Lecture Notes in Computer Science  Springer, 2008, pp. 317–330
  * [Bla09] Eric Blais  “Testing juntas nearly optimally”  In _STOC_ ACM, 2009, pp. 151–158
  * [BLR90] Manuel Blum, Michael Luby and Ronitt Rubinfeld  “Self-Testing/Correcting with Applications to Numerical Problems”  In _STOC_ ACM, 1990, pp. 73–83
  * [BLT20] Guy Blanc, Jane Lange and Li-Yang Tan  “Testing and reconstruction via decision trees”  In _CoRR_ abs/2012.08735, 2020
  * [Bsh19] Nader H. Bshouty  “Almost Optimal Distribution-Free Junta Testing”  In _Computational Complexity Conference_ 137, LIPIcs  Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2019, pp. 2:1–2:13
  * [CG04] Hana Chockler and Dan Gutfreund  “A lower bound for testing juntas”  In _Inf. Process. Lett._ 90.6, 2004, pp. 301–305
  * [Cha+12] Sourav Chakraborty, Eldar Fischer, David Garcia-Soriano and Arie Matsliah  “Junto-Symmetric Functions, Hypergraph Isomorphism and Crunching”  In _Computational Complexity Conference_ IEEE Computer Society, 2012, pp. 148–158
  * [Che+18] Xi Chen, Rocco A. Servedio, Li-Yang Tan, Erik Waingarten and Jinyu Xie  “Settling the Query Complexity of Non-adaptive Junta Testing”  In _J. ACM_ 65.6, 2018, pp. 40:1–40:18
  * [Dia+07] Ilias Diakonikolas et al.  “Testing for Concise Representations”  In _FOCS_ IEEE Computer Society, 2007, pp. 549–558
  * [Dia+08] Ilias Diakonikolas, Homin K. Lee, Kevin Matulef, Rocco A. Servedio and Andrew Wan  “Efficiently Testing Sparse GF(2) Polynomials”  In _ICALP (1)_ 5125, Lecture Notes in Computer Science  Springer, 2008, pp. 502–514
  * [DMN19] Anindya De, Elchanan Mossel and Joe Neeman  “Junta Correlation is Testable”  In _FOCS_ IEEE Computer Society, 2019, pp. 1549–1563
  * [Fis+04] Eldar Fischer, Guy Kindler, Dana Ron, Shmuel Safra and Alex Samorodnitsky  “Testing juntas”  In _J. Comput. Syst. Sci._ 68.4, 2004, pp. 753–787
  * [GGR96] Oded Goldreich, Shafi Goldwasser and Dana Ron  “Property Testing and Its Connection to Learning and Approximation”  In _37th Annual Symposium on Foundations of Computer Science, FOCS ’96, Burlington, Vermont, USA, 14-16 October, 1996_ IEEE Computer Society, 1996, pp. 339–348  DOI: 10.1109/SFCS.1996.548493
  * [Hås01] Johan Håstad  “Some optimal inapproximability results”  In _J. ACM_ 48.4, 2001, pp. 798–859
  * [Kel+20] Esty Kelman, Subhash Khot, Guy Kindler, Dor Minzer and Muli Safra  “Theorems of KKL, Friedgut, and Talagrand via Random Restrictions and Log-Sobolev Inequality”  In _Electron. Colloquium Comput. Complex._ 27, 2020, pp. 9
  * [KKL88] Jeff Kahn, Gil Kalai and Nathan Linial  “The Influence of Variables on Boolean Functions (Extended Abstract)”  In _FOCS_ IEEE Computer Society, 1988, pp. 68–80
  * [KR00] Michael J. Kearns and Dana Ron  “Testing Problems with Sublearning Sample Complexity”  In _J. Comput. Syst. Sci._ 61.3, 2000, pp. 428–456
  * [Liu+19] Zhengyang Liu, Xi Chen, Rocco A. Servedio, Ying Sheng and Jinyu Xie  “Distribution-free Junta Testing”  In _ACM Trans. Algorithms_ 15.1, 2019, pp. 1:1–1:23
  * [LW19] Amit Levi and Erik Waingarten  “Lower Bounds for Tolerant Junta and Unateness Testing via Rejection Sampling of Graphs”  In _ITCS_ 124, LIPIcs  Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2019, pp. 52:1–52:20
  * [MOS03] Elchanan Mossel, Ryan O’Donnell and Rocco A. Servedio  “Learning juntas”  In _STOC_ ACM, 2003, pp. 206–212
  * [ODo14] Ryan O’Donnell  “Analysis of Boolean Functions”  Cambridge University Press, 2014  URL: http://www.cambridge.org/de/academic/subjects/computer-science/algorithmics-complexity-computer-algebra-and-computational-g/analysis-boolean-functions
  * [OW13] Ryan O’Donnell and Karl Wimmer  “Sharpness of KKL on Schreier graphs”  In _Electronic Communications in Probability_ 18.none Institute of Mathematical StatisticsBernoulli Society, 2013, pp. 1–12  DOI: 10.1214/ECP.v18-1961
  * [PRR04] Michal Parnas, Dana Ron and Ronitt Rubinfeld  “Tolerant Property Testing and Distance Approximation”  In _Electron. Colloquium Comput. Complex._ , 2004
  * [PRS01] Michal Parnas, Dana Ron and Alex Samorodnitsky  “Proclaiming Dictators and Juntas or Testing Boolean Formulae”  In _RANDOM-APPROX_ 2129, Lecture Notes in Computer Science  Springer, 2001, pp. 273–284
  * [PRW20] Ramesh Krishnan S. Pallavoor, Sofya Raskhodnikova and Erik Waingarten  “Approximating the Distance to Monotonicity of Boolean Functions”  In _SODA_ SIAM, 2020, pp. 1995–2009
  * [Sag18] Mert Saglam  “Near Log-Convexity of Measured Heat in (Discrete) Time and Consequences”  In _FOCS_ IEEE Computer Society, 2018, pp. 967–978
  * [Ser10] Rocco A. Servedio  “Testing by Implicit Learning: A Brief Survey”  In _Property Testing_ 6390, Lecture Notes in Computer Science  Springer, 2010, pp. 197–210
  * [Tal94] Michel Talagrand  “On Russo’s Approximate Zero-One Law”  In _The Annals of Probability_ 22.3 Institute of Mathematical Statistics, 1994, pp. 1576–1587  DOI: 10.1214/aop/1176988612
  * [Val12] Gregory Valiant  “Finding Correlations in Subquadratic Time, with Applications to Learning Parities and Juntas”  In _FOCS_ IEEE Computer Society, 2012, pp. 11–20
  * [Zha19] Xiaojin Zhang  “Near-Optimal Algorithm for Distribution-Free Junta Testing”  In _CoRR_ abs/1911.10833, 2019

##  Appendix A Maximum $k$-Subset Fourier Mass Approximation

In this section, we sketch a proof of Theorem 5.10, which involves simple
modifications and observations about our algorithm. The main difference is
that we sample from the normalized influence subdistribution at a different
Fourier level – namely, we let $\kappa:=\sqrt{k}$ and
$\alpha=k/\kappa=\sqrt{k}$ in Algorithm 4 and Algorithm 5, respectively
(recall that before, $\kappa=\sqrt{\varepsilon k}$). This improves the query
complexity dependence on $\varepsilon$ in Phase 1.

######  Claim A.1.

The query complexity of phase one of the algorithm for constant $\delta$
(failure probability) is $2^{\widetilde{O}(\sqrt{k}\log(1/\varepsilon))}$.

###### Proof.

The proof is analogous to the proof of 5.5, so we just point out the
differences. We still require our Fourier coefficients to be accurate to
within $1/\poly(k,1/\varepsilon)$, and we require confidence
$1-O(1/\ell)=1-2^{\tilde{\Omega}(\sqrt{k}\log(1/\varepsilon))}$. However, now
our branching process now has depth only $O(\sqrt{k})$, so we need only repeat
this $O(\ell)=2^{\widetilde{O}(\sqrt{k}\log(1/\varepsilon))}$ times, which
yields the improved query complexity. ∎

In Phase 2, we argue that it is not necessary to apply a noise operator in
order to only consider Fourier mass below level $\kappa$ after Phase 1. Recall
that we applied this noise operator in Section 5.2 in order to deal with the
case that a particular $U$ satisfied

$\mathop{\bf E\/}_{z\in\\{\pm
1\\}^{B}}\Bigl{[}\sum_{\begin{subarray}{c}S\subseteq U\setminus B,\\\
|S|>\kappa\end{subarray}}\;\widehat{f_{B\to
z}}(S)^{2}\Bigr{]}>\varepsilon^{2}/4.$ (15)

If this happened, then we could not rule out the possibility that taking the
low-degree part of $f$ within $U$ gives an overestimate to the correlation
with the best $k$-junta. However, now we are not concerned with the junta
correlation, but rather which set has the most mass, so we claim we do not
have to worry about this possibility anymore. To see this, suppose we have
identified $B\subseteq U$, and note that

$\displaystyle\sum_{S\subseteq U}\widehat{f}(S)^{2}$
$\displaystyle=\mathop{\bf E\/}_{x}[f(x)f_{\avg,U}(x)]$
$\displaystyle=\mathop{\bf E\/}_{z\in\\{\pm 1\\}^{B}}\bigg{[}\mathop{\bf
E\/}_{x}[f_{B\to z}(x)(f_{\avg,U})_{B\to z}(x)]\bigg{]}$
$\displaystyle=\mathop{\bf E\/}_{z}\bigg{[}\sum_{\begin{subarray}{c}S\subseteq
U\\\ |S|\leq\kappa\end{subarray}}\widehat{f_{B\to
z}}(S)^{2}+\sum_{\begin{subarray}{c}S\subseteq U\\\
|S|>\kappa\end{subarray}}\widehat{f_{B\to z}}(S)^{2}\bigg{]}$
$\displaystyle\geq\mathop{\bf
E\/}_{z}\bigg{[}\sum_{\begin{subarray}{c}S\subseteq U\\\
|S|\leq\kappa\end{subarray}}\widehat{f_{B\to z}}(S)^{2}\bigg{]}.$

Therefore, we no longer have to apply any noise operator, which negates the
necessity of 5.6 and 5.7. It therefore suffices in Lemma 5.8 to estimate the
mass of each set, rather than the correlation, as

$m_{U}=\mathop{\bf E\/}_{z}\bigg{[}\sum_{\begin{subarray}{c}S\subseteq U\\\
|S|\leq\kappa\end{subarray}}\widehat{f_{B\to z}}(S)^{2}\bigg{]}.$

To do so, as in the proof of Lemma 5.8 we let
$t=O(\log(1/\delta)/\varepsilon^{2})$ be the number of random samples of $z$
we take. Then we estimate all the Fourier coefficients below level $\kappa$.
This requires estimating $\widehat{f}(S)$ for all $S\subseteq\mathcal{S}$ of
size at most $\kappa$ up to additive error
$\varepsilon/\binom{k}{\leq\kappa}=2^{\widetilde{\Omega}(\sqrt{k}\log(1/\varepsilon))}$
with probability $1-\frac{\delta}{t\cdot\binom{k}{\leq\kappa}}$, which is
possible via 2.1 with
$\log(1/\delta)2^{\tilde{O}(\sqrt{k}\log(1/\varepsilon))}$ queries. The rest
of our argument and algorithm is exactly the same as in Section 5.