% This function calculates the components of the heat current vector for thermal conductivity calculations.
function heatcurrent(xn, pn, t, Lo2, m, a, partiden, N, energy_cutoff, heatCurrentFile);

L = zeros(1,3);
q = zeros(1,3);
qconv = zeros(1,3);
qcond = zeros(1,3);
deltaR = zeros(1,3);
vE = zeros(1,3);
vsum = zeros(1,3);
vn = zeros(N,3);
PHI = zeros(1,N);

for i = 1:3 L(1,i) = Lo2(1,i)*2; end;

for i = 1:N
  for j = 1:3 vn(i,j) = pn(i,j)/m(1,i); end;
end;


for i = 1:N
  for j = i+1:N
    for k = 1:3 vsum(1,k) = vn(i,k) + vn(j,k); end;

    for k = 1:3
      deltaR(1,k) = xn(i,k) - xn(j,k);
      if (deltaR(1,k) > Lo2(1,k)) deltaR(1,k) = deltaR(1,k) - L(1,k); end;
      if (deltaR(1,k) < (-1.*Lo2(1,k))) deltaR(1,k) = deltaR(1,k) + L(1,k); end; 
    end;

    if ((deltaR(1,1) < a) && (deltaR(1,2) < a) && (deltaR(1,3) < a)) r = partsep(deltaR);
    else r = 10*a;
    end;

    % Short range contribution.
    if(r<a)
      Fdot_plus = 0.;
	  for k = 1:3 Fdot_plus = Fdot_plus + deltaR(1,k)*vsum(1,k); end;
      Fdot_plus = Fdot_plus * LJForce(r);
      for k = 1:3 qcond(1,k) = qcond(1,k) + deltaR(1,k)*Fdot_plus; end;

      dE = LJEnergy(r, energy_cutoff);
      PHI(1,i) = PHI(1,i) + dE;
      PHI(1,j) = PHI(1,j) + dE;
	end;
  end;
end;

for i= 1:3 qcond(1,i) = qcond(1,i)/2.; end;

for i = 1:N
  PHI(1,i) = PHI(1,i)/2. + 1./(2.*m(1,i))*(pn(i,1)*pn(i,1) + pn(i,2)*pn(i,2) + pn(i,3)*pn(i,3));
  for j = 1:3 qconv(1,j) = qconv(1,j) + PHI(1,i)*vn(i,j); end;
end;

for i = 1:3 q(1,i) = qcond(1,i) + qconv(1,i); end;

averageHeat = sqrt(q(1,1)^2+q(1,2)^2+q(1,3)^2);

fprintf(heatCurrentFile, '%d\t %g\t %g\t %g\t %g\n', t, q(1,1), q(1,2), q(1,3), averageHeat);
