Matlab code to Mathematica, Tensor

Question

I am trying to convert a Matlab code to Mathematica. I have a tensor, RF (a,b,c,d) which is constructed from Hamiltonian (Hs) and 2x2 identity matrix (DS):

RF(a,b,c,d)=Gamma(d,b,a,c)+conj(Gamma(c,a,b,d))-GaP-GaN

where a, b , c and d are the different position of Hs. In Matlab, it is codded as following as:

N=2
%-----------------------------------------------------------
Hs=[-0.3165 , 0.08;0.08,-0.1582]
%---------------------------------------------------------------

%--------------Diagonalization of the Hamiltonian---------------
DS=eye(N,N) ;
DSt=DS' ;
[U,E]=eig(Hs) ;
Ut=U' ;
U_1=inv(U);
%---------------------------------------------------------------

%------------------- Tensor1 -------------------------------------
for a=1:N
    for b=1:N
        for c=1:N
            for d=1:N
                 P=0 ;
                    for n=1:N %Run over the sites n       
                      P=P+Ut(a,:)*DS(:,n)*DSt(n,:)*U(:,b)*Ut(c,:)*DS(:,n)*DSt(n,:)*U(:,d) ;
                    end
                    if  d==c
                        Gamma(a,b,c,d)=P*(E(d,d)-E(c,c));
                    else
                        Gamma(a,b,c,d)=P*(E(d,d)-E(c,c));
                    end;
            end
        end
    end
end

%------------------- Tensor2 -------------------------------------
for a=1:N
    for b=1:N
        for c=1:N
            for d=1:N                
               GaP=0 ;
               if b==d          
                   for i=1:N
                       GaP=GaP+Gamma(a,i,i,c) ;
                   end
               end

               GaN=0 ;
               if a==c
                   for i=1:N
                       GaN=GaN+conj(Gamma(b,i,i,d)) ;
                   end
               end
               RF(a,b,c,d)=Gamma(d,b,a,c)+conj(Gamma(c,a,b,d))-GaP-GaN ;
            end
        end
    end
end    

In Mathematica, I have tried to code this in following way,

1st-- I diagonalize the Hamiltonian and obtained eigen-vectors (U) with corresponding eigen-values in HLS.

2nd-- I defined a function f1[b_, a_] which is equivalent to (E(d,d)-E(c,c)) in Matlab.

3rd-- I defined a function f[a_, b_, c_, d_] which equivalent to Gamma(a,b,c,d) in Matlab.

4rth -- Using Table Function, I evaluated the f[a_, b_, c_, d_] over range on n. At this point I have the values for Gamma(a,b,c,d).

Clear["Global`*"];
n=2;
Hs = {{-0.3165 , 0.08},{0.08,-0.1582}};
{eigsHS, vecsHS} = Eigensystem[Hs];
U = Transpose[vecsHS];
Ut = Transpose[U];
HLS = DiagonalMatrix[eigsHS];
f1[b_, a_] := Indexed[HLS, {b, b}] - Indexed[HLS, {a, a}];
f[a_, b_, c_, d_] := 
  Sum[Indexed[Ut, {a}].Indexed[DS, {m}]*
     Indexed[DSt, {m}].Indexed[Ut, {b}]*
     Indexed[Ut, {c}].Indexed[DS, {m}]*
     Indexed[DSt, {m}].Indexed[Ut, {d}], {m, n}]*(f1[a, b]);
Gamma(a,b,c,d) = Table[f[a, b, c, d], {a, n}, {b, n}, {c, n}, {d, n}];

But, I am unable to translate second part of Matlab code %-----Tensor2, to finally arrive at RF(a,b,c,d). So far, I have tried following option, using If condition in Table, but the output is different.

tesnor2 = 
 Table[If[b == d, f[a, b, c, d], Nothing], {a, n}, {b, n}, {c, n}, {d,
    n}]

I will be really grateful for any help.

Answer 1

Maybe this is a bit closer to what you want to do.

RF = Table[
   Plus[
    f[d, b, a, c],
    Conjugate[f[c, a, b, d]],
    If[b == d, -Sum[f[a, i, i, c], {i, 1, n}], 0],
    If[a == c, -Conjugate[Sum[f[b, i, i, d], {i, 1, n}]], 0]
    ]
   , {a, 1, n}, {b, 1, n}, {c, 1, n}, {d, 1, n}];

Note that I use Plus instead of + only because it is more legible for code that is spread over several lines.

Don't use Nothing in this case. Its effect is that every occurence of it will be removed from an embracing list. So the result will be a ragged array, not a tensor.