% Program "MDargon_NVE"

clc; close all; clear; tic

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% PARAMETRI SIMULACIJE %%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

NParticles=XXX;                % Broj cestica u simulaciji
ro=XXX;                        % Cesticna gustina = NParticles/Lbox^3
T=XXX;                         % Temperatura
dt=XXX;                        % Duzina vremenskog koraka 
Nsteps=XXX;                    % Broj koraka u simulaciji
Nsample=XXX;                   % Broj koraka za racunanje srednje vrednosti

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Lbox=(NParticles/ro)^(1/3); % duzina simulacione kutije
Rcutoff=0.4999*Lbox;  % rastojanje iznad kojeg se ne racuna sila 
Ecutoff=4*(1/(Rcutoff^12)-1/(Rcutoff^6)); % energija na rastojanju Rcutoff
Epot_tail=8/3*pi*ro*(1/(3*Rcutoff^9)-1/(Rcutoff^3)); % korekcija za potencijalnu energiju
P_tail=16/3*pi*ro^2*(2/(3*Rcutoff^9)-1/(Rcutoff^3)); % korekcija za pritisak
maxx_bin_radial=250; % broj kutijica za radijalnu funk.
Delta_g=Lbox/(2*maxx_bin_radial); %  velicina kutijica za radijalnu funk. raspodele
bin_radial=0:maxx_bin_radial-1; % kutijice
g_radial=zeros(1,maxx_bin_radial); % inicijalizacija radijalne funk.
if (Rcutoff<1)
    error('upozorenje:Rcutoff<1')
end
if (Nsteps<Nsample)
    error('upozorenje:Nsteps<Nsample')
end

fprintf('\n')
fprintf('-----------------------------------------------------------------------------------------------\n')
fprintf('\n')
fprintf ( 1, 'Program za proracun termodinamickih osobina argona pomocu molekulske dinamike u NVE ansamblu. \n' );
fprintf ( 1, 'Program koristi redukovane jedinice (r.j.). \n' );
fprintf ( 1, 'Potencijalna energija se racuna koristeci ograniceni i pomereni LJ potencijal. \n' );
fprintf ( 1, 'Korekcija za veca rastojanja se dodaje na kraju. \n' );
fprintf ( 1, 'Parametri LJ potencijala za argon su: %s = %e m, %s = %e J i m = %e kg \n',char(963), 3.405e-10, char(949), 1.654e-21, 6.68e-26);
fprintf('\n')
fprintf('-----------------------------------------------------------------------------------------------\n')
fprintf('\n')
fprintf ( 1, 'Podaci o simulaciji: \n' );
fprintf ( 1, '\n' );
fprintf ( 1, 'Broj atoma u simulaciji je:                      %d \n', NParticles );
fprintf ( 1, 'Gustina atoma je:                                %f r.j. = %f kg/m3 \n', ro, 1688*ro );
fprintf ( 1, 'Temperatura je:                                  %f r.j. = %f K \n', T,T*119.91 );
fprintf ( 1, 'Ukupni broj koraka u simulaciji je:              %d \n', Nsteps );
fprintf ( 1, 'Broj koraka za uravnotezenje je:                 %d \n', Nsteps-Nsample );
fprintf ( 1, 'Broj koraka za uzorkovanje je:                   %d \n', Nsample );
fprintf ( 1, 'Duzina koraka u simulaciji je:                   %f r.j. = %f fs \n', dt, dt*2170 );
fprintf ( 1, 'Duzina sumulacione kutije je:                    %f r.j. = %f nm \n', Lbox, 0.3405*Lbox );
fprintf ( 1, 'Rastojanje iznad kojeg se ne racuna sila je:     %f r.j. = %f %s \n', Rcutoff, 0.3405*Rcutoff/10, char(197) );
fprintf ( 1, 'Potencijalna energija na tom rastojanju je:     %f r.j. \n', Ecutoff );
fprintf ( 1, 'Korekcija za potencijalnu energiju po atomu je: %f r.j. \n', Epot_tail );
fprintf ( 1, 'Korekcija za pritisak je:                       %f r.j. \n', P_tail );

%% definisati nizove koji ce se koristiti
coords=zeros(3,NParticles,Nsteps+1);
Velocities=zeros(3,NParticles,Nsteps);
Ekin=zeros(1,Nsteps);
Epot=zeros(1,Nsteps);
Virial=zeros(1,Nsteps);

% postavi atome na resetku
ni = ceil(NParticles^(1.0/3.0));
ncount = 0;
dx = Lbox/ni;
for ix = 1:ni
 for iy = 1:ni
  for iz = 1:ni
   ncount = ncount + 1;
   if (ncount <= NParticles)
    coords(1,ncount,2) = dx*ix;
    coords(2,ncount,2) = dx*iy;
    coords(3,ncount,2) = dx*iz;
   end
  end
 end
end

% dati mali slucajan pomeraj svakom atomu
coords(:,:,2)=coords(:,:,2)-0.05*(2*rand(3,NParticles)-1);
coords(:,:,2)=coords(:,:,2)-Lbox*(coords(:,:,2)>Lbox);
coords(:,:,2)=coords(:,:,2)+Lbox*(coords(:,:,2)<0);

%%% dati slucajne brzine iz Gausove raspodele 
Velocities(:,:,1)=randn(3,NParticles,1);

% oduzeti ukupni moment(=brzina) 
Velocities(1,:,1)=Velocities(1,:,1)-sum(Velocities(1,:,1))/NParticles;
Velocities(2,:,1)=Velocities(2,:,1)-sum(Velocities(2,:,1))/NParticles;
Velocities(3,:,1)=Velocities(3,:,1)-sum(Velocities(3,:,1))/NParticles;
% skalirati brzine cestica tako da ukupna kineticka energija odgovara
% zadatoj temperaturi
Ukin0=0.5*sum(sum(Velocities(:,:,1).^2));
scale=sqrt(T*(3*NParticles-3)/(2*Ukin0));
Velocities(:,:,1)=Velocities(:,:,1)*scale;
    
% odrediti polozaj u prvom trenutku (uzimamo da se cestice krecu po inerciji bez sila)
coords(:,:,1)=coords(:,:,2)-Velocities(:,:,1)*dt;

if ((Nsteps-Nsample+1)>2)
    poc_uzorkovanje=Nsteps-Nsample+1;
    Nuzoraka=Nsample;
else
    poc_uzorkovanje=2;
    Nuzoraka=Nsteps-1;
end

for m=2:Nsteps
    koordinate=coords(:,:,m);
    starekoordinate=coords(:,:,m-1);
    [PotEnergija,KinEnergija,Virijal,novekoordinate,brzine,radial]=md(m,koordinate,starekoordinate,NParticles,T,Ecutoff,Rcutoff,Lbox,dt,poc_uzorkovanje,bin_radial,maxx_bin_radial);   
    Epot(m)=PotEnergija;
    Ekin(m)=KinEnergija;
    Virial(m)=Virijal;
    g_radial=g_radial+radial;
    coords(:,:,m+1)=novekoordinate;
    Velocities(:,:,m)=brzine;
end

Etot=Epot+Ekin;

if (poc_uzorkovanje>2)
    Pomeraj_energije=sum(abs(((Etot(Nsample+1)-Etot(Nsample+2:Nsteps))/Etot(Nsample+1))))/(Nsteps-Nsample-1);
else
    Pomeraj_energije=sum(abs(((Etot(2)-Etot(3:Nsteps))/Etot(2))))/(Nsteps-2);
end

 figure  
 plot(1:Nsteps,Epot/NParticles,'k-',1:Nsteps,Ekin/NParticles,'r-',1:Nsteps,Etot/NParticles,'b-',...
     'LineWidth',0.6)
 xlabel('\bf \fontsize{15} \fontname{Euclid} Broj koraka');
 ylabel('\bf \fontsize{15} \fontname{Euclid} E_{pot}, E_{kin}, E_{tot} po cestici');
 legend('\fontname{Euclid} E_{pot}/cestici','\fontname{Euclid} E_{kin}/cestici','\fontname{Euclid} E_{tot}/cestici')

Etot_pocestici=sum(Etot(1,poc_uzorkovanje:Nsteps))/(Nuzoraka*NParticles);
Epot_pocestici=sum(Epot(1,poc_uzorkovanje:Nsteps))/(Nuzoraka*NParticles);
Ekin_pocestici=sum(Ekin(1,poc_uzorkovanje:Nsteps))/(Nuzoraka*NParticles);
Av_temperature=2*Ekin_pocestici*NParticles/(3*NParticles-3);
Av_pressure=ro*Av_temperature+(sum(Virial(1,poc_uzorkovanje:Nsteps))/Nuzoraka)/(3*Lbox^3);

fprintf ( 1, '\n' );
fprintf('-----------------------------------------------------------------------------------------------\n')
fprintf ( 1, '\n' );
fprintf ( 1, 'Rezultati simulacije: \n');
fprintf ( 1, '\n' );
fprintf ( 1, 'Srednja korigovana ukupna energija po cestici je:       %f r.j. = %f kJ/mol \n', Etot_pocestici+Epot_tail,(Etot_pocestici+Epot_tail)*0.997);
fprintf ( 1, 'Srednja korigovana potencijalna energija po cestici je: %f r.j. = %f kJ/mol \n', Epot_pocestici+Epot_tail, (Epot_pocestici+Epot_tail)*0.997);
fprintf ( 1, 'Srednja kineticka energija po cestici je:                %f r.j. =  %f kJ/mol \n', Ekin_pocestici, Ekin_pocestici*0.997);
fprintf ( 1, 'Srednja temperatura sistema je:                          %f r.j. = %f K \n', Av_temperature, Av_temperature*120);
fprintf ( 1, 'Srednji korigovani pritisak je:                         %f r.j. = %f MPa \n', Av_pressure+P_tail,(Av_pressure+P_tail)*42);
fprintf ( 1, 'Srednji relativni pomeraj ukupne energije je:            %f r.j. \n', Pomeraj_energije);
fprintf ( 1, 'Ukupni impuls na kraju simulacije je:                    %e r.j., %e r.j., %e r.j. \n', sum(Velocities(1,:,Nsteps)),sum(Velocities(2,:,Nsteps)),sum(Velocities(3,:,Nsteps)));

%%% deo za radijalnu funkciju raspodele
  g_radial(1)=g_radial(1)-Nuzoraka*NParticles;
  rr=(0.5:1:(maxx_bin_radial-0.5))*Delta_g;
  g_rad_final=g_radial./(Nuzoraka*NParticles*(4*pi/3*NParticles/(Lbox^3)*(Delta_g^3)*((1:maxx_bin_radial).^3-(0:maxx_bin_radial-1).^3)));
  
  figure
  plot(rr,g_rad_final,'k-',...
    'LineWidth',1.5);
  xlabel('\bf \fontsize{15} \fontname{Euclid} r');
  ylabel('\bf \fontsize{15} \fontname{Euclid} g(r)');
  legend('\fontname{Euclid} Radijalna funkcija raspodele')
  
  figure
  plot(rr,-Av_temperature*log(g_rad_final),'k-',...
    'LineWidth',1.5);
  xlabel('\bf \fontsize{15} \fontname{Euclid}r');
  ylabel('\bf \fontsize{15} \fontname{Euclid} w(r)');
  legend('\fontname{Euclid} Potencijal srednje sile')

  fprintf ( 1, '\n' );
  fprintf('-----------------------------------------------------------------------------------------------\n')
  fprintf ( 1, '\n' );
  fprintf ( 1, 'Srednja potencijalna energija po cestici izracunata koristeci radijalnu funk. raspodele je: %f r.j. = %f kJ/mol \n',trapz(rr,8*pi*ro*g_rad_final.*rr.^2.*(rr.^(-12)-rr.^(-6))),trapz(rr,8*pi*ro*g_rad_final.*rr.^2.*(rr.^(-12)-rr.^(-6)))*0.997);
  fprintf ( 1, 'Srednji pritisak izracunat koristeci radijalnu funk. raspodele je:                           %f r.j. = %f MPa \n',Av_temperature*ro+trapz(rr,32*pi*(ro^2)*g_rad_final.*(rr.^(-10)-0.5*rr.^(-4))), Av_temperature*ro+trapz(rr,32*pi*(ro^2)*g_rad_final.*(rr.^(-10)-0.5*rr.^(-4)))*42);
  fprintf ( 1, '\n' );
  fprintf('-----------------------------------------------------------------------------------------------\n')
  fprintf ( 1, '\n' );  
% %%% raspodela x-komponente brzine atoma
 figure
 vxbins=-4:0.1:4;
 plot(vxbins,sum(histc(Velocities(1,:,Nsteps-Nsample+1:Nsteps),vxbins),3)/(Nsample*NParticles),'k-',...
    'LineWidth',1.5);
 xlabel('\bf \fontsize{15} \fontname{Euclid} v_{x}');
 ylabel('\bf \fontsize{15} \fontname{Euclid} f(v_{x})');
 legend('\fontname{Euclid} Raspodela x-komponente brzine atoma')

% %%% deo za racunanje autokorelacione funkcije brzine i koeficijenta difuzije
% maxt=750;  % maksimalan broj koraka za autokorelacionu funkciju 
% maxt0=5000;  % maksimalan broj koraka za pocetke uzorkovanja
% freqt0=10; % frekvencija na kojoj se novi pocetak bira
% t0brojac=0; % brojac pocetaka za uzorkovanje autokorelacione funkcije
% t0vreme=zeros(1,maxt0); % niz koji sadrzi vremena pocetka uzorkovanja
% vt0=zeros(3,NParticles,maxt0); % brzine cestica u pocecima uzorkovanja
% r2t0=zeros(3,NParticles,maxt0); % polozaji cestica u pocecima uzorkovanja 
% vacf=zeros(3,NParticles,maxt); % autokorelaciona funkcija brzine
% r2t=zeros(3,NParticles,maxt); % razlike polozaja cestica 
% brojacuzoraka=zeros(1,maxt); % broj uzorkovanja u autokorelacionj funkciji
% for vreme=poc_uzorkovanje:Nsteps;
%   if (mod(vreme,freqt0)==0)
%     t0indeks=mod(t0brojac,maxt0)+1; % indeks trenutnog pocetka t0
%     t0vreme(t0indeks)=vreme;
%     t0brojac=t0brojac+1;
%     vt0(:,:,t0indeks)=Velocities(:,:,vreme);
%     r2t0(:,:,t0indeks)=coords(:,:,vreme); 
%   end
%   for j=1:min(t0brojac,maxt0)
%     korelaciono_vreme=vreme-t0vreme(j)+1;
%     if (korelaciono_vreme<=maxt)
%         brojacuzoraka(korelaciono_vreme)=brojacuzoraka(korelaciono_vreme)+1;
%         vacf(:,:,korelaciono_vreme)=vacf(:,:,korelaciono_vreme)+Velocities(:,:,vreme).*vt0(:,:,j);
%         r2t(:,:,korelaciono_vreme)=r2t(:,:,korelaciono_vreme)+(coords(:,:,vreme)-r2t0(:,:,j)).^2; 
%     end
%   end
% end
% autokor=sum(sum(vacf,1),2)/NParticles; % suma po komponentama i cesticama
% autokor_r2=sum(sum(r2t,1),2)/NParticles; % suma po komponentama i cesticama
% z = size(autokor);
% autokorelaciona=reshape(autokor,[z(2:end) 1]);
% autokorelaciona_r2=reshape(autokor_r2,[z(2:end) 1]);
% autokorelaciona=autokorelaciona./brojacuzoraka;
% autokorelaciona_r2=autokorelaciona_r2./brojacuzoraka;
% figure
% plot(dt*(0:maxt-1),autokorelaciona,dt*(0:maxt-1),0,'-b')
% xlabel('Vreme');
% ylabel('<v(0)v(t)>');
% legend('Autokorelaciona funkcija brzine')
% figure
% plot(dt*(0:maxt-1),autokorelaciona_r2,dt*(0:maxt-1),0,'-b')
% xlabel('Vreme');
% ylabel('<(r(0)-r(t))^2>');
% legend('Srednji kvadratni pomeraj atoma')
% koef_difuzije=trapz(dt*(0:maxt-1),autokorelaciona)/3;
% koef_difuzije_r2=(autokorelaciona_r2(length(autokorelaciona_r2))-autokorelaciona_r2(length(autokorelaciona_r2)/2))/(6*dt*(maxt-1-(maxt-1)/2));
% fprintf ( 1, 'Koeficient difuzije izracunat koriscenjem autokorelacione funkcije brzine: %e= %e cm2/s \n', koef_difuzije,koef_difuzije*(0.3405*1e-7)^2/(2170*1e-15));
% fprintf ( 1, 'Koeficient difuzije izracunat koriscenjem srednjeg kvadratnog pomeraja atoma: %e= %e cm2/s \n', koef_difuzije_r2,koef_difuzije_r2*(0.3405*1e-7)^2/(2170*1e-15));

% %% vratiti cestice u kutiju:
coordsPBC=zeros(3,NParticles,Nsteps+1);
for frame=1:Nsteps+1
coordsPBC(:,:,frame)=coords(:,:,frame)-Lbox*floor(coords(:,:,frame)./Lbox);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Deo za pravljenje filma simulacije (velicina atoma je smanjena da bi se bolje videli; u simulaciji je precnik jednak 1) 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

writerObj = VideoWriter('MD_NVE.avi');
open(writerObj);
figure
for j = poc_uzorkovanje:20:Nsteps
    scatter3(coordsPBC(1,:,j),coordsPBC(2,:,j),coordsPBC(3,:,j),'MarkerEdgeColor','k','MarkerFaceColor',[0 .75 .75])
    axis([0,Lbox,0,Lbox,0,Lbox]); frame = getframe;
     writeVideo(writerObj,frame);
 end
 close(writerObj);
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

toc