Boltzmann distribution in Ising model I've written in Matlab a code for a Ising model in 1 dimension with 40 spins at $k_{B}T=1$.
I record the energy of every step in a Metropolis Monte Carlo algorithm, and then I made an histogram like this.

I want to show the theoretical Boltzmann distribution. What is the exact formula to get this shape? $Ax \exp(-\beta x)$?
 A: I have to make a number of assumptions, as you did not state all the necessary information. So, I am going to assume that you are using periodic boundary conditions, that is, your Hamiltonian is
$$
\mathcal{H}(\sigma) = -\sum_{i=1}^N \sigma_i\sigma_{i+1},
$$
where I have denoted by $N$ the number of spins (that is, $N=40$ in your case) and used the convention that $\sigma_{N+1}=\sigma_1$.
Write $D(\sigma)$ the number of $i\in \{1,\dots,N\}$ such that $\sigma_i\neq \sigma_{i+1}$, still using the convention that $\sigma_{N+1} = \sigma_1$. Note that $D(\sigma)$ is necessarily an even number (because of the periodic boundary conditions).
Then, the total energy can be rewritten as
$$
\mathcal{H}(\sigma) = D(\sigma) - (N-D(\sigma)) = 2D(\sigma) - N \,.
$$
The probability that $D(\sigma)=\delta N$ (with $\delta \in \{0,\frac{2}{N},\frac{4}{N},\dots,\frac{2\lfloor N/2 \rfloor}{N}\}$) is
$$
\mathrm{Prob}(D = \delta N) = \binom{N}{\delta N}\frac{\exp(-2\beta\delta N)}{\frac{1}{2}\bigl(1-\exp(-2\beta)\bigr)^N + \frac{1}{2}\bigl(1+\exp(-2\beta)\bigr)^N},
$$
since there are $\binom{N}{\delta N}$ ways of choosing the $\delta N$ pairs of disagreeing neighbors.
This can be easily reformulated in terms of the energy per spin
$$
\frac{1}{N}\mathcal{H}(\sigma) = \frac{2}{N}D(\sigma) - 1.
$$
Note that the possible values of the latter are $-1$, $-1+\frac{4}{N}$, $-1+\frac{8}{N}$, ... , $-1+\frac{4\lfloor N/2\rfloor}{N}$.
The probability of observing an energy per spin equal to $\epsilon$ is then given by
$$
\mathrm{Prob}(\mathcal{H}=\epsilon N)
=
\binom{N}{\frac{1+\epsilon}{2}N}
\frac{\exp\bigl(-\beta (1+\epsilon) N\bigr)}{\frac{1}{2}\bigl(1-\exp(-2\beta)\bigr)^N + \frac{1}{2}\bigl(1+\exp(-2\beta)\bigr)^N} .
$$
Here is a plot of the distribution for your parameters $N=40$ and $\beta=1$ (only values of $\epsilon$ smaller than $-0.2$ are indicated as higher values have too small probability at this temperature):

(The computation using other boundary conditions are similar.)
