A model for protein in 2D can be formed by adding bonds of fixed length $l\sqrt{2}$ on a square lattice along the diagonal, ie $\hat{\mathbf{b}}_i=\frac{1}{\sqrt{2}}(\pm \hat{\mathbf{x}}\pm \mathbf{y})$. Consider a protein made up of $N$ bonds $\hat{\mathbf{b}}_i=\frac{1}{\sqrt{2}}(x_i,y_i)$ where $x_i=\pm 1, y_i=\pm 1$ and $i \in \{1,\dots,N\}$.
The energy of the underformed protein is given by
$$\mathcal{H}_0=-\sum_{i \neq j}K \hat{\mathbf{b}}_i \cdot \hat{\mathbf{b}}_j=-\sum_{i \neq j}\frac{K}{2}(x_ix_j+y_iy_j)$$
where the sum above is over all pairs $i, j \in \{1,\dots,N \}$. The displacement between the two ends of the protein is $\mathbf{R}=l\sqrt{2}\sum_{i=1}^N \hat{\mathbf{b}}_i$.
Under a tensile force $\hat{\mathbf{f}}=f\hat{\mathbf{n}}$, $f>0$, is the force aligned along the vector $\hat{\mathbf{n}}=(\cos \theta, \sin \theta)$, the system is described by the Hamiltonian $$\mathcal{H}=\mathcal{H}_0-\mathbf{f} \cdot \mathbf{R}=\mathcal{H}_0 -fl \sum_{i=1}^N (x_i \cos \theta + y_i \sin \theta)$$
Finite stiffness. Consider the case $K \neq 0$ for the long protein made of $N \gg 1$ bonds under a tensile force aligned along the x-axis:
It can be shown that the Hamiltonian can be written in the form $$\mathcal{H}(\{x_i,y_i \})=\mathcal{H}_0-\mathbf{f} \cdot \mathbf{R}=-l\sqrt{2}\sum_{i=1}^N \mathbf{f}_i \cdot \hat{\mathbf{b}}_i+E_0$$
where $\displaystyle \mathbf{f}_i=\sum_{j=1}^N \frac{K}{l\sqrt{2}}\hat{\mathbf{b}}_j-\mathbf{f}$
and it is easy to show that $\langle \mathbf{f}_i\rangle=\mathbf{f}+\Delta \mathbf{f}$ and $\Delta \mathbf{f}=KA \langle \mathbf{R} \rangle$ where $A=\frac{1}{2l^2}$.
A mean field Weiss approximation involves replacing the instantaneous force felt by each bond by its average. I am then asked to find the mean field Hamiltonian.
I get the mean field Hamiltonian will be $$\begin{align}\mathcal{H}_{mf} &= -l\sqrt{2}\sum_{i=1}^N (\mathbf{f}+\Delta \mathbf{f})\cdot \hat{\mathbf{b}}_i+E_0 \\ &= -l\sqrt{2} \left(\frac{K}{2l^2}\langle \mathbf{R} \rangle + \mathbf{f}\right)\cdot \sum_{i=1}^N \hat{\mathbf{b}}_i+E_0 \\ \end{align}$$
I then try to calculate the Helmholtz free energy $F=-k_BT \ln Z$
We may drop the $E_0$ so
$$\begin{align}Z &= \sum_{\Gamma}e^{-\beta \mathcal{H}_{MF}(\Gamma)} \\ &= \sum_{\Gamma}e^{-\beta \left(-l\sqrt{2} \left(\frac{K}{2l^2}\langle \mathbf{R} \rangle + \mathbf{f}\right)\cdot \sum_{i=1}^N \hat{\mathbf{b}}_i\right)} \end{align}$$
but this isnt very tractable. Where are am I going wrong?