Hi Miles, I got it, and my program calculates know without error, but I think there is still a mistake, because I get implausible values.
I discovered the correlationmatrix function:
My question: is it possible to use the correlationmatrix function for i,i+1,j,j+1 → in total 4 operators on different sites ?
In my case, ladder look like this (d&p orbitals/sites):
d-p-d-p-d -p -d \\ | \quad \qquad \ |\qquad \quad \ |\qquad \quad \ |\\ p \quad \qquad p \qquad \quad p \qquad \quad p \\ | \quad \qquad \ |\qquad \quad \ |\qquad \quad \ |\\ d-p-d-p-d -p -d
Site numbers:
3-4-8-9-\ \ 13 -14 -18 \\
| \quad \qquad \ |\qquad \quad \ \ \ \ |\qquad \quad \ \ \ \ \ |\\
2 \quad \qquad 7 \qquad \quad \ \ 12 \qquad \quad \ \ 17 \\
| \quad \qquad \ |\qquad \quad \ \ \ \ |\qquad \quad \ \ \ \ \ |\\
1-5-6-10-11 -15 -16
I have to calculate this correlation function P(r)=0.5 (<\Delta^+_i \Delta_j>+<\Delta_i \Delta^+_j>) with \Delta^+_i=2^{-1/2} (d^+_{i,1,\uparrow}d^+_{i,2,\downarrow}-d^+_{i,1,\downarrow}d^+_{i,2,\uparrow}) and \Delta_i=2^{-1/2} (d_{i,2,\downarrow}d_{i,1,\uparrow}-d_{i,2,\uparrow}d_{i,1,\downarrow}).
That’s why I’m asking about the possibility to use the correlationmatrix function for my problem, it would save me time, but I don’t see a way to use it like this 
BR
Gökmen