% Program "MCargon_NVT"

clc; close all; clear; tic

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

NParticles=120;                % Broj cestica u simulacionoj kutiji
T=1.3;                         % Temperatura
ro=0.8;                        % Cesticna gustina = NParticles/Lbox^3
Nsample=1000000;                   % Broj koraka za racunanje srednje vrednosti
Nsteps=2000000;                    % Broj koraka u simulaciji
freq_uzorak=100;               % Frekvencija uzorkovanja osobina koje se racunaju
MaxDisp=0.3;                   % Maksimalni pomeraj atoma

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

Lbox=(NParticles/ro)^(1/3); 
Rcutoff=0.4999*Lbox;  
Ecutoff=4*(1/(Rcutoff^12)-1/(Rcutoff^6));            % PotEnergija za rastojanje 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
if (Rcutoff<1)
    error('upozorenje:Rcutoff<1')
end
if (Nsteps<Nsample)
    error('upozorenje:Nsteps<Nsample')
end

Neq=Nsteps-Nsample;

fprintf('\n')
fprintf('-----------------------------------------------------------------------------------------------\n')
fprintf('\n')
fprintf ( 1, 'Program za proracun termodinamickih osobina argona pomocu Monte Karlo metode u NVT ansamblu. \n' );
fprintf ( 1, 'Program koristi redukovane jedinice. \n' );
fprintf ( 1, 'Potencijalna energija se racuna koristeci ograniceni LJ potencijal. \n' );
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 = %f kg/m3 \n', ro, 1688*ro );
fprintf ( 1, 'Temperatura je:                                 %f = %f K \n', T,T*119.91 );
fprintf ( 1, 'Ukupni broj koraka u simulaciji je:             %d \n', Nsteps );
fprintf ( 1, 'Broj koraka za uravnotezenje sistema je:        %d \n', Neq );
fprintf ( 1, 'Broj koraka za uzorkovanje je:                  %d \n', Nsample );
fprintf ( 1, 'Duzina maksimalnog pomeraja je:                 %f = %f %s \n', MaxDisp, 0.3405*MaxDisp/10, char(197) );
fprintf ( 1, 'Duzina sumulacione kutije je:                   %f = %f %s \n', Lbox, 0.3405*Lbox/10, char(197) );
fprintf ( 1, 'Rastojanje iznad kojeg se ne racuna sila je:    %f = %f %s \n', Rcutoff, 0.3405*Rcutoff/10, char(197) );
fprintf ( 1, 'Potencijalna energija na tom rastojanju je:     %f \n', Ecutoff );
fprintf ( 1, 'Korekcija za potencijalnu energiju po atomu je: %f \n', Epot_tail );
fprintf ( 1, 'Korekcija za pritisak je:                       %f \n', P_tail );
fprintf('\n')
fprintf ( 1, 'Termodinamicke osobine se racunaju na svakih %d koraka. \n', freq_uzorak );

coords=zeros(3,NParticles);

% Postavljanje atoma u 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) = dx*ix;
    coords(2,ncount) = dx*iy;
    coords(3,ncount) = dx*iz;
   end
  end
 end
end

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

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

prihvacen=0; 
vir=zeros(1,(Nsteps-poc_uzorkovanje+1)/freq_uzorak);
ene=zeros(1,(Nsteps-poc_uzorkovanje+1)/freq_uzorak);
brojac_u=1;

for m=1:Nsteps  
    %% Izabrati atom koju treba pomeriti
    x=floor(rand*NParticles)+1;
    %% Izracunati staru potencijalnu energiju za atom koji se pomera
    PotEnergija_stara=PotEne1(coords,x,NParticles,Rcutoff,Lbox);
    %% Izabrati nove koordinate
    nove_x=(rand(3,1)-0.5)*MaxDisp;
    coords_nove=coords;
    coords_nove(:,x)=coords_nove(:,x)+nove_x;
    %% Izracunati novu potencijalnu energiju za atom koji se pomera
    PotEnergija_nova=PotEne1(coords_nove,x,NParticles,Rcutoff,Lbox);
    if (rand<exp(-(PotEnergija_nova-PotEnergija_stara)/T))
      %% Povecati brojac prihvacenih koraka
        coords=coords_nove;
        coords(:,x)=coords(:,x)+Lbox*(coords(:,x)<0)-Lbox*(coords(:,x)>Lbox);
        prihvacen=prihvacen+1; 
    end
    if ((m>=poc_uzorkovanje)&&(mod(m,freq_uzorak)==0))
       [PotE,Vir]=Osobine(coords,NParticles,Rcutoff,Lbox);
       ene(brojac_u)=PotE/NParticles+Epot_tail;
       vir(brojac_u)=Vir;
       brojac_u=brojac_u+1;
    end
end

Epot_poatomu=mean(ene);
Av_pressure=ro*T+mean(vir)/(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 atomu je:       %f = %f kJ/mol \n', Epot_poatomu+1.5*T,(Epot_poatomu+1.5*T)*0.997);
fprintf ( 1, 'Srednja korigovana potencijalna energija po atomu je: %f = %f kJ/mol \n', Epot_poatomu, Epot_poatomu*0.997);
fprintf ( 1, 'Srednji korigovani pritisak je:                       %f = %f MPa \n', Av_pressure+P_tail,(Av_pressure+P_tail)*42);
fprintf ( 1, 'Procenat prihvacenih pomeraja:                        %f %%\n', prihvacen*100/Nsteps);
fprintf ( 1, '\n' );
fprintf('-----------------------------------------------------------------------------------------------\n')
fprintf ( 1, '\n' );

figure
plot((0:length(ene)-1)*freq_uzorak,ene,'k-',...
    'LineWidth',1.5)
xlabel('\bf \fontsize{15} \fontname{Euclid} Vreme u MC koracima');
ylabel('\bf \fontsize{15} \fontname{Euclid} E_{pot}');
legend('\fontname{Euclid} Trenutna vrednost pot. ener. u funkciji MC koraka')


trenutna_sr_potE=zeros(1,length(ene));
for i=1:length(ene) 
trenutna_sr_potE(i)=sum(ene(1:i))/i; 
end

figure
plot((0:length(ene)-1)*freq_uzorak,trenutna_sr_potE,'k-',...
    'LineWidth',1.5)
xlabel('\bf \fontsize{15} \fontname{Euclid} Vreme u MC koracima');
ylabel('\bf \fontsize{15} \fontname{Euclid} < E_{pot} >');
legend('\fontname{Euclid} Trenutna vrednost srednje pot. ener. u funkciji MC koraka')

%%% Proracun autokorelacione funkcije potencijalne energije atoma

% maxt=5000;              % Maksimalan broj koraka za autokorelacionu funkciju 
% maxt0=50000;            % 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
% enet0=zeros(1,maxt0);   % Enerija u pocecima uzorkovanja
% Eacf=zeros(1,maxt);     % Autokorelaciona funkcija
% brojacuzoraka=zeros(1,maxt); % Broj uzorkovanja u autokorelacionj funkciji
% for vreme=1:length(ene)
%   if (mod(vreme,freqt0)==0)
%     t0indeks=mod(t0brojac,maxt0)+1; % Indeks trenutnog pocetka t0
%     t0vreme(t0indeks)=vreme;
%     t0brojac=t0brojac+1;
%     enet0(t0indeks)=ene(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;
%         Eacf(korelaciono_vreme)=Eacf(korelaciono_vreme)+(ene(vreme)-Epot_poatomu)*(enet0(j)-Epot_poatomu);
%     end
%   end
% end
% autokor=Eacf./brojacuzoraka;
% autokor=autokor./autokor(1);
% figure
% plot((0:maxt-1)*freq_uzorak,autokor,(0:maxt-1)*freq_uzorak,0,'-b')
% xlabel('Vreme u MC koracima');
% ylabel('R');
% legend('Normirana autokorelaciona funkcija potencijalne energije po atomu')
% tau=sum(autokor)*freq_uzorak; %% Autokorelaciono vreme
% fprintf ( 1, 'Statisticki nezavisna konfiguracija se javlja posle %d MC pomeraja.\n', round(2*tau));
% fprintf ( 1, 'Srednja korigovana potencijalna energija po atomu je: %f %s %f.\n', Epot_poatomu,char(177),sqrt(var(ene)/length(ene)*2*tau/freq_uzorak));
% fprintf ( 1, 'Toplotni kapacitet po atomu je: %f \n', 1.5+var(ene)*2*tau/freq_uzorak/T^2);

toc