this simulation stops working after the end moves have been done and the corner moves do not work

clc;

close;

int_polymer=[10 10; 9 10; 8 10; 7 10; 7 9; 8 9; 9 9; 10 9; 10 8; 10 7; 9 7; 9 8; 8 8; 8 7; 7 7; 7 8];

new_polymer =zeros(16,2);

% Create video writer object

writer = VideoWriter('polymer_simulation.mp4','MPEG-4');

writer.FrameRate = 10; % Set the frame rate of the video

open(writer);

% Plot initial state of the polymer

plot(int_polymer(:,1), int_polymer(:,2), 'bo-', 'MarkerSize', 3, 'LineWidth', 2);

title(['interaction energy=',num2str(-1.25*9),' no of steps=',num2str(0),' no of cycles=',num2str(0)]);

grid on;

axis([3 14 4 13]);

set(findobj('type','line'),'XData',int_polymer(:,1),'YData',int_polymer(:,2));

set(findobj('type','scatter'),'XData',int_polymer(:,1),'YData',int_polymer(:,2));

int_ener=-1.25;

Ei=9*int_ener;

no_sim=100000;

kt=1;

for i=1:50

for j=1:no_sim

tpoly=temp_poly(int_polymer);

Et=temp_ener(tpoly)*int_ener;

a=mc_sim(Et,Ei);

if a==1

int_polymer=tpoly;

Ei=Et*int_ener;

end

% Plot current state of the polymer

plot(int_polymer(:,1), int_polymer(:,2), 'bo-', 'MarkerSize', 3, 'LineWidth', 2);

title(['interaction energy=',num2str(Ei),' no of steps=',num2str(j),' no of cycles=',num2str(i)]);

grid on;

axis([3 14 4 13]);

set(findobj('type','line'),'XData',int_polymer(:,1),'YData',int_polymer(:,2));

set(findobj('type','scatter'),'XData',int_polymer(:,1),'YData',int_polymer(:,2));

% Write the current frame to the video file

frame = getframe(gcf);

writeVideo(writer,frame); % Add the current frame to the video file

drawnow;

end

end

close(writer); % Close the video writer object

function tpoly=temp_poly(polymer)

tpoly=polymer;

idx = randi(size(polymer,1));

mono_pos=polymer(idx,:);

allow_choice=0;

if idx==1

choices=[mono_pos(1)+1 mono_pos(2)-1;

mono_pos(1)+1 mono_pos(2)+1;

mono_pos(1)-1 mono_pos(2)+1;

mono_pos(1)-1 mono_pos(2)-1];

for i=1:4

if ((choices(i,1)-polymer(idx+1,1))^2+(choices(i,2)-polymer(idx+1,2))^2)==1

a=ismember(choices(i,:),polymer);

if any(~a)

allow_choice=choices(i,:);

end

end

end

elseif idx==16

choices=[mono_pos(1)+1 mono_pos(2)-1;

mono_pos(1)+1 mono_pos(2)+1;

mono_pos(1)-1 mono_pos(2)+1;

mono_pos(1)-1 mono_pos(2)-1];

for i=1:4

if ((choices(i,1)-polymer(idx-1,1))^2+(choices(i,2)-polymer(idx-1,2))^2)==1

a=ismember(choices(i,:),polymer);

if any(~a)

allow_choice=choices(i,:);

end

end

end

else

choices=[mono_pos(1)+1 mono_pos(2)-1;

mono_pos(1)+1 mono_pos(2)+1;

mono_pos(1)-1 mono_pos(2)+1;

mono_pos(1)-1 mono_pos(2)-1];

for i=1:4

if ((choices(i,1)-polymer(idx+1,1))^2+(choices(i,2)-polymer(idx+1,2))^2)==1&&(choices(i,1)-polymer(idx-1,1))^2+((choices(i,2)-polymer(idx-1,2))^2)==1

a=ismember(choices(i,:),polymer);

if any(~a)

allow_choice=choices(i,:);

end

end

end

end

if allow_choice~=0

tpoly(idx,:)=allow_choice;

end

end

function noofinter=temp_ener(polymer)

noofinter=-15;

distances = pdist(polymer);

for i=1:120

if (distances(i)==1)

noofinter=noofinter+1;

end

end

end

function result=mc_sim(Et,Ei)

del_E = Et - Ei;

r=rand;

w = exp(-del_E);

if w >= 1

result=1;

end

if w < 1

if w > r

result=1;

else

result=0;

end

end

end

Answer 1

[Rasmita]

Hi,

As per my understanding, the issue with the provided code is that the algorithm is getting stuck at the end moves of the polymer and the corner moves do not work. This is because there are fewer options to change the position of the monomer. One solution for this issue could be to modify the code to allow the corner moves to take place as well.

Currently, the ‘temp_poly’ function only allows for end moves to be made. To enable corner moves, you could modify this function to consider corner monomers as well. For example, you could add the following code to allow corner moves:

if idx==1 || idx==4 || idx==13 || idx==16 % check if corner monomer 
    choices=[mono_pos(1)+1 mono_pos(2)+1;             mono_pos(1)+1 mono_pos(2)-1;             mono_pos(1)-1 mono_pos(2)+1;             mono_pos(1)-1 mono_pos(2)-1]; 
    for i=1:4 
        if ((choices(i,1)-polymer(idx+1,1))^2+(choices(i,2)-polymer(idx+1,2))^2)==1 && ((choices(i,1)-polymer(idx-1,1))^2+(choices(i,2)-polymer(idx-1,2))^2)==1 
            a=ismember(choices(i,:),polymer); 
            if any(~a)  
                allow_choice=choices(i,:); 
            end 
        end 
    end 
end   

This code checks if the monomer at ‘idx’ is a corner monomer and if so, considers the four choices for the new position. It then checks if the two neighbouring monomers are one unit away from the new position and the new position is not already in the polymer. If there is an allowed choice, ‘tpoly’ is updated accordingly.

With this modification, the polymer simulation should continue to run even after all the end moves have been completed.

Hope this helps!

Regards,

Rasmita