1
$\begingroup$

In my Monte Carlo simulations, hard particles aggregate into sets of triple particles due to the energy minimization. The particles cannot overlap, but they touch at the minimum of the energy. The potential between them is a pair potential which depends on the distance between the pairs and orientation of each particles.

The orientation is measured based on the angle between the orientation vector of each particle and the line which connects the centre of them.

I do understand how the Monte Carlo works, but I would like to use Mathematica to find the minimum energy of only 3 particles besides the relevant parameters like orientation angles. From the Monte Carlo I know they form equilateral triangle so it is possible to fix the distance by putting the particles on the vertices of an equilateral triangle. Then the energy will be minimized by changing the orientation angles.

How can I do such constraint minimization in Mathematica?

I thought something like what happen in Monte Carlo simulation, iterative minimization, should work, but I could not put it together correctly.

Could you please let me know how should I handle it correctly?

Below you can find a minimal version of the problem. The radius of each particle is R=10 which makes the minimum distance, r= 20. The goal is to minimize v12+v13+v23 and finding B1, B2, B3.

Thank you.

    Clear["Global`*"]
(*r12=r13=r23= r ; the distance between the particles. The goal is to minimize v12+v13+v23*)

v12[B1_ , B2_] = 
  10^6 r^-3 E^-r  (-r^2 Cos[B1] Cos[B2] + r Sin[B1] Sin[B2]) ;
v13[B1_, B3_] = 
  10^6 r^-3 E^-r  (-r^2 Cos[B1] Cos[B3] + r Sin[B1] Sin[B3]) ;
v23[B2_, B3_] = 
  10^6 r^-3 E^-r  (-r^2 Cos[B2] Cos[B3] + r Sin[B2] Sin[B3]) ;

R = 10.0;
r = 20.0;
b1 = 0; b2 = 0; b3 = 0; aa0 = 10^4; aa = 0; a1 = 0; a2 = 0; a3 = 0;
{a1, {b1 , b2}} = 
 While[Abs[aa - aa0] > 10^-5, aa0 = aa; 
  NMinimize[
    v12[B1 , 
     B2 ], {{B1, b1 - 0.1, b1 + 0.1}, {B2, b2 - 0.1, b2 + 0.1}}, 
    AccuracyGoal -> 20, PrecisionGoal -> 18, 
    WorkingPrecision -> 10] /. sol : {__Rule} :> Values[sol];
  {a2, {b1 , b3}} = 
   NMinimize[
     v13[B1 , 
      B3], {{B1, b1 - 0.1, b1 + 0.1}, {B3, b3 - 0.1, b3 + 0.1}}, 
     AccuracyGoal -> 20, PrecisionGoal -> 18, 
     WorkingPrecision -> 10] /. sol : {__Rule} :> Values[sol];
  {a3, {b2 , b3}} = 
   NMinimize[
     v23[B2, B3], {{B2, b2 - 0.1, b2 + 0.1}, {B3, b3 - 0.1, 
       b3 + 0.1}}, AccuracyGoal -> 20, PrecisionGoal -> 18, 
     WorkingPrecision -> 10] /. sol : {__Rule} :> Values[sol];
  aa = a1 + a2 + a3;]
Improve this question
$\endgroup$

1 Answer 1

Reset to default
2
$\begingroup$

I am not sure if I understand you correctly. You want to minimize V12+V13+V23 over B1, B2, B3. What is the `While`` loop for? And as V12,V13,V23 are not independent, you must minimize them together, not separately.

Anyway, this can be done by:

r = 20.0;
v12[B1_, B2_] = 
  10^6 r^-3 E^-r (-r^2 Cos[B1] Cos[B2] + r Sin[B1] Sin[B2]);
v13[B1_, B3_] = 
  10^6 r^-3 E^-r (-r^2 Cos[B1] Cos[B3] + r Sin[B1] Sin[B3]);
v23[B2_, B3_] = 
  10^6 r^-3 E^-r (-r^2 Cos[B2] Cos[B3] + r Sin[B2] Sin[B3]);

Minimize[{v12[B1, B2] + v13[B1, B3] + v23[B2, B3], 
  0 <= {B1, B2, B3} < 2 Pi}, {B1, B2, B3}, Reals]

(* {-0.000309173, {B1 -> 3.14159, B2 -> 3.14159, B3 -> 3.14159}} *)
Improve this answer
$\endgroup$
4
  • $\begingroup$ Thank you for your reply. Actually, I have tried that, but it does not work. Specially in this case the answer seems to align particles on a line; head-to-tail orientation, but I was looking for the case where I put them on an equilateral triangle. This simplified version of the potential is similar to the dipole-dipole interaction where the head-to-tail is preferred for free particles. I though by using while loop I can mimic something similar to the Monte Carlo method. $\endgroup$
    – Aa Aa
    Commented Mar 27, 2021 at 16:36
  • $\begingroup$ Why then are you not calculating the sum of the pair potentials for the molecules sitting at the vertices of an equilateral triangle and minimize over B1, B2, B3? $\endgroup$
    – Daniel Huber
    Commented Mar 27, 2021 at 16:48
  • $\begingroup$ Unfortunately, the formalism that I am using has a simple form if I measure the angles relative to the connecting line between the centre of the pair particles. Converting this formalism to a more general case to measure the angles relative to the xyz axes makes the equations very complicated which are not usable in Monte Carlo simulations. It means rewriting the pair interaction based on the coordinates of the particles is not easy. I thought if I iterate over the angles it may help, but I do not have any clue. $\endgroup$
    – Aa Aa
    Commented Mar 27, 2021 at 16:57
  • $\begingroup$ The pair potential can not depend on the absolute orientation. Only the relative position matters. The angles relative to the connecting line will be o.k. $\endgroup$
    – Daniel Huber
    Commented Mar 27, 2021 at 17:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.