  A search for Higgs boson pair production in the $b\bar{b}b\bar{b}$ final
state is carried out with up to 36.1 $\mathrm{fb}^{-1}$ of LHC proton-proton
collision data collected at $\sqrt{s}$ = 13 TeV with the ATLAS detector in 2015
and 2016. Three benchmark signals are studied: a spin-2 graviton decaying into
a Higgs boson pair, a scalar resonance decaying into a Higgs boson pair, and
Standard Model non-resonant Higgs boson pair production. Two analyses are
carried out, each implementing a particular technique for the event
reconstruction that targets Higgs bosons reconstructed as pairs of jets or
single boosted jets. The resonance mass range covered is $260-3000$ GeV. The
analyses are statistically combined and upper limits on the production cross
section of Higgs boson pairs times branching ratio to $b\bar{b}b\bar{b}$ are
set in each model. No significant excess is observed; the largest deviation of
data over prediction is found at a mass of 280 GeV, corresponding to 2.3
standard deviations globally. The observed 95% confidence level upper limit on
the non-resonant production is 13 times the Standard Model prediction.
