Unrecognized function or variable 'calculate_energy' and 'move_monomer'. how to solve the issue?
1 visualizzazione (ultimi 30 giorni)
Mostra commenti meno recenti
% 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 Commento
Askic V
il 22 Feb 2023
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.
Risposte (2)
Anton Kogios
il 22 Feb 2023
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)
0 Commenti
Torsten
il 22 Feb 2023
% 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
0 Commenti
Vedere anche
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!