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

1 次查看（过去 30 天）
Vaibhav
Vaibhav 2023-2-22
评论： Askic V 2023-2-22
% 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'.
  1 个评论
Askic V
Askic V 2023-2-22
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.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

Anton Kogios
Anton Kogios 2023-2-22
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)
Torsten
Torsten 2023-2-22
Ran in:
% 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

类别

在 Help Center 和 File Exchange 中查找有关 Physics 的更多信息

标签

产品


版本

R2022b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by