function [PotEnergija,KinEnergija,Virijal,novekoordinate,brzine,radial]=md(mm,koordinate,starekoordinate,brojcestica,Temperatura,Ecutoff,Rcutoff,Lbox,dt,poc_uzorkovanje,bin_radial,maxx_bin_radial)
    dx=ones(brojcestica,1)*koordinate(1,:)-koordinate(1,:)'*ones(1,brojcestica);
    dy=ones(brojcestica,1)*koordinate(2,:)-koordinate(2,:)'*ones(1,brojcestica);
    dz=ones(brojcestica,1)*koordinate(3,:)-koordinate(3,:)'*ones(1,brojcestica);
    % uslov periodicne granice
    dx=dx-Lbox*round(dx./Lbox);
    dy=dy-Lbox*round(dy./Lbox);
    dz=dz-Lbox*round(dz./Lbox);
    dr2=dx.^2+dy.^2+dz.^2;
    drm2=1./(dr2+eye(brojcestica));
    % Proveriti da li je rastojanje izmedju 
    % cestica vece od granicnog rastojanja; ako je vece zameniti ga sa 0
    drm2(drm2<1/(Rcutoff^2))=0;
    brojnula=sum(sum(drm2==0));
    drm6=drm2.^3;
    drm6(logical(eye(brojcestica,brojcestica)))=0;
    drm8=drm2.^4;
    drm8(logical(eye(brojcestica,brojcestica)))=0;
    ff=48*drm8.*(drm6-0.5*(~eye(brojcestica,brojcestica)));
    F(1,:)=sum(ff.*dx); 
    F(2,:)=sum(ff.*dy);
    F(3,:)=sum(ff.*dz);
    PotEnergija=sum(sum(triu(4*drm6.*(drm6-ones(brojcestica,brojcestica)))))-(brojcestica*(brojcestica-1)/2-brojnula/2)*Ecutoff;
    Virijal=sum(sum(triu(48*drm6.*(drm6-0.5*ones(brojcestica,brojcestica))))); 
    % odrediti nove koordinate i brzine cestica 
    novekoordinate=2*koordinate-starekoordinate+F*(dt^2);
    brzine=(novekoordinate-starekoordinate)/(2*dt);
    KinEnergija=0.5*sum(sum(brzine.^2));
    if (mm<poc_uzorkovanje) 
        scale=sqrt(Temperatura*(3*brojcestica-3)/(2*KinEnergija));
        brzine=brzine*scale;
        novekoordinate=starekoordinate+brzine*2*dt;
        KinEnergija=0.5*sum(sum(brzine.^2));
        radial=zeros(1,maxx_bin_radial);
    else
        %%% deo za radijalnu funkciju raspodele      
        dr=floor(sqrt(dr2)/(Lbox/(2*maxx_bin_radial)));
        radial=sum(histc(dr,bin_radial),2)';
    end
end

