Let be $$H_0= -J \sum_{x} \sigma_{x}\sigma_{x+e_1}+\sigma_{x}\sigma_{x+e_2}$$ with $e_1=(1,0)$ and $e_2=(0,1)$.
The representation of partition function in terms of Grassmann integral is $$Z=\sum_{\{\sigma\}}e^{-\beta H_0}=(cosh(\beta J))^B2^s\int \prod_{x}dH_xd\bar{H}_xdV_xd\bar{V}_xe^S$$ where $S=\sum_{x}(tanh (\beta J))[\bar{H}_x H_{x+e_1}+\bar{V}_x V_{x+e_2}]+\sum_{x}[\bar{H}_x H_x+\bar{V}_x V_x +\bar{V}_x \bar{H}_x+ V_x\bar{H}_x+H_x\bar{V}_x+V_xH_x ]$.
Now if one considerer the coupled terms dependent by coordinate one can say
$$\sigma_x\sigma_{x+e_1}e^{-\beta H_0}=\frac{\partial}{\partial\beta J_{x,x+e_1}}e^{-\beta H_0[J]}\bigg|_{J=1}$$
So one can write
$$ \begin{align*}\sigma_x\sigma_{x+e_1}\sigma_{x+e_1}\sigma_{x+2e_1}e^{-\beta H_0}&=\frac{\partial}{\partial \beta J_{x,x+e_1}}\frac{\partial}{\partial \beta J_{x+e_1,x+2e_1}}e^{-\beta H_0[J]}\bigg|_{J=1} \\ =\frac{\partial}{\partial \beta J_{x,x+e_1}} \frac{1}{cosh^2(\beta J_{x+e_1,x+2e_1})}\bar{H}_{x+e_1}H_{x+2e_1}e^{-\beta H_0[J]}\bigg |_{J=1} \end{align*}$$
My problem is: how obtain the last equality? I am very in trouble with this...