MathGroup Archive 2010

[Date Index] [Thread Index] [Author Index]

Search the Archive

Trouble with coupled quadratic equations where the solutions are degenerate/symmetric

  • To: mathgroup at smc.vnet.net
  • Subject: [mg106476] Trouble with coupled quadratic equations where the solutions are degenerate/symmetric
  • From: Robert Hoy <robert.hoy at yale.edu>
  • Date: Thu, 14 Jan 2010 05:46:47 -0500 (EST)

Hi, all!  I am trying to solve for the structures of small atomic  
clusters using Mathematica.  Using adjacency (contact) matrices as  
input, my the desired output is 3-D Cartesian coordinates.   I  
formulate the problem as coupled quadratic equations and  
inequalities.  For example, for clusters of M = 6 atoms I introduce  
the variables

xvec = {{x1, y1, z1}, {x2, y2, z2}, {x3, y3, z3}, {x4, y4, z4}, {x5, y5,
    z5}, {x6, y6, z6}}

Then, for an adjacency matrix A

{
   {0, 1, 0, 1, 1, 1},
   {1, 0, 1, 1, 0, 1},
   {0, 1, 0, 1, 1, 1},
   {1, 1, 1, 0, 1, 0},
   {1, 0, 1, 1, 0, 1},
   {1, 1, 1, 0, 1, 0}
  }

there are M(M-1)/2 coupled equations and inequalities:

EqnTab = {(x1 - x2)^2 + (y1 - y2)^2 + (z1 - z2)^2 ==
   1, (x1 - x3)^2 + (y1 - y3)^2 + (z1 - z3)^2 >
   1, (x1 - x4)^2 + (y1 - y4)^2 + (z1 - z4)^2 ==
   1, (x1 - x5)^2 + (y1 - y5)^2 + (z1 - z5)^2 ==
   1, (x1 - x6)^2 + (y1 - y6)^2 + (z1 - z6)^2 ==
   1, (x2 - x3)^2 + (y2 - y3)^2 + (z2 - z3)^2 ==
   1, (x2 - x4)^2 + (y2 - y4)^2 + (z2 - z4)^2 ==
   1, (x2 - x5)^2 + (y2 - y5)^2 + (z2 - z5)^2 >
   1, (x2 - x6)^2 + (y2 - y6)^2 + (z2 - z6)^2 ==
   1, (x3 - x4)^2 + (y3 - y4)^2 + (z3 - z4)^2 ==
   1, (x3 - x5)^2 + (y3 - y5)^2 + (z3 - z5)^2 ==
   1, (x3 - x6)^2 + (y3 - y6)^2 + (z3 - z6)^2 ==
   1, (x4 - x5)^2 + (y4 - y5)^2 + (z4 - z5)^2 ==
   1, (x4 - x6)^2 + (y4 - y6)^2 + (z4 - z6)^2 >
   1, (x5 - x6)^2 + (y5 - y6)^2 + (z5 - z6)^2 == 1}


In these equations and inequalities, the left hand sides are  
Dot[xvec[[m]] - xvec[[n]], xvec[[m]] - xvec[[n]] (i. e. the squared  
Euclidean distance between xvec[[m]] and xvec[[n]] and the 'right hand  
sides'  are '==1' if A[[m,n]]==1 and '> 1' if A[[m,n]]=0.  Only terms  
above the diagonal of the adjacency matrix (n > m) are necessary since  
A is symmetric.

I know the solution to these equations; the six elements of xvec form  
the vertices of an octahedron.

Mathematica hangs when I try to solve this set.  This isn't  
surprising, because the above equations alone are invariant under  
translations and rotations, so the solutions are continuously  
degenerate.  To eliminate the continuous degeneracies I add equations  
constraining the first three 'atoms' to a plane (I can do this because  
they are noncollinear).  There is another discrete mirror degeneracy  
(the other atoms can lie on either side of the defined plane), so I  
add another inequality saying the fourth atom is on one side.  This  
adds ten equations/inequalities and gives the system

CompleteEqnTab = {x1 == 0, y1 == 0, z1 > 0, x2 == 0, z2 == 0, y2 > 0,  
y3 == 0, z3 == 0, x3 > 0,
   x4/x3 + y4/y2 + z4/z1 > 1, [the 15 terms of EqnTab]}


Now, the problem is, I have been unable (after a considerable amount  
of effort) to solve this system with Mathematica.   
Reduce[CompleteEqnTab, Flatten[xvec], Reals] fails.  So does  
FindInstance[CompleteEqnTab, Flatten[xvec], Reals].  The memory used  
by the kernel just increases (well past 1GB) until the program hangs.

I think the problem is that the solution for xvec (i. e. the vertices  
of an octahedron) has many rotational symmetries - its point group is  
O_h.  Perhaps these symmetries give Mathematica's solver trouble?

In contrast, the exact same approach for a less symmetric cluster with  
adjacency matrix

(
   {0, 1, 1, 1, 1, 1},
   {1, 0, 1, 1, 1, 1},
   {1, 1, 0, 1, 0, 0},
   {1, 1, 1, 0, 1, 0},
   {1, 1, 0, 1, 0, 1},
   {1, 1, 0, 0, 1, 0}
  )

and equation/inequality set

CompleteEqnTab2 = {x1 == 0, y1 == 0, z1 > 0, x2 == 0, z2 == 0, y2 > 0,  
y3 == 0, z3 == 0, x3 > 0,
   x4/x3 + y4/y2 + z4/z1 > 1, (x1 - x2)^2 + (y1 - y2)^2 + (z1 - z2)^2 ==
   1, (x1 - x3)^2 + (y1 - y3)^2 + (z1 - z3)^2 ==
   1, (x1 - x4)^2 + (y1 - y4)^2 + (z1 - z4)^2 ==
   1, (x1 - x5)^2 + (y1 - y5)^2 + (z1 - z5)^2 ==
   1, (x1 - x6)^2 + (y1 - y6)^2 + (z1 - z6)^2 ==
   1, (x2 - x3)^2 + (y2 - y3)^2 + (z2 - z3)^2 ==
   1, (x2 - x4)^2 + (y2 - y4)^2 + (z2 - z4)^2 ==
   1, (x2 - x5)^2 + (y2 - y5)^2 + (z2 - z5)^2 ==
   1, (x2 - x6)^2 + (y2 - y6)^2 + (z2 - z6)^2 ==
   1, (x3 - x4)^2 + (y3 - y4)^2 + (z3 - z4)^2 ==
   1, (x3 - x5)^2 + (y3 - y5)^2 + (z3 - z5)^2 >
   1, (x3 - x6)^2 + (y3 - y6)^2 + (z3 - z6)^2 >
   1, (x4 - x5)^2 + (y4 - y5)^2 + (z4 - z5)^2 ==
   1, (x4 - x6)^2 + (y4 - y6)^2 + (z4 - z6)^2 >
   1, (x5 - x6)^2 + (y5 - y6)^2 + (z5 - z6)^2 == 1}

gives a quick solution

In[32]:= Timing[Reduce[CompleteEqnTab2, Flatten[xvec], Reals]]

Out[32]= {18.527, x1 == 0 && y1 == 0 && z1 == 1/Sqrt[2] && x2 == 0 &&  
y2 == 1/Sqrt[2] &&
    z2 == 0 && x3 == 1/Sqrt[2] && y3 == 0 && z3 == 0 && x4 == 1/ 
Sqrt[2] &&
   y4 == 1/Sqrt[2] && z4 == 1/Sqrt[2] && x5 == -(1/(3 Sqrt[2])) &&
   y5 == (2 Sqrt[2])/3 && z5 == (2 Sqrt[2])/3 && x6 == -(11/(9  
Sqrt[2])) &&
   y6 == 5/(9 Sqrt[2]) && z6 == 5/(9 Sqrt[2])}

Similar issues hold for larger clusters; Mathematica solves the less- 
symmetric systems and hangs on the more-symmetric ones.  Is there  
anything I can do about this?  I thought FindInstance was designed for  
this sort of thing but it doesn't seem to help.

As you can see from the above, the form of the equations and  
inequalities is simple.  I assume Mathematica can handle the symmetric  
cases and I'm just unaware of the necessary functionality.  Any help  
would be greatly appreciated.  I am using Mathematica 7.0.1.0 on Mac  
OS 10.5.8.

Thanks,
Robert Hoy


  • Prev by Date: Re: More /.{I->-1} craziness
  • Next by Date: Re: Re: More /.{I->-1} craziness
  • Previous by thread: Re: Labeling points on ListPlot - Transform string date element in list of tripes
  • Next by thread: Re: Trouble with coupled quadratic equations where the