Need help with calculating the minimum potential energy of two particles.

Hello, I am attempting to find the minimum potential between particles. One of the particles is fixed at the point (0,0) and the other particle is allowed to move freely and I have set that particle at an arbitrary point, in this case (5,5). I am trying to use the Monte Carlo method of finding a minimum potential energy between the two particles. The free particle is supposed to slightly move around and slightly vary its angle theta and these new values are then used to calculate the new potential energy of the system. If the potential energy of the system is lower then these new values are saved and the calculation is repeated with the new values. However, if the potential energy increases then these new values are discarded and a new calculation is carried out so far. Here is what I have so far,

%Variable declaration

epsilon = (8.85*10^-12); %permittivity of free space
constant = 1/(4*pi*epsilon); %form of the coulomb constant

%position initialization of the fixed particle 

particle.x = 0; 
particle.y = 0; 
dipole.angle.fixed.particle = pi/2;

%position initialization of the free particle 

particle.x1 = 5; %set at arbitrary position
particle.y1 = 5; %set at arbitrary position
particle.theta = pi/6 

%calculation of initial energy
initial.potential.energy = ((-1)*constant*1*1*(cos(dipole.angle.fixed.particle - particle.theta)))./(((particle.x1 - particle.x).^2 + (particle.y1 - particle.y).^2).^(3/2))

for i = 1:100
 deltaX = particle.x1 + 0.1*randn(1) - 0.1; 
 deltaY = particle.y1 + 0.1*randn(1) - 0.1;
 theta = particle.theta + randn(1);
 potential.energy(i) = ((-1)*constant*1*1*(cos(dipole.angle.fixed.particle - particle.theta)))./(((deltaX - particle.x).^2 + (deltaY - particle.y).^2).^(3/2))
 if (potential.energy(i) < initial.potential.energy)
     particle.x1 = deltaX; 
     particle.y1 = deltaY;
     particle.theta = theta;
     initial.potential.energy = potential.energy
 end

end

I believe my problem is that the if condition in my for loop isn't activating and therefore the values that lead to a lower potential energy aren't being saved. Instead the program is doing the calculation with whatever values it wants. I am unsure how to fix this error or where the error in my code is and any help at all would be appreciated to fix this issue and to get this program actually working.

Answer 1

[nl2605]

Hi there! Do you really want structures in your program? If not, then you should remove the punkt in between.

%Variable declaration

epsilon = (8.85*10^-12); %permittivity of free space

constant = 1/(4*pi*epsilon); %form of the coulomb constant

%position initialization of the fixed particle

particle_x = 0;

particle_y = 0;

dipole_angle_fixed_particle = pi/2;

%position initialization of the free particle

particle_x1 = 5; %set at arbitrary position

particle_y1 = 5; %set at arbitrary position

particle_theta = pi/6 ;

%calculation of initial energy

initial_potential_energy = ((-1)*constant*1*1(cos(dipole_angle_fixed_particle - particle_theta)))./(((particle_x1 - particle_x).^2 + (particle_y1 - particle_y).^2).^(3/2));

for i = 1:100

 deltaX = particle_x1 + 0.1*randn(1) - 0.1;

 deltaY = particle_y1 + 0.1*randn(1) - 0.1;

 theta = particle_theta + randn(1);

 potential_energy(i) = ((-1)*constant*1*1*(cos(dipole_angle_fixed_particle - particle_theta)))./(((deltaX - particle_x).^2 + (deltaY - particle_y).^2).^(3/2));

 if (potential_energy(i) < initial_potential_energy)

     particle_x1 = deltaX; 

     particle_y1 = deltaY;

     particle_theta = theta;

     initial_potential_energy = potential_energy;

 end

end

If you do not need structures this should work to give you the potential energy values. I hope I understood it right.

[Gary]

Hello! I appreciate the help but I'm still not getting the correct output. that I want. For some reason the potential energy keeps going up and down and isn't converging to a minimum. The problem is that once the potential energy finds a lower, more stable value, it should only be able to take on even lower values. As it stands right now the potential energy will sometimes increase as well which shouldn't be happening.

[nl2605]

What do you want the program to do when the potential energy is greater than the current potential energy?

[Gary]

Ideally I would like that potential energy value to be discarded or thrown out and a new calculation to be carried out. The only potential energy values that should be allowed should be values that are lower than the previous potential energy.

[nl2605]

Hmmm..I don't really see any problem with the code. Maybe you have used random function in the code that is creating the problem. Because I changed the program with an else statement after the if statement in else I put potential energy = initial potential energy. It still seems to be increasing.

[Gary]

Hello, thank you so much for the help. I found out the problem I was having is that my if condition isn't being executed while the loop iterates which is why I was getting incorrect values. Putting in an else statement and changing the final line to potential_energy(i) = initial_potential_energy; solved the problem.