% Ising Monte Carlo simulation 
% Version 03.09.2009
% 
% Example:
% ising(1.0,20,100,30,0.1,0.1,50,0,0,1,1,0,1)
% arguments            value in the example    meaning
% ---------            --------------------    -------
%       J              1.0                     coupling constant  
%       L              20                      lattice size (L x L)
%       MCSTEPS        100                     # of MC steps
%       THERM          30                      # of thermalization steps
%       T              0.1                     first temperature
%       DT             0.1                     temperature step
%       Tnum           50                      # of temperature values
%       B              0                       first magnetic field value
%       DB             0                       magnetic field step 
%       Bnum           1                       # of magnetic field values
%       plotfreq       1                       plotting frequency
%       pauses         0                       pause between plots
%       init           1                       lattice initialization
% 
% If you have problems, ask Vesa Apaja, room YN212 (Nanoscience Center 2nd floor)


function[]=ising(J,L,MCSTEPS,THERM,T,DT,Tnum,B,DB,Bnum,plotfreq,pauses,init)

    disp('Usage: example   ising(1.0,20,100,50,0.1,0.2,25,0.0,0,1,1,0,1)')
    disp(' ')
    disp('List of parameters:')
    disp('1.  J         Coupling constant (e.g.. 1.0 or 0.0)')
    disp('2.  L         Square lattice size (e.g. 10 means 10x10 lattice)')
    disp('3.  MCSTEPS   # of MC steps (e.g. 100)')
    disp(['4.  THERM     # of thermalisation MC steps (e.g. 50 means no ' ...
          'measurements during first 50 MC steps)'])
    disp('5.  T         First temperature value (e.g. 0.1)')
    disp('6.  DT        Temperature step (e.g. 0.2)')
    disp('7.  Tnum      # of temperatures (e.g. 25)')
    disp('8.  B         First magnetic field value (e.g. 0.0)')
    disp('9.  DB        Magnetic field step (e.g. 0.1)')
    disp('10. Bnum      # of magnetic field values (e.g. 3)')
    disp(['11. plotfreq  Plotting frequency for lattice snaphots (0=never, ' ...
          '1=after each MC step (one sweep), 2=all the time, after ' ...
          'every spin flip)'])
    disp(['12. pause     pause at snapshots (0 or 1 for no or yes; ' ...
          'this lets you save the plot is you wish)'])
    disp('13. init      Initial lattice:')
    disp('                  0=random')
    disp('                  1=spins up, -1=spins down')
    disp('                  2=use previous (start spins up)')
    disp('                 -2=use previous (start spins down)')
    disp('                  3=use previous (start with random spins)')
    %    disp(['                  4=use previous (read first from file ' ...
    %      'lattice.dat)'])
    disp(['                  5=use previous (start with half spins ' ...
          'up, half down)'])
    disp(' ')
    
    disp('You gave parameters:')
    disp(['       J = ',num2str(J)])
    disp(['       L = ',num2str(L)])
    disp([' MCSTEPS = ',num2str(MCSTEPS)])
    disp(['   THERM = ',num2str(THERM)])
    disp(['       T = ',num2str(T)])
    disp(['      DT = ',num2str(DT)])    
    disp(['    Tnum = ',num2str(Tnum)])
    disp(['       B = ',num2str(B)])
    disp(['      DB = ',num2str(DB)])
    disp(['    Bnum = ',num2str(Bnum)])
    disp(['plotfreq = ',num2str(plotfreq)])
    disp(['  pauses = ',num2str(pauses)])
    disp(['    init = ',num2str(init)])
    disp('')
    disp('press a key to continue (or control-C to abort)')
    disp('Info:')
    disp(['   - This program will make plots according to the parameter ' ...
          'plotfreq'])
    disp('   - Some plots may be utterly useless')
    disp(['   - Not all bad parameter choices are trapped in the ' ...
          'beginning, the run may crash'])
    disp(['   - All results except the lattice spins will be saved to an ascii data file ' ...
          'results.dat,']) 
    disp('     which you can use in any plotting program you wish.')
    disp(['   - You may want to backup the the old results.dat file ' ...
          'now, it will be overwritten'])
    pause
    PARAMETERS=[J,L,MCSTEPS,THERM,T,DT,Tnum,B,DB,Bnum,plotfreq,pauses,init];
       
    
   isingJxjY(PARAMETERS);
        
       function[]=isingJxjY(PARAMETERS)
 
 J=PARAMETERS(1);
 L=PARAMETERS(2);
 NINDEP=1 ; % keep fixed
 MCSTEPS=PARAMETERS(3);
 THERM =PARAMETERS(4);
 T=PARAMETERS(5);
 DT=PARAMETERS(6);
 Tnum=PARAMETERS(7); 
 B=PARAMETERS(8);
 DB=PARAMETERS(9);
 Bnum=PARAMETERS(10); 
 plotfreq=PARAMETERS(11);
 pauses=PARAMETERS(12);
 init=PARAMETERS(13);
 %
 % always square lattice
 %
 Lx = L;
 Ly = L;
 %
 % marker size in spy plotting
 %
  if(L<=10)
     msize = 50;
 elseif(L<=20)
     msize = 40;
 elseif(L<40)
     msize = 20;
 else
     msize = 8;
 end
 %
 MN=Bnum;
 MN1=Tnum;
 %
 % spin values in the simulation
 %
 ispin(1)=-1         ;
 ispin(2)=1          ;
 mu = 2 ;
 %
 % sanity tests
 %
 if (THERM>= MCSTEPS)
     disp (' ')
     disp ('Bad parameters')
     disp ('NTHERM >= MCSTEPS : no MC steps to measure')
     disp ('Either increase MCSTEPS or decrease NTHERM')
     return     
 end

 if (init~=0 & init~=1 & init~=-1& init~=2 & init~=-2 & init~=3 ...
     & init~=4 & init~=5 )
     disp('bad init lattice, choose 0, 1, -1, 2, -2, 3, 4 or 5')
     return
 end    
 if(plotfreq ~=0 & plotfreq ~=1 & plotfreq ~=2)
     plotfreq
     disp('bad value for plotfreq')
     return
 end
 if(pauses ~=0 & pauses ~=1)
     disp('bad value for pauses; assuming no pauses')
     pauses = 0;
 end   
 
 
 

 %
 % just a warning
 %
 if(L>100)
     disp ('L>100: Big lattice, this simulation may take a looooong time.')
 end
 
 
 disp('========== INPUT PARAMETERS ==============================================')
 if (abs(J)<0.0001) 
     disp('J=0 : This is a system of IDEAL PARAMAGNETS')
     disp('Program uses parameter x=mu B/(kB T) and keeps T=1 fixed')
     Tnum = 1;
     T = 1;
     dT = 0;  
     %B = 0;
     %DB = 0.2*mu;
     %Bnum = 16;
     if(Bnum<=1) 
         disp('WARNING: Just one x value (taken from B table) will not show much.')
     end
 end
 if (J>0)
     disp(['coupling constant J = ' num2str(J)]);
     disp('J>0: This is a system of FERROMAGNETS')
 end
 if (J<0)
     disp(['coupling constant J = ' num2str(J)]);
     disp('J<0: This is a system of ANTIFERROMAGNETS')
     %disp('fixing B=0')
     %Bnum = 1;
     %B = 0;
     %dP = 0;
 end
     
 disp('Lattice dimensions')
 disp([int2str(Lx) ' x ' int2str(Ly)])
 %
 % tabulate T and B values
 %
 if (Bnum>0)
     for i=1:Bnum
         B1(i)=B+(i-1)*DB;
     end
     B=B1;
 end
if (Tnum>0)
    for i=1:Tnum
        T1(i)=T+(i-1)*DT;
    end
    T=T1;
end

if(abs(J)<0.0001)
    disp('Table of x values') 
    % multiply B with mu
    B = mu*B
    disp(B)
else
    disp('Table of magnetic field values') 
    B
    disp('Table of temperature values') 
    T
end    

if(init==0) 
    disp('init=0: setting spins to random values +/- 1 between MC runs')
end 
if(init==1)
    disp('init=1: resetting spins up between MC runs')
end 
if(init==-1)
    disp('init=-1: resetting spins down between MC runs')
end 
if(init==2)
    disp(['init=2: keeping previous spin configuration between MC ' ...
          'runs (the first one is spins up)'])
end
if(init==-2)
    disp(['init=-2: keeping previous spin configuration between MC ' ...
          'runs (the first one is spins down)'])
end
if(init==3)
    disp(['init=3: keeping previous spin configuration between MC ' ...
          'runs (the first one is random spins)'])
end
%if(init==4)
%    disp(['init=4, trying to read spin configuration from file ' ...
%          'lattice.dat'])
%end
if(init==5)
    disp(['init=5, keeping previous spin configuration between MC ' ...
          'runs (the first one is half up, half down)'])
end


 % set random generator to Mersenne twister and use clock seed
 seed = sum(100*clock);
%  rand('twister',seed)  % version 7 Matlab  only
rand('seed',seed)
 
 
 if(plotfreq == 0)
     disp('plotfreq = 0, not making plots during MC runs')
 elseif(plotfreq==1)
     disp(['making plots after every MC step:  this slows down the ' ...
           'computation but may be useful'])
 elseif(plotfreq==2)
     disp(['making plots after spin flip:  this slows down the ' ...
           'computation but may be useful'])
 else
     disp('Bad plotfreq (not 0,1)')
     return
 end


disp('========== END OF INPUT PARAMETERS =======================================')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% translation vectors for periodic boundary conditions 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for a1=1:Lx
   xright(a1)=a1+1;
   xleft(a1)=a1-1;
end 
xright(Lx)=1;
xleft(1)=Lx;
for a1=1:Ly
   yup(a1)=a1+1;
   ydown(a1)=a1-1;
end 
yup(Ly)=1   ;
ydown(1)=Ly ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nnsum(1)=-2         ;
nnsum(2)=0          ;
nnsum(3)=2          ;
nnsum(4)=-2         ;
nnsum(5)=0          ;
nnsum(6)=2          ;


index=1;
indexT=1;
MCSTEP=1:MCSTEPS;
disp('========== Monte Carlo begins ============================================')
%
% start lattice
%
if(init==1 | init==2) 
    lattice1=ones(Lx,Ly) ; 
    disp('first lattice is spins up')
elseif(init==0 | init==3)
    lattice1=2*round(rand(Lx,Ly))-1; 
    disp('first lattice has random spins')
elseif(init==-1 | init==-2)
     lattice1=-ones(Lx,Ly) ;
     disp('first lattice is spins down')
elseif(init==4)
    disp('first lattice from file')
    eval(['load lattice']);
elseif(init==5)
    disp('first lattice is half up half down')
    lattice1=ones(Lx,Ly) ; 
    lattice1(Lx/2:Lx,:) = -lattice1(Lx/2:Lx,:);
end
%
% Initial spin config
%
figure(1);clf(1);
%spy(lattice1 + ones(size(lattice1)),'r',msize); hold on
%spy(lattice1 - ones(size(lattice1)),'b',msize);
colormap bone
image(-lattice1*128);
line2=['J=',num2str(J)];
title({'snapshot of the spin configuration';line2})
if(pauses) 
    disp('press a key to continue (or use the figure)')
    pause
end

   




if(abs(J)<0.0001);disp('here x=mu B/(kB T)');end
for it=1:Tnum % T loop
   tic
   for ib=1:Bnum  % B loop
        if(abs(J)<0.0001)
            tic
            disp(['x = ',num2str(B(ib)/mu)])
        else            
            disp(['T, B = ' num2str(T(it)) ' ' num2str(B(ib))])
        end
        M=zeros(1,MCSTEPS)              ;
        ABS_M=zeros(1,MCSTEPS)          ;
        E=zeros(1,MCSTEPS)              ;
        BperkT = B(ib)/T(it)    ;
        JxperkT = J/T(it)    ;
        JyperkT = J/T(it)    ;
        line2=['J=',num2str(J),'   T=',num2str(T(it)),'   B=',num2str(B(ib))];
        % 
        % energy from nearest neighbour spin sums
        %
        q_apu=1;
        for q1 =1:3	 % x direction
            for q2=1:3     % y direction
                for q3 =1:2  % spin values up and down
                    w(q_apu,q3)=1;
                    dEperkT(q_apu,q3)=2*JxperkT*nnsum(q1)*ispin(q3)+BperkT*ispin(q3)+2*JyperkT*nnsum(q2)*ispin(q3) ;
                    if (dEperkT(q_apu,q3) > 0) 
                        w(q_apu,q3) = exp(-dEperkT(q_apu,q3));
                    end
                end
                q_apu=1+q_apu;
            end 
        end
        %
        % initial lattice: random, all spins up (1) or keep previous
        %
        if(init==0) 
            lattice1=2*round(rand(Lx,Ly))-1; 
            disp('init lattice to random spins')
        end 
        if(init==1)
            lattice1=ones(Lx,Ly); 
            disp('init lattice to spin ups')
        end
        if(abs(init)>=2)
            % init : do nothing here
            disp('start from previous spin lattice')
        end
        for k=1:NINDEP   % Loop over independent MC runs
	    lattice=lattice1;
            for mcs=1:MCSTEPS   % Loop over MC steps       
                for kk=1:(Lx*Ly) % sweep over the lattice 
                    i = floor(Lx*rand+1)    ;
                    j = floor(Ly*rand+1)     ;
                    spinX=lattice(xleft(i),j)+lattice(xright(i),j);
                    spinY=lattice(i,yup(j))+lattice(i,ydown(j));
                    ijspin=lattice(i,j)  ;
                    
                    if (ijspin==-1)
                        spin=1;
                    else
                        spin=2;
                    end 
                    if (spinX==-2)
                        if (spinY==-2)
                            nnsumma=1;
                        end
                        if (spinY==0)
                            nnsumma=2;
                        end
                        if (spinY==2)
                            nnsumma=3;
                        end
                    end
                    if (spinX==0)
                        if (spinY==-2)
                            nnsumma=4;
                        end
                        if (spinY==0)
                            nnsumma=5;
                        end
                        if (spinY==2)
                            nnsumma=6;
                        end
                    end    
                    if (spinX==2)
                        if (spinY==-2)
                            nnsumma=7;
                        end
                        if (spinY==0)
                            nnsumma=8;
                        end
                        if (spinY==2)
                            nnsumma=9;
                        end
                    end 
                    if (w(nnsumma,spin) > rand)
                        lattice(i,j) = -1*ijspin;
                    end 
                    
                end  % sweep ends
                
                apu=1;
                for i1=1:Ly
                    for i2=1:Lx
                        ENERGY_Y(apu)=lattice(i2,yup(i1))*lattice(i2,i1);
                        ENERGY_X(apu)=lattice(xleft(i2),i1)*lattice(i2,i1);
                        apu=apu+1;
                    end 
                end
                apu1=0;
                apu2=0;
	        apuY=sum(ENERGY_Y);
	        apuX=sum(ENERGY_X);
       	        apu2=sum(sum(lattice));
       		ENERGY(mcs)=-J*apuX-J*apuY-B(ib)*apu2;
	       	Magnetization(mcs)=apu2;
	      	Magnetization1(mcs)=abs(apu2);  
            
            % snaphot between MC steps	    
            if(plotfreq==2)
                figure(1);clf(1);
                %spy(lattice + ones(size(lattice)));
                %spy(lattice + ones(size(lattice)),'r',msize); hold on
                %spy(lattice - ones(size(lattice)),'b',msize);hold on
		colormap bone;
		image(-lattice*128);
                title({'snapshot of the spin configuration';line2})
                drawnow;
		%pause(0.1)
                if(pauses) 
                    disp('press a key to continue (or use the figure)')
                    pause
                end
            end             
            end % loop over MC steps ends
            
	    M=Magnetization+M            ;
	    ABS_M=Magnetization1+ABS_M   ;
	    E=E+ENERGY                 ;
	end    % Loop over independent MC runs ends
		


        if(plotfreq)
            %
            % spin config
            %
            figure(1);clf(1);
            %spy(lattice + ones(size(lattice)));
            %spy(lattice + ones(size(lattice)),'r',msize); hold on
            %spy(lattice - ones(size(lattice)),'b',msize);
	    colormap bone;
	    image(-lattice*128);
            title({'snapshot of the spin configuration';line2})
            drawnow;pause(0.1)
            if(pauses) 
                disp('press a key to continue (or use the figure)')
                pause
            end
            %
            % Plot evolution of |M|
            %
            if( 1==0)   
                figure(5);clf(5);
                if(J>0)
                  plot([THERM THERM],[-.6 1.2],'r');hold on;
                else
                  plot([THERM THERM],[-.1 0.1],'r');hold on;
                end
            
                plot(MCSTEP,ABS_M./(Lx*Ly*k),'g'); hold on
                title(['Stochastic data (the red line shows where thermalization ' ...
                      'ends)'])
                xlabel('MC steps');ylabel('|M|')
                %
                % Plot evolution of M
                %
                figure(6); clf(6)
                if(J>0) 
                    plot([THERM THERM],[-.6 1.2],'r');hold on; 
                    plot(MCSTEP,M./(Lx*Ly*k),'g'); hold on
                else
                plot([THERM THERM],[-.1 0.1],'r');hold on; 
                    plot(MCSTEP,M./(Lx*Ly*k),'g'); hold on
                end 
                title(['Stochastic data (the red line shows where thermalization ends)'])
                xlabel('MC steps');ylabel('M')
            end
        end
        %
        % snapshot of spin configuration
        %
        if(plotfreq>0) 
            figure(1);clf(1);
            %spy(lattice + ones(size(lattice)),'r',msize); hold on
            %spy(lattice - ones(size(lattice)),'b',msize);
	    colormap bone;		 
	    image(-lattice*128);
            title({'snapshot of the spin configuration';line2})
            drawnow;pause(0.1)
            if(pauses) 
                disp('press a key to continue (or use the figure)')
                pause
            end
        end
        norm = Lx*Ly*NINDEP   ;
        M=M./norm           ;
        E=E./norm            ;
        E2=E.^2 ;
        ABS_M=ABS_M./norm    ;
        disp(['<M>=',num2str(mean(M))])
        %
        % save results for later analysis
        %        
        eval(['save tulokset' int2str(index) ,' ABS_M M E E2']);
        eval(['save parametrit' int2str(index) ,[' J L NINDEP MCSTEPS ' ...
                            'it T ib B mu']]);
        index=index+1;
        %
        % timing info for J=0 
        %        
        if(abs(J)<0.0001)
         time = toc;
         time = time*(Bnum-ib);
         if(time>60)
            disp(['Ready in ',num2str(time/60.) ' minutes'])
         else
            disp(['Ready in ',num2str(time) ' seconds'])
         end
        end
        %
        % store current lattice for a run with a new B
        %
        if(abs(init)>1) 
            lattice1 = lattice;
        end
        %
   end % Loop over B ends
   %
   % store current lattice for a run with a new T
   %   
   lattice1 = lattice;  
   %    
   % save lattice to file "lattice"
   %
   if( init == 4) 
       eval(['save lattice',' lattice1']);
   end
   %
   % timing info for J not 0 and plotfreq=0 
   %  
   if(abs(J)>0.0001)
     time = toc;
     time = time*(Tnum-it);
     if(time>60)
       disp(['Ready in ',num2str(time/60.) ' minutes'])
     else
         disp(['Ready in ',num2str(time) ' seconds'])
     end
   end
   
end % Loop over T ends

disp('========== Monte Carlo ends ==============================================')
disp('analysis and plotting the data follows')
analyse(index-1,THERM)




