Unrecognized function or variable 'calculate_energy' and 'move_monomer'. how to solve the issue?

% Define initial energy and energy change when covalent interaction breaks
total_energy = -8.75;
delta_energy = -0.25;

% Define temperature and number of iterations
kT = 1;
nsteps = 100000;

% Define initial state and calculate its energy
state = [10 10; 9 10; 8 10; 7 10; 7 9; 8 9; 9 10; 10 9];
energy = calculate_energy(state, [-0.25 -0.25 -0.25 -0.25 -1 -1 -1 -1]);

% Store accepted moves for movie generation
accepted_moves = zeros(nsteps,2,2);

% Run Metropolis Criterion simulation
for step = 1:nsteps
    % Randomly select a monomer to move
    monomer_index = randi(16);
    
    % Move the selected monomer
    new_state = move_monomer(state, monomer_index, 1, 90);
    
    % Calculate the energy change
    new_energy = calculate_energy(new_state, [-0.25 -0.25 -0.25 -0.25 -1 -1 -1 -1]);
    delta_E = new_energy - energy;
    
    % Accept or reject the move
    if delta_E <= 0 || exp(-delta_E / (kT)) > rand()
        state = new_state;
        energy = new_energy;
        accepted_moves(step,:,:) = [state(monomer_index,:); new_state(monomer_index,:)];
    else
        accepted_moves(step,:,:) = [state(monomer_index,:); state(monomer_index,:)];
    end
end

% Generate movie of accepted moves
figure;
for step = 1:nsteps
    plot(state(:,1), state(:,2), 'bo-', 'MarkerSize', 8, 'LineWidth', 2);
    axis([4 13 4 13]);
    hold on;
    move = squeeze(accepted_moves(step,:,:));
    plot(move(:,1), move(:,2), 'r-', 'LineWidth', 2);
    hold off;
    drawnow;
end




function energy = calculate_energy(coords, energy_vals)
energy = 0;
for i = 1:size(coords,1)-1
    for j = i+1:size(coords,1)
        dist = norm(coords(i,:) - coords(j,:));
        energy = energy + energy_vals(i,j) / dist;
    end
end
end

function new_coords = move_monomer(coords, index, displacement, angle)
% Get the coordinates of the monomer to move and its neighbors
monomer_coords = coords(index,:);
if index == 1
    neighbor_coords = coords(index+1,:);
elseif index == size(coords,1)
    neighbor_coords = coords(index-1,:);
else
    neighbor_coords = [coords(index-1,:); coords(index+1,:)];
end

% Calculate the direction of the displacement
dir_vec = sum(neighbor_coords) / size(neighbor_coords,1) - monomer_coords;
dir_vec = dir_vec / norm(dir_vec);

% Rotate the displacement vector by the specified angle
rot_mat = [cosd(angle) -sind(angle); sind(angle) cosd(angle)];
dir_vec = rot_mat * dir_vec';

% Calculate the new coordinates
new_coords = coords;
new_coords(index,:) = monomer_coords + displacement * dir_vec';
end

Unrecognized function or variable 'calculate_energy'.

[Askic V]

First of all this line will generate an error:

energy = energy + energy_vals(i,j) / dist;

in the function call you send an array with dimensions 1x8, and inside nested loop you're trying to access the element with indices i = 2, j = 3 i.e. energy_vals(i,j). You need to rethink what you really want here.

Answer 1

[Anton Kogios]

Your two functions have to be saved in the same directory as the script you are running, and with the name of the file being the name of the function.

Alternatively, you can have the functions at the very end of the script you are running. (i.e. if I copy-paste the code in your question into one long script, it runs fine, albeit with some index dimensions errors which I think you should be able to solve)

Answer 2

[Torsten]

% Define initial energy and energy change when covalent interaction breaks
total_energy = -8.75;
delta_energy = -0.25;

% Define temperature and number of iterations
kT = 1;
nsteps = 100000;

% Define initial state and calculate its energy
state = [10 10; 9 10; 8 10; 7 10; 7 9; 8 9; 9 10; 10 9];
energy_vals = rand(size(state,1));
energy = calculate_energy(state, energy_vals);%[-0.25 -0.25 -0.25 -0.25 -1 -1 -1 -1]);

% Store accepted moves for movie generation
accepted_moves = zeros(nsteps,2,2);

% Run Metropolis Criterion simulation
for step = 1:nsteps
    % Randomly select a monomer to move
    monomer_index = randi(size(state,1));%randi(16);
    
    % Move the selected monomer
    new_state = move_monomer(state, monomer_index, 1, 90);
    
    % Calculate the energy change
    new_energy = calculate_energy(new_state, energy_vals);%[-0.25 -0.25 -0.25 -0.25 -1 -1 -1 -1]);
    delta_E = new_energy - energy;
    
    % Accept or reject the move
    if delta_E <= 0 || exp(-delta_E / (kT)) > rand()
        state = new_state;
        energy = new_energy;
        accepted_moves(step,:,:) = [state(monomer_index,:); new_state(monomer_index,:)];
    else
        accepted_moves(step,:,:) = [state(monomer_index,:); state(monomer_index,:)];
    end
end

% Generate movie of accepted moves
figure;
for step = 1:nsteps
    plot(state(:,1), state(:,2), 'bo-', 'MarkerSize', 8, 'LineWidth', 2);
    axis([4 13 4 13]);
    hold on;
    move = squeeze(accepted_moves(step,:,:));
    plot(move(:,1), move(:,2), 'r-', 'LineWidth', 2);
    hold off;
    drawnow;
end

function energy = calculate_energy(coords, energy_vals)
energy = 0;
for i = 1:size(coords,1)-1
    for j = i+1:size(coords,1)
        dist = norm(coords(i,:) - coords(j,:));
        energy = energy + energy_vals(i,j) / dist;
    end
end
end

function new_coords = move_monomer(coords, index, displacement, angle)
% Get the coordinates of the monomer to move and its neighbors
monomer_coords = coords(index,:);
if index == 1
    neighbor_coords = coords(index+1,:);
elseif index == size(coords,1)
    neighbor_coords = coords(index-1,:);
else
    neighbor_coords = [coords(index-1,:); coords(index+1,:)];
end

% Calculate the direction of the displacement
dir_vec = sum(neighbor_coords) / size(neighbor_coords,1) - monomer_coords;
dir_vec = dir_vec / norm(dir_vec);

% Rotate the displacement vector by the specified angle
rot_mat = [cosd(angle) -sind(angle); sind(angle) cosd(angle)];
dir_vec = rot_mat * dir_vec';

% Calculate the new coordinates
new_coords = coords;
new_coords(index,:) = monomer_coords + displacement * dir_vec';
end