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

**Follow-Ups**:**Re: SymbolicPolynomialMod***From:*danl@wolfram.com

**Re: Trouble with coupled quadratic equations where the***From:*Daniel Lichtblau <danl@wolfram.com>